Opt.h

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2015 by Laboratoire d'Economie Forestière               *
 *   http://ffsm-project.org                                               *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 3 of the License, or     *
 *   (at your option) any later version, given the compliance with the     *
 *   exceptions listed in the file COPYING that is distribued together     *
 *   with this file.                                                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#ifndef STDOPT_H
#define STDOPT_H


#include "IpTNLP.hpp"
#include <adolc.h>
#include <adolc_sparse.h>

//regmas headers
#include "BaseClass.h"
#include "ThreadManager.h"
#include "ModelData.h"

#define tag_f 1
#define tag_g 2
#define tag_L 3
#define HPOFF 30 //original: 30

/// Class containing the optimization problem (the matrix and its methods)

/**

@author Antonello Lobianco
*/

using namespace Ipopt;

struct constrain;
struct endvar;

class Opt: public BaseClass, public TNLP{
public:
            Opt(ThreadManager* MTHREAD_h); ///< Constructor
             ~Opt();
  /**@name Overloaded from TNLP */
  //@{
  /** Method to return some info about the nlp */
  virtual bool        get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
                   Index& nnz_h_lag, IndexStyleEnum& index_style);
  /** Method to return the bounds for my problem */
  virtual bool        get_bounds_info(Index n, Number* x_l, Number* x_u,
                    Index m, Number* g_l, Number* g_u);

  /** Method to return the starting point for the algorithm */
  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 to return the objective value */
  template<class T> bool eval_obj(Index n, const T *x, T& obj_value);

  /** Template to compute contraints */
  template<class T> bool eval_constraints(Index n, const T *x, Index m, T *g);

  /** Original method from Ipopt to return the objective value */
  /** remains unchanged */
  virtual bool        eval_f(Index n, const Number* x, bool new_x, Number& obj_value);

  /** Original method from Ipopt to return the gradient of the objective */
  /** remains unchanged */
  virtual bool        eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);

  /**  Original method from Ipopt to return the constraint residuals */
  /** remains unchanged */
  virtual bool        eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);

  /** 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 */
  virtual bool        eval_jac_g(Index n, const Number* x, bool new_x,
                   Index m, Index nele_jac, Index* iRow, Index *jCol,
                   Number* values);

  /** 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 */
  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);

  //@}

  /** @name Solution Methods */
  //@{
  /** This method is called when the algorithm is complete so the TNLP can store/write the solution */
  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);
  //@}

    /** Return information on each iteration */
    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);

  //***************    start ADOL-C part ***********************************

  /** Method to generate the required tapes */
  virtual void generate_tapes(Index n, Index m, Index& nnz_jac_g, Index& nnz_h_lag);

  //***************    end   ADOL-C part ***********************************

  // **************   start FFSM part ************************************
  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); ///< return the number of instances of a variable, given his domain type
  ///< 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);

  //template <class T> const T& mymax ( const T& a, const T& b );

  // **************   end FFSM part **************************************


protected:

  // convenient handles to equivalent ModelData functions..
  const double        gpd(const string &type_h, const int& regId_h, const string &prodId_h, const int& year=DATA_NOW, const string &freeDim_h="") const {return MTHREAD->MD->getProdData(type_h, regId_h, prodId_h, year, freeDim_h);};
  const double        gfd(const string &type_h, const int& regId_h, const string &forType_h, const string &diamClass_h, const int& year=DATA_NOW) const {return MTHREAD->MD->getForData(type_h, regId_h, forType_h, diamClass_h, year);};
  void                spd(const double& value_h, const string &type_h, const int& regId_h, const string &prodId_h, const int& year=DATA_NOW, const bool& allowCreate=false, const string &freeDim_h="") const {MTHREAD->MD->setProdData(value_h, type_h, regId_h, prodId_h, year, allowCreate, freeDim_h);};
  void                sfd(const double& value_h, const string &type_h, const int& regId_h, const string &forType_h, const string &diamClass_h, const int& year=DATA_NOW, const bool& allowCreate=false) const {MTHREAD->MD->setForData(value_h, type_h, regId_h, forType_h, diamClass_h, year, allowCreate);};
  bool                app(const string &prod_h, const string &forType_h, const string &dClass_h) const {return MTHREAD->MD->assessProdPossibility(prod_h, forType_h, dClass_h);};
  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
  //vector < vector <int> > createCombinationsVector(const int& nItems); ///< Return a vector containing any possible combination of nItems items (including any possible subset). The returned vector has in each slot the items present in that specific combination.
  void                copyInventoryResourses(); ///< Copy the inventoried resources in the in vector for better performances

  void tempDebug();

  //virtual void        eval_obj (Index n, const T *x, T& obj_value);

  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.
  void                    debugPrintParameters();
  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


  /**@name 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++")
   *
  */
  //@{
  //  MyADOLC_NLP();
  Opt(const Opt&);
  Opt& operator=(const Opt&);
  //@}

  //@{

  double *x_lam;

  //** variables for sparsity exploitation
  unsigned int **HP_t;         /* compressed block row storage */
  unsigned int *rind_g;        /* row indices    */
  unsigned int *cind_g;        /* column indices */
  double *jacval;              /* values         */
  unsigned int *rind_L;        /* row indices    */
  unsigned int *cind_L;        /* column indices */
  unsigned int *rind_L_total;  /* row indices    */
  unsigned int *cind_L_total;  /* column indices */
  double *hessval;             /* values */
  int nnz_jac;
  int nnz_L, nnz_L_total;
  int options_g[4];
  int options_L[4];


  //@}

};

struct constrain{
                        constrain(){comment="";};
  string                                    name;
  string                                 comment;
  int                                     domain;
  int                                  direction;

};

struct endvar{
  string                                    name;
  int                                     domain;
  string                                    desc; ///< Description of the variable
  double                                 l_bound; ///< A fixed numerical lower bound for all the domain
  double                                 u_bound; ///< A fixed numerical upper bound for all the domain
  string                             l_bound_var; ///< A variable giving the lower bound. If present, the value defined in the variable overrides l_bound.
  string                             u_bound_var; ///< A variable giving the upper bound. If present, the value defined in the variable overrides u_bound.
};



#endif