source: trunk/Test/CbcMain.cpp @ 116

Last change on this file since 116 was 116, checked in by forrest, 16 years ago

decided odd hole too slow
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 39.6 KB
Line 
1// copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8#include <cassert>
9#include <typeinfo>
10#include <cstdio>
11#include <cmath>
12#include <cfloat>
13#include <string>
14#include <iostream>
15
16#define CBCVERSION "0.91"
17
18#include "CoinMpsIO.hpp"
19#include "CoinPackedMatrix.hpp"
20#include "CoinPackedVector.hpp"
21#include "CoinWarmStartBasis.hpp"
22#include "CoinTime.hpp"
23
24#include "OsiSolverInterface.hpp"
25#include "OsiCuts.hpp"
26#include "OsiRowCut.hpp"
27#include "OsiColCut.hpp"
28
29#include "CglCutGenerator.hpp"
30#include "CglGomory.hpp"
31#include "CglProbing.hpp"
32#include "CglKnapsackCover.hpp"
33#include "CglOddHole.hpp"
34#include "CglClique.hpp"
35#include "CglFlowCover.hpp"
36#include "CglMixedIntegerRounding.hpp"
37#include "CglTwomir.hpp"
38
39#include "CbcModel.hpp"
40#include "CbcCutGenerator.hpp"
41#include "CbcHeuristic.hpp"
42#include "CbcCompareActual.hpp"
43#include  "CbcParam.hpp"
44
45#ifdef COIN_USE_CLP
46#include "OsiClpSolverInterface.hpp"
47#endif
48#ifdef COIN_USE_DYLP
49#include "OsiDylpSolverInterface.hpp"
50#endif
51
52
53
54
55/* Before first solution do depth first,
56   then it is computed to hit first solution less 2%
57*/
58class CbcCompareUser  : public CbcCompareBase {
59public:
60  // Weight for each infeasibility
61  double weight_;
62  // Number of solutions
63  int numberSolutions_;
64  // Default Constructor
65  CbcCompareUser () : weight_(-1.0), numberSolutions_(0) {test_=this;};
66
67  // Copy constructor
68  CbcCompareUser ( const CbcCompareUser &rhs)
69    : CbcCompareBase(rhs)
70  {
71    weight_=rhs.weight_;
72    numberSolutions_=rhs.numberSolutions_;
73  };
74   
75  // Assignment operator
76  CbcCompareUser & operator=( const CbcCompareUser& rhs)
77  {
78    if (this!=&rhs) { 
79      CbcCompareBase::operator=(rhs);
80      weight_=rhs.weight_;
81      numberSolutions_=rhs.numberSolutions_;
82    }
83    return *this;
84  };
85
86  /// Clone
87  virtual CbcCompareBase * clone() const
88  { 
89    return new CbcCompareUser (*this);
90  };
91
92  ~CbcCompareUser() {};
93
94  /*
95     Return true if y better than x
96     Node y is better than node x if y has fewer unsatisfied (greater depth on tie) or
97     after solution weighted value of y is less than weighted value of x
98  */
99  virtual bool test (CbcNode * x, CbcNode * y) {
100    if (weight_<0.0) {
101      // before solution
102      /* printf("x %d %d %g, y %d %d %g\n",
103             x->numberUnsatisfied(),x->depth(),x->objectiveValue(),
104             y->numberUnsatisfied(),y->depth(),y->objectiveValue()); */
105      if (x->numberUnsatisfied() > y->numberUnsatisfied())
106        return true;
107      else if (x->numberUnsatisfied() < y->numberUnsatisfied())
108        return false;
109      else
110        return x->depth() < y->depth();
111    } else {
112      // after solution
113      return x->objectiveValue()+ weight_*x->numberUnsatisfied() > 
114        y->objectiveValue() + weight_*y->numberUnsatisfied();
115    }
116  }
117  // This allows method to change behavior as it is called
118  // after each solution
119  virtual void newSolution(CbcModel * model,
120                           double objectiveAtContinuous,
121                           int numberInfeasibilitiesAtContinuous) 
122  {
123    if (model->getSolutionCount()==model->getNumberHeuristicSolutions())
124      return; // solution was got by rounding
125    // set to get close to this solution
126    double costPerInteger = 
127      (model->getObjValue()-objectiveAtContinuous)/
128      ((double) numberInfeasibilitiesAtContinuous);
129    weight_ = 0.98*costPerInteger;
130    numberSolutions_++;
131    if (numberSolutions_>5)
132      weight_ =0.0; // this searches on objective
133  }
134  // This allows method to change behavior
135  virtual bool every1000Nodes(CbcModel * model, int numberNodes)
136  {
137    if (numberNodes>10000)
138      weight_ =0.0; // this searches on objective
139    return false;
140  }
141};
142
143
144#define MAXPARAMETERS 100
145
146namespace {
147
148void establishParams (int &numberParameters, CbcParam *const parameters)
149/*
150  Subroutine to establish the cbc parameter array. See the description of
151  class CbcParam for details. Pulled from main() for clarity.
152*/
153{ numberParameters = 0 ;
154
155  parameters[numberParameters++]=
156      CbcParam("?","For help",GENERALQUERY);
157  parameters[numberParameters++]=
158      CbcParam("dualT!olerance",
159               "For an optimal solution no dual infeasibility may "
160               "exceed this value",
161               1.0e-20,1.0e12,DUALTOLERANCE);
162  parameters[numberParameters++]=
163      CbcParam("primalT!olerance",
164               "For an optimal solution no primal infeasibility may "
165               " exceed this value",
166              1.0e-20,1.0e12,PRIMALTOLERANCE);
167    parameters[numberParameters++]=
168      CbcParam("inf!easibilityWeight","Each integer infeasibility is expected \
169to cost this much",
170              0.0,1.0e20,INFEASIBILITYWEIGHT);
171    parameters[numberParameters++]=
172      CbcParam("integerT!olerance","For an optimal solution \
173no integer variable may be this away from an integer value",
174              1.0e-20,0.5,INTEGERTOLERANCE);
175    parameters[numberParameters++]=
176      CbcParam("inc!rement","A valid solution must be at least this \
177much better than last integer solution",
178              -1.0e20,1.0e20,INCREMENT);
179    parameters[numberParameters++]=
180      CbcParam("allow!ableGap","Stop when gap between best possible and \
181best less than this",
182              0.0,1.0e20,ALLOWABLEGAP);
183    parameters[numberParameters++]=
184      CbcParam("ratio!Gap","Stop when gap between best possible and \
185best less than this fraction of larger of two",
186              0.0,1.0e20,GAPRATIO);
187    parameters[numberParameters++]=
188      CbcParam("fix!OnDj","Try heuristic based on fixing variables with \
189reduced costs greater than this",
190              -1.0e20,1.0e20,DJFIX);
191    parameters[numberParameters++]=
192      CbcParam("tighten!Factor","Tighten bounds using this times largest \
193activity at continuous solution",
194              1.0,1.0e20,TIGHTENFACTOR);
195    parameters[numberParameters++]=
196      CbcParam("log!Level","Level of detail in BAB output",
197              0,63,LOGLEVEL);
198    parameters[numberParameters++]=
199      CbcParam("slog!Level","Level of detail in Solver output",
200              0,63,SOLVERLOGLEVEL);
201    parameters[numberParameters++]=
202      CbcParam("maxN!odes","Maximum number of nodes to do",
203              1,999999,MAXNODES);
204    parameters[numberParameters++]=
205      CbcParam("strong!Branching","Number of variables to look at in strong branching",
206              0,999999,STRONGBRANCHING);
207    parameters[numberParameters++]=
208      CbcParam("direction","Minimize or Maximize",
209              "min!imize",DIRECTION);
210    parameters[numberParameters-1].append("max!imize");
211    parameters[numberParameters++]=
212      CbcParam("error!sAllowed","Whether to allow import errors",
213              "off",ERRORSALLOWED);
214    parameters[numberParameters-1].append("on");
215    parameters[numberParameters++]=
216      CbcParam("gomory!Cuts","Whether to use Gomory cuts",
217              "off",GOMORYCUTS);
218    parameters[numberParameters-1].append("on");
219    parameters[numberParameters-1].append("root");
220    parameters[numberParameters++]=
221      CbcParam("probing!Cuts","Whether to use Probing cuts",
222              "off",PROBINGCUTS);
223    parameters[numberParameters-1].append("on");
224    parameters[numberParameters-1].append("root");
225    parameters[numberParameters++]=
226      CbcParam("knapsack!Cuts","Whether to use Knapsack cuts",
227              "off",KNAPSACKCUTS);
228    parameters[numberParameters-1].append("on");
229    parameters[numberParameters-1].append("root");
230    parameters[numberParameters++]=
231      CbcParam("oddhole!Cuts","Whether to use Oddhole cuts",
232              "off",ODDHOLECUTS);
233    parameters[numberParameters-1].append("on");
234    parameters[numberParameters-1].append("root");
235    parameters[numberParameters++]=
236      CbcParam("clique!Cuts","Whether to use Clique cuts",
237              "off",CLIQUECUTS);
238    parameters[numberParameters-1].append("on");
239    parameters[numberParameters-1].append("root");
240    parameters[numberParameters++]=
241      CbcParam("mixed!IntegerRoundingCuts","Whether to use Mixed Integer Rounding cuts",
242              "off",MIXEDCUTS);
243    parameters[numberParameters-1].append("on");
244    parameters[numberParameters-1].append("root");
245    parameters[numberParameters++]=
246      CbcParam("flow!CoverCuts","Whether to use Flow Cover cuts",
247              "off",FLOWCUTS);
248    parameters[numberParameters-1].append("on");
249    parameters[numberParameters-1].append("root");
250    parameters[numberParameters++]=
251      CbcParam("two!MirCuts","Whether to use Two phase Mixed Integer Rounding cuts",
252              "off",TWOMIRCUTS);
253    parameters[numberParameters-1].append("on");
254    parameters[numberParameters-1].append("root");
255    parameters[numberParameters++]=
256      CbcParam("round!ingHeuristic","Whether to use Rounding heuristic",
257              "off",ROUNDING);
258    parameters[numberParameters-1].append("on");
259    parameters[numberParameters++]=
260      CbcParam("cost!Strategy","How to use costs",
261              "off",COSTSTRATEGY);
262    parameters[numberParameters-1].append("pri!orities");
263    parameters[numberParameters-1].append("pseudo!costs(not implemented yet)");
264    parameters[numberParameters++]=
265      CbcParam("keepN!ames","Whether to keep names from import",
266              "on",KEEPNAMES);
267    parameters[numberParameters-1].append("off");
268    parameters[numberParameters++]=
269      CbcParam("scaling","Whether to do scaling",
270              "on",SCALING);
271    parameters[numberParameters-1].append("off");
272    parameters[numberParameters++]=
273      CbcParam("directory","Set Default import directory",
274              DIRECTORY);
275    parameters[numberParameters++]=
276      CbcParam("solver!","Set the solver used by cbc",
277               SOLVER) ;
278    parameters[numberParameters++]=
279      CbcParam("import","Import model from mps file",
280              IMPORT);
281    parameters[numberParameters++]=
282      CbcParam("export","Export model as mps file",
283              EXPORT);
284    parameters[numberParameters++]=
285      CbcParam("save!Model","Save model to binary file",
286              SAVE);
287    parameters[numberParameters++]=
288      CbcParam("restore!Model","Restore model from binary file",
289              RESTORE);
290    parameters[numberParameters++]=
291      CbcParam("presolve","Whether to use integer presolve - be careful",
292              "off",PRESOLVE);
293    parameters[numberParameters-1].append("on");
294    parameters[numberParameters++]=
295      CbcParam("initialS!olve","Solve to continuous",
296              SOLVECONTINUOUS);
297    parameters[numberParameters++]=
298      CbcParam("branch!AndBound","Do Branch and Bound",
299              BAB);
300    parameters[numberParameters++]=
301      CbcParam("sol!ution","Prints solution to file",
302              SOLUTION);
303    parameters[numberParameters++]=
304      CbcParam("max!imize","Set optimization direction to maximize",
305              MAXIMIZE);
306    parameters[numberParameters++]=
307      CbcParam("min!imize","Set optimization direction to minimize",
308              MINIMIZE);
309    parameters[numberParameters++] =
310      CbcParam("time!Limit","Set a time limit for solving this problem",
311              1.0,(double)(60*60*24*365*10),TIMELIMIT) ;
312    parameters[numberParameters++]=
313      CbcParam("exit","Stops cbc execution",
314              EXIT);
315    parameters[numberParameters++]=
316      CbcParam("stop","Stops cbc execution",
317              EXIT);
318    parameters[numberParameters++]=
319      CbcParam("quit","Stops cbc execution",
320              EXIT);
321    parameters[numberParameters++]=
322      CbcParam("-","From stdin",
323              STDIN);
324    parameters[numberParameters++]=
325      CbcParam("stdin","From stdin",
326              STDIN);
327    parameters[numberParameters++]=
328      CbcParam("unitTest","Do unit test",
329              UNITTEST);
330    parameters[numberParameters++]=
331      CbcParam("miplib","Do some of miplib test set",
332              MIPLIB);
333    parameters[numberParameters++]=
334      CbcParam("ver!sion","Print out version",
335              VERSION);
336
337    assert(numberParameters<MAXPARAMETERS);
338
339  return ; }
340
341#ifdef COIN_USE_READLINE     
342#include <readline/readline.h>
343#include <readline/history.h>
344#endif
345
346// Returns next valid field
347
348int read_mode=1;
349char line[1000];
350char * where=NULL;
351
352std::string
353nextField()
354{
355  std::string field;
356  if (!where) {
357    // need new line
358#ifdef COIN_USE_READLINE     
359    // Get a line from the user.
360    where = readline ("Cbc:");
361     
362    // If the line has any text in it, save it on the history.
363    if (where) {
364      if ( *where)
365        add_history (where);
366      strcpy(line,where);
367    }
368#else
369    fprintf(stdout,"Cbc:");
370    fflush(stdout);
371    where = fgets(line,1000,stdin);
372#endif
373    if (!where)
374      return field; // EOF
375    where = line;
376    // clean image
377    char * lastNonBlank = line-1;
378    while ( *where != '\0' ) {
379      if ( *where != '\t' && *where < ' ' ) {
380        break;
381      } else if ( *where != '\t' && *where != ' ') {
382        lastNonBlank = where;
383      }
384      where++;
385    }
386    where=line;
387    *(lastNonBlank+1)='\0';
388  }
389  // munch white space
390  while(*where==' '||*where=='\t')
391    where++;
392  char * saveWhere = where;
393  while (*where!=' '&&*where!='\t'&&*where!='\0')
394    where++;
395  if (where!=saveWhere) {
396    char save = *where;
397    *where='\0';
398    //convert to string
399    field=saveWhere;
400    *where=save;
401  } else {
402    where=NULL;
403    field="EOL";
404  }
405  return field;
406}
407
408std::string
409getCommand(int argc, const char *argv[])
410{
411  std::string field="EOL";
412  while (field=="EOL") {
413    if (read_mode>0) {
414      if (read_mode<argc) {
415        field = argv[read_mode++];
416        if (field=="-") {
417          std::cout<<"Switching to line mode"<<std::endl;
418          read_mode=-1;
419          field=nextField();
420        } else if (field[0]!='-') {
421          if (read_mode!=2) {
422            std::cout<<"skipping non-command "<<field<<std::endl;
423            field="EOL"; // skip
424          } else {
425            // special dispensation - taken as -import name
426            read_mode--;
427            field="import";
428          }
429        } else {
430          if (field!="--") {
431            // take off -
432            field = field.substr(1);
433          } else {
434            // special dispensation - taken as -import --
435            read_mode--;
436            field="import";
437          }
438        }
439      } else {
440        field="";
441      }
442    } else {
443      field=nextField();
444    }
445  }
446  //std::cout<<field<<std::endl;
447  return field;
448}
449std::string
450getString(int argc, const char *argv[])
451{
452  std::string field="EOL";
453  if (read_mode>0) {
454    if (read_mode<argc) {
455      if (argv[read_mode][0]!='-') { 
456        field = argv[read_mode++];
457      } else if (!strcmp(argv[read_mode],"--")) {
458        field = argv[read_mode++];
459        // -- means import from stdin
460        field = "-";
461      }
462    }
463  } else {
464    field=nextField();
465  }
466  //std::cout<<field<<std::endl;
467  return field;
468}
469
470// valid = 0 okay, 1 bad, 2 not there
471int
472getIntField(int argc, const char *argv[],int * valid)
473{
474  std::string field="EOL";
475  if (read_mode>0) {
476    if (read_mode<argc) {
477      // may be negative value so do not check for -
478      field = argv[read_mode++];
479    }
480  } else {
481    field=nextField();
482  }
483  int value=0;
484  //std::cout<<field<<std::endl;
485  if (field!="EOL") {
486    // how do I check valid
487    value =  atoi(field.c_str());
488    *valid=0;
489  } else {
490    *valid=2;
491  }
492  return value;
493}
494
495
496// valid = 0 okay, 1 bad, 2 not there
497double
498getDoubleField(int argc, const char *argv[],int * valid)
499{
500  std::string field="EOL";
501  if (read_mode>0) {
502    if (read_mode<argc) {
503      // may be negative value so do not check for -
504      field = argv[read_mode++];
505    }
506  } else {
507    field=nextField();
508  }
509  double value=0.0;
510  //std::cout<<field<<std::endl;
511  if (field!="EOL") {
512    // how do I check valid
513    value = atof(field.c_str());
514    *valid=0;
515  } else {
516    *valid=2;
517  }
518  return value;
519}
520
521/// For run timing
522
523double totalTime=0.0;
524
525}       /* end unnamed namespace */
526
527
528int main (int argc, const char *argv[])
529{
530  // next {} is just to make sure all memory should be freed - for debug
531  {
532    std::ios::sync_with_stdio() ;
533/*
534  Create a vector of solver prototypes and establish a default solver. After
535  this the code is solver independent.
536
537  Creating multiple solvers is moderately expensive. If you don't want to
538  make use of this feature, best to define just one. The businesss with
539  CBC_DEFAULT_SOLVER will select the first available solver as the default,
540  unless overridden at compile time.
541
542  NOTE that processing of string parameters is case-independent, but maps are
543       case-sensitive. The solver name given here must contain only lower case
544       letters.
545*/
546    typedef std::map<std::string,OsiSolverInterface*> solverMap_t ;
547    typedef solverMap_t::const_iterator solverMapIter_t ;
548
549    solverMap_t solvers ;
550
551#   ifdef COIN_USE_CLP
552#     ifndef CBC_DEFAULT_SOLVER
553#       define CBC_DEFAULT_SOLVER "clp"
554#     endif
555      solvers["clp"] = new OsiClpSolverInterface ;
556#   else
557      solvers["clp"] = 0 ;
558#   endif
559#   ifdef COIN_USE_DYLP
560#     ifndef CBC_DEFAULT_SOLVER
561#       define CBC_DEFAULT_SOLVER "dylp"
562#     endif
563      solvers["dylp"] = new OsiDylpSolverInterface  ;
564#   else
565      solvers["dylp"] = 0 ;
566#   endif
567/*
568  If we don't have a default solver, we're deeply confused.
569*/
570    OsiSolverInterface *dflt_solver = solvers[CBC_DEFAULT_SOLVER] ;
571    if (dflt_solver)
572    { std::cout << "Default solver is " << CBC_DEFAULT_SOLVER << std::endl ; }
573    else
574    { std::cerr << "No solvers! Aborting." << std::endl ;
575      return (1) ; }
576/*
577  For calculating run time.
578*/
579    double time1 = CoinCpuTime() ;
580    double time2;
581/*
582  Establish the command line interface: parameters, with associated info
583  messages, ranges, defaults. See CbcParam for details. Scan the vector of
584  solvers and add the names to the parameter object.
585*/
586    CbcParam parameters[MAXPARAMETERS];
587    int numberParameters ;
588    establishParams(numberParameters,parameters) ;
589
590    { int iSolver = 0 ;
591      for ( ; iSolver < numberParameters ; iSolver++)
592      { int match = parameters[iSolver].matches("solver") ;
593        if (match==1) break ; }
594      for (solverMapIter_t solverIter = solvers.begin() ;
595           solverIter != solvers.end() ;
596           solverIter++)
597      { if (solverIter->second)
598          parameters[iSolver].append(solverIter->first) ; }
599      int iKwd = parameters[iSolver].parameterOption(CBC_DEFAULT_SOLVER) ;
600      parameters[iSolver].setCurrentOption(iKwd) ; }
601/*
602  The rest of the default setup: establish a model, instantiate cut generators
603  and heuristics, set various default values.
604*/
605    dflt_solver->messageHandler()->setLogLevel(0) ;
606    CbcModel *model = new CbcModel(*dflt_solver) ;
607    model->messageHandler()->setLogLevel(1);
608    bool goodModel=false;
609   
610// Set up likely cut generators and defaults
611
612    CglGomory gomoryGen;
613    // try larger limit
614    gomoryGen.setLimit(3000);
615    // set default action (0=off,1=on,2=root)
616    int gomoryAction=1;
617
618    CglProbing probingGen;
619    probingGen.setUsingObjective(true);
620    probingGen.setMaxPass(3);
621    probingGen.setMaxProbe(100);
622    probingGen.setMaxLook(50);
623    probingGen.setRowCuts(3);
624    // set default action (0=off,1=on,2=root)
625    int probingAction=1;
626
627    CglKnapsackCover knapsackGen;
628    // set default action (0=off,1=on,2=root)
629    int knapsackAction=1;
630
631    CglOddHole oddholeGen;
632    oddholeGen.setMinimumViolation(0.005);
633    oddholeGen.setMinimumViolationPer(0.0002);
634    oddholeGen.setMaximumEntries(100);
635    // set default action (0=off,1=on,2=root)
636    int oddholeAction=0;
637
638    CglClique cliqueGen;
639    cliqueGen.setStarCliqueReport(false);
640    cliqueGen.setRowCliqueReport(false);
641    // set default action (0=off,1=on,2=root)
642    int cliqueAction=1;
643
644    CglMixedIntegerRounding mixedGen;
645    // set default action (0=off,1=on,2=root)
646    int mixedAction=1;
647
648    CglFlowCover flowGen;
649    // set default action (0=off,1=on,2=root)
650    int flowAction=1;
651
652    CglTwomir twomirGen;
653    // set default action (0=off,1=on,2=root)
654    int twomirAction=0;
655
656    bool useRounding=true;
657   
658    int allowImportErrors=0;
659    int keepImportNames=1;      // not implemented
660    int doScaling=1;
661    int preSolve=0;
662    double djFix=1.0e100;
663    double gapRatio=1.0e100;
664    double tightenFactor=0.0;
665
666    std::string directory ="./";
667    std::string field;
668/*
669  The main command parsing loop.
670*/
671    // total number of commands read
672    int numberGoodCommands=0;
673   
674    while (1) {
675      // next command
676      field=getCommand(argc,argv);
677     
678      // exit if null or similar
679      if (!field.length()) {
680        if (numberGoodCommands==1&&goodModel) {
681          // we just had file name
682          model->initialSolve();
683          model->solver()->messageHandler()->setLogLevel(0);
684          CbcRounding heuristic1(*model);
685          if (useRounding)
686            model->addHeuristic(&heuristic1) ;
687          // add cut generators if wanted
688          if (probingAction==1)
689            model->addCutGenerator(&probingGen,-1,"Probing");
690          else if (probingAction==2)
691            model->addCutGenerator(&probingGen,-99,"Probing");
692          if (gomoryAction==1)
693            model->addCutGenerator(&gomoryGen,-1,"Gomory");
694          else if (gomoryAction==2)
695            model->addCutGenerator(&gomoryGen,-99,"Gomory");
696          if (knapsackAction==1)
697            model->addCutGenerator(&knapsackGen,-1,"Knapsack");
698          else if (knapsackAction==2)
699            model->addCutGenerator(&knapsackGen,-99,"Knapsack");
700          if (oddholeAction==1)
701            model->addCutGenerator(&oddholeGen,-1,"OddHole");
702          else if (oddholeAction==2)
703            model->addCutGenerator(&oddholeGen,-99,"OddHole");
704          if (cliqueAction==1)
705            model->addCutGenerator(&cliqueGen,-1,"Clique");
706          else if (cliqueAction==2)
707            model->addCutGenerator(&cliqueGen,-99,"Clique");
708          if (mixedAction==1)
709            model->addCutGenerator(&mixedGen,-1,"MixedintegerRounding");
710          else if (mixedAction==2)
711            model->addCutGenerator(&mixedGen,-99,"MixedintegerRounding");
712          if (flowAction==1)
713            model->addCutGenerator(&flowGen,-1,"FlowCover");
714          else if (flowAction==2)
715            model->addCutGenerator(&flowGen,-99,"FlowCover");
716          if (twomirAction==1)
717            model->addCutGenerator(&twomirGen,-1,"TwoMirCuts");
718          else if (twomirAction==2)
719            model->addCutGenerator(&twomirGen,-99,"TwoMirCuts");
720          // Say we want timings
721          int numberGenerators = model->numberCutGenerators();
722          int iGenerator;
723          for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
724            CbcCutGenerator * generator = model->cutGenerator(iGenerator);
725            generator->setTiming(true);
726          }
727          model->branchAndBound();
728          time2 = CoinCpuTime();
729          totalTime += time2-time1;
730          std::cout<<"Result "<<model->getObjValue()<<
731            " iterations "<<model->getIterationCount()<<
732            " nodes "<<model->getNodeCount()<<
733            " took "<<time2-time1<<" seconds - total "<<totalTime<<std::endl;
734        } else if (!numberGoodCommands) {
735          // let's give the sucker a hint
736          std::cout
737            <<"Cbc takes input from arguments ( - switches to stdin)"
738            <<std::endl
739            <<"Enter ? for list of commands or (-)unitTest or -miplib"
740            <<" for tests"<<std::endl;
741        }
742        break;
743      }
744     
745      // see if ? at end
746      int numberQuery=0;
747      if (field!="?") {
748        int length = field.length();
749        int i;
750        for (i=length-1;i>0;i--) {
751          if (field[i]=='?') 
752            numberQuery++;
753          else
754            break;
755        }
756        field=field.substr(0,length-numberQuery);
757      }
758      // find out if valid command
759      int iParam;
760      int numberMatches=0;
761      for ( iParam=0; iParam<numberParameters; iParam++ ) {
762        int match = parameters[iParam].matches(field);
763        if (match==1) {
764          numberMatches = 1;
765          break;
766        } else {
767          numberMatches += match>>1;
768        }
769      }
770      if (iParam<numberParameters&&!numberQuery) {
771        // found
772        CbcParam found = parameters[iParam];
773        CbcParameterType type = found.type();
774        int valid;
775        numberGoodCommands++;
776        if (type==GENERALQUERY) {
777          std::cout<<"In argument list keywords have leading - "
778            ", -stdin or just - switches to stdin"<<std::endl;
779          std::cout<<"One command per line (and no -)"<<std::endl;
780          std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
781          std::cout<<"abcd?? adds explanation, if only one fuller help(LATER)"<<std::endl;
782          std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
783          std::cout<<"abcd value or abcd = value sets value"<<std::endl;
784          std::cout<<"Commands are:"<<std::endl;
785          for ( iParam=0; iParam<numberParameters; iParam+=4 ) {
786            int i;
787            for (i=iParam;i<min(numberParameters,iParam+4);i++) 
788              std::cout<<parameters[i].matchName()<<"  ";
789            std::cout<<std::endl;
790          }
791        } else if (type<81) {
792          // get next field as double
793          double value = getDoubleField(argc,argv,&valid);
794          if (!valid) {
795            parameters[iParam].setDoubleParameter(*model,value);
796          } else if (valid==1) {
797            abort();
798          } else {
799            std::cout<<parameters[iParam].name()<<" has value "<<
800              parameters[iParam].doubleParameter(*model)<<std::endl;
801          }
802        } else if (type<101) {
803          // get next field as double for local use
804          double value = getDoubleField(argc,argv,&valid);
805          if (!valid) {
806            if (!parameters[iParam].checkDoubleParameter(value)) {
807              switch(type) {
808              case DJFIX:
809                djFix=value;
810                preSolve=5;
811                break;
812              case GAPRATIO:
813                gapRatio=value;
814                break;
815              case TIGHTENFACTOR:
816                tightenFactor=value;
817                break;
818              default:
819                abort();
820              }
821            }
822          } else if (valid==1) {
823            abort();
824          } else {
825            switch(type) {
826            case DJFIX:
827              value = djFix ;
828              break;
829            case GAPRATIO:
830              value = gapRatio ;
831              break;
832            case TIGHTENFACTOR:
833              value = tightenFactor ;
834              break;
835            default:
836              abort();
837            }
838            std::cout << parameters[iParam].name() << " has value " <<
839                         value << std::endl ;
840          }
841        } else if (type<201) {
842          // get next field as int
843          int value = getIntField(argc,argv,&valid);
844          if (!valid) {
845            parameters[iParam].setIntParameter(*model,value);
846          } else if (valid==1) {
847            abort();
848          } else {
849            std::cout<<parameters[iParam].name()<<" has value "<<
850              parameters[iParam].intParameter(*model)<<std::endl;
851          }
852        } else if (type<301) {
853          // one of several strings
854          std::string value = getString(argc,argv);
855          int action = parameters[iParam].parameterOption(value);
856          if (action<0) {
857            if (value!="EOL") {
858              // no match
859              parameters[iParam].printOptions();
860            } else {
861              // print current value
862              std::cout<<parameters[iParam].name()<<" has value "<<
863                parameters[iParam].currentOption()<<std::endl;
864            }
865          } else {
866            parameters[iParam].setCurrentOption(action);
867            // for now hard wired
868            switch (type) {
869            case DIRECTION:
870              if (action==0)
871                model->solver()->setObjSense(1);
872              else
873                model->solver()->setObjSense(-1);
874              break;
875            case ERRORSALLOWED:
876              allowImportErrors = action;
877              break;
878            case KEEPNAMES:
879              keepImportNames = 1-action;
880              break;
881            case SCALING:
882              doScaling = 1-action;
883              break;
884            case GOMORYCUTS:
885              gomoryAction = action;
886              break;
887            case PROBINGCUTS:
888              probingAction = action;
889              break;
890            case KNAPSACKCUTS:
891              knapsackAction = action;
892              break;
893            case ODDHOLECUTS:
894              oddholeAction = action;
895              break;
896            case CLIQUECUTS:
897              cliqueAction = action;
898              break;
899            case FLOWCUTS:
900              flowAction = action;
901              break;
902            case MIXEDCUTS:
903              mixedAction = action;
904              break;
905            case TWOMIRCUTS:
906              twomirAction = action;
907              break;
908            case ROUNDING:
909              useRounding = action;
910              break;
911            case COSTSTRATEGY:
912              if (action!=1) {
913                printf("Pseudo costs not implemented yet\n");
914              } else {
915                int numberColumns = model->getNumCols();
916                int * sort = new int[numberColumns];
917                double * dsort = new double[numberColumns];
918                int * priority = new int [numberColumns];
919                const double * objective = model->getObjCoefficients();
920                int iColumn;
921                int n=0;
922                for (iColumn=0;iColumn<numberColumns;iColumn++) {
923                  if (model->isInteger(iColumn)) {
924                    sort[n]=n;
925                    dsort[n++]=-objective[iColumn];
926                  }
927                }
928                CoinSort_2(dsort,dsort+n,sort);
929                int level=0;
930                double last = -1.0e100;
931                for (int i=0;i<n;i++) {
932                  int iPut=sort[i];
933                  if (dsort[i]!=last) {
934                    level++;
935                    last=dsort[i];
936                  }
937                  priority[iPut]=level;
938                }
939                model->passInPriorities( priority,false);
940                delete [] priority;
941                delete [] sort;
942                delete [] dsort;
943              }
944              break;
945            case PRESOLVE:
946              preSolve = action*5;
947              break;
948            case SOLVER:
949            { for (int i = 0 ; i < (int) value.length() ; i++)
950                value[i] = tolower(value[i]) ;
951              OsiSolverInterface *newSolver = solvers[value]->clone() ;
952              model->assignSolver(newSolver) ;
953              std::cout << "Solver set to " << value << "." << std::endl ;
954              break ; }
955            default:
956            { std::cerr << "Unrecognized action. Aborting." << std::endl ;
957              abort(); }
958            }
959          }
960        } else {
961          // action
962          if (type==EXIT)
963            break; // stop all
964          switch (type) {
965          case SOLVECONTINUOUS:
966            if (goodModel) {
967              model->initialSolve();
968              time2 = CoinCpuTime();
969              totalTime += time2-time1;
970              std::cout<<"Result "<<model->solver()->getObjValue()<<
971                " iterations "<<model->solver()->getIterationCount()<<
972                " took "<<time2-time1<<" seconds - total "<<totalTime<<std::endl;
973              time1=time2;
974            } else {
975              std::cout<<"** Current model not valid"<<std::endl;
976            }
977            break;
978/*
979  Run branch-and-cut. First set a few options -- node comparison, scaling. If
980  the solver is Clp, consider running some presolve code (not yet converted
981  this to generic OSI) with branch-and-cut. If presolve is disabled, or the
982  solver is not Clp, simply run branch-and-cut. Print elapsed time at the end.
983*/
984          case BAB: // branchAndBound
985          { if (goodModel)
986            { CbcCompareUser compare; // Definition of node choice
987              model->setNodeComparison(compare);
988              OsiSolverInterface * solver = model->solver();
989              if (!doScaling)
990                solver->setHintParam(OsiDoScale,false,OsiHintTry);
991#ifdef COIN_USE_CLP
992              OsiClpSolverInterface * si =
993                dynamic_cast<OsiClpSolverInterface *>(solver) ;
994              if (preSolve&&si != NULL) {
995                // get clp itself
996                ClpSimplex * modelC = si->getModelPtr();
997                if (si->messageHandler()->logLevel())
998                  si->messageHandler()->setLogLevel(1);
999                if (modelC->tightenPrimalBounds()!=0) {
1000                  std::cout<<"Problem is infeasible!"<<std::endl;
1001                  break;
1002                }
1003                model->initialSolve();
1004                // bounds based on continuous
1005                if (tightenFactor) {
1006                  if (modelC->tightenPrimalBounds(tightenFactor)!=0) {
1007                    std::cout<<"Problem is infeasible!"<<std::endl;
1008                    break;
1009                  }
1010                }
1011                if (gapRatio<1.0e100) {
1012                  double value = si->getObjValue();
1013                  double value2 = gapRatio*(1.0e-5+fabs(value));
1014                  model->setAllowableGap(value2);
1015                  model->setAllowableFractionGap(gapRatio);
1016                  std::cout<<"Continuous "<<value
1017                           <<", so allowable gap set to "<<value2<<std::endl;
1018                }
1019                if (djFix<1.0e20) {
1020                  // do some fixing
1021                  int numberColumns = modelC->numberColumns();
1022                  int i;
1023                  const char * type = modelC->integerInformation();
1024                  double * lower = modelC->columnLower();
1025                  double * upper = modelC->columnUpper();
1026                  double * solution = modelC->primalColumnSolution();
1027                  double * dj = modelC->dualColumnSolution();
1028                  for (i=0;i<numberColumns;i++) {
1029                    if (type[i]) {
1030                      double value = solution[i];
1031                      if (value<lower[i]+1.0e-5&&dj[i]>djFix) {
1032                        solution[i]=lower[i];
1033                        upper[i]=lower[i];
1034                      } else if (value>upper[i]-1.0e-5&&dj[i]<-djFix) {
1035                        solution[i]=upper[i];
1036                        lower[i]=upper[i];
1037                      }
1038                    }
1039                  }
1040                }
1041                {
1042                  // integer presolve
1043                  CbcModel * model2 = model->integerPresolve();
1044                  if (model2) {
1045                    // Do complete search
1046                   
1047                    CbcRounding heuristic1(*model2);
1048                    if (useRounding)
1049                      model2->addHeuristic(&heuristic1);
1050                    model2->branchAndBound();
1051                    // get back solution
1052                    model->originalModel(model2,false);
1053                  } else {
1054                    // infeasible
1055                    exit(1);
1056                  }
1057                }
1058              } else
1059#endif
1060              { if (model->solver()->messageHandler()->logLevel())
1061                  model->solver()->messageHandler()->setLogLevel(1) ;
1062                model->initialSolve() ;
1063                if (gapRatio < 1.0e100)
1064                { double value = model->solver()->getObjValue() ;
1065                  double value2 = gapRatio*(1.0e-5+fabs(value)) ;
1066                  model->setAllowableGap(value2) ;
1067                  std::cout << "Continuous " << value
1068                            << ", so allowable gap set to "
1069                            << value2 << std::endl ; }
1070                CbcRounding heuristic1(*model) ;
1071                if (useRounding)
1072                  model->addHeuristic(&heuristic1) ;
1073                // add cut generators if wanted
1074                if (probingAction==1)
1075                  model->addCutGenerator(&probingGen,-1,"Probing");
1076                else if (probingAction==2)
1077                  model->addCutGenerator(&probingGen,-99,"Probing");
1078                if (gomoryAction==1)
1079                  model->addCutGenerator(&gomoryGen,-1,"Gomory");
1080                else if (gomoryAction==2)
1081                  model->addCutGenerator(&gomoryGen,-99,"Gomory");
1082                if (knapsackAction==1)
1083                  model->addCutGenerator(&knapsackGen,-1,"Knapsack");
1084                else if (knapsackAction==2)
1085                  model->addCutGenerator(&knapsackGen,-99,"Knapsack");
1086                if (oddholeAction==1)
1087                  model->addCutGenerator(&oddholeGen,-1,"OddHole");
1088                else if (oddholeAction==2)
1089                  model->addCutGenerator(&oddholeGen,-99,"OddHole");
1090                if (cliqueAction==1)
1091                  model->addCutGenerator(&cliqueGen,-1,"Clique");
1092                else if (cliqueAction==2)
1093                  model->addCutGenerator(&cliqueGen,-99,"Clique");
1094                if (mixedAction==1)
1095                  model->addCutGenerator(&mixedGen,-1,"MixedintegerRounding");
1096                else if (mixedAction==2)
1097                  model->addCutGenerator(&mixedGen,-99,"MixedintegerRounding");
1098                if (flowAction==1)
1099                  model->addCutGenerator(&flowGen,-1,"FlowCover");
1100                else if (flowAction==2)
1101                  model->addCutGenerator(&flowGen,-99,"FlowCover");
1102                if (twomirAction==1)
1103                  model->addCutGenerator(&twomirGen,-1,"TwoMirCuts");
1104                else if (twomirAction==2)
1105                  model->addCutGenerator(&twomirGen,-99,"TwoMirCuts");
1106                model->branchAndBound() ; }
1107              if (model->bestSolution())
1108              { std::cout << "Optimal solution "
1109                          << model->solver()->getObjValue() << std::endl ; }
1110              else
1111              { std::cout << "No integer solution found." << std::endl ; }
1112                       
1113              time2 = CoinCpuTime() ;
1114              totalTime += time2-time1 ;
1115              std::cout << "Result " << model->solver()->getObjValue()
1116                        << " took " << time2-time1 << " seconds - total "
1117                        << totalTime << std::endl ;
1118              time1 = time2 ; }
1119            else
1120            { std::cout << "** Current model not valid" << std::endl ; }
1121            break ; }
1122          case IMPORT:
1123            {
1124              // get next field
1125              field = getString(argc,argv);
1126              std::string fileName;
1127              bool canOpen=false;
1128              if (field=="-") {
1129                // stdin
1130                canOpen=true;
1131                fileName = "-";
1132              } else {
1133                if (field[0]=='/'||field[0]=='~')
1134                  fileName = field;
1135                else
1136                  fileName = directory+field;
1137                FILE *fp=fopen(fileName.c_str(),"r");
1138                if (fp) {
1139                  // can open - lets go for it
1140                  fclose(fp);
1141                  canOpen=true;
1142                } else {
1143                  std::cout<<"Unable to open file "<<fileName<<std::endl;
1144                }
1145              }
1146              if (canOpen) {
1147                model->gutsOfDestructor();
1148                int status =model->solver()->readMps(fileName.c_str(),"");
1149                if (!status||(status>0&&allowImportErrors)) {
1150                  // I don't think there is any need for this but ..
1151                  //OsiWarmStartBasis allSlack;
1152                  goodModel=true;
1153                  //model->setBasis(allSlack);
1154                  time2 = CoinCpuTime();
1155                  totalTime += time2-time1;
1156                  time1=time2;
1157                } else {
1158                  // errors
1159                  std::cout<<"There were "<<status<<
1160                    " errors on input"<<std::endl;
1161                }
1162              }
1163            }
1164            break;
1165          case EXPORT:
1166            {
1167              // get next field
1168              field = getString(argc,argv);
1169              std::string fileName;
1170              bool canOpen=false;
1171              if (field[0]=='/'||field[0]=='~')
1172                fileName = field;
1173              else
1174                fileName = directory+field;
1175              FILE *fp=fopen(fileName.c_str(),"w");
1176              if (fp) {
1177                // can open - lets go for it
1178                fclose(fp);
1179                canOpen=true;
1180              } else {
1181                std::cout<<"Unable to open file "<<fileName<<std::endl;
1182              }
1183              if (canOpen) {
1184                model->solver()->writeMps(fileName.c_str(),"");
1185                time2 = CoinCpuTime();
1186                totalTime += time2-time1;
1187                time1=time2;
1188              }
1189            }
1190            break;
1191          case MAXIMIZE:
1192            model->solver()->setObjSense(-1);
1193            break;
1194          case MINIMIZE:
1195            model->solver()->setObjSense(1);
1196            break;
1197          case DIRECTORY:
1198          { directory = getString(argc,argv);
1199            if (directory[directory.length()-1] != '/')
1200              directory += '/' ;
1201            break ; }
1202          case STDIN:
1203            read_mode=-1;
1204            break;
1205          case VERSION:
1206            std::cout<<"Coin LP version "<<CBCVERSION
1207                     <<", build "<<__DATE__<<std::endl;
1208            break;
1209          case UNITTEST:
1210            {
1211              // okay so there is not a real unit test
1212
1213              int status =model->solver()->readMps("../Mps/Sample/p0033.mps",
1214                                                   "");
1215              assert(!status);
1216              model->branchAndBound();
1217              model->solver()->resolve();
1218              std::cout<<"Optimal solution "<<model->solver()->getObjValue()<<std::endl;
1219              assert(fabs(model->solver()->getObjValue()-3089.0)<1.0e-5);
1220              fprintf(stderr,"Test was okay\n");
1221              status =model->solver()->readMps("../Mps/Sample/p0033.mps",
1222                                                   "");
1223              assert(!status);
1224              model->setCutoff(1.0e20);
1225              model->setBestObjectiveValue(1.0e20);
1226              model->setMaximumSolutions(1);
1227              model->setSolutionCount(0);
1228              // Switch off strong branching to give better chance of NOT finding best
1229              model->setNumberStrong(0);
1230              // Definition of node choice
1231              CbcCompareDefault compare(100.0);
1232              model->setNodeComparison(compare);
1233              model->solver()->resolve();
1234              model->branchAndBound();
1235              model->solver()->resolve();
1236              std::cout<<"partial solution "<<model->solver()->getObjValue()<<std::endl;
1237              if (model->solver()->getObjValue()<3090.0) {
1238                std::cout<<"Got optimal solution by mistake!"<<std::endl;
1239              }
1240            }
1241            break;
1242          case MIPLIB:
1243            {
1244              int mainTest (int argc, const char *argv[]);
1245              // create fields for test
1246              const char * fields[3];
1247              int nFields=1;
1248              fields[0]="fake main for miplib";
1249              if (directory!="./") {
1250                fields[1]=("-miplibDir="+directory).c_str();
1251                nFields=2;
1252              }
1253              mainTest(nFields,fields);
1254            }
1255            break;
1256          case SOLUTION:
1257            if (goodModel) {
1258              // get next field
1259              field = getString(argc,argv);
1260              std::string fileName;
1261              FILE *fp=NULL;
1262              if (field=="-"||field=="EOL") {
1263                // stdout
1264                fp=stdout;
1265              } else {
1266                if (field[0]=='/'||field[0]=='~')
1267                  fileName = field;
1268                else
1269                  fileName = directory+field;
1270                fp=fopen(fileName.c_str(),"w");
1271              }
1272              if (fp) {
1273                // make fancy later on
1274                int iRow;
1275                int numberRows=model->solver()->getNumRows();
1276                const double * dualRowSolution = model->getRowPrice();
1277                const double * primalRowSolution =  model->getRowActivity();
1278                for (iRow=0;iRow<numberRows;iRow++) {
1279                  fprintf(fp,"%7d ",iRow);
1280                  fprintf(fp,"%15.8g        %15.8g\n",primalRowSolution[iRow],
1281                          dualRowSolution[iRow]);
1282                }
1283                int iColumn;
1284                int numberColumns=model->solver()->getNumCols();
1285                const double * dualColumnSolution = 
1286                  model->getReducedCost();
1287                const double * primalColumnSolution = 
1288                  model->getColSolution();
1289                for (iColumn=0;iColumn<numberColumns;iColumn++) {
1290                  fprintf(fp,"%7d ",iColumn);
1291                  fprintf(fp,"%15.8g        %15.8g\n",
1292                          primalColumnSolution[iColumn],
1293                          dualColumnSolution[iColumn]);
1294                }
1295                if (fp!=stdout)
1296                  fclose(fp);
1297              } else {
1298                std::cout<<"Unable to open file "<<fileName<<std::endl;
1299              }
1300            } else {
1301              std::cout<<"** Current model not valid"<<std::endl;
1302             
1303            }
1304            break;
1305          default:
1306            abort();
1307          }
1308        } 
1309      } else if (!numberMatches) {
1310        std::cout<<"No match for "<<field<<" - ? for list of commands"
1311                 <<std::endl;
1312      } else if (numberMatches==1) {
1313        if (!numberQuery) {
1314          std::cout<<"Short match for "<<field<<" possible completion:"
1315                   <<std::endl;
1316          for ( iParam=0; iParam<numberParameters; iParam++ ) {
1317            int match = parameters[iParam].matches(field);
1318            if (match) 
1319              std::cout<<parameters[iParam].matchName()<<std::endl;
1320          }
1321        } else if (numberQuery) {
1322          std::cout<<"Short match for "<<field<<" completion:"
1323                   <<std::endl;
1324          for ( iParam=0; iParam<numberParameters; iParam++ ) {
1325            int match = parameters[iParam].matches(field);
1326            if (match) {
1327              std::cout<<parameters[iParam].matchName()<<" : ";
1328              std::cout<<parameters[iParam].shortHelp()<<std::endl;
1329            }
1330          }
1331        }
1332      } else {
1333        if (!numberQuery) 
1334          std::cout<<"Multiple matches for "<<field<<" - possible completions:"
1335                   <<std::endl;
1336        else
1337          std::cout<<"Completions of "<<field<<":"<<std::endl;
1338        for ( iParam=0; iParam<numberParameters; iParam++ ) {
1339          int match = parameters[iParam].matches(field);
1340          if (match) {
1341            std::cout<<parameters[iParam].matchName();
1342            if (numberQuery>=2) 
1343              std::cout<<" : "<<parameters[iParam].shortHelp();
1344            std::cout<<std::endl;
1345          }
1346        }
1347      }
1348    }
1349/*
1350  Final cleanup. Delete the model and the vector of available solvers.
1351*/
1352    delete model;
1353
1354    for (solverMapIter_t solverIter = solvers.begin() ;
1355         solverIter != solvers.end() ;
1356         solverIter++)
1357    { if (solverIter->second) delete solverIter->second ; }
1358  }
1359  return 0;
1360}   
Note: See TracBrowser for help on using the repository browser.