#include <Opt.h>

Inheritance diagram for Opt:

[legend]

Collaboration diagram for Opt:

[legend]

Public Member Functions
	Opt (ThreadManager *MTHREAD_h)
	Constructor.

	~Opt ()

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)

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.

void	msgOut (const int &msgCode_h, const int &msg_h, const bool &refreshGUI_h=true) const
	Overloaded function to print the output log.

void	msgOut (const int &msgCode_h, const double &msg_h, const bool &refreshGUI_h=true) const
	Overloaded function to print the output log.

int	s2i (const string &string_h) const
	string to integer conversion

double	s2d (const string &string_h) const
	string to double conversion

double	s2d (const string &string_h, const bool &replaceComma) const
	string to double conversion

bool	s2b (const string &string_h) const
	string to bool conversion

string	i2s (const int &int_h) const
	integer to string conversion

string	d2s (const double &double_h) const
	double to string conversion

string	b2s (const bool &bool_h) const
	bool to string conversion

vector< int >	s2i (const vector< string > &string_h) const
	string to integer conversion (vector)

vector< double >	s2d (const vector< string > &string_h, const bool &replaceComma=false) const
	string to double conversion (vector)

vector< bool >	s2b (const vector< string > &string_h) const
	string to bool conversion (vector)

vector< string >	i2s (const vector< int > &int_h) const
	integer to string conversion (vector)

vector< string >	d2s (const vector< double > &double_h) const
	double to string conversion (vector)

vector< string >	b2s (const vector< bool > &bool_h) const
	bool to string conversion (vector)

int	getType (const string &type_h) const
	Return a type according to enum TYPE_* from a string (eg: "string" -> TYPE_STRING (2))

void	refreshGUI () const
	Ping to periodically return the control to the GUI.

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)

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.

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.

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.

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.

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.

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.

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)

template<typename V >
void	debugMap (const map< iisskey, V > &mymap)
	Debug a map.

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)

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)

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)

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)

template<typename K >
double	getAvg (const vector< K > &v)
	Returns the average of the elements in the vector.

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

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.

template<typename T >
std::string	toString (const T &x) const

Solution Methods
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.

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.

map< string, int >	initPos
	A map that returns the initial index position in the concatenated array for each variable.

map< int, string >	initPos_rev
	A map with the name of the variable keyed by its initial position in the index.

vector< int >	cInitPos
	A vector that returns the initial index position in the concatenated array for each constrain.

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)

map< string, vector< vector< vector< vector< int > > > > >	vpositions
	cached position in the concatenated vector for each variables. Dimensions are l1reg, l2reg, prod, (l2To region).

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

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.

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)

vector< vector< Index > >	nzhelements
	nzero elements for the hessian matrix

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)

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

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.

void	declareConstrains ()
	declare the constrains, their domain, their direction and their associated evaluation function

void	cacheInitialPosition ()
	cache the initial positions of the variables and the constrains

void	calculateNumberVariablesConstrains ()
	calculate the number of variables and constrains

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

int	getDomainElements (int domain)
	return the number of elements of a domain

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

int	getVarInstances (const string &varName)
	build the matrix of the positions for a given variable or contrain

void	calculateSparsityPatternJ ()

void	calculateSparsityPatternH ()

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

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

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.

const int	gip (const int &cn) const
	Return the initial index position of a certain constrain.

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.

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.

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.

const int	gdt (const string &varName)
	Get the domain type of a given variable.

const int	gdt (const int &cn)
	Get the domain type of a given constrain.

int	getConstrainDirectionByIndex (int idx)
	Return the direction of a given constrain.

double	getBoundByIndex (const int &bound_type, const int &idx)
	Return the bound of a given variable (by index)

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().

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.

int	getConNumber (constrain *con)
	Return the position in the cons vector.

void	copyInventoryResourses ()
	Copy the inventoried resources in the in vector for better performances.

void	tempDebug ()

void	debugPrintParameters ()

Methods to block default compiler methods.
The compiler automatically generates the following three methods. Since the default compiler implementation is generally not what you want (for all but the most simple classes), we usually put the declarations of these methods in the private section and never implement them. This prevents the compiler from implementing an incorrect "default" behavior without us knowing. (See Scott Meyers book, "Effective C++")
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 &)

Additional Inherited Members
Protected Attributes inherited from BaseClass
ThreadManager *	MTHREAD
	Pointer to the Thread manager.

Detailed Description

Definition at line 52 of file Opt.h.

Constructor & Destructor Documentation

◆ Opt() [1/2]

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()

~Opt ( )

Definition at line 545 of file Opt.cpp.

         {
 
}

◆ Opt() [2/2]

Opt ( const Opt & )

protected

Member Function Documentation

◆ app()

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

inlineprotected

Definition at line 172 of file Opt.h.

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

ModelData::assessProdPossibility

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

A simple function to assess if a specified product can be made by a certain forest type and diameter ...

Definition ModelData.cpp:441

ThreadManager::MD

ModelData * MD

the model data object

Definition ThreadManager.h:72

Here is the call graph for this function:

◆ buildPositionVector()

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

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;
}

Referenced by cachePositions().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cacheInitialPosition()

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);
  }
}

Referenced by get_nlp_info().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cachePositions()

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));
  }
 
}

Referenced by get_nlp_info().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ calculateNumberVariablesConstrains()

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)");
}

Referenced by get_nlp_info().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ calculateSparsityPatternH()

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);
      }
    }
  }
}

◆ calculateSparsityPatternJ()

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);
    }
  }
}

◆ copyInventoryResourses()

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;
}

Referenced by get_nlp_info().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ debugPrintParameters()

void debugPrintParameters ( )

protected

◆ declareConstrains()

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);
;
 
 
 
}

Referenced by get_nlp_info().

Here is the caller graph for this function:

◆ declareVariable()

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));
}

Referenced by declareVariables().

Here is the caller graph for this function:

◆ declareVariables()

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");
}

Referenced by get_nlp_info().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ eval_constraints()

template<class T >

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;
}

Referenced by eval_g(), and generate_tapes().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ eval_f()

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;
}

Here is the call graph for this function:

◆ eval_g()

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;
}

Here is the call graph for this function:

◆ eval_grad_f()

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;
}

◆ eval_h()

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;
}

◆ eval_jac_g()

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;
}

◆ eval_obj()

template<class T >

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
 
}

Referenced by eval_f(), and generate_tapes().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ finalize_solution()

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);
 
}

Here is the call graph for this function:

◆ gdt() [1/2]

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;
}

◆ gdt() [2/2]

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;
  }
}

Referenced by getVarInstances(), and gix_uncached().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ generate_tapes()

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

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;
 
}

Referenced by get_nlp_info().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ get_bounds_info()

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;
}

Here is the call graph for this function:

◆ get_nlp_info()

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:

◆ get_starting_point()

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;
}

Referenced by generate_tapes().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getBoundByIndex()

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.;
}

Referenced by get_bounds_info().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getConNumber()

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.");
}

Here is the call graph for this function:

◆ getConstrainByIndex()

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!");
}

Here is the call graph for this function:

◆ getConstrainDirectionByIndex()

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!");
}

Referenced by get_bounds_info().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getDetailedBoundByVarAndIndex()

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]);
    }
  }
}

Referenced by getBoundByIndex().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getDomainElements()

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)+")");
  }
}

Referenced by cacheInitialPosition(), calculateNumberVariablesConstrains(), and getVarInstances().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getVarInstances()

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));
}

Here is the call graph for this function:

◆ gfd()

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.

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

ModelData::getForData

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

Definition ModelData.cpp:1281

Here is the call graph for this function:

◆ gip() [1/2]

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);
}

◆ gip() [2/2]

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;
  }
}

Referenced by getConstrainByIndex(), getConstrainDirectionByIndex(), getDetailedBoundByVarAndIndex(), and gix_uncached().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ gix() [1/2]

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];
}

◆ gix() [2/2]

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;
  }
}

Referenced by eval_constraints(), eval_obj(), finalize_solution(), and get_starting_point().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ gix_uncached()

template<class T >

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;
  }
}

Referenced by buildPositionVector().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ gpd()

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.

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

ModelData::getProdData

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

Definition ModelData.cpp:1216

Referenced by eval_constraints(), eval_obj(), get_starting_point(), and getDetailedBoundByVarAndIndex().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ intermediate_callback()

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;
}

Here is the call graph for this function:

◆ mymax() [1/2]

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

Definition at line 1753 of file Opt.cpp.

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

◆ mymax() [2/2]

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

Definition at line 1749 of file Opt.cpp.

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

Referenced by eval_obj().

Here is the caller graph for this function:

◆ operator=()

Opt & operator= ( const Opt & )

protected

◆ sfd()

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.

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

ModelData::setForData

void setForData(const double &value_h, const string &type_h, const int &regId_h, const string &forType_h, const string &freeDim_h, const int &year=DATA_NOW, const bool &allowCreate=false)

Definition ModelData.cpp:1412

Here is the call graph for this function:

◆ spd()

void spd	(	const double &	value_h,
		const string &	type_h,
		const int &	regId_h,
		const string &	prodId_h,
		const int &	year = `DATA_NOW`,
		const bool &	allowCreate = `false`,
		const string &	freeDim_h = `""`
	)		const

inlineprotected

Definition at line 170 of file Opt.h.

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

ModelData::setProdData

void setProdData(const double &value_h, const string &type_h, const int &regId_h, const string &prodId_h, const int &year=DATA_NOW, const bool &allowCreate=false, const string &freeDim_h="")

Definition ModelData.cpp:1352

Referenced by finalize_solution().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ tempDebug()

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);
 
}

◆ unpack()

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");
}

Referenced by getDetailedBoundByVarAndIndex().

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

◆ allPr

vector<string> allPr

protected

Definition at line 195 of file Opt.h.

Referenced by buildPositionVector(), eval_constraints(), eval_obj(), finalize_solution(), get_nlp_info(), get_starting_point(), and getDetailedBoundByVarAndIndex().

◆ cind_g

unsigned int* cind_g

protected

Definition at line 253 of file Opt.h.

Referenced by eval_jac_g(), finalize_solution(), and generate_tapes().

◆ cind_L

unsigned int* cind_L

protected

Definition at line 256 of file Opt.h.

Referenced by eval_h(), finalize_solution(), and generate_tapes().

◆ cind_L_total

unsigned int* cind_L_total

protected

Definition at line 258 of file Opt.h.

Referenced by eval_h(), finalize_solution(), and generate_tapes().

◆ cInitPos

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.

Referenced by cacheInitialPosition(), and gip().

◆ cons

vector<constrain> cons

protected

Definition at line 225 of file Opt.h.

Referenced by cacheInitialPosition(), cachePositions(), calculateNumberVariablesConstrains(), declareConstrains(), gdt(), getConNumber(), getConstrainByIndex(), and getConstrainDirectionByIndex().

◆ cpositions

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.

Referenced by cachePositions(), and gix().

◆ debugRunOnce

bool debugRunOnce

protected

Definition at line 221 of file Opt.h.

Referenced by Opt().

◆ firstYear

int firstYear

protected

Definition at line 218 of file Opt.h.

Referenced by get_nlp_info().

◆ hessval

double* hessval

protected

Definition at line 259 of file Opt.h.

Referenced by eval_h(), finalize_solution(), and generate_tapes().

◆ HP_t

unsigned int** HP_t

protected

Definition at line 251 of file Opt.h.

Referenced by finalize_solution(), and generate_tapes().

◆ initOpt

bool initOpt

protected

Definition at line 224 of file Opt.h.

Referenced by get_nlp_info(), and Opt().

◆ initPos

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.

Referenced by cacheInitialPosition(), and gip().

◆ initPos_rev

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.

Referenced by cacheInitialPosition(), and getBoundByIndex().

◆ ins

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.

Referenced by copyInventoryResourses(), and eval_constraints().

◆ jacval

double* jacval

protected

Definition at line 254 of file Opt.h.

Referenced by eval_jac_g(), finalize_solution(), and generate_tapes().

◆ l2r

vector< vector <int> > l2r

protected

Definition at line 197 of file Opt.h.

Referenced by buildPositionVector(), copyInventoryResourses(), eval_constraints(), eval_obj(), finalize_solution(), get_nlp_info(), get_starting_point(), getDetailedBoundByVarAndIndex(), getDomainElements(), gix_uncached(), and unpack().

◆ nAllPr

int nAllPr

protected

Definition at line 210 of file Opt.h.

Referenced by get_nlp_info(), getDomainElements(), gix_uncached(), and unpack().

◆ nCons

int nCons

protected

Definition at line 213 of file Opt.h.

Referenced by calculateNumberVariablesConstrains(), calculateSparsityPatternH(), calculateSparsityPatternJ(), get_nlp_info(), getConstrainByIndex(), getConstrainDirectionByIndex(), Opt(), and tempDebug().

◆ nEqualityConstrains

int nEqualityConstrains

protected

Definition at line 214 of file Opt.h.

Referenced by calculateNumberVariablesConstrains().

◆ nGreaterEqualZeroConstrains

int nGreaterEqualZeroConstrains

protected

Definition at line 216 of file Opt.h.

Referenced by calculateNumberVariablesConstrains().

◆ nL2r

int nL2r

protected

Definition at line 211 of file Opt.h.

Referenced by get_nlp_info(), and getDomainElements().

◆ nLowerEqualZeroConstrains

int nLowerEqualZeroConstrains

protected

Definition at line 215 of file Opt.h.

Referenced by calculateNumberVariablesConstrains().

◆ nnz_jac

int nnz_jac

protected

Definition at line 260 of file Opt.h.

Referenced by eval_jac_g(), and generate_tapes().

◆ nnz_L

int nnz_L

protected

Definition at line 261 of file Opt.h.

Referenced by eval_h(), and generate_tapes().

◆ nnz_L_total

int nnz_L_total

protected

Definition at line 261 of file Opt.h.

Referenced by eval_h(), and generate_tapes().

◆ nOthPr

int nOthPr

protected

Definition at line 207 of file Opt.h.

Referenced by get_nlp_info().

◆ nPriPr

int nPriPr

protected

Definition at line 206 of file Opt.h.

Referenced by eval_constraints(), eval_obj(), get_nlp_info(), getDomainElements(), gix_uncached(), and unpack().

◆ nPriPrCombs

int nPriPrCombs

protected

Definition at line 208 of file Opt.h.

Referenced by get_nlp_info(), getDomainElements(), and gix_uncached().

◆ nSecPr

int nSecPr

protected

Definition at line 209 of file Opt.h.

Referenced by get_nlp_info(), getDomainElements(), gix_uncached(), and unpack().

◆ nVar

int nVar

protected

Definition at line 212 of file Opt.h.

Referenced by calculateNumberVariablesConstrains(), calculateSparsityPatternH(), calculateSparsityPatternJ(), get_nlp_info(), Opt(), and tempDebug().

◆ nzhelements

vector<vector <Index> > nzhelements

protected

nzero elements for the hessian matrix

Definition at line 227 of file Opt.h.

Referenced by calculateSparsityPatternH(), and tempDebug().

◆ nzjelements

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.

Referenced by calculateSparsityPatternJ().

◆ options_g

int options_g[4]

protected

Definition at line 262 of file Opt.h.

Referenced by eval_jac_g(), and generate_tapes().

◆ options_L

int options_L[4]

protected

Definition at line 263 of file Opt.h.

Referenced by eval_h(), and generate_tapes().

◆ othPr

vector<string> othPr

protected

Definition at line 196 of file Opt.h.

Referenced by eval_constraints(), and get_nlp_info().

◆ overharvestingAllowance

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.

Referenced by eval_constraints(), and get_nlp_info().

◆ previousYear

int previousYear

protected

Definition at line 217 of file Opt.h.

Referenced by eval_constraints(), get_nlp_info(), and get_starting_point().

◆ priPr

vector<string> priPr

protected

Definition at line 193 of file Opt.h.

Referenced by buildPositionVector(), eval_constraints(), eval_obj(), finalize_solution(), get_nlp_info(), and get_starting_point().

◆ priPrCombs

vector< vector <int> > priPrCombs

protected

A vector with all the possible combinations of primary products.

Definition at line 198 of file Opt.h.

Referenced by buildPositionVector(), copyInventoryResourses(), eval_constraints(), and get_nlp_info().

◆ rind_g

unsigned int* rind_g

protected

Definition at line 252 of file Opt.h.

Referenced by eval_jac_g(), finalize_solution(), and generate_tapes().

◆ rind_L

unsigned int* rind_L

protected

Definition at line 255 of file Opt.h.

Referenced by eval_h(), finalize_solution(), and generate_tapes().

◆ rind_L_total

unsigned int* rind_L_total

protected

Definition at line 257 of file Opt.h.

Referenced by eval_h(), finalize_solution(), and generate_tapes().

◆ secondYear

int secondYear

protected

Definition at line 219 of file Opt.h.

Referenced by eval_obj(), and get_nlp_info().

◆ secPr

vector<string> secPr

protected

Definition at line 194 of file Opt.h.

Referenced by buildPositionVector(), eval_constraints(), eval_obj(), finalize_solution(), get_nlp_info(), and get_starting_point().

◆ vars

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.

Referenced by cacheInitialPosition(), cachePositions(), calculateNumberVariablesConstrains(), declareVariable(), finalize_solution(), gdt(), get_starting_point(), and getBoundByIndex().

◆ vpositions

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.

Referenced by cachePositions(), and gix().

◆ worldCodeLev2

int worldCodeLev2

protected

Definition at line 220 of file Opt.h.

Referenced by eval_constraints(), and get_nlp_info().

◆ x_lam

double* x_lam

protected

Definition at line 248 of file Opt.h.

Referenced by eval_h(), finalize_solution(), and generate_tapes().

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

/exports/ffsm/src_latest/src/Opt.h
/exports/ffsm/src_latest/src/Opt.cpp

Public Member Functions

Solution Methods

Methods to block default compiler methods.

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ Opt() [1/2]

◆ ~Opt()

◆ Opt() [2/2]

Member Function Documentation

◆ app()

◆ buildPositionVector()

◆ cacheInitialPosition()

◆ cachePositions()

◆ calculateNumberVariablesConstrains()

◆ calculateSparsityPatternH()

◆ calculateSparsityPatternJ()

◆ copyInventoryResourses()

◆ debugPrintParameters()

◆ declareConstrains()

◆ declareVariable()

◆ declareVariables()

◆ eval_constraints()

◆ eval_f()

◆ eval_g()

◆ eval_grad_f()

◆ eval_h()

◆ eval_jac_g()

◆ eval_obj()

◆ finalize_solution()

◆ gdt() [1/2]

◆ gdt() [2/2]

◆ generate_tapes()

◆ get_bounds_info()

◆ get_nlp_info()

◆ get_starting_point()

◆ getBoundByIndex()

◆ getConNumber()

◆ getConstrainByIndex()

◆ getConstrainDirectionByIndex()

◆ getDetailedBoundByVarAndIndex()

◆ getDomainElements()

◆ getVarInstances()

◆ gfd()

◆ gip() [1/2]

◆ gip() [2/2]

◆ gix() [1/2]

◆ gix() [2/2]

◆ gix_uncached()

◆ gpd()

◆ intermediate_callback()

◆ mymax() [1/2]

◆ mymax() [2/2]

◆ operator=()

◆ sfd()

◆ spd()

◆ tempDebug()

◆ unpack()

Member Data Documentation

◆ allPr

◆ cind_g

◆ cind_L

◆ cind_L_total

◆ cInitPos

◆ cons

◆ cpositions

◆ debugRunOnce

◆ firstYear

◆ hessval

◆ HP_t

◆ initOpt

◆ initPos

◆ initPos_rev

◆ ins

◆ jacval

◆ l2r

◆ nAllPr

◆ nCons

◆ nEqualityConstrains

◆ nGreaterEqualZeroConstrains