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

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

fixing drand48 calls for WIN_

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