FFSM++  1.1.0
French Forest Sector Model ++
Opt.h
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 #ifndef STDOPT_H
23 #define STDOPT_H
24 
25 
26 #include "IpTNLP.hpp"
27 #include <adolc.h>
28 #include <adolc_sparse.h>
29 
30 //regmas headers
31 #include "BaseClass.h"
32 #include "ThreadManager.h"
33 #include "ModelData.h"
34 
35 #define tag_f 1
36 #define tag_g 2
37 #define tag_L 3
38 #define HPOFF 30 //original: 30
39 
40 /// Class containing the optimization problem (the matrix and its methods)
41 
42 /**
43 
44 @author Antonello Lobianco
45 */
46 
47 using namespace Ipopt;
48 
49 struct constrain;
50 struct endvar;
51 
52 class Opt: public BaseClass, public TNLP{
53 public:
54  Opt(ThreadManager* MTHREAD_h); ///< Constructor
55  ~Opt();
56  /**@name Overloaded from TNLP */
57  //@{
58  /** Method to return some info about the nlp */
59  virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
60  Index& nnz_h_lag, IndexStyleEnum& index_style);
61  /** Method to return the bounds for my problem */
62  virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
63  Index m, Number* g_l, Number* g_u);
64 
65  /** Method to return the starting point for the algorithm */
66  virtual bool get_starting_point(Index n, bool init_x, Number* x,
67  bool init_z, Number* z_L, Number* z_U,
68  Index m, bool init_lambda,
69  Number* lambda);
70  /** Template to return the objective value */
71  template<class T> bool eval_obj(Index n, const T *x, T& obj_value);
72 
73  /** Template to compute contraints */
74  template<class T> bool eval_constraints(Index n, const T *x, Index m, T *g);
75 
76  /** Original method from Ipopt to return the objective value */
77  /** remains unchanged */
78  virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
79 
80  /** Original method from Ipopt to return the gradient of the objective */
81  /** remains unchanged */
82  virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
83 
84  /** Original method from Ipopt to return the constraint residuals */
85  /** remains unchanged */
86  virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
87 
88  /** Original method from Ipopt to return:
89  * 1) The structure of the jacobian (if "values" is NULL)
90  * 2) The values of the jacobian (if "values" is not NULL)
91  */
92  /** remains unchanged */
93  virtual bool eval_jac_g(Index n, const Number* x, bool new_x,
94  Index m, Index nele_jac, Index* iRow, Index *jCol,
95  Number* values);
96 
97  /** Original method from Ipopt to return:
98  * 1) The structure of the hessian of the lagrangian (if "values" is NULL)
99  * 2) The values of the hessian of the lagrangian (if "values" is not NULL)
100  */
101  /** remains unchanged */
102  virtual bool eval_h(Index n, const Number* x, bool new_x,
103  Number obj_factor, Index m, const Number* lambda,
104  bool new_lambda, Index nele_hess, Index* iRow,
105  Index* jCol, Number* values);
106 
107  //@}
108 
109  /** @name Solution Methods */
110  //@{
111  /** This method is called when the algorithm is complete so the TNLP can store/write the solution */
112  virtual void finalize_solution(SolverReturn status,
113  Index n, const Number* x, const Number* z_L, const Number* z_U,
114  Index m, const Number* g, const Number* lambda,
115  Number obj_value,
116  const IpoptData* ip_data,
117  IpoptCalculatedQuantities* ip_cq);
118  //@}
119 
120  /** Return information on each iteration */
121  virtual bool intermediate_callback(AlgorithmMode mode,
122  Index iter,
123  Number obj_value,
124  Number inf_pr,
125  Number inf_du,
126  Number mu,
127  Number d_norm,
128  Number regularization_size,
129  Number alpha_du,
130  Number alpha_pr,
131  Index ls_trials,
132  const IpoptData *ip_data,
133  IpoptCalculatedQuantities *ip_cq);
134 
135  //*************** start ADOL-C part ***********************************
136 
137  /** Method to generate the required tapes */
138  virtual void generate_tapes(Index n, Index m, Index& nnz_jac_g, Index& nnz_h_lag);
139 
140  //*************** end ADOL-C part ***********************************
141 
142  // ************** start FFSM part ************************************
143  void declareVariables(); ///< declare the variables, their domains and their bounds
144  void declareVariable(const string &name, const int & domain, const string &desc= "", const double & l_bound=0.0, const double & u_bound=UBOUND_MAX, const string & l_bound_var="", const string & u_bound_var=""); ///< Declare a single variable, its domain and its bounds
145  void declareConstrains(); ///< declare the constrains, their domain, their direction and their associated evaluation function
146  void cacheInitialPosition(); ///< cache the initial positions of the variables and the constrains
147  void calculateNumberVariablesConstrains(); ///< calculate the number of variables and constrains
148  void cachePositions(); ///< cache the exact position index (initial+f(r1,r2,p,r2To) for each variable and constrain
149  int getDomainElements(int domain); ///< return the number of elements of a domain
150 template<class T> vector < vector < vector < vector <int> > > > buildPositionVector(const T &v_or_c, int dType); ///< build the matrix of the positions for a given variable or contrain
151  int getVarInstances(const string& varName); ///< return the number of instances of a variable, given his domain type
152  ///< build the matrix of the positions for a given variable or contrain
153  void calculateSparsityPatternJ();
154  void calculateSparsityPatternH();
155 
156 
157  const Number& mymax(const Number& a, const Number& b);
158  const adouble& mymax(const adouble& a, const adouble& b);
159 
160  //template <class T> const T& mymax ( const T& a, const T& b );
161 
162  // ************** end FFSM part **************************************
163 
164 
165 protected:
166 
167  // convenient handles to equivalent ModelData functions..
168  const double gpd(const string &type_h, const int& regId_h, const string &prodId_h, const int& year=DATA_NOW, const string &freeDim_h="") const {return MTHREAD->MD->getProdData(type_h, regId_h, prodId_h, year, freeDim_h);};
169  const double gfd(const string &type_h, const int& regId_h, const string &forType_h, const string &diamClass_h, const int& year=DATA_NOW) const {return MTHREAD->MD->getForData(type_h, regId_h, forType_h, diamClass_h, year);};
170  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 {MTHREAD->MD->setProdData(value_h, type_h, regId_h, prodId_h, year, allowCreate, freeDim_h);};
171  void sfd(const double& value_h, const string &type_h, const int& regId_h, const string &forType_h, const string &diamClass_h, const int& year=DATA_NOW, const bool& allowCreate=false) const {MTHREAD->MD->setForData(value_h, type_h, regId_h, forType_h, diamClass_h, year, allowCreate);};
172  bool app(const string &prod_h, const string &forType_h, const string &dClass_h) const {return MTHREAD->MD->assessProdPossibility(prod_h, forType_h, dClass_h);};
173  const int gip(const string &varName) const; ///< Get the initial index position of a given variable in the concatenated array
174  const int gip(const int &cn) const; ///< Return the initial index position of a certain constrain
175  template<class T> const int gix_uncached(const T &v_or_c, int r1Ix, int r2Ix, int prIx, int r2IxTo=0); ///< Get the index in the concatenated array gived a certain var name (string) or constrain index (int), the reg lev1 index, the reg lev2 index and the prod. index
176  const int gix(const string &varName, const int& r1Ix, const int& r2Ix, const int& prIx, const int& r2IxTo=0) const; ///< Get the index in the concatenated array gived a certain var name, the reg lev1 index, the reg lev2 index and the prod. index
177  const int gix(const int &cn, const int& r1Ix, const int& r2Ix, const int& prIx, const int& r2IxTo=0) const; ///< Get the index in the concatenated array gived a certain constrain, the reg lev1 index, the reg lev2 index and the prod. index
178  const int gdt(const string &varName); ///< Get the domain type of a given variable
179  const int gdt(const int &cn); ///< Get the domain type of a given constrain
180  int getConstrainDirectionByIndex(int idx); ///< Return the direction of a given constrain
181  double getBoundByIndex(const int & bound_type, const int & idx); ///< Return the bound of a given variable (by index)
182  double getDetailedBoundByVarAndIndex(const endvar & var, const int & idx, const int & bType); ///< Return the bound of a given variable given the variable and the required index. Called by getBoundByIndex().
183  constrain* getConstrainByIndex(int idx);
184  void unpack(int ix_h, int domain, int initial, int &r1_h, int &r2_h, int&p_h, int&r2to_h, bool fullp=false); ///< Return the dimensions given a certain index, domain type and initial position
185  int getConNumber(constrain* con); ///< Return the position in the cons vector
186  //vector < vector <int> > createCombinationsVector(const int& nItems); ///< Return a vector containing any possible combination of nItems items (including any possible subset). The returned vector has in each slot the items present in that specific combination.
187  void copyInventoryResourses(); ///< Copy the inventoried resources in the in vector for better performances
188 
189  void tempDebug();
190 
191  //virtual void eval_obj (Index n, const T *x, T& obj_value);
192 
193  vector<string> priPr;
194  vector<string> secPr;
195  vector<string> allPr;
196  vector<string> othPr;
197  vector < vector <int> > l2r;
198  vector < vector <int> > priPrCombs; ///< A vector with all the possible combinations of primary products
199  vector < vector < vector <double> > > ins; ///< A copy of the inventoried resourses by region and primary product combination. It works also with dynamic loading of the region and the in, but it may be slower.
200  map <string, int> initPos; ///< A map that returns the initial index position in the concatenated array for each variable
201  map <int, string> initPos_rev; ///< A map with the name of the variable keyed by its initial position in the index
202  vector<int> cInitPos; ///< A vector that returns the initial index position in the concatenated array for each constrain
203  map <string, endvar> vars; ///< List of variables in the model and their domain: pr product, sec prod, all products or all products over each subregion pair (exports)
204  map <string, vector < vector < vector < vector <int> > > > > vpositions; ///< cached position in the concatenated vector for each variables. Dimensions are l1reg, l2reg, prod, (l2To region).
205  vector < vector < vector < vector < vector <int> > > > > cpositions; ///< cached position in the concatenated vector for each variables. Dimensions are contrain number, l1reg, l2reg, prod, (l2To region).
206  int nPriPr;
207  int nOthPr;
209  int nSecPr;
210  int nAllPr;
211  int nL2r;
212  int nVar;
213  int nCons;
222  double overharvestingAllowance; ///< Allows to harvest more than the resources available. Useful when resources got completelly exausted and the model refuses to solve.
223  void debugPrintParameters();
224  bool initOpt;
225  vector <constrain> cons;
226  vector <vector <Index> > nzjelements; ///< nzero elements for the jacobian matrix. nzelements[i][0] -> row (constrain), nzelements[i][1] -> column (variable)
227  vector <vector <Index> > nzhelements; ///< nzero elements for the hessian matrix
228 
229 
230  /**@name Methods to block default compiler methods.
231  * The compiler automatically generates the following three methods.
232  * Since the default compiler implementation is generally not what
233  * you want (for all but the most simple classes), we usually
234  * put the declarations of these methods in the private section
235  * and never implement them. This prevents the compiler from
236  * implementing an incorrect "default" behavior without us
237  * knowing. (See Scott Meyers book, "Effective C++")
238  *
239  */
240  //@{
241  // MyADOLC_NLP();
242  Opt(const Opt&);
243  Opt& operator=(const Opt&);
244  //@}
245 
246  //@{
247 
248  double *x_lam;
249 
250  //** variables for sparsity exploitation
251  unsigned int **HP_t; /* compressed block row storage */
252  unsigned int *rind_g; /* row indices */
253  unsigned int *cind_g; /* column indices */
254  double *jacval; /* values */
255  unsigned int *rind_L; /* row indices */
256  unsigned int *cind_L; /* column indices */
257  unsigned int *rind_L_total; /* row indices */
258  unsigned int *cind_L_total; /* column indices */
259  double *hessval; /* values */
260  int nnz_jac;
261  int nnz_L, nnz_L_total;
262  int options_g[4];
263  int options_L[4];
264 
265 
266  //@}
267 
268 };
269 
270 struct constrain{
271  constrain(){comment="";};
272  string name;
273  string comment;
274  int domain;
276 
277 };
278 
279 struct endvar{
280  string name;
281  int domain;
282  string desc; ///< Description of the variable
283  double l_bound; ///< A fixed numerical lower bound for all the domain
284  double u_bound; ///< A fixed numerical upper bound for all the domain
285  string l_bound_var; ///< A variable giving the lower bound. If present, the value defined in the variable overrides l_bound.
286  string u_bound_var; ///< A variable giving the upper bound. If present, the value defined in the variable overrides u_bound.
287 };
288 
289 
290 
291 #endif
292 
unsigned int * cind_L
Definition: Opt.h:256
The required data is for the current year.
Definition: BaseClass.h:73
int direction
Definition: Opt.h:275
double l_bound
A fixed numerical lower bound for all the domain.
Definition: Opt.h:283
unsigned int * rind_L
Definition: Opt.h:255
double u_bound
A fixed numerical upper bound for all the domain.
Definition: Opt.h:284
const double gfd(const string &type_h, const int &regId_h, const string &forType_h, const string &diamClass_h, const int &year=DATA_NOW) const
Definition: Opt.h:169
string desc
Description of the variable.
Definition: Opt.h:282
vector< string > priPr
Definition: Opt.h:193
vector< int > cInitPos
A vector that returns the initial index position in the concatenated array for each constrain...
Definition: Opt.h:202
unsigned int * rind_L_total
Definition: Opt.h:257
int nnz_L_total
Definition: Opt.h:261
string name
Definition: Opt.h:271
double overharvestingAllowance
Allows to harvest more than the resources available. Useful when resources got completelly exausted a...
Definition: Opt.h:222
int nVar
Definition: Opt.h:212
vector< vector< int > > priPrCombs
A vector with all the possible combinations of primary products.
Definition: Opt.h:198
int secondYear
Definition: Opt.h:219
int nL2r
Definition: Opt.h:211
vector< vector< int > > l2r
Definition: Opt.h:197
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: Opt.h:170
Thread manager. Responsable to manage the main thread and "speak" with the GUI.
Definition: ThreadManager.h:65
int domain
Definition: Opt.h:281
int nGreaterEqualZeroConstrains
Definition: Opt.h:216
unsigned int * rind_g
Definition: Opt.h:252
vector< string > secPr
Definition: Opt.h:194
vector< vector< vector< double > > > ins
A copy of the inventoried resourses by region and primary product combination. It works also with dyn...
Definition: Opt.h:199
int nAllPr
Definition: Opt.h:210
bool app(const string &prod_h, const string &forType_h, const string &dClass_h) const
Definition: Opt.h:172
int nOthPr
Definition: Opt.h:207
map< string, endvar > vars
List of variables in the model and their domain: pr product, sec prod, all products or all products o...
Definition: Opt.h:203
bool debugRunOnce
Definition: Opt.h:221
int firstYear
Definition: Opt.h:218
Definition: Opt.h:270
vector< vector< Index > > nzhelements
nzero elements for the hessian matrix
Definition: Opt.h:227
int nnz_jac
Definition: Opt.h:260
string u_bound_var
A variable giving the upper bound. If present, the value defined in the variable overrides u_bound...
Definition: Opt.h:286
unsigned int ** HP_t
Definition: Opt.h:251
int nLowerEqualZeroConstrains
Definition: Opt.h:215
double * hessval
Definition: Opt.h:259
map< int, string > initPos_rev
A map with the name of the variable keyed by its initial position in the index.
Definition: Opt.h:201
double * x_lam
Definition: Opt.h:248
int nPriPr
Definition: Opt.h:206
void sfd(const double &value_h, const string &type_h, const int &regId_h, const string &forType_h, const string &diamClass_h, const int &year=DATA_NOW, const bool &allowCreate=false) const
Definition: Opt.h:171
int nPriPrCombs
Definition: Opt.h:208
vector< constrain > cons
Definition: Opt.h:225
double * jacval
Definition: Opt.h:254
unsigned int * cind_g
Definition: Opt.h:253
Base class for the regmas application.
Definition: BaseClass.h:239
const 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: Opt.h:168
int nEqualityConstrains
Definition: Opt.h:214
string comment
Definition: Opt.h:273
vector< vector< Index > > nzjelements
nzero elements for the jacobian matrix. nzelements[i][0] -> row (constrain), nzelements[i][1] -> colu...
Definition: Opt.h:226
unsigned int * cind_L_total
Definition: Opt.h:258
Definition: Opt.h:52
Definition: Opt.h:279
string name
Definition: Opt.h:280
int worldCodeLev2
Definition: Opt.h:220
int previousYear
Definition: Opt.h:217
string l_bound_var
A variable giving the lower bound. If present, the value defined in the variable overrides l_bound...
Definition: Opt.h:285
#define UBOUND_MAX
Upper bound in optimisation 10^19.
Definition: BaseClass.h:174
int domain
Definition: Opt.h:274
int nCons
Definition: Opt.h:213
map< string, vector< vector< vector< vector< int > > > > > vpositions
cached position in the concatenated vector for each variables. Dimensions are l1reg, l2reg, prod, (l2To region).
Definition: Opt.h:204
int nSecPr
Definition: Opt.h:209
vector< string > othPr
Definition: Opt.h:196
This file is the header of BaseClass and it is included by ALL compiled code.
vector< string > allPr
Definition: Opt.h:195
bool initOpt
Definition: Opt.h:224
map< string, int > initPos
A map that returns the initial index position in the concatenated array for each variable.
Definition: Opt.h:200
vector< vector< vector< vector< vector< int > > > > > cpositions
cached position in the concatenated vector for each variables. Dimensions are contrain number...
Definition: Opt.h:205