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

Last change on this file since 1552 was 1552, checked in by mjs, 10 years ago

Format examples with 'astyle -A4 -p'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.4 KB
Line 
1/* $Id: driver2.cpp 1552 2010-05-25 01:03:59Z mjs $ */
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.