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

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

merging 874 to fix infinite branching problem in ticket #18

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