source: trunk/Test/CbcMain.cpp @ 77

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

get version where I want it

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