source: trunk/Samples/driver2.cpp @ 250

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

For Edwin

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.5 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpSimplex.hpp"
5#include "ClpPresolve.hpp"
6#include "CoinHelperFunctions.hpp"
7#include "CoinTime.hpp"
8#include "ClpDualRowSteepest.hpp"
9#include "ClpPrimalColumnSteepest.hpp"
10#include <iomanip>
11// This driver shows how to trap messages - this is just as in unitTest.cpp
12// ****** THis code is similar to MyMessageHandler.hpp and MyMessagehandler.cpp
13#include "CoinMessageHandler.hpp"
14
15/** This just adds a model to CoinMessage and a void pointer so
16    user can trap messages and do useful stuff. 
17    This is used in Clp/Test/unitTest.cpp
18
19    The file pointer is just there as an example of user stuff.
20
21*/
22class ClpSimplex;
23
24class MyMessageHandler : public CoinMessageHandler {
25 
26public:
27  /**@name Overrides */
28  //@{
29  virtual int print();
30  //@}
31  /**@name set and get */
32  //@{
33  /// Model
34  const ClpSimplex * model() const;
35  void setModel(ClpSimplex * model);
36  //@}
37
38  /**@name Constructors, destructor */
39  //@{
40  /** Default constructor. */
41  MyMessageHandler();
42  /// Constructor with pointer to model
43  MyMessageHandler(ClpSimplex * model,
44                           FILE * userPointer=NULL);
45  /** Destructor */
46  virtual ~MyMessageHandler();
47  //@}
48
49  /**@name Copy method */
50  //@{
51  /** The copy constructor. */
52  MyMessageHandler(const MyMessageHandler&);
53  /** The copy constructor from an CoinSimplexMessageHandler. */
54  MyMessageHandler(const CoinMessageHandler&);
55 
56  MyMessageHandler& operator=(const MyMessageHandler&);
57  /// Clone
58  virtual CoinMessageHandler * clone() const ;
59  //@}
60   
61   
62protected:
63  /**@name Data members
64     The data members are protected to allow access for derived classes. */
65  //@{
66  /// Pointer back to model
67  ClpSimplex * model_;
68  //@}
69};
70
71
72//#############################################################################
73// Constructors / Destructor / Assignment
74//#############################################################################
75
76//-------------------------------------------------------------------
77// Default Constructor
78//-------------------------------------------------------------------
79MyMessageHandler::MyMessageHandler () 
80  : CoinMessageHandler(),
81    model_(NULL)
82{
83}
84
85//-------------------------------------------------------------------
86// Copy constructor
87//-------------------------------------------------------------------
88MyMessageHandler::MyMessageHandler (const MyMessageHandler & rhs) 
89: CoinMessageHandler(rhs),
90    model_(rhs.model_)
91{ 
92}
93
94MyMessageHandler::MyMessageHandler (const CoinMessageHandler & rhs) 
95  : CoinMessageHandler(),
96    model_(NULL)
97{ 
98}
99
100// Constructor with pointer to model
101MyMessageHandler::MyMessageHandler(ClpSimplex * model,
102               FILE * userPointer)
103  : CoinMessageHandler(),
104    model_(model)
105{
106}
107
108//-------------------------------------------------------------------
109// Destructor
110//-------------------------------------------------------------------
111MyMessageHandler::~MyMessageHandler ()
112{
113}
114
115//----------------------------------------------------------------
116// Assignment operator
117//-------------------------------------------------------------------
118MyMessageHandler &
119MyMessageHandler::operator=(const MyMessageHandler& rhs)
120{
121  if (this != &rhs) {
122    CoinMessageHandler::operator=(rhs);
123    model_ = rhs.model_;
124  }
125  return *this;
126}
127//-------------------------------------------------------------------
128// Clone
129//-------------------------------------------------------------------
130CoinMessageHandler * MyMessageHandler::clone() const
131{
132  return new MyMessageHandler(*this);
133}
134// Print out values from first 20 messages
135static int times=0;
136int 
137MyMessageHandler::print()
138{
139  // You could have added a callback flag if you had wanted - see Clp_C_Interface.c
140  times++;
141  if (times<=20) {
142    int messageNumber = currentMessage().externalNumber();
143    if (currentSource()!="Clp")
144      messageNumber += 1000000;
145    int i;
146    int nDouble=numberDoubleFields();
147    printf("%d doubles - ",nDouble);
148    for (i=0;i<nDouble;i++)
149      printf("%g ",doubleValue(i));
150    printf("\n");;
151    int nInt=numberIntFields();
152    printf("%d ints - ",nInt);
153    for (i=0;i<nInt;i++)
154      printf("%d ",intValue(i));
155    printf("\n");;
156    int nString=numberStringFields();
157    printf("%d strings - ",nString);
158    for (i=0;i<nString;i++)
159      printf("%s ",stringValue(i).c_str());
160    printf("\n");;
161  }
162  return CoinMessageHandler::print();
163}
164const ClpSimplex *
165MyMessageHandler::model() const
166{
167  return model_;
168}
169void 
170MyMessageHandler::setModel(ClpSimplex * model)
171{
172  model_ = model;
173}
174
175int main (int argc, const char *argv[])
176{
177  ClpSimplex  model;
178  // Message handler
179  MyMessageHandler messageHandler(&model);
180  std::cout<<"Testing derived message handler"<<std::endl;
181  model.passInMessageHandler(&messageHandler);
182  int status;
183  // Keep names when reading an mps file
184  if (argc<2)
185    status=model.readMps("../../Mps/Sample/p0033.mps",true);
186  else
187    status=model.readMps(argv[1],true);
188
189  if (status) {
190    fprintf(stderr,"Bad readMps %s\n",argv[1]);
191    fprintf(stdout,"Bad readMps %s\n",argv[1]);
192    exit(1);
193  }
194
195  double time1 = CoinCpuTime();
196  /*
197    This driver shows how to do presolve.by hand (rather than with initialSolve)
198  */
199  ClpSimplex * model2;
200  ClpPresolve pinfo;
201  int numberPasses=5; // can change this
202  /* Use a tolerance of 1.0e-8 for feasibility, treat problem as
203     not being integer, do "numberpasses" passes and throw away names
204     in presolved model */
205  model2 = pinfo.presolvedModel(model,1.0e-8,false,numberPasses,false);
206  if (!model2) {
207    fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
208            argv[1],1.0e-8);
209    fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
210            argv[1],1.0e-8);
211    // model was infeasible - maybe try again with looser tolerances
212    model2 = pinfo.presolvedModel(model,1.0e-7,false,numberPasses,false);
213    if (!model2) {
214      fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
215              argv[1],1.0e-7);
216      fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
217              argv[1],1.0e-7);
218      exit(2);
219    }
220  }
221  // change factorization frequency from 200
222  model2->setFactorizationFrequency(100+model2->numberRows()/50);
223  if (argc<3 ||!strstr(argv[2],"primal")) {
224    // Use the dual algorithm unless user said "primal"
225    /* faster if bounds tightened as then dual can flip variables
226        to other bound to stay dual feasible.  We can trash the bounds as
227        this model is going to be thrown away
228    */
229    int numberInfeasibilities = model2->tightenPrimalBounds();
230    if (numberInfeasibilities)
231      std::cout<<"** Analysis indicates model infeasible"
232               <<std::endl;
233    model2->crash(1000.0,2);
234    ClpDualRowSteepest steep(1);
235    model2->setDualRowPivotAlgorithm(steep);
236    model2->dual();
237  } else {
238    ClpPrimalColumnSteepest steep(1);
239    model2->setPrimalColumnPivotAlgorithm(steep);
240    model2->primal();
241  }
242  pinfo.postsolve(true);
243
244  int numberIterations=model2->numberIterations();;
245  delete model2;
246  /* After this postsolve model should be optimal.
247     We can use checkSolution and test feasibility */
248  model.checkSolution();
249  if (model.numberDualInfeasibilities()||
250      model.numberPrimalInfeasibilities()) 
251    printf("%g dual %g(%d) Primal %g(%d)\n",
252           model.objectiveValue(),
253           model.sumDualInfeasibilities(),
254           model.numberDualInfeasibilities(),
255           model.sumPrimalInfeasibilities(),
256           model.numberPrimalInfeasibilities());
257  // But resolve for safety
258  model.primal(1);
259
260  numberIterations += model.numberIterations();;
261  // for running timing tests
262  std::cout<<argv[1]<<" Objective "<<model.objectiveValue()<<" took "<<
263    numberIterations<<" iterations and "<<
264    CoinCpuTime()-time1<<" seconds"<<std::endl;
265
266  std::string modelName;
267  model.getStrParam(ClpProbName,modelName);
268  std::cout<<"Model "<<modelName<<" has "<<model.numberRows()<<" rows and "<<
269    model.numberColumns()<<" columns"<<std::endl;
270
271  // remove this to print solution
272
273  exit(0);
274
275  /*
276    Now to print out solution.  The methods used return modifiable
277    arrays while the alternative names return const pointers -
278    which is of course much more virtuous.
279
280    This version just does non-zero columns
281 
282   */
283#if 0
284  int numberRows = model.numberRows();
285
286  // Alternatively getRowActivity()
287  double * rowPrimal = model.primalRowSolution();
288  // Alternatively getRowPrice()
289  double * rowDual = model.dualRowSolution();
290  // Alternatively getRowLower()
291  double * rowLower = model.rowLower();
292  // Alternatively getRowUpper()
293  double * rowUpper = model.rowUpper();
294  // Alternatively getRowObjCoefficients()
295  double * rowObjective = model.rowObjective();
296   
297  // If we have not kept names (parameter to readMps) this will be 0
298  assert(model.lengthNames());
299
300  // Row names
301  const std::vector<std::string> * rowNames = model.rowNames();
302
303
304  int iRow;
305
306  std::cout<<"                       Primal          Dual         Lower         Upper        (Cost)"
307           <<std::endl;
308
309  for (iRow=0;iRow<numberRows;iRow++) {
310    double value;
311    std::cout<<std::setw(6)<<iRow<<" "<<std::setw(8)<<(*rowNames)[iRow];
312    value = rowPrimal[iRow];
313    if (fabs(value)<1.0e5)
314      std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
315    else
316      std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
317    value = rowDual[iRow];
318    if (fabs(value)<1.0e5)
319      std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
320    else
321      std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
322    value = rowLower[iRow];
323    if (fabs(value)<1.0e5)
324      std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
325    else
326      std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
327    value = rowUpper[iRow];
328    if (fabs(value)<1.0e5)
329      std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
330    else
331      std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
332    if (rowObjective) {
333      value = rowObjective[iRow];
334      if (fabs(value)<1.0e5)
335        std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
336      else
337        std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
338    }
339    std::cout<<std::endl;
340  }
341#endif
342  std::cout<<"--------------------------------------"<<std::endl;
343
344  // Columns
345
346  int numberColumns = model.numberColumns();
347
348  // Alternatively getColSolution()
349  double * columnPrimal = model.primalColumnSolution();
350  // Alternatively getReducedCost()
351  double * columnDual = model.dualColumnSolution();
352  // Alternatively getColLower()
353  double * columnLower = model.columnLower();
354  // Alternatively getColUpper()
355  double * columnUpper = model.columnUpper();
356  // Alternatively getObjCoefficients()
357  double * columnObjective = model.objective();
358   
359  // If we have not kept names (parameter to readMps) this will be 0
360  assert(model.lengthNames());
361
362  // Column names
363  const std::vector<std::string> * columnNames = model.columnNames();
364
365
366  int iColumn;
367 
368  std::cout<<"                       Primal          Dual         Lower         Upper          Cost"
369           <<std::endl;
370
371  for (iColumn=0;iColumn<numberColumns;iColumn++) {
372    double value;
373    value = columnPrimal[iColumn];
374    if (fabs(value)>1.0e-8) {
375      std::cout<<std::setw(6)<<iColumn<<" "<<std::setw(8)<<(*columnNames)[iColumn];
376      if (fabs(value)<1.0e5)
377        std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
378      else
379        std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
380      value = columnDual[iColumn];
381      if (fabs(value)<1.0e5)
382        std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
383      else
384        std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
385      value = columnLower[iColumn];
386      if (fabs(value)<1.0e5)
387        std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
388      else
389        std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
390      value = columnUpper[iColumn];
391      if (fabs(value)<1.0e5)
392        std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
393      else
394        std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
395      value = columnObjective[iColumn];
396      if (fabs(value)<1.0e5)
397        std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
398      else
399        std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
400     
401      std::cout<<std::endl;
402    }
403  }
404  std::cout<<"--------------------------------------"<<std::endl;
405
406  return 0;
407}   
Note: See TracBrowser for help on using the repository browser.