Opt Class Reference

#include <Opt.h>

Inheritance diagram for Opt:

Collaboration diagram for Opt:

Public Member Functions
	Opt (ThreadManager *MTHREAD_h)
	Constructor. More...
	~Opt ()
virtual bool	intermediate_callback (AlgorithmMode mode, Index iter, Number obj_value, Number inf_pr, Number inf_du, Number mu, Number d_norm, Number regularization_size, Number alpha_du, Number alpha_pr, Index ls_trials, const IpoptData *ip_data, IpoptCalculatedQuantities *ip_cq)
virtual void	generate_tapes (Index n, Index m, Index &nnz_jac_g, Index &nnz_h_lag)
void	declareVariables ()
	declare the variables, their domains and their bounds More...
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. More...
void	declareConstrains ()
	declare the constrains, their domain, their direction and their associated evaluation function More...
void	cacheInitialPosition ()
	cache the initial positions of the variables and the constrains More...
void	calculateNumberVariablesConstrains ()
	calculate the number of variables and constrains More...
void	cachePositions ()
	cache the exact position index (initial+f(r1,r2,p,r2To) for each variable and constrain More...
int	getDomainElements (int domain)
	return the number of elements of a domain More...
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 More...
int	getVarInstances (const string &varName)
	build the matrix of the positions for a given variable or contrain More...
void	calculateSparsityPatternJ ()
void	calculateSparsityPatternH ()
const Number &	mymax (const Number &a, const Number &b)
const adouble &	mymax (const adouble &a, const adouble &b)
Overloaded from TNLP
virtual bool	get_nlp_info (Index &n, Index &m, Index &nnz_jac_g, Index &nnz_h_lag, IndexStyleEnum &index_style)
virtual bool	get_bounds_info (Index n, Number *x_l, Number *x_u, Index m, Number *g_l, Number *g_u)
virtual bool	get_starting_point (Index n, bool init_x, Number *x, bool init_z, Number *z_L, Number *z_U, Index m, bool init_lambda, Number *lambda)
template<class T >
bool	eval_obj (Index n, const T *x, T &obj_value)
template<class T >
bool	eval_constraints (Index n, const T *x, Index m, T *g)
virtual bool	eval_f (Index n, const Number *x, bool new_x, Number &obj_value)
virtual bool	eval_grad_f (Index n, const Number *x, bool new_x, Number *grad_f)
virtual bool	eval_g (Index n, const Number *x, bool new_x, Index m, Number *g)
virtual bool	eval_jac_g (Index n, const Number *x, bool new_x, Index m, Index nele_jac, Index *iRow, Index *jCol, Number *values)
virtual bool	eval_h (Index n, const Number *x, bool new_x, Number obj_factor, Index m, const Number *lambda, bool new_lambda, Index nele_hess, Index *iRow, Index *jCol, Number *values)
Solution Methods
virtual void	finalize_solution (SolverReturn status, Index n, const Number *x, const Number *z_L, const Number *z_U, Index m, const Number *g, const Number *lambda, Number obj_value, const IpoptData *ip_data, IpoptCalculatedQuantities *ip_cq)
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

Protected Member Functions
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
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
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 &diamClass_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
const int	gip (const string &varName) const
	Get the initial index position of a given variable in the concatenated array. More...
const int	gip (const int &cn) const
	Return the initial index position of a certain constrain. More...
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. More...
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. More...
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. More...
const int	gdt (const string &varName)
	Get the domain type of a given variable. More...
const int	gdt (const int &cn)
	Get the domain type of a given constrain. More...
int	getConstrainDirectionByIndex (int idx)
	Return the direction of a given constrain. More...
double	getBoundByIndex (const int &bound_type, const int &idx)
	Return the bound of a given variable (by index) More...
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(). More...
constrain *	getConstrainByIndex (int idx)
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. More...
int	getConNumber (constrain *con)
	Return the position in the cons vector. More...
void	copyInventoryResourses ()
	Copy the inventoried resources in the in vector for better performances. More...
void	tempDebug ()
void	debugPrintParameters ()

Protected Attributes
vector< string >	priPr
vector< string >	secPr
vector< string >	allPr
vector< string >	othPr
vector< vector< int > >	l2r
vector< vector< int > >	priPrCombs
	A vector with all the possible combinations of primary products. More...
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. More...
map< string, int >	initPos
	A map that returns the initial index position in the concatenated array for each variable. More...
map< int, string >	initPos_rev
	A map with the name of the variable keyed by its initial position in the index. More...
vector< int >	cInitPos
	A vector that returns the initial index position in the concatenated array for each constrain. More...
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) More...
map< string, vector< vector< vector< vector< int > > > > >	vpositions
	cached position in the concatenated vector for each variables. Dimensions are l1reg, l2reg, prod, (l2To region). More...
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). More...
int	nPriPr
int	nOthPr
int	nPriPrCombs
int	nSecPr
int	nAllPr
int	nL2r
int	nVar
int	nCons
int	nEqualityConstrains
int	nLowerEqualZeroConstrains
int	nGreaterEqualZeroConstrains
int	previousYear
int	firstYear
int	secondYear
int	worldCodeLev2
bool	debugRunOnce
double	overharvestingAllowance
	Allows to harvest more than the resources available. Useful when resources got completelly exausted and the model refuses to solve. More...
bool	initOpt
vector< constrain >	cons
vector< vector< Index > >	nzjelements
	nzero elements for the jacobian matrix. nzelements[i][0] -> row (constrain), nzelements[i][1] -> column (variable) More...
vector< vector< Index > >	nzhelements
	nzero elements for the hessian matrix More...
Protected Attributes inherited from BaseClass
ThreadManager *	MTHREAD
	Pointer to the Thread manager. More...

Methods to block default compiler methods.
double *	x_lam
unsigned int **	HP_t
unsigned int *	rind_g
unsigned int *	cind_g
double *	jacval
unsigned int *	rind_L
unsigned int *	cind_L
unsigned int *	rind_L_total
unsigned int *	cind_L_total
double *	hessval
int	nnz_jac
int	nnz_L
int	nnz_L_total
int	options_g [4]
int	options_L [4]
	Opt (const Opt &)
Opt &	operator= (const Opt &)

Detailed Description

Definition at line 52 of file Opt.h.

Constructor & Destructor Documentation

Opt

(

ThreadManager *

MTHREAD_h

)

Constructor.

Definition at line 537 of file Opt.cpp.

                                {
  MTHREAD = MTHREAD_h;
  nVar    = 0;
  nCons   = 0;
  debugRunOnce = false;
  initOpt = true;
}

~Opt

(

)

Definition at line 545 of file Opt.cpp.

         {

}

Opt ( const Opt &  )

protected

Member Function Documentation

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

inlineprotected

Definition at line 172 of file Opt.h.

{return MTHREAD->MD->assessProdPossibility(prod_h, forType_h, dClass_h);};

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

Definition at line 1401 of file Opt.cpp.

                                                  {
  int pVectorSize;

  switch (dType){
    case DOM_PRI_PR:
      pVectorSize= priPr.size();
      break;
    case DOM_SEC_PR:
      pVectorSize= secPr.size();
      break;
    case DOM_ALL_PR:
      pVectorSize= allPr.size();
      break;
    case DOM_R2_PRI_PR:
      pVectorSize= priPr.size();
      break;
    case DOM_R2_SEC_PR:
      pVectorSize= secPr.size();
      break;
    case DOM_R2_ALL_PR:
      pVectorSize= allPr.size();
      break;
    case DOM_SCALAR:
      pVectorSize= allPr.size(); // it will simply fill the matrix all with the same value (the ip)
      break;
    case DOM_PRI_PR_ALLCOMBS:
      pVectorSize= priPrCombs.size();
      break;
    default:
      msgOut(MSG_CRITICAL_ERROR,"Try to build the position of a variable (or contrain) of unknow type.");
  }


  vector < vector < vector < vector <int> > > > positionsToAdd;
  for(uint r1=0;r1<l2r.size();r1++){
    vector < vector < vector <int> > > dim1;
    for(uint r2=0;r2<l2r[r1].size();r2++){
      vector < vector <int> > dim2;
      for(uint p=0;p<pVectorSize;p++){
        vector <int> dim3;
          for(uint r2To=0;r2To<l2r[r1].size();r2To++){
            dim3.push_back(gix_uncached(v_or_c,r1,r2,p,r2To));
          }
        dim2.push_back(dim3);
      }
      dim1.push_back(dim2);
    }
    positionsToAdd.push_back(dim1);
  }
  return positionsToAdd;
}

void cacheInitialPosition

(

)

cache the initial positions of the variables and the constrains

Definition at line 1370 of file Opt.cpp.

                         {
  int vInitialPosition = 0;
  int cInitialPosition = 0;
  VarMap::iterator viter;
  for (viter = vars.begin(); viter != vars.end(); ++viter) {
    initPos.insert(pair<string, int>(viter->first, vInitialPosition));
    initPos_rev.insert(pair<int, string>(vInitialPosition,viter->first));
    vInitialPosition += getDomainElements(viter->second.domain);
  }
  for (uint i=0;i<cons.size();i++){
    cInitPos.push_back(cInitialPosition);
    cInitialPosition += getDomainElements(cons[i].domain);
  }
}

void cachePositions

(

)

cache the exact position index (initial+f(r1,r2,p,r2To) for each variable and constrain

Definition at line 1386 of file Opt.cpp.

                   {

  // variables..
  VarMap::iterator viter;
  for (viter = vars.begin(); viter != vars.end(); ++viter) {
    vpositions.insert(pair<string, vector < vector < vector < vector <int> > > > >(viter->first,  buildPositionVector(viter->first, viter->second.domain)));
  }
  // constrains..
  for (uint i=0; i<cons.size();i++){
    cpositions.push_back(buildPositionVector(i, cons[i].domain));
  }

}

void calculateNumberVariablesConstrains

(

)

calculate the number of variables and constrains

Definition at line 1455 of file Opt.cpp.

                                       {
    // calculating the number of variables and the initial positions in the concatenated array..
    nVar = 0;
    VarMap::iterator viter;
    for (viter = vars.begin(); viter != vars.end(); ++viter) {
        nVar += getDomainElements(viter->second.domain);
    }

    // calculating the number of constrains..
    nCons = 0;
    nEqualityConstrains = 0;
    nLowerEqualZeroConstrains = 0;
    nGreaterEqualZeroConstrains = 0;
    for(uint i=0;i<cons.size();i++){
      nCons += getDomainElements(cons[i].domain);
      if(cons[i].direction == CONSTR_EQ){
        nEqualityConstrains += getDomainElements(cons[i].domain);
        continue;
      } else if (cons[i].direction == CONSTR_LE0) {
        nLowerEqualZeroConstrains += getDomainElements(cons[i].domain);
        continue;
      } else if (cons[i].direction == CONSTR_GE0) {
        nGreaterEqualZeroConstrains += getDomainElements(cons[i].domain);
        continue;
      } else {
        msgOut(MSG_CRITICAL_ERROR, "Asking for a constrain with unknown direction ("+i2s(cons[i].direction)+")");
      }
    }

    msgOut(MSG_INFO,"The model will work with "+i2s(nVar)+" variables and "+i2s(nCons)+" constrains ("+i2s(nEqualityConstrains)+" equalities, "+i2s(nLowerEqualZeroConstrains)+" lower than 0 and "+i2s(nGreaterEqualZeroConstrains)+" greater than 0)");
}

void calculateSparsityPatternH

(

)

Definition at line 1708 of file Opt.cpp.

                              {

  unsigned int  **hesspat=NULL; // compressed row storage
  int           options_h=0; // options for the hessian patterns
  double        *x;
  int retv_h = -1; // return value

  hesspat = new unsigned int* [(nVar+nCons+1)];
  x = new double[(nVar+nCons+1)];

  retv_h = hess_pat(tag_L,nVar+nCons+1,  x, hesspat, options_h);

  for (int i=0;i<(nVar);i++) {
    for (int j=1;j<=hesspat[i][0];j++){
      if(hesspat[i][j]<=i){
        vector <int> nzhelement;
        nzhelement.push_back(i);
        nzhelement.push_back(hesspat[i][j]);
        nzhelements.push_back(nzhelement);
      }
    }
  }
}

void calculateSparsityPatternJ

(

)

Definition at line 1680 of file Opt.cpp.

                              {

  unsigned int  **jacpat=NULL; // compressed row storage
  int           options_j[3]; // options for the jacobian patterns
  double        *x;
  int retv_j = -1; // return value

  options_j[0] = 0; // index domain propagation
  options_j[1] = 0; // automatic mode choice (ignored here)
  options_j[2] = 0; // safe
  jacpat = new unsigned int* [nCons];
  x = new double[nVar];

  nzjelements.clear();

  retv_j = jac_pat(tag_g, nCons, nVar,  x, jacpat, options_j);

  for (int i=0;i<nCons;i++) {
    for (int j=1;j<=jacpat[i][0];j++){
      vector <int> nzjelement;
      nzjelement.push_back(i);
      nzjelement.push_back(jacpat[i][j]);
      nzjelements.push_back(nzjelement);
    }
  }
}

void copyInventoryResourses ( )

protected

Copy the inventoried resources in the in vector for better performances.

Opt::createCombinationsVector Return a vector containing any possible combination of nItems items (including all subsets).

For example with nItems = 3: 0: []; 1: [0]; 2: [1]; 3: [0,1]; 4: [2]; 5: [0,2]; 6: [1,2]; 7: [0,1,2]

Parameters

nItems

number of items to create p

Returns

A vector with in each slot the items present in that specific combination subset.

Definition at line 1845 of file Opt.cpp.

                           {
  // This function is not really needed, as actually the solver works also picking the region and the in dynamically
  // Caching the inventories in a vector should however be faster.
  // We now need it, as the vector inResByAnyCombination() account for the union between the inv set of the various pp. Also it now include the total mortality (alive plus death, if modelled)
  vector < vector < vector <double> > >  in_temp;
  for (uint r1=0;r1<l2r.size();r1++){
     vector < vector <double> > dim1;
    for (uint r2=0;r2<l2r[r1].size();r2++){
      vector <double> dim2;
      ModelRegion* REG = MTHREAD->MD->getRegion(l2r[r1][r2]);
      for (uint p=0;p<priPrCombs.size();p++){
        double this_in = REG->inResByAnyCombination[p];
        dim2.push_back(this_in);
      }
      dim1.push_back(dim2);
    }
    in_temp.push_back(dim1);
  }
  ins = in_temp;
}

void debugPrintParameters ( )

protected

void declareConstrains

(

)

declare the constrains, their domain, their direction and their associated evaluation function

Declare the constrains and their properties. For the domain type

See also

BaseClass

Definition at line 84 of file Opt.cpp.

                       {
  // domain of constrains variables
  // for domain
  constrain mkeq2;
  mkeq2.name="mkeq2";
  mkeq2.comment="[h1] Conservation of matters of transformed products";
  mkeq2.domain=DOM_SEC_PR;
  mkeq2.direction = CONSTR_EQ;
  //mkeq2.evaluate = Opt::mkteq2f;

  constrain mkeq3;
  mkeq3.name="mkeq3";
  mkeq3.comment="[h2] Conservation of matters of raw products";
  mkeq3.domain=DOM_PRI_PR;
  mkeq3.direction = CONSTR_EQ;
  //mkeq3.evaluate = Opt::mkteq3f;

  constrain mkeq4;
  mkeq4.name="mkeq4";
  mkeq4.comment="[eq 13] Leontief transformation function";
  mkeq4.domain=DOM_PRI_PR;
  mkeq4.direction = CONSTR_EQ;

  constrain mkeq5;
  mkeq5.name="mkeq5";
  mkeq5.comment="[eq 21] Raw product supply function";
  mkeq5.domain=DOM_PRI_PR;
  mkeq5.direction = CONSTR_EQ;

  constrain mkeq6;
  mkeq6.name="mkeq6";
  mkeq6.comment="[eq 20] Trasformed products demand function";
  mkeq6.domain=DOM_SEC_PR;
  mkeq6.direction = CONSTR_EQ;

  constrain mkeq7;
  mkeq7.name="mkeq7";
  mkeq7.comment="[h7 and h3] Transformed products import function";
  mkeq7.domain=DOM_SEC_PR;
  mkeq7.direction = CONSTR_EQ;

  constrain mkeq8;
  mkeq8.name="mkeq8";
  mkeq8.comment="[h8 and h4] Raw products export function";
  mkeq8.domain=DOM_PRI_PR;
  mkeq8.direction = CONSTR_EQ;

  constrain mkeq13;
  mkeq13.name="mkeq13";
  mkeq13.comment="[h9] Calculation of the composite price of transformed products (PPC_Dp)";
  mkeq13.domain=DOM_SEC_PR;
  mkeq13.direction = CONSTR_EQ;

  constrain mkeq14;
  mkeq14.name="mkeq14";
  mkeq14.comment="[h10] Calculation of the composite price of raw products (PPC_Sw)";
  mkeq14.domain=DOM_PRI_PR;
  mkeq14.direction = CONSTR_EQ;

  constrain mkeq17;
  mkeq17.name="mkeq17";
  mkeq17.comment="[h16] Constrain of the transformaton supply (lower than the regional maximal production capacity)";
  mkeq17.domain=DOM_SEC_PR;
  mkeq17.direction = CONSTR_LE0;


  constrain mkeq23;
  mkeq23.name="mkeq23";
  mkeq23.comment="[h3] Composit demand eq. (Dp)";
  mkeq23.domain=DOM_SEC_PR;
  mkeq23.direction = CONSTR_EQ;

  constrain mkeq24;
  mkeq24.name="mkeq24";
  mkeq24.comment="[h4] Composite supply eq. (Sw)";
  mkeq24.domain=DOM_PRI_PR;
  mkeq24.direction = CONSTR_EQ;

  constrain mkeq26;
  mkeq26.name="mkeq26";
  mkeq26.comment="[eq ] Verification of the null transport agents supply";
  mkeq26.domain=DOM_R2_ALL_PR;
  mkeq26.direction = CONSTR_LE0;

  constrain mkeq25;
  mkeq25.name="mkeq25";
  mkeq25.comment="Verification of the null trasformers supply (price of raw product + trasf product > trasf product)";
  mkeq25.domain=DOM_SEC_PR;
  mkeq25.direction = CONSTR_GE0;

  constrain mkeq18;
  mkeq18.name="mkeq18";
  mkeq18.comment="Constrain on raw material supply (lower than inventory)";
  mkeq18.domain=DOM_PRI_PR;
  mkeq18.direction = CONSTR_LE0;

  constrain resbounds;
  resbounds.name="resbounds";
  resbounds.comment="Constrain on raw material supply (lower than inventory, for each possible combination of primary products)";
  resbounds.domain=DOM_PRI_PR_ALLCOMBS;
  resbounds.direction = CONSTR_LE0;



  //constrain steq;
  //steq.name="steq";
  //steq.comment="computation of total supply";
  //steq.domain=DOM_PRI_PR;
  //steq.direction = CONSTR_EQ;

  cons.push_back(mkeq2);
  cons.push_back(mkeq6);
  cons.push_back(mkeq7);
  cons.push_back(mkeq13);
  cons.push_back(mkeq23);
  cons.push_back(mkeq3);
  cons.push_back(mkeq4);
  cons.push_back(mkeq5);
  cons.push_back(mkeq8);
  cons.push_back(mkeq14);
  cons.push_back(mkeq24);
  cons.push_back(mkeq17);
  cons.push_back(mkeq26);
  cons.push_back(mkeq25);
  //cons.push_back(mkeq18);
  cons.push_back(resbounds);
  //cons.push_back(steq);
;



}

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.

Opt::declareVariable Define a single variable together with its domain and optionally its lower and upper bound (default 0.0, +inf)

Parameters

name	var name
domain	domain of the variable
l_bound	lower bound (fixed)
u_bound	upper bound (fixed)
l_bound_var	variable name defining lower bound
u_bound_var	variable name defining upper bound

Definition at line 1791 of file Opt.cpp.

                                                                                                                                                                                      {
   endvar end_var;
   end_var.name = name;
   end_var.domain = domain;
   end_var.l_bound = l_bound;
   end_var.u_bound = u_bound;
   end_var.l_bound_var = l_bound_var;
   end_var.u_bound_var = u_bound_var;
   end_var.desc= desc;
   vars.insert(std::pair<std::string, endvar >(name, end_var));
}

void declareVariables

(

)

declare the variables, their domains and their bounds

Definition at line 59 of file Opt.cpp.

                     {
    // filling the list of variables and their domain and optionally their bonds
    // if you add variables in the model that enter optimisation you'll have to add them here
    // the underlying map goes automatically in alphabetical order
    // original order: pc,pl,dc,dl,da,sc,sl,sa,exp
    // 20140328: if these vars have a lower bound > 0 the model doesn't solve when volumes in a region go to zero !!!

    // syntax: declareVariable("name", domainType, lbound[default=0], ubound[default= +inf], variable defining lower bounds[default=""], variable defining upper bound[default=""])

    // all variables have upper or equal than zero bound:
    declareVariable("da", DOM_SEC_PR,    "Demand from abroad (imports)");
    declareVariable("dc", DOM_SEC_PR,    "Demand, composite");
    declareVariable("dl", DOM_ALL_PR,    "Demand from local");
    declareVariable("pc", DOM_ALL_PR,    "Price, composite");
    declareVariable("pl", DOM_ALL_PR,    "Price, local") ;
    declareVariable("rt", DOM_R2_ALL_PR, "Regional trade"); //it was exp in gams
    declareVariable("sa", DOM_PRI_PR,    "Supply to abroad (exports)");
    declareVariable("sc", DOM_PRI_PR,    "Supply, composite");
    declareVariable("sl", DOM_ALL_PR,    "Supply to locals");
    //declareVariable("st", DOM_PRI_PR,    "Supply, total", 0.0,UBOUND_MAX,"","in");
}

bool eval_constraints	(	Index	n,
		const T *	x,
		Index	m,
		T *	g
	)

Template to compute contraints

Template function to implement (define) the previously declared constains. To the initial macro loop it must be passed the product vector over where to loop (priPr, secPr or allPr) and the order of the constrain has it has been added to the const vector. It could be possible to change this in a map and uses name, but then we would loose control on the constrains order, and we saw that it matters for finding the equilibrium.

Definition at line 311 of file Opt.cpp.

                                                       {

  double a_pr, a, sigma, ff, pol_mktDirInt_s, pol_mktDirInt_d, pol_mktDirInt_d_pSubstituted, pol_mktDirInt_d_1, pol_mktDirInt_d_1_pSubstituted, gg, q1, p1v, t1, r1v, psi, eta, pworld, ct, k, dispor, mv, in, in_1, supCorr, es_d, pc_1, pc_pSubstituted, pc_1_pSubstituted, additionalDemand;
  Index cix = 0;
  Index debug = 0;

  // mkteq2(i,p_tr)..  RVAR('dl',i,p_tr)+sum(j,EXP(i,j,p_tr))  =e=  RVAR('sl',i,p_tr)+ sum(b,EXP(b,i,p_tr)); // h1
  CONSTRAIN_START_LOOP(secPr, 0) // attenction! you have to give the same order number as you inserted in the cons vector
    //g[cix] = x[gix("dl",r1,r2,psec)]-x[gix("sl",r1,r2,psec)]+x[gix("da",r1,r2,p)];
    g[cix] = x[gix("dl",r1,r2,psec)]-x[gix("sl",r1,r2,psec)];
    for (uint r2To=0;r2To<l2r[r1].size();r2To++){
      g[cix] += x[gix("rt",r1,r2,psec,r2To)]-x[gix("rt",r1,r2To,psec,r2)];
    }
  CONSTRAIN_END_LOOP


  // mkteq6(i,p_tr)..  RVAR('dc',i,p_tr)  =e=  GG(i,p_tr)*(RVAR('pc',i,p_tr)**sigma(p_tr)); // eq. 20 20160216: added sustitution elasticity in the demand
  // DEMAND EQUATION of transformed products
  CONSTRAIN_START_LOOP(secPr,1)
    gg         = gpd("gg",l2r[r1][r2],secPr[p]);
    sigma      = gpd("sigma",l2r[r1][r2],secPr[p]);
    pc_1       = gpd("pc",l2r[r1][r2],secPr[p],previousYear);
    pol_mktDirInt_d   = gpd("pol_mktDirInt_d",l2r[r1][r2],secPr[p]);              // market policies (subs or taxes) this year
    pol_mktDirInt_d_1 = gpd("pol_mktDirInt_d",l2r[r1][r2],secPr[p],previousYear); // market policies previous year
    additionalDemand = gpd("additionalDemand",l2r[r1][r2],secPr[p]);
    g[cix]     = - gg*pow(x[gix("pc",r1,r2,psec)]*pol_mktDirInt_d,sigma);  // note that putting 0.0 and 0.000 IS NOT the same thing. With 0.000 I have a seg fault in ADOL-C !!
    vector<string> debug = secPr;
    // Substitution part within forest producs (whose price is endogenous to the model)..
    for (uint p2=0;p2<secPr.size();p2++){
      es_d     = gpd("es_d",l2r[r1][r2],secPr[p],DATA_NOW,secPr[p2]);
      if(es_d){
        pc_1_pSubstituted    = gpd("pc",l2r[r1][r2],secPr[p2],previousYear);
        pol_mktDirInt_d_pSubstituted   = gpd("pol_mktDirInt_d",l2r[r1][r2],secPr[p2]);               // market policies this year for the substitute product
        pol_mktDirInt_d_1_pSubstituted = gpd("pol_mktDirInt_d",l2r[r1][r2],secPr[p2],previousYear);  // market policies last year for the substitute product

        g[cix] *= pow(
                 (
                   ((x[gix("pc",r1,r2,psec)]*pol_mktDirInt_d) / (x[gix("pc",r1,r2,priPr.size()+p2)]*pol_mktDirInt_d_pSubstituted))
                      /
                   ((pc_1*pol_mktDirInt_d_1) / (pc_1_pSubstituted*pol_mktDirInt_d_1_pSubstituted))
                 ), es_d
               );
      }
    }
    // Substitution part with other products with exogenous prices (e.g. fossilFuels)..
    for (uint p3=0;p3<othPr.size();p3++){
      es_d     = gpd("es_d",l2r[r1][r2],secPr[p],DATA_NOW,othPr[p3]);
      if(es_d){
        pc_pSubstituted    = gpd("pl",l2r[r1][r2],othPr[p3],DATA_NOW);
        pc_1_pSubstituted  = gpd("pl",l2r[r1][r2],othPr[p3],previousYear);

        g[cix] *= pow(
                 (
                   (x[gix("pc",r1,r2,psec)]*pol_mktDirInt_d / pc_pSubstituted)
                      /
                   (pc_1*pol_mktDirInt_d_1 / pc_1_pSubstituted)
                 ), es_d
               );
      }
    }
    //g[cix] = x[gix("dc",r1,r2,p)]-gg*pow(x[gix("pc",r1,r2,psec)],sigma); // original without substitution elasticity
    g[cix] += x[gix("dc",r1,r2,p)]-additionalDemand;
  CONSTRAIN_END_LOOP


  // mkteq7(i,p_tr)..  RVAR('da',i,p_tr)/RVAR('dl',i,p_tr)  =e=  ((q1(i,p_tr)*RVAR('pl',i,p_tr))/(p1(i,p_tr)*PT_t(p_tr)))**psi(i,p_tr); // h7 and h3 ?
  CONSTRAIN_START_LOOP(secPr,2)
    q1 = gpd("q1_polCorr",l2r[r1][r2],secPr[p]);
    p1v = 1-q1;
    psi = gpd("psi",l2r[r1][r2],secPr[p]);
    pworld = gpd("pl", worldCodeLev2,secPr[p]);
    g[cix] = x[gix("da",r1,r2,p)]/x[gix("dl",r1,r2,psec)] - pow((q1*x[gix("pl",r1,r2,psec)])/(p1v*pworld),psi);
  CONSTRAIN_END_LOOP

  // mkteq13(i,p_tr).. RVAR('pc',i,p_tr)*RVAR('dc',i,p_tr)  =e=  RVAR('dl',i,p_tr)*RVAR('pl',i,p_tr)+RVAR('da',i,p_tr)*PT_t(p_tr);  // h9
  CONSTRAIN_START_LOOP(secPr,3)
    pworld = gpd("pl", worldCodeLev2,secPr[p]);
    g[cix] = x[gix("pc",r1,r2,psec)]*x[gix("dc",r1,r2,p)]-x[gix("dl",r1,r2,psec)]*x[gix("pl",r1,r2,psec)]-x[gix("da",r1,r2,p)]*pworld;
  CONSTRAIN_END_LOOP

  // mkteq23(i,p_tr).. RVAR('dc',i,p_tr)  =e=  (q1(i,p_tr)*(RVAR('da',i,p_tr)**((psi(i,p_tr)-1)/psi(i,p_tr)))+ p1(i,p_tr)*(RVAR('dl',i,p_tr)**((psi(i,p_tr)-1)/psi(i,p_tr))))**(psi(i,p_tr)/(psi(i,p_tr)-1)); // h3
  CONSTRAIN_START_LOOP(secPr,4)
    q1 = gpd("q1_polCorr",l2r[r1][r2],secPr[p]);
    psi = gpd("psi",l2r[r1][r2],secPr[p]);
    p1v = 1-q1;
    g[cix] = x[gix("dc",r1,r2,p)] -
      pow(
            q1 * pow(x[gix("da",r1,r2,p)],(psi-1)/psi)
          + p1v * pow(x[gix("dl",r1,r2,psec)],(psi-1)/psi),
          psi/(psi-1)
      );
  CONSTRAIN_END_LOOP

  // mkteq3(i,p_pr)..  RVAR('dl',i,p_pr)+sum(j,EXP(i,j,p_pr))   =e=  RVAR('sl',i,p_pr)+ sum(b,EXP(b,i,p_pr))+sum(p_pr2, pres(p_pr2,p_pr)* RVAR('sl',i,p_pr2)); // h2
  CONSTRAIN_START_LOOP(priPr,5)
    //g[cix] = x[gix("dl",r1,r2,p)]-x[gix("sl",r1,r2,p)]-x[gix("sa",r1,r2,p)];
    g[cix] = x[gix("dl",r1,r2,p)]-x[gix("sl",r1,r2,p)];
    for (uint r2To=0;r2To<l2r[r1].size();r2To++){
      g[cix] += x[gix("rt",r1,r2,p,r2To)]-x[gix("rt",r1,r2To,p,r2)];
    }
    for (uint p2=0;p2<priPr.size();p2++){
      a_pr = gpd("a_pr",l2r[r1][r2],priPr[p2],DATA_NOW,priPr[p]);
      g[cix] -= a_pr*x[gix("sl",r1,r2,p2)];
    }
  CONSTRAIN_END_LOOP

  //mkteq4(i,p_pr)..  RVAR('dl',i,p_pr)  =e=  sum(p_tr, a(p_pr,p_tr)*(RVAR('sl',i,p_tr))); // eq. 13
  CONSTRAIN_START_LOOP(priPr,6)
    g[cix] = x[gix("dl",r1,r2,p)];
    for (uint p2=0;p2<secPr.size();p2++){
      a = gpd("a",l2r[r1][r2],priPr[p],DATA_NOW,secPr[p2]);
      g[cix] -= a*x[gix("sl",r1,r2,p2+nPriPr)];
    }
  CONSTRAIN_END_LOOP


  // mkteq5(i,p_pr)..  RVAR('sc',i,p_pr)  =e=  FF(i,p_pr)*(RVAR('pc',i,p_pr)**sigma(p_pr)); // eq. 21
  // SUPPLY EQUATION OF PRIMARY PRODUCTS


  CONSTRAIN_START_LOOP(priPr,7)
    ff = gpd("ff",l2r[r1][r2],priPr[p]);
    pol_mktDirInt_s = gpd("pol_mktDirInt_s",l2r[r1][r2],priPr[p]);
    sigma = gpd("sigma",l2r[r1][r2],priPr[p]);

    // Added for carbon project ----------------------
    double carbonPrice = gpd("carbonPrice",l2r[r1][r2],"",DATA_NOW); //using dummy region until Anto solves issue // Solved and generalised ;-)
    double intRate     = MTHREAD->MD->getDoubleSetting("ir");
    double Pct         = carbonPrice*intRate;
    double omega       = gpd("co2content_products", l2r[r1][r2], priPr[p])/1000;  // Generalised
    // ------------------------------------------------

    g[cix] = x[gix("sc",r1,r2,p)]-ff*pow((x[gix("pc",r1,r2,p)]-omega*Pct)*pol_mktDirInt_s,sigma);
    // -----------------------------------------------

    //g[cix] = x[gix("sc",r1,r2,p)]-ff*pow(x[gix("pc",r1,r2,p)]*pol_mktDirInt_s,sigma);
  CONSTRAIN_END_LOOP


 // mkteq8(i,p_pr)..  RVAR('sa',i,p_pr)/RVAR('sl',i,p_pr)  =e=  ((t1(i,p_pr)*RVAR('pl',i,p_pr))/(r1(i,p_pr)*PT_t(p_pr)))**eta(i,p_pr); // h8 and h4 ?
  CONSTRAIN_START_LOOP(priPr,8)
    t1 = gpd("t1_polCorr",l2r[r1][r2],priPr[p]);
    r1v = 1-t1;
    eta = gpd("eta",l2r[r1][r2],priPr[p]);
    pworld = gpd("pl", worldCodeLev2,priPr[p]);
    g[cix] = x[gix("sa",r1,r2,p)]/x[gix("sl",r1,r2,p)] - pow((t1*x[gix("pl",r1,r2,p)])/(r1v*pworld),eta);
  CONSTRAIN_END_LOOP

  // mkteq14(i,p_pr).. RVAR('pc',i,p_pr)*RVAR('sc',i,p_pr)  =e=  RVAR('sl',i,p_pr)*RVAR('pl',i,p_pr)+RVAR('sa',i,p_pr)*PT_t(p_pr);  // h10
  CONSTRAIN_START_LOOP(priPr,9)
    pworld = gpd("pl", worldCodeLev2,priPr[p]);
    g[cix] = x[gix("pc",r1,r2,p)]*x[gix("sc",r1,r2,p)]-x[gix("sl",r1,r2,p)]*x[gix("pl",r1,r2,p)]-x[gix("sa",r1,r2,p)]*pworld;
  CONSTRAIN_END_LOOP

  //mkteq24(i,p_pr).. RVAR('sc',i,p_pr)  =e=  (t1(i,p_pr)*(RVAR('sa',i,p_pr)**((eta(i,p_pr)-1)/eta(i,p_pr)))+ r1(i,p_pr)*(RVAR('sl',i,p_pr)**((eta(i,p_pr)-1)/eta(i,p_pr))))**(eta(i,p_pr)/(eta(i,p_pr)-1)); // h4
  CONSTRAIN_START_LOOP(priPr,10)
    t1 = gpd("t1_polCorr",l2r[r1][r2],priPr[p]);
    r1v = 1-t1;
    eta = gpd("eta",l2r[r1][r2],priPr[p]);
    g[cix] = x[gix("sc",r1,r2,p)] -
          pow(
             t1 * pow(x[gix("sa",r1,r2,p)],(eta-1)/eta)
          + r1v * pow(x[gix("sl",r1,r2,p)],(eta-1)/eta),
          eta/(eta-1)
        );
  CONSTRAIN_END_LOOP

  // mkteq17(i,p_tr).. RVAR('sl',i,p_tr)  =l=  Kt(i,p_tr);     // h16 in the presentation paper
  CONSTRAIN_START_LOOP(secPr,11)
    k = gpd("k",l2r[r1][r2],secPr[p]);
    g[cix] = x[gix("sl",r1,r2,p+nPriPr)]-k;
  CONSTRAIN_END_LOOP

  // mkeq26(i,prd,j).. RVAR('pl',j,prd)-RVAR('pl',i,prd)-CT(i,j,prd) =l= 0;
  CONSTRAIN_START_LOOP(allPr,12)
    for (uint r2To=0;r2To<l2r[r1].size();r2To++){
      cix = gix(12, r1, r2, p,r2To); // attenction we must redefine it, as we are now in a r2to loop
      ct = (gpd("ct",l2r[r1][r2],allPr[p],DATA_NOW,i2s(l2r[r1][r2To])) *  gpd("pol_tcSub",l2r[r1][r2],allPr[p])) ;
      g[cix] = (x[gix("pl",r1,r2To,p)]-x[gix("pl",r1,r2,p)]-ct);
    }
  CONSTRAIN_END_LOOP

  // mkteq25(i,p_tr).. sum(p_pr, a(p_pr,p_tr)*RVAR('pl',i,p_pr))+m(i,p_tr)  =g=  (RVAR('pl',i,p_tr));        // price of raw products + transf cost > trasf product
  CONSTRAIN_START_LOOP(secPr,13)
    mv = (gpd("m",l2r[r1][r2],secPr[p]) * gpd("pol_trSub",l2r[r1][r2],secPr[p]));
    g[cix] = mv - x[gix("pl",r1,r2,p+nPriPr)];
    for (uint p2=0;p2<priPr.size();p2++){
      a = gpd("a",l2r[r1][r2],priPr[p2],DATA_NOW,secPr[p]);
      g[cix] += a * x[gix("pl",r1,r2,p2)];
    }
  CONSTRAIN_END_LOOP

//  // mkteq18(i,p_pr).. RVAR('sa',i,p_pr)+RVAR('sl',i,p_pr)  =l=  dispor(i,p_pr); // total supply lower than the available stock
//  CONSTRAIN_START_LOOP(priPr,14)
//    in = gpd("in",l2r[r1][r2],priPr[p]);
//    double d1 = gix("sa",r1,r2,p);
//    double d2 = gix("sl",r1,r2,p);
//    g[cix] = x[gix("sa",r1,r2,p)]+x[gix("sl",r1,r2,p)]-in;
//  CONSTRAIN_END_LOOP

  // resbounds(i, p_pr_comb).. RVAR('sa',i,p_pr)+RVAR('sl',i,p_pr)  =l=  dispor(i,p_pr); // total supply lower than the available stock - FOR all combination subsets of ins
  CONSTRAIN_START_LOOP(priPrCombs,14)
    //ModelRegion* REG = MTHREAD->MD->getRegion(l2r[r1][r2]); // possibly slower
    //in = REG->inResByAnyCombination[p];
    in = ins[r1][r2][p];
    //if(p==0){
    //  in = 1.0; // workaround to lead -1<0 rather than 0<0 for the first (empty) subset - notneeded
    //}
    g[cix] = -in;
    for(uint i=0;i<priPrCombs[p].size();i++){
      g[cix] += x[gix("sa",r1,r2,priPrCombs[p][i])]+x[gix("sl",r1,r2,priPrCombs[p][i])];
    }
    g[cix] -= overharvestingAllowance; //0.02 don't work always, expecially intermediate scnearios, 0.1 seems to work but produce a large artefact 20160219: made it a parameter

  CONSTRAIN_END_LOOP

  //CONSTRAIN_START_LOOP(priPr,15)
  //    g[cix] = x[gix("st",r1,r2,p)]-(x[gix("sl",r1,r2,p)]+x[gix("sa",r1,r2,p)]);
  //CONSTRAIN_END_LOOP

  return true;
}

bool eval_f ( Index  n, const Number *  x, bool  new_x, Number &  obj_value  )

virtual

Original method from Ipopt to return the objective value remains unchanged

Definition at line 822 of file Opt.cpp.

                                                                  {
  eval_obj(n,x,obj_value);

  return true;
}

bool eval_g ( Index  n, const Number *  x, bool  new_x, Index  m, Number *  g  )

virtual

Original method from Ipopt to return the constraint residuals remains unchanged

Definition at line 837 of file Opt.cpp.

                                                                   {

  eval_constraints(n,x,m,g);

  return true;
}

bool eval_grad_f ( Index  n, const Number *  x, bool  new_x, Number *  grad_f  )

virtual

Original method from Ipopt to return the gradient of the objective remains unchanged

Definition at line 829 of file Opt.cpp.

                                                                    {

  gradient(tag_f,n,x,grad_f);

  return true;
}

bool eval_h ( Index  n, const Number *  x, bool  new_x, Number  obj_factor, Index  m, const Number *  lambda, bool  new_lambda, Index  nele_hess, Index *  iRow, Index *  jCol, Number *  values  )

virtual

Original method from Ipopt to return: 1) The structure of the hessian of the lagrangian (if "values" is NULL) 2) The values of the hessian of the lagrangian (if "values" is not NULL)remains unchanged

Definition at line 871 of file Opt.cpp.

                                                                                       {

  //return false;

  if (values == NULL) {
    // return the structure. This is a symmetric matrix, fill the lower left
    // triangle only.

    for(Index idx=0; idx<nnz_L; idx++)
      {
  iRow[idx] = rind_L[idx];
  jCol[idx] = cind_L[idx];
      }
  }
  else {
    // return the values. This is a symmetric matrix, fill the lower left
    // triangle only

    for(Index idx = 0; idx<n ; idx++)
      x_lam[idx] = x[idx];
    for(Index idx = 0; idx<m ; idx++)
      x_lam[n+idx] = lambda[idx];
    x_lam[n+m] = obj_factor;

    sparse_hess(tag_L, n+m+1, 1, x_lam, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);

    Index idx = 0;
    for(Index idx_total = 0; idx_total <nnz_L_total ; idx_total++)
      {
  if((rind_L_total[idx_total] < (unsigned int) n) && (cind_L_total[idx_total] < (unsigned int) n))
    {
      values[idx] = hessval[idx_total];
      idx++;
    }
      }
  }

  return true;

  //return false;
}

bool eval_jac_g ( Index  n, const Number *  x, bool  new_x, Index  m, Index  nele_jac, Index *  iRow, Index *  jCol, Number *  values  )

virtual

Original method from Ipopt to return: 1) The structure of the jacobian (if "values" is NULL) 2) The values of the jacobian (if "values" is not NULL)remains unchanged

Definition at line 845 of file Opt.cpp.

                                                         {
  if (values == NULL) {
    // return the structure of the jacobian

    for(Index idx=0; idx<nnz_jac; idx++)
      {
  iRow[idx] = rind_g[idx];
  jCol[idx] = cind_g[idx];
      }
  }
  else {
    // return the values of the jacobian of the constraints

    sparse_jac(tag_g, m, n, 1, x, &nnz_jac, &rind_g, &cind_g, &jacval, options_g);

    for(Index idx=0; idx<nnz_jac; idx++)
      {
  values[idx] = jacval[idx];

      }
  }
  return true;
}

bool eval_obj	(	Index	n,
		const T *	x,
		T &	obj_value
	)

Template to return the objective value

Define the objective function

Definition at line 220 of file Opt.cpp.

                                              {

  double aa, bb, dc0, sigma, a_pr, ct, m, zeromax,supCorr2;
  obj_value = 0.;
  zeromax = 0.;

  for (uint r1=0;r1<l2r.size();r1++){
    for (uint r2=0;r2<l2r[r1].size();r2++){
    //  // consumer's surplus..
    //  sum ((i,p_tr),
    //    AA(i,p_tr)*(RVAR('dc',i,p_tr)**((sigma(p_tr)+1)/sigma(p_tr)))
    //    - AA(i,p_tr)*((0.5*dc0(i,p_tr))**((sigma(p_tr)+1)/sigma(p_tr)))
    //    - RVAR('pc',i,p_tr)*RVAR('dc',i,p_tr)
    //  )
    // 20161003: TODO: check if subsidies should enter also the obj function other than the bounds equations. For the moment, as agreed with Sylvain, they are left outside the obj function, but I am not sure of it.
    // 20170306: Subsides (but not substitute products) added to the aa and bb parameters. They do not change anything.
    // 20170307: obj_value does not account for the minus 0-5dc0 part anymore. It doesn't chang anything if not that the "reported" surplus is negative.
    //           The integral of a power function p=q^a with a < -1 (derived from an q=p^b function with b between 0 and -1) is negative ranging from -inf at q=0
    //           to 0- with q = +inf
      for (uint p=0;p<secPr.size();p++){
        aa = gpd("aa",l2r[r1][r2],secPr[p]);
        sigma = gpd("sigma",l2r[r1][r2],secPr[p]);
        dc0 = gpd("dc",l2r[r1][r2],secPr[p],secondYear);
        //obj_value += aa*pow(mymax(zeromax,x[gix("dc",r1,r2,p)]),(sigma+1)/sigma)-aa*pow(mymax(zeromax,0.5*dc0),(sigma+1)/sigma)-x[gix("pc",r1,r2,p+nPriPr)]*x[gix("dc",r1,r2,p)];
        obj_value += aa*pow(mymax(zeromax,x[gix("dc",r1,r2,p)]),(sigma+1)/sigma)-x[gix("pc",r1,r2,p+nPriPr)]*x[gix("dc",r1,r2,p)];

      }
      //  // producers surplus..
      //  + sum((i,p_pr),
      //    RVAR('pc',i,p_pr)*RVAR('sc',i,p_pr)
      //    - BB(i,p_pr)*(RVAR('sc',i,p_pr)**((sigma(p_pr)+1)/sigma(p_pr)))
      //  )
      for (uint p=0;p<priPr.size();p++){
        bb = gpd("bb",l2r[r1][r2],priPr[p]);
        sigma = gpd("sigma",l2r[r1][r2],priPr[p]);
        //supCorr2 = gpd("supCorr2",l2r[r1][r2],priPr[p]);
        obj_value += x[gix("pc",r1,r2,p)]*x[gix("sc",r1,r2,p)] - bb*pow(mymax(zeromax,x[gix("sc",r1,r2,p)]),((sigma+1)/sigma));
      }
      //  // trasformations between primary products
      //  + sum ((i,p_pr,p_pr2),
      //  +RVAR('pc',i,p_pr2)*pres(p_pr,p_pr2)*RVAR('sc',i,p_pr)
      //  -BB(i,p_pr2)*(pres(p_pr,p_pr2)*RVAR('sc',i,p_pr))**((sigma(p_pr2)+1)/sigma(p_pr2))
      //   )

      for (uint p1=0;p1<priPr.size();p1++){
        for (uint p2=0;p2<priPr.size();p2++){
          a_pr = gpd("a_pr",l2r[r1][r2],priPr[p1],DATA_NOW,priPr[p2]);
          bb = gpd("bb",l2r[r1][r2],priPr[p2]);
          sigma = gpd("sigma",l2r[r1][r2],priPr[p2]);
          obj_value += x[gix("pc",r1,r2,p2)]*a_pr*x[gix("sc",r1,r2,p1)]-bb*pow(mymax(zeromax,a_pr*x[gix("sc",r1,r2,p1)]),(sigma+1)/sigma);
        }
      }
      //  // surplus of transport agents..
      //  + sum((i,j,prd), (RVAR('pl',j,prd)-RVAR('pl',i,prd)-CT(i,j,prd))*EXP(i,j,prd))
      for (uint p=0;p<allPr.size();p++){
        for (uint r2To=0;r2To<l2r[r1].size();r2To++){
          ct = (gpd("ct",l2r[r1][r2],allPr[p],DATA_NOW,i2s(l2r[r1][r2To])) * gpd("pol_tcSub",l2r[r1][r2],allPr[p]));
          obj_value += (x[gix("pl",r1,r2To,p)]-x[gix("pl",r1,r2,p)]-ct)*x[gix("rt",r1,r2,p,r2To)];
        }
      }

      //  // transformers surplus..
      //  + sum((i,p_tr), (RVAR('pl',i,p_tr)-m(i,p_tr))*(RVAR('sl',i,p_tr))) // attenction it's local. if we include w imports or p expports this have to change
      for (uint p=0;p<secPr.size();p++){
        m = (gpd("m",l2r[r1][r2],secPr[p]) * gpd("pol_trSub",l2r[r1][r2],secPr[p]));
        obj_value += (x[gix("pl",r1,r2,p+nPriPr)]-m)*x[gix("sl",r1,r2,p+nPriPr)];
      }
      //  - sum((i,p_pr), RVAR('pl',i,p_pr)*RVAR('dl',i,p_pr))               // to total and an other equation total=local+abroad should be added
      for (uint p=0;p<priPr.size();p++){
        obj_value -= x[gix("pl",r1,r2,p)]*x[gix("dl",r1,r2,p)];
      }
    } // end of each lev2 regions

  } //end of each r1 regions

  //obj_value = -obj_value; // we want maximisation, ipopt minimize! (donei n the options - scaling obj function)

 //exit(0);
  return true;
  // checked 20120802 this function is ok with gams, both in input and in output of the preoptimisation stage

}

void finalize_solution ( SolverReturn  status, Index  n, const Number *  x, const Number *  z_L, const Number *  z_U, Index  m, const Number *  g, const Number *  lambda, Number  obj_value, const IpoptData *  ip_data, IpoptCalculatedQuantities *  ip_cq  )

virtual

This method is called when the algorithm is complete so the TNLP can store/write the solution

Definition at line 732 of file Opt.cpp.

                                                                                                    {

  printf("\n\nObjective value\n");
  printf("f(x*) = %e\n", obj_value);

  // --> here is where to code the assignment of optimal values to to spd()

  VarMap::iterator viter;

  // fixing the starting points for each variable at the level of the previous years
  for (viter = vars.begin(); viter != vars.end(); ++viter) {
    //string debugs = viter->first;
    int vdomtype = viter->second.domain;
    if(vdomtype==DOM_PRI_PR){
      for(uint r1=0;r1<l2r.size();r1++){
        for(uint r2=0;r2<l2r[r1].size();r2++){
          for(uint p=0;p<priPr.size();p++){
            spd(x[gix(viter->first,r1,r2,p)],viter->first,l2r[r1][r2],priPr[p]);
          }
        }
      }
    } else if (vdomtype==DOM_SEC_PR) {
      for(uint r1=0;r1<l2r.size();r1++){
        for(uint r2=0;r2<l2r[r1].size();r2++){
          for(uint p=0;p<secPr.size();p++){
            spd(x[gix(viter->first,r1,r2,p)],viter->first,l2r[r1][r2],secPr[p]);

          }
        }
      }
    } else if (vdomtype==DOM_ALL_PR) {
      for(uint r1=0;r1<l2r.size();r1++){
        for(uint r2=0;r2<l2r[r1].size();r2++){
          for(uint p=0;p<allPr.size();p++){
            spd(x[gix(viter->first,r1,r2,p)],viter->first,l2r[r1][r2],allPr[p]);
          }
        }
      }
    } else if (vdomtype==DOM_R2_ALL_PR) {
      for(uint r1=0;r1<l2r.size();r1++){
        for(uint r2=0;r2<l2r[r1].size();r2++){
          for(uint p=0;p<allPr.size();p++){
            for(uint r2To=0;r2To<l2r[r1].size();r2To++){
              //if(x[gix(viter->first,r1,r2,p,r2To)] > 0){
              //  cout << l2r[r1][r2] << "\t" << allPr[p] << "\t" << l2r[r1][r2To] << "\t" << x[gix(viter->first,r1,r2,p,r2To)] << endl;
              //}
              spd(x[gix(viter->first,r1,r2,p,r2To)],viter->first,l2r[r1][r2],allPr[p],DATA_NOW,false,i2s(l2r[r1][r2To]));
            }
          }
        }
      }
    } else {
      msgOut(MSG_CRITICAL_ERROR,"Try to setting the solved value of a variable of unknow type ("+viter->first+")");
    }
  }

    // memory deallocation of ADOL-C variables
    delete[] x_lam;

    free(rind_g);
    free(cind_g);

    delete[] rind_L;
    delete[] cind_L;

    free(rind_L_total);
    free(cind_L_total);
    free(jacval);
    free(hessval);

    for (int i=0;i<n+m+1;i++) {
        free(HP_t[i]);
    }
    free(HP_t);

}

const int gdt ( const string &  varName )

protected

Get the domain type of a given variable.

Definition at line 1292 of file Opt.cpp.

                             { // get domain type
  VarMap::const_iterator p;
  p=vars.find(varName);
  if(p != vars.end()) {
    return p->second.domain;
  }
  else {
    msgOut(MSG_CRITICAL_ERROR, "Asking the domain type of a variable ("+varName+") that doesn't exist!");
    return 0;
  }
}

const int gdt ( const int &  cn )

protected

Get the domain type of a given constrain.

Definition at line 1305 of file Opt.cpp.

                     { // get domain type
  return cons.at(cn).domain;
}

void generate_tapes ( Index  n, Index  m, Index &  nnz_jac_g, Index &  nnz_h_lag  )

virtual

Method to generate the required tapes

Copied from http://bocop.org/

test avec "<=" (avant on avait "<" : bug, acces memoire non allouee)

valgrind : invalid read

test avec "<=" (pour etre coherent avec la remarque ci dessus, mais pas de cas test, a verifier)

test

test

Definition at line 918 of file Opt.cpp.

                                                                       {
    /// Copied from http://bocop.org/
    Number *xp    = new double[n];
    Number *lamp  = new double[m];
    Number *zl    = new double[m];
    Number *zu    = new double[m];

    adouble *xa   = new adouble[n];
    adouble *g    = new adouble[m];
    adouble *lam  = new adouble[m];
    adouble sig;
    adouble obj_value;

    double dummy;
//  double *jacval;

    int i,j,k,l,ii;

    x_lam   = new double[n+m+1];

//  cout << " Avant get_start" << endl;
    get_starting_point(n, 1, xp, 0, zl, zu, m, 0, lamp);
//   cout << " Apres get_start" << endl;

    //if(initOpt){ // that's funny, if I use this I get it slighly longer times, whatever I then use trace_off() or trace_off(1) (save to disk, seems unnecessary). If I use regenerated tapes I have also slighly inaccurate results.
    trace_on(tag_f);

    for(Index idx=0;idx<n;idx++)
        xa[idx] <<= xp[idx];

    eval_obj(n,xa,obj_value);

    obj_value >>= dummy;

    trace_off();

    trace_on(tag_g);

    for(Index idx=0;idx<n;idx++)
        xa[idx] <<= xp[idx];

    eval_constraints(n,xa,m,g);


    for(Index idx=0;idx<m;idx++)
        g[idx] >>= dummy;

    trace_off();

    trace_on(tag_L);

    for(Index idx=0;idx<n;idx++)
        xa[idx] <<= xp[idx];
    for(Index idx=0;idx<m;idx++)
        lam[idx] <<= 1.0;
    sig <<= 1.0;

    eval_obj(n,xa,obj_value);

    obj_value *= sig;
    eval_constraints(n,xa,m,g);

    for(Index idx=0;idx<m;idx++)
        obj_value += g[idx]*lam[idx];

    obj_value >>= dummy;

    trace_off();
    //} // end of if initOpt()



    rind_g = NULL;
    cind_g = NULL;

    options_g[0] = 0;          /* sparsity pattern by index domains (default) */
    options_g[1] = 0;          /*                         safe mode (default) */
    options_g[2] = -1;         /*                     &jacval is not computed */
    options_g[3] = 0;          /*                column compression (default) */

    jacval=NULL;

    sparse_jac(tag_g, m, n, 0, xp, &nnz_jac, &rind_g, &cind_g, &jacval, options_g);

    options_g[2] = 0;
    nnz_jac_g = nnz_jac;

    unsigned int  **JP_f=NULL;                /* compressed block row storage */
    unsigned int  **JP_g=NULL;                /* compressed block row storage */
    unsigned int  **HP_f=NULL;                /* compressed block row storage */
    unsigned int  **HP_g=NULL;                /* compressed block row storage */
    unsigned int  *HP_length=NULL;            /* length of arrays */
    unsigned int  *temp=NULL;                 /* help array */

    int ctrl_H;

    JP_f = (unsigned int **) malloc(sizeof(unsigned int*));
    JP_g = (unsigned int **) malloc(m*sizeof(unsigned int*));
    HP_f = (unsigned int **) malloc(n*sizeof(unsigned int*));
    HP_g = (unsigned int **) malloc(n*sizeof(unsigned int*));
    HP_t = (unsigned int **) malloc((n+m+1)*sizeof(unsigned int*));
    HP_length = (unsigned int *) malloc((n)*sizeof(unsigned int));
    ctrl_H = 0;

    hess_pat(tag_f, n, xp, HP_f, ctrl_H);

    indopro_forward_safe(tag_f, 1, n, xp, JP_f);
    indopro_forward_safe(tag_g, m, n, xp, JP_g);
    nonl_ind_forward_safe(tag_g, m, n, xp, HP_g);

    for (i=0;i<n;i++)
    {
        if (HP_f[i][0]+HP_g[i][0]!=0)
        {
            if (HP_f[i][0]==0) // number of non zeros in the i-th row
            {
                HP_t[i] = (unsigned int *) malloc((HP_g[i][0]+HPOFF)*sizeof(unsigned int));
                for(j=0;j<=(int) HP_g[i][0];j++)
                {
                    HP_t[i][j] = HP_g[i][j];
                }
                HP_length[i] = HP_g[i][0]+HPOFF;
            }
            else
            {
                if (HP_g[i][0]==0) // number of non zeros in the i-th row
                {
                    HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HPOFF)*sizeof(unsigned int));
                    for(j=0;j<=(int) HP_f[i][0];j++)
                    {
                        HP_t[i][j] = HP_f[i][j];
                    }
                    HP_length[i] = HP_f[i][0]+HPOFF;
                }
                else
                {
                    HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HP_g[i][0]+HPOFF)*sizeof(unsigned int));
                    k = l = j = 1;
                    while ((k<=(int) HP_f[i][0]) && (l <= (int) HP_g[i][0]))
                    {
                        if (HP_f[i][k] < HP_g[i][l])
                        {
                            HP_t[i][j]=HP_f[i][k];
                            j++; k++;
                        }
                        else
                        {
                            if (HP_f[i][k] == HP_g[i][l])
                            {
                                HP_t[i][j]=HP_f[i][k];
                                l++;j++;k++;
                            }
                            else
                            {
                                HP_t[i][j]=HP_g[i][l];
                                j++;l++;
                            }
                        }
                    } // end while

                    // Fill the end of the vector if HP_g[i][0] < HP_f[i][0]
                    for(ii=k;ii<=(int) HP_f[i][0];ii++)
                    {
                        HP_t[i][j] = HP_f[i][ii];
                        j++;
                    }

                    // Fill the end of the vector if HP_f[i][0] < HP_g[i][0]
                    for(ii=l;ii<=(int) HP_g[i][0];ii++)
                    {
                        HP_t[i][j] = HP_g[i][ii];
                        j++;
                    }

                }
            }
            HP_t[i][0]=j-1; // set the first element with the number of non zeros in the i-th line
            HP_length[i] = HP_f[i][0]+HP_g[i][0]+HPOFF; // length of the i-th line
        }
        else
        {
            HP_t[i] = (unsigned int *) malloc((HPOFF+1)*sizeof(unsigned int));
            HP_t[i][0]=0;
            HP_length[i]=HPOFF;
        }

//     if (i==(int)n-1)
//     {
//       cout << " DISPLAY FINAL TIME HP : " << endl;
//       for (ii=0;ii<=(int)HP_length[i];ii++)
//         cout << " -------> HP[last][" << ii << "] = " << HP_t[i][ii] << endl;
//     }
    }

//   cout << " Avant les boucles" << endl;
//   cout << " m = " << m << endl;

    for (i=0;i<m;i++)
    {
//     cout << i << " --> nnz JP_g = " << JP_g[i][0]+1 << " -- ";
        HP_t[n+i] = (unsigned int *) malloc((JP_g[i][0]+1)*sizeof(unsigned int));
        HP_t[n+i][0]=JP_g[i][0];

//     cout << HP_t[n+i][0] << endl;

        for(j=1;j<= (int) JP_g[i][0];j++)
        {
            HP_t[n+i][j]=JP_g[i][j];
//       cout << " ---------> " << HP_t[n+i][j] << endl;
//       cout << " --> HP_length[" << JP_g[i][j] << "] = " << HP_length[JP_g[i][j]] << "  --  HP_t[" << JP_g[i][j] << "][0] = " << HP_t[JP_g[i][j]][0]+1 << endl;
            // We write the rows allocated in the previous "for" loop
            // If the memory allocated for the row is not big enough :
            if (HP_length[JP_g[i][j]] <= HP_t[JP_g[i][j]][0]+1) //! test avec "<=" (avant on avait "<" : bug, acces memoire non allouee)
            {
//         cout << " ---------> WARNING " << endl;
//         cout << " At index " << JP_g[i][j] << endl;


                // save a copy of existing vector elements :
                temp = (unsigned int *) malloc((HP_t[JP_g[i][j]][0]+1)*sizeof(unsigned int));
                for(l=0;l<=(int)HP_t[JP_g[i][j]][0];l++)
                {
                    temp[l] = HP_t[JP_g[i][j]][l]; //! valgrind : invalid read
//           cout << " -------> l = " << l << " -- " << temp[l] << endl;
                }

//         cout << " ----------->  DISPLAY " << endl;
//         for(l=0;l<=(int)HP_t[JP_g[i][j]][0];l++)
//         {
//           temp[l] = HP_t[JP_g[i][j]][l]; //! valgrind : invalid read & write
//           cout << " -------> HP[machin][" << l << "] = " << HP_t[JP_g[i][j]][l] << endl; //! valgrind : invalid read
//         }


                // Free existing row, and allocate more memory for it :
//         cout << " Avant free  --> pointeur = " <<HP_t[JP_g[i][j]]<< endl;
                unsigned int machin = JP_g[i][j];
                free(HP_t[machin]); // !Problem double free or corruption
//         cout << " Apres free --> pointeur = " <<HP_t[JP_g[i][j]]<< endl;

                HP_t[JP_g[i][j]] = (unsigned int *) malloc(2*HP_length[JP_g[i][j]]*sizeof(unsigned int));
                HP_length[JP_g[i][j]] = 2*HP_length[JP_g[i][j]];

                // Put back the values in this bigger vector :
                for(l=0;l<=(int)temp[0];l++)
                    HP_t[JP_g[i][j]][l] =temp[l];
                free(temp);

//         HP_t[JP_g[i][j]] = (unsigned int*) realloc (HP_t[JP_g[i][j]], 2*HP_length[JP_g[i][j]] * sizeof(unsigned int));
//         HP_length[JP_g[i][j]] = 2*HP_length[JP_g[i][j]];
            }
            HP_t[JP_g[i][j]][0] = HP_t[JP_g[i][j]][0]+1; // The size of the row is one greater than before
            HP_t[JP_g[i][j]][HP_t[JP_g[i][j]][0]] = i+n; // Now adding the element at the end //! valgrind : invalid write
        }
    }
//   cout << " Apres les boucles" << endl;

    for(j=1;j<= (int) JP_f[0][0];j++)
    {
        if (HP_length[JP_f[0][j]] <= HP_t[JP_f[0][j]][0]+1) //! test avec "<=" (pour etre coherent avec la remarque ci dessus, mais pas de cas test, a verifier)
        {
            temp = (unsigned int *) malloc((HP_t[JP_f[0][j]][0])*sizeof(unsigned int));
            for(l=0;l<=(int)HP_t[JP_f[0][j]][0];l++)
                temp[l] = HP_t[JP_f[0][j]][l];
            free(HP_t[JP_f[0][j]]);
            HP_t[JP_f[0][j]] = (unsigned int *) malloc(2*HP_length[JP_f[0][j]]*sizeof(unsigned int));
            HP_length[JP_f[0][j]] = 2*HP_length[JP_f[0][j]];
            for(l=0;l<=(int)temp[0];l++)
                HP_t[JP_f[0][j]][l] =temp[l];
            free(temp);
        }
        HP_t[JP_f[0][j]][0] = HP_t[JP_f[0][j]][0]+1;
        HP_t[JP_f[0][j]][HP_t[JP_f[0][j]][0]] = n+m;
    }

    HP_t[n+m] = (unsigned int *) malloc((JP_f[0][0]+2)*sizeof(unsigned int));
    HP_t[n+m][0]=JP_f[0][0]+1;
    for(j=1;j<= (int) JP_f[0][0];j++)
        HP_t[n+m][j]=JP_f[0][j];
    HP_t[n+m][JP_f[0][0]+1]=n+m;

    set_HP(tag_L,n+m+1,HP_t); // set sparsity pattern for the Hessian

    nnz_h_lag = 0;
    for (i=0;i<n;i++)
    {
        for (j=1;j<=(int) HP_t[i][0];j++)
            if ((int) HP_t[i][j] <= i)
                nnz_h_lag++;
        free(HP_f[i]);
        free(HP_g[i]);
    }
    nnz_L = nnz_h_lag;

    options_L[0] = 0;
    options_L[1] = 1;

    rind_L_total = NULL;
    cind_L_total = NULL;
    hessval = NULL;

    sparse_hess(tag_L, n+m+1, -1, xp, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);

    rind_L = new unsigned int[nnz_L];
    cind_L = new unsigned int[nnz_L];
    rind_L_total = (unsigned int*) malloc(nnz_L_total*sizeof(unsigned int)); //! test
    cind_L_total = (unsigned int*) malloc(nnz_L_total*sizeof(unsigned int)); //! test

    unsigned int ind = 0;

    for (int i=0;i<n;i++)
        for (unsigned int j=1;j<=HP_t[i][0];j++)
    {
        if (((int) HP_t[i][j]>=i) &&((int) HP_t[i][j]<n))
        {
            rind_L[ind] = i;
            cind_L[ind++] = HP_t[i][j];
        }
    }

    ind = 0;
    for (int i=0;i<n+m+1;i++)
        for (unsigned int j=1;j<=HP_t[i][0];j++)
    {
        if ((int) HP_t[i][j]>=i)
        {
            rind_L_total[ind] = i;
            cind_L_total[ind++] = HP_t[i][j];
        }
    }

    for (i=0;i<m;i++) {
        free(JP_g[i]);
    }

    free(JP_f[0]);
    free(JP_f);
    free(JP_g);
    free(HP_f);
    free(HP_g);
    free(HP_length);

    delete[] lam;
    delete[] g;
    delete[] xa;
    delete[] zu;
    delete[] zl;
    delete[] lamp;
    delete[] xp;

}

bool get_bounds_info ( Index  n, Number *  x_l, Number *  x_u, Index  m, Number *  g_l, Number *  g_u  )

virtual

Method to return the bounds for my problem

Definition at line 634 of file Opt.cpp.

                                                                                        {

  // Set the bounds for the endogenous variables..
  for (Index i=0; i<n; i++) {
    x_l[i] =  getBoundByIndex(LBOUND,i);
    x_u[i] =  getBoundByIndex(UBOUND,i);
  }

  // Set the bounds for the constraints..
  for (Index i=0; i<m; i++) {
    int direction = getConstrainDirectionByIndex(i);
    switch (direction){
      case CONSTR_EQ:
        g_l[i] = 0.;
        g_u[i] = 0.;
        break;
      case CONSTR_LE0:
        g_l[i] = -2e19;
        g_u[i] = 0.;
        break;
      case CONSTR_GE0:
        g_l[i] = 0.;
        g_u[i] = 2e19;
        break;
    }
  }
  return true;
}

bool get_nlp_info ( Index &  n, Index &  m, Index &  nnz_jac_g, Index &  nnz_h_lag, IndexStyleEnum &  index_style  )

virtual

Method to return some info about the nlp

Definition at line 551 of file Opt.cpp.

                                                                {


  if(initOpt){
    // does this initialisation code only once
    priPr = MTHREAD->MD->getStringVectorSetting("priProducts");
    secPr = MTHREAD->MD->getStringVectorSetting("secProducts");
    othPr = MTHREAD->MD->getStringVectorSetting("othExogenousProducts");
    allPr = priPr;
    allPr.insert( allPr.end(), secPr.begin(), secPr.end() );
    nPriPr = priPr.size();
    nSecPr = secPr.size();
    nAllPr = allPr.size();
    nOthPr = othPr.size();
    std::vector<int> l1regIds =  MTHREAD->MD->getRegionIds(1, true);
    nL2r =  MTHREAD->MD->getRegionIds(2, true).size();
    firstYear = MTHREAD->MD->getIntSetting("initialYear");
    secondYear = firstYear+1;
    worldCodeLev2 = MTHREAD->MD->getIntSetting("worldCodeLev2");

    for(uint i=0;i<l1regIds.size();i++){
      std::vector<int> l2ChildrenIds;
      ModelRegion* l1Region = MTHREAD->MD->getRegion(l1regIds[i]);
      std::vector<ModelRegion*> l2Childrens = l1Region->getChildren(true);
      for(uint j=0;j<l2Childrens.size();j++){
        l2ChildrenIds.push_back(l2Childrens[j]->getRegId());
      }
      if(l2ChildrenIds.size()){
        l2r.push_back(l2ChildrenIds);
      }
    }

    // Create a vector with all possible combinations of primary products
    priPrCombs = MTHREAD->MD->createCombinationsVector(nPriPr);
    nPriPrCombs = priPrCombs.size();

    // put the variables and their domain in the vars map
    declareVariables();

    // declaring the contrains...
    declareConstrains();

    // calculate number of variables and constrains..
    calculateNumberVariablesConstrains();

    // cache initial positions (variables and constrains)..
    cacheInitialPosition();

    // cache initial positions (variables and constrains)..
    cachePositions();

    //tempDebug();

    //debugPrintParameters();

  } // finish initialisation things to be done only the first year

  previousYear = MTHREAD->SCD->getYear()-1; // this has to be done EVERY years !!

  n = nVar; // 300; // nVar;
  m = nCons; // 70; // nCons;

  overharvestingAllowance = MTHREAD->MD->getDoubleSetting("overharvestingAllowance",DATA_NOW);

  copyInventoryResourses();

  generate_tapes(n, m, nnz_jac_g, nnz_h_lag);

  //if(initOpt){
  //  calculateSparsityPatternJ();
  //  calculateSparsityPatternH();
    //tempDebug();
  //}

  // use the C style indexing (0-based)
  index_style = C_STYLE;

  initOpt=false;
  return true;
}

Here is the call graph for this function:

bool get_starting_point ( Index  n, bool  init_x, Number *  x, bool  init_z, Number *  z_L, Number *  z_U, Index  m, bool  init_lambda, Number *  lambda  )

virtual

Method to return the starting point for the algorithm

Definition at line 664 of file Opt.cpp.

                                                                  {

  // function checked on 20120724 on a subset of 3 regions and 4 products. All variables initial values are correctly those outputed by gams in 2006.
  //int thisYear = MTHREAD->SCD->getYear();
  //int initialOptYear = MTHREAD->MD->getIntSetting("initialOptYear");
  //if(thisYear != initialOptYear) return true;

  //msgOut(MSG_DEBUG,"Giving optimising variables previous years value as starting point");
  // Here, we assume we only have starting values for x, if you code
  // your own NLP, you can provide starting values for the others if
  // you wish.
  assert(init_x == true);
  assert(init_z == false);
  assert(init_lambda == false);

  VarMap::iterator viter;

  // fixing the starting points for each variable at the level of the previous years
  for (viter = vars.begin(); viter != vars.end(); ++viter) {
    //string debugs = viter->first;
    int vdomtype = viter->second.domain;
    if(vdomtype==DOM_PRI_PR){
      for(uint r1=0;r1<l2r.size();r1++){
        for(uint r2=0;r2<l2r[r1].size();r2++){
          for(uint p=0;p<priPr.size();p++){
            x[gix(viter->first,r1,r2,p)]= gpd(viter->first,l2r[r1][r2],priPr[p],previousYear);
          }
        }
      }
    } else if (vdomtype==DOM_SEC_PR) {
      for(uint r1=0;r1<l2r.size();r1++){
        for(uint r2=0;r2<l2r[r1].size();r2++){
          for(uint p=0;p<secPr.size();p++){
            x[gix(viter->first,r1,r2,p)]= gpd(viter->first,l2r[r1][r2],secPr[p],previousYear);
          }
        }
      }
    } else if (vdomtype==DOM_ALL_PR) {
      for(uint r1=0;r1<l2r.size();r1++){
        for(uint r2=0;r2<l2r[r1].size();r2++){
          for(uint p=0;p<allPr.size();p++){
            x[gix(viter->first,r1,r2,p)]= gpd(viter->first,l2r[r1][r2],allPr[p],previousYear);
          }
        }
      }
    } else if (vdomtype==DOM_R2_ALL_PR) {
      for(uint r1=0;r1<l2r.size();r1++){
        for(uint r2=0;r2<l2r[r1].size();r2++){
          for(uint p=0;p<allPr.size();p++){
            for(uint r2To=0;r2To<l2r[r1].size();r2To++){
              x[gix(viter->first,r1,r2,p,r2To)]= gpd(viter->first,l2r[r1][r2],allPr[p],previousYear,i2s(l2r[r1][r2To]));
            }
          }
        }
      }
    } else {
      msgOut(MSG_CRITICAL_ERROR,"Try to setting the initial value of a variable of unknow type ("+viter->first+")");
    }
  }

  //msgOut(MSG_DEBUG,"Finisced initial value assignments");

  return true;
}

double getBoundByIndex ( const int &  bound_type, const int &  idx  )

protected

Return the bound of a given variable (by index)

Definition at line 1538 of file Opt.cpp.

                                                           {

  map <int, string>::const_iterator p;
  p=initPos_rev.upper_bound(idx);
  p--;
  VarMap::const_iterator p2;
  p2=vars.find(p->second);
  if(p2 != vars.end()) {
    if (bound_type==LBOUND){
      if (p2->second.l_bound_var == ""){ // this var don't specific a variable as bound
        return p2->second.l_bound;
      } else {
        return getDetailedBoundByVarAndIndex(p2->second,idx,LBOUND);
      }
    } else if (bound_type==UBOUND){
      if (p2->second.u_bound_var == ""){ // this var don't specific a variable as bound
        return p2->second.u_bound;
      } else {
        return getDetailedBoundByVarAndIndex(p2->second,idx,UBOUND);
      }
    } else {
      msgOut(MSG_CRITICAL_ERROR, "Asking the bound with a type ("+i2s(bound_type)+") that I don't know how to handle !");
    }
  }
  else {
    msgOut(MSG_CRITICAL_ERROR, "Asking the bound from a variable ("+p->second+") that doesn't exist!");
  }
  return 0.;
}

int getConNumber ( constrain *  con )

protected

Return the position in the cons vector.

Definition at line 1666 of file Opt.cpp.

                               {
  for(uint i=0;i<cons.size();i++){
    if(   cons[i].name       == con->name
       && cons[i].comment    == con->comment
       && cons[i].domain     == con->domain
       && cons[i].direction  == con->direction){
      return i;
       }
  }
  msgOut(MSG_CRITICAL_ERROR,"Constrain didn't found in list.");
}

constrain * getConstrainByIndex ( int  idx )

protected

Definition at line 1592 of file Opt.cpp.

                               {
  for(uint i=0;i<cons.size();i++){
    if(i!=cons.size()-1){
      if (idx >= gip(i) && idx < gip(i+1)){
        return &cons[i];
      }
    } else {
      if (idx >= gip(i) && idx < nCons){
      return &cons[i];
      }
    }
  }
  msgOut(MSG_CRITICAL_ERROR, "Asking contrain direction for an out of range contrain index!");
}

int getConstrainDirectionByIndex ( int  idx )

protected

Return the direction of a given constrain.

Definition at line 1522 of file Opt.cpp.

                                        {
  for(uint i=0;i<cons.size();i++){
    if(i!=cons.size()-1){
      if (idx >= gip(i) && idx < gip(i+1)){
        return cons[i].direction;
      }
    } else {
      if (idx >= gip(i) && idx < nCons){
      return cons[i].direction;
      }
    }
  }
  msgOut(MSG_CRITICAL_ERROR, "Asking contrain direction for an out of range contrain index!");
}

double getDetailedBoundByVarAndIndex ( const endvar &  var, const int &  idx, const int &  bType  )

protected

Return the bound of a given variable given the variable and the required index. Called by getBoundByIndex().

Definition at line 1569 of file Opt.cpp.

                                                                                        {
  // Tested 2015.01.08 with DOM_ALL_PR, DOM_PRI_PR, DOM_ALL_PR, DOM_R2_ALL_PR.
  int r1,r2,p,r2to;
  unpack(idx,var.domain,gip(var.name),r1,r2,p,r2to,true);
  //cout << "getBoundByVarAndIndex():\t" << var.name << '\t' << idx << '\t' << gip(var.name) << '\t' << r1 << '\t' << r2 << '\t' << p << '\t' << r2to << endl;
  //cout << "  --variables:\t" << var.l_bound_var << '\t' << var.u_bound_var << '\t' << "" << '\t' << l2r[r1][r2] << '\t' << "" << '\t' << allPr[p] << '\t' << l2r[r1][r2to] << endl;
  if(bType==LBOUND){
    if(r2to){
       return gpd(var.l_bound_var,l2r[r1][r2],allPr[p],DATA_NOW,i2s(l2r[r1][r2to]));
    } else {
      return gpd(var.l_bound_var,l2r[r1][r2],allPr[p],DATA_NOW,i2s(l2r[r1][r2to]));
    }
  } else {
    if(r2to){
       return gpd(var.u_bound_var,l2r[r1][r2],allPr[p]);
    } else {
      //cout << gpd(var.u_bound_var,l2r[r1][r2],allPr[p]) << endl;
      return gpd(var.u_bound_var,l2r[r1][r2],allPr[p]);
    }
  }
}

int getDomainElements

(

int

domain

)

return the number of elements of a domain

Definition at line 1488 of file Opt.cpp.

                                {
  int elements = 0;
  switch (domain){
    case DOM_PRI_PR:
      return nL2r*nPriPr;
    case DOM_SEC_PR:
      return nL2r*nSecPr;
    case DOM_ALL_PR:
      return nL2r*nAllPr;
    case DOM_R2_PRI_PR:
      for(uint r1=0;r1<l2r.size();r1++){
          elements += l2r[r1].size()*l2r[r1].size()*nPriPr; // EXP(i,j,p_pr)
      }
      return elements;
    case DOM_R2_SEC_PR:
      for(uint r1=0;r1<l2r.size();r1++){
          elements += l2r[r1].size()*l2r[r1].size()*nSecPr; // EXP(i,j,p_tr)
      }
      return elements;
    case DOM_R2_ALL_PR:
      for(uint r1=0;r1<l2r.size();r1++){
          elements += l2r[r1].size()*l2r[r1].size()*nAllPr; // EXP(i,j,prd)
      }
      return elements;
    case DOM_SCALAR:
      return 1;
    case DOM_PRI_PR_ALLCOMBS:
      return nL2r*nPriPrCombs;
  default:
    msgOut(MSG_CRITICAL_ERROR, "Asking for an unknown domain type ("+i2s(domain)+")");
  }
}

int getVarInstances

(

const string &

varName

)

build the matrix of the positions for a given variable or contrain

return the number of instances of a variable, given his domain type

Definition at line 1768 of file Opt.cpp.

                                         {
  return getDomainElements(gdt(varName));
}

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

inlineprotected

Definition at line 169 of file Opt.h.

{return MTHREAD->MD->getForData(type_h, regId_h, forType_h, diamClass_h, year);};

const int gip ( const string &  varName ) const

protected

Get the initial index position of a given variable in the concatenated array.

Definition at line 1274 of file Opt.cpp.

                                   { // get initial position
  map<string, int>::const_iterator p;
  p=initPos.find(varName);
  if(p != initPos.end()) {
    return p->second;
  }
  else {
    msgOut(MSG_CRITICAL_ERROR, "Asking the initial position in the concatenated array of a variable ("+varName+") that doesn't exist!");
    return 0;
  }
}

const int gip ( const int &  cn ) const

protected

Return the initial index position of a certain constrain.

Definition at line 1287 of file Opt.cpp.

                            { // get initial position
  return cInitPos.at(cn);
}

const int gix ( const string &  varName, const int &  r1Ix, const int &  r2Ix, const int &  prIx, const int &  r2IxTo = 0  ) const

protected

Get the index in the concatenated array gived a certain var name, the reg lev1 index, the reg lev2 index and the prod. index.

Definition at line 1351 of file Opt.cpp.

                                                                                                         {
  // attenction, for computational reasons we are not checking the call is within vectors limits!!!
  map <string, vector < vector < vector < vector <int> > > > >::const_iterator p;
  p=vpositions.find(varName);
  if(p != vpositions.end()) {
    return p->second[r1Ix][r2Ix][prIx][r2IxTo];
  }
  else {
    msgOut(MSG_CRITICAL_ERROR, "Asking the position of a variable ("+varName+") that doesn't exist!");
    return 0;
  }
}

const int gix ( const int &  cn, const int &  r1Ix, const int &  r2Ix, const int &  prIx, const int &  r2IxTo = 0  ) const

protected

Get the index in the concatenated array gived a certain constrain, the reg lev1 index, the reg lev2 index and the prod. index.

Definition at line 1365 of file Opt.cpp.

                                                                                                 {
  return cpositions[cn][r1Ix][r2Ix][prIx][r2IxTo];
}

const int gix_uncached ( const T &  v_or_c, int  r1Ix, int  r2Ix, int  prIx, int  r2IxTo = 0  )

protected

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.

Definition at line 1310 of file Opt.cpp.

                                                                          {

  // attenction, for computational reason we are not checking the call is within vectors limits!!!

  int dType = gdt(v_or_c);
  int othCountriesRegions = 0;
  int othCountriesRegions_r2case = 0;
  for (uint i=0;i<r1Ix;i++){
    othCountriesRegions += l2r[i].size();
  }
  for (uint i=0;i<r1Ix;i++){
    othCountriesRegions_r2case +=l2r[i].size()*l2r[i].size();
  }

  switch (dType){
    case DOM_PRI_PR:
      return gip(v_or_c)+(othCountriesRegions+r2Ix)*nPriPr+prIx;
    case DOM_SEC_PR:
      return gip(v_or_c)+(othCountriesRegions+r2Ix)*nSecPr+prIx;;
    case DOM_ALL_PR:
      return gip(v_or_c)+(othCountriesRegions+r2Ix)*nAllPr+prIx;
    case DOM_R2_PRI_PR:
      return gip(v_or_c)+(othCountriesRegions_r2case)*nAllPr+(r2Ix*nPriPr+prIx)*l2r[r1Ix].size()+r2IxTo;
    case DOM_R2_SEC_PR:
      return gip(v_or_c)+(othCountriesRegions_r2case)*nAllPr+(r2Ix*nSecPr+prIx)*l2r[r1Ix].size()+r2IxTo;
    case DOM_R2_ALL_PR:
      return gip(v_or_c)+(othCountriesRegions_r2case)*nAllPr+(r2Ix*nAllPr+prIx)*l2r[r1Ix].size()+r2IxTo; // new 20120814, looping r1,r2,p,r2to
      // initial position + (other countries region pairs + same country other regions from pair + regions to)* number of all products+product
      //return gip(v_or_c)+(othCountriesRegions_r2case+r2Ix*l2r[r1Ix].size()+r2IxTo)*nAllPr+prIx; // looping r1,r2,r2to,p
    case DOM_SCALAR:
      return gip(v_or_c);
    case DOM_PRI_PR_ALLCOMBS:
      return gip(v_or_c)+(othCountriesRegions+r2Ix)*nPriPrCombs+prIx;
  default:
    msgOut(MSG_CRITICAL_ERROR,"Try to calculate the position of a variable (or constrain) of unknow type.");
    return 0;
  }
}

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

inlineprotected

Definition at line 168 of file Opt.h.

{return MTHREAD->MD->getProdData(type_h, regId_h, prodId_h, year, freeDim_h);};

bool intermediate_callback ( AlgorithmMode  mode, Index  iter, Number  obj_value, Number  inf_pr, Number  inf_du, Number  mu, Number  d_norm, Number  regularization_size, Number  alpha_du, Number  alpha_pr, Index  ls_trials, const IpoptData *  ip_data, IpoptCalculatedQuantities *  ip_cq  )

virtual

Return information on each iteration

Definition at line 1759 of file Opt.cpp.

                                                                                                                                                                                                                                                                             {
    int itnumber = iter;
    if(itnumber%10==0){
        msgOut(MSG_DEBUG,"Running ("+i2s(itnumber)+" iter) ..");
    }
    return true;
}

const Number & mymax	(	const Number &	a,
		const Number &	b
	)

Definition at line 1749 of file Opt.cpp.

                                          {
  return (a<b)?b:a;
}

const adouble & mymax	(	const adouble &	a,
		const adouble &	b
	)

Definition at line 1753 of file Opt.cpp.

                                            {
  return (a<b)?b:a;
}

Opt& operator= ( const Opt &  )

protected

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

inlineprotected

Definition at line 171 of file Opt.h.

{MTHREAD->MD->setForData(value_h, type_h, regId_h, forType_h, diamClass_h, year, allowCreate);};

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

inlineprotected

Definition at line 170 of file Opt.h.

{MTHREAD->MD->setProdData(value_h, type_h, regId_h, prodId_h, year, allowCreate, freeDim_h);};

void tempDebug ( )

protected

Definition at line 1733 of file Opt.cpp.

              {

  cout << "Num of variables: " <<  nVar << " - Num of constrains:" << nCons << endl;
  cout << "IDX;ROW;COL" << endl;
  for(uint i=0;i<nzhelements.size();i++){
    cout << i << ";" << nzhelements[i][0] << ";" << nzhelements[i][1] << endl;
  }

  cout << "Dense jacobian: " << nCons * nVar << " elements" << endl;
  cout << "Dense hessian:  " << nVar*(nVar-1)/2+nVar << " elements" << endl;
  //exit(0);

}

void unpack ( int  ix_h, int  domain, int  initial, int &  r1_h, int &  r2_h, int &  p_h, int &  r2to_h, bool  fullp = false  )

protected

Return the dimensions given a certain index, domain type and initial position.

Definition at line 1609 of file Opt.cpp.

                                                                                                   {
  ix_h = ix_h-initial;
  double ix=0;
  bool r2flag = false;
  int pIndexToAdd = 0;
  int np=0;
  if(domain==DOM_PRI_PR || domain==DOM_R2_PRI_PR) {
    np = nPriPr;
  } else if (domain==DOM_SEC_PR || domain==DOM_R2_SEC_PR) {
    np = nSecPr;
  } else if (domain==DOM_ALL_PR || domain==DOM_R2_ALL_PR) {
    np = nAllPr;
  } else if (domain=DOM_SCALAR){
    r1_h=0;r2_h=0;p_h=0;r2to_h=0;
    return;
  } else {
    msgOut(MSG_CRITICAL_ERROR,"unknow domain ("+i2s(domain)+") in unpack() function.");
  }
  if(domain==DOM_R2_PRI_PR || domain==DOM_R2_SEC_PR ||domain==DOM_R2_ALL_PR){
    r2flag = true;
  }
  if(fullp && (domain==DOM_SEC_PR || domain==DOM_R2_SEC_PR)){ // changed 20140107 (any how, previously the unpack() function was not used!!)
    pIndexToAdd = nPriPr;
    //cout << "pindexToAdd: " << pIndexToAdd << endl;
  }

  for (uint r1=0;r1<l2r.size();r1++){
    for (uint r2=0;r2<l2r[r1].size();r2++){
      for (uint p=0;p<np;p++){
        if(!r2flag){
          if(ix==ix_h){
            r1_h=r1;
            r2_h=r2;
            p_h=p+pIndexToAdd;
            r2to_h=0;
            return;
          }
          ix++;
        } else {
          for (uint r2To=0;r2To<l2r[r1].size();r2To++){
            if(ix==ix_h){
              r1_h=r1;
              r2_h=r2;
              p_h=p+pIndexToAdd;
              r2to_h=r2To;
              return;
            }
            ix++;
          }
        }
      }
    }
  }
  msgOut(MSG_CRITICAL_ERROR, "Error in unpack() function. Ix ("+i2s(ix_h)+") can not be unpacked");
}

Member Data Documentation

vector<string> allPr

protected

Definition at line 195 of file Opt.h.

unsigned int* cind_g

protected

Definition at line 253 of file Opt.h.

unsigned int* cind_L

protected

Definition at line 256 of file Opt.h.

unsigned int* cind_L_total

protected

Definition at line 258 of file Opt.h.

vector<int> cInitPos

protected

A vector that returns the initial index position in the concatenated array for each constrain.

Definition at line 202 of file Opt.h.

vector<constrain> cons

protected

Definition at line 225 of file Opt.h.

vector< vector < vector < vector < vector <int> > > > > cpositions

protected

cached position in the concatenated vector for each variables. Dimensions are contrain number, l1reg, l2reg, prod, (l2To region).

Definition at line 205 of file Opt.h.

bool debugRunOnce

protected

Definition at line 221 of file Opt.h.

int firstYear

protected

Definition at line 218 of file Opt.h.

double* hessval

protected

Definition at line 259 of file Opt.h.

unsigned int** HP_t

protected

Definition at line 251 of file Opt.h.

bool initOpt

protected

Definition at line 224 of file Opt.h.

map<string, int> initPos

protected

A map that returns the initial index position in the concatenated array for each variable.

Definition at line 200 of file Opt.h.

map<int, string> initPos_rev

protected

A map with the name of the variable keyed by its initial position in the index.

Definition at line 201 of file Opt.h.

vector< vector < vector <double> > > ins

protected

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.

Definition at line 199 of file Opt.h.

double* jacval

protected

Definition at line 254 of file Opt.h.

vector< vector <int> > l2r

protected

Definition at line 197 of file Opt.h.

int nAllPr

protected

Definition at line 210 of file Opt.h.

int nCons

protected

Definition at line 213 of file Opt.h.

int nEqualityConstrains

protected

Definition at line 214 of file Opt.h.

int nGreaterEqualZeroConstrains

protected

Definition at line 216 of file Opt.h.

int nL2r

protected

Definition at line 211 of file Opt.h.

int nLowerEqualZeroConstrains

protected

Definition at line 215 of file Opt.h.

int nnz_jac

protected

Definition at line 260 of file Opt.h.

int nnz_L

protected

Definition at line 261 of file Opt.h.

int nnz_L_total

protected

Definition at line 261 of file Opt.h.

int nOthPr

protected

Definition at line 207 of file Opt.h.

int nPriPr

protected

Definition at line 206 of file Opt.h.

int nPriPrCombs

protected

Definition at line 208 of file Opt.h.

int nSecPr

protected

Definition at line 209 of file Opt.h.

int nVar

protected

Definition at line 212 of file Opt.h.

vector<vector <Index> > nzhelements

protected

nzero elements for the hessian matrix

Definition at line 227 of file Opt.h.

vector<vector <Index> > nzjelements

protected

nzero elements for the jacobian matrix. nzelements[i][0] -> row (constrain), nzelements[i][1] -> column (variable)

Definition at line 226 of file Opt.h.

int options_g[4]

protected

Definition at line 262 of file Opt.h.

int options_L[4]

protected

Definition at line 263 of file Opt.h.

vector<string> othPr

protected

Definition at line 196 of file Opt.h.

double overharvestingAllowance

protected

Allows to harvest more than the resources available. Useful when resources got completelly exausted and the model refuses to solve.

Definition at line 222 of file Opt.h.

int previousYear

protected

Definition at line 217 of file Opt.h.

vector<string> priPr

protected

Definition at line 193 of file Opt.h.

vector< vector <int> > priPrCombs

protected

A vector with all the possible combinations of primary products.

Definition at line 198 of file Opt.h.

unsigned int* rind_g

protected

Definition at line 252 of file Opt.h.

unsigned int* rind_L

protected

Definition at line 255 of file Opt.h.

unsigned int* rind_L_total

protected

Definition at line 257 of file Opt.h.

int secondYear

protected

Definition at line 219 of file Opt.h.

vector<string> secPr

protected

Definition at line 194 of file Opt.h.

map<string, endvar> vars

protected

List of variables in the model and their domain: pr product, sec prod, all products or all products over each subregion pair (exports)

Definition at line 203 of file Opt.h.

map<string, vector < vector < vector < vector <int> > > > > vpositions

protected

cached position in the concatenated vector for each variables. Dimensions are l1reg, l2reg, prod, (l2To region).

Definition at line 204 of file Opt.h.

int worldCodeLev2

protected

Definition at line 220 of file Opt.h.

double* x_lam

protected

Definition at line 248 of file Opt.h.

---

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

• /srv/datadisk03/home/ffsm/code/ffsm_pp/src/Opt.h
• /srv/datadisk03/home/ffsm/code/ffsm_pp/src/Opt.cpp