27 #include "IpIpoptApplication.hpp" 28 #include "IpSolveStatistics.hpp" 102 for(uint i=0;i<
regIds2.size();i++){
155 double stvalue = sl+sa;
156 double q1 = 1/ ( 1+pow(dl/da,1/psi)*(pl/pWo) );
159 q1*pow(da,(psi-1)/psi) + p1*pow(dl,(psi-1)/psi)
163 double pc = (da/dc)*pWo
165 double pw = (dl*pl+da*pWo)/(dl+da);
191 double t1 = 1/ ( 1+(pow(sl/sa,1/eta))*(pl/pWo) );
194 t1*pow(sa,(eta-1)/eta) + r1*pow(sl,(eta-1)/eta)
198 double pc = (sa/sc)*pWo+(sl/sc)*pl;
199 double pw = (sl*pl+sa*pWo)/(sl+sa);
214 for(uint r1=0;r1<
l2r.size();r1++){
215 for(uint r2=0;r2<
l2r[r1].size();r2++){
217 for(uint r2To=0;r2To<
l2r[r1].size();r2To++){
219 if(
l2r[r1][r2] ==
l2r[r1][r2To]){
225 for(uint r2To=0;r2To<
l2r[r1].size();r2To++){
226 if(
l2r[r1][r2] ==
l2r[r1][r2To]){
237 static double cumOverHarvesting = 0.0;
248 for(uint i=0;i<
regIds2.size();i++){
256 double pol_mktDirInt_d =
gpd(
"pol_mktDirInt_d",r2,
secProducts[sp]);
257 double pol_mktDirInt_d_1 =
gpd(
"pol_mktDirInt_d",r2,
secProducts[sp],previousYear);
261 double q1_polCorr = min(1.0, (q1-0.5)*pol_mktStr_d+0.5);
265 double k = (1+g1)*k_1;
266 double aa = (sigma/(sigma+1))*pc_1*pow(dc_1,-1/sigma) * pol_mktDirInt_d_1/pol_mktDirInt_d ;
267 double gg = dc_1*pow(pc_1*pol_mktDirInt_d_1,-sigma);
300 double pol_mktDirInt_s =
gpd(
"pol_mktDirInt_s",r2,
priProducts[pp]);
301 double pol_mktDirInt_s_1 =
gpd(
"pol_mktDirInt_s",r2,
priProducts[pp],previousYear);
307 double carbonPrice_1 =
gpd(
"carbonPrice",r2,
"",previousYear);
308 double omega =
gpd(
"co2content_products", r2,
priProducts[pp])/1000;
309 double Pct_1 = carbonPrice_1*intRate;
314 double t1_polCorr = min(1.0, (t1-0.5)*pol_mktStr_s+0.5);
317 double gamma = in>in_1 ? gamma_incr: gamma_decr;
318 double gamma_d = in_d>in_d_1 ? gamma_d_incr: gamma_d_decr;
321 double in_ratio = in/in_1;
322 double in_d_ratio = useDeathTimber ? max(min_in_d,in_d) / max(min_in_d, in_d_1 ) : 1.0 ;
331 if (in_d_ratio < 0.01){
334 if (in_d_ratio > 100){
339 double bb = (sigma/(sigma+1.0))*pc_1*pow(sc_1,-1.0/sigma)*pow(1/in_ratio,gamma/sigma)*pow(1.0,1.0/sigma) * pol_mktDirInt_s_1/pol_mktDirInt_s ;
346 double ff = sc_1*pow((pc_1 - omega * Pct_1)*pol_mktDirInt_s_1,-sigma)*pow(in_ratio,gamma)*pow(in_d_ratio,gamma_d);
365 SmartPtr<IpoptApplication> application =
new IpoptApplication();
367 application->Options()->SetStringValue(
"linear_solver", linearSolver);
370 application->Options()->SetNumericValue(
"obj_scaling_factor", -1);
371 application->Options()->SetNumericValue(
"max_cpu_time", 1800);
372 application->Options()->SetStringValue(
"check_derivatives_for_naninf",
"yes");
386 ApplicationReturnStatus status;
387 status = application->Initialize();
388 if (status != Solve_Succeeded) {
389 printf(
"\n\n*** Error during initialization!\n");
390 msgOut(
MSG_INFO,
"Error during initialization! Do you have the solver compiled for the specified linear solver?");
394 msgOut(
MSG_INFO,
"Running optimisation problem for this year (it may take a few minutes for large models)..");
395 status = application->OptimizeTNLP(
MTHREAD->
OPT);
407 if (status == Solve_Succeeded) {
408 Index iter_count = application->Statistics()->IterationCount();
409 Number final_obj = application->Statistics()->FinalObjective();
410 printf(
"\n*** The problem solved year %d in %d iterations!\n", thisYear, iter_count);
411 printf(
"\n*** The final value of the objective function is %e.\n", final_obj);
412 msgOut(
MSG_INFO,
"The problem solved successfully year "+
i2s(thisYear)+
" in "+
i2s(iter_count)+
" iterations.");
413 int icount = iter_count;
414 double obj = final_obj;
419 cout <<
"***ERROR: MODEL DIDN'T SOLVE FOR THIS YEAR ("<<thisYear<<
")"<< endl;
429 for(uint r2= 0; r2<
regIds2.size();r2++){
466 double pw = dt?(dl*pl+da*pworld)/dt:0.0;
475 double st_fromFor = st *
gpd(
"supply_fromForestShare",regId,
priProducts[p]);
476 double pw = st?(sl*pl+sa*pworld)/st:0.0;
485 int nPriPrCombs = priPrCombs.size();
487 for (uint i=0;i<priPrCombs.size();i++){
488 double stMkMod = 0.0;
491 for (uint p=0;p<priPrCombs[i].size();p++){
502 string pProductsInvolved =
"";
503 for (uint p=0;p<priPrCombs[i].size();p++){
504 pProductsInvolved += (
priProducts[priPrCombs[i][p]]+
"; ");
506 double inV_over_hV_ratio = stMkMod ? sumIn/stMkMod : 0.0;
507 cumOverHarvesting += (stMkMod-sumIn);
508 msgOut(
MSG_DEBUG,
"Overharvesting has happened. Year: "+
i2s(thisYear)+
"Region: "+
i2s(regId)+
"Involved products: "+pProductsInvolved+
". sumIn: "+
d2s(sumIn)+
" stMkMod:" +
d2s(stMkMod) +
" cumOverHarvesting: "+
d2s(cumOverHarvesting));
509 for (uint p=0;p<priPrCombs[i].size();p++){
510 double st_fromFor_orig =
gpd(
"st_fromFor",regId,
priProducts[priPrCombs[i][p]]);
511 spd(st_fromFor_orig*inV_over_hV_ratio,
"st_fromFor",regId,
priProducts[priPrCombs[i][p]]);
523 vector <double> in_deathTimber(
priProducts.size(),0.);
524 vector <double> in_aliveForest (
priProducts.size(),0.);
541 if (cumOverHarvesting>0.0){
542 msgOut(
MSG_DEBUG,
"Overharvesting is present. Year: "+
i2s(thisYear)+
" cumOverHarvesting: "+
d2s(cumOverHarvesting));
568 for(uint i=0;i<
regIds2.size();i++){
574 double shareMortalityUsableTimber = 0.0;
576 for (uint p=0;p<
regPx.size();p++){
579 double pxId = px->
getID();
584 for(uint j=0;j<
fTypes.size();j++){
587 vector <double> hV_byDiam;
588 vector <double> productivity_byDc;
589 vector < vector <double> > hV_byDiamAndPrd;
590 vector <double> hArea_byDc;
591 vector <double> newVol_byDiam;
592 vector <double> vMort_byDc;
593 vector <double> vMortAdd_byDc ;
594 vector <double> areasMovingUp(
dClasses.size(), 0.0);
595 double areaFirstProdClass;
609 int tp_u0 = px->
tp.at(j).at(0);
615 double availableArea = px->
area_l.at(j).at(0);
617 double vHa = px->
vHa.at(j).at(1);
620 areaFirstProdClass = pastRegArea;
622 areaFirstProdClass = min(availableArea, pastRegArea);
624 px->
vReg.push_back(areaFirstProdClass*vHa/1000000.0);
629 if (areaFirstProdClass < 0.0){
632 if ( (availableArea-pastRegArea) < -0.00000001 ){
641 double regRegVolumes =
gfd(
"vReg",r2,ft,
"");
642 double newVReg = regionArea ? regRegVolumes*pxArea/regionArea : 0.0;
643 px->
vReg.push_back(newVReg);
646 areaFirstProdClass = (1.0 / ((double) tp_u0) ) * px->
initialDc0Area.at(j);
656 if (areaFirstProdClass<0.0){
659 if (areaFirstProdClass > px->
area_l.at(j).at(0)){
672 shareMortalityUsableTimber =
gfd(
"shareMortalityUsableTimber",r2,ft,
"");
674 shareMortalityUsableTimber = 0.0;
677 for (uint u=0; u<
dClasses.size(); u++){
683 double pastYearVol = px->
vol_l.at(j).at(u);
684 vector <double> hV_byPrd;
685 vector <double> hr_byPrd;
696 hr_byPrd.push_back( hr_pr);
705 double hr_pr = origHr ? hr_byPrd[pp] * min(1.0,1.0/origHr) : 0.0;
706 hV_byPrd.push_back( hr_pr*pastYearVol*px->
avalCoef.at(j));
709 double hV = hr*pastYearVol*px->
avalCoef.at(j);
712 hV_byDiam.push_back(hV);
713 hV_byDiamAndPrd.push_back(hV_byPrd);
727 double tp = px->
tp.at(j).at(u);
728 double mort = px->
mort.at(j).at(u);
729 double addMort = px->
addMort.at(j).at(u);
730 double vReg = px->
vReg.at(j);
731 double beta = px->
beta.at(j).at(u);
733 double vHa = px->
vHa.at(j).at(u);
734 double finalHarvestFlag =
gfd(
"finalHarvestFlag",regId,ft,dc);
736 double vMort = mort*pastYearVol;
737 double vMortAdd = addMort*pastYearVol;
739 vMort_byDc.push_back(vMort);
740 vMortAdd_byDc.push_back(vMortAdd);
743 iisskey key(thisYear,r2,ft,dc);
751 vol = max(0.0,(1-1/tp-mort-addMort))*pastYearVol+vReg;
754 if ((1-1/tp-mort-addMort)<0.0){
755 msgOut(
MSG_DEBUG,
"The sum of leaving trres and mortality would have lead to nevative volume if we didn't put a max. 1/tp: "+
d2s(1/tp)+
", mort: "+
d2s(mort)+
", total coeff: "+
d2s((1-1/tp-mort))+
" ");
760 double inc = (u==
dClasses.size()-1)?0:1./tp;
761 double tp_1 = px->
tp.at(j).at(u-1);
762 double pastYearVol_1 = px->
vol_l.at(j).at(u-1);
764 vol = max(0.0,(1-inc-mort-addMort)*pastYearVol-hV+(1/tp_1)*beta*pastYearVol_1);
766 if ((1-inc-mort)*pastYearVol-hV+(1/tp_1)*beta*pastYearVol_1 < 0){
767 double realVolumes = (1-inc-mort-addMort)*pastYearVol-hV+(1/tp_1)*beta*pastYearVol_1;
768 msgOut(
MSG_DEBUG,
"Negative real volumes ("+
d2s(realVolumes)+
"), possibly because of little bit larger bounds in the market module to avoid zeros. Volumes in the resource module set back to zero, so it should be ok.");
773 double inc = (u==
dClasses.size()-1)?0:1.0/tp;
774 double volumesMovingUp = inc*pastYearVol;
775 double pastArea = px->
area_l.at(j).at(u);
777 areasMovingUp.at(u) = inc*pastArea;
780 hArea_byDc.push_back(finalHarvestFlag*1000000*hV/vHa);
782 double finalHarvestedVolumes = finalHarvestFlag* hV;
783 double finalHarvestedRate = pastYearVol?finalHarvestedVolumes/pastYearVol:0.0;
785 if (finalHarvestedRate > 1.0){
789 hArea_byDc.push_back(finalHarvestedRate*pastArea);
791 px->
area.at(j).at(u) = max(0.0, px->
area_l.at(j).at(u) - areasMovingUp.at(u) + areasMovingUp.at(u-1) - hArea_byDc.at(u));
793 if ((px->
area_l.at(j).at(u) - areasMovingUp.at(u) + areasMovingUp.at(u-1) - hArea_byDc.at(u))< 0.0){
794 msgOut(
MSG_DEBUG,
"If not for a max, we would have had a negative area ("+
d2s(px->
area_l.at(j).at(u) - areasMovingUp.at(u) + areasMovingUp.at(u-1) - hArea_byDc.at(u))+
" ha).");
798 areasMovingUp.at(u) = areaFirstProdClass;
799 hArea_byDc.push_back(0.);
800 px->
area.at(j).at(u) = px->
area_l.at(j).at(u) - areasMovingUp.at(u) - hArea_byDc.at(u);
805 newVol_byDiam.push_back(vol);
807 if(px->
area.at(j).at(u)< 0.0 || areasMovingUp.at(u) < 0.0 || hArea_byDc.at(u) < 0.0 ){
820 double productivity = hArea_byDc.at(u) ? (1000000.0 * hV_byDiam.at(u))/(hArea_byDc.at(u)*age) : 0.0;
821 productivity_byDc.push_back(productivity);
823 px->
hVol.push_back(hV_byDiam);
825 px->
hArea.push_back(hArea_byDc);
826 px->
vol.push_back(newVol_byDiam);
827 px->
vMort.push_back(vMort_byDc);
828 px->
vMortAdd.push_back(vMortAdd_byDc);
833 for (uint u=1; u<
dClasses.size(); u++){
834 double volMort = vMort_byDc[u];
835 double harvVol = hV_byDiam[u];
836 double vol_new = newVol_byDiam[u];
837 double vol_lagged = px->
vol_l.at(j).at(u);
838 double gain = vol_new - (vol_lagged-harvVol-volMort);
839 if (volMort > vol_lagged){
854 for (uint p=0;p<
regPx.size();p++){
855 for(uint j=0;j<
fTypes.size();j++){
856 for (uint u=0; u<
dClasses.size(); u++){
859 sumHv +=
regPx[p]->hVol_byPrd[j][u][pp];
864 if(abs(sumSt-sumHv) > 0.000001){
878 map<string,double> hAreaByFTypeGroup =
vectorToMap(allFTypes,0.0);
883 for(uint i=0;i<
regIds2.size();i++){
898 double carbonPrice_NOW =
gpd(
"carbonPrice",r2,
"",thisYear+1);
899 vector<double> futureCarbonPrices;
900 if(flagHartman!=
"NONE"){
901 for (
int y=0; y<1000; y++) {
902 futureCarbonPrices.push_back(
gpd(
"carbonPrice",r2,
"",thisYear + y));
906 vector<double> polBal_fiSub (
fTypes.size(),0.0);
909 double fArea_reg = REG->
getArea();
910 double fArea_diff = 0.0;
911 double fArea_reldiff = 0.0;
913 fArea_reldiff =
gfd(
"forestChangeAreaIncrementsRel",r2,
"",
"",
DATA_NOW);
914 fArea_diff = fArea_reg * fArea_reldiff;
916 fArea_diff =
gfd(
"forestChangeAreaIncrementsHa",r2,
"",
"",
DATA_NOW);
917 fArea_reldiff = fArea_diff / fArea_reg;
920 double regHArea = 0.0;
945 vector<vector<double>> futureCarbonRef;
946 vector<vector<double>> futureExtraBiomass_ratio;
948 if(flagHartman!=
"NONE"){
949 for(uint j=0;j<
fTypes.size();j++){
951 vector<double> futureCarbonRef_ft;
952 vector<double>futureExtraBiomass_ratio_ft;
953 for (
int y=0; y<1000; y++) {
954 double carbonRef =
gfd(
"carbonRef",r2,ft,
"", thisYear + y);
955 futureCarbonRef_ft.push_back(carbonRef);
956 double extraBiomass_ratio =
gfd(
"extraBiomass_ratio",regId,ft,
"",
DATA_NOW);
957 futureExtraBiomass_ratio_ft.push_back(extraBiomass_ratio);
959 futureCarbonRef.push_back(futureCarbonRef_ft);
960 futureExtraBiomass_ratio.push_back(futureExtraBiomass_ratio_ft);
965 for (uint p=0;p<
regPx.size();p++){
976 double totalHarvestedAreaOrig =
vSum(px->
hArea);
977 double totalHarvestedArea = 0.0;
978 vector<double> thisYearRegAreas(
fTypes.size(),0.0);
979 vector<double> expectedReturns(
fTypes.size(),0.0);
980 vector<double> expectedReturnsCarbon(
fTypes.size(),0.0);
981 vector<int> rotation_dc(
fTypes.size(),0.0);
982 vector< vector<double> > deltaAreas(
fTypes.size()+1, vector<double>(
fTypes.size()+1,0));
987 double fArea_diff_px = fArea_px * fArea_diff/ fArea_reg;
989 double fArea_incr = max(0.0,fArea_diff_px);
990 double fArea_incr_rel = totalHarvestedAreaOrig?fArea_incr/totalHarvestedAreaOrig:0.0;
991 double fArea_decr = - min(0.0,fArea_diff_px);
992 double fArea_decr_rel = totalHarvestedAreaOrig?min(1.0,fArea_decr/totalHarvestedAreaOrig):0.0;
993 regHArea += totalHarvestedAreaOrig;
994 totalHarvestedArea = totalHarvestedAreaOrig *(1-fArea_decr_rel);
998 for(uint j=0;j<
fTypes.size();j++){
1001 double hAreaThisFt=
vSum(px->
hArea.at(j))*(1-fArea_decr_rel);
1020 for(uint j=0;j<
fTypes.size();j++){
1024 double co2content_inventory =
gfd(
"co2content_inventory",r2,ft,
"",
DATA_NOW)/1000;
1028 double expReturns = std::numeric_limits<double>::lowest();
1029 double expReturns_carbon = std::numeric_limits<double>::lowest();
1030 double expReturns_timber = std::numeric_limits<double>::lowest();
1032 for (uint u=0; u<
dClasses.size(); u++){
1034 double vHa = px->
vHa_exp.at(j).at(u);
1035 double finalHarvestFlag =
gfd(
"finalHarvestFlag",regId,ft,dc);
1036 double cumTp_u = px->
cumTp_exp.at(j).at(u);
1038 double pol_mktDirInt_s =
gpd(
"pol_mktDirInt_s",regId,
priProducts[pp]);
1039 double pol_mktDirInt_s_fut =
gpd(
"pol_mktDirInt_s",regId,
priProducts[pp],thisYear+cumTp_u);
1041 double worldCurPrice =
gpd(
"pl",
WL2,
priProducts[pp])*pol_mktDirInt_s*pol_mktDirInt_s;
1042 double worldFutPrice =
gpd(
"pl",
WL2,
priProducts[pp],thisYear+cumTp_u)*pol_mktDirInt_s_fut;
1046 double raw_amount = finalHarvestFlag*pw_exp*vHa*
app(
priProducts[pp],ft,dc);
1048 double anualised_amount = anualised_amount_timber;
1052 double raw_amount_carbon = 0;
1053 if (flagHartman ==
"EXO") {
1054 double carbonPrice_y = carbonPrice_NOW;
1055 for (
int y=0; y<cumTp_u; y++) {
1057 double carbonRef = futureCarbonRef[j][y];
1058 double expectedVHa_y =
getVHaByYear(px, j, y, futureExtraBiomass_ratio[j][y], regId);
1059 double carbonPayment_y = carbonPrice_y*ir*(co2content_inventory*expectedVHa_y - carbonRef);
1060 raw_amount_carbon += carbonPayment_y*pow(1+ir,cumTp_u - y);
1064 raw_amount_carbon += pickling * carbonPrice_y * finalHarvestFlag * vHa *
app(
priProducts[pp],ft,dc);
1068 anualised_amount += anualised_amount_carbon;
1074 if (anualised_amount>expReturns) {
1075 expReturns=anualised_amount;
1077 expReturns_carbon = anualised_amount_carbon;
1078 expReturns_timber = anualised_amount_timber;
1084 if (flagHartman==
"NONE" || flagHartman==
"EXO") {
1088 px->
optDc.push_back(optDc);
1096 double origExpReturns = expReturns;
1097 expReturns = origExpReturns * (1.0 - ra*cumMort);
1101 if (flagHartman==
"NONE" || flagHartman==
"EXO") {
1105 expectedReturns.at(j) = expReturns;
1106 rotation_dc.at(j) = optDc;
1117 if (flagHartman ==
"ENDO") {
1119 int ft_opt_index =
getMaxPos(expectedReturns);
1120 string ft_opt =
fTypes[ft_opt_index];
1121 int dc_opt_index = rotation_dc[ft_opt_index];
1122 string dc_opt =
dClasses[dc_opt_index];
1123 double cumTP_opt = px->
cumTp_exp.at(ft_opt_index).at(dc_opt_index);
1124 double co2content_inventory_opt =
gfd(
"co2content_inventory",r2,ft_opt,
"",
DATA_NOW)/1000;
1127 vector <double> FutureCarbonContent_opt;
1128 for (
int y=0; y<=cumTP_opt; y++) {
1129 double expectedVHa_opt_y =
getVHaByYear(px, ft_opt_index, y, futureExtraBiomass_ratio[ft_opt_index][y], regId);
1130 FutureCarbonContent_opt.push_back(co2content_inventory_opt*expectedVHa_opt_y);
1133 while (FutureCarbonContent_opt.size()<=2000){
1134 FutureCarbonContent_opt.insert(std::end(FutureCarbonContent_opt), std::begin(FutureCarbonContent_opt), std::end(FutureCarbonContent_opt));
1140 for(uint j=0;j<
fTypes.size();j++){
1142 double co2content_inventory =
gfd(
"co2content_inventory",r2,ft,
"",
DATA_NOW)/1000;
1144 double expReturns_hartman = -10000000000;
1145 double expReturns_timber_hartman;
1146 double expReturns_carbon_hartman;
1147 int optDc_hartman = -1;
1148 for (uint u=1; u<
dClasses.size(); u++){
1150 double cumTp_u = px->
cumTp_exp.at(j).at(u);
1153 double raw_amount_carbon=0;
1154 double carbonPrice_y = carbonPrice_NOW;
1155 for (
int y=0; y<=cumTp_u; y++) {
1156 double expectedVHa_y =
getVHaByYear(px, j, y, futureExtraBiomass_ratio[j][y], regId);
1157 double carbonPayment_y;
1158 carbonPayment_y = carbonPrice_y*ir*(co2content_inventory*expectedVHa_y - FutureCarbonContent_opt[y]);
1170 raw_amount_carbon += carbonPayment_y*pow(1+ir,cumTp_u - y);
1174 double vHa = px->
vHa_exp.at(j).at(u);
1175 double finalHarvestFlag =
gfd(
"finalHarvestFlag",regId,ft,dc);
1183 double raw_amount = finalHarvestFlag*pw_exp*vHa*
app(
priProducts[pp],ft,dc);
1187 raw_amount_carbon += pickling * carbonPrice_y * finalHarvestFlag * vHa *
app(
priProducts[pp],ft,dc);
1192 double anualised_amount = anualised_amount_timber + anualised_amount_carbon;
1195 if (anualised_amount>=expReturns_hartman) {
1196 expReturns_hartman=anualised_amount;
1197 expReturns_carbon_hartman = anualised_amount_carbon;
1198 expReturns_timber_hartman = anualised_amount_timber;
1208 px->
optDc.push_back(optDc_hartman);
1213 double cumMort = 1-px->
cumAlive_exp.at(j).at(optDc_hartman);
1215 double origExpReturns_hartman = expReturns_hartman;
1216 expReturns_hartman = origExpReturns_hartman * (1.0 - ra*cumMort);
1223 expectedReturns.at(j) = expReturns_hartman;
1224 rotation_dc.at(j) = optDc_hartman;
1233 for(uint j=0;j<
fTypes.size();j++){
1237 double harvestedAreaForThisFT =
vSum(px->
hArea.at(j))*(1-fArea_decr_rel);
1238 deltaAreas[j][
fTypes.size()] +=
vSum(px->
hArea.at(j))*(fArea_decr_rel);
1239 vector<double> corrExpectedReturns(
fTypes.size(),0.0);
1240 vector<double> polBal_fiSub_unitvalue(
fTypes.size(),0.0);
1244 for(uint j2=0;j2<
fTypes.size();j2++){
1246 double invTransCost = (
gfd(
"invTransCost",regId,ft,ft2,
DATA_NOW) *
gfd(
"pol_fiSub",regId,ft,ft2,
DATA_NOW) ) ;
1247 polBal_fiSub_unitvalue[j2] =
gfd(
"invTransCost",regId,ft,ft2,
DATA_NOW) * (
gfd(
"pol_fiSub",regId,ft,ft2,
DATA_NOW) - 1) ;
1248 corrExpectedReturns[j2] = (expectedReturns[j2]/
ir)-invTransCost;
1272 int optFtChosen =
getMaxPos(corrExpectedReturns);
1276 thisYearRegAreas[
getMaxPos(corrExpectedReturns)] += harvestedAreaForThisFT*mr;
1277 thisYearRegAreas[
getMaxPos(expectedReturns)] += fArea_incr*mr/((double)
fTypes.size());
1279 deltaAreas[j][
getMaxPos(corrExpectedReturns)] += harvestedAreaForThisFT*mr;
1282 polBal_fiSub[
getMaxPos(corrExpectedReturns)] += harvestedAreaForThisFT*mr*polBal_fiSub_unitvalue[
getMaxPos(corrExpectedReturns)];
1290 double hAreaThisFtGroup =
findMap(hAreaByFTypeGroup,parentFt);
1291 double hRatio = 1.0;
1292 if(hAreaThisFtGroup>0){
1294 hRatio = harvestedAreaForThisFT/hAreaThisFtGroup;
1297 hRatio = 1.0/nFtChilds;
1299 thisYearRegAreas[j] += totalHarvestedArea*(1-mr)*freq*hRatio;
1300 thisYearRegAreas[j] += fArea_incr*(1-mr)*freq*hRatio;
1303 for(uint j2=0;j2<
fTypes.size();j2++){
1304 deltaAreas[j2][j] +=
vSum(px->
hArea.at(j2))*(1-fArea_decr_rel)*(1-mr)*freq*hRatio;
1306 deltaAreas[
fTypes.size()][j] += fArea_incr*(1-mr)*freq*hRatio;
1316 if(mortRatePath > 0){
1320 vector <double> siblingAreas;
1321 for(uint j2=0;j2<siblings.size();j2++){
1322 if(siblings[j2]==ft){
1323 siblingAreas.push_back(0.0);
1328 siblingAreas.push_back(thisSiblingArea);
1331 double areaAllSiblings =
vSum(siblingAreas);
1332 thisYearRegAreas[j] += harvestedAreaForThisFT*(1-mr)*(1-mortRatePath);
1333 deltaAreas[j][j] += harvestedAreaForThisFT*(1-mr)*(1-mortRatePath);
1335 if(areaAllSiblings>0.0){
1336 for(uint j2=0;j2<siblings.size();j2++){
1344 thisYearRegAreas[
getPos(siblings[j2],
fTypes)] += harvestedAreaForThisFT*(1.0-mr)*(mortRatePath)*(siblingAreas[j2]/areaAllSiblings);
1345 deltaAreas[j][
getPos(siblings[j2],
fTypes)] += harvestedAreaForThisFT*(1.0-mr)*(mortRatePath)*(siblingAreas[j2]/areaAllSiblings);
1347 }
else if (siblings.size()>1) {
1348 for(uint j2=0;j2<siblings.size();j2++){
1349 if (siblings[j2] != ft){
1350 thisYearRegAreas[
getPos(siblings[j2],
fTypes)] += harvestedAreaForThisFT*(1.0-mr)*(mortRatePath)* 1.0 / (( (float) siblings.size())-1.0);
1351 deltaAreas[j][
getPos(siblings[j2],
fTypes)] += harvestedAreaForThisFT*(1.0-mr)*(mortRatePath)* 1.0 / (( (float) siblings.size())-1.0);
1356 thisYearRegAreas[j] += harvestedAreaForThisFT*(1.0-mr);
1357 deltaAreas[j][j] += harvestedAreaForThisFT*(1.0-mr);
1362 thisYearRegAreas[j] += newAreaThisFt;
1363 deltaAreas[
fTypes.size()][j] += newAreaThisFt;
1365 if(! (thisYearRegAreas[j] >= 0.0) ){
1366 msgOut(
MSG_ERROR,
"thisYearRegAreas[j] is not non negative (j: "+
i2s(j)+
", thisYearRegAreas[j]: "+
i2s( thisYearRegAreas[j])+
").");
1374 for(uint j=0;j<
fTypes.size();j++){
1375 px->
area.at(j).at(0) += thisYearRegAreas.at(j);
1379 double totalRegArea =
vSum(thisYearRegAreas);
1380 if (! (totalRegArea==0.0 && totalHarvestedArea==0.0)){
1381 double ratio = totalRegArea / totalHarvestedArea ;
1383 msgOut(
MSG_CRITICAL_ERROR,
"Sum of regeneration areas not equal to sum of harvested area in runManagementModel()!");
1392 if (-fArea_diff > regHArea){
1393 msgOut(
MSG_WARNING,
"In region "+
i2s(regId) +
" the exogenous area decrement ("+
d2s(-fArea_diff) +
" ha) is higger than the harvesting ("+
d2s(regHArea) +
" ha). Ratio forced to 1.");
1396 for(uint j=0;j<
fTypes.size();j++){
1397 double debug = polBal_fiSub[j]/1000000;
1398 sfd((polBal_fiSub[j]/1000000.0),
"polBal_fiSub",regId,
fTypes[j],
"",
DATA_NOW,
true);
1427 msgOut(
MSG_INFO,
"Cashing model settings that may change every year..");
1444 for(uint i=0;i<
regIds2.size();i++){
1447 for (uint j=0;j<rpx.size();j++){
1450 rpx[j]->vol.clear();
1451 for(uint y=0;y<
fTypes.size();y++){
1452 vector <double> vol_byu;
1454 for (uint z=0;z<
dClasses.size();z++){
1457 double pxArea = rpx[j]->getDoubleValue(
"forArea_"+
fTypes[y],
true);
1461 double pxVol = regForArea ? regVol * pxArea/regForArea: 0;
1463 vol_byu.push_back(pxVol);
1465 rpx[j]->vol.push_back(vol_byu);
1505 for(uint r=0;r<
regIds2.size();r++){
1509 for(uint f=0;f<
fTypes.size();f++){
1512 double sStDev =
gfd(
"sStDev",
regIds2[r],ft,
"");
1513 vector<double> vols;
1514 vector<double> diffVols;
1515 vector<double> diffVols_rescaled;
1516 double diffVols_rescaled_deviation;
1517 double diffVols_rescaled_deviation_rescaled;
1520 vector<double> multipliers;
1522 double vol_max, rescale_ratio_avg, rescale_ratio_sd;
1523 double diffVols_avg, diffVols_rescaled_avg;
1524 double diffVols_rescaled_sd;
1526 for (uint p=0;p<rpx.size();p++){
1528 vols.push_back(
vSum(px->
vol[f]));
1532 for(uint p=0;p<vols.size();p++){
1533 diffVols.push_back(vol_max-vols[p]);
1536 diffVols_avg =
getAvg(diffVols);
1537 rescale_ratio_avg = (diffVols_avg != 0.0) ? agr/diffVols_avg : 1.0;
1538 for(uint p=0;p<diffVols.size();p++){
1539 diffVols_rescaled.push_back(diffVols[p]*rescale_ratio_avg);
1541 diffVols_rescaled_avg =
getAvg(diffVols_rescaled);
1542 diffVols_rescaled_sd =
getSd(diffVols_rescaled,
false);
1544 rescale_ratio_sd = (diffVols_rescaled_sd != 0.0) ? sStDev/diffVols_rescaled_sd : 1.0;
1545 for(uint p=0;p<diffVols_rescaled.size();p++){
1546 diffVols_rescaled_deviation = diffVols_rescaled[p] - diffVols_rescaled_avg;
1547 diffVols_rescaled_deviation_rescaled = diffVols_rescaled_deviation * rescale_ratio_sd;
1548 final_value = diffVols_rescaled_avg + diffVols_rescaled_deviation_rescaled;
1549 multiplier = (agr != 0.0) ? min(1.6, max(0.4,final_value/agr)) : 1.0;
1554 multipliers.push_back(multiplier);
1559 double avgMultipliers =
getAvg(multipliers);
1560 double sdMultipliers =
getSd(multipliers,
false);
1561 if ( avgMultipliers < 0.9 || avgMultipliers > 1.1){
1564 if ( ( sdMultipliers - (sStDev/agr) ) < -0.5 || ( sdMultipliers - (sStDev/agr) ) > 0.5 ){
1565 double cv = sStDev/agr;
1566 msgOut(
MSG_CRITICAL_ERROR,
"The sd of multipliers of ft "+ ft +
" for the region " +
i2s(rId) +
" is not equal to the spatial cv for the region!");
1581 for(uint i=0;i<
regIds2.size();i++){
1582 vector<double> deathBiomass;
1583 for(uint j=0;j<
fTypes.size();j++){
1585 deathBiomass.push_back(deathBiomass_ft);
1588 vector<double>qProducts;
1592 qProducts.push_back(int_exports);
1597 qProducts.push_back(consumption);
1607 for(
int y=currentYear;y>currentYear-30;y--){
1608 for(uint i=0;i<
regIds2.size();i++){
1609 for(uint j=0;j<
fTypes.size();j++){
1610 for (uint u=0;u<
dClasses.size();u++){
1633 for(uint i=0;i<
regIds2.size();i++){
1636 for (uint p=0;p<rpx.size();p++){
1638 double pxid= px->
getID();
1639 for(uint j=0;j<
fTypes.size();j++){
1641 vector <double> tempAreas;
1642 vector <double> areasByFt;
1644 for (uint u=0;u<
dClasses.size();u++){
1647 double regRegVolumes =
gfd(
"vReg",
regIds2[i],ft,
"");
1648 double newVReg = regionArea ? regRegVolumes*pxArea/regionArea : 0.0;
1649 double tp_u0 = px->
tp.at(j).at(0);
1650 double entryVolHa =
gfd(
"entryVolHa",
regIds2[i],ft,
"");
1651 double tempArea = (newVReg*1000000.0/entryVolHa)*tp_u0;
1652 tempAreas.push_back(tempArea);
1655 double dcVol = px->
vol_l.at(j).at(u)*1000000.0;
1656 double dcVHa = px->
vHa.at(j).at(u);
1658 if(dcVol < 0.0 || dcVHa < 0.0){
1662 double tempArea = dcVHa?dcVol/dcVHa:0;
1663 tempAreas.push_back(tempArea);
1667 double sumTempArea =
vSum(tempAreas);
1671 double normCoef = sumTempArea?pxArea/ sumTempArea:0;
1678 for (uint u=0;u<
dClasses.size();u++){
1679 areasByFt.push_back(tempAreas.at(u)*normCoef);
1683 double ratio =
vSum(areasByFt)/ pxArea;
1684 if(ratio < 0.99999999999 || ratio > 1.00000000001) {
1689 px->
area_l.push_back(areasByFt);
1711 for(uint r2= 0; r2<
regIds2.size();r2++){
1715 for (uint p=0;p<
regPx.size();p++){
1725 for(uint j=0;j<
fTypes.size();j++){
1734 double pathMort_now, pathMort_t0;
1747 vector <double> cumTp_temp;
1748 vector <double> vHa_temp;
1749 vector <double> cumAlive_temp;
1750 vector <double> cumTp_exp_temp;
1751 vector <double> vHa_exp_temp;
1752 vector <double> cumAlive_exp_temp;
1755 for (uint u=0; u<
dClasses.size(); u++){
1757 double cumTp_u, cumTp_u_exp, cumTp_u_noExp, cumTp_u_fullExp;
1758 double tp, tp_exp, tp_noExp, tp_fullExp;
1759 double vHa_u, vHa_u_exp, vHa_u_noExp, vHa_u_fullExp, beta, beta_exp, beta_noExp, beta_fullExp, mort, mort_exp, mort_noExp, mort_fullExp;
1760 double cumAlive_u, cumAlive_exp_u;
1763 double additionalMortality_multiplier_reg_ft_dc =
gfd(
"additionalMortality_multiplier",regId,ft,dc,
DATA_NOW);
1764 double additionalMortality_multiplier_reg_ft_dc_t0 =
gfd(
"additionalMortality_multiplier",regId,ft,dc,
firstYear);
1765 double additionalMortality_now = px->
getSTData(
"additionalMortality", ft,
DATA_NOW, dc, 0.0)*additionalMortality_multiplier_reg_ft_dc;
1766 double additionalMortality_t0 = px->
getSTData(
"additionalMortality", ft,
firstYear, dc, 0.0)*additionalMortality_multiplier_reg_ft_dc_t0;
1776 cumTp_temp.push_back(cumTp_u);
1777 vHa_temp.push_back(vHa_u);
1778 cumTp_exp_temp.push_back(cumTp_u);
1779 vHa_exp_temp.push_back(vHa_u);
1780 cumAlive_temp.push_back(cumAlive_u);
1781 cumAlive_exp_temp.push_back(cumAlive_u);
1786 tp =
gfd(
"tp",regId,ft,
dClasses[u-1],thisYear)*tp_multiplier_now;
1787 cumTp_u = cumTp_temp[u-1] + tp;
1796 vHa_u =
gfd(
"entryVolHa",regId,ft,
"",thisYear);
1799 mort = 1-pow(1-min(
gfd(
"mortCoef",regId,ft,
dClasses[u-1],thisYear)*mortCoef_multiplier_now+pathMort_now+additionalMortality_now,1.0),tp);
1800 beta =
gfd(
"betaCoef",regId,ft,dc, thisYear)*betaCoef_multiplier_now;
1801 vHa_u = vHa_temp[u-1]*beta*(1-mort);
1803 cumAlive_u = max(0.,cumAlive_temp[u-1]*(1-mort));
1804 cumAlive_temp.push_back(cumAlive_u);
1805 cumTp_temp.push_back(cumTp_u);
1806 vHa_temp.push_back(vHa_u);
1816 cumTp_u_exp = cumTp_exp_temp[u-1]+tp_exp;
1817 cumTp_exp_temp.push_back(cumTp_u_exp);
1822 mort_exp = 1-pow(1-min(
gfd(
"mortCoef",regId,ft,
dClasses[u-1],
firstYear)*mortCoef_multiplier_t0+pathMort_t0+additionalMortality_t0,1.0),tp_exp);
1823 beta_exp =
gfd(
"betaCoef",regId,ft,dc,
firstYear)*betaCoef_multiplier_t0;
1824 vHa_u_exp = vHa_exp_temp[u-1]*beta_exp*(1-mort_exp);
1827 double tp_multiplier_dynamic = px->
getMultiplier(
"tp_multiplier",ft,thisYear+ceil(cumTp_exp_temp[u-1]));
1828 tp_noExp =
gfd(
"tp",regId,ft,
dClasses[u-1])*tp_multiplier_now;
1829 cumTp_u_noExp = cumTp_exp_temp[u-1]+tp_noExp;
1830 tp_fullExp =
gfd(
"tp",regId,ft,
dClasses[u-1],thisYear+ceil(cumTp_exp_temp[u-1]))*tp_multiplier_dynamic ;
1831 cumTp_u_fullExp = cumTp_exp_temp[u-1]+tp_fullExp ;
1832 cumTp_u_exp = cumTp_u_fullExp*expType+cumTp_u_noExp*(1-expType);
1833 cumTp_exp_temp.push_back(cumTp_u_exp);
1835 vHa_u_noExp =
gfd(
"entryVolHa",regId,ft,
"",
DATA_NOW);
1836 vHa_u_fullExp =
gfd(
"entryVolHa",regId,ft,
"",thisYear+ceil(cumTp_u));
1837 vHa_u_exp = vHa_u_fullExp*expType+vHa_u_noExp*(1-expType);
1840 mort_noExp = 1-pow(1-min(1.0,
gfd(
"mortCoef",regId,ft,
dClasses[u-1],
DATA_NOW)*mortCoef_multiplier_now+pathMort_now+additionalMortality_now), tp_noExp);
1841 double mortCoef_multiplier_dynamic = px->
getMultiplier(
"mortCoef_multiplier",ft,thisYear+ceil(cumTp_exp_temp[u-1]));
1842 double pathMort_dynamic = px->
getPathMortality(ft,dc,thisYear+ceil(cumTp_exp_temp[u-1]));
1843 double additionalMortality_multiplier_reg_ft_dc_dynamic =
gfd(
"additionalMortality_multiplier",regId,ft,dc, thisYear+ceil(cumTp_exp_temp[u-1]));
1844 double additionalMortality_dynamic = px->
getSTData(
"additionalMortality", ft, thisYear+ceil(cumTp_exp_temp[u-1]), dc, 0.0)*additionalMortality_multiplier_reg_ft_dc_dynamic;
1845 mort_fullExp = 1-pow(1-min(1.0,
gfd(
"mortCoef",regId,ft,
dClasses[u-1],thisYear+ceil(cumTp_exp_temp[u-1]))*mortCoef_multiplier_dynamic+pathMort_dynamic+additionalMortality_dynamic),tp_fullExp);
1854 beta_noExp =
gfd(
"betaCoef",regId,ft,dc,
DATA_NOW)*betaCoef_multiplier_now;
1855 double betaCoef_multiplier_dynamic = px->
getMultiplier(
"betaCoef_multiplier",ft,thisYear+ceil(cumTp_u));
1856 beta_fullExp =
gfd(
"betaCoef",regId,ft,dc, thisYear+ceil(cumTp_u))*betaCoef_multiplier_dynamic;
1857 mort_exp = mort_fullExp*expType+mort_noExp*(1-expType);
1858 beta_exp = beta_fullExp*expType+beta_noExp*(1-expType);
1859 vHa_u_exp = vHa_exp_temp[u-1]*beta_exp*(1-mort_exp);
1862 vHa_exp_temp.push_back(vHa_u_exp);
1863 cumAlive_exp_u = max(0.,cumAlive_exp_temp[u-1]*(1-mort_exp));
1864 cumAlive_exp_temp.push_back(cumAlive_exp_u);
1880 px->
cumTp.push_back(cumTp_temp);
1881 px->
vHa.push_back(vHa_temp);
1882 px->
cumAlive.push_back(cumAlive_temp);
1883 px->
cumTp_exp.push_back(cumTp_exp_temp);
1884 px->
vHa_exp.push_back(vHa_exp_temp);
1906 for(uint r2= 0; r2<
regIds2.size();r2++){
1909 for (uint p=0;p<
regPx.size();p++){
1947 msgOut(
MSG_INFO,
"Starting cashing on pixel spatial-level exogenous data");
1950 for(uint r2= 0; r2<
regIds2.size();r2++){
1953 for (uint p=0;p<
regPx.size();p++){
1961 for(uint j=0;j<
fTypes.size();j++){
1963 vector <double> tp_byu;
1964 vector <double> beta_byu;
1965 vector <double> mort_byu;
1966 vector <double> addMort_byu;
1978 for (uint u=0; u<
dClasses.size(); u++){
1981 double additionalMortality_multiplier_reg_ft_dc =
gfd(
"additionalMortality_multiplier",regId,ft,dc,
DATA_NOW);
1982 double additionalMortality = px->
getSTData(
"additionalMortality", ft,
DATA_NOW, dc, 0.0)*additionalMortality_multiplier_reg_ft_dc;
1983 double tp, beta_real, mort_real;
1993 mort_real = min(u?
gfd(
"mortCoef",regId,ft,
dClasses[u],
DATA_NOW)*mortCoef_multiplier_now+pathMortality :0,1.0);
1994 tp_byu.push_back(tp);
1995 beta_byu.push_back(beta_real);
1996 mort_byu.push_back(mort_real);
1997 addMort_byu.push_back(additionalMortality);
1999 px->
tp.push_back(tp_byu);
2000 px->
beta.push_back(beta_byu);
2001 px->
mort.push_back(mort_byu);
2002 px->
addMort.push_back(addMort_byu);
2003 px->
avalCoef.push_back(avalCoef_ft);
2011 msgOut(
MSG_INFO,
"Starting computing inventary available for this year..");
2016 for(uint i=0;i<
regIds2.size();i++){
2022 vector <double> in_deathTimber_reg(
priProducts.size(),0.);
2023 for (uint p=0;p<
regPx.size();p++){
2030 for(uint ft=0;ft<
fTypes.size();ft++){
2031 for(uint dc=0;dc<
dClasses.size();dc++){
2036 in_reg.at(pp) += in;
2042 vector<string> priProducts_vector;
2046 in_deathTimber_reg.at(pp) += in_deathMortality;
2056 spd(in_deathTimber_reg.at(pp),
"in_deathTimber",r2,priProducts[pp],
DATA_NOW,
true);
2058 if (in_reg.at(pp) < -0.0){
2073 for (uint i=0; i<nbounds; i++){
2074 vector<int> concernedPriProducts = concernedPriProductsTotal[i];
2082 double bound_total = bound_alive + bound_deathTimber;
2096 for(uint i=0;i<
regIds2.size();i++){
2099 for (uint p=0;p<rpx.size();p++){
2101 double pxid= px->
getID();
2102 for(uint j=0;j<
fTypes.size();j++){
2104 double forArea =
vSum(px->
area.at(j));
2116 map<int,double> forestArea;
2117 pair<int,double > forestAreaPair;
2120 int nFTypes = fTypes.size();
2121 int nL2Regions = l2Regions.size();
2122 for(
int i=0;i<nL2Regions;i++){
2123 int regId = l2Regions[i];
2125 for(
int j=0;j<nFTypes;j++){
2126 string ft = fTypes[j];
2133 for(uint z=0;z<rpx.size();z++){
2139 double newValue = oldValue-hArea+regArea;
2140 double areaNetOfRegeneration = oldValue-hArea;
2142 if (areaNetOfRegeneration<0.0){
2149 rpx[z]->changeValue(
"forArea_"+ft, newValue*10000);
2161 int nFTypes = fTypes.size();
2162 int nL2Regions = l2Regions.size();
2163 for(
int i=0;i<nL2Regions;i++){
2164 int regId = l2Regions[i];
2166 for(
int j=0;j<nFTypes;j++){
2167 string ft = fTypes[j];
2168 for(uint z=0;z<rpx.size();z++){
2170 double vol =
vSum(px->
vol.at(j));
2173 rpx[z]->changeValue(
"vol_"+ft, vol);
2176 rpx[z]->changeValue(
"expectedReturns_"+ft, expectedReturns);
2183 for(
int j=0;j<nFTypes;j++){
2184 string ft = fTypes[j];
2201 map<tr1::array<string, 4>,
double> hVol_byPrd;
2204 for(uint r2= 0; r2<
regIds2.size();r2++){
2210 for(uint j=0;j<
fTypes.size();j++){
2213 double regArea = 0.;
2214 double sumAreaByFt = 0.;
2215 double pxForAreaByFt = 0.;
2217 double sumSumHProductivity = 0.;
2218 double sumHArea = 0.;
2220 for (uint u=0; u<
dClasses.size(); u++){
2226 double vMortAdd = 0;
2227 double sumHProductivity = 0.;
2228 for (uint pi=0;pi<
regPx.size();pi++){
2230 vol += px->
vol.at(j).at(u);
2231 hV += px->
hVol.at(j).at(u);
2232 hArea += px->
hArea.at(j).at(u);
2233 vMort += px->
vMort.at(j).at(u);
2234 vMortAdd += px->
vMortAdd.at(j).at(u);
2240 tr1::array<string, 4> hVKey = {
i2s(regId), ft, dc, pr};
2250 double hProductivityByDc = hArea ? sumHProductivity/hArea : 0.;
2251 sumSumHProductivity += (hProductivityByDc*hArea);
2253 sfd(hProductivityByDc,
"hProductivityByDc",regId,ft,dc,
DATA_NOW,
true);
2254 sfd(hArea,
"harvestedArea",regId,ft,dc,
DATA_NOW,
true);
2256 sfd(vMortAdd,
"vMortAdd",regId,ft,dc,
DATA_NOW,
true);
2257 double vol_1 =
gfd(
"vol",regId,ft,dc,currentYear-1);
2266 double hProductivity = sumHArea? sumSumHProductivity / sumHArea : 0.;
2267 sfd(hProductivity,
"hProductivity",regId,ft,
"",
DATA_NOW,
true);
2269 for (uint p=0;p<
regPx.size();p++){
2271 vReg += px->
vReg.at(j);
2275 sumAreaByFt += pxForAreaByFt;
2277 if(! (sumAreaByFt >= 0.0) ){
2282 sfd(regArea,
"regArea",regId,ft,
"",
DATA_NOW,
true);
2283 sfd(sumAreaByFt,
"forArea",regId,ft,
"",
DATA_NOW,
true);
2286 for (uint p=0;p<
regPx.size();p++){
2288 double totPxForArea =
vSum(px->
area);
2291 double totPxForArea_debug = 0.0;
2292 for(uint j=0;j<
fTypes.size();j++){
2294 totPxForArea_debug += (px->
getDoubleValue(
"forArea_"+ft,
true)/10000);
2297 if ( (totPxForArea - totPxForArea_debug) > 0.0001 || (totPxForArea - totPxForArea_debug) < -0.0001 ){
2298 cout <<
"*** ERROR: area discrepance in pixel " << px->
getID() <<
" of " << (totPxForArea - totPxForArea_debug) <<
" ha!" << endl;
2299 msgOut(
MSG_CRITICAL_ERROR,
"Total forest area in pixel do not coincide if token from layer forArea or (pixel) vector area!");
2316 for(uint r2= 0; r2<
regIds2.size();r2++){
2319 double totRegionForArea = 0.;
2320 double totSumExpRet = 0.;
2321 vector <double> totSumExpRet_byFTParent(parentFtypes.size(),0.0);
2322 vector <double> totSumExpRet_byFTypes(
fTypes.size(),0.0);
2325 for (uint p=0;p<
regPx.size();p++){
2329 totRegionForArea += pxForArea;
2331 for(uint i=0;i<parentFtypes.size();i++){
2334 vector<double> pxExpReturnsByChilds(childPos.size(),0.0);
2335 for(uint j=0;j<childPos.size();j++){
2339 pxExpReturnsByChilds.at(j) = (childIds.at(j) ==
"ash") ? 0.0 : pxExpReturn_singleFt;
2341 totSumExpRet_byFTypes.at(childPos[j]) += pxExpReturn_singleFt*pxForArea;
2343 totSumExpRet_byFTParent[i] +=
getMax(pxExpReturnsByChilds)*pxForArea;
2345 totSumExpRet += bestPxExpectedRet * pxForArea;
2349 for(uint i=0;i<parentFtypes.size();i++){
2351 for(uint j=0;j<childPos.size();j++){
2353 sfd(totSumExpRet_byFTypes.at(childPos[j]),
"sumExpReturns",regId,
fTypes.at(childPos[j]),
"",
DATA_NOW,
true);
2354 sfd(totSumExpRet_byFTypes.at(childPos[j])/totRegionForArea,
"expReturns",regId,
fTypes.at(childPos[j]),
"",
DATA_NOW,
true);
2357 sfd(totSumExpRet_byFTParent.at(i),
"sumExpReturns",regId,parentFtypes[i],
"",
DATA_NOW,
true);
2358 sfd(totSumExpRet_byFTParent.at(i)/totRegionForArea,
"expReturns",regId,parentFtypes[i],
"",
DATA_NOW,
true);
2362 sfd(totSumExpRet,
"sumExpReturns",regId,
"",
"",
DATA_NOW,
true);
2363 sfd(totSumExpRet/totRegionForArea,
"expReturns",regId,
"",
"",
DATA_NOW,
true);
2369 for(uint r2= 0; r2<
regIds2.size();r2++){
2372 double totalForArea = 0.0;
2373 double invadedArea = 0.0;
2374 for (uint p=0;p<
regPx.size();p++){
2377 for(uint j=0;j<
fTypes.size();j++){
2378 for (uint u=0; u<
dClasses.size(); u++){
2385 invadedArea +=
vSum(px->
area)*invaded;
2387 sfd(invadedArea/totalForArea,
"totalShareInvadedArea",regId,
"",
"",
DATA_NOW,
true);
2413 for(uint i=0;i<
regIds2.size();i++){
2414 for(uint j=0;j<
fTypes.size();j++){
2434 for (uint r1=0;r1<
l2r.size();r1++){
2435 for (uint r2=0;r2<
l2r[r1].size();r2++){
2436 int rfrom=
l2r[r1][r2];
2437 double distQProd = 0.0;
2438 for (uint r3=0;r3<
l2r[r1].size();r3++){
2439 int rto =
l2r[r1][r3];
2465 double fullExpWPrice = (curLocPrice*(worldFutPrice/worldCurPrice)*sl+worldFutPrice*sa)/(sa+sl);
2466 double curWPrice = (curLocPrice*sl+worldCurPrice*sa)/(sl+sa);
2467 return curWPrice * (1-expCoef) + fullExpWPrice * expCoef;
2484 int nFTypes =
fTypes.size();
2492 pxC +=
regPx.size();
2495 static vector<vector<vector<double>>> original_vols(pxC, vector<vector<double>>(nFTypes, vector<double>(nDC, 0.0)));
2499 for(uint i=0;i<
fTypes.size();i++){
2500 for (uint u=0; u<
dClasses.size(); u++){
2506 for (uint p=0;p<
regPx.size();p++){
2508 original_vols[pxC_loc][i][u] += px->
vol[i][u];
2514 for(uint i=0;i<
fTypes.size();i++){
2516 for(uint o=0;o<
fTypes.size();o++){
2518 for (uint u=1; u<
dClasses.size(); u++){
2519 string layerName =
"spInput#vol#"+fto+
"#"+fti+
"#"+
i2s(u);
2525 for (uint p=0;p<
regPx.size();p++){
2527 double vol_transfer = min(px->
getDoubleValue(layerName,
true)/1000000,px->
vol[i][u]) ;
2528 px->
vol[i][u] -= vol_transfer;
2529 px->
vol[o][u] += vol_transfer;
2559 vector<vector<vector<double>>> original_areas(pxC, vector<vector<double>>(nFTypes, vector<double>(nDC, 0.0)));
2560 for(uint i=0;i<
fTypes.size();i++){
2561 for (uint u=0; u<
dClasses.size(); u++){
2567 for (uint p=0;p<
regPx.size();p++){
2569 original_areas[pxC_loc][i][u] += px->
area_l[i][u];
2578 vector<vector<vector<double>>> transferred_areas(pxC, vector<vector<double>>(nFTypes, vector<double>(nFTypes, 0.0)));
2580 for(uint i=0;i<
fTypes.size();i++){
2582 for(uint o=0;o<
fTypes.size();o++){
2584 for (uint u=1; u<
dClasses.size(); u++){
2585 string layerName =
"spInput#vol#"+fto+
"#"+fti+
"#"+
i2s(u);
2592 for (uint p=0;p<
regPx.size();p++){
2594 double vol_i_orig = original_vols[pxC_loc][i][u];
2595 double vol_transfer = vol_i_orig?px->
getDoubleValue(layerName,
true)/1000000:0.0;
2596 double area_i_orig = original_areas[pxC_loc][i][u];
2597 double area_transfer = vol_i_orig?area_i_orig*vol_transfer/vol_i_orig:0.0;
2598 px->
area_l[i][u] -= area_transfer;
2600 px->
area_l[o][u] += area_transfer;
2602 transferred_areas[pxC_loc][i][o] += area_transfer;
2617 for (uint p=0;p<
regPx.size();p++){
2619 for(uint i=0;i<
fTypes.size();i++){
2620 for(uint o=0;o<
fTypes.size();o++){
2621 double area_i_orig = 0.0;
2622 for (uint u=1; u<
dClasses.size(); u++){
2623 area_i_orig += original_areas[pxC_loc][i][u];
2625 double area_transfer_u0 = area_i_orig?original_areas[pxC_loc][i][0]*(transferred_areas[pxC_loc][i][o]/area_i_orig):0.0;
2626 px->
area_l[i][0] -= area_transfer_u0 ;
2628 px->
area_l[o][0] += area_transfer_u0 ;
2637 for(uint i=0;i<
fTypes.size();i++){
2638 string fti_id =
fTypes[i];
2640 int ft_memType = fti->
memType;
2641 string ft_layerName = fti->
forLayer;
2645 if(ft_memType==3 ||ft_memType==2){
2650 for (uint p=0;p<
regPx.size();p++){
2652 double area_px =
vSum(px->
area[i]);
2664 cout <<
"Printing debug information on initial regional inventories and supplies.." << endl;
2665 cout <<
"Reg\tProduct\t\tInv\tSt\tSa\tSl" << endl;
2666 for(uint r1=0;r1<
l2r.size();r1++){
2667 for(uint r2c=0;r2c<
l2r[r1].size();r2c++){
2669 int r2 =
l2r[r1][r2c];
2674 cout << r2 <<
"\t" <<
priProducts[p] <<
"\t\t" << inv <<
"\t" << st <<
"\t" << sl <<
"\t" << sa << endl;
2702 for(uint y = currentYear-1; y>currentYear-1-maxYears; y--){
2703 for (
int u =
dClasses.size()-1; u>=0; u--){
2705 for (uint f=0; f<
fTypes.size(); f++){
2710 double sum_total_st_allocable = 0;
2712 for(uint ap=0;ap<allocableProducts.size();ap++){
2713 sum_total_st_allocable += total_st.at(allocableProducts[ap]);
2715 for(uint ap=0;ap<allocableProducts.size();ap++){
2716 double allocableShare = sum_total_st_allocable?total_st.at(allocableProducts[ap])/sum_total_st_allocable:0.0;
2717 double allocated = min(total_st[allocableProducts[ap]],deathTimber*allocableShare);
2719 total_st[allocableProducts[ap]] -= allocated;
2739 for(uint r2=0; r2<
regIds2.size();r2++){
2744 double pol_mktDirInt_s =
gpd(
"pol_mktDirInt_s",regId,
priProducts[p]);
2745 double pol_mktDirInt_s_1 =
gpd(
"pol_mktDirInt_s",regId,
priProducts[p],previousYear);
2756 double gamma = in>in_1 ? gamma_incr: gamma_decr;
2757 double polBal_mktDirInt_s = pc * (1-pol_mktDirInt_s) * sc;
2758 double surplus_prod = pc * sc -
2761 * pow(pol_mktDirInt_s,-1)
2762 * pow(sc_1,-1/sigma)
2763 * pow(in_1/in,gamma/sigma)
2765 * pow(sc,(sigma+1)/sigma);
2771 double pol_mktDirInt_d =
gpd(
"pol_mktDirInt_d",regId,
secProducts[p]);
2772 double pol_mktDirInt_d_1 =
gpd(
"pol_mktDirInt_d",regId,
secProducts[p],previousYear);
2783 double polBal_mktDirInt_d = pc * (pol_mktDirInt_d - 1) * dc;
2784 double polBal_trSub = m * (pol_trSub-1) * sl;
2786 double surplus_cons = pc_1
2788 * pow(dc_1,-1/sigma)
2789 * pow(pol_mktDirInt_d,-1)
2791 * (pow(dc,(sigma+1)/sigma) - pow(trs*dc0,(sigma+1)/sigma))
2792 - (dc - trs * dc0) * pc ;
2801 double polBal_tcSub = 0;
2802 vector<ModelRegion*> siblings = reg->
getSiblings();
2803 for(uint rTo=0;rTo<siblings.size();rTo++){
2804 int regId2 = siblings[rTo]->getRegId();
2807 polBal_tcSub += ct*(pol_tcSub-1) * rt;
2821 if(dc ==
dClasses.size()-1) {
return px->
cumTp[ft][dc] * 1.2;}
2822 else {
return (px->
cumTp[ft][dc]+px->
cumTp[ft][dc+1])/2; }
2840 for (
int u = 0; u<
dClasses.size(); u++){
2851 double cumTP_plus = px->
cumTp_exp[ft][u];
2852 double cumTP_minus = px->
cumTp_exp[ft][u-1];
2853 double volHa_plus = px->
vHa_exp[ft][u];
2854 double volHa_minus = px->
vHa_exp[ft][u-1];
2855 double slope = (volHa_plus - volHa_minus)/(cumTP_plus - cumTP_minus);
2857 vHa = (volHa_minus + (year-cumTP_minus)*slope)*(1+extraBiomass_ratio);
2865 double cumTP_max = px-> cumTp_exp[ft][dc];
2866 double time_elapsed = year - cumTP_max;
2867 double volHa_max = px->
cumTp_exp[ft][dc];
2871 vHa = volHa_max*pow(1-mortCoeff,time_elapsed)*(1+extraBiomass_ratio);
void updateImage(string layerName_h)
Add one layer to the system.
vector< vector< double > > deltaArea
void setSpModifier(const double &value, const int &ftindex)
Class to provide a simple integer-integer-string-string key in std maps.
int getPos(const K &element, const vector< K > &v, const int &msgCode_h=MSG_CRITICAL_ERROR)
Print an ERROR message, but don't stop the model.
The required data is for the current year.
void runBiologicalModule()
computes hV, hArea and new vol at end of year
vector< int > getAllocableProductIdsFromDeathTimber(const int ®Id_h, const string &ft, const string &dc, const int &harvesting_year, int request_year=DATA_NOW)
Returns the ids of the primary products that is possible to obtain using the timber recorded death in...
double getArea(const string &fType_h, const string &dClass_h)
Get area by ft and dc (from pixel->area matrix)
void initialiseDeathBiomassStocks(const vector< double > &deathByFt, const int ®Id)
Initialises the stocks of death biomass for the avgLive_* years before the simulation starts...
int getIntSetting(const string &name_h, int position=0, int reg=WORLD) const
vector< vector< int > > l2r
vector< vector< double > > area
vector< string > dClasses
vector< string > getForTypeParents()
vector< vector< double > > cumAlive_exp
This is the expected version of cumAlive, used for calculating profits.
bool getBoolSetting(const string &name_h, int position=0, int reg=WORLD) const
void runManagementModule()
computes regArea and expectedReturns
Do not actually output any message.
void initializePixelVolumes()
distribuite regional exogenous volumes to pixel volumes using corine land cover area as weight ...
vector< vector< int > > createCombinationsVector(const int &nItems)
Return a vector containing any possible combination of nItems items (including any possible subset)...
void sumRegionalForData()
computes vol, hV, harvestedArea, regArea and expReturns at reg level from the pixel level ...
vector< vector< double > > addMort
void spd(const double &value_h, const string &type_h, const int ®Id_h, const string &prodId_h, const int &year=DATA_NOW, const bool &allowCreate=false, const string &freeDim_h="") const
string i2s(const int &int_h) const
integer to string conversion
ThreadManager * MTHREAD
Pointer to the Thread manager.
K getMax(const vector< K > &v)
Returns the value of the maximum element in the vector (the last one in case of multiple equivalent m...
double getAvailableDeathTimber(const vector< string > &primProd_h, int regId_h, int year_h)
Returns the timber available for a given set of primary products as stored in the deathTimberInventor...
bool layerExist(const string &layerName_h, bool exactMatch=true) const
Return a pointer to a layer given its name.
void changeValue(const string &layerName_h, const double &value_h, const bool &setNoValueForZero=false)
Change the value of an existing layerMTHREAD->GIS->pack(parName, forName, dClass, year)...
double getValue(string layerName, int op=OP_SUM)
return the values of its own pixels for the specified layer. Possible operations: OP_SUM or OP_AVG ...
double calculateAnnualisedEquivalent(const double &amount_h, const int &years_h, const double &ir) const
Calculate the annual equivalent flow.
ModelData * MD
the model data object
void resetMapValues(map< K, V > &mymap, const V &value)
Reset all values stored in a map to the specified one.
void loadExogenousForestLayers(const string &what)
Set pixel volumes (what="vol") OR areas (what="area") by specific forest types as defined in gis laye...
double expType
Sampling derived expectation types of this agent (forest bilogical parameters: growth, mortality)
vector< double > avalCoef
Availability (of wood resources) coefficient. A [0,1] coefficient (new: by forest type) that reduces ...
vector< string > priProducts
string forestAreaChangeMethod
double getAvailableAliveTimber(const vector< string > &primProd_h, int regId_h)
Returns the timber available for a given set of primary products as stored in the px->vol_l vector...
vector< vector< double > > mort
void printOptLog(bool optimal, int &nIterations, double &obj)
Scheduler * SCD
the scheduler object (simulation-loops scheduler)
vector< vector< double > > vMort
vector< T > positionsToContent(const vector< T > &vector_h, const vector< int > &positions)
Return a vector of content from a vector and a vector of positions (int)
vector< vector< double > > vol_l
store the volumes of the previous year
void msgOut(const int &msgCode_h, const string &msg_h, const bool &refreshGUI_h=true) const
Overloaded function to print the output log.
vector< double > expectedAnnualIncome_carbon
Gis * GIS
GIS information and methods.
Thread manager. Responsable to manage the main thread and "speak" with the GUI.
void incrMapValue(map< K, V > &mymap, const K &key, const V &value, const int &error_level=MSG_CRITICAL_ERROR)
Increments a value stored in a map of the specified value, given the key.
vector< string > getStringVectorSetting(const string &name_h, int reg=WORLD) const
void resetPixelValues()
swap volumes->lagged_volumes and reset the other pixel vectors
double deathTimberInventory_get(const iisskey &thekey)
vector< vector< double > > vHa_exp
This is the expected version of vHa, used for calculating profits.
vector< string > pDClasses
vector< double > initialDc0Area
vector< vector< double > > cumAlive
Cumulative prob of remaining alive at beginnin of a given diam class.
vector< int > getForTypeChilds_pos(const string &forTypeId_h, bool all=false)
int getMaxPos(const vector< K > &v)
Returns the position of the maximum element in the vector (the last one in case of multiple equivalen...
void computeInventary()
in=f(vol_t-1)
vector< vector< double > > vHa
Volume at hectar by each diameter class [m^3/ha].
vector< double > inResByAnyCombination
Vector of inventory resource for each possible combination of primary products. This store both alive...
string getStringSetting(const string &name_h, int position=0, int reg=WORLD) const
void runMarketModule()
computes st (supply total) and pw (weighted price). Optimisation inside.
void registerTransports(const double &distQ, const int ®Id)
Registers the quantities emitted by transport of wood FROM a given region.
int getNForTypesChilds(const string &forTypeId_h)
void cacheSettings()
just cache exogenous settings from ModelData
Print an error message and stop the model.
void registerDeathBiomass(const double &value, const int ®Id, const string &fType)
Registers the "death" of a given amount of biomass, storing it in the deathBiomass map...
vector< vector< double > > vMortAdd
double getSTData(const string &parName, const string &forName, int year, const string &d2="", double naToReturn=RETNA)
vector< int > optDcChosen
double getDoubleValue(const string &layerName_h, const bool &returnZeroForNoValue=false) const
Return the value for a specific layer.
string d2s(const double &double_h) const
double to string conversion
vector< double > expectedAnnualIncome_timber
const double getForData(const string &type_h, const int ®Id_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW)
void updateOtherMapData()
update (if the layer exists) other gis-based data, as volumes and expected returns, taking them from the data in the px object
int vSum(const vector< int > &vector_h) const
void initialiseDeathTimber()
Set deathTimberInventory to zero for the previous years (under the hipotesis that we don't have advan...
double getPastRegArea(const int &ft_idx, const int &year)
vector< vector< double > > hProductivity
void initialiseProductsStocks(const vector< double > &qByProduct, const int ®Id)
Initialises the stocks of products for the avgLive_* years before the simulation starts.
void setTempBool(bool tempBool_h)
vector< vector< double > > vol
double getVHaByYear(const Pixel *px, const int &ft, const int &year, const double &extraBiomass_ratio, const int ®Id) const
return the Volume/ha in a forest after a given number of year after planting, for a specific forest t...
double getMultiplier(const string &multiplierName, const string &forName, int year=DATA_NOW)
vector< string > allProducts
vector< string > secProducts
vector< vector< double > > cumTp_exp
This is the expected version of cumTp, used for calculating profits.
double expTypePrices
Sampling derived expectation types of this agent (prices)
void registerHarvesting(const double &value, const int ®Id, const string &fType)
Registers the harvesting of trees increasing the value of cumEmittedForOper.
Print a debug message, normally filtered out.
Ipopt::SmartPtr< Ipopt::TNLP > OPT
Market optimisation.
Carbon * CBAL
Module for the Carbon Balance.
vector< vector< vector< double > > > hVol_byPrd
vector< vector< double > > area_l
store the areas of the previous year
vector< Pixel * > getAllPlotsByRegion(ModelRegion ®ion_h, bool shuffle=false)
Return the vector of all plots by a specific region (main region or subregion), optionally shuffled;...
string getForTypeParentId(const string &forTypeId_h)
vector< string > getForTypeIds(bool all=false)
By default it doesn't return forTypes used only as input.
double getSd(const vector< K > &v, bool sample=true)
double gpd(const string &type_h, const int ®Id_h, const string &prodId_h, const int &year=DATA_NOW, const string &freeDim_h="") const
map< K, V > vectorToMap(const vector< K > &keys, const V &value=0.0)
Returns a map built using the given vector and the given (scalar) value as keys/values pairs...
vector< int > optFtChosen
vector< double > expectedReturns
V findMap(const map< K, V > &mymap, const K &key, const int &error_level=MSG_CRITICAL_ERROR, const V ¬FoundValue=numeric_limits< V >::min()) const
Lookup a map for a value. Return the value starting from the key.
vector< vector< double > > tp
double getPathMortality(const string &forType, const string &dC, int year=DATA_NOW)
Return the INCREASED mortality due to pathogen presence for a given ft and dc in a certain year (defa...
ModelCoreSpatial(ThreadManager *MTHREAD_h)
double getAvgAgeByDc(Pixel *px, int ft, int dc)
return the average age of a tree at a given diameter class, using the cumTp vector ...
void swap(const int &swap_what)
Assign to the delayed value the current values, e.g. vol_l = vol.
vector< ModelRegion * > getSiblings(bool excludeResidual=true)
Return a vector of pointers to the siblings regions.
void printDetailedHV(map< tr1::array< string, 4 >, double > hVol_byPrd)
vector< double > expectedReturnsNotCorrByRa
by ft. Attenction, reported expReturns at "forest" level (compared with those at forest type level) d...
vector< vector< double > > hVol
void incrOrAddMapValue(map< K, V > &mymap, const K &key, const V &value)
Increments a value stored in a map of the specified value, given the key.
void initialiseEmissionCounters()
Initialises the emission counters to zero.
void setErrorLevel(int errorLevel_h)
void initMarketModule()
computes st and pw for second year and several needed-only-at-t0-vars for the market module ...
vector< int > getRegionIds(int level_h, bool excludeResidual=true)
double computeExpectedPrice(const double &curLocPrice, const double &worldCurPrice, const double &worldFutPrice, const double &sl, const double &sa, const double &expCoef)
Compute weighted expected price for a given product.
void print(bool earlyPrint=true)
Print output. If earlyPrinting it doesn't print some stuff for which we don't yet have values...
void HWP_eol2energy()
Computes the energy substitution for the quota of HWP that reachs end of life and doesn't go to landf...
void initializePixelArea()
compute px->area for each ft and dc
double getDoubleSetting(const string &name_h, int position=0, int reg=WORLD) const
vector< vector< double > > cumTp
This is time of passage to REACH a diameter class (while the exogenous tp by diameter class is the ti...
vector< vector< double > > hArea
void sfd(const double &value_h, const string &type_h, const int ®Id_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW, const bool &allowCreate=false) const
string getRegSName() const
double gfd(const string &type_h, const int ®Id_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW) const
bool app(const string &prod_h, const string &forType_h, const string &dClass_h) const
vector< vector< double > > beta
vector< Pixel * > getMyPixels()
void assignSpMultiplierPropToVols()
ModelCoreSpatial::assignSpMultiplierPropToVols assigns the spatial multiplier (used in the time of re...
ModelRegion * getRegion(int regId_h)
void cacheDynamicSettings()
cache settings that may change with time
void initialiseCarbonModule()
call initialiseDeathBiomassStocks(), initialiseProductsStocks() and initialiseEmissionCounters() ...
void cachePixelExogenousData()
computes pixel level tp, meta and mort
double getAvg(const vector< K > &v)
Returns the average of the elements in the vector.
#define DIAM_ALL
All diameter classes.
void updateMapAreas()
computes forArea_{ft}
void deathTimberInventory_incrOrAdd(const iisskey &thekey, double value_h)
void registerCarbonEvents()
call registerHarvesting(), registerDeathBiomass(), registerProducts() and registerTransports() ...
vector< string > getForTypeChilds(const string &forTypeId_h)
void computeCumulativeData()
computes cumTp_exp, vHa_exp, vHa
void registerProducts(const double &value, const int ®Id, const string &productName)
Registers the production of a given amount of products, storing it in the products maps...
void printDebugInitRegionalValues()
print initial inv, st, sl and sa in each region
void computeEconomicBalances()
compute the policy balances (- -> costs; + -> rev) and the producer/consumer surpluses ...
forType * getForType(int position)
map< int, vector< double > > regArea
const int getMaxYearUsableDeathTimber(const string &prod_h, const string &forType_h, const string &dClass_h)
vector< double > allocateHarvesting(vector< double > total_st, const int ®Id)
Using the deathTimberInventory map, this function allocate the total st in st from death timber (that...