FFSM++  1.1.0
French Forest Sector Model ++
ModelCoreSpatial.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * Copyright (C) 2015 by Laboratoire d'Economie Forestière *
3  * http://ffsm-project.org *
4  * *
5  * This program is free software; you can redistribute it and/or modify *
6  * it under the terms of the GNU General Public License as published by *
7  * the Free Software Foundation; either version 3 of the License, or *
8  * (at your option) any later version, given the compliance with the *
9  * exceptions listed in the file COPYING that is distribued together *
10  * with this file. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program; if not, write to the *
19  * Free Software Foundation, Inc., *
20  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
21  ***************************************************************************/
22 #include <cmath>
23 #include <algorithm>
24 #include <tr1/array>
25 
26 
27 #include "IpIpoptApplication.hpp"
28 #include "IpSolveStatistics.hpp"
29 
30 #include "ModelCoreSpatial.h"
31 #include "ModelData.h"
32 #include "ThreadManager.h"
33 #include "Opt.h"
34 #include "Scheduler.h"
35 #include "Gis.h"
36 #include "Carbon.h"
37 
38 
39 ModelCoreSpatial::ModelCoreSpatial(ThreadManager *MTHREAD_h){
40  MTHREAD = MTHREAD_h;
41 
42 }
43 
44 ModelCoreSpatial::~ModelCoreSpatial(){
45 
46 }
47 
48 void
49 ModelCoreSpatial::runInitPeriod(){
50  cacheSettings(); ///< cashe things like first year, second year, dClasses...
51  cacheDynamicSettings(); // cache settings that may change every year
52  initializePixelVolumes(); ///< compute px volumes vol for 2005 (including exogenous loaded volumes)
53  assignSpMultiplierPropToVols(); // assign the spatial multiplier (used in the time of return) based no more on a Normal distribution but on the volumes present in the pixel: more volume, more the pixel is fit for the ft
54  initMarketModule(); ///< inside it uses first year, second year
55  initialiseDeathTimber();
56  MTHREAD->DO->print(true);
57  MTHREAD->SCD->advanceYear(); ///< 2005->2006
58  int thisYear = MTHREAD->SCD->getYear(); // for debugging
59  msgOut(MSG_INFO, "### 2013 "+i2s(thisYear)+" year started..");
60  resetPixelValues(); ///< swap volumes->lagged_volumes and reset the other pixel vectors
61  cachePixelExogenousData(); ///< compute pixel tp, meta and mort
62  computeInventary(); ///< in=f(vol_t-1)
63  //printDebugInitRegionalValues();
64  computeCumulativeData(); ///< compute cumTp_exp, vHa_exp, vHa
65  initializePixelArea(); ///< compute px->area for each ft and dc (including exogenous loaded areas)
66  runBiologicalModule();
67  runManagementModule();
68  updateMapAreas(); ///< update the forArea_{ft} layer on each pixel as old value-hArea+regArea
69  updateOtherMapData(); ///< update (if the layer exists) other gis-based data, as volumes and expected returns, taking them from the data in the px object
70  sumRegionalForData(); ///< only for printing stats as forest data is never used at regional level
71  initialiseCarbonModule();
72  //MTHREAD->DO->printDebugPixelValues(); // uncomment to enable pixel-level debugging // moved to output->print
73 
74  MTHREAD->DO->print();
75 }
76 
77 void
78 ModelCoreSpatial::runSimulationYear(){
79  int thisYear = MTHREAD->SCD->getYear(); // for debugging
80  cacheDynamicSettings(); // cache settings that may change every year
81  resetPixelValues(); // swap volumes->lagged_volumes and reset the other pixel vectors
82  cachePixelExogenousData(); // compute pixel tp, meta and mort
83  computeInventary(); // in=f(vol_t-1)
84  runMarketModule(); // RUN THE MARKET OPTIMISATION HERE
85  computeCumulativeData(); // compute cumTp_exp, vHa_exp
86  cachePixelExogenousData();
87  runBiologicalModule();
88  runManagementModule();
89  //MTHREAD->DO->printDebugPixelValues(); // uncomment to enable pixel-level debugging // moved to output->print
90  updateMapAreas();
91  updateOtherMapData(); // update (if the layer exists) other gis-based data, as volumes and expected returns, taking them from the data in the px object
92  sumRegionalForData(); // only for printing stats as forest data is never used at regional level
93  registerCarbonEvents();
94  computeEconomicBalances(); // policy costs and welfare analysis
95  MTHREAD->DO->print(false);
96 }
97 
98 void
99 ModelCoreSpatial::initMarketModule(){
100  msgOut(MSG_INFO, "Starting market module (init stage)..");
101 
102  for(uint i=0;i<regIds2.size();i++){
103  int r2 = regIds2[i];
104  //RPAR('pl',i,p_tr,t-1) = sum(p_pr, a(p_pr,p_tr)*RPAR('pl',i,p_pr,t-1))+m(i,p_tr);
105  for(uint sp=0;sp<secProducts.size();sp++){
106  double value = 0;
107  for (uint pp=0;pp<priProducts.size();pp++){
108  value += gpd("pl",r2,priProducts[pp],secondYear)*
109  gpd("a",r2,priProducts[pp],secondYear,secProducts[sp]);
110  }
111  value += (gpd("m",r2,secProducts[sp],secondYear)* gpd("pol_trSub",r2,secProducts[sp],secondYear));
112  spd(value,"pl",r2,secProducts[sp],secondYear,true);
113  }
114  // RPAR('dl',i,p_pr,t-1) = sum(p_tr, a(p_pr,p_tr)*RPAR('sl',i,p_tr,t-1));
115  for (uint pp=0;pp<priProducts.size();pp++){
116  double value=0;
117  for(uint sp=0;sp<secProducts.size();sp++){
118  value += gpd("sl",r2,secProducts[sp],secondYear)*
119  gpd("a",r2,priProducts[pp],secondYear, secProducts[sp]);
120  }
121  spd(value,"dl",r2,priProducts[pp],secondYear,true);
122  double stvalue = gpd("sl",r2,priProducts[pp],secondYear)
123  + gpd("sa",r2,priProducts[pp],secondYear);
124  spd(stvalue,"st",r2,priProducts[pp],secondYear,true);
125  double supply_fromForestShare = gpd("supply_fromForestShare",r2,priProducts[pp],secondYear);
126  spd(stvalue*supply_fromForestShare,"st_fromHarv",r2,priProducts[pp],secondYear,true); // needed for the biological module in year starty+1
127  }
128  // RPAR('st',i,prd,t-1) = RPAR('sl',i,prd,t-1)+RPAR('sa',i,prd,t-1);
129  // RPAR('dt',i,prd,t-1) = RPAR('dl',i,prd,t-1)+RPAR('da',i,prd,t-1);
130  for (uint ap=0;ap<allProducts.size();ap++){
131  //double debug = gpd("dl",r2,allProducts[ap],secondYear);
132  double dtvalue = gpd("dl",r2,allProducts[ap],secondYear)
133  + gpd("da",r2,allProducts[ap],secondYear);
134 
135  spd(dtvalue,"dt",r2,allProducts[ap],secondYear,true);
136  }
137 
138  // q1(i,p_tr) = 1/(1+((RPAR('dl',i,p_tr,t-1)/RPAR('da',i,p_tr,t-1))**(1/psi(i,p_tr)))*(RPAR('pl',i,p_tr,t-1)/PT(p_tr,t-1)));
139  // p1(i,p_tr) = 1-q1(i,p_tr);
140  // RPAR('dc',i,p_tr,t-1) = (q1(i,p_tr)*RPAR('da',i,p_tr,t-1)**((psi(i,p_tr)-1)/psi(i,p_tr))+ p1(i,p_tr)*RPAR('dl',i,p_tr,t-1)**((psi(i,p_tr)-1)/psi(i,p_tr)))**(psi(i,p_tr)/(psi(i,p_tr)-1));
141  // RPAR('pc',i,p_tr,t-1) = (RPAR('da',i,p_tr,t-1)/RPAR('dc',i,p_tr,t-1))*PT(p_tr,t-1)+(RPAR('dl',i,p_tr,t-1)/RPAR('dc',i,p_tr,t-1))*RPAR('pl',i,p_tr,t-1);
142  // RPAR('pc',i,p_pr,t-1) = (RPAR('sa',i,p_pr,t-1)/RPAR('sc',i,p_pr,t-1))*PT(p_pr,t-1)+(RPAR('sl',i,p_pr,t-1)/RPAR('sc',i,p_pr,t-1))*RPAR('pl',i,p_pr,t-1);
143  // RPAR('pw',i,p_tr,t-1) = (RPAR('dl',i,p_tr,t-1)*RPAR('pl',i,p_tr,t-1)+RPAR('da',i,p_tr,t-1)*PT(p_tr,t-1))/RPAR('dt',i,p_tr,t-1) ; //changed 20120419
144  // K(i,p_tr,t-1) = k1(i,p_tr)*RPAR('sl',i,p_tr,t-1);
145  for(uint sp=0;sp<secProducts.size();sp++){
146  double psi = gpd("psi",r2,secProducts[sp],secondYear);
147  double dl = gpd("dl",r2,secProducts[sp],secondYear);
148  double da = gpd("da",r2,secProducts[sp],secondYear);
149  double pl = gpd("pl",r2,secProducts[sp],secondYear);
150  double sa = gpd("sa",r2,secProducts[sp],secondYear);
151  double sl = gpd("sl",r2,secProducts[sp],secondYear);
152  double k1 = gpd("k1",r2,secProducts[sp],secondYear);
153  double pWo = gpd("pl",WL2,secProducts[sp],secondYear); // World price (local price for region 99999)
154 
155  double stvalue = sl+sa;
156  double q1 = 1/ ( 1+pow(dl/da,1/psi)*(pl/pWo) );
157  double p1 = 1-q1;
158  double dc = pow(
159  q1*pow(da,(psi-1)/psi) + p1*pow(dl,(psi-1)/psi)
160  ,
161  psi/(psi-1)
162  );
163  double pc = (da/dc)*pWo
164  +(dl/dc)*pl;
165  double pw = (dl*pl+da*pWo)/(dl+da);
166  double k = k1*sl;
167 
168  spd(stvalue,"st",r2,secProducts[sp],secondYear,true);
169  spd(q1,"q1",r2,secProducts[sp],firstYear,true);
170  //spd(p1,"p1",r2,secProducts[sp],firstYear,true);
171  spd(dc,"dc",r2,secProducts[sp],secondYear,true);
172  spd(pc,"pc",r2,secProducts[sp],secondYear,true);
173  spd(pw,"pw",r2,secProducts[sp],secondYear,true);
174  spd(k,"k",r2,secProducts[sp],secondYear,true);
175  }
176 
177  // t1(i,p_pr) = 1/(1+((RPAR('sl',i,p_pr,t-1)/RPAR('sa',i,p_pr,t-1))**(1/eta(i,p_pr)))*(RPAR('pl',i,p_pr,t-1)/PT(p_pr,t-1)));
178  // r1(i,p_pr) = 1-t1(i,p_pr);
179  // RPAR('sc',i,p_pr,t-1) = (t1(i,p_pr)*RPAR('sa',i,p_pr,t-1)**((eta(i,p_pr)-1)/eta(i,p_pr))+ r1(i,p_pr)*RPAR('sl',i,p_pr,t-1)**((eta(i,p_pr)-1)/eta(i,p_pr)))**(eta(i,p_pr)/(eta(i,p_pr)-1))
180  // RPAR('pc',i,p_pr,t-1) = (RPAR('sa',i,p_pr,t-1)/RPAR('sc',i,p_pr,t-1))*PT(p_pr,t-1)+(RPAR('sl',i,p_pr,t-1)/RPAR('sc',i,p_pr,t-1))*RPAR('pl',i,p_pr,t-1);
181  // RPAR('pw',i,p_pr,t-1) = (RPAR('sl',i,p_pr,t-1)*RPAR('pl',i,p_pr,t-1)+RPAR('sa',i,p_pr,t-1)*PT(p_pr,t-1))/RPAR('st',i,p_pr,t-1) ; //changed 20120419
182  for(uint pp=0;pp<priProducts.size();pp++){
183 
184  double sl = gpd("sl",r2,priProducts[pp],secondYear);
185  double sa = gpd("sa",r2,priProducts[pp],secondYear);
186  double eta = gpd("eta",r2,priProducts[pp],secondYear);
187  double pl = gpd("pl",r2,priProducts[pp],secondYear);
188  double pWo = gpd("pl",WL2,priProducts[pp],secondYear); // World price (local price for region 99999)
189 
190 
191  double t1 = 1/ ( 1+(pow(sl/sa,1/eta))*(pl/pWo) );
192  double r1 = 1-t1;
193  double sc = pow(
194  t1*pow(sa,(eta-1)/eta) + r1*pow(sl,(eta-1)/eta)
195  ,
196  eta/(eta-1)
197  );
198  double pc = (sa/sc)*pWo+(sl/sc)*pl;
199  double pw = (sl*pl+sa*pWo)/(sl+sa);
200 
201  spd(t1,"t1",r2,priProducts[pp],firstYear,true);
202  //spd(r1,"r1",r2,priProducts[pp],firstYear,true);
203  spd(sc,"sc",r2,priProducts[pp],secondYear,true);
204  spd(pc,"pc",r2,priProducts[pp],secondYear,true);
205  spd(pw,"pw",r2,priProducts[pp],secondYear,true);
206  }
207 
208  // up to here tested with gams output on 20120628, that's fine !!
209  } // end for each region in level 2
210 
211 
212  // initializing the exports to zero quantities
213  // initializing of the transport cost for the same region to one and distance to zero
214  for(uint r1=0;r1<l2r.size();r1++){
215  for(uint r2=0;r2<l2r[r1].size();r2++){
216  for(uint p=0;p<allProducts.size();p++){
217  for(uint r2To=0;r2To<l2r[r1].size();r2To++){
218  spd(0,"rt",l2r[r1][r2],allProducts[p],secondYear,true,i2s(l2r[r1][r2To])); // regional trade, it was exp in gams
219  if(l2r[r1][r2] == l2r[r1][r2To]){
220  spd(1,"ct",l2r[r1][r2],allProducts[p],firstYear,true,i2s(l2r[r1][r2To])); // as long this value is higher than zero, rt within the same region is not choosen by the solver, so the value doesn't really matters. If it is zero, the solver still works and results are the same, but reported rt within the region are crazy high (100000)
221  }
222  }
223  } // end each product
224 
225  for(uint r2To=0;r2To<l2r[r1].size();r2To++){
226  if(l2r[r1][r2] == l2r[r1][r2To]){
227  spd(0,"dist",l2r[r1][r2],"",firstYear,true,i2s(l2r[r1][r2To])); // setting distance zero in code, so no need to put it in the data
228  }
229  }
230  } // end of r2 regions
231  } // end of r1 region
232 }
233 
234 void
235 ModelCoreSpatial::runMarketModule(){
236  msgOut(MSG_INFO, "Starting market module");
237  static double cumOverHarvesting = 0.0;
238  int thisYear = MTHREAD->SCD->getYear();
239  int previousYear = MTHREAD->SCD->getYear()-1;
240  bool useDeathTimber = MTHREAD->MD->getBoolSetting("useDeathTimber",DATA_NOW);
241 
242  // Added for carbon poject -------------------------------
243  double intRate = MTHREAD->MD->getDoubleSetting("ir");
244 
245  // ----------------------------------------------------------------
246 
247  // *** PRE-OPTIMISATION YEARLY OPERATIONS..
248  for(uint i=0;i<regIds2.size();i++){
249  int r2 = regIds2[i];
250  for(uint sp=0;sp<secProducts.size();sp++){
251  double g1 = gpd("g1",r2,secProducts[sp],previousYear);
252  double sigma = gpd("sigma",r2,secProducts[sp]);
253  double pc_1 = gpd("pc",r2,secProducts[sp],previousYear);
254  double dc_1 = gpd("dc",r2,secProducts[sp],previousYear);
255  double k_1 = gpd("k",r2,secProducts[sp],previousYear);
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);
258  double q1 = gpd("q1",r2,secProducts[sp]);
259  double pol_mktStr_d = gpd("pol_mktStr_d",r2,secProducts[sp]);
260 
261  double q1_polCorr = min(1.0, (q1-0.5)*pol_mktStr_d+0.5); // The pol_mktStr_s coefficient makes t1 (the "b" weight in the CES function) converge toward 0.5 as pol_mktStr_s goes to 0. As pol_mktStr_s can be above 1, an upper cap of unity is imposed to t1
262 
263 
264 
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); //alpha
268 
269  spd(k, "k" ,r2,secProducts[sp]);
270  spd(aa,"aa",r2,secProducts[sp],DATA_NOW,true);
271  spd(gg,"gg",r2,secProducts[sp],DATA_NOW,true);
272  spd(q1_polCorr,"q1_polCorr",r2,secProducts[sp],DATA_NOW,true);
273 
274  }
275 
276  // BB(i,p_pr) = (sigma(p_pr)/(sigma(p_pr)+1))*RPAR('pc',i,p_pr,t-1)*(RPAR('sc',i,p_pr,t-1)**(-1/sigma(p_pr)))*(In(i,p_pr,t-1)/In(i,p_pr,t))**(gamma(p_pr)/sigma(p_pr));
277  // FF(i,p_pr) = RPAR('sc',i,p_pr,t-1)*((RPAR('pc',i,p_pr,t-1))**(-sigma(p_pr)))*(In(i,p_pr,t)/In(i,p_pr,t-1))**(gamma(p_pr)); //chi
278  for(uint pp=0;pp<priProducts.size();pp++){
279  double gamma_incr = gpd("gamma_incr",r2,priProducts[pp]); // elast supply to increasing stocks
280  double gamma_decr = gpd("gamma_decr",r2,priProducts[pp]); // elast supply to decreasing stocks
281  double gamma_d_incr = gpd("gamma_d_incr",r2,priProducts[pp]); // elast supply to increasing stocks death timber
282  double gamma_d_decr = gpd("gamma_d_decr",r2,priProducts[pp]); // elast supply to decreasing stocks death timber
283 
284  double sigma = gpd("sigma",r2,priProducts[pp]); // elast supply to price
285  double pc_1 = gpd("pc",r2,priProducts[pp],previousYear);
286  double sc_1 = gpd("sc",r2,priProducts[pp],previousYear);
287  double in = gpd("in",r2,priProducts[pp]);
288  double in_0 = gpd("in",r2,priProducts[pp],MTHREAD->MD->getIntSetting("initialYear"));
289  double in_1 = gpd("in",r2,priProducts[pp],previousYear);
290  double in_d =gpd("in_deathTimber",r2,priProducts[pp]);
291  double in_d_1 ;
292  if (thisYear == MTHREAD->MD->getIntSetting("initialOptYear")){
293  in_d_1 = gpd("in_deathTimber",r2,priProducts[pp],thisYear);
294  } else {
295  in_d_1 = gpd("in_deathTimber",r2,priProducts[pp],previousYear);
296  }
297  double min_in_d = max(0.00001,in_0 * MTHREAD->MD->getDoubleSetting("minShareDeathIn"));
298 
299 
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);
302  double t1 = gpd("t1",r2,priProducts[pp]);
303  double pol_mktStr_s = gpd("pol_mktStr_s",r2,priProducts[pp]);
304 
305 
306  // Added for carbon project ------------------------------------------
307  double carbonPrice_1 = gpd("carbonPrice",r2,"",previousYear); //using dummy region until Anto solves issue // Antonello: solved but also generalised
308  double omega = gpd("co2content_products", r2, priProducts[pp])/1000; //kg to tons
309  double Pct_1 = carbonPrice_1*intRate;
310 
311  // -------------------------------------------------------------------
312 
313 
314  double t1_polCorr = min(1.0, (t1-0.5)*pol_mktStr_s+0.5); // The pol_mktStr_s coefficient makes t1 (the "b" weight in the CES function) converge toward 0.5 as pol_mktStr_s goes to 0. As pol_mktStr_s can be above 1, an upper cap of unity is imposed to t1
315 
316 
317  double gamma = in>in_1 ? gamma_incr: gamma_decr; // If inventories resources are decreasing we use a given elasticity, if they are increasing we use an other one
318  double gamma_d = in_d>in_d_1 ? gamma_d_incr: gamma_d_decr; // If inventories resources are decreasing we use a given elasticity, if they are increasing we use an other one
319 
320 
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 ; // to avoid depending on the ratio of too small variations in mortality
323 
324  //msgOut(MSG_INFO,"in_d: "+d2s(in_d));
325  //msgOut(MSG_INFO,"in_d_1: "+d2s(in_d_1));
326  //msgOut(MSG_INFO,"in: "+d2s(in));
327  //msgOut(MSG_INFO,"in_1: "+d2s(in_1));
328  //msgOut(MSG_INFO,"in_d_ratio: "+d2s(in_d_ratio));
329 
330 #ifdef QT_DEBUG
331 if (in_d_ratio < 0.01){
332  msgOut(MSG_CRITICAL_ERROR,"Very low in_d_ratio ");
333 }
334 if (in_d_ratio > 100){
335  msgOut(MSG_CRITICAL_ERROR,"Very high in_d_ratio ");
336 }
337 #endif
338 
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 ;
340 
341  /* Old FF (without carbon)
342  double ff = sc_1*pow(pc_1*pol_mktDirInt_s_1,-sigma)*pow(in_ratio,gamma); //chi
343  */
344 
345  // Added for carbon project
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);
347  // -----------------------------------------------
348 
349  spd(bb,"bb",r2,priProducts[pp],DATA_NOW,true);
350  spd(ff,"ff",r2,priProducts[pp],DATA_NOW,true);
351  spd(sigma,"sigma",r2,priProducts[pp],DATA_NOW,true);
352  spd(t1_polCorr,"t1_polCorr",r2,priProducts[pp],DATA_NOW,true);
353 
354  }
355 
356  } // end for each region in level 2 (and updating variables)
357 
358  //cout << "### "+i2s(thisYear-2)+"abb: " << d2s( MD->deathTimberInventory_get(iisskey (thisYear-2,11001,"Fut_Feu","30") ) ) << endl;
359 
360  // *** OPTIMISATION....
361 
362  // Create an instance of the IpoptApplication
363  //Opt *OPTa = new Opt(MTHREAD);
364  //SmartPtr<TNLP> OPTa = new Opt(MTHREAD);
365  SmartPtr<IpoptApplication> application = new IpoptApplication();
366  string linearSolver = MTHREAD->MD->getStringSetting("linearSolver",DATA_NOW);
367  application->Options()->SetStringValue("linear_solver", linearSolver); // default in ipopt is ma27
368  //application->Options()->SetStringValue("hessian_approximation", "limited-memory"); // quasi-newton approximation of the hessian
369  //application->Options()->SetIntegerValue("mumps_mem_percent", 100);
370  application->Options()->SetNumericValue("obj_scaling_factor", -1); // maximisation
371  application->Options()->SetNumericValue("max_cpu_time", 1800); // max 1/2 hour to find the optimus for one single year
372  application->Options()->SetStringValue("check_derivatives_for_naninf", "yes");
373  //application->Options()->SetStringValue("halt_on_ampl_error", "yes"); // not a valid option
374  //application->Options()->SetStringValue("print_options_documentation", "yes"); // print the ipopt options (hundreds of pages!)
375  // Relaxing tollerance options... worster effect :-(
376  //application->Options()->SetNumericValue("tol", 1e-07);
377  //application->Options()->SetNumericValue("constr_viol_tol",0.001);
378  //application->Options()->SetNumericValue("compl_inf_tol",0.001);
379  //application->Options()->SetNumericValue("acceptable_tol", 1e-05);
380  //application->Options()->SetNumericValue("acceptable_dual_inf_tol", 1e+9);
381  //application->Options()->SetNumericValue("acceptable_constr_viol_tol", 0.1);
382  //application->Options()->SetNumericValue("acceptable_compl_inf_tol", 0.1);
383 
384 
385  // Initialize the IpoptApplication and process the options
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?");
391  return;
392  }
393 
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);
396 
397  // cout << "### "+i2s(thisYear-2)+"abc: " << d2s( MD->deathTimberInventory_get(iisskey (thisYear-2,11001,"Fut_Feu","30") ) ) << endl;
398  // *** POST OPTIMISATION....
399 
400  // post-equilibrium variables->parameters assignments..
401  // RPAR(type,i,prd,t) = RVAR.l(type,i,prd);
402  // EX(i,j,prd,t) = EXP.l(i,j,prd);
403  // ObjT(t) = Obj.l ;
404  // ==> in Opt::finalize_solution()
405 
406  // Retrieve some statistics about the solve
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;
415  MTHREAD->DO->printOptLog(true, icount, obj);
416  } else {
417  //Number final_obj = application->Statistics()->FinalObjective();
418  int thisYear = MTHREAD->SCD->getYear();
419  cout << "***ERROR: MODEL DIDN'T SOLVE FOR THIS YEAR ("<<thisYear<<")"<< endl;
420  msgOut(MSG_CRITICAL_ERROR, "Model DIDN'T SOLVE for this year ("+i2s(thisYear)+")");
421  // IMPORTANT! Don't place the next two lines above the msgOut() function or it will crash in windows if the user press the stop button
422  //Index iter_count = application->Statistics()->IterationCount(); // syserror if model doesn't solve
423  //Number final_obj = application->Statistics()->FinalObjective();
424  int icount = 0;
425  double obj = 0;
426  MTHREAD->DO->printOptLog(false, icount, obj);
427  }
428 
429  for(uint r2= 0; r2<regIds2.size();r2++){ // you can use r2<=regIds2.size() to try an out-of range memory error that is not detected other than by valgrind (with a message "Invalid read of size 4 in ModelCore::runSimulationYear() in src/ModelCore.cpp:351")
430  int regId = regIds2[r2];
431  ModelRegion* REG = MTHREAD->MD->getRegion(regId);
432 
433  // *** st adjustments:***
434  // st_or[pp] -> original st from the mkt module
435  // st_fromFor[pp] -> cover the volumes from both alive trees and death ones. This should be bounded by inResByAnyCombination, that already includes death resources
436  // st_fromHarv[pp] -> from harvesting. This should be equal to hV and IGN data
437  // st
438  // st_or = st
439  // st_fromFor = st_or * a
440  // st_fromFor > in:
441  // (overharvesting)
442  // overHarv = st_fromFor - in
443  // st_fromFor = in
444  // st = st_fromFor/a
445  // st_fromHarv = distribute(st_fromFor)
446 
447  // // total supply and total demand..
448  // RPAR('st',i,prd,t) = RPAR('sl',i,prd,t)+RPAR('sa',i,prd,t);
449  // RPAR('dt',i,prd,t) = RPAR('dl',i,prd,t)+RPAR('da',i,prd,t);
450  // // weighted prices.. //changed 20120419
451  // RPAR('pw',i,p_tr,t) = (RPAR('dl',i,p_tr,t)*RPAR('pl',i,p_tr,t)+RPAR('da',i,p_tr,t)*PT(p_tr,t))/RPAR('dt',i,p_tr,t) ; //changed 20120419
452  // RPAR('pw',i,p_pr,t) = (RPAR('sl',i,p_pr,t)*RPAR('pl',i,p_pr,t)+RPAR('sa',i,p_pr,t)*PT(p_pr,t))/RPAR('st',i,p_pr,t) ; //changed 20120419
453  for (uint p=0;p<allProducts.size();p++){
454  double st = gpd("sl",regId,allProducts[p])+gpd("sa",regId,allProducts[p]);
455  double dt = gpd("dl",regId,allProducts[p])+gpd("da",regId,allProducts[p]);
456  spd(st,"st",regId,allProducts[p]);
457  spd(st,"st_or",regId,allProducts[p],DATA_NOW,true); // original total supply, not corrected by resetting it to min(st, inv).
458  spd(dt,"dt",regId,allProducts[p]);
459  }
460  for (uint p=0;p<secProducts.size();p++){
461  double dl = gpd("dl",regId,secProducts[p]);
462  double pl = gpd("pl",regId,secProducts[p]);
463  double da = gpd("da",regId,secProducts[p]); // bug corrected 20120913
464  double pworld = gpd("pl", WL2,secProducts[p]);
465  double dt = gpd("dt",regId,secProducts[p]);
466  double pw = dt?(dl*pl+da*pworld)/dt:0.0;
467  spd(pw,"pw",regId,secProducts[p]);
468  }
469  for (uint p=0;p<priProducts.size();p++){
470  double sl = gpd("sl",regId,priProducts[p]);
471  double pl = gpd("pl",regId,priProducts[p]);
472  double sa = gpd("sa",regId,priProducts[p]); // bug corrected 20120913
473  double pworld = gpd("pl", WL2,priProducts[p]);
474  double st = gpd("st",regId,priProducts[p]);
475  double st_fromFor = st * gpd("supply_fromForestShare",regId,priProducts[p]);
476  double pw = st?(sl*pl+sa*pworld)/st:0.0;
477  spd(pw,"pw",regId,priProducts[p]);
478  spd(st_fromFor,"st_fromFor",regId,priProducts[p],DATA_NOW,true);
479  }
480 
481  // Correcting st if this is over the in
482 
483  // Create a vector with all possible combinations of primary products
484  vector<vector<int>> priPrCombs = MTHREAD->MD->createCombinationsVector(priProducts.size());
485  int nPriPrCombs = priPrCombs.size();
486 
487  for (uint i=0;i<priPrCombs.size();i++){
488  double stMkMod = 0.0;
489  double sumIn = REG->inResByAnyCombination[i];
490  // double sumIn2 = 0.0;
491  for (uint p=0;p<priPrCombs[i].size();p++){
492  stMkMod += gpd("st_fromFor",regId,priProducts[priPrCombs[i][p]]);
493  //sumIn2 += gpd("in",regId,priProducts[priPrCombs[i][p]]);
494  }
495 
496  //if(sumIn<=0.00001){
497  // for (uint p=0;p<priPrCombs[i].size();p++){
498  // spd(0.0,"st",regId,priProducts[priPrCombs[i][p]]);
499  // }
500  // } else {
501  if(stMkMod>sumIn){ // if we harvested more than available
502  string pProductsInvolved = "";
503  for (uint p=0;p<priPrCombs[i].size();p++){
504  pProductsInvolved += (priProducts[priPrCombs[i][p]]+"; ");
505  }
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]]);
512  }
513  }
514 
515  //}
516 
517 
518  }
519 
520  //cout << "### "+i2s(thisYear-2)+"abe: " << d2s( MD->deathTimberInventory_get(iisskey (thisYear-2,11001,"Fut_Feu","30") ) ) << endl;
521  // here we create st_fromHarv as st_fromFor - st_from_deathbiomass
522  vector <double> total_st(priProducts.size(),0.);
523  vector <double> in_deathTimber(priProducts.size(),0.);
524  vector <double> in_aliveForest (priProducts.size(),0.);
525  for (uint i=0;i<priProducts.size();i++){
526  total_st[i] = gpd("st_fromFor",regId,priProducts[i]);
527  in_deathTimber[i] = gpd("in_deathTimber",regId,priProducts[i]);
528  in_aliveForest[i] = gpd("in",regId,priProducts[i]);
529  }
530  // cout << "### "+i2s(thisYear-2)+"abf: " << d2s( MD->deathTimberInventory_get(iisskey (thisYear-2,11001,"Fut_Feu","30") ) ) << endl;
531  vector <double> st_fromHarv= allocateHarvesting(total_st, regId);
532  //cout << "### "+i2s(thisYear-2)+"abg: " << d2s( MD->deathTimberInventory_get(iisskey (thisYear-2,11001,"Fut_Feu","30") ) ) << endl;
533  for (uint i=0;i<priProducts.size();i++){
534  spd(st_fromHarv[i],"st_fromHarv",regId,priProducts[i],DATA_NOW,true);
535  }
536 
537  } // end of each region
538 
539  // cout << "### "+i2s(thisYear-2)+"abh: " << d2s( MD->deathTimberInventory_get(iisskey (thisYear-2,11001,"Fut_Feu","30") ) ) << endl;
540 
541  if (cumOverHarvesting>0.0){
542  msgOut(MSG_DEBUG, "Overharvesting is present. Year: "+i2s(thisYear)+" cumOverHarvesting: "+d2s(cumOverHarvesting));
543  }
544 }
545 
546 /**
547  * @brief ModelCoreSpatial::runBiologicalModule
548  *
549  * Changes in Area:
550  * dc area_l area diff
551  * 0 ---------> +regArea -areaFirstProdClass (areaMovingUp_00)
552  * 15 ---------> +areaFirstPrClass -hArea_15 -areaMovingUp_15
553  * 25 ---------> +areaMovingUp15 - hArea_25 - areaMovingUp_25
554  * 35 ---------> +areaMovingUp25 - hArea_35 - areaMovingUp_35
555  * ...
556  * 95 ---------> +areaMovingUp85 - hArea_95 - areaMovingUp_95
557  * 105 ---------> +areaMovingUp95 - hArea_105
558  *
559  * note: regArea is computed in the management module, not here. Further, regArea is already the net one of forest area changes
560  */
561 void
562 ModelCoreSpatial::runBiologicalModule(){
563 
564  msgOut(MSG_INFO, "Starting resource module..");
565  int thisYear = MTHREAD->SCD->getYear();
566  bool useDeathTimber = MD->getBoolSetting("useDeathTimber",DATA_NOW);
567 
568  for(uint i=0;i<regIds2.size();i++){
569  int r2 = regIds2[i];
570  int regId = r2;
571  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
572  //Gis* GIS = MTHREAD->GIS;
573  regPx = REG->getMyPixels();
574  double shareMortalityUsableTimber = 0.0;
575 
576  for (uint p=0;p<regPx.size();p++){
577  Pixel* px = regPx[p];
578 
579  double pxId = px->getID();
580  //if (pxId == 3550.0){
581  // cout << "got the pixel" << endl;
582  //}
583  //px->expectedReturns.clear();
584  for(uint j=0;j<fTypes.size();j++){
585  string ft = fTypes[j];
586  //double pxArea_debug = px->getDoubleValue("forArea_"+ft, true);
587  vector <double> hV_byDiam;
588  vector <double> productivity_byDc; // mc/ha/y
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 ; // fire mortality
594  vector <double> areasMovingUp(dClasses.size(), 0.0);
595  double areaFirstProdClass;
596 
597 
598  // A - COMPUTING THE REGENERATION..
599  // if we are in a year where the time of passage has not yet been reached
600  // for the specific i,e,l then we use the exogenous Vregen, otherwise we
601  // calculate it
602  //if ( not scen("fxVreg") ,
603  // loop( (i,essence,lambda),
604  // if( ord(t)>=(tp_u1(i,essence,lambda)+2),
605  // Vregen(i,lambda,essence,t)=regArea(i,essence,lambda,t-tp_u1(i,essence,lambda))*volHa_u1(i,essence,lambda)/1000000 ;
606  // );
607  // );
608  //);
609  int tp_u0 = px->tp.at(j).at(0); // time of passage to reach the first production diameter class // bug 20140318, added ceil. 20140318 removed it.. model did go crazy with it
610  if(thisYear == secondYear){
611  px->initialDc0Area.push_back(px->area_l.at(j).at(0));
612  }
613  if(regType != "fixed" && (thisYear-secondYear) >= tp_u0 ) { // T.O.D.O to be checked -> 20121109 OK
614  double pastRegArea = px->getPastRegArea(j,thisYear-tp_u0);
615  double availableArea = px->area_l.at(j).at(0);
616  //double entryVolHa = gfd("entryVolHa",regId,ft,"");
617  double vHa = px->vHa.at(j).at(1);
618  //attenction that at times could take the wrong pastRegArea if tp change too suddenly as in some "strange" scenarios
619  if (oldVol2AreaMethod){
620  areaFirstProdClass = pastRegArea;
621  } else {
622  areaFirstProdClass = min(availableArea, pastRegArea); // this is just a start and will need to include the last year area
623  }
624  px->vReg.push_back(areaFirstProdClass*vHa/1000000.0); // TO.DO: check the 1000000. Should be ok, as area in ha vol in Mm^3
625  //if (pxId == 3550.0 && j==3){
626  // cout << "got the pixel" << endl;
627  //}
628  #ifdef QT_DEBUG
629  if (areaFirstProdClass < 0.0){
630  //msgOut(MSG_CRITICAL_ERROR,"Negative regeneration volumes in endogenous regeneration");
631  }
632  if ( (availableArea-pastRegArea) < -0.00000001 ){
633  // in a very rare cases tp change first in a direction and then in the other, so that the wrong past regeneration area
634  // is picken up.
635  //msgOut(MSG_CRITICAL_ERROR,"Upgrading from dc0 more area than the available one in endogenous regeneration");
636  }
637  #endif
638  } else {
639  double regionArea = REG->getValue("forArea_"+ft,OP_SUM);
640  double pxArea = px->getDoubleValue("forArea_"+ft, true); // 20121109 bug solved (add get zero for not data)
641  double regRegVolumes = gfd("vReg",r2,ft,"");
642  double newVReg = regionArea ? regRegVolumes*pxArea/regionArea : 0.0;
643  px->vReg.push_back(newVReg); // 20121108 BUG !!! solved // as now we have the area we could also use here entryVolHa
644  // only a share of the exogenous area goes up, the regeneration one doesn't yet reach tp0:
645  // areaFirstProdClass = (1.0 / px->tp.at(j).at(0) ) * px->area_l.at(j).at(0);
646  areaFirstProdClass = (1.0 / ((double) tp_u0) ) * px->initialDc0Area.at(j);
647  // in the exogenous period we are exogenously upgrading u0->u1 some areas but, as we do not have the regeneration
648  // are corresponding to that we have also to manually add it to u0
649  //px->area_l.at(j).at(0) += areaFirstProdClass;
650  //areaFirstProdClass = entryVolHa ? newVReg*1000000 /entryVolHa:0.0;
651  //if (pxId == 3550.0 && j==3){
652  // cout << "got the pixel" << endl;
653  //}
654 
655  #ifdef QT_DEBUG
656  if (areaFirstProdClass<0.0){
657  // msgOut(MSG_CRITICAL_ERROR,"Negative regeneration volumes in exogenous regeneration");
658  }
659  if (areaFirstProdClass > px->area_l.at(j).at(0)){
660  //msgOut(MSG_CRITICAL_ERROR,"Moving up area higher than available area in exogenous regeneration !");
661  }
662  #endif
663  // vReg and entryVolHa are NOT the same thing. vReg is the yearly regeneration volumes
664  // for the whole region. We can use them when we don't know the harvested area
665  // entryVolHa can lead to vReg calculation only when we know the regeneration area. So in the
666  // first years we use vReg and subsequently the endogenous one.
667  }
668 
669  //double harvestedArea = 0;
670 
671  if(useDeathTimber){
672  shareMortalityUsableTimber = gfd("shareMortalityUsableTimber",r2,ft,"");
673  } else {
674  shareMortalityUsableTimber = 0.0;
675  }
676 
677  for (uint u=0; u<dClasses.size(); u++){
678  string dc = dClasses[u];
679  double hr = 0;
680  double age = getAvgAgeByDc(px, j, u);
681 
682  //double pastYearVol_reg = u ? gfd("vol",r2,ft,dc,thisYear-1): 0;
683  double pastYearVol = px->vol_l.at(j).at(u);
684  vector <double> hV_byPrd;
685  vector <double> hr_byPrd;
686 
687  // harvesting rate & volumes...
688  // hr is by region.. no reasons in one pixel the RATE of harvesting will be different than in an other pixel
689  //hr(u,i,essence,lambda,t) = sum(p_pr, prov(u,essence,lambda,p_pr)*RPAR('st',i,p_pr,t)/In(i,p_pr,t));
690  //hV(u,i,essence,lambda,t) = hr(u,i,essence,lambda,t) * V(u,i,lambda,essence,t-1);
691  //hV_byPrd(u,i,essence,lambda,p_pr,t) = prov(u,essence,lambda,p_pr)*(RPAR('st',i,p_pr,t)/In(i,p_pr,t))*V(u,i,lambda,essence,t-1);
692  for(uint pp=0;pp<priProducts.size();pp++){
693  double st = gpd("st_fromHarv",r2,priProducts[pp]);
694  double in = gpd("in",r2,priProducts[pp]);
695  double hr_pr = in ? app(priProducts[pp],ft,dc)*st/in : 0.0;
696  hr_byPrd.push_back( hr_pr);
697  hr += hr_pr;
698  }
699 
700  // adjusting for overharvesting..
701  // 20160204: inserted to account that we let supply to be marginally higher than in in the mamarket module, to let the solver solving
702  double origHr = hr;
703  hr = min(1.0,hr);
704  for(uint pp=0;pp<priProducts.size();pp++){
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));
707  }
708 
709  double hV = hr*pastYearVol*px->avalCoef.at(j);
710 
711 
712  hV_byDiam.push_back(hV);
713  hV_byDiamAndPrd.push_back(hV_byPrd);
714 
715  // post harvesting remained volumes computation..
716  // loop(u$(ord(u)=1),
717  // first diameter class, no harvesting and fixed regenaration..
718  // V(u,i,lambda,essence,t)=(1-1/(tp(u,i,lambda,essence))-mort(u,i,lambda,essence) )*V(u,i,lambda,essence,t-1)
719  // +Vregen(i,lambda,essence,t);
720  // );
721  // loop(u$(ord(u)>1),
722  // generic case..
723  // V(u,i,lambda,essence,t)=((1-1/(tp(u,i,lambda,essence))
724  // -mort(u,i,lambda,essence) - hr(u,i,essence,lambda,t))*V(u,i,lambda,essence,t-1)
725  // +(1/(tp(u-1,i,lambda,essence)))*beta(u,i,lambda,essence)*V(u-1,i,lambda,essence,t-1));
726  double vol;
727  double tp = px->tp.at(j).at(u); //gfd("tp",regId,ft,dc);
728  double mort = px->mort.at(j).at(u); //gfd("mortCoef",regId,ft,dc);
729  double addMort = px->addMort.at(j).at(u); // fire mortality
730  double vReg = px->vReg.at(j); //gfd("vReg",regId,ft,""); // Taking it from the memory database as we could be in a fixed vReg scenario and not having calculated it from above!
731  double beta = px->beta.at(j).at(u); //gfd("betaCoef",regId,ft,dc);
732  //double hv2fa = gfd("hv2fa",regId,ft,dc);
733  double vHa = px->vHa.at(j).at(u); //gfd("vHa",regId,ft,dc);
734  double finalHarvestFlag = gfd("finalHarvestFlag",regId,ft,dc);
735 
736  double vMort = mort*pastYearVol;
737  double vMortAdd = addMort*pastYearVol;
738 
739  vMort_byDc.push_back(vMort);
740  vMortAdd_byDc.push_back(vMortAdd); // fire mortality
741 
742  if(useDeathTimber){
743  iisskey key(thisYear,r2,ft,dc);
744 
745  MD->deathTimberInventory_incrOrAdd(key,(vMort+vMortAdd)*shareMortalityUsableTimber);
746  }
747 
748  if(u==0){
749  vol = 0.0;
750  }else if(u==1){
751  vol = max(0.0,(1-1/tp-mort-addMort))*pastYearVol+vReg; //Antonello, "bug" fixed 20160203: In case of very strong mortality this quantity (that doesn't include harvesting) could be negative!
752  //double debug = vol;
753  #ifdef QT_DEBUG
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))+" ");
756  }
757  #endif
758  } else {
759  // time of passage and volume of smaller diameter class
760  double inc = (u==dClasses.size()-1)?0:1./tp; // we exclude the possibility for trees in the last diameter class to move to an upper class
761  double tp_1 = px->tp.at(j).at(u-1); //gfd("tp",regId,ft,dClasses[u-1]);
762  double pastYearVol_1 = px->vol_l.at(j).at(u-1); //gfd("vol",regId,ft,dClasses[u-1],thisYear-1);
763  //vol = max(0.0,(1-inc-mort-hr)*pastYearVol+(1/tp_1)*beta*pastYearVol_1);
764  vol = max(0.0,(1-inc-mort-addMort)*pastYearVol-hV+(1/tp_1)*beta*pastYearVol_1); // I can't use any more hr as it is the harvesting rate over the available volumes, not the whole ones
765  #ifdef QT_DEBUG
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.");
769  }
770  #endif
771  }
772  if(u != 0){ // this if is required to avoid a 0/0 and na error that then propage also in vSum()
773  double inc = (u==dClasses.size()-1)?0:1.0/tp; // we exclude the possibility for trees in the last diameter class to move to an upper class
774  double volumesMovingUp = inc*pastYearVol;
775  double pastArea = px->area_l.at(j).at(u);
776 
777  areasMovingUp.at(u) = inc*pastArea;
778 
779  if(oldVol2AreaMethod) {
780  hArea_byDc.push_back(finalHarvestFlag*1000000*hV/vHa); // volumes are in Mm^3, area in ha, vHa in m^3/ha
781  } else {
782  double finalHarvestedVolumes = finalHarvestFlag* hV;
783  double finalHarvestedRate = pastYearVol?finalHarvestedVolumes/pastYearVol:0.0; // Here we want the harvested rate over the whole volumes, not just the available ones, so we don't need to multiply to px->avalCoef.at(j)
784  #ifdef QT_DEBUG
785  if (finalHarvestedRate > 1.0){
786  msgOut(MSG_CRITICAL_ERROR,"Negative final harvested rate.");
787  }
788  #endif
789  hArea_byDc.push_back(finalHarvestedRate*pastArea); // volumes are in Mm^3, area in ha, vHa in m^3/ha
790  }
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));
792  #ifdef QT_DEBUG
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).");
795  }
796  #endif
797  } else {
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);
801  //if (pxId == 3550.0 && j==3){
802  // cout << "got the pixel" << endl;
803  //}
804  }
805  newVol_byDiam.push_back(vol);
806  #ifdef QT_DEBUG
807  if(px->area.at(j).at(u)< 0.0 || areasMovingUp.at(u) < 0.0 || hArea_byDc.at(u) < 0.0 ){
808  msgOut(MSG_CRITICAL_ERROR, "Negative values in runBiologicalModel");
809  }
810  #endif
811 
812  //double debug = hv2fa*hr*pastYearVol*100;
813  //cout << "regId|ft|dc| debug | freeArea: " << r2 << "|"<<ft<<"|"<<dc<<"| "<< debug << " | " << freeArea_byU << endl;
814 
815  //sfd(hr,"hr",regId,ft,dc);
816  //sfd(hV,"hV",regId,ft,dc);
817  //sfd(vol,"vol",regId,ft,dc);
818 
819  //sfd(freeArea_byU,"harvestedArea",regId,ft,dc,DATA_NOW,true);
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);
822  } // end foreach diameter classes
823  px->hVol.push_back(hV_byDiam);
824  px->hVol_byPrd.push_back(hV_byDiamAndPrd);
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); // fire mortality
829  px->hProductivity.push_back(productivity_byDc);
830 
831 
832  #ifdef QT_DEBUG
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){
840  msgOut(MSG_CRITICAL_ERROR,"mort vol > lagged volumes ?");
841  }
842  }
843  #endif
844  } // end of each forest type
845  } // end of each pixel
846 
847  #ifdef QT_DEBUG
848  // checking that in a region the total hVol is equal to the st for each products. 20150122 Test passed with the new availCoef
849  double sumSt = 0.0;
850  double sumHv = 0.0;
851  for(uint pp=0;pp<priProducts.size();pp++){
852  sumSt += gpd("st_fromHarv",r2,priProducts[pp]);
853  }
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++){
857  for(uint pp=0;pp<priProducts.size();pp++){
858  // by ft, dc, pp
859  sumHv += regPx[p]->hVol_byPrd[j][u][pp];
860  }
861  }
862  }
863  }
864  if(abs(sumSt-sumHv) > 0.000001){
865  msgOut(MSG_DEBUG, "St and harvested volumes diverge in region "+REG->getRegSName()+". St: "+d2s(sumSt)+" hV: "+d2s(sumHv));
866  }
867  #endif
868  } // end of each region
869 
870 }
871 
872 
873 
874 void
875 ModelCoreSpatial::runManagementModule(){
876  msgOut(MSG_INFO, "Starting management module..");
877  vector<string> allFTypes = MTHREAD->MD->getForTypeIds(true);
878  map<string,double> hAreaByFTypeGroup = vectorToMap(allFTypes,0.0);
879  int thisYear = MTHREAD->SCD->getYear();
880  double ir = MTHREAD->MD->getDoubleSetting("ir");
881 
882  // Post optimisation management module..
883  for(uint i=0;i<regIds2.size();i++){
884  int r2 = regIds2[i];
885  int regId = r2;
886  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
887  regPx = REG->getMyPixels();
888 
889  string flagHartman = MTHREAD->MD->getStringSetting("flagHartman",DATA_NOW,r2);
890  double pickling = MTHREAD->MD->getDoubleSetting("pickling"); // pickling factor used to factor in sequestration ex-situ and substitution effects
891 
892  //double exogenousReference = MTHREAD->MD->getBoolSetting("exogenousReference"); // new setting to distinguish between exogenous and endogenous carbon reference levels
893 
894 
895 
896 
897  // Caching futureCarbonPrices..
898  double carbonPrice_NOW = gpd("carbonPrice",r2,"",thisYear+1); // first carbon payments take place the year after replanting
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));
903  }
904  }
905 
906  vector<double> polBal_fiSub (fTypes.size(),0.0); // Expenses for forest investment subsides by region and destination ft €
907 
908  // Dealing with area change..
909  double fArea_reg = REG->getArea();
910  double fArea_diff = 0.0;
911  double fArea_reldiff = 0.0;
912  if(forestAreaChangeMethod=="relative"){
913  fArea_reldiff = gfd("forestChangeAreaIncrementsRel",r2,"","",DATA_NOW);
914  fArea_diff = fArea_reg * fArea_reldiff;
915  } else if (forestAreaChangeMethod=="absolute"){
916  fArea_diff = gfd("forestChangeAreaIncrementsHa",r2,"","",DATA_NOW);
917  fArea_reldiff = fArea_diff / fArea_reg;
918  }
919 
920  double regHArea = 0.0; // for the warning
921 
922  // Management rate
923  // 20170123: if mr is defined at regional level in the forData table, good, otherwise we look at settings table.
924  // If it isn't there as well, we rise a critical error.
925  // 20170302: removed the (false positive) error message raised if we put mr in setting table
926  // 20180626: mr is now, like all the other settings, powered with a regional dimension, so this looking in the
927  // forest data is no longer needed, but we keep it as some old projects have mr in the forData sheet.
928  double mr;
929  int origErrorLevel = MD->getErrorLevel();
930  MTHREAD->MD->setTempBool(true); // this should not be needed as it is set positive in getForData() map soubroutine..
931  MD->setErrorLevel(MSG_NO_MSG);
932  mr = MD->getForData("mr",regId,"","");
933  if(!MTHREAD->MD->getTempBool()) { // mr not found in the forest data
934  MD->setErrorLevel(MSG_CRITICAL_ERROR);
935  mr = MD->getDoubleSetting("mr",DATA_NOW,regId);
936  }
937  //try{
938  // mr = MD->getForData("mr",regId,"","");
939  //} catch (...){
940  // mr = MD->getDoubleSetting("mr");
941  //}
942  MD->setErrorLevel(origErrorLevel);
943 
944  // Caching carbon reference levels by ftypes
945  vector<vector<double>> futureCarbonRef;
946  vector<vector<double>> futureExtraBiomass_ratio;
947 
948  if(flagHartman!="NONE"){
949  for(uint j=0;j<fTypes.size();j++){
950  string ft = fTypes[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);
958  }
959  futureCarbonRef.push_back(futureCarbonRef_ft);
960  futureExtraBiomass_ratio.push_back(futureExtraBiomass_ratio_ft);
961  }
962  }
963 
964 
965  for (uint p=0;p<regPx.size();p++){
966  Pixel* px = regPx[p];
967  px->expectedReturns.clear();
968  px->expectedReturnsNotCorrByRa.clear(); // BUG discovered 20160825
969  //px->expectedReturnsCarbon.clear(); TODO
970  px->optDc.clear();
971  px->optFtChosen.clear();
972  px->optDcChosen.clear();
973  px->expectedAnnualIncome_carbon.clear();
974  px->expectedAnnualIncome_timber.clear();
975  resetMapValues(hAreaByFTypeGroup,0.0);
976  double totalHarvestedAreaOrig = vSum(px->hArea); // still need to remove the forest decrease areas..
977  double totalHarvestedArea = 0.0; // account for (eventual) forest decreases area)
978  vector<double> thisYearRegAreas(fTypes.size(),0.0); // initialize a vector of fTypes.size() zeros.
979  vector<double> expectedReturns(fTypes.size(),0.0); // uncorrected expected returns (without considering transaction costs). These are in form of eai
980  vector<double> expectedReturnsCarbon(fTypes.size(),0.0); // expected return of the carbon only
981  vector<int> rotation_dc(fTypes.size(),0.0); //vector to store optimal dc for each ftype (to be reused when flagHartman == "BAU")
982  vector< vector<double> > deltaAreas(fTypes.size()+1, vector<double>(fTypes.size()+1,0)); // matrix of regeneration areas ft_from/ft_to
983 
984 
985 
986  double fArea_px = vSum(px->area);
987  double fArea_diff_px = fArea_px * fArea_diff/ fArea_reg;
988 
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; // we can "remove" at maximum what has been harvested
993  regHArea += totalHarvestedAreaOrig;
994  totalHarvestedArea = totalHarvestedAreaOrig *(1-fArea_decr_rel);
995 
996 
997  // A - Computing the harvestingArea by parent ft group (for the allocation according to the prob of presence):
998  for(uint j=0;j<fTypes.size();j++){
999  string ft = fTypes[j];
1000  string parentFt = MTHREAD->MD->getForTypeParentId(ft);
1001  double hAreaThisFt=vSum(px->hArea.at(j))*(1-fArea_decr_rel);
1002  incrMapValue(hAreaByFTypeGroup,parentFt,hAreaThisFt); // increment the parent ft of the harvested area, neeed for assigning the frequences (prob. of presence)
1003  }
1004 
1005  // B - Computing the uncorrected expected returns (without considering transaction costs)
1006  // 20120910, Antonello: changed.. calculating the expected returns also for fixed and fromHrLevel regeneration (then not used but gives indication)
1007  // calculating the expected returns..
1008  // loop ( (u,i,essence,lambda,p_pr),
1009  // if (sum(u2, hV(u2,i,essence,lambda,t))= 0,
1010  // expRetPondCoef(u,i,essence,lambda,p_pr) = 0;
1011  // else
1012  // expRetPondCoef(u,i,essence,lambda,p_pr) = hV_byPrd(u,i,essence,lambda,p_pr,t)/ sum(u2, hV(u2,i,essence,lambda,t));
1013  // );
1014  // );
1015  // expReturns(i,essence,lambda) = sum( (u,p_pr),
1016  // RPAR("pl",i,p_pr,t)*hv2fa(i,essence,lambda,u)*(1/df_byFT(u,i,lambda,essence))* // df_byFT(u,i,lambda,essence)
1017  // expRetPondCoef(u,i,essence,lambda,p_pr)
1018  // );
1019 
1020  for(uint j=0;j<fTypes.size();j++){
1021  string ft = fTypes[j];
1022 
1023  // Added for carbon project
1024  double co2content_inventory = gfd("co2content_inventory",r2,ft,"",DATA_NOW)/1000; // kilos to tons
1025  // End of added
1026 
1027 
1028  double expReturns = std::numeric_limits<double>::lowest(); // these are reset at the beginning of each forest type, and need to be declared before the next loop
1029  double expReturns_carbon = std::numeric_limits<double>::lowest();
1030  double expReturns_timber = std::numeric_limits<double>::lowest();
1031  int optDc = 0; // "optimal diameter class", the one on which the expected returns are computed
1032  for (uint u=0; u<dClasses.size(); u++){
1033  string dc = dClasses[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);
1037  for (uint pp=0;pp<priProducts.size();pp++){
1038  double pol_mktDirInt_s = gpd("pol_mktDirInt_s",regId,priProducts[pp]); // This is valid also for the exported timber
1039  double pol_mktDirInt_s_fut = gpd("pol_mktDirInt_s",regId,priProducts[pp],thisYear+cumTp_u);
1040  double pl = gpd("pl",regId,priProducts[pp])*pol_mktDirInt_s; // note that this is the OBSERVED price. If we call it at current year+cumTp_u we would have the expected price. But we would first have to compute it, as pw is weigthed price world-local and we don't have local price for the future. DONE 20141202 ;-)
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;
1043  double sl = gpd("sl",regId,priProducts[pp]);
1044  double sa = gpd("sa",regId,priProducts[pp]);
1045  double pw_exp = computeExpectedPrice(pl, worldCurPrice, worldFutPrice, sl, sa, px->expTypePrices); //20141030: added the expected price!
1046  double raw_amount = finalHarvestFlag*pw_exp*vHa*app(priProducts[pp],ft,dc); // B.U.G. 20121126, it was missing app(pp,ft,dc) !!
1047  double anualised_amount_timber = MD->calculateAnnualisedEquivalent(raw_amount,cumTp_u, ir); // expected revenues from timber
1048  double anualised_amount = anualised_amount_timber; // total expected revenues (for now there is no carbon payment yet)
1049 
1050  // If carbon payments are included and the reference is an exogenous fixed value (e.g. year of reference), we compute revenues from carbon after retrieving reference from input files
1051  // Added for carbon project-------------------------------------------------------------------------
1052  double raw_amount_carbon = 0; // declared and initialized before IF statement and FOR loop
1053  if (flagHartman =="EXO") {
1054  double carbonPrice_y = carbonPrice_NOW;
1055  for (int y=0; y<cumTp_u; y++) {
1056  //double carbonPrice_y = futureCarbonPrices[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); // raw amount is calculated iteratively
1061  }
1062  //Now we add the final payment corresponding to carbon sequestered in long-lived products
1063  //double CarbonPrice_cumTp = futureCarbonPrices[floor(cumTp_u)]; INSTEAD WE USE CURRENT CARBON PRICE
1064  raw_amount_carbon += pickling * carbonPrice_y * finalHarvestFlag * vHa * app(priProducts[pp],ft,dc); // we add the final payment to raw amount for the final year, outside the FOR loop
1065  }
1066  // And we finally annualise everything
1067  double anualised_amount_carbon = MD->calculateAnnualisedEquivalent(raw_amount_carbon,cumTp_u, ir);
1068  anualised_amount += anualised_amount_carbon; // anualised_amount already had the timber revenues in it, now we add carbon
1069 
1070  // End of added-----------------------------------------------------------------------------------
1071 
1072 
1073 
1074  if (anualised_amount>expReturns) {
1075  expReturns=anualised_amount;
1076  optDc = u;// this is the index of the optimal dc, not the dc itself
1077  expReturns_carbon = anualised_amount_carbon;
1078  expReturns_timber = anualised_amount_timber;
1079  }
1080  } // end for each product
1081  } // end for each diameter class
1082 
1083  //results are pushed back to pixel only if they are final
1084  if (flagHartman=="NONE" || flagHartman=="EXO") {
1085  px->expectedAnnualIncome_timber.push_back(expReturns_timber);
1086  px->expectedAnnualIncome_carbon.push_back(expReturns_carbon);
1087  px->expectedReturnsNotCorrByRa.push_back(expReturns);
1088  px->optDc.push_back(optDc);
1089  }
1090 
1091  //calculation of expected returns with risk aversion
1092  if(MD->getBoolSetting("heterogeneousRiskAversion",DATA_NOW)){
1093  double ra = px->getDoubleValue("ra");
1094  double cumMort = 1-px->cumAlive_exp.at(j).at(optDc);
1095  //cout << px->getID() << "\t" << ft << "\t\t" << "optDc" << optDc << "\t" << cumMort << endl;
1096  double origExpReturns = expReturns;
1097  expReturns = origExpReturns * (1.0 - ra*cumMort);
1098  }
1099 
1100  // results are pushed back to pixel only if they are final
1101  if (flagHartman=="NONE" || flagHartman=="EXO") {
1102  px->expectedReturns.push_back(expReturns); // if there is no risk version, then expectedReturns and expectedReturnsNotCorrByRa will have the same value
1103  }
1104 
1105  expectedReturns.at(j) = expReturns;
1106  rotation_dc.at(j) = optDc; // here we store the index of optimal dc for the forest type considered
1107 
1108  } // end foreach forest type
1109 
1110 
1111  // Looping again by forest types, here by ft in input (those harvested)
1112  // Added for second carbon project-------------------------------------------------------------------------------------------------------------------------------------
1113  // If carbon payments are included and the reference is BAU (ie, endogenous, we need a new loop-------------------------------------------------------------------
1114  // In this case, the reference is what happens without policy, i.e. the result of the previous loop when (flagHartman =="REFERENCE") is FALSE
1115  // So we just redo another loop with carbon payments using the previous result as reference
1116 
1117  if (flagHartman =="ENDO") { // if carbon payments take place and reference is endogenous (is BAU)
1118  // STEP 1 - retrieve optimal regeneration choice without carbon payments
1119  int ft_opt_index = getMaxPos(expectedReturns); // Faustmannian optimal forest type index
1120  string ft_opt = fTypes[ft_opt_index]; // Faustmannian optimal forest type name
1121  int dc_opt_index = rotation_dc[ft_opt_index]; // Faustmannian optimal diameter class index
1122  string dc_opt = dClasses[dc_opt_index]; // Faustmannian optimal diameter class name
1123  double cumTP_opt = px->cumTp_exp.at(ft_opt_index).at(dc_opt_index); // Faustmannian optimal rotation length
1124  double co2content_inventory_opt = gfd("co2content_inventory",r2,ft_opt,"",DATA_NOW)/1000; // kilos to tons
1125 
1126  // STEP 2 - calculate expected carbon contents on the optimal choice without carbon payments (Faustmannian optimum), to use as reference
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);
1131  }
1132  // the Hartmanian choice will usually be longer, so the Faustmannian rotation needs to be replicated several times (at least 2000 years, should be enough)
1133  while (FutureCarbonContent_opt.size()<=2000){
1134  FutureCarbonContent_opt.insert(std::end(FutureCarbonContent_opt), std::begin(FutureCarbonContent_opt), std::end(FutureCarbonContent_opt));
1135  }
1136  //int cumTP_opt_floor = FutureCarbonContent_opt.size(); // need to get an integer value
1137 
1138  // STEP 3 - new loop to find optimal regeneration choice WITH carbon payments, using the choice WITHOUT carbon payments as a reference
1139 
1140  for(uint j=0;j<fTypes.size();j++){
1141  string ft = fTypes[j];
1142  double co2content_inventory = gfd("co2content_inventory",r2,ft,"",DATA_NOW)/1000; // kilos to tons
1143 
1144  double expReturns_hartman = -10000000000; // need to initialize because we use it for a comparison
1145  double expReturns_timber_hartman; // no need to give an initial value because they are not used before being given a value
1146  double expReturns_carbon_hartman;
1147  int optDc_hartman = -1; // "optimal diameter class", the one on which the expected returns are computed
1148  for (uint u=1; u<dClasses.size(); u++){
1149  string dc = dClasses[u];
1150  double cumTp_u = px->cumTp_exp.at(j).at(u);
1151 
1152  // First we calculate expected revenues from carbon (they used to be inside the products loop, but they don't need to be since they don't depend on products to be sold)
1153  double raw_amount_carbon=0;
1154  double carbonPrice_y = carbonPrice_NOW; //carbon price is by default current price, ie it is assumed to be constant from owner's perspective
1155  for (int y=0; y<=cumTp_u; y++) { // compute yearly carbon payments
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]);
1159 
1160  //r if (y==0) { // we use euclidian division, so an exception must be made when denominator is equal to 0
1161  //r carbonPayment_y = 0; // in the first year, no forest can reach a non-zero volume so payments are 0
1162  //r} else {
1163  //r int factor = y/cumTP_opt_floor; // euclidian division of the two lengths we compare
1164  //r if (cumTP_opt_floor%y >0) { // if the year is NOT a multiple of optimal rotation time, we need to take the year-th value of vector modulus how many rotations we have passed
1165  //r carbonPayment_y = carbonPrice_y*ir*(co2content_inventory*expectedVHa_y - FutureCarbonContent_opt[y - factor*cumTP_opt_floor]);
1166  //r } else { // if the year is EXACTLY a multiple of optimal rotation time, we need to take the last value of vector, not the first one
1167  //r carbonPayment_y = carbonPrice_y*ir*(co2content_inventory*expectedVHa_y - FutureCarbonContent_opt[cumTP_opt_floor]);
1168  //r }
1169  //r}
1170  raw_amount_carbon += carbonPayment_y*pow(1+ir,cumTp_u - y); // all payments compounded to end of rotation
1171  }
1172 
1173  // And now we compute revenues from timber (we try several products, as in the Faustmannian rotation)
1174  double vHa = px->vHa_exp.at(j).at(u);
1175  double finalHarvestFlag = gfd("finalHarvestFlag",regId,ft,dc);
1176  for (uint pp=0;pp<priProducts.size();pp++){
1177  double pl = gpd("pl",regId,priProducts[pp]); // note that this is the OBSERVED price. If we call it at current year+cumTp_u we would have the expected price. But we would first have to compute it, as pw is weigthed price world-local and we don't have local price for the future. DONE 20141202 ;-)
1178  double worldCurPrice = gpd("pl",WL2,priProducts[pp]);
1179  double worldFutPrice = gpd("pl",WL2,priProducts[pp],thisYear+cumTp_u);
1180  double sl = gpd("sl",regId,priProducts[pp]);
1181  double sa = gpd("sa",regId,priProducts[pp]);
1182  double pw_exp = computeExpectedPrice(pl, worldCurPrice, worldFutPrice, sl, sa, px->expTypePrices); //20141030: added the expected price!
1183  double raw_amount = finalHarvestFlag*pw_exp*vHa*app(priProducts[pp],ft,dc); // B.U.G. 20121126, it was missing app(pp,ft,dc) !!
1184  double anualised_amount_timber = MD->calculateAnnualisedEquivalent(raw_amount,cumTp_u, ir);
1185 
1186  // Now we add the final carbon payment corresponding to carbon sequestered in long-lived products
1187  raw_amount_carbon += pickling * carbonPrice_y * finalHarvestFlag * vHa * app(priProducts[pp],ft,dc);
1188  // And we annualise carbon revenues and add them to annualised timber revenues
1189  double anualised_amount_carbon = MD->calculateAnnualisedEquivalent(raw_amount_carbon,cumTp_u, ir);
1190  //px->expectedReturnsCarbon = anualised_amount_carbon; TODO
1191 
1192  double anualised_amount = anualised_amount_timber + anualised_amount_carbon;
1193 
1194  // If the new expected revenues is higher than previously, we store the maximum value and corresponding dc.
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;
1199  optDc_hartman = u;
1200  }
1201  } // end foreach primary product
1202  } // end foreach DC
1203 
1204  // Here we send results to pixel because they are final
1205  px->expectedAnnualIncome_timber.push_back(expReturns_timber_hartman);
1206  px->expectedAnnualIncome_carbon.push_back(expReturns_carbon_hartman);
1207  px->expectedReturnsNotCorrByRa.push_back(expReturns_hartman);
1208  px->optDc.push_back(optDc_hartman);
1209 
1210  // And with risk aversion
1211  if(MD->getBoolSetting("heterogeneousRiskAversion",DATA_NOW)){
1212  double ra = px->getDoubleValue("ra");
1213  double cumMort = 1-px->cumAlive_exp.at(j).at(optDc_hartman);
1214  //cout << px->getID() << "\t" << ft << "\t\t" << "optDc" << optDc << "\t" << cumMort << endl;
1215  double origExpReturns_hartman = expReturns_hartman;
1216  expReturns_hartman = origExpReturns_hartman * (1.0 - ra*cumMort);
1217  }
1218 
1219  px->expectedReturns.push_back(expReturns_hartman);
1220 
1221  // Now we replace results from the previous loop with the new ones
1222 
1223  expectedReturns.at(j) = expReturns_hartman;
1224  rotation_dc.at(j) = optDc_hartman;
1225 
1226 
1227  } // end foreach forest type
1228 
1229  } // end IF Hartman reference is BAU (no policy)----------------------------------------------------------------
1230 
1231  // End of added second carbon project---------------------------------------------------------------------------
1232 
1233  for(uint j=0;j<fTypes.size();j++){
1234  string ft = fTypes[j];
1235  forType* thisFt = MTHREAD->MD->getForType(ft);
1236 
1237  double harvestedAreaForThisFT = vSum(px->hArea.at(j))*(1-fArea_decr_rel); // gfd("harvestedArea",regId,ft,DIAM_ALL);
1238  deltaAreas[j][fTypes.size()] += vSum(px->hArea.at(j))*(fArea_decr_rel);
1239  vector<double> corrExpectedReturns(fTypes.size(),0.0); // corrected expected returns (considering transaction costs). These are in form of NPV
1240  vector<double> polBal_fiSub_unitvalue(fTypes.size(),0.0); // subside (tax) per ha per destination ft €/ha
1241 
1242  // C - Computing the corrected expected returns including transaction costs
1243 
1244  for(uint j2=0;j2<fTypes.size();j2++){
1245  string ft2 = fTypes[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; // changed 20150718: npv = eai/ir + tr. cost // HUGE BUG 20151202: transaction costs should be REDUCED, not added to the npv...
1249  }
1250 
1251  //int highestReturnFtIndex = getMaxPos(corrExpectedReturns);
1252 
1253  // D - Assigning the Managed area
1254  // calculating freeArea at the end of the year and choosing the new regeneration area..
1255  //freeArea(i,essence,lambda) = sum(u, hv2fa(i,essence,lambda,u)*hr(u,i,essence,lambda,t)*V(u,i,lambda,essence,t-1)*100);
1256  //if(scen("endVreg") ,
1257  // regArea(i,essence,lambda,t) = freeArea(i,essence, lambda); // here we could introduce in/out area from other land usages
1258  //else
1259  // loop (i,
1260  // loop( (essence,lambda),
1261  // if ( expReturns(i,essence,lambda) = smax( (essence2,lambda2),expReturns(i,essence2,lambda2) ),
1262  // regArea (i,essence,lambda,t) = sum( (essence2, lambda2), freeArea(i,essence2, lambda2) ) * mr;
1263  // );
1264  // );
1265  // regArea(i,essence,lambda,t) = freeArea(i,essence, lambda)*(1-mr); // here we could introduce in/out area from other land usages
1266  // );
1267  //if (j==highestReturnFtIndex){
1268  // thisYearRegAreas[j] += totalHarvestedArea*mr;
1269  //}
1270  // If I Implement this I'll have a minimal diff in total area.. why ?????
1271 
1272  int optFtChosen = getMaxPos(corrExpectedReturns);
1273  px->optFtChosen.push_back(optFtChosen);
1274  px->optDcChosen.push_back(px->optDc.at(optFtChosen));
1275 
1276  thisYearRegAreas[getMaxPos(corrExpectedReturns)] += harvestedAreaForThisFT*mr;
1277  thisYearRegAreas[getMaxPos(expectedReturns)] += fArea_incr*mr/((double)fTypes.size()); // mr quota of new forest area assigned to highest expected returns ft (not considering transaction costs). Done for each forest types
1278 
1279  deltaAreas[j][getMaxPos(corrExpectedReturns)] += harvestedAreaForThisFT*mr;
1280  deltaAreas[fTypes.size()][getMaxPos(expectedReturns)] += fArea_incr*mr/((double)fTypes.size());
1281 
1282  polBal_fiSub[getMaxPos(corrExpectedReturns)] += harvestedAreaForThisFT*mr*polBal_fiSub_unitvalue[getMaxPos(corrExpectedReturns)];
1283 
1284  // E - Assigning unmanaged area
1285  //for(uint j2=0;j2<fTypes.size();j2++){
1286  if(natRegAllocation=="pp"){ // according to prob presence
1287  //string ft2 = fTypes[j2];
1288  string parentFt = MTHREAD->MD->getForTypeParentId(ft);
1289  double freq = rescaleFrequencies ? gfd("freq_norm",regId,parentFt,""):gfd("freq",regId,parentFt,""); // "probability of presence" for unmanaged forest, added 20140318
1290  double hAreaThisFtGroup = findMap(hAreaByFTypeGroup,parentFt);
1291  double hRatio = 1.0;
1292  if(hAreaThisFtGroup>0){
1293  //double harvestedAreaForThisFT2 = vSum(px->hArea.at(j2));
1294  hRatio = harvestedAreaForThisFT/hAreaThisFtGroup;
1295  } else {
1296  int nFtChilds = MTHREAD->MD->getNForTypesChilds(parentFt);
1297  hRatio = 1.0/nFtChilds;
1298  }
1299  thisYearRegAreas[j] += totalHarvestedArea*(1-mr)*freq*hRatio;
1300  thisYearRegAreas[j] += fArea_incr*(1-mr)*freq*hRatio; // non-managed quota of new forest area assigning proportionally on pp at sp group level
1301  //thisYearRegAreas[j2] += harvestedAreaForThisFT*(1-mr)*freq*hRatio;
1302  // Here is the problem with building a transaction matrix ft_from/ft_to: we don't explicitly assign from/to for unmanaged areas !
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;
1305  }
1306  deltaAreas[fTypes.size()][j] += fArea_incr*(1-mr)*freq*hRatio;
1307 
1308 
1309 
1310  } else { // prob presence not used..
1311 
1312  // Accounting for mortality arising from pathogens. Assigning the area to siblings according to area..
1313  // TODO: chek that mortRatePath doesn't involve adding forest area to deltaAreas[j][""]
1314 
1315  double mortRatePath = px->getPathMortality(ft, "0");
1316  if(mortRatePath > 0){
1317 
1318  string parentFt = MTHREAD->MD->getForTypeParentId(ft);
1319  vector <string> siblings = MTHREAD->MD->getForTypeChilds(parentFt);
1320  vector <double> siblingAreas;
1321  for(uint j2=0;j2<siblings.size();j2++){
1322  if(siblings[j2]==ft){
1323  siblingAreas.push_back(0.0);
1324  } else {
1325  //string debug_sibling_ft = siblings[j2];
1326  //int debug_positin = getPos(debug_sibling_ft,fTypes);
1327  double thisSiblingArea = vSum(px->area.at(getPos(siblings[j2],fTypes)));
1328  siblingAreas.push_back(thisSiblingArea);
1329  }
1330  }
1331  double areaAllSiblings = vSum(siblingAreas);
1332  thisYearRegAreas[j] += harvestedAreaForThisFT*(1-mr)*(1-mortRatePath);
1333  deltaAreas[j][j] += harvestedAreaForThisFT*(1-mr)*(1-mortRatePath);
1334 
1335  if(areaAllSiblings>0.0){ // area of siblings is >0: we attribute the area from the pathogen induced mortality to the siblings proportionally to area..
1336  for(uint j2=0;j2<siblings.size();j2++){
1337 // int debug1 = getPos(siblings[j2],fTypes);
1338 // double debug2= harvestedAreaForThisFT;
1339 // double debug3 = 1.0-mr;
1340 // double debug4 = mortRatePath;
1341 // double debug5 = siblingAreas[j2];
1342 // double debug6 = areaAllSiblings;
1343 // double debug7 = harvestedAreaForThisFT*(1.0-mr)*(mortRatePath)*(siblingAreas[j2]/areaAllSiblings);
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);
1346  }
1347  } else if (siblings.size()>1) { // area of all siblings is 0, we just give them the mortality area in equal parts..
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);
1352  }
1353  }
1354  }
1355  } else { // mortRatePath == 0
1356  thisYearRegAreas[j] += harvestedAreaForThisFT*(1.0-mr);
1357  deltaAreas[j][j] += harvestedAreaForThisFT*(1.0-mr);
1358  }
1359 
1360  // Allocating non-managed quota of new forest area to ft proportionally to the current area share by ft
1361  double newAreaThisFt = vSum(px->area) ? fArea_incr*(1-mr)*vSum(px->area.at(j))/vSum(px->area): 0.0;
1362  thisYearRegAreas[j] += newAreaThisFt;
1363  deltaAreas[fTypes.size()][j] += newAreaThisFt;
1364 
1365  if(! (thisYearRegAreas[j] >= 0.0) ){
1366  msgOut(MSG_ERROR,"thisYearRegAreas[j] is not non negative (j: "+i2s(j)+", thisYearRegAreas[j]: "+i2s( thisYearRegAreas[j])+").");
1367  }
1368  //thisYearRegAreas[j2] += harvestedAreaForThisFT*(1-mr);
1369  }
1370  //}
1371  } // end for each forest type
1372 
1373  // adding regeneration area to the fist (00) diameter class
1374  for(uint j=0;j<fTypes.size();j++){
1375  px->area.at(j).at(0) += thisYearRegAreas.at(j);
1376  }
1377 
1378  #ifdef QT_DEBUG
1379  double totalRegArea = vSum(thisYearRegAreas);
1380  if (! (totalRegArea==0.0 && totalHarvestedArea==0.0)){
1381  double ratio = totalRegArea / totalHarvestedArea ;
1382  if(rescaleFrequencies && (ratio < 0.99999999999 || ratio > 1.00000000001) && forestAreaChangeMethod == "fixed") {
1383  msgOut(MSG_CRITICAL_ERROR, "Sum of regeneration areas not equal to sum of harvested area in runManagementModel()!");
1384  }
1385  }
1386  #endif
1387 
1388  px->regArea.insert(pair <int, vector<double> > (MTHREAD->SCD->getYear(), thisYearRegAreas));
1389  px->deltaArea = deltaAreas;
1390 
1391  } // end of each pixel
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.");
1394  }
1395 
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); //M€
1399 
1400  }
1401  } // end of each region
1402 }
1403 
1404 void
1405 ModelCoreSpatial::cacheSettings(){
1406  msgOut(MSG_INFO, "Cashing initial model settings..");
1407  MD = MTHREAD->MD;
1408  firstYear = MD->getIntSetting("initialYear");
1409  secondYear = firstYear+1;
1410  thirdYear = firstYear+2;
1411  WL2 = MD->getIntSetting("worldCodeLev2");
1412  regIds2 = MD->getRegionIds(2);
1413  priProducts = MD->getStringVectorSetting("priProducts");
1414  secProducts = MD->getStringVectorSetting("secProducts");
1415  allProducts = priProducts;
1416  allProducts.insert( allProducts.end(), secProducts.begin(), secProducts.end() );
1417  dClasses = MD->getStringVectorSetting("dClasses");
1418  pDClasses; // production diameter classes: exclude the fist diameter class below 15 cm
1419  pDClasses.insert(pDClasses.end(), dClasses.begin()+1, dClasses.end() );
1420  fTypes= MD->getForTypeIds();
1421  l2r = MD->getRegionIds();
1422 }
1423 
1424 void
1425 ModelCoreSpatial::cacheDynamicSettings(){
1426  // Settings needs explicitly DATA_NOW to have a temporal dymension..
1427  msgOut(MSG_INFO, "Cashing model settings that may change every year..");
1428  regType = MTHREAD->MD->getStringSetting("regType",DATA_NOW); // how the regeneration should be computed (exogenous, from hr, from allocation choises)
1429  natRegAllocation = MTHREAD->MD->getStringSetting("natRegAllocation",DATA_NOW); // how to allocate natural regeneration
1430  rescaleFrequencies = MD->getBoolSetting("rescaleFrequencies",DATA_NOW);
1431  oldVol2AreaMethod = MD->getBoolSetting("oldVol2AreaMethod",DATA_NOW);
1432  //mr = MD->getDoubleSetting("mr");
1433  forestAreaChangeMethod = MTHREAD->MD->getStringSetting("forestAreaChangeMethod",DATA_NOW);
1434  ir = MD->getDoubleSetting("ir",DATA_NOW);
1435 }
1436 
1437 void
1438 ModelCoreSpatial::initializePixelVolumes(){
1439  msgOut(MSG_INFO, "Starting initializing pixel-level values");
1440 
1441  // pxVol = regVol * pxArea/regForArea
1442  // this function can be done only at the beginning of the model, as it assume that the distribution of volumes by diameter class in the pixels within a certain region is homogeneous, but as the model progress along the time dimension this is no longer true.
1443  if(!MD->getBoolSetting("usePixelData")) return;
1444  for(uint i=0;i<regIds2.size();i++){
1445  ModelRegion* reg = MD->getRegion(regIds2[i]);
1446  vector <Pixel*> rpx = MTHREAD->GIS->getAllPlotsByRegion(regIds2[i]);
1447  for (uint j=0;j<rpx.size();j++){
1448  //int debugPx = rpx[j]->getID();
1449  //int debug2 = debugPx;
1450  rpx[j]->vol.clear(); // not actually necessary
1451  for(uint y=0;y<fTypes.size();y++){
1452  vector <double> vol_byu;
1453  double regForArea = reg->getValue("forArea_"+fTypes[y]);
1454  for (uint z=0;z<dClasses.size();z++){
1455  double regVol;
1456  regVol = z ? gfd("vol",regIds2[i],fTypes[y],dClasses[z],firstYear) : 0 ; // if z=0-> regVol= gfd(), otherwise regVol=0;
1457  double pxArea = rpx[j]->getDoubleValue("forArea_"+fTypes[y], true); // bug solved 20121109. get zero for not data
1458  if (pxArea<0.0){
1459  msgOut(MSG_CRITICAL_ERROR,"Error in initializePixelVolumes, negative pxArea!");
1460  }
1461  double pxVol = regForArea ? regVol * pxArea/regForArea: 0; // if we introduce new forest types without initial area we must avoid a 0/0 division
1462  //rpx[j]->changeValue(pxVol,"vol",fTypes[y],dClasses[z],firstYear);
1463  vol_byu.push_back(pxVol);
1464  }
1465  rpx[j]->vol.push_back(vol_byu);
1466  }
1467  }
1468  }
1469  loadExogenousForestLayers("vol");
1470 }
1471 
1472 /**
1473  * @brief ModelCoreSpatial::assignSpMultiplierPropToVols assigns the spatial multiplier (used in the time of return) based no more on a Normal distribution but on the volumes present in the pixel: more volume, more the pixel is fit for the ft
1474  *
1475  * This function apply to the pixel a multiplier of time of passage that is inversely proportional to the volumes of that forest type present in the pixel.
1476  * The idea is that in the spots where we observe more of a given forest type are probabily the most suited ones to it.
1477  *
1478  * The overall multipliers **of time of passage** (that is, the one returned by Pixel::getMultiplier("tp_multiplier") ) will then be the product of this multiplier that account for spatial heterogeneity and of an eventual exogenous
1479  * multiplier that accounts for different scenarios among the spatio-temporal dimensions.
1480  *
1481  * Given that (forest type index omitted):
1482  * - \f$V_{p}\f$ = volume of a given ft in each pixel (p)
1483  * - \f$\bar{g}\f$ and \f$\sigma_{g}\f$ = regional average and standard deviation of the growth rate
1484  * - \f$m_{p}\f$ = multiplier of time of passage
1485  *
1486  * This multiplier is computed as:
1487  * - \f$ v_{p} = max(V) - V_{p}~~ \f$ A diff from the max volume is computed in each pixel
1488  * - \f$ vr_{p} = v_{p} * \bar{g}/\bar{v}~~ \f$ The volume diff is rescaled to match the regional growth rate
1489  * - \f$ vrd_{p} = vr_{p} - \bar{vr}~~ \f$ Deviation of the rescaled volumes are computed
1490  * - \f$ vrdr_{p} = vrd_{p} * \sigma_{g}/\sigma_{vr}~~ \f$ The deviations are then rescaled to match the standard deviations of the regional growth rate
1491  * - \f$ m_{p} = (vrdr_{p} + \bar{vr}) / \bar{g}~~ \f$ The multiplier is computed from the ratio of the average rescaled volumes plus rescaled deviation over the average growth rate.
1492  *
1493  * And it has the following properties:
1494  * - \f$\bar{m} = 1\f$
1495  * - \f$\sigma_{m} = cv_{g}\f$
1496  * - \f$m_{p} = V_{p}*\alpha+\beta\f$
1497  * - \f$m_{\bar{V}} = 1\f$
1498  *
1499  * For spreadsheet "proof" see the file computation_of_growth_multipliers_from_know_avg_sd_and_proportional_to_share_of_area_in_each_pixel.ods
1500  */
1501 void
1502 ModelCoreSpatial::assignSpMultiplierPropToVols(){
1503 
1504  if(!MTHREAD->MD->getBoolSetting("useSpatialVarPropToVol")){return;}
1505  for(uint r=0;r<regIds2.size();r++){
1506  int rId = regIds2[r];
1507  ModelRegion* reg = MD->getRegion(regIds2[r]);
1508  vector <Pixel*> rpx = MTHREAD->GIS->getAllPlotsByRegion(regIds2[r]);
1509  for(uint f=0;f<fTypes.size();f++){
1510  string ft = fTypes[f];
1511  double agr = gfd("agr",regIds2[r],ft,"");
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;
1518  double final_value;
1519  double multiplier;
1520  vector<double> multipliers; // for tests
1521 
1522  double vol_max, rescale_ratio_avg, rescale_ratio_sd;
1523  double diffVols_avg, diffVols_rescaled_avg;
1524  double diffVols_rescaled_sd;
1525 
1526  for (uint p=0;p<rpx.size();p++){
1527  Pixel* px = rpx[p];
1528  vols.push_back(vSum(px->vol[f]));
1529  } // end for each pixel
1530  vol_max=getMax(vols);
1531 
1532  for(uint p=0;p<vols.size();p++){
1533  diffVols.push_back(vol_max-vols[p]);
1534  }
1535 
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);
1540  }
1541  diffVols_rescaled_avg = getAvg(diffVols_rescaled);
1542  diffVols_rescaled_sd = getSd(diffVols_rescaled,false);
1543 
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; //20151130: added bounds for extreme cases. Same bonds as in Gis::applySpatialStochasticValues()
1550  // multiplier = 1.0;
1551 
1552  Pixel* px = rpx[p];
1553  px->setSpModifier(multiplier,f);
1554  multipliers.push_back(multiplier);
1555  }
1556 
1557  #ifdef QT_DEBUG
1558  // Check relaxed as we introduced bounds that may change slighly the avg and sd...
1559  double avgMultipliers = getAvg(multipliers);
1560  double sdMultipliers = getSd(multipliers,false);
1561  if ( avgMultipliers < 0.9 || avgMultipliers > 1.1){
1562  msgOut(MSG_CRITICAL_ERROR, "The average of multipliers of ft "+ ft +" for the region " + i2s(rId) + " is not 1!");
1563  }
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!");
1567  }
1568  #endif
1569  } // end for each ft
1570  } // end for each region
1571 }
1572 
1573 
1574 
1575 void
1576 ModelCoreSpatial::initialiseCarbonModule(){
1577 
1578  ///< call initialiseDeathBiomassStocks(), initialiseProductsStocks() and initialiseEmissionCounters()
1579  MTHREAD->CBAL->initialiseEmissionCounters();
1580 
1581  for(uint i=0;i<regIds2.size();i++){
1582  vector<double> deathBiomass;
1583  for(uint j=0;j<fTypes.size();j++){
1584  double deathBiomass_ft = gfd("vMort",regIds2[i],fTypes[j],DIAM_ALL,DATA_NOW)+gfd("vMortAdd",regIds2[i],fTypes[j],DIAM_ALL,DATA_NOW);
1585  deathBiomass.push_back(deathBiomass_ft);
1586  }
1587  MTHREAD->CBAL->initialiseDeathBiomassStocks(deathBiomass, regIds2[i]);
1588  vector<double>qProducts;
1589  for(int p=0;p<priProducts.size();p++){
1590  // for the primary products we consider only the exports as the domestic consumption is entirelly transformed in secondary products
1591  double int_exports = gpd("sa",regIds2[i],priProducts[p],DATA_NOW);
1592  qProducts.push_back(int_exports);
1593  }
1594  for(int p=0;p<secProducts.size();p++){
1595  // for the tranformed product we skip those that are imported, hence derived from other forest systems
1596  double consumption = gpd("dl",regIds2[i],secProducts[p],DATA_NOW); // dl = sl + net regional imports
1597  qProducts.push_back(consumption);
1598  }
1599  MTHREAD->CBAL->initialiseProductsStocks(qProducts, regIds2[i]);
1600 
1601  }
1602 }
1603 
1604 void
1605 ModelCoreSpatial::initialiseDeathTimber(){
1606  int currentYear = MTHREAD->SCD->getYear();
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++){
1611  iisskey key(y,regIds2[i],fTypes[j],dClasses[u]);
1612  MD->deathTimberInventory_incrOrAdd(key,0.0);
1613  }
1614  }
1615  }
1616  }
1617 }
1618 
1619 /**
1620  * @brief ModelCoreSpatial::initializePixelArea
1621  *
1622  * This function compute the initial area by ft and dc. It requires vHa computed in computeCumulativeData, this is why it is
1623  * separated form the other initialisedPixelValues().
1624  * As the sum of area computed using vHa may differ from the one memorised in forArea_* layer, all values are scaled to match
1625  * it before being memorised.
1626  * Also assign area = area_l
1627  */
1628 
1629 void
1630 ModelCoreSpatial::initializePixelArea(){
1631  msgOut(MSG_INFO, "Starting initializing pixel-level area");
1632  if(!MD->getBoolSetting("usePixelData")) return;
1633  for(uint i=0;i<regIds2.size();i++){
1634  ModelRegion* reg = MD->getRegion(regIds2[i]);
1635  vector <Pixel*> rpx = MTHREAD->GIS->getAllPlotsByRegion(regIds2[i]);
1636  for (uint p=0;p<rpx.size();p++){
1637  Pixel* px = rpx[p];
1638  double pxid= px->getID();
1639  for(uint j=0;j<fTypes.size();j++){
1640  string ft = fTypes[j];
1641  vector <double> tempAreas;
1642  vector <double> areasByFt;
1643  double pxArea = px->getDoubleValue("forArea_"+ft,true)/10000.0; //ha
1644  for (uint u=0;u<dClasses.size();u++){
1645  if(u==0){
1646  double regionArea = reg->getValue("forArea_"+ft,OP_SUM)/10000.0; //ha
1647  double regRegVolumes = gfd("vReg",regIds2[i],ft,""); // regional regeneration volumes.. ugly name !!
1648  double newVReg = regionArea ? regRegVolumes*pxArea/regionArea : 0.0;
1649  double tp_u0 = px->tp.at(j).at(0); // time of passage to reach the first production diameter class
1650  double entryVolHa = gfd("entryVolHa",regIds2[i],ft,"");
1651  double tempArea = (newVReg*1000000.0/entryVolHa)*tp_u0;
1652  tempAreas.push_back(tempArea);
1653  } else {
1654  string dc = dClasses[u];
1655  double dcVol = px->vol_l.at(j).at(u)*1000000.0; // m^3
1656  double dcVHa = px->vHa.at(j).at(u); // m^3/ha
1657  #ifdef QT_DEBUG
1658  if(dcVol < 0.0 || dcVHa < 0.0){
1659  msgOut(MSG_CRITICAL_ERROR, "Negative volumes or density in initializePixelArea");
1660  }
1661  #endif
1662  double tempArea = dcVHa?dcVol/dcVHa:0;
1663  tempAreas.push_back(tempArea);
1664  }
1665 
1666  } // end dc
1667  double sumTempArea = vSum(tempAreas);
1668  // double sharedc0 = 5.0/90.0; // an arbitrary share of total area allocated to first diameter class
1669  //tempAreas.at(0) = sumTempArea * sharedc0;
1670  //sumTempArea = vSum(tempAreas);
1671  double normCoef = sumTempArea?pxArea/ sumTempArea:0;
1672  //cout << i << '\t' << pxid << '\t' << ft << '\t' << normCoef << endl;
1673  #ifdef QT_DEBUG
1674  if(normCoef < 0.0){
1675  msgOut(MSG_CRITICAL_ERROR, "Negative normCoef in initializePixelArea");
1676  }
1677  #endif
1678  for (uint u=0;u<dClasses.size();u++){
1679  areasByFt.push_back(tempAreas.at(u)*normCoef); //manca la costruzione originale del vettore
1680  }
1681  #ifdef QT_DEBUG
1682  if (pxArea != 0.0){
1683  double ratio = vSum(areasByFt)/ pxArea; // vSum(areasByFt) should be equal to pxArea
1684  if(ratio < 0.99999999999 || ratio > 1.00000000001) {
1685  msgOut(MSG_CRITICAL_ERROR, "pxArea is not equal to vSum(areasByFt) in initializePixelArea");
1686  }
1687  }
1688  #endif
1689  px->area_l.push_back(areasByFt);
1690  /// \todo here I have finally also area_ft_dc_px and I can implement the new one I am in 2006
1691  } // end ft
1692  px->area = px->area_l; //Assigning initial value of area to the area of the old year
1693  } // end px
1694  } // end region
1695  loadExogenousForestLayers("area");
1696  /// \todo: also update area_l
1697 }
1698 
1699 void
1700 ModelCoreSpatial::computeCumulativeData(){
1701 
1702  msgOut(MSG_INFO, "Starting computing some cumulative values..");
1703  int thisYear = MTHREAD->SCD->getYear();
1704 
1705 // double sumCumTP=0;
1706 // double sumVHa = 0;
1707 // double count = 0;
1708 // double avg_sumCumTp;
1709 // double avg_sumVHa;
1710 
1711  for(uint r2= 0; r2<regIds2.size();r2++){
1712  int regId = regIds2[r2];
1713  regPx = MTHREAD->MD->getRegion(regId)->getMyPixels();
1714 
1715  for (uint p=0;p<regPx.size();p++){
1716  Pixel* px = regPx[p];
1717  px->cumTp.clear();
1718  px->cumTp_exp.clear();
1719  px->vHa_exp.clear();
1720  px->vHa.clear();
1721  px->cumAlive.clear();
1722  px->cumAlive_exp.clear();
1723  double expType = px->expType;
1724 
1725  for(uint j=0;j<fTypes.size();j++){
1726  string ft = fTypes[j];
1727 
1728  double tp_multiplier_now = px->getMultiplier("tp_multiplier",ft,DATA_NOW);
1729  double tp_multiplier_t0 = px->getMultiplier("tp_multiplier",ft,firstYear);
1730  double mortCoef_multiplier_now = px->getMultiplier("mortCoef_multiplier",ft,DATA_NOW);
1731  double mortCoef_multiplier_t0 = px->getMultiplier("mortCoef_multiplier",ft,firstYear);
1732  double betaCoef_multiplier_now = px->getMultiplier("betaCoef_multiplier",ft,DATA_NOW);
1733  double betaCoef_multiplier_t0 = px->getMultiplier("betaCoef_multiplier",ft,firstYear);
1734  double pathMort_now, pathMort_t0;
1735 
1736  // calculating the cumulative time of passage and the (cumulativelly generated) vHa for each diameter class (depending on forest owners diam growth expectations)
1737  //loop(u$(ord(u)=1),
1738  // cumTp(u,i,lambda,essence) = tp_u1(i,essence,lambda);
1739  //);
1740  //loop(u$(ord(u)>1),
1741  // cumTp(u,i,lambda,essence) = cumTp(u-1,i,lambda,essence)+tp(u-1,i,lambda,essence);
1742  //);
1743  ////ceil(x) DNLP returns the smallest integer number greater than or equal to x
1744  //loop( (u,i,lambda,essence),
1745  // cumTp(u,i,lambda,essence) = ceil(cumTp(u,i,lambda,essence));
1746  //);
1747  vector <double> cumTp_temp; // cumulative time of passage to REACH a diameter class (tp is to LEAVE to the next one)
1748  vector <double> vHa_temp; // volume at hectar by each diameter class [m^3/ha]
1749  vector <double> cumAlive_temp; // cumulated alive rate to reach a given diameter class
1750  vector <double> cumTp_exp_temp; // expected version of cumTp_temp
1751  vector <double> vHa_exp_temp; // expected version of vHa_temp
1752  vector <double> cumAlive_exp_temp; // "expected" version of cumMort
1753 
1754  MD->setErrorLevel(MSG_NO_MSG); // as otherwise on 2007 otherwise sfd() will complain that is filling multiple years (2006 and 2007)
1755  for (uint u=0; u<dClasses.size(); u++){
1756  string dc = dClasses[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;
1761  pathMort_now = px->getPathMortality(ft,dc,DATA_NOW);
1762  pathMort_t0 = px->getPathMortality(ft,dc,firstYear);
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;
1767 
1768  // only cumTp is depending for the expectations, as it is what it is used by owner to calculate return of investments.
1769  // the tp, beta and mort coefficients instead are the "real" ones as predicted by scientist for that specific time
1770 
1771  if(u==0) {
1772  // first diameter class.. expected and real values are the same (0)
1773  cumTp_u = 0.;
1774  vHa_u = 0.;
1775  cumAlive_u = 1.;
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);
1782  } else {
1783  // other diameter classes.. first dealing with real values and then with expected ones..
1784  // real values..
1785  // real values..
1786  tp = gfd("tp",regId,ft,dClasses[u-1],thisYear)*tp_multiplier_now;
1787  cumTp_u = cumTp_temp[u-1] + tp;
1788  if (u==1){
1789  /**
1790  Note on the effect of mortality modifiers on the entryVolHa.
1791  Unfortunatly for how it is defined the mortality multiplier (the ratio with the new mortality rate over the old one) we can't
1792  compute a entryVolHa based on it. It is NOT infact just like: vHa_adjusted = vHa_orig / mort_multiplier.
1793  The effect of mortality on the vHa of the first diameter class is unknow, and so we can't compute the effect of a relative
1794  increase.
1795  */
1796  vHa_u = gfd("entryVolHa",regId,ft,"",thisYear);
1797  mort = 0.; // not info about mortality first diameter class ("00")
1798  } else {
1799  mort = 1-pow(1-min(gfd("mortCoef",regId,ft,dClasses[u-1],thisYear)*mortCoef_multiplier_now+pathMort_now+additionalMortality_now,1.0),tp); // mortality of the previous diameter class
1800  beta = gfd("betaCoef",regId,ft,dc, thisYear)*betaCoef_multiplier_now;
1801  vHa_u = vHa_temp[u-1]*beta*(1-mort);
1802  }
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);
1807  // expected values..
1808  /**
1809  param expType Specify how the forest owners (those that make the investments) behave will be the time of passage in the future in order to calculate the cumulative time of passage in turn used to discount future revenues.
1810  Will forest owners behave adaptively believing the time of passage between diameter classes will be like the observed one at time they make decision (0) or they will have full expectations believing forecasts (1) or something in the middle ?
1811  For compatibility with the GAMS code, a -1 value means using initial simulation tp values (fixed cumTp)."
1812  */
1813  if (expType == -1){
1814  tp_exp = gfd("tp",regId,ft,dClasses[u-1],firstYear)*tp_multiplier_t0;
1815  //tp = px->tp.at(u); no. not possible, tp stored at pixel level is the current year one
1816  cumTp_u_exp = cumTp_exp_temp[u-1]+tp_exp;
1817  cumTp_exp_temp.push_back(cumTp_u_exp);
1818  if(u==1) {
1819  vHa_u_exp = gfd("entryVolHa",regId,ft,"",firstYear);
1820  mort_exp = 0.; // not info about mortality first diameter class ("00")
1821  } else {
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); // mortality rate of previous diameter class
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);
1825  }
1826  } else {
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 ; // time of passage that there should be to reach this diameter class in the year where the previous diameter class will be reached
1831  cumTp_u_fullExp = cumTp_exp_temp[u-1]+tp_fullExp ; // it adds to the time of passage to reach the previous diameter class the time of passage that there should be to reach this diameter class in the year where the previous diameter class will be reached
1832  cumTp_u_exp = cumTp_u_fullExp*expType+cumTp_u_noExp*(1-expType); // 20121108: it's math the same as cumTp_exp_temp[u-1] + tp
1833  cumTp_exp_temp.push_back(cumTp_u_exp);
1834  if(u==1) {
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);
1838  mort_exp = 0.; // not info about mortality first diameter class ("00")
1839  } else {
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); // mortCoef is a yearly value. Mort coeff between class is 1-(1-mortCoeff)^tp
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); // mortality of the previous diameter class
1846  //double debug1 = gfd("mortCoef",regId,ft,dClasses[u-1],thisYear+ceil(cumTp_exp_temp[u-1]));
1847  //double debug2 = debug1*mortCoef_multiplier_dynamic+pathMort_dynamic;
1848  //double debug3 = min(1.0,debug2);
1849  //double debug4 = 1.0-debug3;
1850  //double debug5 = pow(debug4,tp_fullExp);
1851  //double debug6 = 1.0-debug5;
1852 
1853 
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); // BUG !!! mort is yearly value, not between diameter class. SOLVED 20121108
1860  }
1861  }
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);
1865 
1866  //cout << "********" << endl;
1867  //cout << "dc;mort;cumAlive;cumAlive_exp "<< endl ;
1868  //cout << dClasses[u] << ";"<< mort << ";" << cumAlive_u << ";" << cumAlive_exp_u << endl;
1869 
1870  }
1871  // debug stuff on vHa
1872  //double vHa_new = gfd("vHa",regId,ft,dc,DATA_NOW);
1873  //double hv2fa_old = gfd("hv2fa",regId,ft,dc,DATA_NOW);
1874  //cout << "Reg|Ft|dc|vHa (new)|1/hv2fa (old): " << regId << " | " << ft;
1875  //cout << " | " << dc << " | " << vHa_new << " | " << 1/hv2fa_old << endl;
1876 
1877  } // end of each diam
1878  //double pixID = px->getID();
1879  //cout << thisYear << ";"<< regIds2[r2] << ";" << pixID << ";" << ft << ";" << cumTp_exp_temp[3] << ";" << vHa_exp_temp[3] << endl;
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);
1885  px->cumAlive_exp.push_back(cumAlive_exp_temp);
1886 
1887  //sumCumTP += cumTp_exp_temp[3];
1888  //sumVHa += vHa_exp_temp[3];
1889  //count ++;
1890 
1891 
1892  } // end of each ft
1893  //double debug = 0.0;
1894  } // end of each pixel
1895  } // end of each region
1896  MD->setErrorLevel(MSG_ERROR);
1897  //avg_sumCumTp = sumCumTP/ count;
1898  //avg_sumVHa = sumVHa / count;
1899  //cout << "Avg sumCumTp_35 and sumVha_35: " << avg_sumCumTp << " and " << avg_sumVHa << " (" << count << ")" << endl;
1900  //exit(0);
1901 }
1902 
1903 void
1904 ModelCoreSpatial::resetPixelValues(){
1905  msgOut(MSG_INFO, "Starting resetting pixel level values");
1906  for(uint r2= 0; r2<regIds2.size();r2++){
1907  int regId = regIds2[r2];
1908  regPx = MTHREAD->MD->getRegion(regId)->getMyPixels();
1909  for (uint p=0;p<regPx.size();p++){
1910  Pixel* px = regPx[p];
1911  px->swap(VAR_VOL); // vol_l = vol
1912  px->swap(VAR_AREA); // area_l = area
1913  // 20121108 BUG! Solved, used empty (just return true if the vector is empty) instead of clear (it actually clears the vector)
1914  px->vol.clear(); // by ft,dc
1915  px->area = px->area_l; // ATTENCTION, DIFFERENT FROM THE OTHERS. Here it is not cleared, it is assigned the previous year as default
1916  /*px->area.clear(); // by ft,dc*/
1917  px->hArea.clear(); // by ft, dc
1918  px->deltaArea.clear();
1919  //px->regArea.clear(); // by year, ft NO, this one is a map, it doesn't need to be changed
1920  px->hVol.clear(); // by ft, dc
1921  px->hVol_byPrd.clear(); // by ft, dc, pp
1922  //px->in.clear(); // by pp
1923  //px->hr.clear(); // by pp
1924  px->vReg.clear(); // by ft
1925  px->expectedReturns.clear(); // by ft
1926  px->beta.clear();
1927  px->mort.clear();
1928  px->addMort.clear(); // fire mortality
1929  px->tp.clear();
1930  px->cumTp.clear();
1931  px->vHa.clear();
1932  px->cumTp_exp.clear();
1933  px->vHa_exp.clear();
1934  px->cumAlive.clear();
1935  px->cumAlive_exp.clear();
1936  px->vMort.clear();
1937  px->vMortAdd.clear(); // fire mortality
1938  //std::fill(rpx[j]->vMort.begin(), rpx[j]->vMort.end(), 0.0);
1939 
1940  }
1941  }
1942 }
1943 
1944 void
1945 ModelCoreSpatial::cachePixelExogenousData(){
1946 
1947  msgOut(MSG_INFO, "Starting cashing on pixel spatial-level exogenous data");
1948  bool applyAvalCoef = MTHREAD->MD->getBoolSetting("applyAvalCoef");
1949 
1950  for(uint r2= 0; r2<regIds2.size();r2++){
1951  int regId = regIds2[r2];
1952  regPx = MTHREAD->MD->getRegion(regId)->getMyPixels();
1953  for (uint p=0;p<regPx.size();p++){
1954  Pixel* px = regPx[p];
1955  px->tp.clear();
1956  px->beta.clear();
1957  px->mort.clear();
1958  px->addMort.clear(); // fire mortality
1959  px->avalCoef.clear();
1960 
1961  for(uint j=0;j<fTypes.size();j++){
1962  string ft = fTypes[j];
1963  vector <double> tp_byu;
1964  vector <double> beta_byu;
1965  vector <double> mort_byu;
1966  vector <double> addMort_byu; // fire mortality
1967 
1968  double tp_multiplier_now = px->getMultiplier("tp_multiplier",ft,DATA_NOW);
1969  double mortCoef_multiplier_now = px->getMultiplier("mortCoef_multiplier",ft,DATA_NOW);
1970 
1971  double betaCoef_multiplier_now = px->getMultiplier("betaCoef_multiplier",ft,DATA_NOW);
1972  double avalCoef_ft = applyAvalCoef?px->getSTData("avalCoef",ft,DATA_NOW):1.0;
1973 
1974  if(avalCoef_ft == MTHREAD->MD->getDoubleSetting("noValue")){
1975  msgOut(MSG_CRITICAL_ERROR,"applyAvalCoef has benn asked for, but I can't find the avalCoef in the data!");
1976  }
1977 
1978  for (uint u=0; u<dClasses.size(); u++){
1979  string dc = dClasses[u];
1980  double pathMortality = px->getPathMortality(ft,dc,DATA_NOW);
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; // these are volumes not rates !
1983  double tp, beta_real, mort_real;
1984  if (u==0){
1985  // tp of first diameter class not making it change across the time dimension, otherwise problems in getting the rigth past
1986  // regenerations. BUT good, px->tp.at(0) is used only to pick up the right regeneration, so the remaining of the model
1987  // uses the getMultiplier version and cumTp consider the dynamic effects also in the first dc.
1988  tp = gfd("tp",regId,ft,dClasses[u],firstYear)*px->getMultiplier("tp_multiplier",ft,firstYear); // tp is defined also in the first diameter class, as it is the years to reach the NEXT diameter class
1989  } else {
1990  tp = gfd("tp",regId,ft,dClasses[u],DATA_NOW)*tp_multiplier_now; // tp is defined also in the first diameter class, as it is the years to reach the NEXT diameter class
1991  }
1992  beta_real = u?gfd("betaCoef",regId,ft,dClasses[u],DATA_NOW)*betaCoef_multiplier_now:0;
1993  mort_real = min(u?gfd("mortCoef",regId,ft,dClasses[u],DATA_NOW)*mortCoef_multiplier_now+pathMortality :0,1.0); // 20230428 bugfix Antonello, removed additionalMortality //Antonello, bug fixed 20160203: In any case, natural plus pathogen mortality can not be larger than 1!
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); // fire mortality
1998  } // end of each dc
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); // fire mortality
2003  px->avalCoef.push_back(avalCoef_ft);
2004  } // end of each ft
2005  } // end of each pixel
2006  } // end of each region
2007 }
2008 
2009 void
2010 ModelCoreSpatial::computeInventary(){ // in=f(vol_t-1)
2011  msgOut(MSG_INFO, "Starting computing inventary available for this year..");
2012  int nbounds = pow(2,priProducts.size());
2013  vector<vector<int>> concernedPriProductsTotal = MTHREAD->MD->createCombinationsVector(priProducts.size());
2014  int currentYear = MTHREAD->SCD->getYear();
2015 
2016  for(uint i=0;i<regIds2.size();i++){
2017  int r2 = regIds2[i];
2018  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
2019  //Gis* GIS = MTHREAD->GIS;
2020  regPx = REG->getMyPixels();
2021  vector <double> in_reg(priProducts.size(),0.); // should have ceated a vector of size priProducts.size(), all filled with zeros
2022  vector <double> in_deathTimber_reg(priProducts.size(),0.); // should have ceated a vector of size priProducts.size(), all filled with zeros
2023  for (uint p=0;p<regPx.size();p++){
2024  Pixel* px = regPx[p];
2025  //int debugPx = px->getID();
2026  //int debug2 = debugPx;
2027  //px->in.clear();
2028  for(uint pp=0;pp<priProducts.size();pp++){
2029  double in = 0;
2030  for(uint ft=0;ft<fTypes.size();ft++){
2031  for(uint dc=0;dc<dClasses.size();dc++){
2032  in += app(priProducts[pp],fTypes[ft],dClasses[dc])*px->vol_l.at(ft).at(dc)*px->avalCoef.at(ft);
2033  }
2034  }
2035  //px->in.push_back(in);
2036  in_reg.at(pp) += in;
2037  } // end of each priProduct
2038  } // end each pixel
2039 
2040 
2041  for(uint pp=0;pp<priProducts.size();pp++){
2042  vector<string> priProducts_vector;
2043  priProducts_vector.push_back(priProducts[pp]);
2044 
2045  double in_deathMortality = MD->getAvailableDeathTimber(priProducts_vector,r2,currentYear-1);
2046  in_deathTimber_reg.at(pp) += in_deathMortality;
2047 
2048  // Even if I fixed all the lower bounds to zero in Opt::get_bounds_info still the model
2049  // doesn't solve with no-forest in a region.
2050  // Even with 0.0001 doesn't solve !!
2051  // With 0.001 some scenarios doesn't solve in 2093
2052  // With 0.003 vRegFixed doesn't solve in 2096
2053  // Tried with 0.2 but no changes, so put it back on 0.003
2054  //spd(max(0.001,in_reg.at(pp)),"in",r2,priProducts[pp],DATA_NOW,true);
2055  spd(in_reg.at(pp),"in",r2,priProducts[pp],DATA_NOW,true);
2056  spd(in_deathTimber_reg.at(pp),"in_deathTimber",r2,priProducts[pp],DATA_NOW,true);
2057  #ifdef QT_DEBUG
2058  if (in_reg.at(pp) < -0.0){
2059  msgOut(MSG_CRITICAL_ERROR,"Negative inventory");
2060  }
2061  #endif
2062  }
2063 
2064  // ##### Now creating a set of bonds for the optimisation that account of the fact that the same ft,dc can be used for multiple products:
2065 
2066  // 20160928: Solved a big bug: for each combination instead of taking the UNION of the various priProduct inventory sets I was taking the sum
2067  // Now both the alive and the death timber are made from the union
2068  // 20150116: As the same (ft,dc) can be used in more than one product knowing -and bounding the supply in the optimisation- each single
2069  // in(pp) is NOT enought.
2070  // We need to bound the supply for each possible combination, that is for 2^(number of prim.pr)
2071  // Here we compute the detailed inventory. TO.DO: Create the pounds in Opt. done
2072  // 20160209: Rewritten and corrected a bug that was not giving enought inv to multiproduct combinations
2073  for (uint i=0; i<nbounds; i++){
2074  vector<int> concernedPriProducts = concernedPriProductsTotal[i];
2075  vector<string> concernedPriProducts_ids = positionsToContent(priProducts,concernedPriProducts);
2076  //double debug = 0.0;
2077  //for(uint z=0;z<concernedPriProducts.size();z++){
2078  // debug += gpd("in",r2,priProducts[concernedPriProducts[z]]); // to.do: this will need to be rewritten checked!
2079  //}
2080  double bound_alive = MD->getAvailableAliveTimber(concernedPriProducts_ids,r2); // From px->vol_l, as in "in"
2081  double bound_deathTimber = MD->getAvailableDeathTimber(concernedPriProducts_ids,r2,currentYear-1); // From deathTimberInventory map
2082  double bound_total = bound_alive + bound_deathTimber;
2083 
2084  REG->inResByAnyCombination[i] = bound_total;
2085  //REG->inResByAnyCombination_deathTimber[i] = bound_deathTimber;
2086  } // end for each bond
2087  } // end each region
2088 }
2089 
2090 void
2091 ModelCoreSpatial::updateMapAreas(){
2092  msgOut(MSG_INFO, "Updating map areas..");
2093 
2094  if (!oldVol2AreaMethod){
2095  if(!MD->getBoolSetting("usePixelData")) return;
2096  for(uint i=0;i<regIds2.size();i++){
2097  ModelRegion* reg = MD->getRegion(regIds2[i]);
2098  vector <Pixel*> rpx = MTHREAD->GIS->getAllPlotsByRegion(regIds2[i]);
2099  for (uint p=0;p<rpx.size();p++){
2100  Pixel* px = rpx[p];
2101  double pxid= px->getID();
2102  for(uint j=0;j<fTypes.size();j++){
2103  string ft = fTypes[j];
2104  double forArea = vSum(px->area.at(j));
2105  #ifdef QT_DEBUG
2106  if(forArea < 0.0 ){
2107  msgOut(MSG_CRITICAL_ERROR, "Negative forArea in updateMapAreas");
2108  }
2109  #endif
2110  px->changeValue("forArea_"+ft, forArea*10000);
2111  } // end ft
2112  } // end px
2113  } // end region
2114  } else {
2115  int currentYear = MTHREAD->SCD->getYear();
2116  map<int,double> forestArea; // foresta area by each region
2117  pair<int,double > forestAreaPair;
2118  vector<int> l2Regions = MTHREAD->MD->getRegionIds(2, true);
2119  vector <string> fTypes = MTHREAD->MD->getForTypeIds();
2120  int nFTypes = fTypes.size();
2121  int nL2Regions = l2Regions.size();
2122  for(int i=0;i<nL2Regions;i++){
2123  int regId = l2Regions[i];
2124  vector<Pixel*> rpx = MTHREAD->GIS->getAllPlotsByRegion(regId);
2125  for(int j=0;j<nFTypes;j++){
2126  string ft = fTypes[j];
2127  //double regForArea = reg->getValue("forArea_"+ft);
2128  //double harvestedArea = gfd("harvestedArea",regId,ft,DIAM_ALL);
2129  //double regArea = gfd("regArea",regId,ft,DIAM_ALL);
2130  //cout << "Regid/ft/area/harvested/regeneration: " <<regId<<";"<<ft<<";"<<regForArea<<";"<<harvestedArea<<";" <<regArea<<endl;
2131  //double newAreaNet = regArea-harvestedArea;
2132  //double newAreaRatio = newAreaNet / regForArea;
2133  for(uint z=0;z<rpx.size();z++){
2134  Pixel* px = rpx[z];
2135  double oldValue = px->getDoubleValue("forArea_"+ft,true)/10000;
2136  double hArea = vSum(px->hArea.at(j)); //bug 20140205 areas in the model are in ha, in the layer in m^2
2137  double regArea = findMap(px->regArea,currentYear).at(j); //bug 20140205 areas in the model are in ha, in the layer in m^2
2138  //double newValue = oldValue*(1. + newAreaRatio);
2139  double newValue = oldValue-hArea+regArea;
2140  double areaNetOfRegeneration = oldValue-hArea;
2141  #ifdef QT_DEBUG
2142  if (areaNetOfRegeneration<0.0){
2143  msgOut(MSG_CRITICAL_ERROR,"areaNetOfRegeneration negative in updateMapAreas");
2144  }
2145  if (newValue<0.0){
2146  msgOut(MSG_CRITICAL_ERROR,"for area negative in updateMapAreas");
2147  }
2148  #endif
2149  rpx[z]->changeValue("forArea_"+ft, newValue*10000);
2150  }
2151  }
2152  }
2153  }
2154 }
2155 
2156 void
2157 ModelCoreSpatial::updateOtherMapData(){
2158 
2159 vector<int> l2Regions = MTHREAD->MD->getRegionIds(2, true);
2160 vector <string> fTypes = MTHREAD->MD->getForTypeIds();
2161 int nFTypes = fTypes.size();
2162 int nL2Regions = l2Regions.size();
2163 for(int i=0;i<nL2Regions;i++){
2164  int regId = l2Regions[i];
2165  vector<Pixel*> rpx = MTHREAD->GIS->getAllPlotsByRegion(regId);
2166  for(int j=0;j<nFTypes;j++){
2167  string ft = fTypes[j];
2168  for(uint z=0;z<rpx.size();z++){
2169  Pixel* px = rpx[z];
2170  double vol = vSum(px->vol.at(j));
2171  double expectedReturns = px->expectedReturns.at(j);
2172  if(MTHREAD->GIS->layerExist("vol_"+ft)){
2173  rpx[z]->changeValue("vol_"+ft, vol);
2174  }
2175  if(MTHREAD->GIS->layerExist("expectedReturns_"+ft)){
2176  rpx[z]->changeValue("expectedReturns_"+ft, expectedReturns);
2177  }
2178  }
2179  }
2180 }
2181 
2182 // update GUI image..
2183 for(int j=0;j<nFTypes;j++){
2184  string ft = fTypes[j];
2185  MTHREAD->GIS->updateImage("vol_"+ft);
2186  MTHREAD->GIS->updateImage("expectedReturns_"+ft);
2187 }
2188 
2189 
2190 }
2191 
2192 
2193 void
2194 ModelCoreSpatial::sumRegionalForData(){
2195 
2196  msgOut(MSG_INFO, "Summing data pixels->region..");
2197  //vector <string> outForVariables = MTHREAD->MD->getStringVectorSetting("outForVariables");
2198  int currentYear = MTHREAD->SCD->getYear();
2199 
2200  //map<vector<string>,double> hVol_byPrd; // by regid, ft, dc, pp
2201  map<tr1::array<string, 4>,double> hVol_byPrd; // by regid, ft, dc, pp
2202 
2203  // OLD CODE TO
2204  for(uint r2= 0; r2<regIds2.size();r2++){
2205  int regId = regIds2[r2];
2206  regPx = MTHREAD->MD->getRegion(regId)->getMyPixels(); // by regId, ft, dc, pp
2207 
2208  // vector < vector <vector <double> > >hVol_byPrd; // by ft, dc, pp
2209 
2210  for(uint j=0;j<fTypes.size();j++){
2211  string ft = fTypes[j];
2212 
2213  double regArea = 0.;
2214  double sumAreaByFt = 0.;
2215  double pxForAreaByFt = 0.;
2216  double vReg = 0.;
2217  double sumSumHProductivity = 0.;
2218  double sumHArea = 0.;
2219 
2220  for (uint u=0; u<dClasses.size(); u++){
2221  string dc = dClasses[u];
2222  double vol =0.;
2223  double hV = 0.;
2224  double hArea = 0.;
2225  double vMort = 0.;
2226  double vMortAdd = 0; // fire mortality
2227  double sumHProductivity = 0.; // productivity * area
2228  for (uint pi=0;pi<regPx.size();pi++){
2229  Pixel* px = regPx[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); // fire mortality
2235  sumHProductivity += (px->hProductivity.at(j).at(u) * px->hArea.at(j).at(u));
2236  if(MD->getBoolSetting("outDetailedHv",DATA_NOW)){
2237  for(uint pp=0;pp<priProducts.size();pp++){
2238  string pr = priProducts[pp];
2239  // vector <string> hVKey = {i2s(regId), ft, dc, pr};
2240  tr1::array<string, 4> hVKey = {i2s(regId), ft, dc, pr};
2241  //double debug = px->hVol_byPrd.at(j).at(u).at(pp);
2242  //cout << px->getID() << '\t' << j << '\t' << u << '\t' << pp << '\t' << debug << endl;
2243  incrOrAddMapValue(hVol_byPrd, hVKey, px->hVol_byPrd.at(j).at(u).at(pp));
2244  }
2245  }
2246  }
2247  if(u){
2248  sfd(vol,"vol",regId,ft,dc,DATA_NOW);
2249  sfd(hV,"hV",regId,ft,dc,DATA_NOW,true);
2250  double hProductivityByDc = hArea ? sumHProductivity/hArea : 0.;
2251  sumSumHProductivity += (hProductivityByDc*hArea);
2252  sumHArea += hArea;
2253  sfd(hProductivityByDc,"hProductivityByDc",regId,ft,dc,DATA_NOW,true);
2254  sfd(hArea,"harvestedArea",regId,ft,dc,DATA_NOW, true);
2255  sfd(vMort,"vMort",regId,ft,dc,DATA_NOW,true);
2256  sfd(vMortAdd,"vMortAdd",regId,ft,dc,DATA_NOW,true); // fire mortality
2257  double vol_1 = gfd("vol",regId,ft,dc,currentYear-1);
2258  if(vol_1){
2259  sfd(hV/vol_1,"hr",regId,ft,dc,DATA_NOW, true);
2260  } else {
2261  sfd(0.,"hr",regId,ft,dc,DATA_NOW, true);
2262  }
2263 
2264  }
2265  } // end foreach dc
2266  double hProductivity = sumHArea? sumSumHProductivity / sumHArea : 0.;
2267  sfd(hProductivity,"hProductivity",regId,ft,"",DATA_NOW,true);
2268 
2269  for (uint p=0;p<regPx.size();p++){
2270  Pixel* px = regPx[p];
2271  vReg += px->vReg.at(j);
2272  regArea += findMap(px->regArea,currentYear).at(j);
2273  pxForAreaByFt = (px->getDoubleValue("forArea_"+ft,true)/10000);
2274 
2275  sumAreaByFt += pxForAreaByFt;
2276  //double debug1 = sumAreaByFt;
2277  if(! (sumAreaByFt >= 0.0) ){
2278  msgOut(MSG_CRITICAL_ERROR,"sumAreaByFt is not non negative.");
2279  }
2280  }
2281  sfd(vReg,"vReg",regId,ft,"",DATA_NOW, true);
2282  sfd(regArea,"regArea",regId,ft,"",DATA_NOW, true);
2283  sfd(sumAreaByFt,"forArea",regId,ft,"",DATA_NOW, true);
2284  } // end of for each ft
2285 
2286  for (uint p=0;p<regPx.size();p++){
2287  Pixel* px = regPx[p];
2288  double totPxForArea = vSum(px->area);
2289 
2290 #ifdef QT_DEBUG
2291  double totPxForArea_debug = 0.0;
2292  for(uint j=0;j<fTypes.size();j++){
2293  string ft = fTypes[j];
2294  totPxForArea_debug += (px->getDoubleValue("forArea_"+ft,true)/10000);
2295  }
2296 
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!");
2300  }
2301 #endif
2302  } // end of each pixel
2303 
2304 
2305  } // end each region
2306  //cout << "done" << endl;
2307  MTHREAD->DO->printDetailedHV(hVol_byPrd);
2308 
2309 
2310  // Taking care of expected returns here..
2311  // (Changed 25/08/2016 afternoon: expRet{ft,r} are now sum{px}{expRet{ft,px}*fArea_{px}}/fArea{r} and no longer sum{px}{expRet{ft,px}*fArea_{px,ft}}/fArea{r,ft} )
2312  // Also now we report the expReturns by group and by forest, each of which is made only with the best ones within their group
2313 
2314  vector<string> parentFtypes = MTHREAD->MD->getForTypeParents();
2315 
2316  for(uint r2= 0; r2<regIds2.size();r2++){
2317  int regId = regIds2[r2];
2318  regPx = MTHREAD->MD->getRegion(regId)->getMyPixels();
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);
2323 
2324  // First computing the sumExpectedReturns..
2325  for (uint p=0;p<regPx.size();p++){
2326  Pixel* px = regPx[p];
2327  //int debug_pxid = px->getID();
2328  double pxForArea = vSum(px->area);
2329  totRegionForArea += pxForArea;
2330  double bestPxExpectedRet = getMax(px->expectedReturnsNotCorrByRa);
2331  for(uint i=0;i<parentFtypes.size();i++){
2332  vector <string> childIds = MTHREAD->MD->getForTypeChilds(parentFtypes[i]);
2333  vector <int> childPos = MTHREAD->MD->getForTypeChilds_pos(parentFtypes[i]);
2334  vector<double> pxExpReturnsByChilds(childPos.size(),0.0);
2335  for(uint j=0;j<childPos.size();j++){
2336  double pxExpReturn_singleFt = px->expectedReturns.at(childPos[j]);
2337  // Manual fix to not have the expected returns of ash within the general "broadL" expected returns.
2338  // To do: remove it after we work on the ash project.. I don't like manual fixes !!!
2339  pxExpReturnsByChilds.at(j) = (childIds.at(j) == "ash") ? 0.0 : pxExpReturn_singleFt;
2340  //pxExpReturnsByChilds.at(j) = pxExpReturn_singleFt;
2341  totSumExpRet_byFTypes.at(childPos[j]) += pxExpReturn_singleFt*pxForArea;
2342  } // end of each ft
2343  totSumExpRet_byFTParent[i] += getMax(pxExpReturnsByChilds)*pxForArea;
2344  } // end for each partentFt
2345  totSumExpRet += bestPxExpectedRet * pxForArea;
2346  } // end for each px
2347 
2348  // ..and now computing the expReturns and storing them
2349  for(uint i=0;i<parentFtypes.size();i++){
2350  vector <int> childPos = MTHREAD->MD->getForTypeChilds_pos(parentFtypes[i]);
2351  for(uint j=0;j<childPos.size();j++){
2352  //double debug1 = totSumExpRet_byFTypes.at(childPos[j])/totRegionForArea;
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);
2355  } // end of each ft
2356  //double debug2 = totSumExpRet_byFTParent.at(i)/totRegionForArea;
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);
2359 
2360  } // end for each partentFt
2361  //double debug3 = totSumExpRet/totRegionForArea;
2362  sfd(totSumExpRet,"sumExpReturns",regId,"","",DATA_NOW, true);
2363  sfd(totSumExpRet/totRegionForArea,"expReturns",regId,"","",DATA_NOW, true);
2364 
2365  } // end for each region
2366 
2367  // Computing pathogens share of forest invasion
2368  if(MD->getBoolSetting("usePathogenModule")){
2369  for(uint r2= 0; r2<regIds2.size();r2++){
2370  int regId = regIds2[r2];
2371  regPx = MTHREAD->MD->getRegion(regId)->getMyPixels();
2372  double totalForArea = 0.0;
2373  double invadedArea = 0.0;
2374  for (uint p=0;p<regPx.size();p++){
2375  Pixel* px = regPx[p];
2376  int invaded = 0.0;
2377  for(uint j=0;j<fTypes.size();j++){
2378  for (uint u=0; u<dClasses.size(); u++){
2379  if(px->getPathMortality(fTypes[j],dClasses[u]) > 0){
2380  invaded = 1.0;
2381  }
2382  }
2383  }
2384  totalForArea += vSum(px->area);
2385  invadedArea += vSum(px->area)*invaded;
2386  }
2387  sfd(invadedArea/totalForArea,"totalShareInvadedArea",regId,"","",DATA_NOW, true);
2388  }
2389  } // end we are using path model
2390 
2391 
2392 
2393 
2394 }
2395 /**
2396  * This function call registerHarvesting() (accounts for emissions from for. operations), registerDeathBiomass() (registers new stocks of death biomass),
2397  * registerProducts() (registers new stock of products) and registerTransports() (accounts for emissions from transportation).
2398  *
2399  * It pass to registerProducts():
2400  * - for primary products, the primary products exported out of the country, but not those exported to other regions or used in the region as
2401  * these are assumed to be totally transformed to secondary products;
2402  * - for secondary products, those produced in the region from locally or regionally imported primary product plus those secondary products
2403  * imported from other regions, less those exported to other regions. It doesn't include the secondary products imported from abroad the country.
2404  */
2405 void
2406 ModelCoreSpatial::registerCarbonEvents(){
2407 
2408  //void registerHarvesting(const int & regId, const string & fType, const double & value); ///< register the harvesting of trees -> cumEmittedForOper
2409  //void registerDeathBiomass(const double &value, const int & regId, const string &fType);
2410  //void registerProducts(const double &value, const int & regId, const string &productName);
2411  //void registerTransports(const double &distQ, const int & regId);
2412 
2413  for(uint i=0;i<regIds2.size();i++){
2414  for(uint j=0;j<fTypes.size();j++){
2415  double deathBiomass = gfd("vMort",regIds2[i],fTypes[j],DIAM_ALL,DATA_NOW)+ gfd("vMortAdd",regIds2[i],fTypes[j],DIAM_ALL,DATA_NOW);
2416  double harvesting = gfd("hV",regIds2[i],fTypes[j],DIAM_ALL,DATA_NOW);
2417  MTHREAD->CBAL->registerDeathBiomass(deathBiomass, regIds2[i], fTypes[j]); // register new stock
2418  MTHREAD->CBAL->registerHarvesting(harvesting, regIds2[i], fTypes[j]); // account for emissions. Added 201500715: it also moves the extra biomass to the death biomass pool
2419  }
2420 
2421  for(uint p=0;p<priProducts.size();p++){
2422  // for the primary products we consider only the exports as the domestic consumption is entirelly transformed in secondary products
2423  double int_exports = gpd("sa",regIds2[i],priProducts[p],DATA_NOW);
2424  MTHREAD->CBAL->registerProducts(int_exports, regIds2[i], priProducts[p]); // register new stock
2425  }
2426  for(uint p=0;p<secProducts.size();p++){
2427  // for the tranformed product we skip those that are imported, hence derived from other forest systems
2428  // but we consider those coming from other regions
2429  double consumption = gpd("dl",regIds2[i],secProducts[p],DATA_NOW); // dl = sl + net regional imports
2430  MTHREAD->CBAL->registerProducts(consumption, regIds2[i], secProducts[p]); // register new stock
2431  }
2432 
2433  }
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];
2440  double dist = gpd("dist",rfrom,"",DATA_NOW,i2s(rto)); //km
2441  for(uint p=0;p<allProducts.size();p++){
2442  distQProd += dist*gpd("rt",rfrom,allProducts[p],DATA_NOW,i2s(rto)); //km*Mm^3
2443  }
2444  }
2445  MTHREAD->CBAL->registerTransports(distQProd, rfrom);
2446  }
2447  }
2448  MTHREAD->CBAL->HWP_eol2energy(); // used to compute the energy substitution from hwp that reach the end of life and doesn't go to landfil. Previously the energy substitution was computed in registerProducts(), that is at the time when the product was produced.
2449 
2450 }
2451 
2452 /**
2453  * Compute the expectation weighted price based on the ratio of the international (world) price between the future and now.
2454  *
2455  * @param curLocPrice The local current price
2456  * @param worldCurPrice The world current price
2457  * @param worldFutPrice The world future price
2458  * @param sl Supply local
2459  * @param sa Supply abroad
2460  * @param expCoef The expectation coefficient for prices for the agent [0,1]
2461  * @return The expType-averaged local (or weighter) price
2462  */
2463 double
2464 ModelCoreSpatial::computeExpectedPrice(const double & curLocPrice, const double & worldCurPrice, const double & worldFutPrice, const double & sl, const double & sa, const double & expCoef){
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;
2468 }
2469 
2470 /**
2471  * It uses volumes from gis data to "move" volumes from one forest type to the other (when called with what="vol"). Then it moves areas
2472  * proportionally and, as dc0 volumes are not defined but area it is, compute, again proportionally, area in destination forest times for dc=0
2473  * It acts on the pix->vol, pix->area and pix->area_l vectors. It also create/update the px->values layer map for the area, but it doesn't cash the
2474  * results in forDataMap.
2475  *
2476  * It is called first with parameter what="vol" in initializePixelVolumes() and then with what="area" in initializePixelAreas().
2477  * As we need the original volumes in the area allocation, original_vols is set as a static variable.
2478  *
2479  */
2480 void
2481 ModelCoreSpatial::loadExogenousForestLayers(const string & what){
2482  if(!MD->getBoolSetting("useSpExplicitForestTypes")) return;
2483 
2484  int nFTypes = fTypes.size();
2485  int nDC = dClasses.size();
2486  int pxC = 0;
2487 
2488  for(uint ir=0;ir<regIds2.size();ir++){
2489  int r2 = regIds2[ir];
2490  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
2491  regPx = REG->getMyPixels();
2492  pxC += regPx.size();
2493  }
2494 
2495  static vector<vector<vector<double>>> original_vols(pxC, vector<vector<double>>(nFTypes, vector<double>(nDC, 0.0))); // by px counter, ftype, dc
2496 
2497  if(what=="vol"){
2498  // first, before transfering volumes, saving the original ones..
2499  for(uint i=0;i<fTypes.size();i++){
2500  for (uint u=0; u<dClasses.size(); u++){
2501  int pxC_loc = 0;
2502  for(uint ir=0;ir<regIds2.size();ir++){
2503  int r2 = regIds2[ir];
2504  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
2505  regPx = REG->getMyPixels();
2506  for (uint p=0;p<regPx.size();p++){
2507  Pixel* px = regPx[p];
2508  original_vols[pxC_loc][i][u] += px->vol[i][u];
2509  pxC_loc ++;
2510  }
2511  }
2512  }
2513  }
2514  for(uint i=0;i<fTypes.size();i++){
2515  string fti = fTypes[i];
2516  for(uint o=0;o<fTypes.size();o++){
2517  string fto = fTypes[o];
2518  for (uint u=1; u<dClasses.size(); u++){ // first diameter class volumes are computed from the model..
2519  string layerName = "spInput#vol#"+fto+"#"+fti+"#"+i2s(u);
2520  if (MTHREAD->GIS->layerExist(layerName)){
2521  for(uint ir=0;ir<regIds2.size();ir++){
2522  int r2 = regIds2[ir];
2523  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
2524  regPx = REG->getMyPixels();
2525  for (uint p=0;p<regPx.size();p++){
2526  Pixel* px = regPx[p];
2527  double vol_transfer = min(px->getDoubleValue(layerName,true)/1000000,px->vol[i][u]) ; // Vol in the layer are in m^3, in the model in Mm^3
2528  px->vol[i][u] -= vol_transfer;
2529  px->vol[o][u] += vol_transfer;
2530  }
2531  }
2532  }
2533  }
2534  }
2535  }
2536  }
2537 
2538  if(what=="area"){
2539  /**
2540  Allocate area proportionally to volumes (see file test_proportional_computation_of_areas_from_volumes.ods)
2541  Example:
2542  FtIn FtOut Vtrasfer
2543  con ash 0.2
2544  brHf ash 0.1
2545  brCopp ash 0.3
2546  con oak 0.3
2547  brHf oak 0.2
2548  brCopp oak 0.1
2549 
2550  Vorig Aorig Vnew Anew
2551  con 10 30 9.5 28.5 Aorig-Aorig*(Vtrasfer1/Vorig)-Aorig(Vtrasfer2/Vorig)
2552  brHf 5 20 4.7 18.8
2553  brCopp 2 20 1.6 16
2554  ash 0 0 0.6 4 Aorig1*Vtrasfer1/(Vorig1)+Aorig2*Vtrasfer2/(Vorig2)+...
2555  oak 0 0 0.6 2.7
2556  70 70
2557  */
2558  // first, before transfering areas, saving the original ones (we already saved the vols in the what="vol" section, that is called before this one)..
2559  vector<vector<vector<double>>> original_areas(pxC, vector<vector<double>>(nFTypes, vector<double>(nDC, 0.0))); // by px counter, ftype, dc
2560  for(uint i=0;i<fTypes.size();i++){
2561  for (uint u=0; u<dClasses.size(); u++){
2562  int pxC_loc = 0;
2563  for(uint ir=0;ir<regIds2.size();ir++){
2564  int r2 = regIds2[ir];
2565  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
2566  regPx = REG->getMyPixels();
2567  for (uint p=0;p<regPx.size();p++){
2568  Pixel* px = regPx[p];
2569  original_areas[pxC_loc][i][u] += px->area_l[i][u];
2570  pxC_loc ++;
2571  }
2572  }
2573  }
2574  }
2575 
2576 
2577  // transferred areas ordered by pxcounter, i and then o ftype. Used to then repart the 0 diameter class..
2578  vector<vector<vector<double>>> transferred_areas(pxC, vector<vector<double>>(nFTypes, vector<double>(nFTypes, 0.0))); // initialize a 3d vector of nFTypes zeros.
2579 
2580  for(uint i=0;i<fTypes.size();i++){
2581  string fti = fTypes[i];
2582  for(uint o=0;o<fTypes.size();o++){
2583  string fto = fTypes[o];
2584  for (uint u=1; u<dClasses.size(); u++){ // first diameter class area is comuted proportionally..
2585  string layerName = "spInput#vol#"+fto+"#"+fti+"#"+i2s(u);
2586  if (MTHREAD->GIS->layerExist(layerName)){
2587  int pxC_loc = 0;
2588  for(uint ir=0;ir<regIds2.size();ir++){
2589  int r2 = regIds2[ir];
2590  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
2591  regPx = REG->getMyPixels();
2592  for (uint p=0;p<regPx.size();p++){
2593  Pixel* px = regPx[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; // Vol in the layer are in m^3, in the model in Mm^3
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;
2599  px->area[i][u] = px->area_l[i][u];
2600  px->area_l[o][u] += area_transfer;
2601  px->area[o][u] = px->area_l[o][u];
2602  transferred_areas[pxC_loc][i][o] += area_transfer;
2603  pxC_loc ++;
2604  }
2605  }
2606  }
2607  }
2608  }
2609  }
2610 
2611  // Moving the area in the 0 diameter class, for which no info is normally available..
2612  double pxC_loc = 0;
2613  for(uint ir=0;ir<regIds2.size();ir++){
2614  int r2 = regIds2[ir];
2615  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
2616  regPx = REG->getMyPixels();
2617  for (uint p=0;p<regPx.size();p++){
2618  Pixel* px = regPx[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++){ // we want to skip the 0 diameter class, this is why we don't simply use vSum()..
2623  area_i_orig += original_areas[pxC_loc][i][u];
2624  }
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 ;
2627  px->area[i][0] = px->area_l[i][0];
2628  px->area_l[o][0] += area_transfer_u0 ; // bug corrected 20151130: it was 0 instead of o (output) !!
2629  px->area[o][0] = px->area_l[0][0]; // bug corrected 20151130: it was 0 instead of o (output) !!
2630  }
2631  }
2632  pxC_loc++;
2633  }
2634  }
2635 
2636  // Alligning the area memorised in the px layers to the new areas of the ft..
2637  for(uint i=0;i<fTypes.size();i++){
2638  string fti_id = fTypes[i];
2639  forType* fti = MTHREAD->MD->getForType(fti_id);
2640  int ft_memType = fti->memType;
2641  string ft_layerName = fti->forLayer;
2642  //if(ft_memType==3){
2643  // MTHREAD->GIS->addLayer(ft_layerName,ft_layerName,false,true); //20151130: no needed as we already added it in applyForestReclassification (yes, odd, as memory type 3 layerrs do not have any reclassification rule associated, but if I don't add the layer at that time I got other errors)
2644  // }
2645  if(ft_memType==3 ||ft_memType==2){
2646  for(uint ir=0;ir<regIds2.size();ir++){
2647  int r2 = regIds2[ir];
2648  ModelRegion* REG = MTHREAD->MD->getRegion(r2);
2649  regPx = REG->getMyPixels();
2650  for (uint p=0;p<regPx.size();p++){
2651  Pixel* px = regPx[p];
2652  double area_px = vSum(px->area[i]);
2653  px->changeValue(ft_layerName,area_px*10000);
2654  }
2655  }
2656  }
2657  }
2658  } // end if what is area
2659 }
2660 
2661 void
2662 ModelCoreSpatial::printDebugInitRegionalValues(){
2663  // Print debug stats on inventory and supplies in each region..
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++){
2668  for(uint p=0;p<priProducts.size();p++){
2669  int r2 = l2r[r1][r2c];
2670  double inv = gpd("in",r2,priProducts[p],secondYear);
2671  double st = gpd("st",r2,priProducts[p],secondYear);
2672  double sl = gpd("sl",r2,priProducts[p],secondYear);
2673  double sa = gpd("sa",r2,priProducts[p],secondYear);
2674  cout << r2 << "\t" << priProducts[p] << "\t\t" << inv << "\t" << st << "\t" << sl << "\t" << sa << endl;
2675  }
2676  }
2677  } // end of r1 region
2678  exit(0);
2679 
2680 }
2681 
2682 /**
2683  * @brief ModelCoreSpatial::allocateHarvesting
2684  * @param total_st vector of total supply by primary products
2685  * @return a vector of the remaining supply that goes allocated to alive timber (that is, to harvesting)
2686  *
2687  * The algorithm is such that is loops the deathTimberInventory map for each year (newer to older), dc (higher to smaller) and ft.
2688  * It compute the primary products allocable from that combination and allocate the cell amount to decrease the total_st of that products
2689  * in a proportional way to what still remain of the allocable products.
2690  *
2691  * It is called in the runMarketModule() function.
2692  *
2693  */
2694 
2695 vector <double>
2696 ModelCoreSpatial::allocateHarvesting(vector<double> total_st, const int &regId){
2697  if(!MD->getBoolSetting("useDeathTimber",DATA_NOW)) return total_st;
2698  //vector <double> st_fromHarv(priProducts.size(),0.);
2699  //map<iisskey, double > * deathTimberInventory= MD->getDeathTimberInventory();
2700  int maxYears = MD->getMaxYearUsableDeathTimber();
2701  int currentYear = MTHREAD->SCD->getYear();
2702  for(uint y = currentYear-1; y>currentYear-1-maxYears; y--){
2703  for (int u = dClasses.size()-1; u>=0; u--){ // I need to specify u as an integer !
2704  string dc = dClasses.at(u);
2705  for (uint f=0; f<fTypes.size(); f++){
2706  string ft = fTypes[f];
2707  vector<int>allocableProducts = MD->getAllocableProductIdsFromDeathTimber(regId, ft, dc, y, currentYear-1);
2708  iisskey key(y,regId,ft,dc);
2709  double deathTimber = MD->deathTimberInventory_get(key);
2710  double sum_total_st_allocable = 0;
2711  // Computing shares/weights or remaining st to allcate
2712  for(uint ap=0;ap<allocableProducts.size();ap++){
2713  sum_total_st_allocable += total_st.at(allocableProducts[ap]);
2714  }
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);
2718  MD->deathTimberInventory_incrOrAdd(key,-allocated);
2719  total_st[allocableProducts[ap]] -= allocated;
2720  }
2721  }
2722  }
2723  }
2724  return total_st;
2725 }
2726 
2727 
2728 /**
2729  * This function has two tasks:
2730  * 1) compute the policy balances (- -> public cost,subside; + -> public revenue, tax) for pol_mktDirInt_s, pol_mktDirInt_d, pol_trSub and pol_tcSub
2731  * (pol_fiSub are done directly in the runManagementModule() function)
2732  * 2) compute the producer/consumer surpluses
2733  */
2734 void
2735 ModelCoreSpatial::computeEconomicBalances(){
2736  //int thisYear = MTHREAD->SCD->getYear();
2737  int previousYear = MTHREAD->SCD->getYear()-1;
2738 
2739  for(uint r2=0; r2<regIds2.size();r2++){
2740  int regId = regIds2[r2];
2741  ModelRegion* reg = MD->getRegion(regId);
2742 
2743  for(uint p=0;p<priProducts.size();p++){
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);
2746  double sigma = gpd("sigma",regId,priProducts[p]);
2747  double in = gpd("in",regId,priProducts[p]);
2748  double in_1 = gpd("in",regId,priProducts[p], previousYear);
2749  double pc = gpd("pc",regId,priProducts[p]);
2750  double pc_1 = gpd("pc",regId,priProducts[p], previousYear);
2751  double sc = gpd("sc",regId,priProducts[p]);
2752  double sc_1 = gpd("sc",regId,priProducts[p],previousYear);
2753  double gamma_incr = gpd("gamma_incr",regId,priProducts[p]); // elast supply to increasing stocks
2754  double gamma_decr = gpd("gamma_decr",regId,priProducts[p]); // elast supply to decreasing stocks
2755 
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 -
2759  pc_1
2760  * pol_mktDirInt_s_1
2761  * pow(pol_mktDirInt_s,-1)
2762  * pow(sc_1,-1/sigma)
2763  * pow(in_1/in,gamma/sigma)
2764  * (sigma/(sigma+1))
2765  * pow(sc,(sigma+1)/sigma);
2766 
2767  spd(polBal_mktDirInt_s,"polBal_mktDirInt_s",regId,priProducts[p],DATA_NOW,true); // M€
2768  spd(surplus_prod,"surplus_prod",regId,priProducts[p],DATA_NOW,true);
2769  }
2770  for(uint p=0;p<secProducts.size();p++){
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);
2773  double sigma = gpd("sigma",regId,secProducts[p]);
2774  double pol_trSub = gpd("pol_trSub",regId,secProducts[p]);
2775  double pc = gpd("pc",regId,secProducts[p]);
2776  double pc_1 = gpd("pc",regId,secProducts[p],previousYear);
2777  double dc = gpd("dc",regId,secProducts[p]);
2778  double dc_1 = gpd("dc",regId,secProducts[p],previousYear);
2779  double dc0 = gpd("dc",regId,secProducts[p],secondYear);
2780  double m = gpd("m",regId,secProducts[p]);
2781  double sl = gpd("sl",regId,secProducts[p]); // sl because dl include also the production from byproduct that has not subsides in the model
2782 
2783  double polBal_mktDirInt_d = pc * (pol_mktDirInt_d - 1) * dc;
2784  double polBal_trSub = m * (pol_trSub-1) * sl;
2785  double trs = 0.5; // 0.1: we consider 90 percent of the demand integral
2786  double surplus_cons = pc_1
2787  * pol_mktDirInt_d_1
2788  * pow(dc_1,-1/sigma)
2789  * pow(pol_mktDirInt_d,-1)
2790  * (sigma/(sigma+1))
2791  * (pow(dc,(sigma+1)/sigma) - pow(trs*dc0,(sigma+1)/sigma))
2792  - (dc - trs * dc0) * pc ;
2793 
2794  spd(polBal_mktDirInt_d,"polBal_mktDirInt_d",regId,secProducts[p],DATA_NOW,true); // M€
2795  spd(polBal_trSub,"polBal_trSub",regId,secProducts[p],DATA_NOW,true); // M€
2796  spd(surplus_cons,"surplus_cons",regId,secProducts[p],DATA_NOW,true); // M€
2797  }
2798 
2799  for(uint p=0;p<allProducts.size();p++){
2800  double pol_tcSub = gpd("pol_tcSub",regId,allProducts[p]);
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();
2805  double ct = gpd("ct",regId,allProducts[p],DATA_NOW,i2s(regId2));
2806  double rt = gpd("rt",regId,allProducts[p],DATA_NOW,i2s(regId2));
2807  polBal_tcSub += ct*(pol_tcSub-1) * rt;
2808  }
2809  spd(polBal_tcSub,"polBal_tcSub",regId,allProducts[p],DATA_NOW,true); // M€
2810  }
2811 
2812  }
2813 
2814 }
2815 
2816 /**
2817  * Return the average age using cumTp vector
2818  */
2819 double
2820 ModelCoreSpatial::getAvgAgeByDc(Pixel* px, int ft, int dc){
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; }
2823 }
2824 
2825 
2826 
2827 
2828 // added for carbon project
2829 
2830 double
2831 ModelCoreSpatial::getVHaByYear(const Pixel* px, const int& ft, const int& year, const double& extraBiomass_ratio, const int& regId ) const {
2832  // This function computes and send back the expected volume/ha in a forest type after a specific amount of years
2833  // It uses expected volumes/ha at each diameter class, and makes a linear interpolation between the corresponding times of passage
2834  // Ancient version was a "stairs" function of time, now it is a succession of linear approximations
2835 
2836  double vHa = 0;
2837  int thisYear = MTHREAD->SCD->getYear();
2838 
2839  // This loop goes through all diameter classes to find the last one that was reached at year tested
2840  for (int u = 0; u<dClasses.size(); u++){
2841  // We retrieve the expected time of passage for the current diameter class
2842  // and we compare it to the year tested
2843  double cumTP = px->cumTp_exp[ft][u];
2844  if (year>=cumTP) {
2845  // Forest has already reached this class,
2846  // We do nothing and go to the next one
2847  }
2848  else {// Forest has not reached this class yet, we stop here
2849  // Assumption : Volume growth is a linear approximation between volume at previous DC and at current DC
2850  // There would be a problem if "year" was negative, but it should not happen since year starts at 0
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);
2856 
2857  vHa = (volHa_minus + (year-cumTP_minus)*slope)*(1+extraBiomass_ratio);
2858  return vHa;
2859  }
2860  }
2861  // if we arrive here, we have exited the loop
2862  // and the forest has reached the last diameter class
2863  int dc = dClasses.size()-1; // index of the last diameter class
2864  double mortCoeff = gfd("mort",regId,fTypes[ft],dClasses[dc],thisYear); // mortality rate of the last diameter class
2865  double cumTP_max = px-> cumTp_exp[ft][dc]; // cumulative TP to reach last diameter class
2866  double time_elapsed = year - cumTP_max; // time since the last class was reached
2867  double volHa_max = px->cumTp_exp[ft][dc]; // volume/ha when the last class was reached
2868 
2869  // Now we compute the volume/ha at the year tested
2870  // Assumption: trees don't grow anymore but slowly die according to mortality rate
2871  vHa = volHa_max*pow(1-mortCoeff,time_elapsed)*(1+extraBiomass_ratio);
2872  return vHa;
2873 } // end of function
2874 
2875 
2876 
2877 /** ANCIENT VERSION OF FUNCTION. VOLUME WAS A STAIRS FUNCTION OF TIME
2878  * Return the Volume/ha in a forest after a given number of year after planting, for a specific forest type
2879 
2880 double
2881 ModelCoreSpatial::getVHaByYear(const Pixel* px, const int& ft, const int& year, const double& extraBiomass_ratio) const {
2882  double vHa = 0;
2883  //int regId = px->getMyRegion()->getRegId();
2884  //double extraBiomass_ratio = gfd("extraBiomass_ratio",regId,fTypes[ft],"",DATA_NOW);
2885 
2886  for (int u = 0; u<dClasses.size(); u++){
2887  double cumTP_exp = px->cumTp_exp[ft][u];
2888  if (year>=cumTP_exp){
2889  vHa = px->vHa_exp[ft][u]*(1+extraBiomass_ratio);
2890  } else {
2891  return vHa;
2892  }
2893  }
2894  return vHa;
2895 }
2896 */
Gis::updateImage
void updateImage(string layerName_h)
Add one layer to the system.
Definition: Gis.cpp:700
Pixel::deltaArea
vector< vector< double > > deltaArea
Definition: Pixel.h:108
Pixel::setSpModifier
void setSpModifier(const double &value, const int &ftindex)
Definition: Pixel.h:87
iisskey
Class to provide a simple integer-integer-string-string key in std maps.
Definition: BaseClass.h:213
Scheduler::advanceYear
void advanceYear()
Definition: Scheduler.h:51
BaseClass::getPos
int getPos(const K &element, const vector< K > &v, const int &msgCode_h=MSG_CRITICAL_ERROR)
Definition: BaseClass.h:421
MSG_ERROR
Print an ERROR message, but don&#39;t stop the model.
Definition: BaseClass.h:61
DATA_NOW
The required data is for the current year.
Definition: BaseClass.h:73
ModelCoreSpatial::runBiologicalModule
void runBiologicalModule()
computes hV, hArea and new vol at end of year
Definition: ModelCoreSpatial.cpp:562
ModelData::getAllocableProductIdsFromDeathTimber
vector< int > getAllocableProductIdsFromDeathTimber(const int &regId_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...
Definition: ModelData.cpp:2105
ModelRegion::getArea
double getArea(const string &fType_h, const string &dClass_h)
Get area by ft and dc (from pixel->area matrix)
Definition: ModelRegion.cpp:125
Carbon::initialiseDeathBiomassStocks
void initialiseDeathBiomassStocks(const vector< double > &deathByFt, const int &regId)
Initialises the stocks of death biomass for the avgLive_* years before the simulation starts...
Definition: Carbon.cpp:169
ModelData::getIntSetting
int getIntSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1105
ModelCoreSpatial::l2r
vector< vector< int > > l2r
Definition: ModelCoreSpatial.h:142
Pixel::area
vector< vector< double > > area
Definition: Pixel.h:107
ModelCoreSpatial::dClasses
vector< string > dClasses
Definition: ModelCoreSpatial.h:139
ModelData::getForTypeParents
vector< string > getForTypeParents()
Definition: ModelData.cpp:122
Pixel::cumAlive_exp
vector< vector< double > > cumAlive_exp
This is the expected version of cumAlive, used for calculating profits.
Definition: Pixel.h:141
ModelData::getBoolSetting
bool getBoolSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1117
ModelCoreSpatial::runManagementModule
void runManagementModule()
computes regArea and expectedReturns
Definition: ModelCoreSpatial.cpp:875
MSG_NO_MSG
Do not actually output any message.
Definition: BaseClass.h:57
ModelCoreSpatial::initializePixelVolumes
void initializePixelVolumes()
distribuite regional exogenous volumes to pixel volumes using corine land cover area as weight ...
Definition: ModelCoreSpatial.cpp:1438
ModelData::createCombinationsVector
vector< vector< int > > createCombinationsVector(const int &nItems)
Return a vector containing any possible combination of nItems items (including any possible subset)...
Definition: ModelData.cpp:2055
ModelCoreSpatial::regPx
vector< Pixel * > regPx
Definition: ModelCoreSpatial.h:146
ModelCoreSpatial::sumRegionalForData
void sumRegionalForData()
computes vol, hV, harvestedArea, regArea and expReturns at reg level from the pixel level ...
Definition: ModelCoreSpatial.cpp:2194
Pixel::addMort
vector< vector< double > > addMort
Definition: Pixel.h:134
forType
Forest types (struct)
Definition: ModelData.h:293
ModelCoreSpatial::spd
void spd(const double &value_h, const string &type_h, const int &regId_h, const string &prodId_h, const int &year=DATA_NOW, const bool &allowCreate=false, const string &freeDim_h="") const
Definition: ModelCoreSpatial.h:125
BaseClass::i2s
string i2s(const int &int_h) const
integer to string conversion
Definition: BaseClass.cpp:219
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
BaseClass::getMax
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...
Definition: BaseClass.h:391
ModelData::getAvailableDeathTimber
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...
Definition: ModelData.cpp:2076
Gis::layerExist
bool layerExist(const string &layerName_h, bool exactMatch=true) const
Return a pointer to a layer given its name.
Definition: Gis.cpp:552
Pixel::changeValue
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)...
Definition: Pixel.cpp:135
ModelRegion::getValue
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 ...
Definition: ModelRegion.cpp:191
ModelData::calculateAnnualisedEquivalent
double calculateAnnualisedEquivalent(const double &amount_h, const int &years_h, const double &ir) const
Calculate the annual equivalent flow.
Definition: ModelData.cpp:1961
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
BaseClass::resetMapValues
void resetMapValues(map< K, V > &mymap, const V &value)
Reset all values stored in a map to the specified one.
Definition: BaseClass.h:341
ModelCoreSpatial::loadExogenousForestLayers
void loadExogenousForestLayers(const string &what)
Set pixel volumes (what="vol") OR areas (what="area") by specific forest types as defined in gis laye...
Definition: ModelCoreSpatial.cpp:2481
forType::memType
int memType
Definition: ModelData.h:296
Pixel::expType
double expType
Sampling derived expectation types of this agent (forest bilogical parameters: growth, mortality)
Definition: Pixel.h:145
Pixel::avalCoef
vector< double > avalCoef
Availability (of wood resources) coefficient. A [0,1] coefficient (new: by forest type) that reduces ...
Definition: Pixel.h:148
ModelCoreSpatial::priProducts
vector< string > priProducts
Definition: ModelCoreSpatial.h:136
ModelData::getAvailableAliveTimber
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...
Definition: ModelData.cpp:2124
Pixel::mort
vector< vector< double > > mort
Definition: Pixel.h:133
Output::printOptLog
void printOptLog(bool optimal, int &nIterations, double &obj)
Definition: Output.cpp:827
ThreadManager::SCD
Scheduler * SCD
the scheduler object (simulation-loops scheduler)
Definition: ThreadManager.h:75
Pixel::vMort
vector< vector< double > > vMort
Definition: Pixel.h:118
BaseClass::positionsToContent
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)
Definition: BaseClass.h:359
Pixel::vol_l
vector< vector< double > > vol_l
store the volumes of the previous year
Definition: Pixel.h:129
BaseClass::msgOut
void msgOut(const int &msgCode_h, const string &msg_h, const bool &refreshGUI_h=true) const
Overloaded function to print the output log.
Definition: BaseClass.cpp:50
Pixel::expectedAnnualIncome_carbon
vector< double > expectedAnnualIncome_carbon
Definition: Pixel.h:124
ThreadManager::GIS
Gis * GIS
GIS information and methods.
Definition: ThreadManager.h:73
ThreadManager
Thread manager. Responsable to manage the main thread and "speak" with the GUI.
Definition: ThreadManager.h:65
BaseClass::incrMapValue
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.
Definition: BaseClass.h:312
ModelData::getStringVectorSetting
vector< string > getStringVectorSetting(const string &name_h, int reg=WORLD) const
Definition: ModelData.cpp:1129
OP_SUM
Perform a SUM operation.
Definition: BaseClass.h:77
ModelCoreSpatial::resetPixelValues
void resetPixelValues()
swap volumes->lagged_volumes and reset the other pixel vectors
Definition: ModelCoreSpatial.cpp:1904
ModelData::deathTimberInventory_get
double deathTimberInventory_get(const iisskey &thekey)
Definition: ModelData.h:197
Pixel::vHa_exp
vector< vector< double > > vHa_exp
This is the expected version of vHa, used for calculating profits.
Definition: Pixel.h:140
ModelCoreSpatial::pDClasses
vector< string > pDClasses
Definition: ModelCoreSpatial.h:140
Pixel::initialDc0Area
vector< double > initialDc0Area
Definition: Pixel.h:109
Pixel::cumAlive
vector< vector< double > > cumAlive
Cumulative prob of remaining alive at beginnin of a given diam class.
Definition: Pixel.h:138
ModelData::getForTypeChilds_pos
vector< int > getForTypeChilds_pos(const string &forTypeId_h, bool all=false)
Definition: ModelData.cpp:109
BaseClass::getMaxPos
int getMaxPos(const vector< K > &v)
Returns the position of the maximum element in the vector (the last one in case of multiple equivalen...
Definition: BaseClass.h:383
ModelCoreSpatial::computeInventary
void computeInventary()
in=f(vol_t-1)
Definition: ModelCoreSpatial.cpp:2010
Pixel::vHa
vector< vector< double > > vHa
Volume at hectar by each diameter class [m^3/ha].
Definition: Pixel.h:137
ModelRegion::inResByAnyCombination
vector< double > inResByAnyCombination
Vector of inventory resource for each possible combination of primary products. This store both alive...
Definition: ModelRegion.h:86
ModelData::getStringSetting
string getStringSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1113
ModelCoreSpatial::runMarketModule
void runMarketModule()
computes st (supply total) and pw (weighted price). Optimisation inside.
Definition: ModelCoreSpatial.cpp:235
Carbon::registerTransports
void registerTransports(const double &distQ, const int &regId)
Registers the quantities emitted by transport of wood FROM a given region.
Definition: Carbon.cpp:277
ModelData::getNForTypesChilds
int getNForTypesChilds(const string &forTypeId_h)
Definition: ModelData.cpp:135
ModelCoreSpatial::cacheSettings
void cacheSettings()
just cache exogenous settings from ModelData
Definition: ModelCoreSpatial.cpp:1405
MSG_CRITICAL_ERROR
Print an error message and stop the model.
Definition: BaseClass.h:62
Carbon::registerDeathBiomass
void registerDeathBiomass(const double &value, const int &regId, const string &fType)
Registers the "death" of a given amount of biomass, storing it in the deathBiomass map...
Definition: Carbon.cpp:243
Pixel::vMortAdd
vector< vector< double > > vMortAdd
Definition: Pixel.h:119
Pixel::getSTData
double getSTData(const string &parName, const string &forName, int year, const string &d2="", double naToReturn=RETNA)
Definition: Pixel.cpp:246
ThreadManager::DO
Output * DO
data output
Definition: ThreadManager.h:76
Pixel::optDcChosen
vector< int > optDcChosen
Definition: Pixel.h:123
Pixel::getDoubleValue
double getDoubleValue(const string &layerName_h, const bool &returnZeroForNoValue=false) const
Return the value for a specific layer.
Definition: Pixel.cpp:158
BaseClass::d2s
string d2s(const double &double_h) const
double to string conversion
Definition: BaseClass.cpp:229
Pixel::expectedAnnualIncome_timber
vector< double > expectedAnnualIncome_timber
Definition: Pixel.h:125
ModelData::getForData
const double getForData(const string &type_h, const int &regId_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW)
Definition: ModelData.cpp:1281
ModelCoreSpatial::updateOtherMapData
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
Definition: ModelCoreSpatial.cpp:2157
BaseClass::vSum
int vSum(const vector< int > &vector_h) const
Definition: BaseClass.h:276
ModelData::getTempBool
bool getTempBool()
Definition: ModelData.h:145
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
ModelCoreSpatial::initialiseDeathTimber
void initialiseDeathTimber()
Set deathTimberInventory to zero for the previous years (under the hipotesis that we don&#39;t have advan...
Definition: ModelCoreSpatial.cpp:1605
Pixel::getPastRegArea
double getPastRegArea(const int &ft_idx, const int &year)
Definition: Pixel.cpp:390
ModelCoreSpatial::fTypes
vector< string > fTypes
Definition: ModelCoreSpatial.h:141
Pixel
Pixel-level class.
Definition: Pixel.h:47
Pixel::hProductivity
vector< vector< double > > hProductivity
Definition: Pixel.h:112
Carbon::initialiseProductsStocks
void initialiseProductsStocks(const vector< double > &qByProduct, const int &regId)
Initialises the stocks of products for the avgLive_* years before the simulation starts.
Definition: Carbon.cpp:203
ModelData::setTempBool
void setTempBool(bool tempBool_h)
Definition: ModelData.h:146
Pixel::vol
vector< vector< double > > vol
Definition: Pixel.h:90
ModelCoreSpatial::getVHaByYear
double getVHaByYear(const Pixel *px, const int &ft, const int &year, const double &extraBiomass_ratio, const int &regId) const
return the Volume/ha in a forest after a given number of year after planting, for a specific forest t...
Definition: ModelCoreSpatial.cpp:2831
Pixel::getMultiplier
double getMultiplier(const string &multiplierName, const string &forName, int year=DATA_NOW)
Definition: Pixel.cpp:184
ModelCoreSpatial::allProducts
vector< string > allProducts
Definition: ModelCoreSpatial.h:138
ModelCoreSpatial::secProducts
vector< string > secProducts
Definition: ModelCoreSpatial.h:137
Pixel::cumTp_exp
vector< vector< double > > cumTp_exp
This is the expected version of cumTp, used for calculating profits.
Definition: Pixel.h:139
Pixel::expTypePrices
double expTypePrices
Sampling derived expectation types of this agent (prices)
Definition: Pixel.h:146
Carbon::registerHarvesting
void registerHarvesting(const double &value, const int &regId, const string &fType)
Registers the harvesting of trees increasing the value of cumEmittedForOper.
Definition: Carbon.cpp:225
MSG_DEBUG
Print a debug message, normally filtered out.
Definition: BaseClass.h:58
ThreadManager::OPT
Ipopt::SmartPtr< Ipopt::TNLP > OPT
Market optimisation.
Definition: ThreadManager.h:81
ThreadManager::CBAL
Carbon * CBAL
Module for the Carbon Balance.
Definition: ThreadManager.h:79
Pixel::hVol_byPrd
vector< vector< vector< double > > > hVol_byPrd
Definition: Pixel.h:113
MSG_WARNING
Print a WARNING message.
Definition: BaseClass.h:60
Pixel::area_l
vector< vector< double > > area_l
store the areas of the previous year
Definition: Pixel.h:130
Gis::getAllPlotsByRegion
vector< Pixel * > getAllPlotsByRegion(ModelRegion &region_h, bool shuffle=false)
Return the vector of all plots by a specific region (main region or subregion), optionally shuffled;...
Definition: Gis.cpp:855
ModelData::getForTypeParentId
string getForTypeParentId(const string &forTypeId_h)
Definition: ModelData.cpp:90
ModelData::getForTypeIds
vector< string > getForTypeIds(bool all=false)
By default it doesn&#39;t return forTypes used only as input.
Definition: ModelData.cpp:430
BaseClass::getSd
double getSd(const vector< K > &v, bool sample=true)
Definition: BaseClass.h:408
ModelCoreSpatial::gpd
double gpd(const string &type_h, const int &regId_h, const string &prodId_h, const int &year=DATA_NOW, const string &freeDim_h="") const
Definition: ModelCoreSpatial.h:123
BaseClass::vectorToMap
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...
Definition: BaseClass.h:349
Pixel::optFtChosen
vector< int > optFtChosen
Definition: Pixel.h:122
Pixel::optDc
vector< int > optDc
Definition: Pixel.h:121
Pixel::expectedReturns
vector< double > expectedReturns
Definition: Pixel.h:120
BaseClass::findMap
V findMap(const map< K, V > &mymap, const K &key, const int &error_level=MSG_CRITICAL_ERROR, const V &notFoundValue=numeric_limits< V >::min()) const
Lookup a map for a value. Return the value starting from the key.
Definition: BaseClass.h:286
Pixel::tp
vector< vector< double > > tp
Definition: Pixel.h:135
Pixel::getPathMortality
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...
Definition: Pixel.cpp:329
ModelCoreSpatial::ModelCoreSpatial
ModelCoreSpatial(ThreadManager *MTHREAD_h)
Definition: ModelCoreSpatial.cpp:39
ModelCoreSpatial::getAvgAgeByDc
double getAvgAgeByDc(Pixel *px, int ft, int dc)
return the average age of a tree at a given diameter class, using the cumTp vector ...
Definition: ModelCoreSpatial.cpp:2820
ModelData::getErrorLevel
int getErrorLevel()
Definition: ModelData.h:144
Pixel::swap
void swap(const int &swap_what)
Assign to the delayed value the current values, e.g. vol_l = vol.
Definition: Pixel.cpp:422
forType::forLayer
string forLayer
Definition: ModelData.h:297
ModelRegion::getSiblings
vector< ModelRegion * > getSiblings(bool excludeResidual=true)
Return a vector of pointers to the siblings regions.
Definition: ModelRegion.cpp:85
Output::printDetailedHV
void printDetailedHV(map< tr1::array< string, 4 >, double > hVol_byPrd)
Definition: Output.cpp:1173
Pixel::expectedReturnsNotCorrByRa
vector< double > expectedReturnsNotCorrByRa
by ft. Attenction, reported expReturns at "forest" level (compared with those at forest type level) d...
Definition: Pixel.h:127
Pixel::hVol
vector< vector< double > > hVol
Definition: Pixel.h:111
BaseClass::incrOrAddMapValue
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.
Definition: BaseClass.h:325
Carbon::initialiseEmissionCounters
void initialiseEmissionCounters()
Initialises the emission counters to zero.
Definition: Carbon.cpp:158
ModelData::setErrorLevel
void setErrorLevel(int errorLevel_h)
Definition: ModelData.h:143
ModelCoreSpatial::initMarketModule
void initMarketModule()
computes st and pw for second year and several needed-only-at-t0-vars for the market module ...
Definition: ModelCoreSpatial.cpp:99
ModelData::getRegionIds
vector< int > getRegionIds(int level_h, bool excludeResidual=true)
Definition: ModelData.cpp:366
ModelCoreSpatial::computeExpectedPrice
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.
Definition: ModelCoreSpatial.cpp:2464
Output::print
void print(bool earlyPrint=true)
Print output. If earlyPrinting it doesn&#39;t print some stuff for which we don&#39;t yet have values...
Definition: Output.cpp:426
Carbon::HWP_eol2energy
void HWP_eol2energy()
Computes the energy substitution for the quota of HWP that reachs end of life and doesn&#39;t go to landf...
Definition: Carbon.cpp:288
ModelCoreSpatial::initializePixelArea
void initializePixelArea()
compute px->area for each ft and dc
Definition: ModelCoreSpatial.cpp:1630
ModelCoreSpatial::regIds2
vector< int > regIds2
Definition: ModelCoreSpatial.h:135
ModelData::getDoubleSetting
double getDoubleSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1109
Pixel::cumTp
vector< vector< double > > cumTp
This is time of passage to REACH a diameter class (while the exogenous tp by diameter class is the ti...
Definition: Pixel.h:136
Pixel::hArea
vector< vector< double > > hArea
Definition: Pixel.h:110
ModelCoreSpatial::sfd
void sfd(const double &value_h, const string &type_h, const int &regId_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW, const bool &allowCreate=false) const
Definition: ModelCoreSpatial.h:126
ModelRegion::getRegSName
string getRegSName() const
Definition: ModelRegion.h:67
Pixel::vReg
vector< double > vReg
Definition: Pixel.h:117
ModelCoreSpatial::gfd
double gfd(const string &type_h, const int &regId_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW) const
Definition: ModelCoreSpatial.h:124
ModelCoreSpatial::app
bool app(const string &prod_h, const string &forType_h, const string &dClass_h) const
Definition: ModelCoreSpatial.h:127
Pixel::beta
vector< vector< double > > beta
Definition: Pixel.h:132
ModelRegion::getMyPixels
vector< Pixel * > getMyPixels()
Definition: ModelRegion.h:86
MSG_INFO
Print an INFO message.
Definition: BaseClass.h:59
ModelCoreSpatial::assignSpMultiplierPropToVols
void assignSpMultiplierPropToVols()
ModelCoreSpatial::assignSpMultiplierPropToVols assigns the spatial multiplier (used in the time of re...
Definition: ModelCoreSpatial.cpp:1502
ModelData::getRegion
ModelRegion * getRegion(int regId_h)
Definition: ModelData.cpp:346
ModelCoreSpatial::cacheDynamicSettings
void cacheDynamicSettings()
cache settings that may change with time
Definition: ModelCoreSpatial.cpp:1425
ModelCoreSpatial::initialiseCarbonModule
void initialiseCarbonModule()
call initialiseDeathBiomassStocks(), initialiseProductsStocks() and initialiseEmissionCounters() ...
Definition: ModelCoreSpatial.cpp:1576
ModelCoreSpatial::cachePixelExogenousData
void cachePixelExogenousData()
computes pixel level tp, meta and mort
Definition: ModelCoreSpatial.cpp:1945
BaseClass::getAvg
double getAvg(const vector< K > &v)
Returns the average of the elements in the vector.
Definition: BaseClass.h:399
DIAM_ALL
#define DIAM_ALL
All diameter classes.
Definition: BaseClass.h:157
ModelCoreSpatial::updateMapAreas
void updateMapAreas()
computes forArea_{ft}
Definition: ModelCoreSpatial.cpp:2091
ModelData::deathTimberInventory_incrOrAdd
void deathTimberInventory_incrOrAdd(const iisskey &thekey, double value_h)
Definition: ModelData.h:195
ModelCoreSpatial::registerCarbonEvents
void registerCarbonEvents()
call registerHarvesting(), registerDeathBiomass(), registerProducts() and registerTransports() ...
Definition: ModelCoreSpatial.cpp:2406
ModelData::getForTypeChilds
vector< string > getForTypeChilds(const string &forTypeId_h)
Definition: ModelData.cpp:98
ModelCoreSpatial::computeCumulativeData
void computeCumulativeData()
computes cumTp_exp, vHa_exp, vHa
Definition: ModelCoreSpatial.cpp:1700
Carbon::registerProducts
void registerProducts(const double &value, const int &regId, const string &productName)
Registers the production of a given amount of products, storing it in the products maps...
Definition: Carbon.cpp:257
Pixel::getID
double getID() const
Definition: Pixel.h:67
ModelCoreSpatial::printDebugInitRegionalValues
void printDebugInitRegionalValues()
print initial inv, st, sl and sa in each region
Definition: ModelCoreSpatial.cpp:2662
ModelCoreSpatial::computeEconomicBalances
void computeEconomicBalances()
compute the policy balances (- -> costs; + -> rev) and the producer/consumer surpluses ...
Definition: ModelCoreSpatial.cpp:2735
ModelData::getForType
forType * getForType(int position)
Definition: ModelData.h:128
Pixel::regArea
map< int, vector< double > > regArea
Definition: Pixel.h:114
ModelData::getMaxYearUsableDeathTimber
const int getMaxYearUsableDeathTimber(const string &prod_h, const string &forType_h, const string &dClass_h)
Definition: ModelData.cpp:468
ModelCoreSpatial::allocateHarvesting
vector< double > allocateHarvesting(vector< double > total_st, const int &regId)
Using the deathTimberInventory map, this function allocate the total st in st from death timber (that...
Definition: ModelCoreSpatial.cpp:2696