00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "IpTNLP.hpp"
00011 #include "IpIpoptApplication.hpp"
00012
00013 #include "CouenneSparseMatrix.hpp"
00014 #include "CouenneProblem.hpp"
00015 #include "CouenneProblemElem.hpp"
00016 #include "CouenneExprVar.hpp"
00017 #include "CouenneExprJac.hpp"
00018 #include "CouenneExprHess.hpp"
00019 #include "CouenneTNLP.hpp"
00020
00021 #include <stdio.h>
00022
00023 #include "CoinHelperFunctions.hpp"
00024 #include "CoinFinite.hpp"
00025
00026
00027
00028 using namespace Ipopt;
00029 using namespace Couenne;
00030
00032 CouenneTNLP::CouenneTNLP ():
00033
00034 problem_ (NULL),
00035 sol0_ (NULL),
00036 sol_ (NULL),
00037 HLa_ (NULL),
00038
00039 optHessian_ (NULL),
00040 saveOptHessian_ (false) {}
00041
00042
00044 CouenneTNLP::~CouenneTNLP () {
00045
00046 if (sol0_) delete [] sol0_;
00047 if (sol_) delete [] sol_;
00048 if (HLa_) delete HLa_;
00049 if (optHessian_) delete optHessian_;
00050
00051 for (std::vector <std::pair <int, expression *> >::iterator i = gradient_. begin ();
00052 i != gradient_. end (); ++i)
00053 delete (*i). second;
00054 }
00055
00056
00058 CouenneTNLP::CouenneTNLP (CouenneProblem *p):
00059
00060 problem_ (p),
00061 sol0_ (NULL),
00062 sol_ (NULL),
00063 bestZ_ (COIN_DBL_MAX),
00064 Jac_ (p),
00065 HLa_ (new ExprHess (p)),
00066 optHessian_ (NULL),
00067 saveOptHessian_ (false) {
00068
00069 std::set <int> objDep;
00070
00071 expression *obj = problem_ -> Obj (0) -> Body ();
00072
00073
00074
00075 obj -> DepList (objDep, STOP_AT_AUX);
00076
00077 for (std::set <int>::iterator i = objDep.begin (); i != objDep. end (); ++i) {
00078
00079 expression *gradcomp = obj -> differentiate (*i);
00080 gradcomp -> realign (problem_);
00081 gradient_ . push_back (std::pair <int, expression *> (*i, gradcomp));
00082 }
00083
00084
00085
00086
00087
00088
00089 for (int i = 0; i < problem_ -> nCons (); i++) {
00090
00091 expression *e = problem_ -> Con (i) -> Body ();
00092
00093
00094
00095
00096 if (e -> Type () == AUX ||
00097 e -> Type () == VAR ||
00098 e -> Linearity () <= LINEAR)
00099 continue;
00100
00101
00102
00103
00104 e -> DepList (nonLinVars_, STOP_AT_AUX);
00105 }
00106
00107
00108
00109 for (int i = 0; i < problem_ -> nVars (); i++) {
00110
00111 exprVar *e = problem_ -> Var (i);
00112
00113 if ((e -> Type () != AUX) ||
00114 (e -> Multiplicity () <= 0) ||
00115 (e -> Linearity () <= LINEAR))
00116 continue;
00117
00118 e -> Image () -> DepList (nonLinVars_, STOP_AT_AUX);
00119 }
00120 }
00121
00122
00124 CouenneTNLP::CouenneTNLP (const CouenneTNLP &rhs)
00125 {operator= (rhs);}
00126
00127
00129 CouenneTNLP &CouenneTNLP::operator= (const CouenneTNLP &rhs) {
00130
00131 problem_ = rhs.problem_;
00132
00133 sol0_ = rhs.sol0_ && problem_ ? CoinCopyOfArray (rhs.sol0_, problem_ -> nVars ()) : NULL;
00134 sol_ = rhs.sol_ && problem_ ? CoinCopyOfArray (rhs.sol_, problem_ -> nVars ()) : NULL;
00135
00136 bestZ_ = rhs.bestZ_;
00137 gradient_ = rhs.gradient_;
00138 nonLinVars_ = rhs.nonLinVars_;
00139
00140 Jac_ = rhs.Jac_;
00141 HLa_ = rhs.HLa_ ? new ExprHess (*(rhs.HLa_)) : NULL;
00142
00143 optHessian_ = rhs.optHessian_ ? new CouenneSparseMatrix (*(rhs.optHessian_)) : NULL;
00144 saveOptHessian_ = rhs.saveOptHessian_;
00145
00146 return *this;
00147 }
00148
00150 CouenneTNLP *CouenneTNLP::clone ()
00151 {return new CouenneTNLP (*this);}
00152
00153
00154
00155
00156
00157
00158
00159 bool CouenneTNLP::get_nlp_info (Index& n,
00160 Index& m,
00161 Index& nnz_jac_g,
00162 Index& nnz_h_lag,
00163 IndexStyleEnum& index_style) {
00164 n = problem_ -> nVars ();
00165 m = Jac_. nRows ();
00166
00167 nnz_jac_g = Jac_ . nnz ();
00168 nnz_h_lag = HLa_ -> nnz ();
00169
00170 index_style = C_STYLE;
00171
00172 return true;
00173 }
00174
00175
00177 void CouenneTNLP::setInitSol (double *sol) {
00178
00179 if (sol) {
00180 if (!sol0_)
00181 sol0_ = new CouNumber [problem_ -> nVars ()];
00182 CoinCopyN (sol, problem_ -> nVars (), sol0_);
00183 }
00184 }
00185
00186
00187
00188
00189
00190
00191
00192 bool CouenneTNLP::get_bounds_info (Index n, Number* x_l, Number* x_u,
00193 Index m, Number* g_l, Number* g_u) {
00194
00195
00196
00197 #ifdef DEBUG
00198 printf ("get_bounds_info on %d cons, %d vars\n", m, n);
00199 #endif
00200
00201 for (int i = 0; i < problem_ -> nCons (); i++) {
00202
00203 CouenneConstraint *c = problem_ -> Con (i);
00204
00205 if (c -> Body () -> Type () == AUX ||
00206 c -> Body () -> Type () == VAR)
00207 continue;
00208
00209 CouNumber
00210 clb = (*c -> Lb ()) (),
00211 cub = (*c -> Ub ()) ();
00212
00213
00214 if (clb <= cub) {*g_l++ = clb; *g_u++ = cub;}
00215 else {*g_l++ = cub; *g_u++ = clb;}
00216 }
00217
00218
00219
00220 for (int i = 0; i < problem_ -> nVars (); i++) {
00221
00222 exprVar *e = problem_ -> Var (i);
00223
00224 CouNumber
00225 lb = e -> lb (),
00226 ub = e -> ub ();
00227
00228
00229 if (lb <= ub) {*x_l++ = lb; *x_u++ = ub;}
00230 else {*x_l++ = ub; *x_u++ = lb;}
00231
00232 if ((e -> Type () != AUX) ||
00233 (e -> Multiplicity () <= 0))
00234 continue;
00235
00236 *g_l = (e -> sign () != expression::AUX_GEQ) ? 0. : -COIN_DBL_MAX;
00237 *g_u = (e -> sign () != expression::AUX_LEQ) ? 0. : COIN_DBL_MAX;
00238
00239 ++g_l;
00240 ++g_u;
00241 }
00242
00243 return true;
00244 }
00245
00246
00247
00248
00249
00250
00251 bool CouenneTNLP::get_variables_linearity (Index n, Ipopt::TNLP::LinearityType* var_types) {
00252
00253 CoinFillN (var_types, n, Ipopt::TNLP::LINEAR);
00254
00255 for (std::set <int>:: iterator i = nonLinVars_. begin (); i != nonLinVars_. end (); ++i)
00256 var_types [*i] = Ipopt::TNLP::NON_LINEAR;
00257
00258 return true;
00259 }
00260
00261
00262
00263
00264 bool CouenneTNLP::get_constraints_linearity (Index m, Ipopt::TNLP::LinearityType* const_types) {
00265
00266
00267
00268 for (int i = 0; i < problem_ -> nCons (); i++) {
00269
00270 expression *b = problem_ -> Con (i) -> Body ();
00271
00272 if (b -> Type () == AUX ||
00273 b -> Type () == VAR)
00274 continue;
00275
00276 *const_types++ =
00277 (b -> Linearity () > LINEAR) ?
00278 Ipopt::TNLP::NON_LINEAR :
00279 Ipopt::TNLP::LINEAR;
00280 }
00281
00282
00283
00284 for (int i = 0; i < problem_ -> nVars (); i++) {
00285
00286 exprVar *e = problem_ -> Var (i);
00287
00288 if ((e -> Type () != AUX) ||
00289 (e -> Multiplicity () <= 0))
00290 continue;
00291
00292 *const_types++ =
00293 (e -> Image () -> Linearity () > LINEAR) ?
00294 Ipopt::TNLP::NON_LINEAR :
00295 Ipopt::TNLP::LINEAR;
00296 }
00297
00298 return true;
00299 }
00300
00301
00302
00303
00304
00305
00306
00307
00308 bool CouenneTNLP::get_starting_point (Index n,
00309 bool init_x, Number* x,
00310 bool init_z, Number* z_L, Number* z_U,
00311 Index m,
00312 bool init_lambda, Number* lambda) {
00313 if (init_x)
00314 CoinCopyN (sol0_, n, x);
00315
00316 assert (!init_z);
00317 assert (!init_lambda);
00318
00319 return true;
00320 }
00321
00322
00323
00324 bool CouenneTNLP::eval_f (Index n, const Number* x, bool new_x,
00325 Number& obj_value) {
00326 if (new_x)
00327 CoinCopyN (x, n, problem_ -> X ());
00328
00329
00330 obj_value = (*(problem_ -> Obj (0) -> Body ())) ();
00331 return true;
00332 }
00333
00334
00335
00336
00337 bool CouenneTNLP::eval_grad_f (Index n, const Number* x, bool new_x,
00338 Number* grad_f) {
00339
00340 #ifdef DEBUG
00341 printf ("eval_grad_f: [");
00342 for (int i=0; i<n; i++)
00343 printf ("%.2g ", x [i]);
00344 printf ("] --> [");
00345 #endif
00346
00347 if (new_x)
00348 CoinCopyN (x, n, problem_ -> X ());
00349
00350
00351 CoinFillN (grad_f, n, 0.);
00352
00353 for (std::vector <std::pair <int, expression *> >::iterator i = gradient_. begin ();
00354 i != gradient_. end (); ++i)
00355 grad_f [i -> first] = (*(i -> second)) ();
00356
00357 #ifdef DEBUG
00358 for (int i=0; i<n; i++)
00359 printf ("%.2g ", grad_f [i]);
00360 printf ("]\n");
00361 #endif
00362
00363 return true;
00364 }
00365
00366
00367
00368 bool CouenneTNLP::eval_g (Index n, const Number* x, bool new_x,
00369 Index m, Number* g) {
00370
00371 if (new_x)
00372 CoinCopyN (x, n, problem_ -> X ());
00373
00374
00375 #ifdef DEBUG
00376 if (x) {
00377 printf ("eval_g: [");
00378 for (int i=0; i<n; i++)
00379 printf ("%.2g ", x [i]);
00380 printf ("] --> [");
00381 }
00382 #endif
00383
00384 int nEntries = 0;
00385
00386 for (int i = 0; i < problem_ -> nCons (); i++) {
00387
00388 expression *b = problem_ -> Con (i) -> Body ();
00389
00390 if (b -> Type () == AUX ||
00391 b -> Type () == VAR)
00392 continue;
00393
00394 nEntries ++;
00395
00396 *g++ = (*b) ();
00397 }
00398
00399
00400
00401 assert (n == problem_ -> nVars ());
00402
00403 for (int i = 0; i < problem_ -> nVars (); i++) {
00404
00405 exprVar *e = problem_ -> Var (i);
00406
00407 if ((e -> Type () != AUX) ||
00408 (e -> Multiplicity () <= 0))
00409 continue;
00410
00411 *g++ = (*(e -> Image ())) () - (*e) ();
00412
00413 nEntries ++;
00414 }
00415
00416 #ifdef DEBUG
00417 if (x) {
00418 for (int i=0; i<nEntries; i++)
00419 printf ("%.2g ", *(g - nEntries + i));
00420 printf ("]\n");
00421 }
00422 #endif
00423
00424 return true;
00425 }
00426
00427
00428
00429
00430
00431
00432
00433 bool CouenneTNLP::eval_jac_g (Index n, const Number* x, bool new_x,
00434 Index m, Index nele_jac, Index* iRow,
00435 Index *jCol, Number* values) {
00436 if (new_x)
00437 CoinCopyN (x, n, problem_ -> X ());
00438
00439
00440 #ifdef DEBUG
00441 if (x) {
00442 printf ("eval_jac_g: ["); fflush (stdout);
00443 for (int i=0; i<n; i++)
00444 {printf ("%.2g ", x [i]); fflush (stdout);}
00445 printf ("] --> ["); fflush (stdout);
00446 }
00447 #endif
00448
00449 if (values == NULL &&
00450 iRow != NULL &&
00451 jCol != NULL) {
00452
00453
00454
00455
00456 CoinCopyN (Jac_.iRow (), nele_jac, iRow);
00457 CoinCopyN (Jac_.jCol (), nele_jac, jCol);
00458
00459 } else {
00460
00461
00462
00463
00464 register expression **e = Jac_. expr ();
00465
00466 for (register int i=nele_jac; i--;)
00467 *values++ = (**(e++)) ();
00468 }
00469
00470 #ifdef DEBUG
00471 if (values) {
00472 for (int i=0; i<nele_jac; i++)
00473 {printf ("%.2g ", *(values - nele_jac + i)); fflush (stdout);}
00474 printf ("]\n");
00475 } else printf ("empty\n");
00476 #endif
00477
00478 return true;
00479 }
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493 bool CouenneTNLP::eval_h (Index n, const Number* x, bool new_x, Number obj_factor,
00494 Index m, const Number* lambda, bool new_lambda,
00495 Index nele_hess,
00496 Index* iRow, Index* jCol, Number* values) {
00497
00498 if (new_x)
00499 CoinCopyN (x, n, problem_ -> X ());
00500
00501
00502 #ifdef DEBUG
00503 if (x) {
00504 printf ("eval_h: ["); fflush (stdout);
00505 for (int i=0; i<n; i++)
00506 {printf ("%.2g ", x [i]); fflush (stdout);}
00507 printf ("], lambda: ["); fflush (stdout);
00508 for (int i=0; i<m; i++)
00509 {printf ("%.2g ", lambda [i]); fflush (stdout);}
00510 printf ("] --> ["); fflush (stdout);
00511 }
00512 #endif
00513
00514 if (values == NULL &&
00515 iRow != NULL &&
00516 jCol != NULL) {
00517
00519
00520 CoinCopyN (HLa_ -> iRow (), nele_hess, iRow);
00521 CoinCopyN (HLa_ -> jCol (), nele_hess, jCol);
00522
00523 } else {
00524
00527
00528 CoinZeroN (values, nele_hess);
00529
00530 for (int i=0; i<nele_hess; i++, values++) {
00531
00532 int
00533 numL = HLa_ -> numL () [i],
00534 *lamI = HLa_ -> lamI () [i];
00535
00536 expression
00537 **expr = HLa_ -> expr () [i];
00538
00539 #ifdef DEBUG
00540 printf ("[%d %d] %d lambdas: ", HLa_ -> iRow () [i], HLa_ -> jCol () [i], numL); fflush (stdout);
00541 for (int k=0; k<numL; k++) {
00542 printf ("%d ", lamI [k]);
00543 fflush (stdout);
00544 expr [k] -> print ();
00545 printf ("\n");
00546 }
00547 #endif
00548
00549
00550
00551 if (0 == *lamI) {*values += obj_factor * (*(*expr++)) (); --numL; ++lamI;}
00552 while (numL--) *values += lambda [*lamI++ - 1] * (*(*expr++)) ();
00553 }
00554 }
00555
00556
00557 #ifdef DEBUG
00558 if (values) {
00559 for (int i=0; i<nele_hess; i++)
00560 {printf ("%.2g ", *(values - nele_hess + i)); fflush (stdout);}
00561 printf ("]\n");
00562 } else printf ("empty\n");
00563 #endif
00564
00565 return true;
00566 }
00567
00568
00569
00570 void CouenneTNLP::setObjective (expression *newObj) {
00571
00572 if (HLa_)
00573 delete HLa_;
00574
00575
00576
00577 HLa_ = new ExprHess (problem_);
00578
00579 std::set <int> objDep;
00580
00581
00582
00583 newObj -> DepList (objDep, STOP_AT_AUX);
00584
00585 for (std::vector <std::pair <int, expression *> >::iterator i = gradient_. begin ();
00586 i != gradient_. end (); ++i)
00587 delete (*i). second;
00588
00589 gradient_ . erase (gradient_ . begin (), gradient_ . end ());
00590
00591 for (std::set <int>::iterator i = objDep.begin (); i != objDep. end (); ++i) {
00592
00593 expression
00594 *gradcomp = newObj -> differentiate (*i),
00595 *gsimp = gradcomp -> simplify ();
00596
00597 if (gsimp) {
00598 delete gradcomp;
00599 gradcomp = gsimp;
00600 }
00601
00602 gradcomp -> realign (problem_);
00603 gradient_ . push_back (std::pair <int, expression *> (*i, gradcomp));
00604 }
00605 }
00606
00607
00608
00609 void CouenneTNLP::finalize_solution (SolverReturn status,
00610 Index n, const Number* x, const Number* z_L, const Number* z_U,
00611 Index m, const Number* g, const Number* lambda,
00612 Number obj_value,
00613 const IpoptData* ip_data,
00614 IpoptCalculatedQuantities* ip_cq) {
00615
00616
00617
00618 bestZ_ = obj_value;
00619
00620 if (sol_) CoinCopyN (x, n, sol_);
00621 else sol_ = CoinCopyOfArray (x, n);
00622
00623
00624
00625
00626 if (saveOptHessian_) {
00627
00628 if (!optHessian_)
00629 optHessian_ = new CouenneSparseMatrix;
00630
00631 problem_ -> domain () -> push (n, x, NULL, NULL);
00632
00633 int nnz = HLa_ -> nnz ();
00634
00635
00636
00637 double *&optHessianVal = optHessian_ -> val ();
00638 int *&optHessianRow = optHessian_ -> row ();
00639 int *&optHessianCol = optHessian_ -> col ();
00640
00641 int &optHessianNum = optHessian_ -> num ();
00642
00643 optHessianVal = (double *) realloc (optHessianVal, nnz * sizeof (double));
00644 optHessianRow = (int *) realloc (optHessianRow, nnz * sizeof (int));
00645 optHessianCol = (int *) realloc (optHessianCol, nnz * sizeof (int));
00646
00647 optHessianNum = 0;
00648
00649 for (int i=0; i < HLa_ -> nnz (); ++i) {
00650
00651 double hessMember = 0.;
00652 expression **elist = HLa_ -> expr () [i];
00653
00654 for (int j=0; j < HLa_ -> numL () [i]; ++j) {
00655
00656 int indLam = HLa_ -> lamI () [i][j];
00657
00658 hessMember += (indLam == 0) ?
00659 (*(elist [j])) () :
00660 (*(elist [j])) () * lambda [indLam-1];
00661 }
00662
00663 if (hessMember != 0.) {
00664
00665
00666
00667
00668
00669 optHessianVal [optHessianNum] = hessMember;
00670 optHessianRow [optHessianNum] = HLa_ -> iRow () [i];
00671 optHessianCol [optHessianNum++] = HLa_ -> jCol () [i];
00672 }
00673 }
00674
00675 optHessianVal = (double *) realloc (optHessianVal, optHessianNum * sizeof (double));
00676 optHessianRow = (int *) realloc (optHessianRow, optHessianNum * sizeof (int));
00677 optHessianCol = (int *) realloc (optHessianCol, optHessianNum * sizeof (int));
00678
00679 problem_ -> domain () -> pop ();
00680 }
00681 }
00682
00683
00684
00685
00686 bool CouenneTNLP::intermediate_callback (AlgorithmMode mode,
00687 Index iter, Number obj_value,
00688 Number inf_pr, Number inf_du,
00689 Number mu, Number d_norm,
00690 Number regularization_size,
00691 Number alpha_du, Number alpha_pr,
00692 Index ls_trials,
00693 const IpoptData* ip_data,
00694 IpoptCalculatedQuantities* ip_cq) {
00695
00696
00697 return true;
00698 }
00699
00700
00701
00702
00703
00704
00705
00706
00707 Index CouenneTNLP::get_number_of_nonlinear_variables ()
00708 {return nonLinVars_. size ();}
00709
00710
00711
00712
00713
00714
00715
00716
00717 bool CouenneTNLP::get_list_of_nonlinear_variables (Index num_nonlin_vars,
00718 Index* pos_nonlin_vars) {
00719
00720 for (std::set <int>:: iterator i = nonLinVars_. begin (); i != nonLinVars_. end (); ++i)
00721 *pos_nonlin_vars++ = *i;
00722
00723 return true;
00724 }