main.cc 5.63 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
/*===========================================================================*\
 *                                                                           *
 *                               CoMISo                                      *
 *      Copyright (C) 2008-2019 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 <iostream>

#include <CoMISo/NSolver/NewtonSolver.hh>
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <vector>

32
#include <CoMISo/Utils/ExactConstraintSatisfaction.hh>
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
//------------------------------------------------------------------------------------------------------

class SmallNProblem : public COMISO::NProblemInterface
{
public:
  // specify a function which has several local minima
  // f(x,y) = x^4 + y^4

  // 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] = 4.0;
    _x[1] = 2.0;
  }

  // function evaluation at location _x
  virtual double eval_f       ( const double* _x               )
  {
    return std::pow(_x[0], 4) + std::pow(_x[1], 4);
  }

  // gradient evaluation at location _x
  virtual void   eval_gradient( const double* _x, double*    _g)
  {
    _g[0] =  4.0*std::pow(_x[0], 3);
    _g[1] =  4.0*std::pow(_x[1], 3);
   }

  // 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) =  12.0*std::pow(_x[0], 2);
    _H.coeffRef(1,0) =  0.0;
    _H.coeffRef(0,1) =  0.0;
    _H.coeffRef(1,1) =  12.0*std::pow(_x[1], 2);
  }

  // print result
  virtual void   store_result ( const double* _x               )
  {
    solution.resize(n_unknowns());
    for (int i = 0; i < n_unknowns(); ++i)
      solution[i] = _x[i];
  }

  // advanced properties
Robin Brost's avatar
Robin Brost committed
88
  virtual bool   constant_hessian() const { return false; }
89
90
91
92
93
94
95

  std::vector<double> solution;
};

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

99
  std::cout << "---------- 2) Set up constraints..." << std::endl;
100

101
102
103
104
105

  int n_constraints = 2; // there will be two constraints
  Eigen::VectorXi b;
  COMISO::ExactConstraintSatisfaction::SP_Matrix_R A(n_constraints, problem.n_unknowns());
  b.setZero(n_constraints);
106

Robin Brost's avatar
Robin Brost committed
107
  // first constraint: first variable equals three times second
Robin Brost's avatar
Robin Brost committed
108
  //different number of constraints :
109
110
  if(n_constraints == 2)
  {
111
    A.coeffRef(0,0) =  2;
112
    A.coeffRef(0,1) =  -1;
113
114
    b.coeffRef(0)   =  0;
    A.coeffRef(1,0) =  -2;
115
    A.coeffRef(1,1) =  8;
116
    b.coeffRef(1)   =  21;
Robin Brost's avatar
Robin Brost committed
117
  }
118

119
  std::cout << "Constraints: Ax = b with A = \n" << A << "and b = \n" << b << std::endl;
120
121
122

  std::cout << "---------- 3) Solve with Newton Solver..." << std::endl;
  COMISO::NewtonSolver nsolver;
123
  nsolver.set_verbosity(15);
124
125
126
  Eigen::SparseMatrix<double> Ad = A.cast<double>();
  Eigen::VectorXd bd = b.cast<double>();
  nsolver.solve(&problem, Ad, bd);
127
128
129
130
131
132
133
134
135
136
137
138
139

  std::cout << "---------- 4) Print solution..." << std::endl;
  std::cout << std::setprecision(100);
  for (unsigned int i = 0; i < problem.n_unknowns(); ++i)
    std::cout << "x[" << i << "] = " << problem.solution[i] << std::endl;


  std::cout << "---------- 5) Check constraint violation..." << std::endl;
  Eigen::VectorXd x;
  x.resize(problem.n_unknowns());
  for (unsigned int i = 0; i < problem.n_unknowns(); ++i)
    x.coeffRef(i) = problem.solution[i];

140
  std::cout << "Constraint violation: " << (A.cast<double>() *x - b.cast<double>()).squaredNorm() << std::endl;
141

142
  std::cout << "---------- 6) Try to exactly fulfill constraints..." << std::endl;
143

144
  COMISO::ExactConstraintSatisfaction satisfy;
Robin Brost's avatar
Robin Brost committed
145
//  satisfy.print_matrix(A);
146
  satisfy.evaluation(A, b, x);
147
148


149
150
  for (unsigned int i = 0; i < problem.n_unknowns(); ++i)
    std::cout << "x[" << i << "] = " << x[i] << std::endl;
151

152
  std::cout << "---------- 7) Check constraint violation again..." << std::endl;
153

154
  std::cout << "Constraint violation: " << (A.cast<double>() *x - b.cast<double>()).squaredNorm() << std::endl;
155
156
157
158


  return 0;
}