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

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

only CouenneObjects? allowed, no integer objects (eliminated upon first call to setuplist). Cleanup

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