source: trunk/Clp/examples/driver2.cpp

Last change on this file was 1662, checked in by lou, 8 years ago

Add EPL license notice in examples.

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