source: trunk/Clp/examples/driver2.cpp @ 1370

Last change on this file since 1370 was 1370, checked in by forrest, 11 years ago

add ids

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