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

Last change on this file since 940 was 940, checked in by pbelotti, 8 years ago

applied patch by Stefan on signed powers. Dormant for now, not verified if stable yet. Some code cleanup in checkNLP (debug output). Revert to SCIP/trunk

  • Property svn:keywords set to Author Date Id Revision
File size: 28.8 KB
Line 
1/* $Id: checkNLP.cpp 940 2013-01-13 19:49:02Z pbelotti $
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        v = (*(variables_ [i])) (), 
422        f = (*(variables_ [i] -> Image ())) ();
423
424      printf ("[%g,%g]\n", v, f);
425
426      printf ("checking aux -- %+.12e = %+.12e [%+.12e] ", v, f, v-f);
427      variables_ [i] -> print ();             printf (" := ");
428      variables_ [i] -> Image () -> print (); printf ("\n");
429    }
430   
431    double 
432      vval = (*v) (),
433      fval = (*(v -> Image ())) (),
434      denom  = CoinMax (1., v -> Image () -> gradientNorm (X ()));
435   
436    // check if fval is a number (happens with e.g. w13 = w12/w5 and w5=0, see test/harker.nl)
437    if (CoinIsnan (fval)) {
438      fval = vval + 1.;
439      denom = 1.;
440    }
441   
442    if (fabs (fval) > COUENNE_INFINITY)
443      fval = COUENNE_INFINITY;
444
445    double
446      delta = 
447      ((v -> sign () == expression::AUX_GEQ) && (vval >= fval)) ? 0. : 
448      ((v -> sign () == expression::AUX_LEQ) && (vval <= fval)) ? 0. : fabs (vval - fval),
449     
450      ratio = (CoinMax (1., fabs (vval)) / 
451               CoinMax (1., fabs (fval)));
452   
453    //printf ("checkinf --> v=%e f=%e den=%e ret=%e ratio=%e\n", vval, fval, denom, retval, ratio);
454   
455    double deldenom = delta/denom;
456    if (maxViol <= deldenom) maxViol =  deldenom;
457
458    if ((delta > 0.) &&
459        (((ratio > 2.)  ||  // check delta > 0 to take into account semi-auxs
460          (ratio <  .5)) ||
461         ((delta /= denom) > CoinMin (COUENNE_EPS, feas_tolerance_)))) {
462
463      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);
464     
465      isFeas = false;
466
467      Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkAux(): auxiliary %d violates tolerance %g by %g (deldenom: %g  ratio %g  COUENNE_EPS: %g)\n", 
468                          i, feas_tolerance_, delta, deldenom, ratio, COUENNE_EPS);
469
470      if (stopAtFirstViol)
471        break;
472    }
473  }
474
475  return isFeas;
476}
477
478
479/************************************************************************/
480bool CouenneProblem::checkCons(const CouNumber *sol,
481                               const bool stopAtFirstViol,
482                               const double precision, double &maxViol) const {
483
484  bool isFeas = true;
485  for (int i=0; i<nCons(); i++) {
486   
487    CouenneConstraint *c = Con(i);
488   
489    CouNumber
490      body = (*(c -> Body ())) (),
491      lhs  = (*(c -> Lb   ())) (),
492      rhs  = (*(c -> Ub   ())) ();
493   
494    double denomUb = 1 + CoinMax (fabs (body), fabs (rhs));
495    double denomLb = 1 + CoinMax (fabs (body), fabs (lhs));
496    double violUb = 0, violRelUb = 0, violAbsUb = 0;
497    if(rhs < COUENNE_INFINITY) {
498      violAbsUb = body - rhs;
499      violRelUb = violAbsUb / denomUb;
500      violUb = violAbsUb - precision * denomUb;
501
502#ifdef FM_USE_ABS_VIOL_CONS
503      if (maxViol <= violAbsUb) maxViol = violAbsUb;
504#else
505      if (maxViol <= violRelUb) maxViol = violRelUb;
506#endif
507
508      Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "violAbsUb: %12.10f  violRelUb: %12.10f  violUb: %12.10f maxViol: %12.10f\n", violAbsUb, violRelUb, violUb, maxViol);
509    }
510
511    double violLb = 0, violRelLb = 0, violAbsLb = 0;
512    if(lhs > -COUENNE_INFINITY) {
513      violAbsLb = - body + lhs;
514      violRelLb = violAbsLb / denomLb;
515      violLb = violAbsLb - precision * denomLb;
516
517#ifdef FM_USE_ABS_VIOL_CONS
518      if (maxViol <= violAbsLb) maxViol = violAbsLb;
519#else
520      if (maxViol <= violRelLb) maxViol = violRelLb;
521#endif
522
523      Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "violAbsLb: %12.10f  violRelLb: %12.10f  violLb: %12.10f maxViol: %12.10f\n", violAbsLb, violRelLb, violLb, maxViol);
524    }
525
526#ifdef FM_USE_ABS_VIOL_CONS
527    if((violAbsUb > precision) || (violAbsLb > precision)) {
528      if (Jnlst()->ProduceOutput(Ipopt::J_MOREVECTOR, J_PROBLEM)) {
529       
530        printf ("CouenneProblem::checkCons(): constraint %d violated (lhs=%+e body=%+e rhs=%+e, absolute violation: %g)\n", i, lhs, body, rhs, CoinMax(violAbsUb, violAbsLb));
531        c -> print ();
532      }
533
534      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));
535
536      isFeas = false;
537      if(stopAtFirstViol) {
538        break;
539      }
540    }
541#else /* not FM_USE_ABS_VIOL_CONS */
542    if((violUb > 0) || (violLb > 0)) {
543      if (Jnlst()->ProduceOutput(Ipopt::J_MOREVECTOR, J_PROBLEM)) {
544       
545        printf ("CouenneProblem::checkCons(): constraint %d violated (lhs=%+e body=%+e rhs=%+e, relative violation: %g)\n", i, lhs, body, rhs, CoinMax(violRelUb, violRelLb));
546        c -> print ();
547      }
548
549      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));
550
551      isFeas = false;
552      if(stopAtFirstViol) {
553        break;
554      }
555    }
556#endif /* not FM_USE_ABS_VIOL_CONS */
557  }
558  return(isFeas);
559}
560
561/************************************************************************/
562bool CouenneProblem::checkNLP2(const double *solution, 
563                               const double obj, const bool careAboutObj,
564                               const bool stopAtFirstViol,
565                               const bool checkAll,
566                               const double precision) const {
567
568  if (careAboutObj && stopAtFirstViol) {
569    printf("CouenneProblem::checkNLP2(): ### ERROR: careAboutObj: true and stopAtFirstViol: true are incompatible\n");
570    exit(1);
571  }
572
573  const std::vector<int> listInt = getRecordBestSol()->getListInt();
574  bool isFeas = false;
575  double maxViolCouSol = 0;
576  double maxViolRecSol = 0;
577
578  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
579
580    const bool *initIsInt = getRecordBestSol()->getInitIsInt();
581    printf("Integrality:\n");
582    for (register int i=0; i<nVars(); i++) {
583
584      if (variables_ [i] -> Multiplicity () <= 0) 
585        continue;
586
587      if(initIsInt[i]) {
588        printf(" %d", i);
589      }
590    }
591    printf("\n");
592
593    printf("VAR:\n");
594    for (register int i=0; i<nVars(); i++) {
595
596      if (variables_ [i] -> Multiplicity () <= 0) 
597        continue;
598      exprVar *v = variables_ [i];
599      if(       (v -> Type () == VAR)) {
600        printf(" %d", i);
601      }
602    }
603    printf("\n");
604
605    printf("AUX:\n");
606    for (register int i=0; i<nVars(); i++) {
607
608      if (variables_ [i] -> Multiplicity () <= 0) 
609        continue;
610      exprVar *v = variables_ [i];
611      if(       (v -> Type () == AUX)) {
612        printf(" %d", i);
613      }
614    }
615    printf("\n");
616
617    printf("mult 0:\n");
618    for (register int i=0; i<nVars(); i++) {
619
620      if (variables_ [i] -> Multiplicity () <= 0) { 
621        printf(" %d", i);
622      }
623    }
624    printf("\n");
625  }
626
627  if ((Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM))) {
628    printf ("checking solution:\n");
629    for (int i=0; i<nOrigVars_ - ndefined_; i++)
630      printf ("%.12e ", solution [i]);
631    printf ("\nCouenneProblem::checkNLP2(): Start checking computed_solution\n");
632  }
633
634  // check integrality of integer constrained variables
635
636  int from = 0;
637  bool isFeasRec = checkInt(solution, from, nOrigVars_ - ndefined_, listInt,
638                            false, stopAtFirstViol,
639                            precision, maxViolRecSol);
640  bool isFeasCou = isFeasRec;
641  maxViolCouSol = maxViolRecSol;
642
643  if (stopAtFirstViol && !isFeasRec && Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) 
644    printf("CouenneProblem::checkNLP2(): modified_solution is infeasible (some orig vars not integer feasible; violation: %12.10g)\n", maxViolRecSol);
645
646#ifdef CHECK
647  if(getRecordBestSol()->getCardInitDom() != nVars()) {
648    printf("CouenneProblem::checkNLP2(): ### ERROR: cardInitDom: %d  nVars: %d\n", getRecordBestSol()->getCardInitDom(), nVars());
649    exit(1);
650  }
651  if(getInitDomLb() == NULL) {
652    printf("CouenneProblem::checkNLP2(): ### WARNING: initDomLb == NULL\n");
653  }
654  if(getInitDomUb() == NULL) {
655    printf("CouenneProblem::checkNLP2(): ### WARNING: initDomUb == NULL\n");
656  }
657#endif
658
659  // install NL solution candidate and original bounds in evaluation structure
660  // bounds are important so that getAuxs below works properly
661  domain_.push(nVars(), domain_.x(), getRecordBestSol()->getInitDomLb(), 
662               getRecordBestSol()->getInitDomUb(), false);
663
664  CouNumber *couRecSol = new CouNumber[nVars()];
665  CoinCopyN (solution, nOrigVars_ - ndefined_, couRecSol);
666  getAuxs(couRecSol);
667  //CoinCopyN (solution, nOrigVars_, couRecSol);
668
669  domain_.pop (); // getting rid of current domain now as won't be used again
670
671  // install couRecSol in evaluation structure
672  domain_.push(nVars(), couRecSol, 
673               getRecordBestSol()->getInitDomLb(), 
674               getRecordBestSol()->getInitDomUb(), false);
675
676  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
677    printf ("checkNLP2(): recomputed_solution: %d vars -------------------\n", domain_.current () -> Dimension ());
678    double maxDelta = 0;
679    for (int i=0; i<nVars (); i++) {
680      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));
681      if (maxDelta < fabs(solution[i] - domain_.x(i))) maxDelta = fabs (solution[i] - domain_.x(i));
682    }
683    printf ("maxDelta: %.12g\n", maxDelta);
684  }
685
686  if(checkAll) {
687    if(!stopAtFirstViol || isFeasRec) {
688      bool isFeasInt = checkInt(couRecSol, 0, nVars(), listInt,
689                                false, stopAtFirstViol,
690                                precision, maxViolRecSol);
691
692      if (!isFeasInt) {
693        Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): recomputed_solution is infeasible (some aux vars not integer feasible; violation: %12.10g)\n", maxViolRecSol);
694       
695        isFeasRec = false;
696      }
697    }
698  }
699
700  double objRecSol = checkObj(couRecSol, precision);
701  double objCouSol = 0;
702
703  if(checkAll) {
704    if(!stopAtFirstViol || isFeasRec) {
705      bool isFeasBnd = checkBounds(couRecSol, stopAtFirstViol, 
706                                   precision, maxViolRecSol);
707
708      if(!isFeasBnd) {
709       
710        Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): recomputed_solution is infeasible (violated bounds; violation: %12.10g)\n", maxViolRecSol);
711
712        isFeasRec = false;
713      }
714    }
715
716    if(!stopAtFirstViol || isFeasRec) {
717      bool isFeasAux = checkAux(couRecSol, stopAtFirstViol, 
718                                precision, maxViolRecSol);
719
720      if(!isFeasAux) {
721       
722        Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): recomputed_solution is infeasible (violated Aux; violation: %12.10g)\n", maxViolRecSol);
723       
724        isFeasRec = false;
725      }
726    }
727  }
728
729  if(!stopAtFirstViol || isFeasRec) {
730    bool isFeasCons = checkCons(couRecSol, stopAtFirstViol, 
731                                precision, maxViolRecSol);
732
733    if(!isFeasCons) {
734
735      Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): recomputed_solution is infeasible (violated constraint; violation: %12.10g)\n", maxViolRecSol);
736
737      isFeasRec = false;
738    }
739  }
740
741  Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): end check recomputed_solution (maxViol: %12.10g)\n", maxViolRecSol);
742
743  double objErrorRecSol = objRecSol - obj;
744  if(!careAboutObj)
745    objErrorRecSol = 0;
746
747  CouNumber *couSol = new CouNumber[nVars()];
748  bool useRecSol = false;
749  if(isFeasRec && (objErrorRecSol < precision)) {
750    useRecSol = true;
751    isFeas = true;
752  }
753  else {
754
755    if(checkAll) { // otherwise duplicates above calculations
756
757      Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): Start checking solution (maxViol: %g)\n", maxViolCouSol);
758
759      CoinCopyN(solution, nVars(), couSol);
760      restoreUnusedOriginals(couSol);
761      domain_.push(nVars(), couSol, getRecordBestSol()->getInitDomLb(), 
762                   getRecordBestSol()->getInitDomUb(), false);
763
764      if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
765        printf ("checkNLP2(): solution: %d vars -------------------\n", domain_.current () -> Dimension ());   
766        double maxDelta = 0;
767        for (int i=0; i<domain_.current()->Dimension(); i++) {
768          //printf ("%4d %.12g %.12g [%.12g %.12g]\n", i, solution[i], domain_.x(i), domain_.lb(i), domain_.ub(i));
769          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));     
770          if (fabs (solution[i] - domain_.x(i)) > maxDelta)
771            maxDelta = fabs(solution[i] - domain_.x(i));
772        }
773        printf("maxDelta: %.12g\n", maxDelta);
774      }
775
776      if(!stopAtFirstViol || isFeasCou) {
777        bool isFeasInt = checkInt(couSol, nOrigVars_ - ndefined_, nVars(), listInt,
778                                  false, stopAtFirstViol,
779                                  precision, maxViolCouSol);
780        if(!isFeasInt) {
781
782          Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): solution is infeasible (some aux vars not integer feasible; violation: %12.10g)\n", maxViolCouSol);
783         
784          isFeasCou = false;
785        }
786      }
787     
788      objCouSol = checkObj(couSol, precision);
789     
790      if(!stopAtFirstViol || isFeasCou) {
791        bool isFeasCouBnd = checkBounds(couSol, stopAtFirstViol, 
792                                        precision, maxViolCouSol);
793        if(!isFeasCouBnd) {
794         
795          Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): solution is infeasible (some bounds are violated; violation: %12.10g)\n", maxViolCouSol);
796         
797          isFeasCou = false;
798        }
799      }
800     
801      if(!stopAtFirstViol || isFeasCou) {
802        bool isFeasCouAux = checkAux(couSol, stopAtFirstViol, 
803                                     precision, maxViolCouSol);
804        if(!isFeasCouAux) {
805         
806          Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): solution is infeasible (violated Aux; violation: %12.10g)\n", maxViolCouSol);
807         
808          isFeasCou = false;
809        }
810      }
811     
812      if(!stopAtFirstViol || isFeasCou) {
813        bool isFeasCouCons = checkCons(couSol, stopAtFirstViol, 
814                                       precision, maxViolCouSol);
815        if(!isFeasCouCons) {
816         
817          Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): solution is infeasible (violated constraint; violation: %12.10g)\n", maxViolCouSol);
818         
819          isFeasCou = false;
820        }
821      }
822   
823      Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): end check solution (maxViol: %12.10g)\n", maxViolCouSol);
824   
825      double objErrorCouSol = objCouSol - obj;
826      if(!careAboutObj) {
827        objErrorCouSol = 0;
828      }
829      double delObj = objErrorCouSol - objErrorRecSol;
830      double delViol = maxViolCouSol - maxViolRecSol;
831
832      if(isFeasRec) {
833        if(isFeasCou) {
834          // careAboutObj must be true
835          if(delObj > 0) {
836            useRecSol = true;
837          }
838          else {
839            useRecSol = false;           
840          }
841          isFeas = true;
842        }
843        else { /* isFeasRec == true and isFeasCou == false */
844          useRecSol = true;
845          isFeas = true;
846        }
847      }
848      else { /* isFeasRec == false */
849        if(isFeasCou) {
850          useRecSol = false;           
851          isFeas = true;
852        }
853        else { /* isFeasRec == false and isFeasCou == false */
854          isFeas = false;
855          if(careAboutObj) {
856            if(fabs(delViol) < 10 * precision) {
857              useRecSol = (delObj <= 0 ? false : true);
858            }
859            else {
860              if(fabs(delObj) < 10 * precision) {
861                useRecSol = (delViol > 0 ? false : true);
862              }
863              else {
864                double ratObj = fabs(delObj)/(1+fabs(obj));
865                if(ratObj < fabs(delViol)) {
866                  useRecSol = (delViol > 0 ? false : true);
867                }
868                else {
869                  useRecSol = (delObj > 0 ? false : true);
870                }
871              }
872            }
873          }
874          else {
875            if(delViol < 0) {
876              useRecSol = false;
877            }
878            else {
879              useRecSol = true;
880            }
881          }
882          useRecSol = true;
883        }
884      }
885      domain_.pop (); // pop couSol
886    } 
887  }
888 
889  double maxViol = 0;
890 
891  if(!stopAtFirstViol || isFeas) {
892    if(useRecSol) {
893      recBSol->setModSol(couRecSol, nVars(), objRecSol, maxViolRecSol);
894      maxViol = maxViolRecSol;
895     
896      Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): select recomputed_solution (maxViol: %12.10g)\n", maxViol);     
897    }
898    else {
899      recBSol -> setModSol(couSol, nVars(), objCouSol, maxViolCouSol);
900      maxViol = maxViolCouSol;
901     
902      Jnlst () -> Printf (Ipopt::J_ALL, J_PROBLEM, "CouenneProblem::checkNLP2(): select solution (maxViol: %12.10g)\n", maxViol);     
903    }
904  }
905
906  if (Jnlst () -> ProduceOutput (Ipopt::J_ALL, J_PROBLEM)) {
907    if(isFeas) {
908
909      // printf ("Solution: [");
910      // for (int i = 0; i < nVars (); ++i)
911      //   printf ("%g ", couSol [i]);
912      // printf ("]\n");
913
914      printf ("checkNLP2(): RETURN: selected solution is feasible (maxViol: %g)\n", maxViol);
915    }
916    else {
917      printf ("checkNLP2(): RETURN: modified_solution and solution are infeasible\n");
918      if(!stopAtFirstViol) {
919        printf("(maxViol: %g)\n", maxViol);
920      }
921      else {
922        printf("\n");
923      }
924    }
925  }
926
927  delete[] couSol;
928  delete[] couRecSol;
929
930  domain_.pop (); // pop bounds
931   
932  return isFeas;
933}
934
935// comprehensive method to call one of the two variants
936bool CouenneProblem::checkNLP0 (const double *solution, 
937                               double &obj,
938                               bool recompute_obj, 
939                               const bool careAboutObj,
940                               const bool stopAtFirstViol,
941                               const bool checkAll,
942                               const double precision) const {
943
944  bool retval;
945
946#ifdef FM_CHECKNLP2
947
948  retval = checkNLP2 (solution,
949                      obj,
950                      careAboutObj,
951                      stopAtFirstViol,
952                      checkAll,
953                      (precision < 0.) ? feas_tolerance_ : precision);
954
955  if (retval)
956    obj = getRecordBestSol () -> getModSolVal ();
957
958#else
959
960  retval = checkNLP1 (solution, obj, recompute_obj);
961
962#endif
963
964  return retval;
965}
Note: See TracBrowser for help on using the repository browser.