main.cc 5.79 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
/*===========================================================================*\
 *                                                                           *
 *                               CoMISo                                      *
 *      Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen      *
 *                           www.rwth-graphics.de                            *
 *                                                                           *
 *---------------------------------------------------------------------------* 
 *  This file is part of CoMISo.                                             *
 *                                                                           *
 *  CoMISo 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.                                      *
 *                                                                           *
 *  CoMISo 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 CoMISo.  If not, see <http://www.gnu.org/licenses/>.          *
 *                                                                           *
\*===========================================================================*/ 

#include <CoMISo/Config/config.hh>

#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <CoMISo/NSolver/NPDerivativeChecker.hh>
#include <CoMISo/NSolver/GUROBISolver.hh>
#include <CoMISo/NSolver/CPLEXSolver.hh>
33
#include <CoMISo/NSolver/COMISOSolver.hh>
34
#include <CoMISo/NSolver/LinearConstraint.hh>
35
#include <CoMISo/NSolver/VariableType.hh>
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50


// generate an instance of a nonlinear problem by deriving from base class NProblemInterface
// implement all virtual functions in order to solve this problem by any of the solvers located
// in CoMISo/NSolver

class SmallNProblem : public COMISO::NProblemInterface
{
public:

  // Sparse Matrix Type
  //  typedef Eigen::DynamicSparseMatrix<double,Eigen::ColMajor> SMatrixNP;


  // specify a function which has several local minima
David Bommes's avatar
David Bommes committed
51
  // f(x,y)=(x-2y+1)^2 + (x-5)^2
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69

  // number of unknown variables, here x and y = 2
  virtual int    n_unknowns   (                                )
  {
    return 2;
  }

  // initial value where the optimization should start from
  virtual void   initial_x    (       double* _x               )
  {
    _x[0] = 0.0;
    _x[1] = 0.0;
  }

  // function evaluation at location _x
  virtual double eval_f       ( const double* _x               )
  {
    double term  = _x[0] - 2.0*_x[1] + 1.0;
David Bommes's avatar
David Bommes committed
70
    double term2 = _x[0] - 5.0;
71 72 73 74 75 76 77 78

    return term*term + term2*term2;
  }

  // gradient evaluation at location _x
  virtual void   eval_gradient( const double* _x, double*    _g)
  {
    double term  = _x[0] - 2.0*_x[1] + 1.0;
David Bommes's avatar
David Bommes committed
79
    double term2 = _x[0] - 5.0;
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104

    _g[0] =  2.0*term + 2.0*term2;
    _g[1] = -4.0*term;
   }

  // hessian matrix evaluation at location _x
  virtual void   eval_hessian ( const double* _x, SMatrixNP& _H)
  {
    _H.resize(n_unknowns(), n_unknowns());
    _H.setZero();

    _H.coeffRef(0,0) =  4.0;
    _H.coeffRef(1,0) = -4.0;
    _H.coeffRef(0,1) = -4.0;
    _H.coeffRef(1,1) =  8.0;
  }

  // print result
  virtual void   store_result ( const double* _x               )
  {
    std::cerr << "Energy: " << eval_f(_x) << std::endl;
    std::cerr << "(x,y) = (" << _x[0] << "," << _x[1] << ")" << std::endl;
  }

  // advanced properties
105
  virtual bool   constant_hessian() const { return true; }
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
};


//------------------------------------------------------------------------------------------------------

// Example main
int main(void)
{
  std::cout << "---------- 1) Get an instance of a NProblem..." << std::endl;
  SmallNProblem snp;

  std::cout << "---------- 2) (optional for debugging) Check derivatives of problem..." << std::endl;
  COMISO::NPDerivativeChecker npd;
  npd.check_all(&snp);

  std::cout << "---------- 3) setup list of integer variables..." << std::endl;
  std::vector<COMISO::PairIndexVtype> discrete_variables;
  discrete_variables.push_back( COMISO::PairIndexVtype(0,COMISO::Integer) );

  std::cout << "---------- 4) setup constraints..." << std::endl;
  std::vector<COMISO::NConstraintInterface*> constraints;
David Bommes's avatar
David Bommes committed
127
  // setup constraint x+y <= 6.5
128 129 130
  COMISO::LinearConstraint::SVectorNC coeffs(2);
  coeffs.coeffRef(0) = 1.0;
  coeffs.coeffRef(1) = 1.0;
David Bommes's avatar
David Bommes committed
131
  COMISO::LinearConstraint lc(coeffs, -6.5, COMISO::LinearConstraint::NC_LESS_EQUAL);
132
  constraints.push_back(&lc);
133 134 135 136 137 138 139


// check if IPOPT solver available in current configuration
#if( COMISO_GUROBI_AVAILABLE)
  std::cout << "---------- 5) Get GUROBI solver... " << std::endl;
  COMISO::GUROBISolver gsol;

140
  std::cout << "---------- 6) Solve..." << std::endl;
141 142 143 144 145 146

  gsol.solve(&snp, constraints, discrete_variables);
#endif

  // check if TAO solver available in current configuration
#if( COMISO_CPLEX_AVAILABLE)
147
  std::cout << "---------- 7) Solve with CPLEX solver... " << std::endl;
148 149
  COMISO::CPLEXSolver csol;

150
  std::cout << "---------- 8) Solve..." << std::endl;
151 152 153 154

  csol.solve(&snp, constraints, discrete_variables);
#endif

155 156 157 158 159 160 161
  std::cout << "---------- 9) Approximate with COMISO solver... " << std::endl;
  COMISO::COMISOSolver cosol;

  std::cout << "---------- 10) Solve..." << std::endl;

  cosol.solve(&snp, constraints, discrete_variables);

162 163 164
  return 0;
}