Homec4science
Diffusion SRR - Source receptor relationship (master) Browse

srr/training.cc
900b64533547master

Last Change
Raw File

training.cc
View Options

#include "training.hh"
#include <math.h>
#include <limits>
#include <algorithm>
#include "dlib-18.12/dlib/optimization.h"
#include "dlib-18.12/dlib/optimization/optimization.h"
#include "dlib-18.12/dlib/matrix/matrix_abstract.h"
#include "dlib-18.12/dlib/statistics/statistics_abstract.h"
#include "dlib-18.12/dlib/statistics.h"
#include "dlib-18.12/dlib/clustering.h"
Training::Training(ReadScenarios * input_inf) : input_info(input_inf){
SingleVectorD dim_PreAlpha(this->input_info->get_nPre(),0.00001);
SingleVectorD dim_PreOmega(this->input_info->get_nPre(),-2.0);
DoubleVectorD AlphaCoeff_Step1_(this->input_info->get_nPol(),dim_PreAlpha);
DoubleVectorD OmegaCoeff_Step1_(this->input_info->get_nPol(),dim_PreOmega);
this->input_info->set_AlphaCoeff_Step1(AlphaCoeff_Step1_);
this->input_info->set_OmegaCoeff_Step1(OmegaCoeff_Step1_);
SingleVectorD dim_Ny_Alpha(this->input_info->get_ny(),0.00001);
SingleVectorD dim_Ny_Omega(this->input_info->get_ny(),-2);
DoubleVectorD dim_Nx_Alpha(this->input_info->get_nx(),dim_Ny_Alpha);
DoubleVectorD dim_Nx_Omega(this->input_info->get_nx(),dim_Ny_Omega);
ThreeVectorD dim_nPre_Alpha(this->input_info->get_nPre(), dim_Nx_Alpha);
ThreeVectorD dim_nPre_Omega(this->input_info->get_nPre(), dim_Nx_Omega);
FourVectorD AlphaCoeff_Step2_(this->input_info->get_nPol(), dim_nPre_Alpha); // Vector containing the Alpha coefficient resulting from Step2 (Ordered by precursos, size = (nPol, nPre, nx, ny)
FourVectorD OmegaCoeff_Step2_(this->input_info->get_nPol(), dim_nPre_Omega); // Vector containing the Omega coefficient resulting from Step2 (Ordered by precursos, size = (nPol, nPre, nx, ny)
this->input_info->set_AlphaCoeff_Step2(AlphaCoeff_Step2_);
this->input_info->set_OmegaCoeff_Step2(OmegaCoeff_Step2_);
SingleVectorD d1(this->input_info->get_nPre(),1.);
DoubleVectorD dist_influence(this->input_info->get_nPol(),d1); //(nPoll, nPre)
this->Dist_influence = dist_influence;
}
void Training::ClassifyScenarios_Step1(){
UInt Nx = input_info->get_nx();
UInt Ny = input_info->get_ny();
SingleVectorD dim1(Ny,0);
DoubleVectorD dim2(Nx,dim1);
// Verifying if exist scenario available by precursor
std::cout << " > Verifying the existence of a scenario by precursor (emissions reductions)" << std::endl;
FourVectorD DETra_ = input_info->get_DETra(); //(nTra, nPre, nx, ny)
UInt nPrec_ = input_info->get_nPre();
SingleVectorS namePrec_ = input_info->get_NamePrecursors();
std::vector<int> validScn_Prec(nPrec_,-1);
for (UInt p=0; p<nPrec_; ++p){
int ok_p=-1;
for (UInt s=0; s<DETra_.size(); ++s){
int ok_p1=-1;
for (UInt p1=0; p1<nPrec_; ++p1){
if (p==p1){
if (dim2 != DETra_[s][p]){ok_p=1;}
}
else{
if (dim2 != DETra_[s][p1]){ok_p1=1;}
}
}
if (ok_p==1 && ok_p1==-1) {
validScn_Prec[p]=s;
}
}
}
for (UInt i=0; i<validScn_Prec.size(); ++i){
if (validScn_Prec[i] > 0) std::cout << " For " << namePrec_[i] << ": scenario " << validScn_Prec[i] << " is ok" <<std::endl;
else {std::cout << " For " << namePrec_[i] << ": there is no adequate scenario" << std::endl;
}
}
this->validScn_ByPrec = validScn_Prec;
}
double r;
double Lat_resolution;
double getR(){return r;}
void setR(double r_){r = r_;}
void set_LatResol(double lat_resol_){Lat_resolution = lat_resol_;}
static bool IsInBounds (double value){
double r = getR();
double low = r-Lat_resolution/2;
double high = r+Lat_resolution/2;
return !(value < low) && (value < high);
}
void Training::DefiningDistInfluent_Step1(){
std::cout << " > Defining distance of influence by precursor" << std::endl;
SingleVectorD Lat_ = this->input_info->get_Lat_BC();
SingleVectorD Lon_ = this->input_info->get_Lon_BC();
UInt nPrec_ = input_info->get_nPre();
UInt nPol_ = input_info->get_nPol();
//UInt Nx_ = input_info->get_nx();
//UInt Ny_ = input_info->get_ny();
SingleVectorS namePrec_ = input_info->get_NamePrecursors();
SingleVectorS namePoll_ = input_info->get_NamePollutants();
FourVectorD DETra_ = input_info->get_DETra(); //(nTra, nPre, nx, ny)
FourVectorD DCTra_ = input_info->get_DCTra(); //(nTra, nPol, nx, ny)
std::vector<int> validScn_Prec = this->validScn_ByPrec;
double Lat_resolution = sqrt(pow(Lat_[1]-Lat_[0],2));
double Lon_resolution = sqrt(pow(Lon_[1]-Lon_[0],2));
set_LatResol(Lat_resolution);
double dist_resol = (Lat_resolution + Lon_resolution)/2;
double Lat_min = Lat_[0];
double Lat_max = Lat_[Lat_.size()-1];
double maxDist_Lat = sqrt(pow(Lat_[0]-Lat_[Lat_.size()-1],2));
double maxDist_Lon = sqrt(pow(Lon_[0]-Lon_[Lon_.size()-1],2));
double maxDist = sqrt(pow(maxDist_Lon,2)+pow(maxDist_Lat,2));
double ddd =0;
UInt ii=0;
SingleVectorD dist_vector;
while (ddd<maxDist){
ddd = dist_resol *pow(ii,2);
dist_vector.push_back(ddd);
ii=ii+2;
}
SingleVectorD d1(nPrec_,0);
DoubleVectorD dist_influence(nPol_,d1); //(nPoll, nPre)
//std::cout << "aqui " << validScn_Prec.size() << std::endl;
for (UInt poll=0; poll<nPol_; ++poll){
for (UInt pre=0; pre<nPrec_; ++pre){
int scn = validScn_Prec[pre];
std::cout << " For Poll= " << namePoll_[poll] << " and Prec= "<< namePrec_[pre] << "; dist_influenc = " << std::endl;;
if (scn >=0){
DoubleVectorD DEScn_prec = DETra_[scn][pre]; //this vector contains the deltaEmissions of the correct scenario respect to a given precursor
DoubleVectorD DCScn_prec = DCTra_[scn][poll]; //this vector contains the deltaConcentrations of the correct scenario for a given pollutant
double dist_stop;
double corr_coef = 1.;
for (UInt d=0; d<dist_vector.size(); ++d){
if (corr_coef > 0.09){
if (d>0) dist_stop = dist_vector[d-1];
else dist_stop = dist_vector[1];
double dist = dist_vector[d];
SingleVectorD DE_dist;
SingleVectorD DC_dist;
for (UInt yc=0; yc<Lat_.size(); ++yc){
for (UInt xc=0; xc<Lon_.size(); ++xc){
for (UInt xe=0; xe<Lon_.size();++xe){
double a_=1;
double b_= -2*Lat_[yc];
double c_= pow(Lon_[xc],2)+pow(Lon_[xe],2)+pow(Lat_[yc],2)-pow(dist,2)-2*Lon_[xc]*Lon_[xe];
double determinant = pow(b_,2)-4*a_*c_;
float r1=-1,r2=-1;
UInt ye;
if (determinant >=0){
r1= ((-b_+sqrt(determinant))/2*a_);
r2= ((-b_-sqrt(determinant))/2*a_);
if (r1>=Lat_min && r1<=Lat_max){
setR(r1);
ye = std::find_if(Lat_.begin(),Lat_.end(),IsInBounds)-Lat_.begin();
if (DEScn_prec[xe][ye] !=0.){
DE_dist.push_back(DEScn_prec[xe][ye]);
DC_dist.push_back(DCScn_prec[xc][yc]);
}
}
if ((r2>=Lat_min && r2<=Lat_max) && r2 != r1){
setR(r2);
ye = std::find_if(Lat_.begin(),Lat_.end(),IsInBounds) -Lat_.begin();
if (DEScn_prec[xe][ye] !=0.){
DE_dist.push_back(DEScn_prec[xe][ye]);
DC_dist.push_back(DCScn_prec[xc][yc]);
}
}
}
}
}
}
//std::cout << dist << " " << DE_dist.size() << " " << DC_dist.size() << std::endl;
corr_coef = dlib::correlation(DE_dist,DC_dist);
//std::cout << corr_coef << std::endl;
}
}
dist_influence[poll][pre] = dist_stop;
std::cout << dist_stop << std::endl;
}
else{
//dist_influence[poll][pre] = maxDist;
//std::cout << " (maximum) " << maxDist << std::endl;
dist_influence[poll][pre] = 3.4;
std::cout << " (ojo fixing) " << 3.4 << std::endl;
}
}
}
this->Dist_influence = dist_influence;
}
void Training::ComputingLocalIMorans_Step1(){
std::cout << " > Computing spatial autocorrelation index (Morans's I), weighting criterion is distance influence" << std::endl;
SingleVectorD Lat_ = this->input_info->get_Lat_BC();
SingleVectorD Lon_ = this->input_info->get_Lon_BC();
UInt nPrec_ = input_info->get_nPre();
UInt nPol_ = input_info->get_nPol();
SingleVectorS namePrec_ = input_info->get_NamePrecursors();
SingleVectorS namePoll_ = input_info->get_NamePollutants();
std::vector<int> validScn_Prec = this->validScn_ByPrec;
FourVectorD DETra_ = input_info->get_DETra(); //(nTra, nPre, nx, ny)
UInt Nx_ = input_info->get_nx();
UInt Ny_ = input_info->get_ny();
int Nx__ = Nx_;
int Ny__ = Ny_;
DoubleVectorD Dist_influence = this->Dist_influence;
double Lat_resolution = sqrt(pow(Lat_[1]-Lat_[0],2));
double Lon_resolution = sqrt(pow(Lon_[1]-Lon_[0],2));
for (UInt poll=0; poll<nPol_; ++poll){
ThreeVectorD MoransI_byPrec;
for (UInt pre=0; pre<nPrec_; ++pre){
SingleVectorD DimIy(Ny_,0);
DoubleVectorD I_Morans(Nx_,DimIy);
int scn = validScn_Prec[pre];
if (scn >=0){
std::cout << " For "<< namePrec_[pre] << "; I ranges = ";
DoubleVectorD DEScn_prec = DETra_[scn][pre]; //this vector contains the deltaEmissions of the correct scenario respect to a given precursor
dlib::matrix<double> Ii; // Local Morans' I
Ii.set_size(Nx_,Ny_);
double dist_inf = Dist_influence[poll][pre]; // Weighting criteria
UInt x_cells = round(dist_inf/Lon_resolution);
UInt y_cells = round(dist_inf/Lat_resolution);
for (UInt xei=0; xei<Nx_; ++xei){
for (UInt yei=0; yei< Ny_; ++yei){
UInt Ncells=0;
dlib::matrix<double>DEScn_i;
DEScn_i.set_size(x_cells*2+1,y_cells*2+1);
DEScn_i=0;
dlib::matrix<double> Wi; // Weighting matrix
Wi.set_size(x_cells*2+1,y_cells*2+1);
Wi =0;
dlib::matrix<double> Zi; // deviation from the mean
Zi.set_size(x_cells*2+1,y_cells*2+1);
Zi=0;
int start_x = xei-x_cells; if (start_x<0) start_x=0;
int start_y = yei-y_cells; if (start_y<0) start_y=0;
int end_x = xei+x_cells; if (end_x>Nx__) end_x=Nx_-1;
int end_y = yei+y_cells; if (end_y>Ny__) end_y=Ny_-1;
for (int xek=start_x; xek<=end_x; ++xek){
for (int yek=start_y; yek<=end_y; ++yek){
double dist_loc = sqrt(pow(Lon_[xei]-Lon_[xek],2)+pow(Lat_[yei]-Lat_[yek],2));
if (dist_loc<=dist_inf){
DEScn_i(xek-start_x,yek-start_y) = DEScn_prec[xek][yek];
Wi(xek-start_x,yek-start_y) = 1;
++Ncells;
}
}
}
double ww= 1./Ncells;
Wi = ww*Wi;
double mean_ = dlib::mean(DEScn_i);
double variance_ = dlib::variance(DEScn_i);
Zi = DEScn_i-mean_;
double sum_ =(sum(dlib::pointwise_multiply(Zi,Wi))-(Zi(x_cells+1,y_cells+1)*Wi(x_cells+1,y_cells+1)));
if (variance_ !=0){
Ii(xei,yei) = Zi(x_cells+1,y_cells+1)/variance_ * sum_;
}
else Ii(xei,yei) = 0; // If variace equal to cero, No spatial dependence
I_Morans[xei][yei] = Ii(xei,yei);
}
}
std::cout << dlib::min(Ii) << "-" << dlib::max(Ii) << std::endl;
}
MoransI_byPrec.push_back(I_Morans);
}
this->MoransI.push_back(MoransI_byPrec);
}
}
double residual( const std::pair<double,DoubleVectorD> Pair_CellDC_DistEmi,
const dlib::matrix<double> & Alpha_Omega){
double res=0;
double DC_cell = Pair_CellDC_DistEmi.first;
SingleVectorD Dist_= (Pair_CellDC_DistEmi.second)[0];
SingleVectorD Emis_= (Pair_CellDC_DistEmi.second)[1];
double DC_cell_calculated=0;
for (UInt i=0; i<Dist_.size(); ++i){
DC_cell_calculated += Alpha_Omega(0) * Emis_[i] * (pow(1+Dist_[i],Alpha_Omega(1)));
}
res = pow((DC_cell_calculated-DC_cell),2);
//cout << res << " ";
return res;
}
dlib::matrix<double> residual_derivative( const std::pair<double,DoubleVectorD> Pair_CellDC_DistEmi,
const dlib::matrix<double> & Alpha_Omega){
dlib::matrix<double> der;
der.set_size(2,1);
der=0;
double DC_cell = Pair_CellDC_DistEmi.first;
SingleVectorD Dist_= (Pair_CellDC_DistEmi.second)[0];
SingleVectorD Emis_= (Pair_CellDC_DistEmi.second)[1];
double DC_cell_calculated=0;
for (UInt i=0; i<Dist_.size(); ++i){
DC_cell_calculated += Alpha_Omega(0) * Emis_[i] * (pow(1+Dist_[i],Alpha_Omega(1)));
}
double der_general = 2*(DC_cell-DC_cell_calculated);
double der_alpha=0;
for (UInt i=0; i<Dist_.size(); ++i){
der_alpha += (-( Emis_[i] * (pow(1+Dist_[i],Alpha_Omega(1)))));
}
der(0)= der_general* der_alpha;
double der_omega;
for (UInt i=0; i<Dist_.size(); ++i){
der_omega += Alpha_Omega(0) * Emis_[i] * (pow(1+Dist_[i],Alpha_Omega(1))) * log (1+Dist_[i]);
}
der(1)= der_general* der_omega;
//std::cout << "derivate " << der << std::endl;
return der;
}
class Residual{
public:
double operator() ( const std::pair<double,DoubleVectorD> Pair_CellDC_DistEmi,
const dlib::matrix<double> & Alpha_Omega
) const {
return residual (Pair_CellDC_DistEmi, Alpha_Omega);
};
};
void Training::BellShape_Step1_withMoransICriteria(){
std::cout << " > BellShape by precursor with Morans' I criteria " << std::endl;
FourVectorD ETra_ = input_info->get_ETra(); //(nTra, nPre, nx, ny)
FourVectorD DETra_ = input_info->get_DETra(); //(nTra, nPre, nx, ny)
FourVectorD DCTra_ = input_info->get_DCTra(); //(nTra, nPol, nx, ny)
SingleVectorS namePrec_ = input_info->get_NamePrecursors();
SingleVectorS namePoll_ = input_info->get_NamePollutants();
UInt Nx = input_info->get_nx();
UInt Ny = input_info->get_ny();
UInt nPrec_ = input_info->get_nPre();
UInt nPol_ = input_info->get_nPol();
SingleVectorD Lat_ = this->input_info->get_Lat_BC();
SingleVectorD Lon_ = this->input_info->get_Lon_BC();
std::vector<int> validScn_Prec = this->validScn_ByPrec;
double lat_resolution = sqrt(pow(Lat_[1]-Lat_[0],2));
double lon_resolution = sqrt(pow(Lon_[1]-Lon_[0],2));
SingleVectorD dim1(Ny,0);
DoubleVectorD dim2(Nx,dim1);
// Computing fitting.
std::cout << " - computing coefficients Alpha and Omega by precursor" << std::endl;
//double Alpha=1., Omega=-2.;
DoubleVectorD AlphaCoeff_Step1_ = this->input_info->get_AlphaCoeff_Step1();
DoubleVectorD OmegaCoeff_Step1_ = this->input_info->get_OmegaCoeff_Step1();
ThreeVectorD dim3(nPrec_,dim2);
FourVectorD AlphaCoeff_Step2_(nPol_,dim3); // Vector containing the Alpha coefficient resulting from Step2 (Ordered by precursos, size = (nPol, nPre, nx, ny)
FourVectorD OmegaCoeff_Step2_(nPol_,dim3); // Vector containing the Omega coefficient resulting from Step2 (Ordered by precursos, size = (nPol, nPre, nx, ny)
SingleVectorD Clust; Clust.push_back(-1000000); Clust.push_back(-1); Clust.push_back(-0.5); Clust.push_back(-0.25);
Clust.push_back(-0.00001); Clust.push_back(0.00001);Clust.push_back(0.25); Clust.push_back(0.5); Clust.push_back(1);Clust.push_back(1000000);
for (UInt poll=0; poll<nPol_; ++poll){
for (UInt pre=0; pre<nPrec_; ++pre){
int scn = validScn_Prec[pre];
if (scn >=0){
DoubleVectorD DEScn_prec = DETra_[scn][pre]; //this vector contains the deltaEmissions of the correct scenario respect to precursor
DoubleVectorD DCScn_prec = DCTra_[scn][poll]; //this vector contains the deltaEmissions of the correct scenario for a given pollutant
std::cout << " Coefficients for pollutant " << namePoll_[poll] << " from emissions of precursor " << namePrec_[pre] << ": " << std::endl;
// For computing a "Unique" by Precursor
/* std::vector<std::pair<double,DoubleVectorD>> Pair_CellDC_DistEmi;
for (UInt xc=0; xc<Nx; ++xc){
for (UInt yc=0; yc<Ny; ++yc){
SingleVectorD DistCell;
SingleVectorD Emiss;
DoubleVectorD DistEmi;
for (UInt xe=0; xe<Nx; ++xe){
for (UInt ye=0; ye<Ny; ++ye){
Emiss.push_back(DETra_[scn][pre][xe][ye]);
//DistCell.push_back(this->Dist[xc][yc][xe][ye]);
DistCell.push_back(sqrt(pow(Lon_[xe]-Lon_[xc],2) + pow(Lat_[ye]-Lat_[yc],2)));
}
}
DistEmi.push_back(DistCell);
DistEmi.push_back(Emiss);
Pair_CellDC_DistEmi.push_back(std::make_pair(DCTra_[scn][poll][xc][yc],DistEmi));
}
}
param_vector Alpha_Omega;
Alpha_Omega.set_size(2,1);
Alpha_Omega(0)=AlphaCoeff_Step1_[poll][pre];
Alpha_Omega(1)=OmegaCoeff_Step1_[poll][pre];
std::cout << " Coefficients for pollutant " << namePoll_[poll] << " from emissions of precursor " << namePrec_[pre] << ": " << std::endl;
dlib::solve_least_squares(dlib::objective_delta_stop_strategy(1e-7).be_verbose(),
residual,
residual_derivative,
Pair_CellDC_DistEmi,
Alpha_Omega);
AlphaCoeff_Step1_[poll][pre] = Alpha_Omega(0);
OmegaCoeff_Step1_[poll][pre] = Alpha_Omega(1);
std::cout << " Alpha:" << AlphaCoeff_Step1_[poll][pre] << " Omega:" << OmegaCoeff_Step1_[poll][pre] << std::endl;
for (UInt xc=0; xc<Nx; ++xc){
for (UInt yc=0; yc<Ny; ++yc){
AlphaCoeff_Step2_[poll][pre][xc][yc] = AlphaCoeff_Step1_[poll][pre];
OmegaCoeff_Step2_[poll][pre][xc][yc] = OmegaCoeff_Step1_[poll][pre];
}
}
*/
//COmputing Omega by precursor and by cluster according to Morans'I coefficient
int nX_ = Nx; int nY_= Ny;
double dist_inf = Dist_influence[poll][pre]; // Weighting criteria
UInt x_cells = round(dist_inf/lon_resolution);
UInt y_cells = round(dist_inf/lat_resolution);
for (UInt c=0; c<Clust.size()-1; ++c){
std::vector<std::pair<double,DoubleVectorD>> Pair_CellDC_DistEmi;
for (int xc=0; xc<nX_; ++xc){
for (int yc=0; yc<nY_; ++yc){
if (MoransI[poll][pre][xc][yc]>Clust[c] && MoransI[poll][pre][xc][yc]<=Clust[c+1]){
SingleVectorD DistCell;
SingleVectorD Emiss;
DoubleVectorD DistEmi;
int start_x = xc-x_cells; if (start_x<0) start_x=0;
int start_y = yc-y_cells; if (start_y<0) start_y=0;
int end_x = xc+x_cells; if (end_x>nX_) end_x=nX_-1;
int end_y = yc+y_cells; if (end_y>nY_) end_y=nY_-1;
for (int xc_t=start_x; xc_t<=end_x; ++xc_t){
for (int yc_t=start_y; yc_t<=end_y; ++yc_t){
double dist_loc = sqrt(pow(Lon_[xc]-Lon_[xc_t],2)+pow(Lat_[yc]-Lat_[yc_t],2));
if (dist_loc<=dist_inf){
Emiss.push_back(DETra_[scn][pre][xc_t][yc_t]);
DistCell.push_back(dist_loc);
}
}
}
DistEmi.push_back(DistCell);
DistEmi.push_back(Emiss);
Pair_CellDC_DistEmi.push_back(std::make_pair(DCTra_[scn][poll][xc][yc],DistEmi));
}
}
}
if (Pair_CellDC_DistEmi.size()>0){
param_vector Alpha_Omega;
Alpha_Omega.set_size(2,1);
Alpha_Omega(0)=AlphaCoeff_Step1_[poll][pre];
Alpha_Omega(1)=OmegaCoeff_Step1_[poll][pre];
std::cout << " Coefficients for pollutant " << namePoll_[poll] << " from emissions of precursor " << namePrec_[pre] << " in cluster I" << Clust[c] << "-" << Clust[c+1] <<" : " << std::endl;
std::cout << "sample size = " << Pair_CellDC_DistEmi.size() << std::endl;
//dlib::solve_least_squares(dlib::objective_delta_stop_strategy(1e-7).be_verbose(),
dlib::solve_least_squares(dlib::objective_delta_stop_strategy(1e-7),
residual,
residual_derivative,
Pair_CellDC_DistEmi,
Alpha_Omega);
for (UInt xc=0; xc<Nx; ++xc){
for (UInt yc=0; yc<Ny; ++yc){
if (this->MoransI[poll][pre][xc][yc]>Clust[c] && this->MoransI[poll][pre][xc][yc]<=Clust[c+1]){
AlphaCoeff_Step2_[poll][pre][xc][yc] = Alpha_Omega(0);
OmegaCoeff_Step2_[poll][pre][xc][yc] = Alpha_Omega(1);
}
}
}
std::cout << " Alpha:" << Alpha_Omega(0) << " Omega:" << Alpha_Omega(1) << std::endl;
}
}
//For computing an Omega by cell by precursor (take a lot of time ~ 2 seconds by cell by precursor)
/* for (int xc=0; xc<nX_; ++xc){
for (int yc=0; yc<nY_; ++yc){
std::vector<std::pair<double,DoubleVectorD>> Pair_CellDC_DistEmi;
for (int xc_t=xc-10; xc_t<xc+10; ++xc_t){ // t : it is mean "around"
if (xc_t>=0 && xc_t<nX_){
for (int yc_t=yc-10; yc_t<yc+10; ++yc_t){
if (yc_t>=0 && yc_t<nY_){
SingleVectorD DistCell;
SingleVectorD Emiss;
DoubleVectorD DistEmi;
for (int xe=xc_t-10; xe<xc_t+10; ++xe){
for (int ye=yc_t-10; ye<yc_t+10; ++ye){
if (xe>=0 && xe<nX_ && ye>=0 && ye<nY_){
Emiss.push_back(DETra_[scn][pre][xe][ye]);
DistCell.push_back(sqrt(pow(Lon_[xe]-Lon_[xc_t],2) + pow(Lat_[ye]-Lat_[yc_t],2)));
}
}
}
DistEmi.push_back(DistCell);
DistEmi.push_back(Emiss);
Pair_CellDC_DistEmi.push_back(std::make_pair(DCTra_[scn][poll][xc_t][yc_t],DistEmi));
}
}
}
}
param_vector Alpha_Omega;
Alpha_Omega.set_size(2,1);
Alpha_Omega(0)=AlphaCoeff_Step1_[poll][pre];
Alpha_Omega(1)=OmegaCoeff_Step1_[poll][pre];
//std::cout << " Coefficients for pollutant " << namePoll_[poll] << " from emissions of precursor " << namePrec_[pre] << ": " << std::endl;
//dlib::solve_least_squares(dlib::objective_delta_stop_strategy(1e-7).be_verbose(),
dlib::solve_least_squares(dlib::objective_delta_stop_strategy(1e-7),
residual,
residual_derivative,
Pair_CellDC_DistEmi,
Alpha_Omega);
AlphaCoeff_Step2_[poll][pre][xc][yc] = Alpha_Omega(0);
OmegaCoeff_Step2_[poll][pre][xc][yc] = Alpha_Omega(1);
//nPol, nPre, nx, ny
//std::cout << " Alpha:" << AlphaCoeff_Step2_[poll][pre][xc][yc] << " Omega:" << OmegaCoeff_Step2_[poll][pre][xc][yc] << std::endl;
}
}
*/
}
}
}
this->input_info->set_AlphaCoeff_Step2(AlphaCoeff_Step2_);
this->input_info->set_OmegaCoeff_Step2(OmegaCoeff_Step2_);
}
void Training::ComputingPotency_Step1(){
//UInt nTra_ = this->input_info->get_nTra();
UInt nPol_ = this->input_info->get_nPol();
UInt nPre_ = this->input_info->get_nPre();
UInt nX_ = this->input_info->get_nx();
UInt nY_ = this->input_info->get_ny();
int nX__ = nX_;
int nY__ = nY_;
std::vector<int> validScn_Prec = this->validScn_ByPrec;
FourVectorD DETra_ = input_info->get_DETra(); //(nTra, nPre, nx, ny)
FourVectorD DCTra_ = input_info->get_DCTra(); //(nTra, nPol, nx, ny)
SingleVectorD Lat_ = this->input_info->get_Lat_BC();
SingleVectorD Lon_ = this->input_info->get_Lon_BC();
FourVectorD OmegaCoeff_Step2_ = this->input_info->get_OmegaCoeff_Step2();
SingleVectorD dimY(nY_,0);
DoubleVectorD dimYX(nX_,dimY);
double lat_resolution = sqrt(pow(Lat_[1]-Lat_[0],2));
double lon_resolution = sqrt(pow(Lon_[1]-Lon_[0],2));
for (UInt poll=0; poll<nPol_; ++poll){
ThreeVectorD Potency_byPrec;
for (UInt pre=0; pre<nPre_; ++pre){
DoubleVectorD Potency = dimYX;
int scn = validScn_Prec[pre];
if (scn >=0){
double dist_inf = Dist_influence[poll][pre]; // Weighting criteria
UInt x_cells = round(dist_inf/lon_resolution);
UInt y_cells = round(dist_inf/lat_resolution);
//DoubleVectorD DETra_aggregated = dimYX;
for (UInt xc=0; xc<nX_; ++xc){
for (UInt yc=0; yc<nY_; ++yc){
double DETra_aggregated=0;
int start_x = xc-x_cells; if (start_x<0) start_x=0;
int start_y = yc-y_cells; if (start_y<0) start_y=0;
int end_x = xc+x_cells; if (end_x>nX__) end_x=nX_-1;
int end_y = yc+y_cells; if (end_y>nY__) end_y=nY_-1;
for (int xe_t=start_x; xe_t<=end_x; ++xe_t){
for (int ye_t=start_y; ye_t<=end_y; ++ye_t){
double dist_loc = sqrt(pow(Lon_[xc]-Lon_[xe_t],2)+pow(Lat_[yc]-Lat_[ye_t],2));
if (dist_loc<=dist_inf){
DETra_aggregated += DETra_[scn][pre][xe_t][ye_t]*(pow(1+dist_loc,OmegaCoeff_Step2_[poll][pre][xc][yc]));
}
}
}
Potency[xc][yc] = DCTra_[scn][poll][xc][yc]/DETra_aggregated;
//if (pre==2 && DCTra_[scn][poll][xc][yc]<0)
//std::cout << DCTra_[scn][poll][xc][yc] << " " << DETra_aggregated << " " << Potency[xc][yc] << " ";
}
}
}
Potency_byPrec.push_back(Potency);
}
this->Potency.push_back(Potency_byPrec);
}
}
void Training::BellShape_Step1_withPotencyCriteria(){
std::cout << " > BellShape by precursor with Morans' I criteria " << std::endl;
//UInt nTra_ = this->input_info->get_nTra();
UInt nPol_ = this->input_info->get_nPol();
UInt nPre_ = this->input_info->get_nPre();
UInt nX_ = this->input_info->get_nx();
UInt nY_ = this->input_info->get_ny();
int nX__ = nX_;
int nY__ = nY_;
std::vector<int> validScn_Prec = this->validScn_ByPrec;
FourVectorD DETra_ = input_info->get_DETra(); //(nTra, nPre, nx, ny)
FourVectorD DCTra_ = input_info->get_DCTra(); //(nTra, nPol, nx, ny)
SingleVectorD Lat_ = this->input_info->get_Lat_BC();
SingleVectorD Lon_ = this->input_info->get_Lon_BC();
FourVectorD OmegaCoeff_Step2_ = this->input_info->get_OmegaCoeff_Step2();
double lat_resolution = sqrt(pow(Lat_[1]-Lat_[0],2));
double lon_resolution = sqrt(pow(Lon_[1]-Lon_[0],2));
//SingleVectorD dim1(Ny,0);
//DoubleVectorD dim2(Nx,dim1);
SingleVectorD dimY(nY_,0);
DoubleVectorD dimYX(nX_,dimY);
// Computing fitting.
std::cout << " - computing coefficients Alpha and Omega by precursor" << std::endl;
//double Alpha=1., Omega=-2.;
DoubleVectorD AlphaCoeff_Step1_ = this->input_info->get_AlphaCoeff_Step1();
DoubleVectorD OmegaCoeff_Step1_ = this->input_info->get_OmegaCoeff_Step1();
SingleVectorD Clust; Clust.push_back(-1000000); Clust.push_back(-1); Clust.push_back(-0.5); Clust.push_back(-0.25);
Clust.push_back(-0.00001); Clust.push_back(0.00001);Clust.push_back(0.25); Clust.push_back(0.5); Clust.push_back(1);Clust.push_back(1000000);
for (UInt poll=0; poll<nPol_; ++poll){
for (UInt pre=0; pre<nPre_; ++pre){
int scn = validScn_Prec[pre];
if (scn >=0){
double dist_inf = Dist_influence[poll][pre]; // Weighting criteria
UInt x_cells = round(dist_inf/lon_resolution);
UInt y_cells = round(dist_inf/lat_resolution);
std::vector<dlib::sample_pair> Pair_ConcEmi;
//DoubleVectorD DETra_aggregated = dimYX;
for (UInt xc=0; xc<nX_; ++xc){
for (UInt yc=0; yc<nY_; ++yc){
double DETra_aggregated=0;
int start_x = xc-x_cells; if (start_x<0) start_x=0;
int start_y = yc-y_cells; if (start_y<0) start_y=0;
int end_x = xc+x_cells; if (end_x>nX__) end_x=nX_-1;
int end_y = yc+y_cells; if (end_y>nY__) end_y=nY_-1;
for (int xe_t=start_x; xe_t<=end_x; ++xe_t){
for (int ye_t=start_y; ye_t<=end_y; ++ye_t){
double dist_loc = sqrt(pow(Lon_[xc]-Lon_[xe_t],2)+pow(Lat_[yc]-Lat_[ye_t],2));
if (dist_loc<=dist_inf){
DETra_aggregated += DETra_[scn][pre][xe_t][ye_t]*(pow(1+dist_loc,OmegaCoeff_Step2_[poll][pre][xc][yc]));
}
}
}
dlib::sample_pair aa(DCTra_[scn][poll][xc][yc],DETra_aggregated);
//Potency[xc][yc] = DCTra_[scn][poll][xc][yc]/DETra_aggregated;
Pair_ConcEmi.push_back(aa);
}
}
std::vector<unsigned long> labels;
int N_Clusters = dlib::newman_cluster(Pair_ConcEmi,labels,1e-4,2000);
std::cout << N_Clusters << " " << labels.size() << std::endl;
}
}
}
/*
for (UInt poll=0; poll<nPol_; ++poll){
for (UInt pre=0; pre<nPrec_; ++pre){
DoubleVectorD Potency_ = this->Potency[poll][pre];
std::vector<std::pair<double,double>> Pair_ConcEmi;
dlib::sample_pair aa(1.,2.);
const std::vector<dlib::sample_pair> edges;
std::vector<unsigned long> labels;
dlib::newman_cluster(edges,labels,1e-4,2000);
double max = std::numeric_limits<double>::lowest();
double min = std::numeric_limits<double>::max();
for (const auto& v : Potency_){
double current_max = *std::max_element(v.cbegin(), v.cend());
double current_min = *std::min_element(v.cbegin(), v.cend());
max = max < current_max ? current_max : max; // max = std::max(current_max, max);
min = min > current_min ? current_min : min; // max = std::max(current_max, max);
}
std::cout << min << " " << max << std::endl;
}
}
*/
}
void Training::BellShape_Step2_TraditionalRegression(int SelectWindows){
std::cout << " > BellShape by precursor and Receptor cell" << std::endl;
UInt nTra_ = this->input_info->get_nTra();
UInt nPol_ = this->input_info->get_nPol();
UInt nPre_ = this->input_info->get_nPre();
UInt nX_ = this->input_info->get_nx();
UInt nY_ = this->input_info->get_ny();
FourVectorD DCTra_ = input_info->get_DCTra(); //(nTra, nPol, nx, ny)
FourVectorD DETra_ = input_info->get_DETra(); //(nTra, nPre, nx, ny)
SingleVectorD Lon_ = this->input_info->get_Lon_BC();
SingleVectorD Lat_ = this->input_info->get_Lat_BC();
//DoubleVectorD AlphaCoeff_Step1_ = this->input_info->get_AlphaCoeff_Step1(); //(nPol, nPre)
//DoubleVectorD OmegaCoeff_Step1_ = this->input_info->get_OmegaCoeff_Step1(); //(nPol, nPre)
std::vector<int> validScn_Prec = this->validScn_ByPrec;
SingleVectorD dim1(nY_,0);
DoubleVectorD dim2(nX_,dim1);
ThreeVectorD dim3(nPre_,dim2);
FourVectorD AlphaCoeff_Step2_(nPol_,dim3); //(nPol, nPre, nx, ny)
//FourVectorD OmegaCoeff_Step2_(nPol_,dim3); //(nPol, nPre, nx, ny)
FourVectorD OmegaCoeff_Step2_ = this->input_info->get_OmegaCoeff_Step2();
for (UInt poll=0; poll<nPol_; ++poll){
// Compute the aggregations of emissions with distances & omega values
FourVectorD DETra_AggVRAI(nTra_,dim3); //(nTra, nPre, nx, ny) initialized with 0 values
for (UInt scn=0; scn<nTra_; ++scn){
for (UInt pre=0; pre<nPre_; ++pre){
for (UInt xc=0; xc<nX_; ++xc){
for (UInt yc=0; yc<nY_; ++yc){
for (UInt xe=0; xe<nX_; ++xe){
for (UInt ye=0; ye<nY_; ++ye){
double dist_c_e = sqrt(pow(Lon_[xe]-Lon_[xc],2) + pow(Lat_[ye]-Lat_[yc],2));
DETra_AggVRAI[scn][pre][xc][yc] += DETra_[scn][pre][xe][ye]* pow((1+dist_c_e),OmegaCoeff_Step2_[poll][pre][xc][yc]);
}
}
/*for (int xe=xc-10; xe<xc+10; ++xe){
for (int ye=yc-10; ye<yc+10; ++ye){
if (xe>=0 && xe<nX_ && ye>=0 && ye<nY_){
double dist_c_e = sqrt(pow(Lon_[xe]-Lon_[xc],2) + pow(Lat_[ye]-Lat_[yc],2));
DETra_AggVRAI[scn][pre][xc][yc] += DETra_[scn][pre][xe][ye]* pow((1+dist_c_e),OmegaCoeff_Step2_[poll][pre][xc][yc]);
}
}
}*/
}
}
}
}
// Computing Alpha by cell
int validScn_Prec_size = validScn_Prec.size();
int nTra_int = nTra_;
for (UInt xc=0; xc<nX_; ++xc){
for (UInt yc=0; yc<nY_; ++yc){
SingleVectorD All_poll_cell;
DoubleVectorD All_prec_cell;
for (int xw=-SelectWindows; xw<=SelectWindows; ++xw){
for (int yw=-SelectWindows; yw<=SelectWindows; ++yw){
if ((xc+xw)>=0 && (yc+yw)>=0 && (xc+xw)<nX_ && (yc+yw)<nY_){
for (int scn=0; scn<nTra_int; ++scn){
UInt ok_=1; // this avoid the use of scenarios used on step1
for (int scnNO=0; scnNO<validScn_Prec_size; ++scnNO){
if (scn == validScn_Prec[scnNO]) ok_=0;
}
if (ok_==1){
double poll_cell = DCTra_[scn][poll][xc+xw][yc+yw];
All_poll_cell.push_back(poll_cell);
SingleVectorD precs_ByPollCell;
for (UInt pre=0; pre<nPre_; ++pre){
double prec_i = DETra_AggVRAI[scn][pre][xc+xw][yc+yw];
precs_ByPollCell.push_back(prec_i);
}
All_prec_cell.push_back(precs_ByPollCell);
}
}
}
}
}
dlib::matrix<double,0,1> B_concCells;
B_concCells.set_size(All_poll_cell.size());
//B_concCells =0;
dlib::matrix<double> A_precursors_agg;
A_precursors_agg.set_size(All_prec_cell.size(),All_prec_cell[0].size());
//A_precursors_agg =0;
for (UInt r=0; r<All_poll_cell.size(); ++r){
B_concCells(r) = All_poll_cell[r];
for (UInt c=0; c<All_prec_cell[r].size(); ++c){
A_precursors_agg(r,c) = All_prec_cell[r][c];
}
}
dlib::qr_decomposition < dlib::matrix<double> > qr(A_precursors_agg);
//dlib::matrix<double> R_precursors_agg = qr.get_r();
//dlib::matrix<double> Q_precursors_agg = qr.get_q();
dlib::matrix<double> Alpha_coeff = qr.solve(B_concCells);
/*
if (qr.is_full_rank()) std::cout <<""; //std::cout << " ok: input matrix of precursors is full of rank " << std::endl;
else {
std::cout << " error: input matrix is not full of rank, that is mean the first columns are not linear independent " << std::endl;
std::cout << "xc:" << xc << " yc:" << yc << std::endl;
std::cout << "nEqu: " << B_concCells.nr() << std::endl;
std::cout << "xc:" << Lon_[xc] << " yc:" << Lat_[yc] << std::endl << "coeff: \n" << trans(Alpha_coeff) << std::endl;
//std::cout << B_concCells << std::endl;
//std::cout << A_precursors_agg << std::endl;
}
*/
//std::cout << "xc:" << Lon_[xc] << " yc:" << Lat_[yc] << std::endl << "coeff: \n" << trans(Alpha_coeff) << std::endl;
for (UInt pre=0; pre<nPre_; ++pre){
AlphaCoeff_Step2_[poll][pre][xc][yc] = Alpha_coeff(pre);
//OmegaCoeff_Step2_[poll][pre][xc][yc] = OmegaCoeff_Step1_[poll][pre];
}
//Compute Determination Coefficient R
//dlib::matrix <double> Rxy = dlib::xcorr(B_concCells,A_precursors_agg);
//std::cout << Rxy << std::endl;
}
}
}
this->input_info->set_AlphaCoeff_Step2(AlphaCoeff_Step2_);
//this->input_info->set_OmegaCoeff_Step2(OmegaCoeff_Step2_);
}
bool myfunction (int i,int j) {return (i>j);}
void Training::BellShape_Step2_PCRegression(int SelectWindows){
std::cout << " > BellShape by precursor and Receptor cell" << std::endl;
UInt nTra_ = this->input_info->get_nTra();
UInt nPol_ = this->input_info->get_nPol();
UInt nPre_ = this->input_info->get_nPre();
UInt nX_ = this->input_info->get_nx();
UInt nY_ = this->input_info->get_ny();
SingleVectorD Dist_Agg_ = this->input_info->get_DistAgg(); // (nAgg)
FiveVectorD DETraAgg_ = this->input_info->get_DEtraAgg(); //(nTra, nPre, nx, ny, nAgg)
FourVectorD DCTra_ = input_info->get_DCTra(); //(nTra, nPol, nx, ny)
FourVectorD DETra_ = input_info->get_DETra(); //(nTra, nPre, nx, ny)
SingleVectorD Lon_ = this->input_info->get_Lon_BC();
SingleVectorD Lat_ = this->input_info->get_Lat_BC();
DoubleVectorD AlphaCoeff_Step1_ = this->input_info->get_AlphaCoeff_Step1(); //(nPol, nPre)
DoubleVectorD OmegaCoeff_Step1_ = this->input_info->get_OmegaCoeff_Step1(); //(nPol, nPre)
std::vector<int> validScn_Prec = this->validScn_ByPrec;
SingleVectorD dim1(nY_,0);
DoubleVectorD dim2(nX_,dim1);
ThreeVectorD dim3(nPre_,dim2);
FourVectorD AlphaCoeff_Step2_(nPol_,dim3); //(nPol, nPre, nx, ny)
//FourVectorD OmegaCoeff_Step2_(nPol_,dim3); //(nPol, nPre, nx, ny)
FourVectorD OmegaCoeff_Step2_ = this->input_info->get_OmegaCoeff_Step2();
int validScn_Prec_size = validScn_Prec.size();
int nTra_int = nTra_;
for (UInt poll=0; poll<nPol_; ++poll){
/*
// Complete the aggregations with distances & omega values
FourVectorD DETra_AggVRAI(nTra_,dim3); //(nTra, nPre, nx, ny) initialized with 0 values
for (UInt scn=0; scn<nTra_; ++scn){
for (UInt pre=0; pre<nPre_; ++pre){
for (UInt xe=0; xe<nX_; ++xe){
for (UInt ye=0; ye<nY_; ++ye){
for (UInt agg=0; agg<DETraAgg_[scn][pre][xe][ye].size(); ++agg){
//DETra_AggVRAI[scn][pre][xe][ye] += DETraAgg_[scn][pre][xe][ye][agg]* pow((1+Dist_Agg_[agg]),OmegaCoeff_Step1_[poll][pre]);
DETra_AggVRAI[scn][pre][xe][ye] += DETraAgg_[scn][pre][xe][ye][agg]* pow((1+Dist_Agg_[agg]),OmegaCoeff_Step2_[poll][pre][xe][ye]);
}
}
}
}
}
*/
// Compute the aggregations of emissions with distances & omega values
FourVectorD DETra_AggVRAI(nTra_,dim3); //(nTra, nPre, nx, ny) initialized with 0 values
for (UInt scn=0; scn<nTra_; ++scn){
for (UInt pre=0; pre<nPre_; ++pre){
for (UInt xc=0; xc<nX_; ++xc){
for (UInt yc=0; yc<nY_; ++yc){
for (UInt xe=0; xe<nX_; ++xe){
for (UInt ye=0; ye<nY_; ++ye){
double dist_c_e = sqrt(pow(Lon_[xe]-Lon_[xc],2) + pow(Lat_[ye]-Lat_[yc],2));
DETra_AggVRAI[scn][pre][xc][yc] += DETra_[scn][pre][xe][ye]* pow((1+dist_c_e),OmegaCoeff_Step2_[poll][pre][xc][yc]);
}
}
/*for (int xe=xc-10; xe<xc+10; ++xe){
for (int ye=yc-10; ye<yc+10; ++ye){
if (xe>=0 && xe<nX_ && ye>=0 && ye<nY_){
double dist_c_e = sqrt(pow(Lon_[xe]-Lon_[xc],2) + pow(Lat_[ye]-Lat_[yc],2));
DETra_AggVRAI[scn][pre][xc][yc] += DETra_[scn][pre][xe][ye]* pow((1+dist_c_e),OmegaCoeff_Step2_[poll][pre][xc][yc]);
}
}
}*/
}
}
}
}
// Computing Alpha by cell
for (UInt xc=0; xc<nX_; ++xc){
for (UInt yc=0; yc<nY_; ++yc){
SingleVectorD All_poll_cell;
DoubleVectorD All_prec_cell;
for (int xw=-SelectWindows; xw<=SelectWindows; ++xw){
for (int yw=-SelectWindows; yw<=SelectWindows; ++yw){
if ((xc+xw)>=0 && (yc+yw)>=0 && (xc+xw)<nX_ && (yc+yw)<nY_){
for (int scn=0; scn<nTra_int; ++scn){
UInt ok_=1; // this avoid the use of scenarios used on step1
for (int scnNO=0; scnNO<validScn_Prec_size; ++scnNO){
if (scn == validScn_Prec[scnNO]) ok_=0;
}
if (ok_==1){
double poll_cell = DCTra_[scn][poll][xc+xw][yc+yw];
All_poll_cell.push_back(poll_cell);
SingleVectorD precs_ByPollCell;
for (UInt pre=0; pre<nPre_; ++pre){
double prec_i = DETra_AggVRAI[scn][pre][xc+xw][yc+yw];
precs_ByPollCell.push_back(prec_i);
}
All_prec_cell.push_back(precs_ByPollCell);
}
}
}
}
}
dlib::matrix<double,0,1> B_concCells;
B_concCells.set_size(All_poll_cell.size());
//B_concCells =0;
dlib::matrix<double> A_precursors_agg;
A_precursors_agg.set_size(All_prec_cell.size(),All_prec_cell[0].size());
//A_precursors_agg =0;
for (UInt r=0; r<All_poll_cell.size(); ++r){
B_concCells(r) = All_poll_cell[r];
for (UInt c=0; c<All_prec_cell[r].size(); ++c){
A_precursors_agg(r,c) = All_prec_cell[r][c];
}
}
dlib::matrix <double> U, D, V; // D is W
dlib::svd(A_precursors_agg,U,D,V);
//Ordering vectors according to the values of D
SingleVectorD DiagValues;
for (UInt r=0; r<D.nr(); ++r){
for (UInt c=0; c<D.nc(); ++c){
if (r==c) DiagValues.push_back(D(r,c));
}
}
SingleVectorD DiagValues_ordered = DiagValues;
SingleVectorD position_v(DiagValues_ordered.size(),0);
std::sort(DiagValues_ordered.begin(),DiagValues_ordered.end(),myfunction);
for (UInt v=0; v<DiagValues_ordered.size(); ++v){
for (UInt v_=0; v_<DiagValues.size(); ++v_){
if (DiagValues_ordered[v] == DiagValues[v_])
position_v[v] = v_;
}
}
UInt remained_dimensions;
if (A_precursors_agg.nc() > 2) remained_dimensions = A_precursors_agg.nc() * 0.7;
else remained_dimensions = A_precursors_agg.nc();
dlib::matrix<double> V_reduced;
V_reduced.set_size(V.nr(), remained_dimensions);
for (UInt r=0; r<V_reduced.nr(); ++r){
for (UInt c=0; c<V_reduced.nc(); c++){
V_reduced(r,c) = V(r,position_v[c]);
}
}
dlib::matrix<double> W = A_precursors_agg*V_reduced;
dlib::matrix<double> Ro = inv(trans(W)*W) * trans(W)*B_concCells;
dlib::matrix<double> Alpha_coeff = V_reduced * Ro;
//std::cout << trans(Alpha_coeff) << std::endl;
//std::cout << Alpha_coeff.size() << " " << W.nr() << " " << W.nc() << " " << B_concCells.nr() << std::endl;
for (UInt pre=0; pre<nPre_; ++pre){
AlphaCoeff_Step2_[poll][pre][xc][yc] = Alpha_coeff(pre);
//OmegaCoeff_Step2_[poll][pre][xc][yc] = OmegaCoeff_Step1_[poll][pre];
}
}
}
}
this->input_info->set_AlphaCoeff_Step2(AlphaCoeff_Step2_);
//this->input_info->set_OmegaCoeff_Step2(OmegaCoeff_Step2_);
}
c4science · Help