source: trunk/Couenne/src/main/BonCouenne.cpp @ 972

Last change on this file since 972 was 966, checked in by pbelotti, 7 years ago

show true cutoff when it appears. fix minor leak in FP.

  • Property svn:keywords set to Author Date Id Revision
File size: 15.5 KB
Line 
1// $Id: BonCouenne.cpp 966 2013-05-30 19:41:07Z tkr $
2//
3// (C) Copyright International Business Machines Corporation and Carnegie Mellon University 2006, 2007
4// All Rights Reserved.
5// This code is published under the Eclipse Public License (EPL).
6//
7// Authors :
8// Pietro Belotti, Carnegie Mellon University,
9// Pierre Bonami, International Business Machines Corporation
10//
11// Date : 12/19/2006
12
13#include <iomanip>
14#include <fstream>
15
16#include <stdlib.h>
17
18#include "CoinPragma.hpp"
19#include "CoinTime.hpp"
20#include "CoinError.hpp"
21#include "BonCouenneInterface.hpp"
22
23#include "BonCouenneSetup.hpp"
24
25#include "BonCbc.hpp"
26#include "CouenneBab.hpp"
27
28#include "CbcCutGenerator.hpp"
29#include "CouenneProblem.hpp"
30#include "CouenneCutGenerator.hpp"
31
32#include "CouenneRecordBestSol.hpp"
33
34using namespace Couenne;
35
36// the maximum difference between a printed optimum and a CouNumber
37#define PRINTED_PRECISION 1e-5
38
39#include "CouenneExprVar.hpp"
40#include "CouenneExprConst.hpp"
41#include "CouenneExprSum.hpp"
42#include "CouenneExprClone.hpp"
43#include "CouenneProblemElem.hpp"
44#include "CouenneProblem.hpp"
45#include "CouenneJournalist.hpp"
46
47#ifdef COIN_HAS_NTY
48#include "Nauty.h"
49
50// FIXME: horrible global variables. Brrr.
51#include "CouenneBranchingObject.hpp"
52#endif
53
54#ifdef COIN_HAS_SCIP
55#include "lpiswitch.h"
56#endif
57
58#include "CoinSignal.hpp"
59
60#if 0
61extern "C" {
62
63  static int nInterrupts = 0;
64  static void signal_handler (int sig) {
65
66    if (!nInterrupts) {
67      std::cerr << "[BREAK]" << std::endl;
68      abort ();
69    }
70    return;
71  }
72}
73#endif
74
75//#define FM_FRES
76
77// restored working version
78
79int main (int argc, char *argv[]) {
80
81#ifdef WIN_
82  // FIXME really want Couenne runs to be nondeterministic?
83  srand ((long) time (NULL));
84#endif
85
86  //CoinSighandler_t saveSignal = SIG_DFL;
87  //saveSignal = signal (SIGINT, signal_handler);
88
89    printf ("Couenne %s --  an Open-Source solver for Mixed Integer Nonlinear Optimization\n\
90Mailing list: couenne@list.coin-or.org\n\
91Instructions: http://www.coin-or.org/Couenne\n", 
92            strcmp (COUENNE_VERSION, "trunk") ? COUENNE_VERSION : "");
93
94  WindowsErrorPopupBlocker();
95  using namespace Ipopt;
96
97#ifdef COIN_HAS_SCIP
98  //SCIPlpiSwitchSetDefaultSolver();
99  SCIPlpiSwitchSetSolver(SCIP_LPISW_CLP);
100#endif
101
102  char * pbName = NULL;
103
104  bool infeasible = false;
105
106  try {
107
108    CouenneBab bb;
109
110    CouenneProblem *p = NULL;
111    CouenneInterface *ci = NULL;
112
113#if 0
114    //ci = new CouenneInterface;
115    p = new CouenneProblem;
116
117    p -> addVariable (false, p -> domain ());
118    p -> addVariable (false, p -> domain ());
119    p -> addVariable (false, p -> domain ());
120    p -> addVariable (false, p -> domain ());
121
122    p -> addObjective    (new exprSum (new exprClone (p->Var (1)), new exprClone (p->Var (2))), "min");
123    p -> addLEConstraint (new exprSum (new exprClone (p->Var (0)), new exprClone (p->Var (2))), new exprConst (1));
124    p -> addEQConstraint (new exprSum (new exprClone (p->Var (1)), new exprClone (p->Var (2))), new exprConst (1));
125    p -> addEQConstraint (new exprSum (new exprClone (p->Var (1)), new exprClone (p->Var (3))), new exprConst (1));
126    p -> addEQConstraint (new exprSum (new exprClone (p->Var (2)), new exprClone (p->Var (3))), new exprConst (1));
127#endif
128
129    CouenneSetup couenne;
130    CouenneCutGenerator *cg = NULL;
131    ConstJnlstPtr jnlst;
132    CouenneProblem *prob = NULL;
133
134    infeasible = !(couenne.InitializeCouenne (argv, p, NULL, ci, &bb));
135
136    // there is only one CouenneCutGenerator object; scan array until
137    // dynamic_cast returns non-NULL
138
139    if (couenne. cutGenerators () . size () > 0) {
140
141      for (std::list <Bonmin::BabSetupBase::CuttingMethod>::iterator
142             i  = couenne.cutGenerators () . begin ();
143           !cg && (i != couenne.cutGenerators () . end ()); 
144           ++i) 
145
146        cg = dynamic_cast <CouenneCutGenerator *> (i -> cgl);
147    }
148
149    // This assumes that first cut generator is CouenneCutGenerator
150    // CouenneCutGenerator *cg = dynamic_cast<CouenneCutGenerator *>
151    //   (couenne.cutGenerators().begin()->cgl);
152
153    if (cg)
154
155      cg -> setBabPtr (&bb);
156
157    else if (!infeasible) {
158
159      printf ("main(): ### ERROR: Can not get CouenneCutGenerator\n");
160      exit (-1);
161    }
162
163    // initial printout
164
165    jnlst = couenne. couennePtr () -> Jnlst ();
166    prob  = couenne. couennePtr () -> Problem ();
167
168    bb. setProblem (prob);
169
170    int retcomp = 2; // unspecified.
171
172    jnlst -> Printf (J_ERROR, J_COUENNE, "\
173Loaded instance \"%s\"\n\
174Constraints:     %8d\n\
175Variables:       %8d (%d integer)\n\
176Auxiliaries:     %8d (%d integer)\n\n",
177                     prob -> problemName ().c_str (),
178                     prob -> nOrigCons (),
179                     prob -> nOrigVars (),
180                     prob -> nOrigIntVars (),
181                     prob -> nVars    () - prob -> nOrigVars (),
182                     CoinMax (0, prob -> nIntVars () - prob -> nOrigIntVars ()));
183
184    double time_start = CoinCpuTime();
185
186#if 0
187    CouenneFeasibility feasibility;
188    bb.model().setProblemFeasibility (feasibility);
189#endif
190
191    /// update time limit (read/preprocessing might have taken some)
192    double timeLimit = 0;
193    couenne.options () -> GetNumericValue ("time_limit", timeLimit, "couenne.");
194    couenne.setDoubleParameter (Bonmin::BabSetupBase::MaxTime, 
195                                CoinMax (1., timeLimit - time_start));
196
197    //jnlst -> Printf (J_ERROR, J_COUENNE, "Starting branch-and-bound\n");
198
199    //////////////////////////////////
200
201#ifdef COIN_HAS_NTY
202    double symmGroupSize = prob -> orbitalBranching () ? prob -> getNtyInfo () -> getGroupSize () : -1;
203#endif
204
205    if (!infeasible)
206      bb (couenne); // do branch and bound
207
208#ifdef COIN_HAS_NTY
209    if (CouenneBranchingObject::nOrbBr)
210      printf ("%d orbital nontrivial branchings\n", CouenneBranchingObject::nOrbBr);
211#endif
212
213    std::cout.precision (10);
214
215    int nr=-1, nt=-1;
216    double st=-1;
217
218    if (cg) cg -> getStats (nr, nt, st);
219    else printf ("Warning, could not get pointer to CouenneCutGenerator\n");
220
221    CouenneProblem *cp = cg ? cg -> Problem () : NULL;
222
223#if defined (FM_TRACE_OPTSOL) || defined (FM_FRES)
224    double cbcLb = (infeasible ? -COIN_DBL_MAX : bb.model (). getBestPossibleObjValue ());
225    double printObj = 0;
226    bool foundSol = false;
227#endif
228
229#ifdef FM_TRACE_OPTSOL
230
231    //FILE *fSol = fopen ("bidon.sol", "w");
232
233    // if(fSol == NULL) {
234    //   printf("### ERROR: can not open bidon.sol\n");
235    //   exit(1);
236    // }
237
238    if(cp != NULL) {
239      double cbcObjVal = infeasible ? COIN_DBL_MAX : bb.model().getObjValue();
240      int modelNvars = prob -> nVars ();//bb.model().getNumCols();
241
242      CouenneRecordBestSol *rs = cp->getRecordBestSol(); 
243      const double *cbcSol = infeasible ? NULL : bb.model().getColSolution();
244      double *modCbcSol = new double[modelNvars];
245      double modCbcSolVal= 1e100, modCbcSolMaxViol = 0;
246      bool cbcSolIsFeas = false;
247
248      if(modelNvars != cp->nVars()) {
249        printf("### ERROR: modelNvars: %d nVars: %d\n", 
250               modelNvars, cp->nVars());
251        exit(1);
252      }
253
254      // round cbcSol's integer coordinates
255
256      // for (int i=0; i<modelNvars; i++)
257      //        if (cp -> Var (i) -> isDefinedInteger ())
258      //          cbcSol [i] = COUENNE_round (cbcSol [i]);
259
260      if (cbcObjVal < 1e49 && !infeasible) {
261
262#ifdef FM_CHECKNLP2
263        int cMS = rs->getCardModSol();
264        cbcSolIsFeas = cp->checkNLP2(cbcSol, 0, false, // do not care about obj
265                                     false, // do not stop at first viol
266                                     true, // checkAll
267                                     cp->getFeasTol());
268        CoinCopyN(rs->getModSol(cMS), cMS, modCbcSol);
269        modCbcSolVal = rs->getModSolVal();
270        modCbcSolMaxViol = rs->getModSolMaxViol();
271#else /* not FM_CHECKNLP2 */
272        int cMS = cp->nVars();
273        cbcSolIsFeas = cp->checkNLP(cbcSol, modCbcSolVal, true);
274        CoinCopyN(cbcSol, cMS, modCbcSol);
275        modCbcSolMaxViol = cp->getFeasTol();
276#endif /* not FM_CHECKNLP2 */
277        foundSol = true;
278      }
279
280      const double *couenneSol = rs->getSol();
281      double *modCouenneSol = new double[modelNvars];
282      double modCouenneSolVal= 1e100, modCouenneSolMaxViol = 0;
283      bool couenneSolIsFeas = false;
284
285      // round couenneSol's integer coordinates
286
287      // for (int i=0; i<modelNvars; i++)
288      //        if (cp -> Var (i) -> isDefinedInteger ())
289      //          couenneSol [i] = COUENNE_round (couenneSol [i]);
290
291      if(couenneSol != NULL) {
292#ifdef FM_CHECKNLP2
293        int cMS = rs->getCardModSol();
294        couenneSolIsFeas = cp->checkNLP2(couenneSol, 0, false, 
295                                         false, true, 
296                                         cp->getFeasTol());
297        CoinCopyN(rs->getModSol(cMS), cMS, modCouenneSol);
298        modCouenneSolVal = rs->getModSolVal();
299        modCouenneSolMaxViol = rs->getModSolMaxViol();
300#else /* not FM_CHECKNLP2 */
301        int cMS = cp->nVars();
302        couenneSolIsFeas = cp->checkNLP(couenneSol, modCouenneSolVal, true);
303        CoinCopyN(couenneSol, cMS, modCouenneSol);
304        modCouenneSolMaxViol = cp->getFeasTol();
305#endif /* not FM_CHECKNLP2 */
306        foundSol = true;
307      }
308
309      retcomp = rs -> compareAndSave (modCbcSol,     modCbcSolVal,     modCbcSolMaxViol,     cbcSolIsFeas,
310                                      modCouenneSol, modCouenneSolVal, modCouenneSolMaxViol, couenneSolIsFeas, 
311                                      modelNvars, cp->getFeasTol());
312
313      // switch (retcomp) {
314      // case -1: printf("No solution found\n"); break;
315      // case 0: printf("Best solution found by Cbc. Value: %10.4f. Tolerance: %10g\n", modCbcSolVal, modCbcSolMaxViol); break;
316      // case 1: //printf("Best solution found by Couenne  Value: %10.4f  Tolerance: %10g\n", modCouenneSolVal, modCouenneSolMaxViol); break;
317      // default: break; // never happens
318      // }
319
320      if(rs->getHasSol()) {
321        if(cbcLb > rs->getVal()) { // Best sol found by Couenne and not
322                                   // transmitted to Cbc
323          cbcLb = rs->getVal();
324        }
325        printObj = rs->getVal();
326        //rs->printSol(fSol);
327      }
328      delete[] modCbcSol;
329      delete[] modCouenneSol;
330    }
331    //fclose(fSol);
332#endif /* FM_TRACE_OPTSOL */
333
334#ifdef FM_FRES
335    if(cp != NULL) {
336      FILE *f_res = NULL;
337      f_res = fopen("fres.xxx", "r");
338      if(f_res == NULL) {
339        f_res = fopen("fres.xxx", "w");
340        fprintf(f_res, "END_OF_HEADER\n");
341      }
342      else {
343        fclose(f_res);
344        f_res = fopen("fres.xxx", "a");
345      }
346      char *pbName, shortName[256]; 
347     
348      pbName = strdup(cp -> problemName ().c_str ());
349      char *f_name_pos = strrchr(pbName, '/');
350      if(f_name_pos != NULL) {
351        strcpy(shortName, &(f_name_pos[1]));
352      }
353      else {
354        strcpy(shortName, pbName);
355      }
356     
357      fprintf(f_res, "%20s %10.4f", shortName, cbcLb);
358      if(foundSol) {
359        fprintf(f_res, " %10.4f", printObj);
360      }
361      else {
362        fprintf(f_res, "         *");
363      }
364      fprintf(f_res, " %10d %10.4f\n", infeasible ? 0 : bb.numNodes (),
365              CoinCpuTime () - time_start);
366      fclose(f_res);
367    }
368#endif
369
370    // retrieve test value to check
371    double global_opt;
372    couenne.options () -> GetNumericValue ("couenne_check", global_opt, "couenne.");
373
374    double 
375      ub = infeasible ?  COIN_DBL_MAX : bb. model (). getObjValue (),
376      lb = infeasible ? -COIN_DBL_MAX : bb. model (). getBestPossibleObjValue ();
377
378    if (cp -> getRecordBestSol () &&
379        cp       -> getRecordBestSol () -> getHasSol () &&
380        (ub > cp -> getRecordBestSol () -> getVal    ()))   
381      ub = cp -> getRecordBestSol () -> getVal ();
382
383    if (false || //(fabs (lb) > COUENNE_INFINITY / 1e4) ||
384        (lb > ub))
385      lb = ub;
386
387    char 
388      *gapstr = new char [40],
389      *lbstr  = new char [40],
390      *ubstr  = new char [40];
391
392    // CouenneSolverInterface <OsiClpSolverInterface> *csi = dynamic_cast <CouenneSolverInterface *> (couenne.continuousSolver ());
393
394    // double rootLB = csi -> rootLB ();
395
396    sprintf (lbstr,  "%10g",     lb);
397    sprintf (ubstr,  "%10g",     ub);
398    sprintf (gapstr, "%.2f%%", fabs (100. * (ub - lb) / (1. + fabs (lb))));
399
400    if (!infeasible)
401      jnlst -> Printf (J_ERROR, J_COUENNE, "\n\
402Linearization cuts added at root node:   %8d\n\
403Linearization cuts added in total:       %8d  (separation time: %gs)\n",
404                       nr, nt, st);
405
406    else jnlst -> Printf (J_ERROR, J_COUENNE, "Problem infeasible\n");
407
408    jnlst -> Printf (J_NONE, J_COUENNE, "\
409Total time:                              %8gs (%gs in branch-and-bound)\n\
410Lower bound:                           %s\n\
411Upper bound:                           %s  (gap: %s)\n\
412Branch-and-bound nodes:                  %8d\n",
413                     CoinCpuTime () - time_start,
414                     cg ? (CoinCpuTime () - CoinMax (time_start, cg -> rootTime ())) : CoinCpuTime () - time_start,
415                     ((lb <= -8.9999e12) ||
416                                       infeasible ||          (fabs (lb)             > COUENNE_INFINITY/1e4)) ? "      -inf" : lbstr,
417                     ((retcomp < 0) || infeasible ||                           (ub   > COUENNE_INFINITY/1e4)) ? "       inf" : ubstr,
418                     (                 infeasible || (CoinMax (fabs (lb), fabs (ub)) > COUENNE_INFINITY/1e4)) ? "--"         : gapstr,
419                     infeasible ? 0 : bb.numNodes ());
420
421    // if (fabs (ub - bb. model (). getObjValue ()) > COUENNE_EPS * fabs (ub))
422    //   jnlst -> Printf (J_ERROR, J_COUENNE,
423    //                 "Warning: upper bounds differ between Couenne and Cbc. Saving Couenne's (more reliable).\n");
424
425    delete [] lbstr;
426    delete [] ubstr;
427    delete [] gapstr;
428
429    if (global_opt < COUENNE_INFINITY) { // some value found in couenne.opt
430
431      double opt = infeasible ? -COIN_DBL_MAX : bb.model (). getBestPossibleObjValue ();
432
433      printf ("Global Optimum Test on %-40s %s\n", 
434              cp ? cp -> problemName ().c_str () : "unknown", 
435              (fabs (opt - global_opt) / 
436               (1. + CoinMax (fabs (opt), fabs (global_opt))) < PRINTED_PRECISION) ? 
437              (const char *) "OK" : (const char *) "FAILED");
438              //opt, global_opt,
439              //fabs (opt - global_opt));
440
441    } else // good old statistics
442
443    if (couenne.displayStats ()) { // print statistics
444
445      if (cg && !cp) printf ("Warning, could not get pointer to problem\n");
446      else
447        printf ("Stats: %-15s %4d [var] %4d [int] %4d [con] %4d [aux] "
448                "%6d [root] %8d [tot] %6g [sep] %8g [time] %8g [bb] "
449#ifdef COIN_HAS_NTY
450                "%20e [lower] %20e [upper] %7d [nodes] %.0g [sg] %d [sgc]\n",
451#else
452                "%20e [lower] %20e [upper] %7d [nodes]\n",
453#endif
454                cp ? cp -> problemName (). c_str () : "unknown",
455                (cp) ? cp -> nOrigVars     () : -1, 
456                (cp) ? cp -> nOrigIntVars  () : -1, 
457                (cp) ? cp -> nOrigCons     () : -1,
458                (cp) ? (cp -> nVars     () - 
459                        cp -> nOrigVars ()): -1,
460                nr, nt, st, 
461                CoinCpuTime () - time_start,
462                cg ? (CoinCpuTime () - cg -> rootTime ()) : CoinCpuTime (),
463                lb, //bb.model (). getBestPossibleObjValue (),
464                ub, //bb.model (). getObjValue (),
465                //bb.bestBound (),
466                //bb.bestObj (),
467                infeasible ? 0 : bb.numNodes ()
468#ifdef COIN_HAS_NTY
469                ,symmGroupSize
470                ,CouenneBranchingObject::nSGcomputations
471#endif
472);
473                //bb.iterationCount ());
474                //status.c_str (), message.c_str ());
475    }
476
477//    nlp_and_solver -> writeAmplSolFile (message, bb.bestSolution (), NULL);
478  }
479  catch(Bonmin::TNLPSolver::UnsolvedError *E) {
480     E->writeDiffFiles();
481     E->printError(std::cerr);
482    //There has been a failure to solve a problem with Ipopt.
483    //And we will output file with information on what has been changed in the problem to make it fail.
484    //Now depending on what algorithm has been called (B-BB or other) the failed problem may be at different place.
485    //    const OsiSolverInterface &si1 = (algo > 0) ? nlpSolver : *model.solver();
486  }
487  catch (Bonmin::OsiTMINLPInterface::SimpleError &E) {
488    std::cerr<<E.className()<<"::"<<E.methodName()
489             <<std::endl
490             <<E.message()<<std::endl;
491  }
492  catch (CoinError &E) {
493    std::cerr<<E.className()<<"::"<<E.methodName()
494             <<std::endl
495             <<E.message()<<std::endl;
496  }
497  catch (Ipopt::OPTION_INVALID &E)
498  {
499   std::cerr<<"Ipopt exception : "<<E.Message()<<std::endl;
500  }
501  // catch (int generic_error) {
502  //   // if (generic_error == infeasible)
503  //   //   printf ("problem infeasible\n");
504  // }
505
506//  catch(...) {
507//    std::cerr<<pbName<<" unrecognized excpetion"<<std::endl;
508//    std::cerr<<pbName<<"\t Finished \t exception"<<std::endl;
509//    throw;
510//  }
511
512  delete [] pbName;
513  return 0;
514}
Note: See TracBrowser for help on using the repository browser.