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

Last change on this file since 698 was 698, checked in by pbelotti, 9 years ago

FP: avoid passing depth as an argument -- retrieved from model_. Optimized some code in pool.findClosestAndReplace(). Removed global variable currentmilpmethod, replaced with static int within findSolution() --- though still not so thread safe. Diversified power operators to make it easy to insert new patch by Stefan Vigerske.

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