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

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

replaced some iterators with const_iterator for compiling issues in MSVS, as well as drand48().

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