checkNLP.cpp
Go to the documentation of this file.
1 /* $Id: checkNLP.cpp 967 2013-06-24 10:13:05Z fmargot $
2  *
3  * Name: checkNLP.cpp
4  * Author: Pietro Belotti
5  * Francois Margot
6  * Purpose: check NLP feasibility of incumbent integer solution
7  *
8  * (C) Carnegie-Mellon University, 2006-11.
9  * This file is licensed under the Eclipse Public License (EPL)
10  */
11 
12 #include "CouenneProblem.hpp"
13 #include "CouenneProblemElem.hpp"
14 #include "CouenneExprVar.hpp"
15 #include "CoinHelperFunctions.hpp"
16 #include "CoinFinite.hpp"
17 
18 #include "CouenneRecordBestSol.hpp"
19 
20 using namespace Couenne;
21 
22 // check if solution is MINLP feasible
23 bool CouenneProblem::checkNLP (const double *solution, double &obj, bool recompute) const {
24 
25  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
26 
27  printf ("Checking solution: %.12e (", obj);
28 
29  for (int i=0; i<nOrigVars_ - ndefined_; i++)
30  printf ("%g ", solution [i]);
31  printf (")\n");
32  }
33 
34  // pre-check on original variables --- this is done after every LP,
35  // and should be efficient
36  for (register int i=0; i < nOrigVars_ - ndefined_; i++) {
37 
38  if (variables_ [i] -> Multiplicity () <= 0)
39  continue;
40 
41  CouNumber val = solution [i];
42 
43  // check (original and auxiliary) variables' integrality
44 
45  exprVar *v = variables_ [i];
46 
47  if ((v -> Type () == VAR) &&
48  (v -> isDefinedInteger ()) &&
49  (v -> Multiplicity () > 0) &&
50  (fabs (val - COUENNE_round (val)) > feas_tolerance_)) {
51 
52  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
53  "checkNLP: integrality %d violated: %.6f [%g,%g]\n",
54  i, val, domain_.lb (i), domain_.ub (i));
55 
56  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "Done: (0)\n");
57 
58  return false;
59  }
60  }
61 
62  const int infeasible = 1;
63  const int wrong_obj = 2;
64 
65  // copy solution, evaluate the corresponding aux, and then replace
66  // the original variables again for checking
67  CouNumber *sol = new CouNumber [nVars ()];
68  CoinZeroN (sol + nOrigVars_ - ndefined_, nVars() - (nOrigVars_ - ndefined_));
69  CoinCopyN (solution, nOrigVars_ - ndefined_, sol);
70  getAuxs (sol);
71  CoinCopyN (solution, nOrigVars_ - ndefined_, sol);
72 
73  // install NL solution candidate in evaluation structure
74  domain_.push (nVars (), sol, domain_.lb (), domain_.ub (), false);
75 
76  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
77  printf (" checkNLP: %d vars -------------------\n ", domain_.current () -> Dimension ());
78  for (int i=0; i<domain_.current () -> Dimension (); i++) {
79  if (i && !(i%5)) printf ("\n ");
80  printf ("%4d %16g [%16e %16e] ", i, domain_.x (i), domain_.lb (i), domain_.ub (i));
81  }
82  }
83 
84  expression *objBody = Obj (0) -> Body ();
85 
86  // BUG: if Ipopt solution violates bounds of original variables and
87  // objective depends on originals, we may have a "computed object"
88  // out of bounds
89 
90  //CouNumber realobj = (*(objBody -> Image () ? objBody -> Image () : objBody)) ();
91  CouNumber realobj = obj;
92 
93  if (objBody)
94  realobj =
95  (objBody -> Index () >= 0) ?
96  sol [objBody -> Index ()] :
97  (*(objBody -> Image () ? objBody -> Image () : objBody)) ();
98 
99  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, " Objective: %.12e %.12e %.12e\n",
100  realobj, objBody -> Index () >= 0 ? sol [objBody -> Index ()] : objBody -> Value (),
101  (*(objBody -> Image () ? objBody -> Image () : objBody)) ());
102 
103  bool retval = true;
104 
105  try {
106 
107  // check if objective corresponds
108 
109  if (fabs (realobj - obj) / (1. + fabs (realobj)) > feas_tolerance_) {
110 
111  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
112  " checkNLP, false objective: computed %g != %g xQ (diff. %g)\n",
113  realobj, obj, realobj - obj);
114 
115  if (!recompute)
116  throw wrong_obj;
117  }
118 
119  if (recompute)
120  obj = realobj;
121 
122  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, " recomputed: %.12e\n", obj);
123 
124  for (int i=0; i < nOrigVars_ - ndefined_; i++) {
125 
126  if (variables_ [i] -> Multiplicity () <= 0)
127  continue;
128 
129  CouNumber val = domain_.x (i);
130 
131  // check bounds
132 
133  if ((val > domain_.ub (i) + feas_tolerance_) ||
134  (val < domain_.lb (i) - feas_tolerance_)) {
135 
136  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
137  " checkNLP: variable %d out of bounds: %.6f [%g,%g] (diff %g)\n",
138  i, val, domain_.lb (i), domain_.ub (i),
139  CoinMax (fabs (val - domain_.lb (i)),
140  fabs (val - domain_.ub (i))));
141  throw infeasible;
142  }
143 
144  // check (original and auxiliary) variables' integrality
145 
146  if (variables_ [i] -> isDefinedInteger () &&
147  (fabs (val - COUENNE_round (val)) > feas_tolerance_)) {
148 
149  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
150  " checkNLP: integrality %d violated: %.6f [%g,%g]\n",
151  i, val, domain_.lb (i), domain_.ub (i));
152 
153  throw infeasible;
154  }
155  }
156 
157  // check ALL auxs
158 
159  for (int i=0; i < nVars (); i++) {
160 
161  exprVar *v = variables_ [i];
162 
163  if ((v -> Type () != AUX) ||
164  (v -> Multiplicity () <= 0))
165  continue;
166 
167  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
168  printf (" "); v -> print ();
169  CouNumber
170  val = (*(v)) (),
171  img = (*(v -> Image ())) (),
172  diff = fabs (val - img);
173  printf (": val = %15g; img = %-15g ", val, img);
174  if (diff > 1e-9)
175  printf ("[diff %12e] ", diff);
176  //for (int j=0; j<nVars (); j++) printf ("%.12e ", (*(variables_ [j])) ());
177  v -> Image () -> print ();
178  printf ("\n");
179  }
180 
181  // check if auxiliary has zero infeasibility
182 
183  // same as in CouenneObject::checkInfeasibility -- main difference is use of gradientNorm()
184 
185  double
186  vval = (*v) (),
187  fval = (*(v -> Image ())) (),
188  denom = CoinMax (1., v -> Image () -> gradientNorm (X ()));
189 
190  // check if fval is a number (happens with e.g. w13 = w12/w5 and w5=0, see test/harker.nl)
191  if (CoinIsnan (fval)) {
192  fval = vval + 1.;
193  denom = 1.;
194  }
195 
196  if (fabs (fval) > COUENNE_INFINITY)
197  fval = COUENNE_INFINITY;
198 
199  double
200  delta =
201  ((v -> sign () == expression::AUX_GEQ) && (vval >= fval)) ? 0. :
202  ((v -> sign () == expression::AUX_LEQ) && (vval <= fval)) ? 0. : fabs (vval - fval),
203 
204  ratio = (CoinMax (1., fabs (vval)) /
205  CoinMax (1., fabs (fval)));
206 
207  // printf ("checkinf --> v=%e f=%e den=%e ret=%d ratio=%e delta=%e, delta/denom=%e, thres=%e [",
208  // vval, fval, denom, retval, ratio, delta, delta/denom, CoinMin (COUENNE_EPS, feas_tolerance_));
209  // v -> print ();
210  // printf (" %c= ", v -> sign () == expression::AUX_LEQ ? '<' :
211  // v -> sign () == expression::AUX_GEQ ? '>' : ':');
212  // v -> Image () -> print ();
213  // printf ("]\n");
214 
215  if ((delta > 0.) &&
216  (((ratio > 2.) || // check delta > 0 to take into account semi-auxs
217  (ratio < .5)) ||
218  ((delta /= denom) > CoinMin (COUENNE_EPS, feas_tolerance_)))) {
219 
220  Jnlst () -> Printf (Ipopt::J_MOREVECTOR, J_PROBLEM,
221  " checkNLP: auxiliary %d violates tolerance %g by %g/%g = %g\n",
222  i, CoinMin (COUENNE_EPS, feas_tolerance_), delta*denom, denom, delta);
223 
224  throw infeasible;
225  }
226  }
227 
228  // check constraints
229 
230  for (int i=0; i < nCons (); i++) {
231 
232  CouenneConstraint *c = Con (i);
233 
234  CouNumber
235  body = (*(c -> Body ())) (),
236  lhs = (*(c -> Lb ())) (),
237  rhs = (*(c -> Ub ())) ();
238 
239  if (((rhs < COUENNE_INFINITY) && (body > rhs + feas_tolerance_ * (1. + CoinMax (fabs (body), fabs (rhs))))) ||
240  ((lhs > -COUENNE_INFINITY) && (body < lhs - feas_tolerance_ * (1. + CoinMax (fabs (body), fabs (lhs)))))) {
241 
242  if (Jnlst () -> ProduceOutput (Ipopt::J_MOREVECTOR, J_PROBLEM)) {
243 
244  printf (" checkNLP: constraint %d violated (lhs=%+e body=%+e rhs=%+e, violation %g): ",
245  i, lhs, body, rhs, CoinMax (lhs-body, body-rhs));
246 
247  c -> print ();
248  }
249 
250  throw infeasible;
251  }
252  }
253  }
254 
255  catch (int exception) {
256 
257  switch (exception) {
258 
259  case wrong_obj:
260  retval = false;
261  break;
262 
263  case infeasible:
264  default:
265  retval = false;
266  break;
267  }
268  }
269 
270  delete [] sol;
271  domain_.pop ();
272 
273  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "Done: %d\n", retval);
274 
275  return retval;
276 }
277 
278 /************************************************************************/
279 // Recompute objective value for sol
280 double CouenneProblem::checkObj(const CouNumber *sol, const double &precision)
281  const {
282 
283  expression *objBody = Obj(0)->Body();
284 
285  // BUG: if Ipopt couSol violates bounds of original variables and
286  // objective depends on originals, we may have a "computed object"
287  // out of bounds
288 
289  //CouNumber realobj = (*(objBody -> Image () ? objBody -> Image () : objBody)) ();
290  CouNumber realObj = 0;
291 
292  if (objBody) {
293  realObj =
294  (objBody ->Index() >= 0) ?
295  sol[objBody->Index()] :
296  (*(objBody->Image() ? objBody->Image() : objBody)) ();
297 
298  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM,
299  "%.12e %.12e %.12e ------------------------------\n",
300  realObj, objBody -> Index () >= 0 ? sol[objBody -> Index ()] : 0.,
301  (*(objBody -> Image () ? objBody -> Image () : objBody)) ());
302  }
303  else {
304  printf("### ERROR: CouenneProblem::checkObj(): no objective body\n");
305  exit(1);
306  }
307  return realObj;
308 }
309 
310 /************************************************************************/
311 // check integrality of original vars in sol; return true if all
312 // original integer vars are within precision of an integer value
313 bool CouenneProblem::checkInt(const CouNumber *sol,
314  const int from, const int upto,
315  const std::vector<int> listInt,
316  const bool origVarOnly,
317  const bool stopAtFirstViol,
318  const double precision, double &maxViol) const {
319 
320  bool isFeas = true;
321 
322  for(unsigned int i=0; i<listInt.size(); i++) {
323 
324  int ind = listInt[i];
325 
326  if((ind < from) || (variables_ [ind] -> Multiplicity () <= 0))
327  continue;
328 
329  if(ind >= upto)
330  break;
331 
332  CouNumber val = sol[ind];
333  exprVar *v = variables_ [ind];
334 
335  if ((!origVarOnly) || (v -> Type () == VAR)) {
336 
337  double viol = fabs (val - COUENNE_round (val));
338  if (viol > maxViol) maxViol = viol;
339  if (viol > precision) {
340 
341  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM,
342  "checkInt(): integrality %d violated: %.6f [%g,%g]: integer distance %e > %e (by %e)\n",
343  i, val, domain_.lb (i), domain_.ub (i),
344  fabs (val - COUENNE_round (val)), feas_tolerance_,
345  fabs (val - COUENNE_round (val)) - feas_tolerance_);
346 
347  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkInt(): integrality %d violated: %.6f [%g,%g]\n", ind, val, domain_.lb (ind), domain_.ub (ind));
348 
349  isFeas = false;
350 
351  if (stopAtFirstViol)
352  break;
353  }
354  }
355  }
356  return(isFeas);
357 }
358 
359 /************************************************************************/
360 // Check bounds; returns true iff feasible for given precision
361 bool CouenneProblem::checkBounds(const CouNumber *sol,
362  const bool stopAtFirstViol,
363  const double precision, double &maxViol) const {
364 
365  bool isFeas = true;
366  for(int i=0; i<nOrigVars_ - ndefined_; i++) {
367 
368  if (variables_[i]-> Multiplicity () <= 0)
369  continue;
370 
371  CouNumber val = domain_.x (i);
372  double viol = 0;
373  double violUb = val - domain_.ub (i);
374  double violLb = domain_.lb (i) - val;
375 
376  if (viol < violUb) viol = violUb;
377  if (viol < violLb) viol = violLb;
378 
379  maxViol = (maxViol > viol ? maxViol : viol);
380 
381  if (viol > precision) {
382 
383  Jnlst()->Printf(Ipopt::J_MOREVECTOR, J_PROBLEM, "checkBounds(): variable %d out of bounds: %.6f [%g,%g] (diff %g)\n",
384  i, val, domain_.lb (i), domain_.ub (i),
385  CoinMax (fabs (val - domain_.lb (i)),
386  fabs (val - domain_.ub (i))));
387 
388  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkBounds: variable %d out of bounds: %.6f [%g,%g] (diff %g)\n",
389  i, val, domain_.lb (i), domain_.ub (i),
390  CoinMax (fabs (val - domain_.lb (i)),
391  fabs (val - domain_.ub (i))));
392 
393  isFeas = false;
394  if(stopAtFirstViol)
395  break;
396  }
397  }
398  return isFeas;
399 }
400 
401 /************************************************************************/
402 // returns true iff value of all auxiliaries are within bounds
403 bool CouenneProblem::checkAux(const CouNumber *sol,
404  const bool stopAtFirstViol,
405  const double precision, double &maxViol) const {
406 
407  bool isFeas = true;
408  for (int i=0; i<nVars(); i++) {
409 
410  exprVar *v = variables_ [i];
411 
412  if ((v -> Type () != AUX) ||
413  (v -> Multiplicity () <= 0))
414  continue;
415 
416  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
417 
418  printf ("before check\n");
419 
420  double
421  vdb = (*(variables_ [i])) (),
422  fdb = (*(variables_ [i] -> Image ())) ();
423 
424  double
425  del =
426  ((v -> sign () == expression::AUX_GEQ) && (vdb >= fdb)) ? 0. :
427  ((v -> sign () == expression::AUX_LEQ) && (vdb <= fdb)) ? 0. :
428  fabs (vdb - fdb);
429 
430  printf ("[%g,%g]\n", vdb, fdb);
431 
432  printf ("checking aux -- %+.12e = %+.12e [%+.12e] ", vdb, fdb, del);
433  variables_ [i] -> print ();
434  if(v->sign()== expression::AUX_GEQ) printf(" >= ");
435  if(v->sign()== expression::AUX_LEQ) printf(" <= ");
436  if(v->sign()== expression::AUX_EQ) printf(" := ");
437  variables_ [i] -> Image () -> print (); printf ("\n");
438  }
439 
440  double
441  vval = (*v) (),
442  fval = (*(v -> Image ())) (),
443  denom = CoinMax (1., v -> Image () -> gradientNorm (X ()));
444 
445  // check if fval is a number (happens with e.g. w13 = w12/w5 and w5=0, see test/harker.nl)
446  if (CoinIsnan (fval)) {
447  fval = vval + 1.;
448  denom = 1.;
449  }
450 
451  if (fabs (fval) > COUENNE_INFINITY)
452  fval = COUENNE_INFINITY;
453 
454  double
455  delta =
456  ((v -> sign () == expression::AUX_GEQ) && (vval >= fval)) ? 0. :
457  ((v -> sign () == expression::AUX_LEQ) && (vval <= fval)) ? 0. : fabs (vval - fval),
458 
459  ratio = (CoinMax (1., fabs (vval)) /
460  CoinMax (1., fabs (fval)));
461 
462  //printf ("checkinf --> v=%e f=%e den=%e ret=%e ratio=%e\n", vval, fval, denom, retval, ratio);
463 
464  double deldenom = delta/denom;
465  if (maxViol <= deldenom) maxViol = deldenom;
466 
467  if ((delta > 0.) &&
468  (((ratio > 2.) || // check delta > 0 to take into account semi-auxs
469  (ratio < .5)) ||
470  ((delta /= denom) > CoinMin (COUENNE_EPS, feas_tolerance_)))) {
471 
472  Jnlst () -> Printf (Ipopt::J_MOREVECTOR, J_PROBLEM, "checkAux(): auxiliary %d violates tolerance %g by %g (deldenom: %g ratio %g)\n", i, feas_tolerance_, delta, deldenom, ratio);
473 
474  isFeas = false;
475 
476  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkAux(): auxiliary %d violates tolerance %g by %g (deldenom: %g ratio %g COUENNE_EPS: %g)\n",
477  i, feas_tolerance_, delta, deldenom, ratio, COUENNE_EPS);
478 
479  if (stopAtFirstViol)
480  break;
481  }
482  }
483 
484  return isFeas;
485 }
486 
487 
488 /************************************************************************/
489 bool CouenneProblem::checkCons(const CouNumber *sol,
490  const bool stopAtFirstViol,
491  const double precision, double &maxViol) const {
492 
493  bool isFeas = true;
494  for (int i=0; i<nCons(); i++) {
495 
496  CouenneConstraint *c = Con(i);
497 
498  CouNumber
499  body = (*(c -> Body ())) (),
500  lhs = (*(c -> Lb ())) (),
501  rhs = (*(c -> Ub ())) ();
502 
503  double denomUb = 1 + CoinMax (fabs (body), fabs (rhs));
504  double denomLb = 1 + CoinMax (fabs (body), fabs (lhs));
505  double violUb = 0, violRelUb = 0, violAbsUb = 0;
506  if(rhs < COUENNE_INFINITY) {
507  violAbsUb = body - rhs;
508  violRelUb = violAbsUb / denomUb;
509  violUb = violAbsUb - precision * denomUb;
510 
511 #ifdef FM_USE_ABS_VIOL_CONS
512  if (maxViol <= violAbsUb) maxViol = violAbsUb;
513 #else
514  if (maxViol <= violRelUb) maxViol = violRelUb;
515 #endif
516 
517  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "violAbsUb: %12.10f violRelUb: %12.10f violUb: %12.10f maxViol: %12.10f\n", violAbsUb, violRelUb, violUb, maxViol);
518  }
519 
520  double violLb = 0, violRelLb = 0, violAbsLb = 0;
521  if(lhs > -COUENNE_INFINITY) {
522  violAbsLb = - body + lhs;
523  violRelLb = violAbsLb / denomLb;
524  violLb = violAbsLb - precision * denomLb;
525 
526 #ifdef FM_USE_ABS_VIOL_CONS
527  if (maxViol <= violAbsLb) maxViol = violAbsLb;
528 #else
529  if (maxViol <= violRelLb) maxViol = violRelLb;
530 #endif
531 
532  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "violAbsLb: %12.10f violRelLb: %12.10f violLb: %12.10f maxViol: %12.10f\n", violAbsLb, violRelLb, violLb, maxViol);
533  }
534 
535 #ifdef FM_USE_ABS_VIOL_CONS
536  if((violAbsUb > precision) || (violAbsLb > precision)) {
537  if (Jnlst()->ProduceOutput(Ipopt::J_MOREVECTOR, J_PROBLEM)) {
538 
539  printf ("CouenneProblem::checkCons(): constraint %d violated (lhs=%+e body=%+e rhs=%+e, absolute violation: %g)\n", i, lhs, body, rhs, CoinMax(violAbsUb, violAbsLb));
540  c -> print ();
541  }
542 
543  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkCons(): constraint %d violated (lhs=%+e body=%+e rhs=%+e, absolute violation: %g)\n", i, lhs, body, rhs, CoinMax (violAbsUb, violAbsLb));
544 
545  isFeas = false;
546  if(stopAtFirstViol) {
547  break;
548  }
549  }
550 #else /* not FM_USE_ABS_VIOL_CONS */
551  if((violUb > 0) || (violLb > 0)) {
552  if (Jnlst()->ProduceOutput(Ipopt::J_MOREVECTOR, J_PROBLEM)) {
553 
554  printf ("CouenneProblem::checkCons(): constraint %d violated (lhs=%+e body=%+e rhs=%+e, relative violation: %g)\n", i, lhs, body, rhs, CoinMax(violRelUb, violRelLb));
555  c -> print ();
556  }
557 
558  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkCons(): constraint %d violated (lhs=%+e body=%+e rhs=%+e, relative violation: %g)\n", i, lhs, body, rhs, CoinMax (violRelUb, violRelLb));
559 
560  isFeas = false;
561  if(stopAtFirstViol) {
562  break;
563  }
564  }
565 #endif /* not FM_USE_ABS_VIOL_CONS */
566  }
567  return(isFeas);
568 }
569 
570 /************************************************************************/
571 bool CouenneProblem::checkNLP2(const double *solution,
572  const double obj, const bool careAboutObj,
573  const bool stopAtFirstViol,
574  const bool checkAll,
575  const double precision) const {
576 
577  if (careAboutObj && stopAtFirstViol) {
578  printf("CouenneProblem::checkNLP2(): ### ERROR: careAboutObj: true and stopAtFirstViol: true are incompatible\n");
579  exit(1);
580  }
581 
582  const std::vector<int> listInt = getRecordBestSol()->getListInt();
583  bool isFeas = false;
584  double maxViolCouSol = 0;
585  double maxViolRecSol = 0;
586 
587  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
588 
589  const bool *initIsInt = getRecordBestSol()->getInitIsInt();
590  printf("Integrality:\n");
591  for (register int i=0; i<nVars(); i++) {
592 
593  if (variables_ [i] -> Multiplicity () <= 0)
594  continue;
595 
596  if(initIsInt[i]) {
597  printf(" %d", i);
598  }
599  }
600  printf("\n");
601 
602  printf("VAR:\n");
603  for (register int i=0; i<nVars(); i++) {
604 
605  if (variables_ [i] -> Multiplicity () <= 0)
606  continue;
607  exprVar *v = variables_ [i];
608  if( (v -> Type () == VAR)) {
609  printf(" %d", i);
610  }
611  }
612  printf("\n");
613 
614  printf("AUX:\n");
615  for (register int i=0; i<nVars(); i++) {
616 
617  if (variables_ [i] -> Multiplicity () <= 0)
618  continue;
619  exprVar *v = variables_ [i];
620  if( (v -> Type () == AUX)) {
621  printf(" %d", i);
622  }
623  }
624  printf("\n");
625 
626  printf("mult 0:\n");
627  for (register int i=0; i<nVars(); i++) {
628 
629  if (variables_ [i] -> Multiplicity () <= 0) {
630  printf(" %d", i);
631  }
632  }
633  printf("\n");
634  }
635 
636  if ((Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM))) {
637  printf ("checking solution:\n");
638  for (int i=0; i<nOrigVars_ - ndefined_; i++)
639  printf ("%.12e ", solution [i]);
640  printf ("\nCouenneProblem::checkNLP2(): Start checking recomputed_solution\n");
641  }
642 
643  // check integrality of integer constrained variables
644 
645  int from = 0;
646  bool isFeasRec = checkInt(solution, from, nOrigVars_ - ndefined_, listInt,
647  false, stopAtFirstViol,
648  precision, maxViolRecSol);
649  bool isFeasCou = isFeasRec;
650  maxViolCouSol = maxViolRecSol;
651 
652  if (stopAtFirstViol && !isFeasRec && Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM))
653  printf("CouenneProblem::checkNLP2(): recomputed_solution is infeasible (some orig vars not integer feasible; violation: %12.10g)\n", maxViolRecSol);
654 
655 #ifdef CHECK
656  if(getRecordBestSol()->getCardInitDom() != nVars()) {
657  printf("CouenneProblem::checkNLP2(): ### ERROR: cardInitDom: %d nVars: %d\n", getRecordBestSol()->getCardInitDom(), nVars());
658  exit(1);
659  }
660  if(getInitDomLb() == NULL) {
661  printf("CouenneProblem::checkNLP2(): ### WARNING: initDomLb == NULL\n");
662  }
663  if(getInitDomUb() == NULL) {
664  printf("CouenneProblem::checkNLP2(): ### WARNING: initDomUb == NULL\n");
665  }
666 #endif
667 
668  // install NL solution candidate and original bounds in evaluation structure
669  // bounds are important so that getAuxs below works properly
670  domain_.push(nVars(), solution, getRecordBestSol()->getInitDomLb(),
671  getRecordBestSol()->getInitDomUb(), false);
672 
673  CouNumber *couRecSol = new CouNumber[nVars()];
674  CoinCopyN (solution, nOrigVars_ - ndefined_, couRecSol);
675  getAuxs(couRecSol);
676  //CoinCopyN (solution, nOrigVars_, couRecSol);
677 
678  domain_.pop (); // getting rid of current domain now as won't be used again
679 
680  // install couRecSol in evaluation structure
681  domain_.push(nVars(), couRecSol,
682  getRecordBestSol()->getInitDomLb(),
683  getRecordBestSol()->getInitDomUb(), false);
684 
685  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
686  printf ("checkNLP2(): recomputed_solution: %d vars -------------------\n", domain_.current () -> Dimension ());
687  double maxDelta = 0;
688  for (int i=0; i<nVars (); i++) {
689  exprVar *v = variables_ [i];
690  if (v -> Multiplicity () <= 0)
691  continue;
692  if(i < nOrigVars_ - ndefined_) {
693  double soli = solution[i];
694  double domi = domain_.x (i);
695  double domlbi = domain_.lb (i);
696  double domubi = domain_.ub (i);
697  printf ("%4d %+e %+e [%+e %+e] %+e\n", i, solution[i], domain_.x (i), domain_.lb (i), domain_.ub (i), solution[i] - domain_.x (i));
698  if (maxDelta < fabs(solution[i] - domain_.x(i))) maxDelta = fabs (solution[i] - domain_.x(i));
699  }
700  else {
701  if(v -> Type() == AUX) {
702  printf ("%4d ------ %+e [%+e %+e]\n", i, domain_.x (i), domain_.lb (i), domain_.ub (i));
703  }
704  }
705  printf ("maxDelta: %.12g\n", maxDelta);
706  }
707  }
708 
709  if(checkAll) {
710  if(!stopAtFirstViol || isFeasRec) {
711  bool isFeasInt = checkInt(couRecSol, 0, nVars(), listInt,
712  false, stopAtFirstViol,
713  precision, maxViolRecSol);
714 
715  if (!isFeasInt) {
716  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): recomputed_solution is infeasible (some aux vars not integer feasible; violation: %12.10g)\n", maxViolRecSol);
717 
718  isFeasRec = false;
719  }
720  }
721  }
722 
723  double objRecSol = checkObj(couRecSol, precision);
724  double objCouSol = 0;
725 
726  if(checkAll) {
727  if(!stopAtFirstViol || isFeasRec) {
728  bool isFeasBnd = checkBounds(couRecSol, stopAtFirstViol,
729  precision, maxViolRecSol);
730 
731  if(!isFeasBnd) {
732 
733  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): recomputed_solution is infeasible (violated bounds; violation: %12.10g)\n", maxViolRecSol);
734 
735  isFeasRec = false;
736  }
737  }
738 
739  if(!stopAtFirstViol || isFeasRec) {
740  bool isFeasAux = checkAux(couRecSol, stopAtFirstViol,
741  precision, maxViolRecSol);
742 
743  if(!isFeasAux) {
744 
745  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): recomputed_solution is infeasible (violated Aux; violation: %12.10g)\n", maxViolRecSol);
746 
747  isFeasRec = false;
748  }
749  }
750  }
751 
752  if(!stopAtFirstViol || isFeasRec) {
753  bool isFeasCons = checkCons(couRecSol, stopAtFirstViol,
754  precision, maxViolRecSol);
755 
756  if(!isFeasCons) {
757 
758  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): recomputed_solution is infeasible (violated constraint; violation: %12.10g)\n", maxViolRecSol);
759 
760  isFeasRec = false;
761  }
762  }
763 
764  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): end check recomputed_solution (maxViol: %12.10g)\n", maxViolRecSol);
765 
766  double objErrorRecSol = objRecSol - obj;
767  if(!careAboutObj)
768  objErrorRecSol = 0;
769 
770  CouNumber *couSol = new CouNumber[nVars()];
771  bool useRecSol = false;
772  if(isFeasRec && (objErrorRecSol < precision)) {
773  useRecSol = true;
774  isFeas = true;
775  }
776  else {
777 
778  if(checkAll) { // otherwise duplicates above calculations
779 
780  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): Start checking solution (maxViol: %g)\n", maxViolCouSol);
781 
782  CoinCopyN(solution, nVars(), couSol);
783  restoreUnusedOriginals(couSol);
784  domain_.push(nVars(), couSol, getRecordBestSol()->getInitDomLb(),
785  getRecordBestSol()->getInitDomUb(), false);
786 
787  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
788  printf ("checkNLP2(): solution: %d vars -------------------\n", domain_.current () -> Dimension ());
789  double maxDelta = 0;
790  for (int i=0; i<domain_.current()->Dimension(); i++) {
791  //printf ("%4d %.12g %.12g [%.12g %.12g]\n", i, solution[i], domain_.x(i), domain_.lb(i), domain_.ub(i));
792  printf ("%4d %+e %+e [%+e %+e] %+e\n", i, solution[i], domain_.x (i), domain_.lb (i), domain_.ub (i), solution[i] - domain_.x (i));
793  if (fabs (solution[i] - domain_.x(i)) > maxDelta)
794  maxDelta = fabs(solution[i] - domain_.x(i));
795  }
796  printf("maxDelta: %.12g\n", maxDelta);
797  }
798 
799  if(!stopAtFirstViol || isFeasCou) {
800  bool isFeasInt = checkInt(couSol, nOrigVars_ - ndefined_, nVars(), listInt,
801  false, stopAtFirstViol,
802  precision, maxViolCouSol);
803  if(!isFeasInt) {
804 
805  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): solution is infeasible (some aux vars not integer feasible; violation: %12.10g)\n", maxViolCouSol);
806 
807  isFeasCou = false;
808  }
809  }
810 
811  objCouSol = checkObj(couSol, precision);
812 
813  if(!stopAtFirstViol || isFeasCou) {
814  bool isFeasCouBnd = checkBounds(couSol, stopAtFirstViol,
815  precision, maxViolCouSol);
816  if(!isFeasCouBnd) {
817 
818  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): solution is infeasible (some bounds are violated; violation: %12.10g)\n", maxViolCouSol);
819 
820  isFeasCou = false;
821  }
822  }
823 
824  if(!stopAtFirstViol || isFeasCou) {
825  bool isFeasCouAux = checkAux(couSol, stopAtFirstViol,
826  precision, maxViolCouSol);
827  if(!isFeasCouAux) {
828 
829  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): solution is infeasible (violated Aux; violation: %12.10g)\n", maxViolCouSol);
830 
831  isFeasCou = false;
832  }
833  }
834 
835  if(!stopAtFirstViol || isFeasCou) {
836  bool isFeasCouCons = checkCons(couSol, stopAtFirstViol,
837  precision, maxViolCouSol);
838  if(!isFeasCouCons) {
839 
840  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): solution is infeasible (violated constraint; violation: %12.10g)\n", maxViolCouSol);
841 
842  isFeasCou = false;
843  }
844  }
845 
846  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): end check solution (maxViol: %12.10g)\n", maxViolCouSol);
847 
848  double objErrorCouSol = objCouSol - obj;
849  if(!careAboutObj) {
850  objErrorCouSol = 0;
851  }
852  double delObj = objErrorCouSol - objErrorRecSol;
853  double delViol = maxViolCouSol - maxViolRecSol;
854 
855  if(isFeasRec) {
856  if(isFeasCou) {
857  // careAboutObj must be true
858  if(delObj > 0) {
859  useRecSol = true;
860  }
861  else {
862  useRecSol = false;
863  }
864  isFeas = true;
865  }
866  else { /* isFeasRec == true and isFeasCou == false */
867  useRecSol = true;
868  isFeas = true;
869  }
870  }
871  else { /* isFeasRec == false */
872  if(isFeasCou) {
873  useRecSol = false;
874  isFeas = true;
875  }
876  else { /* isFeasRec == false and isFeasCou == false */
877  isFeas = false;
878  if(careAboutObj) {
879  if(fabs(delViol) < 10 * precision) {
880  useRecSol = (delObj <= 0 ? false : true);
881  }
882  else {
883  if(fabs(delObj) < 10 * precision) {
884  useRecSol = (delViol > 0 ? false : true);
885  }
886  else {
887  double ratObj = fabs(delObj)/(1+fabs(obj));
888  if(ratObj < fabs(delViol)) {
889  useRecSol = (delViol > 0 ? false : true);
890  }
891  else {
892  useRecSol = (delObj > 0 ? false : true);
893  }
894  }
895  }
896  }
897  else {
898  if(delViol < 0) {
899  useRecSol = false;
900  }
901  else {
902  useRecSol = true;
903  }
904  }
905  useRecSol = true;
906  }
907  }
908  domain_.pop (); // pop couSol
909  }
910  }
911 
912  double maxViol = 0;
913 
914  if(!stopAtFirstViol || isFeas) {
915  if(useRecSol) {
916  recBSol->setModSol(couRecSol, nVars(), objRecSol, maxViolRecSol);
917  maxViol = maxViolRecSol;
918 
919  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): select recomputed_solution (maxViol: %12.10g)\n", maxViol);
920  }
921  else {
922  recBSol -> setModSol(couSol, nVars(), objCouSol, maxViolCouSol);
923  maxViol = maxViolCouSol;
924 
925  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): select solution (maxViol: %12.10g)\n", maxViol);
926  }
927  }
928 
929  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
930  if(isFeas) {
931 
932  // printf ("Solution: [");
933  // for (int i = 0; i < nVars (); ++i)
934  // printf ("%g ", couSol [i]);
935  // printf ("]\n");
936 
937  printf ("checkNLP2(): RETURN: feasible (maxViol: %g)\n", maxViol);
938  }
939  else {
940  printf ("checkNLP2(): RETURN: recomputed_solution and solution are infeasible\n");
941  if(!stopAtFirstViol) {
942  printf("(maxViol: %g)\n", maxViol);
943  }
944  else {
945  printf("\n");
946  }
947  }
948  }
949 
950  delete[] couSol;
951  delete[] couRecSol;
952 
953  domain_.pop (); // pop bounds
954 
955  return isFeas;
956 }
957 
958 // comprehensive method to call one of the two variants
959 bool CouenneProblem::checkNLP0 (const double *solution,
960  double &obj,
961  bool recompute_obj,
962  const bool careAboutObj,
963  const bool stopAtFirstViol,
964  const bool checkAll,
965  const double precision) const {
966 
967  bool retval;
968 
969 #ifdef FM_CHECKNLP2
970 
971  retval = checkNLP2 (solution,
972  obj,
973  careAboutObj,
974  stopAtFirstViol,
975  checkAll,
976  (precision < 0.) ? feas_tolerance_ : precision);
977 
978  if (retval)
979  obj = getRecordBestSol () -> getModSolVal ();
980 
981 #else
982 
983  retval = checkNLP1 (solution, obj, recompute_obj);
984 
985 #endif
986 
987  return retval;
988 }
COUENNE_round
#define COUENNE_round(x)
Definition: CouennePrecisions.hpp:42
Couenne::CouenneProblem::checkInt
bool checkInt(const CouNumber *sol, const int from, const int upto, const std::vector< int > listInt, const bool origVarOnly, const bool stopAtFirstViol, const double precision, double &maxViol) const
check integrality of vars in sol with index between from and upto (original vars only if origVarOnly ...
Definition: checkNLP.cpp:313
Couenne::CouenneProblem::nVars
int nVars() const
Total number of variables.
Definition: CouenneProblem.hpp:380
diff
static Bigint * diff(Bigint *a, Bigint *b)
Definition: OSdtoa.cpp:1120
Couenne::Domain::ub
CouNumber & ub(register int index)
current upper bound
Definition: CouenneDomain.hpp:158
Couenne::Domain::lb
CouNumber & lb(register int index)
current lower bound
Definition: CouenneDomain.hpp:157
Couenne::CouenneProblem::print
void print(std::ostream &=std::cout)
Display current representation of problem: objective, linear and nonlinear constraints, and auxiliary variables.
Definition: problemIO.cpp:24
Couenne::CouenneObjective::Body
expression * Body() const
get body
Definition: CouenneProblemElem.hpp:135
Couenne::CouenneProblem::checkNLP0
bool checkNLP0(const double *solution, double &obj, bool recompute_obj=false, const bool careAboutObj=false, const bool stopAtFirstViol=true, const bool checkAll=false, const double precision=-1) const
And finally a method to get both.
Definition: checkNLP.cpp:959
Couenne::CouenneProblem::checkAux
bool checkAux(const CouNumber *sol, const bool stopAtFirstViol, const double precision, double &maxViol) const
returns true iff value of all auxiliaries are within bounds
Definition: checkNLP.cpp:403
Couenne::CouenneProblem::checkNLP2
bool checkNLP2(const double *solution, const double obj, const bool careAboutObj, const bool stopAtFirstViol, const bool checkAll, const double precision) const
Return true if either solution or recomputed_solution obtained using getAuxs() from the original vari...
Definition: checkNLP.cpp:571
Couenne::CouenneProblem::checkObj
double checkObj(const CouNumber *sol, const double &precision) const
Recompute objective value for sol.
Definition: checkNLP.cpp:280
Couenne::DomainPoint::Dimension
int Dimension()
return dimension_
Definition: CouenneDomain.hpp:82
Couenne::expression::Image
virtual expression * Image() const
return pointer to corresponding expression (for auxiliary variables only)
Definition: CouenneExpression.hpp:100
e
void fint fint fint real fint real real real real real real real real real * e
Definition: BonBqpdSolver.cpp:27
Couenne::CouenneConstraint
Class to represent nonlinear constraints.
Definition: CouenneProblemElem.hpp:39
Couenne::CouenneProblem::restoreUnusedOriginals
void restoreUnusedOriginals(CouNumber *=NULL) const
Some originals may be unused due to their zero multiplicity (that happens when they are duplicates)...
Definition: CouenneRestoreUnused.cpp:54
Couenne::CouenneProblem::checkBounds
bool checkBounds(const CouNumber *sol, const bool stopAtFirstViol, const double precision, double &maxViol) const
Check bounds; returns true iff feasible for given precision.
Definition: checkNLP.cpp:361
Couenne::CouenneRecordBestSol::getListInt
std::vector< int > getListInt() const
Definition: CouenneRecordBestSol.hpp:63
ratio
static double ratio(Bigint *a, Bigint *b)
Definition: OSdtoa.cpp:1460
Couenne::CouenneProblem::getAuxs
void getAuxs(CouNumber *) const
Get auxiliary variables from original variables.
Definition: problem.cpp:162
Couenne::CouenneProblem::ndefined_
int ndefined_
Number of &quot;defined variables&quot; (aka &quot;common expressions&quot;)
Definition: CouenneProblem.hpp:228
Couenne::expression::Index
virtual int Index() const
Return index of variable (only valid for exprVar and exprAux)
Definition: CouenneExpression.hpp:73
Couenne::CouenneProblem::getRecordBestSol
CouenneRecordBestSol * getRecordBestSol() const
returns recorded best solution
Definition: CouenneProblem.hpp:821
Couenne::CouenneProblem::Jnlst
ConstJnlstPtr Jnlst() const
Provide Journalist.
Definition: CouenneProblem.cpp:288
Couenne::CouenneProblem::Con
CouenneConstraint * Con(int i) const
i-th constraint
Definition: CouenneProblem.hpp:415
Couenne::Domain::pop
void pop()
restore previous point
Definition: domain.cpp:255
COUENNE_EPS
#define COUENNE_EPS
Definition: CouennePrecisions.hpp:19
Couenne::CouenneProblem::variables_
std::vector< exprVar * > variables_
Variables (original, auxiliary, and defined)
Definition: CouenneProblem.hpp:195
Couenne::CouenneProblem::recBSol
CouenneRecordBestSol * recBSol
Definition: CouenneProblem.hpp:310
Couenne::CouenneProblem::domain_
Domain domain_
current point and bounds;
Definition: CouenneProblem.hpp:202
Couenne::CouenneProblem::Ub
CouNumber * Ub() const
Return vector of upper bounds.
Definition: CouenneProblem.hpp:448
infeasible
static const int infeasible
Definition: Couenne.cpp:68
Couenne::CouNumber
double CouNumber
main number type in Couenne
Definition: CouenneTypes.hpp:100
Couenne::CouenneProblem::nOrigVars_
int nOrigVars_
Number of original variables.
Definition: CouenneProblem.hpp:237
Couenne::CouenneProblem::feas_tolerance_
CouNumber feas_tolerance_
feasibility tolerance (to be used in checkNLP)
Definition: CouenneProblem.hpp:270
COUENNE_INFINITY
#define COUENNE_INFINITY
Definition: CouennePrecisions.hpp:32
Couenne::CouenneProblem::X
CouNumber * X() const
Return vector of variables.
Definition: CouenneProblem.hpp:446
Couenne::Domain::x
CouNumber & x(register int index)
current variable
Definition: CouenneDomain.hpp:156
Couenne::exprVar::sign
virtual enum auxSign sign() const
return its sign in the definition constraint
Definition: CouenneExprVar.hpp:217
Couenne::Domain::push
void push(int dim, CouNumber *x, CouNumber *lb, CouNumber *ub, bool copy=true)
save current point and start using another
Definition: domain.cpp:166
Couenne::exprVar
variable-type operator
Definition: CouenneExprVar.hpp:45
Couenne::CouenneProblem::nCons
int nCons() const
Get number of constraints.
Definition: CouenneProblem.hpp:373
Couenne::CouenneProblem::checkCons
bool checkCons(const CouNumber *sol, const bool stopAtFirstViol, const double precision, double &maxViol) const
returns true iff value of all auxiliaries are within bounds
Definition: checkNLP.cpp:489
CouenneExprVar.hpp
Couenne::CouenneRecordBestSol::setModSol
void setModSol(const double *givenModSol, const int givenModCard, const double givenModVal, const double givenModMaxViol)
Definition: CouenneRecordBestSol.cpp:333
Couenne::CouenneProblem::checkNLP
bool checkNLP(const double *solution, double &obj, bool recompute=false) const
Check if solution is MINLP feasible.
Definition: checkNLP.cpp:23
Couenne::expression
Expression base class.
Definition: CouenneExpression.hpp:48
Couenne::CouenneProblem::Obj
CouenneObjective * Obj(int i) const
i-th objective
Definition: CouenneProblem.hpp:416
CouenneRecordBestSol.hpp
CouenneProblemElem.hpp
Couenne::Domain::current
DomainPoint * current()
return current point
Definition: CouenneDomain.hpp:154
CouenneProblem.hpp
Couenne::J_PROBLEM
const Ipopt::EJournalCategory J_PROBLEM(Ipopt::J_USER4)
c
real c
Definition: BonBqpdSolver.cpp:98
Couenne::CouenneProblem::Lb
CouNumber * Lb() const
Return vector of lower bounds.
Definition: CouenneProblem.hpp:447
Generated by   doxygen 1.8.5