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

Last change on this file since 967 was 967, checked in by fmargot, 7 years ago

fix bug in checkNLP2 and handling of non-var objects in strong branching

  • Property svn:keywords set to Author Date Id Revision
File size: 29.5 KB
Line 
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
20using namespace Couenne;
21
22// check if solution is MINLP feasible
23bool 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
280double 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
313bool 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
361bool 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
403bool 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/************************************************************************/
489bool 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/************************************************************************/
571bool 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
959bool 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}
Note: See TracBrowser for help on using the repository browser.