FFSM++  1.1.0
French Forest Sector Model ++
Public Member Functions | Private Member Functions | Private Attributes | List of all members
ModelCore Class Reference

#include <ModelCore.h>

Inheritance diagram for ModelCore:
Inheritance graph
[legend]
Collaboration diagram for ModelCore:
Collaboration graph
[legend]

Public Member Functions

 ModelCore (ThreadManager *MTHREAD_h)
 
 ~ModelCore ()
 
void runInitPeriod ()
 
void runSimulationYear ()
 
void initMarketModule ()
 computes st and pw for second year and several needed-only-at-t0-vars for the market module More...
 
void runMarketModule ()
 computes st (supply total) and pw (weighted price). Optimisation inside. More...
 
void runBiologicalModule ()
 computes hV, hArea and new vol at end of year More...
 
void runManagementModule ()
 computes regArea and expectedReturns More...
 
void cacheSettings ()
 just cache exogenous settings from ModelData More...
 
void cachePixelExogenousData ()
 computes pixel level tp, meta and mort More...
 
void computeInventary ()
 in=f(vol_t-1) More...
 
void computeCumulativeData ()
 computes cumTp, vHa, cumTp_exp, vHa_exp, More...
 
void updateMapAreas ()
 computes forArea_{ft} More...
 
- Public Member Functions inherited from BaseClass
 BaseClass ()
 
 ~BaseClass ()
 
void msgOut (const int &msgCode_h, const string &msg_h, const bool &refreshGUI_h=true) const
 Overloaded function to print the output log. More...
 
void msgOut (const int &msgCode_h, const int &msg_h, const bool &refreshGUI_h=true) const
 Overloaded function to print the output log. More...
 
void msgOut (const int &msgCode_h, const double &msg_h, const bool &refreshGUI_h=true) const
 Overloaded function to print the output log. More...
 
int s2i (const string &string_h) const
 string to integer conversion More...
 
double s2d (const string &string_h) const
 string to double conversion More...
 
double s2d (const string &string_h, const bool &replaceComma) const
 string to double conversion More...
 
bool s2b (const string &string_h) const
 string to bool conversion More...
 
string i2s (const int &int_h) const
 integer to string conversion More...
 
string d2s (const double &double_h) const
 double to string conversion More...
 
string b2s (const bool &bool_h) const
 bool to string conversion More...
 
vector< int > s2i (const vector< string > &string_h) const
 string to integer conversion (vector) More...
 
vector< double > s2d (const vector< string > &string_h, const bool &replaceComma=false) const
 string to double conversion (vector) More...
 
vector< bool > s2b (const vector< string > &string_h) const
 string to bool conversion (vector) More...
 
vector< string > i2s (const vector< int > &int_h) const
 integer to string conversion (vector) More...
 
vector< string > d2s (const vector< double > &double_h) const
 double to string conversion (vector) More...
 
vector< string > b2s (const vector< bool > &bool_h) const
 bool to string conversion (vector) More...
 
int getType (const string &type_h) const
 Return a type according to enum TYPE_* from a string (eg: "string" -> TYPE_STRING (2)) More...
 
void refreshGUI () const
 Ping to periodically return the control to the GUI. More...
 
template<typename T >
string toString (const T &x) const
 
template<typename T >
T stringTo (const std::string &s) const
 
int vSum (const vector< int > &vector_h) const
 
double vSum (const vector< double > &vector_h) const
 
int vSum (const vector< vector< int > > &vector_h) const
 
double vSum (const vector< vector< double > > &vector_h) const
 
void tokenize (const string &str, vector< string > &tokens, const string &delimiter=" ") const
 Tokenize a string using a delimiter (default is space) More...
 
void untokenize (string &str, vector< string > &tokens, const string &delimiter=" ") const
 
template<typename K , typename V >
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. More...
 
template<typename K , typename V >
void changeMapValue (map< K, V > &mymap, const K &key, const V &value, const int &error_level=MSG_CRITICAL_ERROR)
 Change the value stored in a map given the key and the new value. More...
 
template<typename K , typename V >
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. More...
 
template<typename K , typename V >
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. More...
 
template<typename K , typename V >
void resetMapValues (map< K, V > &mymap, const V &value)
 Reset all values stored in a map to the specified one. More...
 
template<typename K , typename V >
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. More...
 
template<typename T >
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) More...
 
template<typename V >
void debugMap (const map< iisskey, V > &mymap)
 Debug a map. More...
 
template<typename K , typename V >
void debugMap (const map< K, V > &mymap, const K &key)
 
template<typename K >
int getMaxPos (const vector< K > &v)
 Returns the position of the maximum element in the vector (the last one in case of multiple equivalent maxima) More...
 
template<typename K >
int getMinPos (const vector< K > &v)
 Returns the position of the minimum element in the vector (the first one in case of multiple equivalent minima) More...
 
template<typename K >
K getMax (const vector< K > &v)
 Returns the value of the maximum element in the vector (the last one in case of multiple equivalent maxima) More...
 
template<typename K >
K getMin (const vector< K > &v)
 Returns the value of the minimum element in the vector (the first one in case of multiple equivalent minima) More...
 
template<typename K >
double getAvg (const vector< K > &v)
 Returns the average of the elements in the vector. More...
 
template<typename K >
double getSd (const vector< K > &v, bool sample=true)
 
template<typename K >
int getPos (const K &element, const vector< K > &v, const int &msgCode_h=MSG_CRITICAL_ERROR)
 
template<typename K >
bool inVector (const K &element, const vector< K > &v)
 
double normSample (const double &avg, const double &stdev, const double &minval=NULL, const double &maxval=NULL) const
 Sample from a normal distribution with bounds. Slower (double time, but still you see the diff only after milion of loops). More...
 
template<typename K >
K normSample (normal_distribution< K > &d, std::mt19937 &gen, const K &minval=NULL, const K &maxval=NULL) const
 Sample from a normal distribution with bounds. Faster (half time) as the normal_distribution is made only once. More...
 
template<typename T >
std::string toString (const T &x) const
 

Private Member Functions

double gpd (const string &type_h, const int &regId_h, const string &prodId_h, const int &year=DATA_NOW, const string &freeDim_h="") const
 
double gfd (const string &type_h, const int &regId_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW) const
 
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
 
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
 
bool app (const string &prod_h, const string &forType_h, const string &dClass_h) const
 

Private Attributes

ModelData * MD
 
int firstYear
 
int secondYear
 
int thirdYear
 
int WL2
 
vector< int > regIds2
 
vector< string > priProducts
 
vector< string > secProducts
 
vector< string > allProducts
 
vector< string > dClasses
 
vector< string > pDClasses
 
vector< string > fTypes
 
vector< vector< int > > l2r
 
string regType
 
double expType
 
double mr
 
vector< vector< vector< vector< double > > > > hV_byPrd
 
bool rescaleFrequencies
 

Additional Inherited Members

- Protected Attributes inherited from BaseClass
ThreadManager * MTHREAD
 Pointer to the Thread manager. More...
 

Detailed Description

Definition at line 43 of file ModelCore.h.

Constructor & Destructor Documentation

ModelCore ( ThreadManager *  MTHREAD_h)

Definition at line 37 of file ModelCore.cpp.

37  {
38  MTHREAD = MTHREAD_h;
39 }
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
~ModelCore ( )

Definition at line 41 of file ModelCore.cpp.

41  {
42 
43 }

Member Function Documentation

bool app ( const string &  prod_h,
const string &  forType_h,
const string &  dClass_h 
) const
inlineprivate

Definition at line 70 of file ModelCore.h.

Referenced by computeInventary(), runBiologicalModule(), and runManagementModule().

70 {return MTHREAD->MD->assessProdPossibility(prod_h, forType_h, dClass_h);};
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
ModelData::assessProdPossibility
const bool assessProdPossibility(const string &prod_h, const string &forType_h, const string &dClass_h)
A simple function to assess if a specified product can be made by a certain forest type and diameter ...
Definition: ModelData.cpp:441

Here is the caller graph for this function:

void cachePixelExogenousData ( )

computes pixel level tp, meta and mort

void cacheSettings ( )

just cache exogenous settings from ModelData

Definition at line 693 of file ModelCore.cpp.

Referenced by runInitPeriod().

693  {
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 }
ModelCore::thirdYear
int thirdYear
Definition: ModelCore.h:78
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
ModelData::getBoolSetting
bool getBoolSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1117
ModelCore::regIds2
vector< int > regIds2
Definition: ModelCore.h:80
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ModelCore::rescaleFrequencies
bool rescaleFrequencies
Definition: ModelCore.h:93
ModelCore::dClasses
vector< string > dClasses
Definition: ModelCore.h:84
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
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
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
ModelData::getStringSetting
string getStringSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1113
MSG_CRITICAL_ERROR
Print an error message and stop the model.
Definition: BaseClass.h:62
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
ModelCore::expType
double expType
Definition: ModelCore.h:89
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::pDClasses
vector< string > pDClasses
Definition: ModelCore.h:85
ModelCore::mr
double mr
Definition: ModelCore.h:90
ModelData::getRegionIds
vector< int > getRegionIds(int level_h, bool excludeResidual=true)
Definition: ModelData.cpp:366
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
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::MD
ModelData * MD
Definition: ModelCore.h:70
ModelCore::regType
string regType
Definition: ModelCore.h:88

Here is the call graph for this function:

Here is the caller graph for this function:

void computeCumulativeData ( )

computes cumTp, vHa, cumTp_exp, vHa_exp,

Computing some fully exogenous parameters that require complex operations, e.g. cumulative time of passage or volume per hectare. This happen at the very beginning of the init period and after each simulated year

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. 

    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.
    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 ?

For compatibility with the GAMS code, a -1 value means using initial simulation tp values (fixed cumTp)."

Definition at line 728 of file ModelCore.cpp.

Referenced by runInitPeriod(), and runSimulationYear().

728  {
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 }
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
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
MSG_NO_MSG
Do not actually output any message.
Definition: BaseClass.h:57
ModelCore::regIds2
vector< int > regIds2
Definition: ModelCore.h:80
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ModelCore::dClasses
vector< string > dClasses
Definition: ModelCore.h:84
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
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
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
ModelCore::expType
double expType
Definition: ModelCore.h:89
ModelData::setErrorLevel
void setErrorLevel(int errorLevel_h)
Definition: ModelData.h:143
ModelCore::firstYear
int firstYear
Definition: ModelCore.h:76
ModelCore::fTypes
vector< string > fTypes
Definition: ModelCore.h:86
MSG_INFO
Print an INFO message.
Definition: BaseClass.h:59
ModelCore::MD
ModelData * MD
Definition: ModelCore.h:70

Here is the call graph for this function:

Here is the caller graph for this function:

void computeInventary ( )

in=f(vol_t-1)

Definition at line 671 of file ModelCore.cpp.

Referenced by runInitPeriod(), and runSimulationYear().

671  {
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 }
ModelCore::priProducts
vector< string > priProducts
Definition: ModelCore.h:81
ModelCore::regIds2
vector< int > regIds2
Definition: ModelCore.h:80
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ModelCore::dClasses
vector< string > dClasses
Definition: ModelCore.h:84
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
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
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
ModelCore::app
bool app(const string &prod_h, const string &forType_h, const string &dClass_h) const
Definition: ModelCore.h:70
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
ModelCore::fTypes
vector< string > fTypes
Definition: ModelCore.h:86
MSG_INFO
Print an INFO message.
Definition: BaseClass.h:59

Here is the call graph for this function:

Here is the caller graph for this function:

double gfd ( const string &  type_h,
const int &  regId_h,
const string &  forType_h,
const string &  freeDim_h,
const int &  year = DATA_NOW 
) const
inlineprivate

Definition at line 67 of file ModelCore.h.

Referenced by computeCumulativeData(), computeInventary(), runBiologicalModule(), runManagementModule(), and updateMapAreas().

67 {return MTHREAD->MD->getForData(type_h, regId_h, forType_h, freeDim_h, year);};
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
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

Here is the call graph for this function:

Here is the caller graph for this function:

double gpd ( const string &  type_h,
const int &  regId_h,
const string &  prodId_h,
const int &  year = DATA_NOW,
const string &  freeDim_h = "" 
) const
inlineprivate

Definition at line 66 of file ModelCore.h.

Referenced by initMarketModule(), runBiologicalModule(), runManagementModule(), and runMarketModule().

66 {return MTHREAD->MD->getProdData(type_h, regId_h, prodId_h, year, freeDim_h);};
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
ModelData::getProdData
const double getProdData(const string &type_h, const int &regId_h, const string &prodId_h, const int &year=DATA_NOW, const string &freeDim_h="")
Definition: ModelData.cpp:1216

Here is the call graph for this function:

Here is the caller graph for this function:

void initMarketModule ( )

computes st and pw for second year and several needed-only-at-t0-vars for the market module

Definition at line 93 of file ModelCore.cpp.

Referenced by runInitPeriod().

93  {
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 }
ModelCore::priProducts
vector< string > priProducts
Definition: ModelCore.h:81
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
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
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
ModelCore::secondYear
int secondYear
Definition: ModelCore.h:77
ModelCore::l2r
vector< vector< int > > l2r
Definition: ModelCore.h:87
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
ModelCore::firstYear
int firstYear
Definition: ModelCore.h:76
ModelCore::secProducts
vector< string > secProducts
Definition: ModelCore.h:82
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

Here is the call graph for this function:

Here is the caller graph for this function:

void runBiologicalModule ( )

computes hV, hArea and new vol at end of year

Todo:
Harvest volumes from death trees

Definition at line 363 of file ModelCore.cpp.

Referenced by runInitPeriod(), and runSimulationYear().

363  {
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 }
DATA_NOW
The required data is for the current year.
Definition: BaseClass.h:73
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
ModelCore::regIds2
vector< int > regIds2
Definition: ModelCore.h:80
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ModelCore::dClasses
vector< string > dClasses
Definition: ModelCore.h:84
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
ModelCore::secondYear
int secondYear
Definition: ModelCore.h:77
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
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
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
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
ModelCore::fTypes
vector< string > fTypes
Definition: ModelCore.h:86
MSG_INFO
Print an INFO message.
Definition: BaseClass.h:59
ModelCore::regType
string regType
Definition: ModelCore.h:88

Here is the call graph for this function:

Here is the caller graph for this function:

void runInitPeriod ( )

IMPORTANT NOTE: Volumes in Mm^3, Areas in the model in Ha (10000 m^2), in the layers in m^2

Some importan notes: V (volumes) -> at the end of the year In (inventary) -> at the beginning of the year Volumes are in Mm^3, Areas in the model in Ha (10000 m^2), in the layers in m^2

Definition at line 50 of file ModelCore.cpp.

Referenced by Init::setInitLevel3().

50  {
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 }
Scheduler::advanceYear
void advanceYear()
Definition: Scheduler.h:51
ModelCore::computeCumulativeData
void computeCumulativeData()
computes cumTp, vHa, cumTp_exp, vHa_exp,
Definition: ModelCore.cpp:728
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ModelCore::runManagementModule
void runManagementModule()
computes regArea and expectedReturns
Definition: ModelCore.cpp:484
ThreadManager::SCD
Scheduler * SCD
the scheduler object (simulation-loops scheduler)
Definition: ThreadManager.h:75
ThreadManager::DO
Output * DO
data output
Definition: ThreadManager.h:76
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::runBiologicalModule
void runBiologicalModule()
computes hV, hArea and new vol at end of year
Definition: ModelCore.cpp:363
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
ModelCore::updateMapAreas
void updateMapAreas()
computes forArea_{ft}
Definition: ModelCore.cpp:896
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

Here is the call graph for this function:

Here is the caller graph for this function:

void runManagementModule ( )

computes regArea and expectedReturns

see ::calculateAnnualisedEquivalent

Definition at line 484 of file ModelCore.cpp.

Referenced by runInitPeriod(), and runSimulationYear().

484  {
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 }
DATA_NOW
The required data is for the current year.
Definition: BaseClass.h:73
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
ModelCore::regIds2
vector< int > regIds2
Definition: ModelCore.h:80
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ModelCore::rescaleFrequencies
bool rescaleFrequencies
Definition: ModelCore.h:93
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
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
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
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
ThreadManager::DO
Output * DO
data output
Definition: ThreadManager.h:76
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
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
FT_ALL
#define FT_ALL
All forest types.
Definition: BaseClass.h:166
ModelCore::mr
double mr
Definition: ModelCore.h:90
ModelData::getDoubleSetting
double getDoubleSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1109
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
MSG_INFO
Print an INFO message.
Definition: BaseClass.h:59
DIAM_ALL
#define DIAM_ALL
All diameter classes.
Definition: BaseClass.h:157
ModelCore::MD
ModelData * MD
Definition: ModelCore.h:70
ModelCore::regType
string regType
Definition: ModelCore.h:88

Here is the call graph for this function:

Here is the caller graph for this function:

void runMarketModule ( )

computes st (supply total) and pw (weighted price). Optimisation inside.

Definition at line 211 of file ModelCore.cpp.

Referenced by runSimulationYear().

211  {
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 }
DATA_NOW
The required data is for the current year.
Definition: BaseClass.h:73
ModelCore::priProducts
vector< string > priProducts
Definition: ModelCore.h:81
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
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
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
ModelData::getStringSetting
string getStringSetting(const string &name_h, int position=0, int reg=WORLD) const
Definition: ModelData.cpp:1113
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::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
ModelCore::secProducts
vector< string > secProducts
Definition: ModelCore.h:82
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

Here is the call graph for this function:

Here is the caller graph for this function:

void runSimulationYear ( )

Definition at line 70 of file ModelCore.cpp.

Referenced by Scheduler::run().

70  {
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 }
ModelCore::computeCumulativeData
void computeCumulativeData()
computes cumTp, vHa, cumTp_exp, vHa_exp,
Definition: ModelCore.cpp:728
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ModelCore::runManagementModule
void runManagementModule()
computes regArea and expectedReturns
Definition: ModelCore.cpp:484
ThreadManager::SCD
Scheduler * SCD
the scheduler object (simulation-loops scheduler)
Definition: ThreadManager.h:75
ThreadManager::DO
Output * DO
data output
Definition: ThreadManager.h:76
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
ModelCore::runBiologicalModule
void runBiologicalModule()
computes hV, hArea and new vol at end of year
Definition: ModelCore.cpp:363
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
ModelCore::updateMapAreas
void updateMapAreas()
computes forArea_{ft}
Definition: ModelCore.cpp:896
ModelCore::runMarketModule
void runMarketModule()
computes st (supply total) and pw (weighted price). Optimisation inside.
Definition: ModelCore.cpp:211
ModelCore::computeInventary
void computeInventary()
in=f(vol_t-1)
Definition: ModelCore.cpp:671

Here is the call graph for this function:

Here is the caller graph for this function:

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
inlineprivate

Definition at line 69 of file ModelCore.h.

Referenced by computeCumulativeData(), runBiologicalModule(), runManagementModule(), and updateMapAreas().

69 {MTHREAD->MD->setForData(value_h, type_h, regId_h, forType_h, freeDim_h, year, allowCreate);};
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
ModelData::setForData
void setForData(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)
Definition: ModelData.cpp:1412

Here is the call graph for this function:

Here is the caller graph for this function:

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
inlineprivate

Definition at line 68 of file ModelCore.h.

Referenced by computeInventary(), initMarketModule(), and runMarketModule().

68 {MTHREAD->MD->setProdData(value_h, type_h, regId_h, prodId_h, year, allowCreate, freeDim_h);};
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
ModelData::setProdData
void setProdData(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="")
Definition: ModelData.cpp:1352

Here is the call graph for this function:

Here is the caller graph for this function:

void updateMapAreas ( )

computes forArea_{ft}

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.

Definition at line 896 of file ModelCore.cpp.

Referenced by runInitPeriod(), and runSimulationYear().

896  {
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 }
DATA_NOW
The required data is for the current year.
Definition: BaseClass.h:73
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
BaseClass::MTHREAD
ThreadManager * MTHREAD
Pointer to the Thread manager.
Definition: BaseClass.h:467
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
ThreadManager::MD
ModelData * MD
the model data object
Definition: ThreadManager.h:72
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
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
Scheduler::getYear
int getYear()
Definition: Scheduler.h:49
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
ModelData::getRegionIds
vector< int > getRegionIds(int level_h, bool excludeResidual=true)
Definition: ModelData.cpp:366
ModelCore::firstYear
int firstYear
Definition: ModelCore.h:76
ModelCore::fTypes
vector< string > fTypes
Definition: ModelCore.h:86
MSG_INFO
Print an INFO message.
Definition: BaseClass.h:59
ModelData::getRegion
ModelRegion * getRegion(int regId_h)
Definition: ModelData.cpp:346
DIAM_ALL
#define DIAM_ALL
All diameter classes.
Definition: BaseClass.h:157

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

vector<string> allProducts
private

Definition at line 83 of file ModelCore.h.

Referenced by cacheSettings(), initMarketModule(), and runMarketModule().

vector<string> dClasses
private

Definition at line 84 of file ModelCore.h.

Referenced by cacheSettings(), computeCumulativeData(), computeInventary(), runBiologicalModule(), and runManagementModule().

double expType
private

Definition at line 89 of file ModelCore.h.

Referenced by cacheSettings(), and computeCumulativeData().

int firstYear
private

Definition at line 76 of file ModelCore.h.

Referenced by cacheSettings(), computeCumulativeData(), initMarketModule(), and updateMapAreas().

vector<string> fTypes
private

Definition at line 86 of file ModelCore.h.

Referenced by cacheSettings(), computeCumulativeData(), computeInventary(), runBiologicalModule(), runManagementModule(), and updateMapAreas().

vector< vector < vector < vector <double> > > > hV_byPrd
private

Definition at line 91 of file ModelCore.h.

Referenced by runBiologicalModule(), and runManagementModule().

vector<vector <int> > l2r
private

Definition at line 87 of file ModelCore.h.

Referenced by cacheSettings(), and initMarketModule().

ModelData* MD
private

Definition at line 70 of file ModelCore.h.

Referenced by cacheSettings(), computeCumulativeData(), and runManagementModule().

double mr
private

Definition at line 90 of file ModelCore.h.

Referenced by cacheSettings(), and runManagementModule().

vector<string> pDClasses
private

Definition at line 85 of file ModelCore.h.

Referenced by cacheSettings().

vector<string> priProducts
private

Definition at line 81 of file ModelCore.h.

Referenced by cacheSettings(), computeInventary(), initMarketModule(), runBiologicalModule(), runManagementModule(), and runMarketModule().

vector<int> regIds2
private

Definition at line 80 of file ModelCore.h.

Referenced by cacheSettings(), computeCumulativeData(), computeInventary(), initMarketModule(), runBiologicalModule(), runManagementModule(), and runMarketModule().

string regType
private

Definition at line 88 of file ModelCore.h.

Referenced by cacheSettings(), runBiologicalModule(), and runManagementModule().

bool rescaleFrequencies
private

Definition at line 93 of file ModelCore.h.

Referenced by cacheSettings(), and runManagementModule().

int secondYear
private

Definition at line 77 of file ModelCore.h.

Referenced by cacheSettings(), initMarketModule(), and runBiologicalModule().

vector<string> secProducts
private

Definition at line 82 of file ModelCore.h.

Referenced by cacheSettings(), initMarketModule(), and runMarketModule().

int thirdYear
private

Definition at line 78 of file ModelCore.h.

Referenced by cacheSettings().

int WL2
private

Definition at line 79 of file ModelCore.h.

Referenced by cacheSettings(), initMarketModule(), and runMarketModule().

The documentation for this class was generated from the following files: