1 /****************************************************************************/
2 /* This file is part of FreeFEM. */
3 /* */
4 /* FreeFEM is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU Lesser General Public License as */
6 /* published by the Free Software Foundation, either version 3 of */
7 /* the License, or (at your option) any later version. */
8 /* */
9 /* FreeFEM is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU Lesser General Public License for more details. */
13 /* */
14 /* You should have received a copy of the GNU Lesser General Public License */
15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>. */
16 /****************************************************************************/
17 // SUMMARY : ...
18 // LICENSE : LGPLv3
19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
20 // AUTHORS : Sylvain Auliac
21 // E-MAIL : auliac@ann.jussieu.fr
22
23 /* clang-format off */
24 //ff-c++-LIBRARY-dep: ipopt mumps_seq blas libseq fc
25 //ff-c++-cpp-dep:
26 /* clang-format on */
27
28 // using namespace std;
29
30 // TODO: remove this block as soon as autoconf is removed from FreeFem++
31 #ifndef CMAKE
32 #include "../../config.h"
33 #endif
34
35 #include "coin/IpTNLP.hpp"
36 #include "coin/IpIpoptApplication.hpp"
37 #include "ff++.hpp"
38
39 extern Block *currentblock;
40
41 typedef double R;
42 typedef KN_< R > Rn_;
43 typedef KN< R > Rn;
44 typedef KNM_< R > Rnm_;
45 typedef KNM< R > Rnm;
46
47 /*****************************************************************************************************************************
48 * Some misc. function usefull later...
49 *****************************************************************************************************************************/
50
51 // A variadic function to add an undefinite number of elements to a set of short int
52 // This is used to define the set of named parameter which are not used when some assumptions
53 // upon the optimization poblem functions are met
AddElements(std::set<unsigned short> & _set,int amount,int first,...)54 void AddElements(std::set< unsigned short > &_set, int amount, int first, ...) {
55 int elem = 0;
56 va_list vl;
57
58 va_start(vl, first);
59 _set.insert(first);
60
61 for (int i = 1; i < amount; i++) {
62 elem = va_arg(vl, int);
63 _set.insert(elem);
64 }
65
66 va_end(vl);
67 }
68
69 // A raw pointer cleaner
70 template< class T >
clean(T * p)71 inline void clean(T *p) {
72 if (p) {
73 delete p;
74 p = 0;
75 }
76 }
77
78 // Pair compare (certainly already implemented in the STL with KeyLess...)
operator <=(const std::pair<int,int> & l,const std::pair<int,int> & r)79 inline bool operator<=(const std::pair< int, int > &l, const std::pair< int, int > &r) {
80 return (l.first < r.first) || (l.first == r.first && l.second <= r.second);
81 }
82
83 // Some logical operators (exclussive or and its negation)
XOR(bool a,bool b)84 inline bool XOR(bool a, bool b) { return (!a && b) || (a && !b); }
85
NXOR(bool a,bool b)86 inline bool NXOR(bool a, bool b) { return !XOR(a, b); }
87
88 // A debug tool
89 #ifdef DEBUG
SONDE()90 inline void SONDE( ) {
91 static int i = 1;
92 cout << "SONDE " << i << endl;
93 ++i;
94 }
95
96 #else
SONDE()97 inline void SONDE( ) {}
98
99 #endif
100
101 /*****************************************************************************************************************************
102 * FreeFem function callers
103 * ffcalfunc : template abstract mother class with a pointer to the freefem stack and the J virtual
104 *method which computes the function
105 *****************************************************************************************************************************/
106 template< class K >
107 class ffcalfunc {
108 public:
109 Stack stack;
ffcalfunc(const ffcalfunc & f)110 ffcalfunc(const ffcalfunc &f) : stack(f.stack) {}
111
ffcalfunc(Stack _stack)112 ffcalfunc(Stack _stack) : stack(_stack) {}
113
114 virtual K J(Rn_) const = 0;
~ffcalfunc()115 virtual ~ffcalfunc( ) {}
116 };
117
118 /*****************************************************************************************************************************
119 * GeneralFunc : Most general case (specialized for sparse matrix returning functions, because
120 *IPOPT need the hessian func to take some additional parameters).
121 * @theparame: ff expression of the parameter of the ff function, computing J(x) need the
122 *associated KN to be set to the values of x
123 * @JJ : ff expression of the function
124 *****************************************************************************************************************************/
125 template< class K >
126 class GeneralFunc : public ffcalfunc< K > {
127 public:
128 Expression JJ, theparame;
GeneralFunc(const GeneralFunc & f)129 GeneralFunc(const GeneralFunc &f) : ffcalfunc< K >(f), JJ(f.JJ), theparame(f.theparame) {}
130
GeneralFunc(Stack s,Expression JJJ,Expression epar)131 GeneralFunc(Stack s, Expression JJJ, Expression epar)
132 : ffcalfunc< K >(s), JJ(JJJ), theparame(epar) {}
133
J(Rn_ x) const134 K J(Rn_ x) const {
135 KN< double > *p = GetAny< KN< double > * >((*theparame)(this->stack));
136 *p = x;
137 K ret = GetAny< K >((*JJ)(this->stack));
138 WhereStackOfPtr2Free(this->stack)->clean( );
139 return ret;
140 }
141 };
142
143 /*****************************************************************************************************************************
144 * P2ScalarFunc: encapsulate a function which is the sum of a bilinear and a linear form (no
145 *constant part since it will be used as fitness function). It also handles the case of pure
146 *quadratic or linear forms.
147 * @vf : If true J will compute 0.5xMx - bx (x is the solution of Mx = b in the
148 *unconstrained optimization process) if false xMx + bx is returned
149 * @M : FF expression of the matrix of the bilinear form (null pointer for linear form
150 *case)
151 * @b : FF expression of the vector representation of the linear part (null for pure
152 *quadratic case)
153 *****************************************************************************************************************************/
154 class P2ScalarFunc : public ffcalfunc< R > {
155 public:
156 const bool vf;
157 Expression M, b; // Matrix of the quadratic part, vector of the linear part
P2ScalarFunc(const P2ScalarFunc & f)158 P2ScalarFunc(const P2ScalarFunc &f) : ffcalfunc< R >(f), M(f.M), b(f.b), vf(f.vf) {}
159
P2ScalarFunc(Stack s,Expression _M,Expression _b,bool _vf=false)160 P2ScalarFunc(Stack s, Expression _M, Expression _b, bool _vf = false)
161 : ffcalfunc< R >(s), M(_M), b(_b), vf(_vf) {}
162
J(Rn_ x) const163 R J(Rn_ x) const {
164 Rn tmp(x.N( ), 0.);
165
166 if (M) {
167 Matrice_Creuse< R > *a = GetAny< Matrice_Creuse< R > * >((*M)(stack));
168 MatriceMorse< R > *A = dynamic_cast< MatriceMorse< R > * >(&(*a->A));
169 assert(A);
170 tmp = (*A) * x;
171 if (vf) {
172 tmp /= 2.;
173 }
174 }
175
176 if (b) {
177 Rn *B = GetAny< Rn * >((*b)(stack));
178 tmp += *B;
179 }
180
181 R res = 0.;
182
183 for (int i = 0; i < x.N( ); ++i) {
184 res += x[i] * tmp[i];
185 }
186
187 return res;
188 }
189 };
190
191 /*****************************************************************************************************************************
192 * P1VectorFunc: encapsulate a function which is the sum of a linear part and a constant,
193 *mostly used for affine/linear constraints, or for P2 fitness function gradient
194 * @vf : Set to true if this is expected the gradient of a P2 scalar function associated to
195 *Ax=b linear system J will then return Ax - b. Otherwize Ax+b is returned.
196 * @M : FF expression of the matrix of the linear part
197 * @b : FF expression of the vector representation of the constant part
198 *****************************************************************************************************************************/
199 class P1VectorFunc : public ffcalfunc< Rn > {
200 public:
201 const bool vf;
202 Expression M, b;
P1VectorFunc(const P1VectorFunc & f)203 P1VectorFunc(const P1VectorFunc &f) : ffcalfunc< Rn >(f), M(f.M), b(f.b), vf(f.vf) {}
204
P1VectorFunc(Stack s,Expression _M,Expression _b,bool _vf=false)205 P1VectorFunc(Stack s, Expression _M, Expression _b, bool _vf = false)
206 : ffcalfunc< Rn >(s), M(_M), b(_b), vf(_vf) {}
207
J(Rn_ x) const208 Rn J(Rn_ x) const {
209 Rn tmp(0);
210
211 if (M) {
212 Matrice_Creuse< R > *a = GetAny< Matrice_Creuse< R > * >((*M)(stack));
213 MatriceMorse< R > *A = dynamic_cast< MatriceMorse< R > * >(&(*a->A));
214 assert(A);
215 if (tmp.N( ) != A->n) {
216 tmp.resize(A->n);
217 tmp = 0.;
218 }
219
220 tmp = (*A) * x;
221 }
222
223 if (b) {
224 Rn *B = GetAny< Rn * >((*b)(stack));
225 if (tmp.N( ) != B->N( )) {
226 tmp.resize(B->N( ));
227 tmp = 0.;
228 }
229
230 tmp += *B;
231 }
232
233 return tmp;
234 }
235 };
236
237 /*****************************************************************************************************************************
238 * ffcalfunc<Matrice_Creuse>R>*>: specialization for sparse matrix returning function. When it
239 *encapsulates the hessian function of the lagragian, non-linear constraints will need the
240 *additional obj_factor and lagrange multiplier parameters. The one parameter version of J is called
241 *if there is no non-linear constraints or if the objects represents the jacobian of the
242 *constraints.
243 *****************************************************************************************************************************/
244 template<>
245 class ffcalfunc< Matrice_Creuse< R > * > {
246 public:
247 typedef Matrice_Creuse< R > *K;
248 Stack stack;
ffcalfunc(const ffcalfunc & f)249 ffcalfunc(const ffcalfunc &f) : stack(f.stack) {}
250
ffcalfunc(Stack s)251 ffcalfunc(Stack s) : stack(s) {}
252
253 virtual K J(Rn_) const = 0;
254 virtual K J(Rn_, double, Rn_) const = 0;
255 virtual bool NLCHPEnabled( ) const = 0; // Non Linear Constraints Hessian Prototype
~ffcalfunc()256 virtual ~ffcalfunc( ) {}
257 };
258
259 /*****************************************************************************************************************************
260 * GeneralSparseMatFunc: general case of sparse matrix returning function. Members datas added
261 *are ff expression of the scalar objective factor and vectorial lagrange multipliers.
262 *****************************************************************************************************************************/
263 class GeneralSparseMatFunc : public ffcalfunc< Matrice_Creuse< R > * > {
264 private:
265 typedef ffcalfunc< Matrice_Creuse< R > * > FFF;
266
267 public:
268 Expression JJ, param, paramlm, paramof;
GeneralSparseMatFunc(const GeneralSparseMatFunc & f)269 GeneralSparseMatFunc(const GeneralSparseMatFunc &f)
270 : FFF(f), JJ(f.JJ), param(f.param), paramlm(f.paramlm), paramof(f.paramof){};
GeneralSparseMatFunc(Stack s,Expression JJJ,Expression epar,Expression eparof=0,Expression eparlm=0)271 GeneralSparseMatFunc(Stack s, Expression JJJ, Expression epar, Expression eparof = 0,
272 Expression eparlm = 0)
273 : FFF(s), JJ(JJJ), param(epar), paramlm(eparlm), paramof(eparof) {
274 ffassert(NXOR(paramlm, paramof));
275 }
276
NLCHPEnabled() const277 bool NLCHPEnabled( ) const { return paramlm && paramof; }
278
J(Rn_ x) const279 K J(Rn_ x) const {
280 KN< double > *p = GetAny< KN< double > * >((*param)(stack));
281 *p = x;
282 K ret = GetAny< K >((*JJ)(stack));
283 // cout << "call to ffcalfunc.J with " << *p << " and ret=" << ret << endl;
284 WhereStackOfPtr2Free(stack)->clean( );
285 return ret;
286 }
287
J(Rn_ x,double of,Rn_ lm) const288 K J(Rn_ x, double of, Rn_ lm) const {
289 if (paramlm && paramof) {
290 KN< double > *p = GetAny< KN< double > * >((*param)(stack));
291 double *pof = GetAny< double * >((*paramof)(stack));
292 KN< double > *plm = GetAny< KN< double > * >((*paramlm)(stack));
293 *p = x;
294 *pof = of;
295 int m = lm.N( ), mm = plm->N( );
296 if ((m != mm) && mm) {
297 cout << " ff-ipopt H : big bug int size ???" << m << " != " << mm << endl;
298 abort( );
299 }
300
301 ;
302 *plm = lm;
303 K ret = GetAny< K >((*JJ)(stack));
304 // cout << "call to ffcalfunc.J with " << *p << " and ret=" << ret << endl;
305 WhereStackOfPtr2Free(stack)->clean( );
306 return ret;
307 } else {
308 return J(x);
309 }
310 }
311 };
312
313 /*****************************************************************************************************************************
314 * ConstantSparseMatFunc: Encapsulate a constant matrix returning function. Just contains the
315 *ff expression of the matrix (and stack inherited from mother class), this matrix is returned
316 *regardless of x.
317 *****************************************************************************************************************************/
318 class ConstantSparseMatFunc : public ffcalfunc< Matrice_Creuse< R > * > {
319 private:
320 typedef ffcalfunc< Matrice_Creuse< R > * > FFF;
321
322 public:
323 Expression M; // Expression of the matrix
ConstantSparseMatFunc(const ConstantSparseMatFunc & f)324 ConstantSparseMatFunc(const ConstantSparseMatFunc &f) : FFF(f), M(f.M) {}
325
ConstantSparseMatFunc(Stack s,Expression _M)326 ConstantSparseMatFunc(Stack s, Expression _M) : FFF(s), M(_M) {}
327
NLCHPEnabled() const328 bool NLCHPEnabled( ) const { return false; }
329
J(Rn_) const330 K J(Rn_) const {
331 K ret = M ? GetAny< K >((*M)(stack)) : 0;
332
333 WhereStackOfPtr2Free(stack)->clean( );
334 return ret;
335 }
336
J(Rn_ x,double,Rn_) const337 K J(Rn_ x, double, Rn_) const { return J(x); }
338 };
339
340 typedef ffcalfunc< double > ScalarFunc;
341 typedef ffcalfunc< Rn > VectorFunc;
342 typedef ffcalfunc< Rnm > FullMatrixFunc;
343 typedef ffcalfunc< Matrice_Creuse< R > * > SparseMatFunc;
344
345 /*****************************************************************************************************************************
346 * SparseMatStructure: a class for sparse matrix structure management (mostly merging). The
347 *most interesting methods in this class are : AddMatrix : merge the structure of the given matrix
348 *to the structure of current object AddArrays : merge structure in arrays form to the current
349 *object ToKn : allocate the raws and cols pointers and fill them with the std::set<Z2> form of
350 *the structure structure is then emptied if this method is passed a true value
351 * ==> update 28/03/2012, autostruct proved useless since the structure merging can be done with
352 *operator + (I did not no whether nullify coefficients where removed from the result but it
353 *actually doesn't so the structure of the lagrangian hessian can be guessed exactly by evaluating
354 *on a point yeilding the biggest fitness function hessian along with a dual vector filled with 1).
355 *****************************************************************************************************************************/
356 class SparseMatStructure {
357 public:
358 typedef std::pair< int, int > Z2;
359 typedef std::set< Z2 > Structure;
360 typedef std::pair< KN< int >, KN< int > > Zn2;
361 typedef Structure::const_iterator const_iterator;
362 typedef Structure::iterator iterator;
363
SparseMatStructure(bool _sym=0)364 SparseMatStructure(bool _sym = 0) : structure( ), sym(_sym), n(0), m(0), raws(0), cols(0) {}
365
SparseMatStructure(Matrice_Creuse<R> * M,bool _sym=0)366 SparseMatStructure(Matrice_Creuse< R > *M, bool _sym = 0)
367 : structure( ), sym(_sym), n(M->N( )), m(M->M( )), raws(0), cols(0) {
368 this->AddMatrix(M);
369 }
370
371 template< class INT >
SparseMatStructure(const KN<INT> & I,const KN<INT> & J,bool _sym=0)372 SparseMatStructure(const KN< INT > &I, const KN< INT > &J, bool _sym = 0)
373 : structure( ), sym(_sym), n(I.max( )), m(J.max( )), raws(0), cols(0) {
374 this->AddArrays(I, J);
375 }
376
~SparseMatStructure()377 ~SparseMatStructure( ) {
378 if (raws) {
379 delete raws;
380 }
381
382 if (cols) {
383 delete cols;
384 }
385 }
386
begin() const387 const_iterator begin( ) const { return structure.begin( ); }
388
begin()389 iterator begin( ) { return structure.begin( ); }
390
end() const391 const_iterator end( ) const { return structure.end( ); }
392
end()393 iterator end( ) { return structure.end( ); }
394
395 // Structure& operator()() {return structure;}
396 // const Structure& operator()() const {return structure;}
empty() const397 bool empty( ) const { return structure.empty( ) && !raws && !cols; }
398
N() const399 int N( ) const { return n; }
400
M() const401 int M( ) const { return m; }
402
clear()403 SparseMatStructure &clear( ) {
404 structure.clear( );
405 if (raws) {
406 delete raws;
407 }
408
409 if (cols) {
410 delete cols;
411 }
412
413 sym = false;
414 n = 0;
415 m = 0;
416 return *this;
417 }
418
size() const419 int size( ) const { return structure.size( ) ? structure.size( ) : (raws ? raws->N( ) : 0); }
420
421 SparseMatStructure &AddMatrix(Matrice_Creuse< R > *);
422 template< class INT >
423 SparseMatStructure &AddArrays(const KN< INT > &, const KN< INT > &);
424 SparseMatStructure &ToKn(bool emptystruct = false);
425
Raws()426 KN< int > &Raws( ) { return *raws; }
427
Raws() const428 KN< int > const &Raws( ) const { return *raws; }
429
Cols()430 KN< int > &Cols( ) { return *cols; }
431
Cols() const432 KN< int > const &Cols( ) const { return *cols; }
433
434 private:
435 int n, m;
436 Structure structure;
437 bool sym;
438 KN< int > *raws, *cols;
439 };
440
ToKn(bool emptystruct)441 SparseMatStructure &SparseMatStructure::ToKn(bool emptystruct) {
442 if (raws) {
443 delete raws;
444 }
445
446 if (cols) {
447 delete cols;
448 }
449
450 raws = new KN< int >(structure.size( ));
451 cols = new KN< int >(structure.size( ));
452 int k = 0;
453
454 for (const_iterator i = begin( ); i != end( ); ++i) {
455 (*raws)[k] = i->first;
456 (*cols)[k] = i->second;
457 ++k;
458 }
459
460 if (emptystruct) {
461 structure.clear( );
462 }
463
464 return *this;
465 }
466
AddMatrix(Matrice_Creuse<R> * const _M)467 SparseMatStructure &SparseMatStructure::AddMatrix(Matrice_Creuse< R > *const _M) {
468 n = n > _M->N( ) ? n : _M->N( );
469 m = m > _M->M( ) ? m : _M->M( );
470 MatriceMorse< R > *M = _M->pHM( );
471 if (!M) {
472 cerr << " Err= "
473 << " Matrix is not morse or CSR " << &(*_M->A) << endl;
474 ffassert(M);
475 }
476 M->CSR( );
477 {
478 if (!sym || (sym && M->half)) {
479 for (int i = 0; i < M->N; ++i) {
480 for (int k = M->p[i]; k < M->p[i + 1]; ++k) {
481 structure.insert(Z2(i, M->j[k]));
482 }
483 }
484 } else { // sym && !M->symetrique
485 for (int i = 0; i < M->N; ++i) {
486 for (int k = M->p[i]; k < M->p[i + 1]; ++k) {
487 if (i >= M->j[k]) {
488 structure.insert(Z2(i, M->j[k]));
489 }
490 }
491 }
492 }
493 }
494 return *this;
495 }
496
497 template< class INT >
AddArrays(const KN<INT> & I,const KN<INT> & J)498 SparseMatStructure &SparseMatStructure::AddArrays(const KN< INT > &I, const KN< INT > &J) {
499 ffassert(I.N( ) == J.N( ));
500 n = n > I.max( ) + 1 ? n : I.max( ) + 1;
501 m = m > J.max( ) + 1 ? m : J.max( ) + 1;
502 if (!sym) {
503 for (int k = 0; k < I.N( ); ++k) {
504 structure.insert(Z2(I[k], J[k]));
505 }
506 } else {
507 for (int k = 0; k < I.N( ); ++k) {
508 if (I[k] >= J[k]) {
509 structure.insert(Z2(I[k], J[k]));
510 }
511 }
512 }
513
514 return *this;
515 }
516
517 /*****************************************************************************************************************************
518 * ffNLP : Derived from the TNLP non-linear problem wrapper class of Ipopt. Virtual methods are
519 *defined as explain in the IPOPT documentation. Some of them are tricky because the sparse matrix
520 *format in freefem is CRS, whereas IPOPT use COO storage. It is even more tricky because most of
521 *time, FreeFem will remove null coefficient from the structure, leading to non constant indexing of
522 *the coefficient through the algorithm in case of very non linear functions. As IPOPT need a
523 *constant structure, a FindIndex method involving a dichotomic search has been implemented to
524 *prevent the errors related to that.
525 *****************************************************************************************************************************/
526 using namespace Ipopt;
527
528 class ffNLP : public TNLP {
529 public:
ffNLP()530 ffNLP( ) : xstart(0) {}
531
532 ffNLP(Rn &, const Rn &, const Rn &, const Rn &, const Rn &, ScalarFunc *, VectorFunc *,
533 SparseMatFunc *, VectorFunc *, SparseMatFunc *);
534 ffNLP(Rn &, const Rn &, const Rn &, const Rn &, const Rn &, ScalarFunc *, VectorFunc *,
535 SparseMatFunc *, VectorFunc *, SparseMatFunc *, int, int, int);
536 virtual ~ffNLP( );
537
538 bool get_nlp_info(Index &, Index &, Index &, Index &, IndexStyleEnum &); // the IPOPT methods
539 bool get_bounds_info(Index, Number *, Number *, Index, Number *, Number *);
540 bool get_starting_point(Index, bool, Number *, bool, Number *, Number *, Index, bool, Number *);
541 bool eval_f(Index, const Number *, bool, Number &);
542 bool eval_grad_f(Index, const Number *, bool, Number *);
543 bool eval_g(Index, const Number *, bool, Index, Number *);
544 bool eval_jac_g(Index, const Number *, bool, Index, Index, Index *, Index *, Number *);
545 bool eval_h(Index, const Number *, bool, Number, Index, const Number *, bool, Index, Index *,
546 Index *, Number *);
547 void finalize_solution(SolverReturn, Index, const Number *, const Number *, const Number *, Index,
548 const Number *, const Number *, Number, const IpoptData *ip_data,
549 IpoptCalculatedQuantities *ip_cq);
550
551 template< class INT >
552 ffNLP &SetHessianStructure(const KN< INT > &, const KN< INT > &, bool reset = 0);
553 template< class INT >
554 ffNLP &SetJacobianStructure(const KN< INT > &, const KN< INT > &, bool reset = 0);
555 enum Level { do_nothing, user_defined, one_evaluation, basis_analysis };
556 ffNLP &BuildMatrixStructures(Level, Level, int);
EnableCheckStruct()557 ffNLP &EnableCheckStruct( ) {
558 checkstruct = true;
559 return *this;
560 }
561
DisableCheckStruct()562 ffNLP &DisableCheckStruct( ) {
563 checkstruct = false;
564 return *this;
565 }
566
567 Rn lambda_start, x_start, uz_start, lz_start;
568 double sigma_start;
569 double final_value;
570
571 private:
572 // algorithm datas
573 Rn *xstart, xl, xu, gl, gu;
574 ScalarFunc *fitness; // Pointers to functions wrappers
575 VectorFunc *dfitness, *constraints;
576 SparseMatFunc *hessian, *dconstraints;
577 int mm, nnz_jac, nnz_h; // duplicated datas? did not seems to be reachable in the base class
578 // bool sym;
579 bool checkstruct;
580 SparseMatStructure HesStruct, JacStruct;
581
582 // some static functions...
583 template< class A, class B >
KnToPtr(const KN<A> & a,B * b)584 static void KnToPtr(const KN< A > &a, B *b) {
585 for (int i = 0; i < a.N( ); ++i) {
586 b[i] = a[i];
587 }
588 } // Fill a pointer with a KN
589
590 template< class A, class B >
KnFromPtr(KN<A> & a,B const * b)591 static void KnFromPtr(KN< A > &a, B const *b) {
592 for (int i = 0; i < a.N( ); ++i) {
593 a[i] = b[i];
594 }
595 } // Fill a KN with a pointer <-- to avoid the use of const_cast
596
597 static int FindIndex(const KN< int > &irow, const KN< int > &jrow, int i, int j, int kmin,
598 int kmax);
599 };
600
ffNLP(Rn & x,const Rn & _xl,const Rn & _xu,const Rn & _gl,const Rn & _gu,ScalarFunc * _fitness,VectorFunc * _dfitness,SparseMatFunc * _hessian,VectorFunc * _constraints,SparseMatFunc * _dconstraints)601 ffNLP::ffNLP(Rn &x, const Rn &_xl, const Rn &_xu, const Rn &_gl, const Rn &_gu,
602 ScalarFunc *_fitness, VectorFunc *_dfitness, SparseMatFunc *_hessian,
603 VectorFunc *_constraints, SparseMatFunc *_dconstraints)
604 : xstart(&x), xl(_xl), xu(_xu), gl(_gl), gu(_gu),
605 final_value(299792458.), // sym(0),unsymind(),
606 fitness(_fitness), dfitness(_dfitness), constraints(_constraints), uz_start( ), lz_start( ),
607 hessian(_hessian), dconstraints(_dconstraints), mm(-1), nnz_jac(-1), nnz_h(-1), HesStruct(true),
608 JacStruct(false), sigma_start(1.), lambda_start( ), x_start(x), checkstruct(1) {}
609
ffNLP(Rn & x,const Rn & _xl,const Rn & _xu,const Rn & _gl,const Rn & _gu,ScalarFunc * _fitness,VectorFunc * _dfitness,SparseMatFunc * _hessian,VectorFunc * _constraints,SparseMatFunc * _dconstraints,int _mm,int _nnz_jac,int _nnz_h)610 ffNLP::ffNLP(Rn &x, const Rn &_xl, const Rn &_xu, const Rn &_gl, const Rn &_gu,
611 ScalarFunc *_fitness, VectorFunc *_dfitness, SparseMatFunc *_hessian,
612 VectorFunc *_constraints, SparseMatFunc *_dconstraints, int _mm, int _nnz_jac,
613 int _nnz_h)
614 : xstart(&x), xl(_xl), xu(_xu), gl(_gl), gu(_gu), hessian(_hessian),
615 final_value(299792458.), // sym(0),unsymind(),
616 fitness(_fitness), dfitness(_dfitness), constraints(_constraints), dconstraints(_dconstraints),
617 uz_start( ), lz_start( ), mm(_mm), nnz_jac(_nnz_jac), nnz_h(_nnz_h), HesStruct(true),
618 JacStruct(false), sigma_start(1.), lambda_start( ), x_start(x), checkstruct(1) {}
619
~ffNLP()620 ffNLP::~ffNLP( ) {
621 /*
622 * clean(fitness);
623 * clean(dfitness);
624 * clean(constraints);
625 * clean(hessian);
626 * clean(dconstraints);
627 */
628 }
629
630 template< class INT >
SetHessianStructure(const KN<INT> & I,const KN<INT> & J,bool reset)631 ffNLP &ffNLP::SetHessianStructure(const KN< INT > &I, const KN< INT > &J, bool reset) {
632 if (reset) {
633 HesStruct.clear( );
634 }
635
636 HesStruct.AddArrays(I, J);
637 return *this;
638 }
639
640 template< class INT >
SetJacobianStructure(const KN<INT> & I,const KN<INT> & J,bool reset)641 ffNLP &ffNLP::SetJacobianStructure(const KN< INT > &I, const KN< INT > &J, bool reset) {
642 if (reset) {
643 JacStruct.clear( );
644 }
645
646 JacStruct.AddArrays(I, J);
647 return *this;
648 }
649
BuildMatrixStructures(Level hlvl,Level jlvl,int _mm)650 ffNLP &ffNLP::BuildMatrixStructures(Level hlvl, Level jlvl, int _mm) {
651 if (jlvl != do_nothing && dconstraints) {
652 if (jlvl == user_defined) {
653 ffassert(JacStruct.size( ));
654 } else if ((jlvl == one_evaluation || jlvl == basis_analysis) && dconstraints) {
655 JacStruct.AddMatrix(dconstraints->J(x_start));
656 }
657 }
658
659 if (hlvl != do_nothing && hessian) {
660 if (hlvl == user_defined) {
661 ffassert(HesStruct.size( ));
662 } else if (hlvl == one_evaluation || !hessian->NLCHPEnabled( )) {
663 Rn lms = lambda_start;
664 lms = 1.;
665 HesStruct.AddMatrix(hessian->J(x_start, sigma_start, lms));
666 } else if (hlvl == basis_analysis) {
667 {
668 Rn lambda(_mm, 0.);
669 HesStruct.AddMatrix(hessian->J(x_start, 1., lambda));
670 }
671
672 for (int i = 0; i < _mm; ++i) {
673 Rn lambda(_mm, 0.);
674 lambda[i] = 1.;
675 HesStruct.AddMatrix(hessian->J(x_start, 0., lambda));
676 lambda[i] = 0.;
677 }
678 }
679 }
680
681 JacStruct.ToKn( );
682 HesStruct.ToKn( );
683 return *this;
684 }
685
FindIndex(const KN<int> & irow,const KN<int> & jcol,int i,int j,int kmin,int kmax)686 int ffNLP::FindIndex(const KN< int > &irow, const KN< int > &jcol, int i, int j, int kmin,
687 int kmax) {
688 typedef std::pair< int, int > Z2;
689 Z2 ij(i, j), ijmin(irow[kmin], jcol[kmin]), ijmax(irow[kmax], jcol[kmax]);
690 if (abs(kmin - kmax) <= 1) {
691 if (ij == ijmin) {
692 return kmin;
693 } else if (ij == ijmax) {
694 return kmax;
695 } else {
696 return -1;
697 }
698 } else {
699 int knew = (kmin + kmax) / 2;
700 Z2 ijnew(irow[knew], jcol[knew]);
701 if (ij <= ijnew) {
702 return FindIndex(irow, jcol, i, j, kmin, knew);
703 } else {
704 return FindIndex(irow, jcol, i, j, knew, kmax);
705 }
706 }
707 }
708
get_nlp_info(Index & n,Index & m,Index & nnz_jac_g,Index & nnz_h_lag,IndexStyleEnum & index_style)709 bool ffNLP::get_nlp_info(Index &n, Index &m, Index &nnz_jac_g, Index &nnz_h_lag,
710 IndexStyleEnum &index_style) {
711 bool ret = true;
712
713 n = xstart ? xstart->N( ) : (ret = 0);
714 mm = m = constraints ? JacStruct.N( ) : 0;
715 nnz_jac = nnz_jac_g = constraints ? JacStruct.size( ) : 0;
716 nnz_h = nnz_h_lag = HesStruct.size( );
717 index_style = TNLP::C_STYLE;
718 return ret;
719 }
720
get_bounds_info(Index n,Number * x_l,Number * x_u,Index m,Number * g_l,Number * g_u)721 bool ffNLP::get_bounds_info(Index n, Number *x_l, Number *x_u, Index m, Number *g_l, Number *g_u) {
722 KnToPtr(xl, x_l);
723 KnToPtr(xu, x_u);
724 if (mm) {
725 KnToPtr(gl, g_l);
726 }
727
728 if (mm) {
729 KnToPtr(gu, g_u);
730 }
731
732 /* DEBUG
733 * cout << "constraints lower bound = (";
734 * for(int i=0;i<m;++i) cout << g_l[i] << (i<m-1 ? ',':')');
735 * cout << endl << "constraints upper bound = (";
736 * for(int i=0;i<m;++i) cout << g_u[i] << (i<m-1 ? ',':')');
737 * cout << endl;*/
738 return true;
739 }
740
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)741 bool ffNLP::get_starting_point(Index n, bool init_x, Number *x, bool init_z, Number *z_L,
742 Number *z_U, Index m, bool init_lambda, Number *lambda) {
743 assert(init_x == true);
744 assert(xstart->N( ) == n);
745 KnToPtr(*xstart, x);
746 if (init_z) {
747 if (uz_start.N( ) != n) {
748 if (xu.min( ) < 1.e19) {
749 cout << "ff-IPOPT warm start : upper simple bounds start multipliers array doesn't have "
750 "the expected size ("
751 << uz_start.N( ) << "!=" << n << ")." << endl;
752 cout << " ";
753 if (uz_start.N( ) == 0) {
754 cout << "maybe because no upper bounds multiplier has been given. " << endl;
755 }
756
757 cout << " Initializing them to 1..." << endl;
758 }
759
760 uz_start.resize(n);
761 uz_start = 1.;
762 }
763
764 if (lz_start.N( ) != n) {
765 if (xl.max( ) > -1e19) {
766 cout << "ff-IPOPT warm start : lower simple bounds start multipliers array doesn't have "
767 "the expected size ("
768 << lz_start.N( ) << "!=" << n << ")." << endl;
769 cout << " ";
770 if (lz_start.N( ) == 0) {
771 cout << "maybe because no lower bounds multiplier has been given. " << endl;
772 }
773
774 cout << " Initializing them to 1..." << endl;
775 }
776
777 lz_start.resize(n);
778 lz_start = 1.;
779 }
780
781 KnToPtr(uz_start, z_U);
782 KnToPtr(lz_start, z_L);
783 }
784
785 if (init_lambda) {
786 if (lambda_start.N( ) != m) {
787 cout << "ff-IPOPT warm start : constraints start multipliers array doesn't have the expected "
788 "size ("
789 << lambda_start.N( ) << "!=" << m << ")." << endl;
790 cout << " ";
791 if (lambda_start.N( ) == 0) {
792 cout << "maybe because no constraints multiplier has been given. " << endl;
793 }
794
795 cout << " Initializing them to 1..." << endl;
796 lambda_start.resize(m);
797 lambda_start = 1.;
798 }
799
800 KnToPtr(lambda_start, lambda);
801 }
802
803 return true;
804 }
805
eval_f(Index n,const Number * x,bool new_x,Number & obj_value)806 bool ffNLP::eval_f(Index n, const Number *x, bool new_x, Number &obj_value) {
807 assert(n == xstart->N( ));
808 Rn X(n);
809 KnFromPtr(X, x);
810 obj_value = fitness->J(X);
811 return true;
812 }
813
eval_grad_f(Index n,const Number * x,bool new_x,Number * grad_f)814 bool ffNLP::eval_grad_f(Index n, const Number *x, bool new_x, Number *grad_f) {
815 assert(n == xstart->N( ));
816 Rn X(n);
817 KnFromPtr(X, x);
818 Rn _grad_f = dfitness->J(X);
819 KnToPtr(_grad_f, grad_f);
820 return true;
821 }
822
eval_g(Index n,const Number * x,bool new_x,Index m,Number * g)823 bool ffNLP::eval_g(Index n, const Number *x, bool new_x, Index m, Number *g) {
824 Rn X(n);
825
826 KnFromPtr(X, x);
827 if (constraints) {
828 Rn _g = constraints->J(X);
829 KnToPtr(_g, g);
830 }
831
832 return true;
833 }
834
eval_jac_g(Index n,const Number * x,bool new_x,Index m,Index nele_jac,Index * iRow,Index * jCol,Number * values)835 bool ffNLP::eval_jac_g(Index n, const Number *x, bool new_x, Index m, Index nele_jac, Index *iRow,
836 Index *jCol, Number *values) {
837 assert(n == xstart->N( ));
838 Rn X(n);
839 if (x) {
840 KnFromPtr(X, x);
841 } else {
842 X = *xstart;
843 }
844
845 if (values == 0) {
846 int k = 0;
847
848 for (SparseMatStructure::const_iterator i = JacStruct.begin( ); i != JacStruct.end( ); ++i) {
849 iRow[k] = i->first;
850 jCol[k] = i->second;
851 ++k;
852 }
853 } else if (dconstraints) {
854 Matrice_Creuse< R > *M = dconstraints->J(X);
855 MatriceMorse< R > *MM = dynamic_cast< MatriceMorse< R > * >(&(*M->A)); // ugly!
856 MM->CSR( );
857 for (int i = 0; i < MM->N; ++i) {
858 for (int k = MM->p[i]; k < MM->p[i + 1]; ++k) {
859 if (checkstruct) {
860 int kipopt =
861 FindIndex(JacStruct.Raws( ), JacStruct.Cols( ), i, MM->j[k], 0, nele_jac - 1);
862 if (kipopt >= 0) {
863 values[kipopt] = MM->aij[k];
864 }
865 } else {
866 values[k] = MM->aij[k];
867 }
868 }
869 }
870 }
871
872 return true;
873 }
874
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)875 bool ffNLP::eval_h(Index n, const Number *x, bool new_x, Number obj_factor, Index m,
876 const Number *lambda, bool new_lambda, Index nele_hess, Index *iRow, Index *jCol,
877 Number *values) {
878 Rn X(n), L(m);
879
880 if (x) {
881 KnFromPtr(X, x);
882 } else {
883 X = *xstart;
884 }
885
886 if (lambda) {
887 KnFromPtr(L, lambda);
888 } else {
889 L = 0.;
890 }
891
892 bool NLCHPE = hessian->NLCHPEnabled( );
893 Number _obj_factor = NLCHPE ? 1. : obj_factor;
894 if (values == 0) {
895 int k = 0;
896
897 for (SparseMatStructure::const_iterator i = HesStruct.begin( ); i != HesStruct.end( ); ++i) {
898 iRow[k] = i->first;
899 jCol[k] = i->second;
900 ++k;
901 }
902 } else {
903 Matrice_Creuse< R > *M = 0;
904 if (NLCHPE) {
905 M = hessian->J(X, obj_factor, L);
906 } else {
907 M = hessian->J(X);
908 }
909
910 MatriceMorse< R > *MM = dynamic_cast< MatriceMorse< R > * >(&(*M->A)); // ugly!
911 MM->CSR( );
912 if (MM) {
913 if (checkstruct) {
914 for (int i = 0; i < MM->N; ++i) {
915 for (int k = MM->p[i]; k < MM->p[i + 1]; ++k) {
916 int kipopt =
917 FindIndex(HesStruct.Raws( ), HesStruct.Cols( ), i, MM->j[k], 0, nele_hess - 1);
918 if (kipopt >= 0) {
919 values[kipopt] = _obj_factor * (MM->aij[k]);
920 }
921 }
922 }
923 } else if (!MM->half) {
924 for (int i = 0, kipopt = 0; i < MM->N; ++i) {
925 for (int k = MM->p[i]; k < MM->p[i + 1]; ++k) {
926 if (i >= MM->j[k]) {
927 values[kipopt] = _obj_factor * (MM->aij[k]);
928 ++kipopt;
929 }
930 }
931 }
932 } else {
933 for (int i = 0; i < MM->N; ++i) {
934 for (int k = MM->p[i]; k < MM->p[i + 1]; ++k) {
935 values[k] = _obj_factor * (MM->aij[k]);
936 }
937 }
938 }
939 }
940 }
941
942 return true;
943 }
944
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)945 void ffNLP::finalize_solution(SolverReturn status, Index n, const Number *x, const Number *z_L,
946 const Number *z_U, Index m, const Number *g, const Number *lambda,
947 Number obj_value, const IpoptData *ip_data,
948 IpoptCalculatedQuantities *ip_cq) {
949 KnFromPtr(*xstart, x);
950 KnFromPtr(lambda_start, lambda);
951 KnFromPtr(lz_start, z_L);
952 KnFromPtr(uz_start, z_U);
953 final_value = obj_value;
954 }
955
956 /*****************************************************************************************************************************
957 * Assumptions : these are tags used as template parameters for case specific function wrapping
958 *or warning message in the interface. Some case can be added here (but some class has to be
959 *specialized for the new cases) AssumptionF : undeff --> undefined case (not used)
960 * no_assumption_f --> most general case when the fitness function and all its
961 *derivative are coded in the freefem script with the func keyword (type Polymorphic in c++). These
962 *functions are then wrapped in GeneralFunc objects. P2_f --> no longer used (it was
963 *used for fitness and its gradient defined as func in freefem script, while the hessian is a
964 *constant matrix directly given to the interface, but it leads to ambiguities). unavailable_hessian
965 *--> fitness function and its gradients coded with func in the ff script, wrapped into GeneralFunc
966 *objects, without second order derivative function. Enables the BFGS option of IPOPT. mv_P2_f -->
967 *fitness function is a P2 function which will be defined by a [matrix,vector] array. The functions
968 *are passed to the ffNLP object as P2ScalarFunc, P1VectorFunc and ConstantSparseMatFunc
969 *respectively for the fitness function, its gradient and its hessian (with all vf=1). quadratic_f
970 *--> f is a pure quadratic fonction, defined by a single matrix. Same type as mv_P2_f for function
971 *wrappers with a vf=0 tag. linear_f --> f is a linear form, defined by a single vector.
972 *Same type as mv_P2_f for function wrappers with a vf=0 tag. AssumptionG : undeff -->
973 *undefined case (not used) no_assumption_f --> most general case when the constraint functions
974 *and all its derivative are coded in the freefem script with the func keyword (type Polymorphic in
975 *c++). These functions are then wrapped in GeneralFunc objects. P1_g --> no longer
976 *used (it was used for constraints defined as func in freefem script , while the jacobian is a
977 *constant matrix directly given to the interface, but it leads to ambiguities). mv_P1_g -->
978 *Constraints function is a P1 function which will be defined by a [matrix,vector] array. The
979 *functions are passed to the ffNLP object as P1VectorFunc and ConstantSparseMatFunc respectively
980 *for the constraints and its jacobian (with all vf=0). linear_g --> Constraints are
981 *linear, defined by a single matrix. Same type as mv_P1_g for function wrappers with a vf=0 tag.
982 * Case : templatized with a pair of AssumptionF and AssumptionG, is used to build different
983 *constructor for the interface class in order to overload the freefem function which will call
984 *IPOPT
985 *****************************************************************************************************************************/
986
987 enum AssumptionF {
988 undeff,
989 no_assumption_f,
990 P2_f,
991 unavailable_hessian,
992 mv_P2_f,
993 quadratic_f,
994 linear_f
995 };
996 enum AssumptionG { undefg, without_constraints, no_assumption_g, P1_g, mv_P1_g, linear_g };
997
998 template< AssumptionF AF, AssumptionG AG >
999 struct Case {
CaseCase1000 Case( ) {}
1001
1002 static const AssumptionF af = AF;
1003 static const AssumptionG ag = AG;
1004 };
1005
1006 /*****************************************************************************************************************************
1007 * CheckMatrixVectorPair : Small function taking an E_Array and check whether the type of the 2
1008 *objects contained in the array are matrix and vector. Returns false if types are not
1009 *matrix/vector. order is modified to know whether the matrix is in first position or not.
1010 *****************************************************************************************************************************/
CheckMatrixVectorPair(const E_Array * mv,bool & order)1011 bool CheckMatrixVectorPair(const E_Array *mv, bool &order) {
1012 const aType t1 = (*mv)[0].left( ), t2 = (*mv)[1].left( );
1013
1014 if (NXOR(t1 == atype< Matrice_Creuse< R > * >( ), t2 == atype< Matrice_Creuse< R > * >( ))) {
1015 return false;
1016 } else if (NXOR(t1 == atype< Rn * >( ), t2 == atype< Rn * >( ))) {
1017 return false;
1018 } else {
1019 order = (t1 == atype< Matrice_Creuse< R > * >( ));
1020 return true;
1021 }
1022 }
1023
1024 /*****************************************************************************************************************************
1025 * The following class offers a polymorphic way to build the function wrappers to pass to the
1026 *ffNLP object Each element of the assumption enum define a "FunctionDatas" class in which the
1027 *constructor and the operator() makes case specific task. If some new value in the Assumption enums
1028 *are to be added, the FitnessFunctionDatas and/or ConstraintFunctionDatas with the new value as
1029 *template parameter has to be specialized. What should the method do is (exemple at the end of the
1030 *file with already coded cases): Constructor : define the Expression members using the arguments
1031 *passed to the IPOPT function in the script operator() : allocate with appropriate dynamic type
1032 *the ScalarFunc, VectorFunc, SparseMatFunc pointers, and display some case dependant errors or
1033 *warnings (note that there is no ScalarFunc ptr to allocate for constraints)
1034 *****************************************************************************************************************************/
1035 class GenericFitnessFunctionDatas {
1036 public:
1037 static GenericFitnessFunctionDatas *New(AssumptionF, const basicAC_F0 &, Expression const *,
1038 const C_F0 &, const C_F0 &, const C_F0 &);
1039 bool CompletelyNonLinearConstraints;
1040 Expression JJ, GradJ, Hessian;
GenericFitnessFunctionDatas()1041 GenericFitnessFunctionDatas( )
1042 : CompletelyNonLinearConstraints(true), JJ(0), GradJ(0), Hessian(0) {}
1043
A() const1044 virtual const AssumptionF A( ) const { return undeff; }
1045
1046 virtual void operator( )(Stack, const C_F0 &, const C_F0 &, const C_F0 &, Expression const *,
1047 ScalarFunc *&, VectorFunc *&, SparseMatFunc *&,
1048 bool) const = 0; // Build the functions
~GenericFitnessFunctionDatas()1049 virtual ~GenericFitnessFunctionDatas( ) {}
1050 };
1051
1052 template< AssumptionF AF >
1053 class FitnessFunctionDatas
1054 : public GenericFitnessFunctionDatas // not really a template, since most of the methods of all
1055 // cases have to be specialized
1056 {
1057 public:
1058 FitnessFunctionDatas(const basicAC_F0 &, Expression const *, const C_F0 &, const C_F0 &,
1059 const C_F0 &);
A() const1060 const AssumptionF A( ) const { return AF; }
1061
1062 void operator( )(Stack, const C_F0 &, const C_F0 &, const C_F0 &, Expression const *,
1063 ScalarFunc *&, VectorFunc *&, SparseMatFunc *&, bool) const;
1064 };
1065
1066 class GenericConstraintFunctionDatas {
1067 public:
1068 static GenericConstraintFunctionDatas *New(AssumptionG, const basicAC_F0 &, Expression const *,
1069 const C_F0 &);
1070 Expression Constraints, GradConstraints;
GenericConstraintFunctionDatas()1071 GenericConstraintFunctionDatas( ) : Constraints(0), GradConstraints(0) {}
1072
A() const1073 virtual const AssumptionG A( ) const { return undefg; }
1074
1075 virtual const bool WC( ) const = 0; // with constraints
1076 virtual void operator( )(Stack, const C_F0 &, Expression const *, VectorFunc *&, SparseMatFunc *&,
1077 bool) const = 0; // build the functions`
~GenericConstraintFunctionDatas()1078 virtual ~GenericConstraintFunctionDatas( ) {}
1079 };
1080
1081 template< AssumptionG AG >
1082 class ConstraintFunctionDatas : public GenericConstraintFunctionDatas {
1083 public:
1084 ConstraintFunctionDatas(const basicAC_F0 &, Expression const *, const C_F0 &);
A() const1085 const AssumptionG A( ) const { return AG; }
1086
WC() const1087 const bool WC( ) const { return AG != without_constraints; }
1088
1089 void operator( )(Stack, const C_F0 &, Expression const *, VectorFunc *&, SparseMatFunc *&,
1090 bool) const;
1091 };
1092
1093 /*****************************************************************************************************************************
1094 * OptimIpopt & OptimIpopt::E_Ipopt - The interface class
1095 * Do the link beetween freefem and Ipopt
1096 *****************************************************************************************************************************/
1097 class OptimIpopt : public OneOperator {
1098 public:
1099 const AssumptionF AF;
1100 const AssumptionG AG;
1101 class E_Ipopt : public E_F0mps {
1102 private:
1103 bool spurious_cases;
1104
1105 public:
1106 const AssumptionF AF;
1107 const AssumptionG AG;
1108 const bool WC;
1109 std::set< unsigned short > unused_name_param; // In some case, some parameter are usless,
1110 // this is the list of their index in nargs
1111 void InitUNP( ); // Initialize unusued_name_param at freefem compile time
1112 static basicAC_F0::name_and_type name_param[];
1113 static const int n_name_param = 29;
1114 Expression nargs[n_name_param];
1115 Expression X;
1116 mutable Rn lm;
1117 C_F0 L_m;
1118 C_F0 inittheparam, theparam, closetheparam;
1119 C_F0 initobjfact, objfact;
1120 GenericFitnessFunctionDatas *fitness_datas;
1121 GenericConstraintFunctionDatas *constraints_datas;
arg(int i,Stack stack,bool a) const1122 bool arg(int i, Stack stack, bool a) const {
1123 return nargs[i] ? GetAny< bool >((*nargs[i])(stack)) : a;
1124 }
1125
arg(int i,Stack stack,long a) const1126 long arg(int i, Stack stack, long a) const {
1127 return nargs[i] ? GetAny< long >((*nargs[i])(stack)) : a;
1128 }
1129
arg(int i,Stack stack,R a) const1130 R arg(int i, Stack stack, R a) const { return nargs[i] ? GetAny< R >((*nargs[i])(stack)) : a; }
1131
arg(int i,Stack stack,Rn_ a) const1132 Rn_ arg(int i, Stack stack, Rn_ a) const {
1133 return nargs[i] ? GetAny< Rn_ >((*nargs[i])(stack)) : a;
1134 }
1135
1136 template< typename T >
Arg(int i,Stack s) const1137 T Arg(int i, Stack s) const {
1138 return GetAny< T >((*nargs[i])(s));
1139 }
1140
E_Ipopt(const basicAC_F0 & args,AssumptionF af,AssumptionG ag)1141 E_Ipopt(const basicAC_F0 &args, AssumptionF af, AssumptionG ag)
1142 : lm( ), L_m(CPValue(lm)), AF(af), AG(ag), WC(ag != without_constraints),
1143 unused_name_param( ), spurious_cases(false), fitness_datas(0), constraints_datas(0) {
1144 InitUNP( );
1145 int nbj = args.size( ) - 1;
1146 Block::open(currentblock); // make a new block to
1147 X = to< Rn * >(args[nbj]);
1148 C_F0 X_n(args[nbj], "n");
1149 // the expression to init the theparam of all
1150 inittheparam =
1151 currentblock->NewVar< LocalVariable >("the parameter", atype< KN< R > * >( ), X_n);
1152 initobjfact = currentblock->NewVar< LocalVariable >("objective factor", atype< double * >( ));
1153 theparam = currentblock->Find("the parameter"); // the expression for the parameter
1154 objfact = currentblock->Find("objective factor");
1155 args.SetNameParam(n_name_param, name_param, nargs);
1156 fitness_datas = GenericFitnessFunctionDatas::New(
1157 AF, args, nargs, theparam, objfact, L_m); // Creates links to the freefem objects
1158 constraints_datas =
1159 GenericConstraintFunctionDatas::New(AG, args, nargs, theparam); // defining the functions
1160 spurious_cases = AG == no_assumption_g &&
1161 (AF == P2_f || AF == mv_P2_f || AF == quadratic_f || AF == linear_f);
1162 closetheparam = C_F0((Expression)Block::snewclose(currentblock), atype< void >( ));
1163 }
1164
~E_Ipopt()1165 ~E_Ipopt( ) {
1166 if (fitness_datas) {
1167 delete fitness_datas;
1168 }
1169
1170 if (constraints_datas) {
1171 delete constraints_datas;
1172 }
1173 }
1174
operator ( )(Stack stack) const1175 virtual AnyType operator( )(Stack stack) const {
1176 double cost = nan("");
1177
1178 WhereStackOfPtr2Free(stack) = new StackOfPtr2Free(stack); // FH mars 2005
1179 Rn &x = *GetAny< Rn * >((*X)(stack));
1180 {
1181 Expression test(theparam); // in some case the KN object associated to the param is never
1182 // initialized, leading to failed assertion in KN::destroy
1183 Rn *tt = GetAny< Rn * >((*test)(stack)); // this lines prevent this to happen
1184 *tt = x;
1185 }
1186 long n = x.N( );
1187 bool warned = false;
1188 cout << endl;
1189 if (spurious_cases) {
1190 cout << "ff-IPOPT Spurious case detected : the hessian is defined as a constant matrix but "
1191 "constraints are given in function form."
1192 << endl;
1193 cout << "If they are not affine, the optimization is likely to fail. In this case, try one "
1194 "of the following suggestions:"
1195 << endl;
1196 cout << " - if constraints have computable hessians, use function form for the fitness "
1197 "function and all its derivatives"
1198 << endl;
1199 cout
1200 << " and check the documentation to know how to express the whole lagrangian hessian."
1201 << endl;
1202 cout << " - if constraints hessians are difficult to obtain, force the BFGS mode using "
1203 "named parameter "
1204 << name_param[12].name << '.' << endl;
1205 cout << "Do not worry about this message if you know all your constraints has a constant "
1206 "null hessian."
1207 << endl
1208 << endl;
1209 }
1210
1211 if (nargs[7]) {
1212 cout << "ff-IPOPT : the named parameter autostruct is no longer used in this version of "
1213 "the interface."
1214 << endl;
1215 }
1216
1217 // Detection of mixed case dependant warnings or error
1218 for (int i = 0; i < n_name_param; ++i) {
1219 if (nargs[i] && unused_name_param.find(i) != unused_name_param.end( )) {
1220 cout << "ff-IPOPT Warning: named parameter " << name_param[i].name
1221 << " is useless for the problem you have set." << endl;
1222 warned = true;
1223 }
1224 }
1225
1226 if (nargs[4] && nargs[5] && nargs[7]) {
1227 cout << "ff-IPOPT Warning: both " << name_param[4].name << " and " << name_param[5].name
1228 << " has been defined, so " << name_param[7].name;
1229 cout << " will be ignored." << endl;
1230 }
1231
1232 if (warned) {
1233 if (!WC && AF == unavailable_hessian && nargs[8]) {
1234 cout << " ==> " << name_param[8].name
1235 << " is useless because there should not be any function returning matrix in your "
1236 "problem,"
1237 << endl;
1238 cout << " (2 functions can only be J and dJ). You may as well have forgotten one "
1239 "function (IPOPT will certainly crash if so)."
1240 << endl;
1241 }
1242
1243 if (AF != no_assumption_f && AF != unavailable_hessian && AG != no_assumption_g &&
1244 nargs[5]) {
1245 cout << " ==> your lagrangian hessian is a constant matrix, so there is no need to "
1246 "specify its structure with "
1247 << name_param[5].name << endl;
1248 cout << " since it is contained in the matrix object." << endl;
1249 }
1250
1251 if (AF != no_assumption_f && AF != unavailable_hessian && AG != no_assumption_g &&
1252 nargs[7]) {
1253 cout << " ==> " << name_param[7].name
1254 << " will be ignored since all matrices are constants and constraints do not"
1255 << endl;
1256 cout << " contribute to the hessian, matrix structure determination is trivial."
1257 << endl;
1258 }
1259
1260 if (AF == unavailable_hessian && AG != no_assumption_g && (nargs[7] || nargs[8])) {
1261 cout << " ==> " << name_param[7].name << " or " << name_param[8].name
1262 << " will be ignored since the only matrix you have passed is constant. " << endl;
1263 cout << " Or maybe did you forget to pass a function (IPOPT will certainly crash if "
1264 "so)."
1265 << endl;
1266 }
1267
1268 if (AF != no_assumption_f && AF != unavailable_hessian && AG != no_assumption_g &&
1269 nargs[8]) {
1270 cout << " ==> no need to use " << name_param[8].name
1271 << " since all matrices are constant, structures won't change through the algorithm,"
1272 << endl;
1273 cout << " it is automatically set to the default disabling value." << endl;
1274 }
1275 }
1276
1277 long iprint = verbosity;
1278 ScalarFunc *ffJ = 0;
1279 VectorFunc *ffdJ = 0, *ffC = 0;
1280 SparseMatFunc *ffH = 0, *ffdC = 0;
1281
1282 (*fitness_datas)(stack, theparam, objfact, L_m, nargs, ffJ, ffdJ, ffH,
1283 warned); // Fill the functions
1284 (*constraints_datas)(stack, theparam, nargs, ffC, ffdC, warned);
1285
1286 Rn xl(n), xu(n), gl(nargs[2] ? Arg< Rn_ >(2, stack).N( ) : 0),
1287 gu(nargs[3] ? Arg< Rn_ >(3, stack).N( ) : 0);
1288 int mmm = 0;
1289 if (WC && (gl.N( ) + gu.N( )) == 0) {
1290 cout << "IPOPT Warning : constrained problem without constraints bounds." << endl;
1291 mmm = ffC->J(x).N( );
1292 } else {
1293 mmm = gl.N( ) > gu.N( ) ? gl.N( ) : gu.N( );
1294 }
1295
1296 Rn_ *lag_mul = 0, *l_z = 0, *u_z = 0; // Rn(mmm,1.);
1297 // int niter=arg(6,stack,100L);
1298 int autostructmode = ffNLP::one_evaluation;
1299 bool checkindex =
1300 (AF != no_assumption_f && AG != no_assumption_g) ? false : arg(8, stack, true),
1301 cberror = false;
1302
1303 if (nargs[0]) {
1304 xl = Arg< Rn_ >(0, stack);
1305 } else {
1306 xl = -1.e19;
1307 }
1308
1309 if (nargs[1]) {
1310 xu = Arg< Rn_ >(1, stack);
1311 } else {
1312 xu = 1.e19;
1313 }
1314
1315 if (nargs[2]) {
1316 gl = Arg< Rn_ >(2, stack);
1317 } else {
1318 gl.resize(mmm);
1319 gl = -1.e19;
1320 }
1321
1322 if (nargs[3]) {
1323 gu = Arg< Rn_ >(3, stack);
1324 } else {
1325 gu.resize(mmm);
1326 gu = 1.e19;
1327 }
1328
1329 const E_Array *ejacstruct = (WC && AF == no_assumption_f && AG == no_assumption_g && nargs[4])
1330 ? dynamic_cast< const E_Array * >(nargs[4])
1331 : 0,
1332 *ehesstruct = (AF == no_assumption_f && nargs[5])
1333 ? dynamic_cast< const E_Array * >(nargs[5])
1334 : 0;
1335
1336 if (nargs[6] && WC) {
1337 lag_mul = new Rn_(GetAny< Rn_ >((*nargs[6])(stack)));
1338 }
1339
1340 if (nargs[21]) {
1341 l_z = new Rn_(GetAny< Rn_ >((*nargs[21])(stack)));
1342 }
1343
1344 if (nargs[20]) {
1345 u_z = new Rn_(GetAny< Rn_ >((*nargs[20])(stack)));
1346 }
1347
1348 SmartPtr< TNLP > optim = new ffNLP(x, xl, xu, gl, gu, ffJ, ffdJ, ffH, ffC, ffdC);
1349 ffNLP *_optim = dynamic_cast< ffNLP * >(&(*optim));
1350 assert(_optim);
1351 if (WC && nargs[6]) {
1352 _optim->lambda_start = *lag_mul;
1353 } else if (WC) {
1354 _optim->lambda_start.resize(mmm);
1355 _optim->lambda_start = 1.;
1356 }
1357
1358 _optim->sigma_start = 1.;
1359 if (nargs[21] && nargs[0]) {
1360 _optim->lz_start = *l_z;
1361 } else if (nargs[0]) {
1362 _optim->lz_start.resize(xl.N( ));
1363 _optim->lz_start = 1.;
1364 }
1365
1366 if (nargs[20] && nargs[1]) {
1367 _optim->uz_start = *u_z;
1368 } else if (nargs[1]) {
1369 _optim->uz_start.resize(xu.N( ));
1370 _optim->uz_start = 1.;
1371 }
1372
1373 if (ejacstruct) {
1374 if (ejacstruct->nbitem( ) != 2) {
1375 ExecError(
1376 "\nSorry, we were expecting an array with two componants in structjac=[iraw,jcol]");
1377 }
1378
1379 if ((*ejacstruct)[0].left( ) != atype< KN< long > * >( )) {
1380 CompileError("Sorry, array componants in structjac=[iraw,jcol] must be integer arrays");
1381 }
1382
1383 if ((*ejacstruct)[1].left( ) != atype< KN< long > * >( )) {
1384 CompileError("Sorry, array componants in structjac=[iraw,jcol] must be integer arrays");
1385 }
1386
1387 Expression raws = (*ejacstruct)[0], cols = (*ejacstruct)[1];
1388 _optim->SetJacobianStructure(*GetAny< KN< long > * >((*raws)(stack)),
1389 *GetAny< KN< long > * >((*cols)(stack)), true);
1390 }
1391
1392 if (ehesstruct) {
1393 if (ehesstruct->nbitem( ) != 2) {
1394 ExecError(
1395 "\nSorry, we were expecting an array with two componants in structhess=[iraw,jcol]");
1396 }
1397
1398 if ((*ehesstruct)[0].left( ) != atype< KN< long > * >( )) {
1399 CompileError("Sorry, array componants in structhess=[iraw,jcol] must be integer arrays");
1400 }
1401
1402 if ((*ehesstruct)[1].left( ) != atype< KN< long > * >( )) {
1403 CompileError("Sorry, array componants in structhess=[iraw,jcol] must be integer arrays");
1404 }
1405
1406 Expression raws = (*ehesstruct)[0], cols = (*ehesstruct)[1];
1407 _optim->SetHessianStructure(*GetAny< KN< long > * >((*raws)(stack)),
1408 *GetAny< KN< long > * >((*cols)(stack)), true);
1409 }
1410
1411 ffNLP::Level lh = ehesstruct ? ffNLP::user_defined : ffNLP::Level(autostructmode),
1412 lj = ejacstruct ? ffNLP::user_defined : ffNLP::Level(autostructmode);
1413 if (AF == unavailable_hessian) {
1414 lh = ffNLP::do_nothing;
1415 }
1416
1417 _optim->BuildMatrixStructures(lh, lj, mmm);
1418 if (checkindex) {
1419 _optim->EnableCheckStruct( );
1420 }
1421
1422 SmartPtr< IpoptApplication > app = new IpoptApplication( );
1423
1424 if (nargs[9]) {
1425 app->Options( )->SetNumericValue("tol", GetAny< double >((*nargs[9])(stack)));
1426 }
1427
1428 if (nargs[10]) {
1429 app->Options( )->SetIntegerValue("max_iter", GetAny< long >((*nargs[10])(stack)));
1430 }
1431
1432 if (nargs[11]) {
1433 app->Options( )->SetNumericValue("max_cpu_time", GetAny< double >((*nargs[11])(stack)));
1434 }
1435
1436 bool bfgs = nargs[12] ? GetAny< bool >((*nargs[12])(stack)) : false;
1437 if (AF == unavailable_hessian || bfgs) {
1438 if (AF == unavailable_hessian && !bfgs) {
1439 cout << "IPOPT Note : No hessian given ==> LBFGS hessian approximation enabled" << endl;
1440 }
1441
1442 app->Options( )->SetStringValue("hessian_approximation", "limited-memory");
1443 }
1444
1445 if (nargs[13]) {
1446 string derivative_test = *GetAny< string * >((*nargs[13])(stack));
1447 app->Options( )->SetStringValue("derivative_test", derivative_test.c_str( ));
1448 }
1449
1450 if (nargs[14]) {
1451 string options_file = *GetAny< string * >((*nargs[14])(stack));
1452 app->Options( )->SetStringValue("option_file_name", options_file.c_str( ));
1453 }
1454
1455 if (nargs[15]) {
1456 app->Options( )->SetIntegerValue("print_level", GetAny< long >((*nargs[15])(stack)));
1457 }
1458
1459 if (AG == without_constraints || AG == mv_P1_g || AG == linear_g) {
1460 app->Options( )->SetStringValue("jac_c_constant", "yes");
1461 app->Options( )->SetStringValue("jac_d_constant", "yes");
1462 }
1463
1464 if (AF == mv_P2_f || AF == quadratic_f || AF == linear_f) {
1465 app->Options( )->SetStringValue("hessian_constant", "yes");
1466 }
1467
1468 if (nargs[16]) {
1469 app->Options( )->SetNumericValue("derivative_test_perturbation",
1470 GetAny< double >((*nargs[16])(stack)));
1471 }
1472
1473 if (nargs[17]) {
1474 app->Options( )->SetNumericValue("derivative_test_tol",
1475 GetAny< double >((*nargs[17])(stack)));
1476 }
1477
1478 if (nargs[18]) {
1479 app->Options( )->SetStringValue("fixed_variable_treatment",
1480 GetAny< string * >((*nargs[18])(stack))->c_str( ));
1481 }
1482
1483 if (nargs[19]) {
1484 app->Options( )->SetStringValue("warm_start_init_point", "yes");
1485 if (WC && !nargs[6]) {
1486 cout << "ff-IPOPT Warning : warm start for constrained problem without initial "
1487 "constraints dual variables ("
1488 << name_param[6].name << " parameter)." << endl;
1489 cout << " ==> Starting with " << name_param[6].name << "=(1,1,...,1)."
1490 << endl;
1491 }
1492
1493 if (nargs[0] && !nargs[21]) {
1494 cout << "ff-IPOPT Warning : warm start with simple lower bounds without initial lower "
1495 "bounds dual variables ("
1496 << name_param[21].name << " parameter)." << endl;
1497 cout << " ==> Starting with " << name_param[21].name << "=(1,1,...,1)."
1498 << endl;
1499 }
1500
1501 if (nargs[1] && !nargs[20]) {
1502 cout << "ff-IPOPT Warning : warm start with simple upper bounds without initial upper "
1503 "bounds dual variables ("
1504 << name_param[20].name << " parameter)." << endl;
1505 cout << " ==> Starting with " << name_param[20].name << "=(1,1,...,1)."
1506 << endl;
1507 }
1508
1509 if (l_z) {
1510 _optim->lz_start = *l_z;
1511 }
1512
1513 if (u_z) {
1514 _optim->uz_start = *u_z;
1515 }
1516
1517 if (lag_mul) {
1518 _optim->lambda_start = *lag_mul;
1519 }
1520 }
1521
1522 if (nargs[22]) {
1523 app->Options( )->SetNumericValue("mu_init", GetAny< double >((*nargs[22])(stack)));
1524 } else {
1525 app->Options( )->SetStringValue("mu_strategy", "adaptive");
1526 }
1527
1528 if (nargs[23]) {
1529 app->Options( )->SetNumericValue("mumps_pivtol", GetAny< double >((*nargs[23])(stack)));
1530 }
1531
1532 if (nargs[24]) {
1533 app->Options( )->SetNumericValue("bound_relax_factor",
1534 GetAny< double >((*nargs[24])(stack)));
1535 }
1536
1537 if (nargs[25]) {
1538 app->Options( )->SetStringValue("mu_strategy",
1539 GetAny< string * >((*nargs[25])(stack))->c_str( ));
1540 }
1541
1542 if (nargs[27]) {
1543 app->Options( )->SetNumericValue("mu_min", GetAny< double >((*nargs[27])(stack)));
1544 }
1545
1546 if (nargs[28]) {
1547 if (!GetAny< bool >((*nargs[28])(stack))) {
1548 app->Options( )->SetStringValue("accept_every_trial_step", "yes");
1549 }
1550 }
1551
1552 if (verbosity > 1) {
1553 app->Options( )->SetStringValue("print_user_options", "yes");
1554 }
1555
1556 app->Options( )->SetStringValue("output_file", "ipopt.out");
1557 if (AF != no_assumption_f && AF != unavailable_hessian && AG != no_assumption_g) {
1558 app->Options( )->SetStringValue("mehrotra_algorithm", "yes");
1559 }
1560
1561 ApplicationReturnStatus status;
1562 app->Initialize( );
1563
1564 // Ask Ipopt to solve the problem
1565 status = app->OptimizeTNLP(optim);
1566
1567 if (lag_mul) {
1568 *lag_mul = _optim->lambda_start;
1569 }
1570
1571 if (l_z) {
1572 *l_z = _optim->lz_start;
1573 }
1574
1575 if (u_z) {
1576 *u_z = _optim->uz_start;
1577 }
1578
1579 cost = _optim->final_value;
1580
1581 if (nargs[26]) {
1582 double *pfv = GetAny< double * >((*nargs[26])(stack));
1583 *pfv = cost;
1584 }
1585
1586 if (verbosity) {
1587 if (status == Solve_Succeeded) {
1588 printf("\n\n*** Ipopt succeeded \n");
1589 } else if (static_cast< int >(status) < 0) {
1590 printf("\n\n*** Ipopt failure!\n");
1591 } else {
1592 printf("\n\n*** Ipopt mixed results.\n");
1593 }
1594 }
1595
1596 clean(lag_mul);
1597 clean(l_z);
1598 clean(u_z);
1599 clean(ffJ);
1600 clean(ffdJ);
1601 clean(ffH);
1602 clean(ffC);
1603 clean(ffdC);
1604 if (lm) {
1605 lm.destroy( ); // clean memory of LM
1606 }
1607
1608 closetheparam.eval(stack); // clean memory
1609 WhereStackOfPtr2Free(stack)->clean( ); // FH mars 2005
1610 return SetAny< long >(static_cast< long >(
1611 static_cast< int >(status))); // SetAny<long>(0); Modif FH july 2005
1612 }
1613
operator aType() const1614 operator aType( ) const { return atype< long >( ); }
1615 };
1616
code(const basicAC_F0 & args) const1617 E_F0 *code(const basicAC_F0 &args) const { return new E_Ipopt(args, AF, AG); }
1618
1619 // Constructors - they define the different prototype of the overloaded IPOPT function reachable
1620 // in freefem scripts
1621
OptimIpopt(Case<no_assumption_f,no_assumption_g>)1622 OptimIpopt(Case< no_assumption_f, no_assumption_g >)
1623 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1624 atype< Polymorphic * >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1625 atype< KN< R > * >( )),
1626 AF(no_assumption_f), AG(no_assumption_g) {}
1627
OptimIpopt(Case<no_assumption_f,without_constraints>)1628 OptimIpopt(Case< no_assumption_f, without_constraints >)
1629 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1630 atype< Polymorphic * >( ), atype< KN< R > * >( )),
1631 AF(no_assumption_f), AG(without_constraints) {}
1632
OptimIpopt(Case<no_assumption_f,P1_g>)1633 OptimIpopt(Case< no_assumption_f, P1_g >)
1634 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1635 atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1636 atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1637 AF(no_assumption_f), AG(P1_g) {}
1638
OptimIpopt(Case<no_assumption_f,mv_P1_g>)1639 OptimIpopt(Case< no_assumption_f, mv_P1_g >)
1640 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1641 atype< Polymorphic * >( ), atype< E_Array >( ), atype< KN< R > * >( )),
1642 AF(no_assumption_f), AG(mv_P1_g) {}
1643
OptimIpopt(Case<no_assumption_f,linear_g>)1644 OptimIpopt(Case< no_assumption_f, linear_g >)
1645 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1646 atype< Polymorphic * >( ), atype< Matrice_Creuse< R > * >( ),
1647 atype< KN< R > * >( )),
1648 AF(no_assumption_f), AG(linear_g) {}
1649
OptimIpopt(Case<P2_f,P1_g>)1650 OptimIpopt(Case< P2_f, P1_g >)
1651 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1652 atype< Matrice_Creuse< R > * >( ), atype< Polymorphic * >( ),
1653 atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1654 AF(P2_f), AG(P1_g) {}
1655
OptimIpopt(Case<P2_f,without_constraints>)1656 OptimIpopt(Case< P2_f, without_constraints >)
1657 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1658 atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1659 AF(P2_f), AG(without_constraints) {}
1660
OptimIpopt(Case<P2_f,no_assumption_g>)1661 OptimIpopt(Case< P2_f, no_assumption_g >)
1662 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1663 atype< Matrice_Creuse< R > * >( ), atype< Polymorphic * >( ),
1664 atype< Polymorphic * >( ), atype< KN< R > * >( )),
1665 AF(P2_f), AG(no_assumption_g) {}
1666
OptimIpopt(Case<P2_f,mv_P1_g>)1667 OptimIpopt(Case< P2_f, mv_P1_g >)
1668 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1669 atype< Matrice_Creuse< R > * >( ), atype< E_Array >( ), atype< KN< R > * >( )),
1670 AF(P2_f), AG(mv_P1_g) {}
1671
OptimIpopt(Case<P2_f,linear_g>)1672 OptimIpopt(Case< P2_f, linear_g >)
1673 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1674 atype< Matrice_Creuse< R > * >( ), atype< Matrice_Creuse< R > * >( ),
1675 atype< KN< R > * >( )),
1676 AF(P2_f), AG(linear_g) {}
1677
OptimIpopt(Case<unavailable_hessian,no_assumption_g>)1678 OptimIpopt(Case< unavailable_hessian, no_assumption_g >)
1679 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1680 atype< Polymorphic * >( ), atype< Polymorphic * >( ), atype< KN< R > * >( )),
1681 AF(unavailable_hessian), AG(no_assumption_g) {}
1682
OptimIpopt(Case<unavailable_hessian,without_constraints>)1683 OptimIpopt(Case< unavailable_hessian, without_constraints >)
1684 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1685 atype< KN< R > * >( )),
1686 AF(unavailable_hessian), AG(without_constraints) {}
1687
OptimIpopt(Case<unavailable_hessian,P1_g>)1688 OptimIpopt(Case< unavailable_hessian, P1_g >)
1689 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1690 atype< Polymorphic * >( ), atype< Matrice_Creuse< R > * >( ),
1691 atype< KN< R > * >( )),
1692 AF(unavailable_hessian), AG(P1_g) {}
1693
OptimIpopt(Case<unavailable_hessian,mv_P1_g>)1694 OptimIpopt(Case< unavailable_hessian, mv_P1_g >)
1695 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1696 atype< E_Array >( ), atype< KN< R > * >( )),
1697 AF(unavailable_hessian), AG(mv_P1_g) {}
1698
OptimIpopt(Case<unavailable_hessian,linear_g>)1699 OptimIpopt(Case< unavailable_hessian, linear_g >)
1700 : OneOperator(atype< long >( ), atype< Polymorphic * >( ), atype< Polymorphic * >( ),
1701 atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1702 AF(unavailable_hessian), AG(linear_g) {}
1703
OptimIpopt(Case<mv_P2_f,no_assumption_g>)1704 OptimIpopt(Case< mv_P2_f, no_assumption_g >)
1705 : OneOperator(atype< long >( ), atype< E_Array >( ), atype< Polymorphic * >( ),
1706 atype< Polymorphic * >( ), atype< KN< R > * >( )),
1707 AF(mv_P2_f), AG(no_assumption_g) {}
1708
OptimIpopt(Case<mv_P2_f,without_constraints>)1709 OptimIpopt(Case< mv_P2_f, without_constraints >)
1710 : OneOperator(atype< long >( ), atype< E_Array >( ), atype< KN< R > * >( )), AF(mv_P2_f),
1711 AG(without_constraints) {}
1712
OptimIpopt(Case<mv_P2_f,P1_g>)1713 OptimIpopt(Case< mv_P2_f, P1_g >)
1714 : OneOperator(atype< long >( ), atype< E_Array >( ), atype< Polymorphic * >( ),
1715 atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1716 AF(mv_P2_f), AG(P1_g) {}
1717
OptimIpopt(Case<mv_P2_f,mv_P1_g>)1718 OptimIpopt(Case< mv_P2_f, mv_P1_g >)
1719 : OneOperator(atype< long >( ), atype< E_Array >( ), atype< E_Array >( ),
1720 atype< KN< R > * >( )),
1721 AF(mv_P2_f), AG(mv_P1_g) {}
1722
OptimIpopt(Case<mv_P2_f,linear_g>)1723 OptimIpopt(Case< mv_P2_f, linear_g >)
1724 : OneOperator(atype< long >( ), atype< E_Array >( ), atype< Matrice_Creuse< R > * >( ),
1725 atype< KN< R > * >( )),
1726 AF(mv_P2_f), AG(linear_g) {}
1727
OptimIpopt(Case<quadratic_f,no_assumption_g>)1728 OptimIpopt(Case< quadratic_f, no_assumption_g >)
1729 : OneOperator(atype< long >( ), atype< Matrice_Creuse< R > * >( ), atype< Polymorphic * >( ),
1730 atype< Polymorphic * >( ), atype< KN< R > * >( )),
1731 AF(quadratic_f), AG(no_assumption_g) {}
1732
OptimIpopt(Case<quadratic_f,without_constraints>)1733 OptimIpopt(Case< quadratic_f, without_constraints >)
1734 : OneOperator(atype< long >( ), atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1735 AF(quadratic_f), AG(without_constraints) {}
1736
OptimIpopt(Case<quadratic_f,P1_g>)1737 OptimIpopt(Case< quadratic_f, P1_g >)
1738 : OneOperator(atype< long >( ), atype< Matrice_Creuse< R > * >( ), atype< Polymorphic * >( ),
1739 atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1740 AF(quadratic_f), AG(P1_g) {}
1741
OptimIpopt(Case<quadratic_f,mv_P1_g>)1742 OptimIpopt(Case< quadratic_f, mv_P1_g >)
1743 : OneOperator(atype< long >( ), atype< Matrice_Creuse< R > * >( ), atype< E_Array >( ),
1744 atype< KN< R > * >( )),
1745 AF(quadratic_f), AG(mv_P1_g) {}
1746
OptimIpopt(Case<quadratic_f,linear_g>)1747 OptimIpopt(Case< quadratic_f, linear_g >)
1748 : OneOperator(atype< long >( ), atype< Matrice_Creuse< R > * >( ),
1749 atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1750 AF(quadratic_f), AG(linear_g) {}
1751
OptimIpopt(Case<linear_f,no_assumption_g>)1752 OptimIpopt(Case< linear_f, no_assumption_g >)
1753 : OneOperator(atype< long >( ), atype< KN< R > * >( ), atype< Polymorphic * >( ),
1754 atype< Polymorphic * >( ), atype< KN< R > * >( )),
1755 AF(linear_f), AG(no_assumption_g) {}
1756
OptimIpopt(Case<linear_f,without_constraints>)1757 OptimIpopt(Case< linear_f, without_constraints >)
1758 : OneOperator(atype< long >( ), atype< KN< R > * >( ), atype< KN< R > * >( )), AF(linear_f),
1759 AG(without_constraints) {}
1760
OptimIpopt(Case<linear_f,P1_g>)1761 OptimIpopt(Case< linear_f, P1_g >)
1762 : OneOperator(atype< long >( ), atype< KN< R > * >( ), atype< Polymorphic * >( ),
1763 atype< Matrice_Creuse< R > * >( ), atype< KN< R > * >( )),
1764 AF(linear_f), AG(P1_g) {}
1765
OptimIpopt(Case<linear_f,mv_P1_g>)1766 OptimIpopt(Case< linear_f, mv_P1_g >)
1767 : OneOperator(atype< long >( ), atype< KN< R > * >( ), atype< E_Array >( ),
1768 atype< KN< R > * >( )),
1769 AF(linear_f), AG(mv_P1_g) {}
1770
OptimIpopt(Case<linear_f,linear_g>)1771 OptimIpopt(Case< linear_f, linear_g >)
1772 : OneOperator(atype< long >( ), atype< KN< R > * >( ), atype< Matrice_Creuse< R > * >( ),
1773 atype< KN< R > * >( )),
1774 AF(linear_f), AG(linear_g) {}
1775 };
1776
1777 /*
1778 * enum AssumptionF {no_assumption_f, P2_f, unavailable_hessian, mv_P2_f, quadratic_f,linear_f};
1779 * enum AssumptionG {without_constraints, no_assumption_g, P1_g, mv_P1_g, linear_g};
1780 */
1781
1782 // Case dependant initialization of unused_name_param
1783 // exemple : AF==no_assumption_f && AG==without_constraints ==> no constraint related named
1784 // parameter should be used (index 2,3,4,6 - first integer is how many are not used)
InitUNP()1785 void OptimIpopt::E_Ipopt::InitUNP( ) {
1786 if (AF == no_assumption_f && AG == no_assumption_g) {
1787 } // no unused named parameter
1788
1789 if (AF == no_assumption_f && AG == without_constraints) {
1790 AddElements(unused_name_param, 4, 2, 3, 4, 6);
1791 }
1792
1793 if (AF == no_assumption_f && AG == P1_g) {
1794 AddElements(unused_name_param, 1, 4);
1795 }
1796
1797 if (AF == no_assumption_f && AG == mv_P1_g) {
1798 AddElements(unused_name_param, 1, 4);
1799 }
1800
1801 if (AF == no_assumption_f && AG == linear_g) {
1802 AddElements(unused_name_param, 1, 4);
1803 }
1804
1805 if (AF == P2_f && AG == P1_g) {
1806 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1807 }
1808
1809 if (AF == P2_f && AG == without_constraints) {
1810 AddElements(unused_name_param, 8, 2, 3, 4, 5, 6, 7, 8, 12);
1811 }
1812
1813 if (AF == P2_f && AG == no_assumption_g) {
1814 AddElements(unused_name_param, 1, 5);
1815 }
1816
1817 if (AF == P2_f && AG == mv_P1_g) {
1818 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1819 }
1820
1821 if (AF == P2_f && AG == linear_g) {
1822 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1823 }
1824
1825 if (AF == unavailable_hessian && AG == no_assumption_g) {
1826 AddElements(unused_name_param, 1, 5);
1827 }
1828
1829 if (AF == unavailable_hessian && AG == without_constraints) {
1830 AddElements(unused_name_param, 7, 2, 3, 4, 5, 6, 7, 8);
1831 }
1832
1833 if (AF == unavailable_hessian && AG == P1_g) {
1834 AddElements(unused_name_param, 4, 4, 5, 7, 8);
1835 }
1836
1837 if (AF == unavailable_hessian && AG == mv_P1_g) {
1838 AddElements(unused_name_param, 4, 4, 5, 7, 8);
1839 }
1840
1841 if (AF == unavailable_hessian && AG == linear_g) {
1842 AddElements(unused_name_param, 4, 4, 5, 7, 8);
1843 }
1844
1845 if (AF == mv_P2_f && AG == without_constraints) {
1846 AddElements(unused_name_param, 8, 2, 3, 4, 5, 6, 7, 8, 12);
1847 }
1848
1849 if (AF == mv_P2_f && AG == no_assumption_g) {
1850 AddElements(unused_name_param, 1, 5);
1851 }
1852
1853 if (AF == mv_P2_f && AG == P1_g) {
1854 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1855 }
1856
1857 if (AF == mv_P2_f && AG == mv_P1_g) {
1858 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1859 }
1860
1861 if (AF == mv_P2_f && AG == linear_g) {
1862 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1863 }
1864
1865 if (AF == quadratic_f && AG == without_constraints) {
1866 AddElements(unused_name_param, 8, 2, 3, 4, 5, 6, 7, 8, 12);
1867 }
1868
1869 if (AF == quadratic_f && AG == no_assumption_g) {
1870 AddElements(unused_name_param, 1, 5);
1871 }
1872
1873 if (AF == quadratic_f && AG == P1_g) {
1874 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1875 }
1876
1877 if (AF == quadratic_f && AG == mv_P1_g) {
1878 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1879 }
1880
1881 if (AF == quadratic_f && AG == linear_g) {
1882 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1883 }
1884
1885 if (AF == linear_f && AG == without_constraints) {
1886 AddElements(unused_name_param, 8, 2, 3, 4, 5, 6, 7, 8, 12);
1887 }
1888
1889 if (AF == linear_f && AG == no_assumption_g) {
1890 AddElements(unused_name_param, 1, 5);
1891 }
1892
1893 if (AF == linear_f && AG == P1_g) {
1894 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1895 }
1896
1897 if (AF == linear_f && AG == mv_P1_g) {
1898 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1899 }
1900
1901 if (AF == linear_f && AG == linear_g) {
1902 AddElements(unused_name_param, 5, 4, 5, 7, 8, 12);
1903 }
1904 }
1905
1906 basicAC_F0::name_and_type OptimIpopt::E_Ipopt::name_param[] = {
1907 // DONT CHANGE THE ORDER!!!! If some parameters need to be added, add them after the last one
1908 // otherwize warning message of this interface will be a mess
1909 {"lb", &typeid(KN_< double >)}, // 0 - lower bound on optimization parameter X
1910 {"ub", &typeid(KN_< double >)}, // 1 - upper bound on optimization parameter X
1911 {"clb", &typeid(KN_< double >)}, // 2 - constraints lower bounds
1912 {"cub", &typeid(KN_< double >)}, // 3 - constraints upper bounds
1913 {"structjacc", &typeid(E_Array)}, // 4 - constraints jacobian structure
1914 {"structhess", &typeid(E_Array)}, // 5 - lagrangian hessian structure
1915 {"lm", &typeid(KN_< double >)}, // 6 - lagrange multiplier (for autostruct or to get their
1916 // value at the end of the algorithm)
1917 {"autostruct", &typeid(long)}, // 7 - automatic structure determination
1918 {"checkindex", &typeid(bool)}, // 8 - whether to use the FindIndex function or not
1919 {"tol", &typeid(double)}, // 9 - stopping criteria tol
1920 {"maxiter", &typeid(long)}, // 10 - stopping criteria : maximum number of iterations
1921 {"maxcputime", &typeid(double)}, // 11 - stopping criteria : maximum CPU time
1922 {"bfgs", &typeid(bool)}, // 12 - force the bfgs hessian approximation
1923 {"derivativetest", &typeid(string *)}, // 13 - call the derivative checker
1924 {"optfile", &typeid(string *)}, // 14 - set the ipopt option file name (default is ipopt.opt)
1925 {"printlevel", &typeid(long)}, // 15 - controls IPOPT print level output
1926 {"dth", &typeid(double)}, // 16 - perturbation for finite difference derivative test
1927 {"dttol", &typeid(double)}, // 17 - relative tolerence for the derivative test error detection
1928 {"fixedvar", &typeid(string *)}, // 18 - remove the equality simple bounds from problem
1929 {"warmstart", &typeid(bool)}, // 19 - do we initialize multipliers with given values
1930 {"uz", &typeid(KN_< double >)}, // 20 - simple upper bounds dual variable
1931 {"lz", &typeid(KN_< double >)}, // 21 - simple lower bounds dual variable
1932 {"muinit", &typeid(double)}, // 22 - barrier parameter initialization
1933 {"pivtol", &typeid(double)}, // 23 - pivot tolerance for the linear solver
1934 {"brf", &typeid(double)}, // 24 - bounds relax factor
1935 {"mustrategy", &typeid(string *)}, // 25 - strategy for barrier parameter update
1936 {"objvalue", &typeid(double *)}, // 26 - to get the last objective function value
1937 {"mumin", &typeid(double)}, // 27 - minimal value for the barrier parameter
1938 {"linesearch",
1939 &typeid(bool)} // 28 - use the line search or not (if no, the usual Newton step is kept)
1940 // {"osf", &typeid(double) }
1941 // //26 - objective function scalling factor
1942 };
1943
Load_Init()1944 static void Load_Init( ) {
1945 Global.Add("IPOPT", "(", new OptimIpopt(Case< no_assumption_f, no_assumption_g >( )));
1946 Global.Add("IPOPT", "(", new OptimIpopt(Case< no_assumption_f, without_constraints >( )));
1947 Global.Add("IPOPT", "(", new OptimIpopt(Case< no_assumption_f, mv_P1_g >( )));
1948 Global.Add("IPOPT", "(", new OptimIpopt(Case< no_assumption_f, linear_g >( )));
1949 Global.Add("IPOPT", "(", new OptimIpopt(Case< unavailable_hessian, no_assumption_g >( )));
1950 Global.Add("IPOPT", "(", new OptimIpopt(Case< unavailable_hessian, without_constraints >( )));
1951 Global.Add("IPOPT", "(", new OptimIpopt(Case< unavailable_hessian, mv_P1_g >( )));
1952 Global.Add("IPOPT", "(", new OptimIpopt(Case< unavailable_hessian, linear_g >( )));
1953 Global.Add("IPOPT", "(", new OptimIpopt(Case< mv_P2_f, no_assumption_g >( )));
1954 Global.Add("IPOPT", "(", new OptimIpopt(Case< mv_P2_f, without_constraints >( )));
1955 Global.Add("IPOPT", "(", new OptimIpopt(Case< mv_P2_f, mv_P1_g >( )));
1956 Global.Add("IPOPT", "(", new OptimIpopt(Case< mv_P2_f, linear_g >( )));
1957 Global.Add("IPOPT", "(", new OptimIpopt(Case< quadratic_f, no_assumption_g >( )));
1958 Global.Add("IPOPT", "(", new OptimIpopt(Case< quadratic_f, without_constraints >( )));
1959 Global.Add("IPOPT", "(", new OptimIpopt(Case< quadratic_f, mv_P1_g >( )));
1960 Global.Add("IPOPT", "(", new OptimIpopt(Case< quadratic_f, linear_g >( )));
1961 Global.Add("IPOPT", "(", new OptimIpopt(Case< linear_f, no_assumption_g >( )));
1962 Global.Add("IPOPT", "(", new OptimIpopt(Case< linear_f, without_constraints >( )));
1963 Global.Add("IPOPT", "(", new OptimIpopt(Case< linear_f, mv_P1_g >( )));
1964 Global.Add("IPOPT", "(", new OptimIpopt(Case< linear_f, linear_g >( )));
1965 }
1966
1967 /*****************************************************************************************************************************
1968 * Specialization of functions builders' constructor and operator()
1969 *****************************************************************************************************************************/
1970
1971 template<>
FitnessFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m)1972 FitnessFunctionDatas< no_assumption_f >::FitnessFunctionDatas(const basicAC_F0 &args,
1973 Expression const *nargs,
1974 const C_F0 &theparam,
1975 const C_F0 &objfact, const C_F0 &L_m)
1976 : GenericFitnessFunctionDatas( ) {
1977 const Polymorphic *opJ = dynamic_cast< const Polymorphic * >(args[0].LeftValue( )),
1978 *opdJ = dynamic_cast< const Polymorphic * >(args[1].LeftValue( )),
1979 *opH = dynamic_cast< const Polymorphic * >(args[2].LeftValue( ));
1980 ArrayOfaType hprototype2(atype< KN< R > * >( ), atype< double >( ), atype< KN< R > * >( )),
1981 hprototype1(atype< KN< R > * >( ));
1982
1983 JJ = to< R >(C_F0(opJ, "(", theparam));
1984 GradJ = to< Rn_ >(C_F0(opdJ, "(", theparam));
1985 if (opH->Find("(", hprototype2)) {
1986 CompletelyNonLinearConstraints = true;
1987 Hessian = to< Matrice_Creuse< R > * >(C_F0(opH, "(", theparam, objfact, L_m));
1988 } else if (opH->Find("(", hprototype1)) {
1989 CompletelyNonLinearConstraints =
1990 false; // When constraints are affine, lagrange multipliers are not used in the hessian,
1991 // obj_factor is also hidden to the user
1992 Hessian = to< Matrice_Creuse< R > * >(C_F0(opH, "(", theparam));
1993 } else {
1994 CompileError(
1995 "Error, wrong hessian function prototype. Must be either (real[int] &) or (real[int] "
1996 "&,real,real[int] &)");
1997 }
1998 }
1999
2000 template<>
operator ( )(Stack stack,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m,Expression const * nargs,ScalarFunc * & ffJ,VectorFunc * & ffdJ,SparseMatFunc * & ffH,bool warned) const2001 void FitnessFunctionDatas< no_assumption_f >::operator( )(Stack stack, const C_F0 &theparam,
2002 const C_F0 &objfact, const C_F0 &L_m,
2003 Expression const *nargs, ScalarFunc *&ffJ,
2004 VectorFunc *&ffdJ, SparseMatFunc *&ffH,
2005 bool warned) const {
2006 ffJ = new GeneralFunc< R >(stack, JJ, theparam);
2007 ffdJ = new GeneralFunc< Rn >(stack, GradJ, theparam);
2008 if (CompletelyNonLinearConstraints) {
2009 ffH = new GeneralSparseMatFunc(stack, Hessian, theparam, objfact, L_m);
2010 } else {
2011 ffH = new GeneralSparseMatFunc(stack, Hessian, theparam);
2012 }
2013 }
2014
2015 template<>
FitnessFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m)2016 FitnessFunctionDatas< P2_f >::FitnessFunctionDatas(const basicAC_F0 &args, Expression const *nargs,
2017 const C_F0 &theparam, const C_F0 &objfact,
2018 const C_F0 &L_m)
2019 : GenericFitnessFunctionDatas( ) {
2020 CompletelyNonLinearConstraints = false;
2021 const Polymorphic *opJ = dynamic_cast< const Polymorphic * >(args[0].LeftValue( )),
2022 *opdJ = dynamic_cast< const Polymorphic * >(args[1].LeftValue( ));
2023 JJ = to< R >(C_F0(opJ, "(", theparam));
2024 GradJ = to< Rn_ >(C_F0(opdJ, "(", theparam));
2025 Hessian = to< Matrice_Creuse< R > * >(args[2]);
2026 }
2027
2028 template<>
operator ( )(Stack stack,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m,Expression const * nargs,ScalarFunc * & ffJ,VectorFunc * & ffdJ,SparseMatFunc * & ffH,bool warned) const2029 void FitnessFunctionDatas< P2_f >::operator( )(Stack stack, const C_F0 &theparam,
2030 const C_F0 &objfact, const C_F0 &L_m,
2031 Expression const *nargs, ScalarFunc *&ffJ,
2032 VectorFunc *&ffdJ, SparseMatFunc *&ffH,
2033 bool warned) const {
2034 if (warned && nargs[5]) {
2035 cout << " ==> your lagrangian hessian is a constant matrix, so there is no need to specify "
2036 "its structure with ";
2037 cout << OptimIpopt::E_Ipopt::name_param[5].name << endl;
2038 cout << " since it is contained in the matrix object." << endl;
2039 }
2040
2041 ffJ = new GeneralFunc< R >(stack, JJ, theparam);
2042 ffdJ = new GeneralFunc< Rn >(stack, GradJ, theparam);
2043 ffH = new ConstantSparseMatFunc(stack, Hessian);
2044 }
2045
2046 template<>
FitnessFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m)2047 FitnessFunctionDatas< unavailable_hessian >::FitnessFunctionDatas(const basicAC_F0 &args,
2048 Expression const *nargs,
2049 const C_F0 &theparam,
2050 const C_F0 &objfact,
2051 const C_F0 &L_m)
2052 : GenericFitnessFunctionDatas( ) {
2053 CompletelyNonLinearConstraints = false;
2054 const Polymorphic *opJ = dynamic_cast< const Polymorphic * >(args[0].LeftValue( )),
2055 *opdJ = dynamic_cast< const Polymorphic * >(args[1].LeftValue( ));
2056 JJ = to< R >(C_F0(opJ, "(", theparam));
2057 GradJ = to< Rn_ >(C_F0(opdJ, "(", theparam));
2058 }
2059
2060 template<>
operator ( )(Stack stack,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m,Expression const * nargs,ScalarFunc * & ffJ,VectorFunc * & ffdJ,SparseMatFunc * & ffH,bool warned) const2061 void FitnessFunctionDatas< unavailable_hessian >::operator( )(
2062 Stack stack, const C_F0 &theparam, const C_F0 &objfact, const C_F0 &L_m, Expression const *nargs,
2063 ScalarFunc *&ffJ, VectorFunc *&ffdJ, SparseMatFunc *&ffH, bool warned) const {
2064 if (warned && nargs[5]) {
2065 cout << " ==> no hessian has been given, the LBFGS mode has been enabled, thus making ";
2066 cout << OptimIpopt::E_Ipopt::name_param[5].name << " useless. You may also" << endl
2067 << " have forgoten a function (IPOPT will certainly crash if so)." << endl;
2068 }
2069
2070 ffJ = new GeneralFunc< R >(stack, JJ, theparam);
2071 ffdJ = new GeneralFunc< Rn >(stack, GradJ, theparam);
2072 ffH = 0;
2073 }
2074
2075 template<>
FitnessFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m)2076 FitnessFunctionDatas< mv_P2_f >::FitnessFunctionDatas(const basicAC_F0 &args,
2077 Expression const *nargs, const C_F0 &theparam,
2078 const C_F0 &objfact, const C_F0 &L_m)
2079 : GenericFitnessFunctionDatas( ) {
2080 const E_Array *M_b = dynamic_cast< const E_Array * >(args[0].LeftValue( ));
2081
2082 if (M_b->nbitem( ) != 2) {
2083 CompileError(
2084 "\nSorry, we were expecting an array with two componants, either [M,b] or [b,M] for the "
2085 "affine constraints expression.");
2086 }
2087
2088 bool order = true;
2089 if (CheckMatrixVectorPair(M_b, order)) {
2090 Hessian = to< Matrice_Creuse< R > * >((*M_b)[order ? 0 : 1]);
2091 GradJ = to< Rn * >((*M_b)[order ? 1 : 0]); // This is gradJ evaluated on x=0
2092 }
2093 }
2094
2095 template<>
operator ( )(Stack stack,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m,Expression const * nargs,ScalarFunc * & ffJ,VectorFunc * & ffdJ,SparseMatFunc * & ffH,bool warned) const2096 void FitnessFunctionDatas< mv_P2_f >::operator( )(Stack stack, const C_F0 &theparam,
2097 const C_F0 &objfact, const C_F0 &L_m,
2098 Expression const *nargs, ScalarFunc *&ffJ,
2099 VectorFunc *&ffdJ, SparseMatFunc *&ffH,
2100 bool warned) const {
2101 if (warned && nargs[5]) {
2102 cout << " ==> your lagrangian hessian is a constant matrix, so there is no need to specify "
2103 "its structure with ";
2104 cout << OptimIpopt::E_Ipopt::name_param[5].name << endl;
2105 cout << " since it is contained in the matrix object." << endl;
2106 }
2107
2108 ffJ = new P2ScalarFunc(stack, Hessian, GradJ, true);
2109 ffdJ = new P1VectorFunc(stack, Hessian, GradJ, true);
2110 ffH = new ConstantSparseMatFunc(stack, Hessian);
2111 }
2112
2113 template<>
FitnessFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m)2114 FitnessFunctionDatas< quadratic_f >::FitnessFunctionDatas(const basicAC_F0 &args,
2115 Expression const *nargs,
2116 const C_F0 &theparam, const C_F0 &objfact,
2117 const C_F0 &L_m)
2118 : GenericFitnessFunctionDatas( ) {
2119 Hessian = to< Matrice_Creuse< R > * >(args[0]);
2120 }
2121
2122 template<>
operator ( )(Stack stack,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m,Expression const * nargs,ScalarFunc * & ffJ,VectorFunc * & ffdJ,SparseMatFunc * & ffH,bool warned) const2123 void FitnessFunctionDatas< quadratic_f >::operator( )(Stack stack, const C_F0 &theparam,
2124 const C_F0 &objfact, const C_F0 &L_m,
2125 Expression const *nargs, ScalarFunc *&ffJ,
2126 VectorFunc *&ffdJ, SparseMatFunc *&ffH,
2127 bool warned) const {
2128 if (warned && nargs[5]) {
2129 cout << " ==> your lagrangian hessian is a constant matrix, so there is no need to specify "
2130 "its structure with ";
2131 cout << OptimIpopt::E_Ipopt::name_param[5].name << endl;
2132 cout << " since it is contained in the matrix object." << endl;
2133 }
2134
2135 ffJ = new P2ScalarFunc(stack, Hessian, 0, true);
2136 ffdJ = new P1VectorFunc(stack, Hessian, 0, true);
2137 ffH = new ConstantSparseMatFunc(stack, Hessian);
2138 }
2139
2140 template<>
FitnessFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m)2141 FitnessFunctionDatas< linear_f >::FitnessFunctionDatas(const basicAC_F0 &args,
2142 Expression const *nargs,
2143 const C_F0 &theparam, const C_F0 &objfact,
2144 const C_F0 &L_m)
2145 : GenericFitnessFunctionDatas( ) {
2146 GradJ = to< Rn * >(args[0]);
2147 }
2148
2149 template<>
operator ( )(Stack stack,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & L_m,Expression const * nargs,ScalarFunc * & ffJ,VectorFunc * & ffdJ,SparseMatFunc * & ffH,bool warned) const2150 void FitnessFunctionDatas< linear_f >::operator( )(Stack stack, const C_F0 &theparam,
2151 const C_F0 &objfact, const C_F0 &L_m,
2152 Expression const *nargs, ScalarFunc *&ffJ,
2153 VectorFunc *&ffdJ, SparseMatFunc *&ffH,
2154 bool warned) const {
2155 if (warned && nargs[5]) {
2156 cout << " ==> your lagrangian hessian is a null matrix, so there is no need to specify its "
2157 "structure with ";
2158 cout << OptimIpopt::E_Ipopt::name_param[5].name << endl;
2159 cout << " since it is empty." << endl;
2160 }
2161
2162 ffJ = new P2ScalarFunc(stack, 0, GradJ);
2163 ffdJ = new P1VectorFunc(stack, 0, GradJ);
2164 ffH = 0;
2165 }
2166
2167 template<>
ConstraintFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam)2168 ConstraintFunctionDatas< without_constraints >::ConstraintFunctionDatas(const basicAC_F0 &args,
2169 Expression const *nargs,
2170 const C_F0 &theparam)
2171 : GenericConstraintFunctionDatas( ) {}
2172
2173 template<>
operator ( )(Stack stack,const C_F0 & theparam,Expression const * nargs,VectorFunc * & ffC,SparseMatFunc * & ffdC,bool warned) const2174 void ConstraintFunctionDatas< without_constraints >::operator( )(Stack stack, const C_F0 &theparam,
2175 Expression const *nargs,
2176 VectorFunc *&ffC,
2177 SparseMatFunc *&ffdC,
2178 bool warned) const {
2179 if (warned) {
2180 if (nargs[2] || nargs[3]) {
2181 cout << " ==> Some constraints bounds have been defined while no constraints function has "
2182 "been passed."
2183 << endl;
2184 }
2185
2186 if (nargs[4]) {
2187 cout << " ==> A structure has been provided for the constraints jacobian but there is no "
2188 "constraint function."
2189 << endl;
2190 }
2191
2192 if (nargs[6]) {
2193 cout << " ==> Unconstrained problem make the use of "
2194 << OptimIpopt::E_Ipopt::name_param[6].name
2195 << " pointless (see the documentation for more details)." << endl;
2196 }
2197 }
2198
2199 ffC = 0;
2200 ffdC = 0;
2201 }
2202
2203 template<>
ConstraintFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam)2204 ConstraintFunctionDatas< no_assumption_g >::ConstraintFunctionDatas(const basicAC_F0 &args,
2205 Expression const *nargs,
2206 const C_F0 &theparam)
2207 : GenericConstraintFunctionDatas( ) {
2208 int nbj = args.size( ) - 1;
2209 const Polymorphic *opG = dynamic_cast< const Polymorphic * >(args[nbj - 2].LeftValue( )),
2210 *opjG = dynamic_cast< const Polymorphic * >(args[nbj - 1].LeftValue( ));
2211
2212 Constraints = to< Rn_ >(C_F0(opG, "(", theparam));
2213 GradConstraints = to< Matrice_Creuse< R > * >(C_F0(opjG, "(", theparam));
2214 }
2215
2216 template<>
operator ( )(Stack stack,const C_F0 & theparam,Expression const * nargs,VectorFunc * & ffC,SparseMatFunc * & ffdC,bool) const2217 void ConstraintFunctionDatas< no_assumption_g >::operator( )(Stack stack, const C_F0 &theparam,
2218 Expression const *nargs,
2219 VectorFunc *&ffC, SparseMatFunc *&ffdC,
2220 bool) const {
2221 ffC = new GeneralFunc< Rn >(stack, Constraints, theparam);
2222 ffdC = new GeneralSparseMatFunc(stack, GradConstraints, theparam);
2223 }
2224
2225 template<>
ConstraintFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam)2226 ConstraintFunctionDatas< P1_g >::ConstraintFunctionDatas(const basicAC_F0 &args,
2227 Expression const *nargs,
2228 const C_F0 &theparam)
2229 : GenericConstraintFunctionDatas( ) {
2230 int nbj = args.size( ) - 1;
2231 const Polymorphic *opG = dynamic_cast< const Polymorphic * >(args[nbj - 2].LeftValue( ));
2232
2233 Constraints = to< Rn_ >(C_F0(opG, "(", theparam));
2234 GradConstraints = to< Matrice_Creuse< R > * >(args[nbj - 1]);
2235 }
2236
2237 template<>
operator ( )(Stack stack,const C_F0 & theparam,Expression const * nargs,VectorFunc * & ffC,SparseMatFunc * & ffdC,bool warned) const2238 void ConstraintFunctionDatas< P1_g >::operator( )(Stack stack, const C_F0 &theparam,
2239 Expression const *nargs, VectorFunc *&ffC,
2240 SparseMatFunc *&ffdC, bool warned) const {
2241 if (warned && nargs[4]) {
2242 cout << " ==> your constraints jacobian is a constant matrix, there is no need to specify its "
2243 "structure with "
2244 << OptimIpopt::E_Ipopt::name_param[4].name << endl;
2245 cout << " since it is contained in the matrix object." << endl;
2246 }
2247
2248 ffC = new GeneralFunc< Rn >(stack, Constraints, theparam);
2249 ffdC = new ConstantSparseMatFunc(stack, GradConstraints);
2250 }
2251
2252 template<>
ConstraintFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam)2253 ConstraintFunctionDatas< mv_P1_g >::ConstraintFunctionDatas(const basicAC_F0 &args,
2254 Expression const *nargs,
2255 const C_F0 &theparam)
2256 : GenericConstraintFunctionDatas( ) {
2257 int nbj = args.size( ) - 1;
2258 const E_Array *M_b = dynamic_cast< const E_Array * >(args[nbj - 1].LeftValue( ));
2259
2260 if (M_b->nbitem( ) != 2) {
2261 CompileError(
2262 "\nSorry, we were expecting an array with two componants, either [M,b] or [b,M] for the "
2263 "affine constraints expression.");
2264 }
2265
2266 bool order = true;
2267 if (CheckMatrixVectorPair(M_b, order)) {
2268 GradConstraints = to< Matrice_Creuse< R > * >((*M_b)[order ? 0 : 1]);
2269 Constraints = to< Rn * >((*M_b)[order ? 1 : 0]); // Constraint on x=0
2270 } else {
2271 CompileError(
2272 "\nWrong types in the constraints [matrix,vector] pair, expecting a sparse matrix and "
2273 "real[int].");
2274 }
2275 }
2276
2277 template<>
operator ( )(Stack stack,const C_F0 & theparam,Expression const * nargs,VectorFunc * & ffC,SparseMatFunc * & ffdC,bool warned) const2278 void ConstraintFunctionDatas< mv_P1_g >::operator( )(Stack stack, const C_F0 &theparam,
2279 Expression const *nargs, VectorFunc *&ffC,
2280 SparseMatFunc *&ffdC, bool warned) const {
2281 if (warned && nargs[4]) {
2282 cout << " ==> your constraints jacobian is a constant matrix, there is no need to specify its "
2283 "structure with "
2284 << OptimIpopt::E_Ipopt::name_param[4].name << endl;
2285 cout << " since it is contained in the matrix object." << endl;
2286 }
2287
2288 ffC = new P1VectorFunc(stack, GradConstraints, Constraints);
2289 ffdC = new ConstantSparseMatFunc(stack, GradConstraints);
2290 }
2291
2292 template<>
ConstraintFunctionDatas(const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam)2293 ConstraintFunctionDatas< linear_g >::ConstraintFunctionDatas(const basicAC_F0 &args,
2294 Expression const *nargs,
2295 const C_F0 &theparam)
2296 : GenericConstraintFunctionDatas( ) {
2297 int nbj = args.size( ) - 1;
2298
2299 GradConstraints = to< Matrice_Creuse< R > * >(args[nbj - 1]);
2300 }
2301
2302 template<>
operator ( )(Stack stack,const C_F0 & theparam,Expression const * nargs,VectorFunc * & ffC,SparseMatFunc * & ffdC,bool warned) const2303 void ConstraintFunctionDatas< linear_g >::operator( )(Stack stack, const C_F0 &theparam,
2304 Expression const *nargs, VectorFunc *&ffC,
2305 SparseMatFunc *&ffdC, bool warned) const {
2306 if (warned && nargs[4]) {
2307 cout << " ==> your constraints jacobian is a constant matrix, there is no need to specify its "
2308 "structure with "
2309 << OptimIpopt::E_Ipopt::name_param[4].name << endl;
2310 cout << " since it is contained in the matrix object." << endl;
2311 }
2312
2313 ffC = new P1VectorFunc(stack, GradConstraints, 0);
2314 ffdC = new ConstantSparseMatFunc(stack, GradConstraints);
2315 }
2316
New(AssumptionF AF,const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam,const C_F0 & objfact,const C_F0 & lm)2317 GenericFitnessFunctionDatas *GenericFitnessFunctionDatas::New(AssumptionF AF,
2318 const basicAC_F0 &args,
2319 Expression const *nargs,
2320 const C_F0 &theparam,
2321 const C_F0 &objfact, const C_F0 &lm) {
2322 switch (AF) {
2323 case no_assumption_f:
2324 return new FitnessFunctionDatas< no_assumption_f >(args, nargs, theparam, objfact, lm);
2325 break;
2326 case P2_f:
2327 return new FitnessFunctionDatas< P2_f >(args, nargs, theparam, objfact, lm);
2328 break;
2329 case unavailable_hessian:
2330 return new FitnessFunctionDatas< unavailable_hessian >(args, nargs, theparam, objfact, lm);
2331 break;
2332 case mv_P2_f:
2333 return new FitnessFunctionDatas< mv_P2_f >(args, nargs, theparam, objfact, lm);
2334 break;
2335 case quadratic_f:
2336 return new FitnessFunctionDatas< quadratic_f >(args, nargs, theparam, objfact, lm);
2337 break;
2338 case linear_f:
2339 return new FitnessFunctionDatas< linear_f >(args, nargs, theparam, objfact, lm);
2340 break;
2341 default:
2342 return 0;
2343 break;
2344 }
2345 }
2346
New(AssumptionG AG,const basicAC_F0 & args,Expression const * nargs,const C_F0 & theparam)2347 GenericConstraintFunctionDatas *GenericConstraintFunctionDatas::New(AssumptionG AG,
2348 const basicAC_F0 &args,
2349 Expression const *nargs,
2350 const C_F0 &theparam) {
2351 switch (AG) {
2352 case no_assumption_g:
2353 return new ConstraintFunctionDatas< no_assumption_g >(args, nargs, theparam);
2354 break;
2355 case without_constraints:
2356 return new ConstraintFunctionDatas< without_constraints >(args, nargs, theparam);
2357 break;
2358 case P1_g:
2359 return new ConstraintFunctionDatas< P1_g >(args, nargs, theparam);
2360 break;
2361 case mv_P1_g:
2362 return new ConstraintFunctionDatas< mv_P1_g >(args, nargs, theparam);
2363 break;
2364 case linear_g:
2365 return new ConstraintFunctionDatas< linear_g >(args, nargs, theparam);
2366 break;
2367 default:
2368 return 0;
2369 break;
2370 }
2371 }
2372
2373 LOADFUNC(Load_Init)
2374