source: trunk/Couenne/src/problem/checkNLP.cpp @ 576

Last change on this file since 576 was 576, checked in by fmargot, 9 years ago

Fix bug in CouenneProblem::checkAux(), cosmetic changes otherwise

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