FFSM++  1.1.0
French Forest Sector Model ++
ModelCore.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 
25 #include "IpIpoptApplication.hpp"
26 #include "IpSolveStatistics.hpp"
27 
28 #include "ModelCore.h"
29 #include "ModelCoreSpatial.h"
30 #include "ModelData.h"
31 #include "ThreadManager.h"
32 #include "Opt.h"
33 #include "Scheduler.h"
34 #include "Gis.h"
35 
36 
37 ModelCore::ModelCore(ThreadManager* MTHREAD_h){
38  MTHREAD = MTHREAD_h;
39 }
40 
41 ModelCore::~ModelCore(){
42 
43 }
44 
45 
46 /**
47 IMPORTANT NOTE: Volumes in Mm^3, Areas in the model in Ha (10000 m^2), in the layers in m^2
48 */
49 void
50 ModelCore::runInitPeriod(){
51  /**
52  Some importan notes:
53  V (volumes) -> at the end of the year
54  In (inventary) -> at the beginning of the year
55  Volumes are in Mm^3, Areas in the model in Ha (10000 m^2), in the layers in m^2
56  */
57  cacheSettings(); // cashe things like first year, second year, dClasses...
58  initMarketModule(); // inside it uses first year, second year
59  MTHREAD->DO->print();
60  MTHREAD->SCD->advanceYear(); // 2005->2006
61  computeInventary(); // in=f(vol_t-1)
62  computeCumulativeData(); // compute cumTp_exp, vHa_exp, vHa
63  runBiologicalModule(); //
64  runManagementModule();
65  updateMapAreas(); // update the forArea_{ft} layer on each pixel as old value-hArea+regArea
66  MTHREAD->DO->print();
67 }
68 
69 void
70 ModelCore::runSimulationYear(){
71  int thisYear = MTHREAD->SCD->getYear();
72  computeInventary(); // in=f(vol_t-1)
73  runMarketModule();
74  computeCumulativeData(); // compute cumTp_exp, vHa_exp
75  runBiologicalModule();
76 
77  /*double sl = gpd("sl",11041,'softWRoundW');
78  double pl = gpd("pl",11041,'softWRoundW');
79  double sa = gpd("sa",11041,'softWRoundW');
80  double pworld = gpd("pl", WL2,'softWRoundW');
81  double st = gpd("st",11041,'softWRoundW');
82  double pw = (sl*pl+sa*pworld)/st;
83  cout << thisYear <<
84  */
85 
86  runManagementModule();
87  updateMapAreas();
88  MTHREAD->DO->print();
89 }
90 
91 
92 void
93 ModelCore::initMarketModule(){
94  msgOut(MSG_INFO, "Starting market module (init stage)..");
95  for(uint i=0;i<regIds2.size();i++){
96  int r2 = regIds2[i];
97 
98  //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);
99  for(uint sp=0;sp<secProducts.size();sp++){
100  double value = 0;
101  for (uint pp=0;pp<priProducts.size();pp++){
102  value += gpd("pl",r2,priProducts[pp],secondYear)*
103  gpd("a",r2,priProducts[pp],secondYear,secProducts[sp]);
104  }
105  value += gpd("m",r2,secProducts[sp],secondYear);
106  spd(value,"pl",r2,secProducts[sp],secondYear);
107  }
108  // RPAR('dl',i,p_pr,t-1) = sum(p_tr, a(p_pr,p_tr)*RPAR('sl',i,p_tr,t-1));
109  for (uint pp=0;pp<priProducts.size();pp++){
110  double value=0;
111  for(uint sp=0;sp<secProducts.size();sp++){
112  value += gpd("sl",r2,secProducts[sp],secondYear)*
113  gpd("a",r2,priProducts[pp],secondYear, secProducts[sp]);
114  }
115  spd(value,"dl",r2,priProducts[pp],secondYear,true);
116  }
117  // RPAR('st',i,prd,t-1) = RPAR('sl',i,prd,t-1)+RPAR('sa',i,prd,t-1);
118  // RPAR('dt',i,prd,t-1) = RPAR('dl',i,prd,t-1)+RPAR('da',i,prd,t-1);
119  for (uint ap=0;ap<allProducts.size();ap++){
120  double stvalue = gpd("sl",r2,allProducts[ap],secondYear)
121  + gpd("sa",r2,allProducts[ap],secondYear);
122  double dtvalue = gpd("dl",r2,allProducts[ap],secondYear)
123  + gpd("da",r2,allProducts[ap],secondYear);
124  spd(stvalue,"st",r2,allProducts[ap],secondYear,true);
125  spd(dtvalue,"dt",r2,allProducts[ap],secondYear,true);
126  }
127 
128  // 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)));
129  // p1(i,p_tr) = 1-q1(i,p_tr);
130  // 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));
131  // 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);
132  // 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);
133  // 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
134  // K(i,p_tr,t-1) = k1(i,p_tr)*RPAR('sl',i,p_tr,t-1);
135  for(uint sp=0;sp<secProducts.size();sp++){
136  double psi = gpd("psi",r2,secProducts[sp],secondYear);
137  double dl = gpd("dl",r2,secProducts[sp],secondYear);
138  double da = gpd("da",r2,secProducts[sp],secondYear);
139  double pl = gpd("pl",r2,secProducts[sp],secondYear);
140  double sl = gpd("sl",r2,secProducts[sp],secondYear);
141  double k1 = gpd("k1",r2,secProducts[sp],secondYear);
142  double pWo = gpd("pl",WL2,secProducts[sp],secondYear); // World price (local price for region 99999)
143 
144 
145  double q1 = 1/ ( 1+pow(dl/da,1/psi)*(pl/pWo) );
146  double p1 = 1-q1;
147  double dc = pow(
148  q1*pow(da,(psi-1)/psi) + p1*pow(dl,(psi-1)/psi)
149  ,
150  psi/(psi-1)
151  );
152  double pc = (da/dc)*pWo
153  +(dl/dc)*pl;
154  double pw = (dl*pl+da*pWo)/(dl+da);
155  double k = k1*sl;
156 
157  spd(q1,"q1",r2,secProducts[sp],firstYear,true);
158  spd(p1,"p1",r2,secProducts[sp],firstYear,true);
159  spd(dc,"dc",r2,secProducts[sp],secondYear,true);
160  spd(pc,"pc",r2,secProducts[sp],secondYear,true);
161  spd(pw,"pw",r2,secProducts[sp],secondYear,true);
162  spd(k,"k",r2,secProducts[sp],secondYear,true);
163  }
164 
165  // 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)));
166  // r1(i,p_pr) = 1-t1(i,p_pr);
167  // 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))
168  // 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);
169  // 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
170  for(uint pp=0;pp<priProducts.size();pp++){
171 
172  double sl = gpd("sl",r2,priProducts[pp],secondYear);
173  double sa = gpd("sa",r2,priProducts[pp],secondYear);
174  double eta = gpd("eta",r2,priProducts[pp],secondYear);
175  double pl = gpd("pl",r2,priProducts[pp],secondYear);
176  double pWo = gpd("pl",WL2,priProducts[pp],secondYear); // World price (local price for region 99999)
177 
178 
179  double t1 = 1/ ( 1+(pow(sl/sa,1/eta))*(pl/pWo) );
180  double r1 = 1-t1;
181  double sc = pow(
182  t1*pow(sa,(eta-1)/eta) + r1*pow(sl,(eta-1)/eta)
183  ,
184  eta/(eta-1)
185  );
186  double pc = (sa/sc)*pWo+(sl/sc)*pl;
187  double pw = (sl*pl+sa*pWo)/(sl+sa);
188 
189  spd(t1,"t1",r2,priProducts[pp],firstYear,true);
190  spd(r1,"r1",r2,priProducts[pp],firstYear,true);
191  spd(sc,"sc",r2,priProducts[pp],secondYear,true);
192  spd(pc,"pc",r2,priProducts[pp],secondYear,true);
193  spd(pw,"pw",r2,priProducts[pp],secondYear,true);
194  }
195  } // end of each region
196 
197 
198  // initializing the exports to zero quantities
199  for(uint r1=0;r1<l2r.size();r1++){
200  for(uint r2=0;r2<l2r[r1].size();r2++){
201  for(uint p=0;p<allProducts.size();p++){
202  for(uint r2To=0;r2To<l2r[r1].size();r2To++){
203  spd(0,"rt",l2r[r1][r2],allProducts[p],secondYear,true,i2s(l2r[r1][r2To])); // regional trade, it was exp in gams
204  }
205  }
206  }
207  } // end of r1 region
208 }
209 
210 void
211 ModelCore::runMarketModule(){
212 
213  // *** PRE-OPTIMISATION YEARLY OPERATIONS..
214 
215  // Updating variables
216  // notes:
217  // - dispo_sup: not actually entering anywhere, forgiving it for now..
218  // - dc0: not needed, it is just the first year demand composite
219  int thisYear = MTHREAD->SCD->getYear();
220  int previousYear = thisYear-1;
221 
222  for(uint i=0;i<regIds2.size();i++){
223  int r2 = regIds2[i];
224 
225  // K(i,p_tr,t) = (1+g1(i,p_tr))*K(i,p_tr,t-1);
226  // AA(i,p_tr) = (sigma(p_tr)/(sigma(p_tr)+1))*RPAR('pc',i,p_tr,t-1)*(RPAR('dc',i,p_tr,t-1)**(-1/sigma(p_tr)));
227  // GG(i,p_tr) = RPAR('dc',i,p_tr,t-1)*((RPAR('pc',i,p_tr,t-1))**(-sigma(p_tr))); //alpha
228  for(uint sp=0;sp<secProducts.size();sp++){
229  double g1 = gpd("g1",r2,secProducts[sp],previousYear);
230  double sigma = gpd("sigma",r2,secProducts[sp]);
231  double pc_1 = gpd("pc",r2,secProducts[sp],previousYear);
232  double dc_1 = gpd("dc",r2,secProducts[sp],previousYear);
233  double k_1 = gpd("k",r2,secProducts[sp],previousYear);
234 
235  double k = (1+g1)*k_1;
236  double aa = (sigma/(sigma+1))*pc_1*pow(dc_1,-1/sigma);
237  double gg = dc_1*pow(pc_1,-sigma); //alpha
238 
239  spd(k, "k" ,r2,secProducts[sp]);
240  spd(aa,"aa",r2,secProducts[sp],DATA_NOW,true);
241  spd(gg,"gg",r2,secProducts[sp],DATA_NOW,true);
242  }
243 
244  // 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));
245  // 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
246  for(uint pp=0;pp<priProducts.size();pp++){
247  double gamma = gpd("gamma",r2,priProducts[pp]);
248  double sigma = gpd("sigma",r2,priProducts[pp]);
249  double pc_1 = gpd("pc",r2,priProducts[pp],previousYear);
250  double sc_1 = gpd("sc",r2,priProducts[pp],previousYear);
251  double in = gpd("in",r2,priProducts[pp]);
252  double in_1 = gpd("in",r2,priProducts[pp],previousYear);
253 
254  double bb = (sigma/(sigma+1))*pc_1*pow(sc_1,-1/sigma)*pow(in_1/in,gamma/sigma);
255  double ff = sc_1*pow(pc_1,-sigma)*pow(in/in_1,gamma); //chi
256 
257  spd(bb,"bb",r2,priProducts[pp],DATA_NOW,true);
258  spd(ff,"ff",r2,priProducts[pp],DATA_NOW,true);
259  }
260 
261  } // end for each region in level 2 (and updating variables)
262 
263  // *** OPTIMISATION....
264 
265  // Create an instance of the IpoptApplication
266  //Opt *OPTa = new Opt(MTHREAD);
267  //SmartPtr<TNLP> OPTa = new Opt(MTHREAD);
268  //int initialOptYear = MTHREAD->MD->getIntSetting("initialOptYear");
269  SmartPtr<IpoptApplication> application = new IpoptApplication();
270  //if(thisYear == initialOptYear){
271  //application = new IpoptApplication();
272  //} else {
273  // application->Options()->SetStringValue("warm_start_init_point", "yes");
274  //}
275  string linearSolver = MTHREAD->MD->getStringSetting("linearSolver",DATA_NOW);
276  application->Options()->SetStringValue("linear_solver", linearSolver); // default in ipopt is ma27
277  //application->Options()->SetStringValue("hessian_approximation", "limited-memory"); // quasi-newton approximation of the hessian
278  //application->Options()->SetIntegerValue("mumps_mem_percent", 100);
279  application->Options()->SetNumericValue("obj_scaling_factor", -1); // maximisation
280  application->Options()->SetNumericValue("max_cpu_time", 1800); // max 1/2 hour to find the optimus for one single year
281  application->Options()->SetStringValue("check_derivatives_for_naninf", "yes"); // detect error but may crash the application.. TO.DO catch this error!
282  //application->Options()->SetStringValue("nlp_scaling_method", "equilibration-based"); // much worster
283  // Initialize the IpoptApplication and process the options
284  ApplicationReturnStatus status;
285  status = application->Initialize();
286  if (status != Solve_Succeeded) {
287  printf("\n\n*** Error during initialization!\n");
288  msgOut(MSG_INFO,"Error during initialization! Do you have the solver compiled for the specified linear solver?");
289  return;
290  }
291 
292 
293  msgOut(MSG_INFO,"Running optimisation problem for this year (it may take a few minutes for large models)..");
294  status = application->OptimizeTNLP(MTHREAD->OPT);
295 
296  // *** POST OPTIMISATION....
297 
298  // post-equilibrium variables->parameters assignments..
299  // RPAR(type,i,prd,t) = RVAR.l(type,i,prd);
300  // EX(i,j,prd,t) = EXP.l(i,j,prd);
301  // ObjT(t) = Obj.l ;
302  // ==> in Opt::finalize_solution()
303 
304  // Retrieve some statistics about the solve
305  if (status == Solve_Succeeded) {
306  Index iter_count = application->Statistics()->IterationCount();
307  Number final_obj = application->Statistics()->FinalObjective();
308  printf("\n*** The problem solved in %d iterations!\n", iter_count);
309  printf("\n*** The final value of the objective function is %e.\n", final_obj);
310  msgOut(MSG_INFO, "The problem solved successfully in "+i2s(iter_count)+" iterations.");
311  int icount = iter_count;
312  double obj = final_obj;
313  MTHREAD->DO->printOptLog(true, icount, obj);
314  } else {
315  //Number final_obj = application->Statistics()->FinalObjective();
316  cout << "***ERROR: MODEL DIDN'T SOLVE FOR THIS YEAR" << endl;
317  msgOut(MSG_CRITICAL_ERROR, "Model DIDN'T SOLVE for this year");
318  // 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
319  //Index iter_count = application->Statistics()->IterationCount(); // sys error if model didn't solve
320  //Number final_obj = application->Statistics()->FinalObjective();
321  int icount = 0;
322  double obj = 0;
323  MTHREAD->DO->printOptLog(false, icount, obj);
324  }
325 
326  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")
327  int regId = regIds2[r2];
328 
329  // // total supply and total demand..
330  // RPAR('st',i,prd,t) = RPAR('sl',i,prd,t)+RPAR('sa',i,prd,t);
331  // RPAR('dt',i,prd,t) = RPAR('dl',i,prd,t)+RPAR('da',i,prd,t);
332  // // weighted prices.. //changed 20120419
333  // 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
334  // 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
335  for (uint p=0;p<allProducts.size();p++){
336  double st = gpd("sl",regId,allProducts[p])+gpd("sa",regId,allProducts[p]);
337  double dt = gpd("dl",regId,allProducts[p])+gpd("da",regId,allProducts[p]);
338  spd(st,"st",regId,allProducts[p]);
339  spd(dt,"dt",regId,allProducts[p]);
340  }
341  for (uint p=0;p<secProducts.size();p++){
342  double dl = gpd("dl",regId,secProducts[p]);
343  double pl = gpd("pl",regId,secProducts[p]);
344  double da = gpd("da",regId,secProducts[p]); // bug corrected 20120913
345  double pworld = gpd("pl", WL2,secProducts[p]);
346  double dt = gpd("dt",regId,secProducts[p]);
347  double pw = (dl*pl+da*pworld)/dt;
348  spd(pw,"pw",regId,secProducts[p]);
349  }
350  for (uint p=0;p<priProducts.size();p++){
351  double sl = gpd("sl",regId,priProducts[p]);
352  double pl = gpd("pl",regId,priProducts[p]);
353  double sa = gpd("sa",regId,priProducts[p]); // bug corrected 20120913
354  double pworld = gpd("pl", WL2,priProducts[p]);
355  double st = gpd("st",regId,priProducts[p]);
356  double pw = (sl*pl+sa*pworld)/st;
357  spd(pw,"pw",regId,priProducts[p]);
358  }
359  } // end of foreach region
360 }
361 
362 void
363 ModelCore::runBiologicalModule(){
364 
365  msgOut(MSG_INFO, "Starting resource module..");
366  hV_byPrd.clear();
367  int thisYear = MTHREAD->SCD->getYear();
368  int previousYear = thisYear-1;
369 
370  for(uint i=0;i<regIds2.size();i++){
371 
372  int r2 = regIds2[i];
373  int regId = r2;
374  // Post optimisation biological mobule..
375  vector < vector < vector <double> > > hV_byPrd_regional;
376  for(uint j=0;j<fTypes.size();j++){
377  string ft = fTypes[j];
378  vector < vector <double> > hV_byPrd_ft;
379 
380  // calculating the regeneration..
381  // if we are in a year where the time of passage has not yet been reached
382  // for the specific i,e,l then we use the exogenous Vregen, otherwise we
383  // calculate it
384  //if ( not scen("fxVreg") ,
385  // loop( (i,essence,lambda),
386  // if( ord(t)>=(tp_u1(i,essence,lambda)+2),
387  // Vregen(i,lambda,essence,t)=regArea(i,essence,lambda,t-tp_u1(i,essence,lambda))*volHa_u1(i,essence,lambda)/1000000 ;
388  // );
389  // );
390  //);
391  int tp_u0 = gfd("tp",regId,ft,dClasses[0]); // time of passage to reach the first diameter class // bug 20140318, added ceil
392  if(regType != "fixed" && (thisYear-secondYear) >= tp_u0 ) { // T.O.D.O to be checked -> 20121109 OK
393  double pastRegArea = gfd("regArea",regId,ft,"",thisYear-tp_u0);
394  double vHa = gfd("vHa",regId,ft,dClasses[1]);
395  //cout << "vHa - entryVolHa: " << vHa << " / " << entryVolHa << endl;
396  double vReg = pastRegArea*vHa/1000000; // T.O.D.O: check the 1000000. -> Should be ok, as area in ha vol in Mm^3
397  sfd(vReg,"vReg",regId,ft,"");
398  }
399 
400  for (uint u=0; u<dClasses.size(); u++){
401  string dc = dClasses[u];
402  double hr = 0;
403  double pastYearVol = u?gfd("vol",regId,ft,dc,previousYear):0.;
404  double hV_mort = 0.; /// \todo Harvest volumes from death trees
405  vector <double> hV_byPrd_dc;
406 
407  // harvesting rate & volumes...
408  //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));
409  //hV(u,i,essence,lambda,t) = hr(u,i,essence,lambda,t) * V(u,i,lambda,essence,t-1);
410  //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);
411  //double debug =0;
412  for(uint pp=0;pp<priProducts.size();pp++){
413  double st = gpd("st",regId,priProducts[pp]);
414  double in = gpd("in",regId,priProducts[pp]);
415  double hr_pr = u?app(priProducts[pp],ft,dc)*st/ in:0;
416  double hV_byPrd_dc_prd = hr_pr*pastYearVol;
417  hr += hr_pr;
418  hV_byPrd_dc.push_back( hV_byPrd_dc_prd);
419  //debug += hV_byPrd_dc_prd;
420  }
421  double hV = hr*pastYearVol;
422  //double debug2 = debug-hV;
423 
424  // test passed 20131203
425  //if(debug2 < -0.000000000001 || debug2 > 0.000000000001){
426  // cout << "Problems!" << endl;
427  //}
428 
429  // post harvesting remained volumes computation..
430  // loop(u$(ord(u)=1),
431  // first diameter class, no harvesting and fixed regenaration..
432  // 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)
433  // +Vregen(i,lambda,essence,t);
434  // );
435  // loop(u$(ord(u)>1),
436  // generic case..
437  // V(u,i,lambda,essence,t)=((1-1/(tp(u,i,lambda,essence))
438  // -mort(u,i,lambda,essence) - hr(u,i,essence,lambda,t))*V(u,i,lambda,essence,t-1)
439  // +(1/(tp(u-1,i,lambda,essence)))*beta(u,i,lambda,essence)*V(u-1,i,lambda,essence,t-1));
440  double vol;
441  double newMortVol; // new mortality volumes
442  double tp = gfd("tp",regId,ft,dc);
443  double mort = u?gfd("mortCoef",regId,ft,dc):0.;
444  double vReg = 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!
445  double beta = u?gfd("betaCoef",regId,ft,dc):0.;
446  //double hv2fa = gfd("hv2fa",regId,ft,dc);
447  double vHa = gfd("vHa",regId,ft,dc);
448  double finalHarvestFlag = gfd("finalHarvestFlag",regId,ft,dc);
449 
450  if(u==0){
451  vol = 0.;
452  } else if(u==1){
453  vol = (1-1/tp-mort)*pastYearVol+vReg;
454  newMortVol = mort*pastYearVol;
455  double debug = vol;
456  } else {
457  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
458  double tp_1 = gfd("tp",regId,ft,dClasses[u-1]);
459  double pastYearVol_1 = gfd("vol",regId,ft,dClasses[u-1],previousYear);
460  vol = (1-inc-mort-hr)*pastYearVol+(1/tp_1)*beta*pastYearVol_1;
461  newMortVol = mort*pastYearVol;
462  double debug = vol;
463  }
464  double freeArea_byU = u?finalHarvestFlag*1000000*hV/vHa:0; // volumes are in Mm^3, area in ha, vHa in m^3/ha
465  //double debug = hv2fa*hr*pastYearVol*100;
466  //cout << "regId|ft|dc| debug | freeArea: " << r2 << "|"<<ft<<"|"<<dc<<"| "<< debug << " | " << freeArea_byU << endl;
467 
468  sfd(hr,"hr",regId,ft,dc,DATA_NOW,true);
469  sfd(hV,"hV",regId,ft,dc,DATA_NOW,true);
470  sfd(vol,"vol",regId,ft,dc,DATA_NOW,true); // allowCreate needed for u==0
471  sfd(newMortVol,"mortV",regId,ft,dc,DATA_NOW,true);
472 
473  sfd(freeArea_byU,"harvestedArea",regId,ft,dc,DATA_NOW, true);
474  hV_byPrd_ft.push_back(hV_byPrd_dc);
475  } // end foreach diameter classes
476  hV_byPrd_regional.push_back(hV_byPrd_ft);
477  } // end of each forest type
478  hV_byPrd.push_back(hV_byPrd_regional);
479  } // end of for each region
480 
481 }
482 
483 void
484 ModelCore::runManagementModule(){
485 
486  msgOut(MSG_INFO, "Starting management module..");
487  //int thisYear = MTHREAD->SCD->getYear();
488  //int previousYear = thisYear-1;
489  MTHREAD->DO->expReturnsDebug.clear();
490  int outputLevel = MTHREAD->MD->getIntSetting("outputLevel");
491  bool weightedAverageExpectedReturns = MTHREAD->MD->getBoolSetting("weightedAverageExpectedReturns");
492 
493  //vector <vector < vector <vector <vector <double> > > > > expReturnsDebug; ///< l2_region, for type, d.c., pr prod, variable name
494  //cout << "year/dc/pp/eai/cumTp/vHa/pw" << endl;
495 
496  int thisYear = MTHREAD->SCD->getYear();
497  double ir = MTHREAD->MD->getDoubleSetting("ir");
498 
499  for(uint i=0;i<regIds2.size();i++){
500  vector < vector <vector <vector <double> > > > expReturnsDebug_region;
501 
502  int r2 = regIds2[i];
503  int regId = r2;
504  vector <double> cachedExpectedReturnsByFt;
505 
506  // PART 1: COMPUTING THE EXPECTED RETURNS FOR EACH FOREST TYPE
507 
508  for(uint j=0;j<fTypes.size();j++){
509  string ft = fTypes[j];
510  vector <vector <vector <double> > > expReturnsDebug_ft;
511  // Post optimisation management mobule..
512 
513  //if(regType != "fixed" && regType != "fromHrLevel"){
514  // 20120910, Antonello: changed.. calculating the expected returns also for fixed and fromHrLevel regeneration (then not used but gives indication)
515  // calculating the expected returns..
516  // loop ( (u,i,essence,lambda,p_pr),
517  // if (sum(u2, hV(u2,i,essence,lambda,t))= 0,
518  // expRetPondCoef(u,i,essence,lambda,p_pr) = 0;
519  // else
520  // expRetPondCoef(u,i,essence,lambda,p_pr) = hV_byPrd(u,i,essence,lambda,p_pr,t)/ sum(u2, hV(u2,i,essence,lambda,t));
521  // );
522  // );
523  // expReturns(i,essence,lambda) = sum( (u,p_pr),
524  // RPAR("pl",i,p_pr,t)*hv2fa(i,essence,lambda,u)*(1/df_byFT(u,i,lambda,essence))* // df_byFT(u,i,lambda,essence)
525  // expRetPondCoef(u,i,essence,lambda,p_pr)
526  // );
527  double hV_byFT = 0.; // gfd("hV",regId,ft,DIAM_PROD); // it must include only final harvested products in order to act as weightering agent
528  double expReturns = 0;
529 
530 
531  for (uint u=0; u<dClasses.size(); u++){
532  string dc = dClasses[u];
533  double finalHarvestFlag = gfd("finalHarvestFlag",regId,ft,dc);
534  double hV = gfd("hV",regId,ft,dc);
535  hV_byFT += finalHarvestFlag * hV;
536  }
537 
538  if(hV_byFT==0. || !weightedAverageExpectedReturns){ // nothing has been harvested in this pixel for this specific forest type. Let's calculate the combination product/diameter class with the highest expected return
539  for (uint u=0; u<dClasses.size(); u++){
540  vector <vector <double> > expReturnsDebug_dc;
541  string dc = dClasses[u];
542  double vHa = gfd("vHa_exp",regId,ft,dc);
543  double finalHarvestFlag = gfd("finalHarvestFlag",regId,ft,dc);
544  double cumTp_u = gfd("cumTp_exp",regId,ft,dc);
545  for (uint pp=0;pp<priProducts.size();pp++){
546  vector <double> expReturnsDebug_pp;
547  double pw = gpd("pw",regId,priProducts[pp]);
548  double raw_amount = finalHarvestFlag*pw*vHa*app(priProducts[pp],ft,dc); // B.U.G. 20121126, it was missing app(pp,ft,dc) !!
549  double anualised_amount = MD->calculateAnnualisedEquivalent(raw_amount,cumTp_u, ir);
550  if (anualised_amount>expReturns) {
551  expReturns=anualised_amount;
552  // if (ft == "con_highF" && regId == 11041){
553  // cout << thisYear << "/" << dc << "/" << priProducts[pp] << "/" << anualised_amount << "/" << cumTp_u << "/" << vHa << "/" << pw << endl;
554  // }
555  }
556  if(outputLevel >= OUTVL_ALL){
557  expReturnsDebug_pp.push_back(0.0);
558  expReturnsDebug_pp.push_back(hV_byFT);
559  expReturnsDebug_pp.push_back(finalHarvestFlag);
560  expReturnsDebug_pp.push_back(0.0);
561  expReturnsDebug_pp.push_back(pw);
562  expReturnsDebug_pp.push_back(cumTp_u);
563  expReturnsDebug_pp.push_back(vHa);
564  expReturnsDebug_pp.push_back(anualised_amount);
565  expReturnsDebug_pp.push_back(0);
566  }
567  expReturnsDebug_dc.push_back(expReturnsDebug_pp);
568  } // end each pp
569  expReturnsDebug_ft.push_back(expReturnsDebug_dc);
570  } // end dc
571  } else {
572  for (uint u=0; u<dClasses.size(); u++){
573  vector <vector <double> > expReturnsDebug_dc;
574  string dc = dClasses[u];
575  double vHa = gfd("vHa_exp",regId,ft,dc);
576  double finalHarvestFlag = gfd("finalHarvestFlag",regId,ft,dc);
577  double cumTp_u = gfd("cumTp_exp",regId,ft,dc);
578 
579  for (uint pp=0;pp<priProducts.size();pp++){
580  vector <double> expReturnsDebug_pp;
581  double pw = gpd("pw",regId,priProducts[pp]);
582  double pl = gpd("pl",regId,priProducts[pp]);
583  double pwor = gpd("pl",99999,priProducts[pp]);
584 
585  double hVol_byUPp = hV_byPrd.at(i).at(j).at(u).at(pp); // harvested volumes for this product, diameter class on this region and forest type
586 
587  //double raw_amount_old = pw*hv2fa* hVol_byUPp/hV_byFT; // old and wrong. it was in €/m^4
588  double raw_amount = finalHarvestFlag*pw*vHa* hVol_byUPp/hV_byFT; // now in €/ha if there is also thining volumes in hV_byFT, I underestimate expected returns !! TO.DO change it !! DONE, DONE...
589  /**
590  see @ModelData::calculateAnnualisedEquivalent
591  */
592  double anualised_amount = MD->calculateAnnualisedEquivalent(raw_amount,cumTp_u, ir); //comTp is on diamClasses, u here is on pDiamClasses
593  //cout << "reg|ft|dc|prd|raw amount|ann.amount|tp|hV|hVTot|pw|pl|pW|vHa|fHFlag;";
594  //cout << regId <<";"<< ft <<";"<< dc <<";" << priProducts[pp] <<";" << raw_amount <<";"<< anualised_amount<<";";
595  //cout << cumTp_u <<";"<< hVol_byUPp << ";" << hV_byFT << ";" << pw << ";" << pl << ";" << pwor << ";" << vHa << ";" << finalHarvestFlag << endl;
596  expReturns += anualised_amount;
597 
598  if(outputLevel >= OUTVL_ALL){
599  expReturnsDebug_pp.push_back(hVol_byUPp);
600  expReturnsDebug_pp.push_back(hV_byFT);
601  expReturnsDebug_pp.push_back(finalHarvestFlag);
602  expReturnsDebug_pp.push_back(finalHarvestFlag*hVol_byUPp/hV_byFT);
603  expReturnsDebug_pp.push_back(pw);
604  expReturnsDebug_pp.push_back(cumTp_u);
605  expReturnsDebug_pp.push_back(vHa);
606  expReturnsDebug_pp.push_back(MD->calculateAnnualisedEquivalent(finalHarvestFlag*pw*vHa,cumTp_u,ir));
607  expReturnsDebug_pp.push_back(1);
608  }
609  expReturnsDebug_dc.push_back(expReturnsDebug_pp);
610  } // end for each priProducts
611 
612  expReturnsDebug_ft.push_back(expReturnsDebug_dc);
613  //expReturnsPondCoef.push_back(expReturnsPondCoef_byPrd);
614  } // end for each dc
615  } // ending "it has been harvested" condition
616  double debug = expReturns;
617  sfd(expReturns,"expReturns",regId, ft,"",DATA_NOW,true);
618  cachedExpectedReturnsByFt.push_back(expReturns);
619  expReturnsDebug_region.push_back(expReturnsDebug_ft);
620  } // end foreach forest
621  MTHREAD->DO->expReturnsDebug.push_back(expReturnsDebug_region);
622 
623 
624  // PART 2: ALLOCATING THE HARVESTED AREA TO REGENERATION AREA BASED ON EXPECTED RETURNS
625 
626  // calculating freeArea at the end of the year and choosing the new regeneration area..
627  //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);
628  //if(scen("endVreg") ,
629  // regArea(i,essence,lambda,t) = freeArea(i,essence, lambda); // here we could introduce in/out area from other land usages
630  //else
631  // loop (i,
632  // loop( (essence,lambda),
633  // if ( expReturns(i,essence,lambda) = smax( (essence2,lambda2),expReturns(i,essence2,lambda2) ),
634  // regArea (i,essence,lambda,t) = sum( (essence2, lambda2), freeArea(i,essence2, lambda2) ) * mr;
635  // );
636  // );
637  // regArea(i,essence,lambda,t) = freeArea(i,essence, lambda)*(1-mr); // here we could introduce in/out area from other land usages
638  // );
639  double totalHarvestedArea = gfd("harvestedArea",regId,FT_ALL,DIAM_ALL);
640 
641  double maxExpReturns = *( max_element( cachedExpectedReturnsByFt.begin(), cachedExpectedReturnsByFt.end() ) );
642  bool foundMaxExpReturns = false;
643  for(uint j=0;j<fTypes.size();j++){
644  string ft = fTypes[j];
645  double harvestedAreaForThisFT = gfd("harvestedArea",regId,ft,DIAM_ALL);
646  if(regType == "fixed" || regType == "fromHrLevel"){
647  // regeneration area is the harvested area..
648  double harvestedArea = harvestedAreaForThisFT;
649  sfd(harvestedArea,"regArea",regId,ft,"",DATA_NOW,true);
650  } else {
651  // regeneration area is a mix between harvested area and the harvested area of te most profitting forest type..
652  double regArea = 0;
653  if (!foundMaxExpReturns && cachedExpectedReturnsByFt[j] == maxExpReturns){
654  // I use the foundMaxExpReturns for the unlikely event that two forest types have the same expected return to avoid overcounting of the area
655  regArea += totalHarvestedArea*mr;
656  foundMaxExpReturns = true;
657  }
658  double freq = rescaleFrequencies ? gfd("freq_norm",regId,ft,""):gfd("freq",regId,ft,""); // "probability of presence" for unmanaged forest, added 20140318
659  regArea += harvestedAreaForThisFT*(1-mr)*freq;
660  sfd(regArea,"regArea",regId,ft,"",DATA_NOW,true);
661  }
662  }// end of foreach forest type
663  double totalRegArea = gfd("regArea",regId,FT_ALL,DIAM_ALL);
664  } // end of each r2
665  //vector <vector < vector <vector <vector <double> > > > > expReturnsDebug = MTHREAD->DO->expReturnsDebug;
666  //cout << "bla" << endl;
667 
668 }
669 
670 void
671 ModelCore::computeInventary(){
672  msgOut(MSG_INFO, "Starting computing inventary available for this year..");
673  int thisYear = MTHREAD->SCD->getYear();
674 
675  // In(i,p_pr,t) = sum((u,lambda,essence),prov(u,essence,lambda,p_pr)*V(u,i,lambda,essence,t-1));
676  for(uint i=0;i<regIds2.size();i++){
677  int r2 = regIds2[i];
678  for(uint pp=0;pp<priProducts.size();pp++){
679  double in = 0;
680  for(uint ft=0;ft<fTypes.size();ft++){
681  for(uint dc=0;dc<dClasses.size();dc++){
682  double vol = dc?gfd("vol",r2,fTypes[ft],dClasses[dc],thisYear-1):0.;
683  in += app(priProducts[pp],fTypes[ft],dClasses[dc])*vol;
684  }
685  }
686  spd(in,"in",r2,priProducts[pp],thisYear,true);
687 
688  }
689  } // end of for each region
690 }
691 
692 void
693 ModelCore::cacheSettings(){
694  msgOut(MSG_INFO, "Cashing initial model settings..");
695  int currentYear = MTHREAD->SCD->getYear();
696 
697  MD = MTHREAD->MD;
698  firstYear = MD->getIntSetting("initialYear");
699  secondYear = firstYear+1;
700  thirdYear = firstYear+2;
701  WL2 = MD->getIntSetting("worldCodeLev2");
702  regIds2 = MD->getRegionIds(2);
703  priProducts = MD->getStringVectorSetting("priProducts");
704  secProducts = MD->getStringVectorSetting("secProducts");
705  allProducts = priProducts;
706  allProducts.insert( allProducts.end(), secProducts.begin(), secProducts.end() );
707  dClasses = MD->getStringVectorSetting("dClasses");
708  pDClasses; // production diameter classes: exclude the fist diameter class below 15 cm
709  pDClasses.insert(pDClasses.end(), dClasses.begin()+1, dClasses.end() );
710  fTypes= MD->getForTypeIds();
711  l2r = MD->getRegionIds();
712  regType = MTHREAD->MD->getStringSetting("regType"); // how the regeneration should be computed (exogenous, from hr, from allocation choises)
713  expType = MD->getDoubleSetting("expType");
714  rescaleFrequencies = MD->getBoolSetting("rescaleFrequencies");
715  if((expType<0 || expType>1) && expType != -1){
716  msgOut(MSG_CRITICAL_ERROR, "expType parameter must be between 1 (expectations) and 0 (adaptative) or -1 (fixed).");
717  }
718  mr = MD->getDoubleSetting("mr");
719 }
720 
721 /**
722 * Computing some fully exogenous parameters that require complex operations, e.g. cumulative time of passage or volume per hectare.
723 * This happen at the very beginning of the init period and after each simulated year
724 *
725 * It doesn't include tp and mort multipliers, but this could be added as now there is a regional versiopn of them and not just a pixel version.
726 */
727 void
728 ModelCore::computeCumulativeData(){
729 
730  msgOut(MSG_INFO, "Starting computing some cumulative values..");
731  int thisYear = MTHREAD->SCD->getYear();
732 
733  // debug
734  //cout << "cumTp and vHa by dc:" << endl;
735  //cout << "regId|ft|varName|0|15|25|35|45|55|65|75|85|95|150|" << endl;
736 
737  for(uint r2= 0; r2<regIds2.size();r2++){
738  int regId = regIds2[r2];
739  for(uint j=0;j<fTypes.size();j++){
740  string ft = fTypes[j];
741  // calculating the cumulative time of passage and the (cumulativelly generated) vHa for each diameter class (depending on forest owners diam growth expectations)
742  //loop(u$(ord(u)=1),
743  // cumTp(u,i,lambda,essence) = tp_u1(i,essence,lambda);
744  //);
745  //loop(u$(ord(u)>1),
746  // cumTp(u,i,lambda,essence) = cumTp(u-1,i,lambda,essence)+tp(u-1,i,lambda,essence);
747  //);
748  ////ceil(x) DNLP returns the smallest integer number greater than or equal to x
749  //loop( (u,i,lambda,essence),
750  // cumTp(u,i,lambda,essence) = ceil(cumTp(u,i,lambda,essence));
751  //);
752  /**
753  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.
754  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 ?
755 For compatibility with the GAMS code, a -1 value means using initial simulation tp values (fixed cumTp)."
756  */
757  vector <double> cumTp_temp; // cumulative time of passage to REACH a diameter class (tp is to LEAVE to the next one)
758  vector <double> vHa_temp; // volume at hectar by each diameter class [m^3/ha]
759  vector <double> cumAlive_temp; // cumulated alive rate to reach a given diameter class
760  vector <double> cumTp_exp_temp; // "expected" version of cumTp
761  vector <double> vHa_exp_temp; // "expected" version of vHa
762  vector <double> cumAlive_exp_temp; // "expected" version of cumMort
763 
764  MD->setErrorLevel(MSG_NO_MSG); // as otherwise on 2007 otherwise sfd() will complain that is filling multiple years (2006 and 2007)
765  for (uint u=0; u<dClasses.size(); u++){
766  string dc = dClasses[u];
767  double cumTp_u, cumTp_u_exp, cumTp_u_noExp, cumTp_u_fullExp;
768  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;
769  double tp_u, tp_exp;
770  double cumAlive_u, cumAlive_exp_u;
771 
772  if(u==0) {
773  // first diameter class.. expected and real values are the same (0)
774  cumTp_u = 0.;
775  vHa_u = 0.;
776  cumAlive_u = 1.;
777  cumTp_temp.push_back(cumTp_u);
778  cumTp_exp_temp.push_back(cumTp_u);
779  vHa_temp.push_back(vHa_u);
780  vHa_exp_temp.push_back(vHa_u);
781  cumAlive_temp.push_back(cumAlive_u);
782  cumAlive_exp_temp.push_back(cumAlive_u);
783  sfd(cumTp_u,"cumTp",regId,ft,dc,DATA_NOW,true);
784  sfd(cumTp_u,"cumTp_exp",regId,ft,dc,DATA_NOW,true);
785  sfd(vHa_u, "vHa",regId,ft,dc,DATA_NOW,true);
786  sfd(vHa_u, "vHa_exp",regId,ft,dc,DATA_NOW,true);
787  sfd(cumAlive_u,"cumAlive",regId,ft,dc,DATA_NOW,true);
788  sfd(cumAlive_u,"cumAlive_exp",regId,ft,dc,DATA_NOW,true);
789  } else {
790  // other diameter classes.. first dealing with real values and then with expected ones..
791  // real values..
792  cumTp_u = cumTp_temp[u-1] + gfd("tp",regId,ft,dClasses[u-1],thisYear); // 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
793  if (u==1){
794  vHa_u = gfd("entryVolHa",regId,ft,"",thisYear);
795  mort = 0.; // not info about mortality first diameter class ("00")
796  } else {
797  mort = 1-pow(1-gfd("mortCoef",regId,ft,dClasses[u-1],thisYear),gfd("tp",regId,ft,dClasses[u-1],thisYear)); // mortality of the previous diameter class
798  beta = gfd("betaCoef",regId,ft,dc, thisYear);
799  vHa_u = vHa_temp[u-1]*beta*(1-mort);
800  }
801  cumAlive_u = max(0.,cumAlive_temp[u-1]*(1-mort));
802  cumAlive_temp.push_back(cumAlive_u);
803  cumTp_temp.push_back(cumTp_u);
804  vHa_temp.push_back(vHa_u);
805  sfd(cumTp_u,"cumTp",regId,ft,dc,DATA_NOW,true);
806  sfd(vHa_u,"vHa",regId,ft,dc,DATA_NOW,true);
807  sfd(cumAlive_u,"cumAlive",regId,ft,dc,DATA_NOW,true);
808 
809  // expected values..
810  if (expType == -1){
811  cumTp_u_exp = cumTp_exp_temp[u-1]+gfd("tp",regId,ft,dClasses[u-1],firstYear); // 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
812  cumTp_exp_temp.push_back(cumTp_u_exp);
813  if(u==1) {
814  vHa_u_exp = gfd("entryVolHa",regId,ft,"",firstYear);
815  mort_exp = 0.; // not info about mortality first diameter class ("00")
816  } else {
817  mort_exp = 1-pow(1-gfd("mortCoef",regId,ft,dClasses[u-1], firstYear),gfd("tp",regId,ft,dClasses[u-1],firstYear)) ; // mortality rate of previous diameter class
818  beta_exp = gfd("betaCoef",regId,ft,dc, firstYear);
819  vHa_u_exp = vHa_exp_temp[u-1]*beta_exp*(1-mort_exp);
820  }
821  } else {
822  cumTp_u_noExp = cumTp_exp_temp[u-1]+gfd("tp",regId,ft,dClasses[u-1]);
823  cumTp_u_fullExp = cumTp_exp_temp[u-1]+gfd("tp",regId,ft,dClasses[u-1],thisYear+ceil(cumTp_exp_temp[u-1])); // 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
824  cumTp_u_exp = cumTp_u_fullExp*expType+cumTp_u_noExp*(1-expType);
825  cumTp_exp_temp.push_back(cumTp_u_exp);
826  if(u==1) {
827  vHa_u_noExp = gfd("entryVolHa",regId,ft,"",DATA_NOW);
828  vHa_u_fullExp = gfd("entryVolHa",regId,ft,"",thisYear+ceil(cumTp_u));
829  vHa_u_exp = vHa_u_fullExp*expType+vHa_u_noExp*(1-expType);
830  mort_exp = 0. ;
831  } else {
832  mort_noExp = 1-pow(1-gfd("mortCoef",regId,ft,dClasses[u-1], DATA_NOW),cumTp_exp_temp[u]-cumTp_exp_temp[u-1]);
833  mort_fullExp = 1-pow(1-gfd("mortCoef",regId,ft,dClasses[u-1],thisYear+ceil(cumTp_temp[u-1])),cumTp_exp_temp[u]-cumTp_exp_temp[u-1]); // mortality of the previous diameter class
834  beta_noExp = gfd("betaCoef",regId,ft,dc, DATA_NOW);
835  beta_fullExp = gfd("betaCoef",regId,ft,dc, thisYear+ceil(cumTp_u));
836  mort_exp = mort_fullExp*expType+mort_noExp*(1-expType);
837  beta_exp = beta_fullExp*expType+beta_noExp*(1-expType);
838  vHa_u_exp = vHa_exp_temp[u-1]*beta_exp*(1-mort_exp);
839  }
840  }
841  vHa_exp_temp.push_back(vHa_u_exp);
842  cumAlive_exp_u = max(0.,cumAlive_exp_temp[u-1]*(1-mort_exp));
843  cumAlive_exp_temp.push_back(cumAlive_exp_u);
844  sfd(cumTp_u_exp,"cumTp_exp",regId,ft,dc,DATA_NOW,true);
845  sfd(vHa_u_exp, "vHa_exp",regId,ft,dc,DATA_NOW,true);
846  sfd(cumAlive_exp_u,"cumAlive_exp",regId,ft,dc,DATA_NOW,true);
847  //sfd(cumMort_u_exp, "cumMort_exp",regId,ft,dc,DATA_NOW,true);
848 
849  //cout << "********" << endl;
850  //cout << "dc: " << dClasses[u] << endl ;
851  //cout << "mort: " << mort << endl;
852  //cout << "mort_exp: " << mort_exp << endl;
853  //cout << "cumAlive: " << cumAlive_u << endl;
854  //cout << "cumAlive_exp: " << cumAlive_exp_u << endl;
855 
856 
857  }
858 
859  } // end of each diam class
860 
861 
862  // // debug stuff on vHa
863  // cout << regId << "|" << ft << "|cumTp_temp|";
864  // for (uint u=0; u<dClasses.size(); u++){
865  // cout << cumTp_temp.at(u)<<"|";
866  // }
867  // cout << endl;
868  // cout << regId << "|" << ft << "|cumTp_exp|";
869  // for (uint u=0; u<dClasses.size(); u++){
870  // cout << cumTp_exp_temp.at(u)<<"|";
871  // }
872  // cout << endl;
873  // cout << regId << "|" << ft << "|vHa_temp|";
874  // for (uint u=0; u<dClasses.size(); u++){
875  // cout << vHa_temp.at(u)<<"|";
876  // }
877  // cout << endl;
878  // cout << regId << "|" << ft << "|vHa_exp|";
879  // for (uint u=0; u<dClasses.size(); u++){
880  // cout << vHa_exp_temp.at(u)<<"|";
881  // }
882  // cout << endl;
883 
884 
885  } // end of each ft
886  } // end of each region
887  MD->setErrorLevel(MSG_ERROR);
888 }
889 
890 
891 
892 /**
893 This function take for each region the difference for each forest type between the harvested area and the new regeneration one and apply such delta to each pixel of the region proportionally to the area that it already hosts.
894 */
895 void
896 ModelCore::updateMapAreas(){
897 
898  msgOut(MSG_INFO, "Updating map areas..");
899  map<int,double> forestArea; // foresta area by each region
900  pair<int,double > forestAreaPair;
901  int thisYear = MTHREAD->SCD->getYear();
902  vector<int> l2Regions = MTHREAD->MD->getRegionIds(2, true);
903  vector <string> fTypes = MTHREAD->MD->getForTypeIds();
904  int nFTypes = fTypes.size();
905  int nL2Regions = l2Regions.size();
906  for(uint i=0;i<nL2Regions;i++){
907  int regId = l2Regions[i];
908  vector<Pixel*> rpx = MTHREAD->GIS->getAllPlotsByRegion(regId);
909  ModelRegion* reg = MTHREAD->MD->getRegion(regId);
910  for(uint j=0;j<nFTypes;j++){
911  string ft = fTypes[j];
912  double oldRegioForArea;
913  if(thisYear <= firstYear+1) {
914  oldRegioForArea = reg->getValue("forArea_"+ft)/10000;
915  } else {
916  oldRegioForArea = gfd("forArea",regId,ft,DIAM_ALL,thisYear-1);
917  }
918  //oldRegioForArea = reg->getValue("forArea_"+ft)/10000;
919  //double debug = gfd("forArea",regId,ft,DIAM_ALL,thisYear-1);
920  //double debug_diff = oldRegioForArea - debug;
921  //cout << thisYear << ";" << regId << ";" << ft << ";" << oldRegioForArea << ";" << debug << ";" << debug_diff << endl;
922  double harvestedArea = gfd("harvestedArea",regId,ft,DIAM_ALL); //20140206
923  double regArea = gfd("regArea",regId,ft,DIAM_ALL); //20140206
924  double newRegioForArea = oldRegioForArea + regArea - harvestedArea;
925  sfd(newRegioForArea,"forArea",regId,ft,"",DATA_NOW, true);
926  for(uint z=0;z<rpx.size();z++){
927  double oldValue = rpx[z]->getDoubleValue("forArea_"+ft,true);
928  double ratio = newRegioForArea/oldRegioForArea;
929  double newValue = oldValue*ratio;
930  rpx[z]->changeValue("forArea_"+ft, newValue);
931  }
932 
933  }
934  }
935 }
936 
937 
938 
Scheduler::advanceYear
void advanceYear()
Definition: Scheduler.h:51
ModelCore::thirdYear
int thirdYear
Definition: ModelCore.h:78
MSG_ERROR
Print an ERROR message, but don&#39;t stop the model.
Definition: BaseClass.h:61
ModelCore::runInitPeriod
void runInitPeriod()
Definition: ModelCore.cpp:50
DATA_NOW
The required data is for the current year.
Definition: BaseClass.h:73
ModelCore::computeCumulativeData
void computeCumulativeData()
computes cumTp, vHa, cumTp_exp, vHa_exp,
Definition: ModelCore.cpp:728
ModelData::getIntSetting
int getIntSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1105
ModelCore::priProducts
vector< string > priProducts
Definition: ModelCore.h:81
ModelCore::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: ModelCore.h:69
ModelData::getBoolSetting
bool getBoolSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1117
MSG_NO_MSG
Do not actually output any message.
Definition: BaseClass.h:57
ModelCore::regIds2
vector< int > regIds2
Definition: ModelCore.h:80
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
ModelCore::rescaleFrequencies
bool rescaleFrequencies
Definition: ModelCore.h:93
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
ModelCore::dClasses
vector< string > dClasses
Definition: ModelCore.h:84
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
OUTVL_ALL
Output verbosity level print everything.
Definition: BaseClass.h:89
ModelCore::runManagementModule
void runManagementModule()
computes regArea and expectedReturns
Definition: ModelCore.cpp:484
ModelCore::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: ModelCore.h:68
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
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
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
ModelData::getStringVectorSetting
vector< string > getStringVectorSetting(const string &name_h, int reg=WORLD) const
Definition: ModelData.cpp:1129
ModelCore::secondYear
int secondYear
Definition: ModelCore.h:77
ModelCore::l2r
vector< vector< int > > l2r
Definition: ModelCore.h:87
ModelCore::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: ModelCore.h:67
ModelData::getStringSetting
string getStringSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1113
ModelCore::hV_byPrd
vector< vector< vector< vector< double > > > > hV_byPrd
Definition: ModelCore.h:91
ModelCore::app
bool app(const string &prod_h, const string &forType_h, const string &dClass_h) const
Definition: ModelCore.h:70
MSG_CRITICAL_ERROR
Print an error message and stop the model.
Definition: BaseClass.h:62
ThreadManager::DO
Output * DO
data output
Definition: ThreadManager.h:76
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
ModelCore::expType
double expType
Definition: ModelCore.h:89
ModelCore::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: ModelCore.h:66
ThreadManager::OPT
Ipopt::SmartPtr< Ipopt::TNLP > OPT
Market optimisation.
Definition: ThreadManager.h:81
FT_ALL
#define FT_ALL
All forest types.
Definition: BaseClass.h:166
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::getForTypeIds
vector< string > getForTypeIds(bool all=false)
By default it doesn&#39;t return forTypes used only as input.
Definition: ModelData.cpp:430
ModelCore::initMarketModule
void initMarketModule()
computes st and pw for second year and several needed-only-at-t0-vars for the market module ...
Definition: ModelCore.cpp:93
ModelCore::pDClasses
vector< string > pDClasses
Definition: ModelCore.h:85
ModelCore::runBiologicalModule
void runBiologicalModule()
computes hV, hArea and new vol at end of year
Definition: ModelCore.cpp:363
ModelData::setErrorLevel
void setErrorLevel(int errorLevel_h)
Definition: ModelData.h:143
ModelCore::mr
double mr
Definition: ModelCore.h:90
ModelData::getRegionIds
vector< int > getRegionIds(int level_h, bool excludeResidual=true)
Definition: ModelData.cpp:366
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
ModelData::getDoubleSetting
double getDoubleSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1109
ModelCore::firstYear
int firstYear
Definition: ModelCore.h:76
ModelCore::secProducts
vector< string > secProducts
Definition: ModelCore.h:82
ModelCore::fTypes
vector< string > fTypes
Definition: ModelCore.h:86
Output::expReturnsDebug
vector< vector< vector< vector< vector< double > > > > > expReturnsDebug
l2_region, for type, d.c., pr prod, variable name
Definition: Output.h:75
ModelCore::WL2
int WL2
Definition: ModelCore.h:79
ModelCore::allProducts
vector< string > allProducts
Definition: ModelCore.h:83
MSG_INFO
Print an INFO message.
Definition: BaseClass.h:59
ModelCore::updateMapAreas
void updateMapAreas()
computes forArea_{ft}
Definition: ModelCore.cpp:896
ModelData::getRegion
ModelRegion * getRegion(int regId_h)
Definition: ModelData.cpp:346
ModelCore::runMarketModule
void runMarketModule()
computes st (supply total) and pw (weighted price). Optimisation inside.
Definition: ModelCore.cpp:211
DIAM_ALL
#define DIAM_ALL
All diameter classes.
Definition: BaseClass.h:157
ModelCore::MD
ModelData * MD
Definition: ModelCore.h:70
ModelCore::runSimulationYear
void runSimulationYear()
Definition: ModelCore.cpp:70
ModelCore::ModelCore
ModelCore(ThreadManager *MTHREAD_h)
Definition: ModelCore.cpp:37
ModelCore::cacheSettings
void cacheSettings()
just cache exogenous settings from ModelData
Definition: ModelCore.cpp:693
ModelCore::computeInventary
void computeInventary()
in=f(vol_t-1)
Definition: ModelCore.cpp:671
ModelCore::regType
string regType
Definition: ModelCore.h:88