source: trunk/include/ClpSimplex.hpp @ 451

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

Changes to make a difficult problem faster (dual) + quadratic

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5   Authors
6 
7   John Forrest
8
9 */
10#ifndef ClpSimplex_H
11#define ClpSimplex_H
12
13#include <iostream>
14#include <cfloat>
15#include "ClpModel.hpp"
16#include "ClpMatrixBase.hpp"
17#include "ClpSolve.hpp"
18class ClpDualRowPivot;
19class ClpPrimalColumnPivot;
20class ClpFactorization;
21class CoinIndexedVector;
22class ClpNonLinearCost;
23class ClpSimplexProgress;
24
25/** This solves LPs using the simplex method
26
27    It inherits from ClpModel and all its arrays are created at
28    algorithm time. Originally I tried to work with model arrays
29    but for simplicity of coding I changed to single arrays with
30    structural variables then row variables.  Some coding is still
31    based on old style and needs cleaning up.
32
33    For a description of algorithms:
34
35    for dual see ClpSimplexDual.hpp and at top of ClpSimplexDual.cpp
36    for primal see ClpSimplexPrimal.hpp and at top of ClpSimplexPrimal.cpp
37
38    There is an algorithm data member.  + for primal variations
39    and - for dual variations
40
41    This file also includes (at end) a very simple class ClpSimplexProgress
42    which is where anti-looping stuff should migrate to
43
44*/
45
46class ClpSimplex : public ClpModel {
47   friend void ClpSimplexUnitTest(const std::string & mpsDir,
48                                  const std::string & netlibDir);
49
50public:
51  /** enums for status of various sorts.
52      First 4 match CoinWarmStartBasis,
53      isFixed means fixed at lower bound and out of basis
54  */
55  enum Status {
56    isFree = 0x00,
57    basic = 0x01,
58    atUpperBound = 0x02,
59    atLowerBound = 0x03,
60    superBasic = 0x04,
61    isFixed = 0x05
62  };
63  // For Dual
64  enum FakeBound {
65    noFake = 0x00,
66    bothFake = 0x01,
67    upperFake = 0x02,
68    lowerFake = 0x03
69  };
70
71  /**@name Constructors and destructor and copy */
72  //@{
73  /// Default constructor
74    ClpSimplex (  );
75
76  /** Copy constructor. May scale depending on mode
77      -1 leave mode as is
78      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
79  */
80  ClpSimplex(const ClpSimplex & rhs, int scalingMode =-1);
81  /** Copy constructor from model. May scale depending on mode
82      -1 leave mode as is
83      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
84  */
85  ClpSimplex(const ClpModel & rhs, int scalingMode=-1);
86  /** Subproblem constructor.  A subset of whole model is created from the
87      row and column lists given.  The new order is given by list order and
88      duplicates are allowed.  Name and integer information can be dropped
89  */
90  ClpSimplex (const ClpModel * wholeModel,
91              int numberRows, const int * whichRows,
92              int numberColumns, const int * whichColumns,
93              bool dropNames=true, bool dropIntegers=true);
94  /** This constructor modifies original ClpSimplex and stores
95      original stuff in created ClpSimplex.  It is only to be used in
96      conjunction with originalModel */
97  ClpSimplex (ClpSimplex * wholeModel,
98              int numberColumns, const int * whichColumns);
99  /** This copies back stuff from miniModel and then deletes miniModel.
100      Only to be used with mini constructor */
101  void originalModel(ClpSimplex * miniModel);
102  /// Assignment operator. This copies the data
103    ClpSimplex & operator=(const ClpSimplex & rhs);
104  /// Destructor
105   ~ClpSimplex (  );
106  // Ones below are just ClpModel with some changes
107  /** Loads a problem (the constraints on the
108        rows are given by lower and upper bounds). If a pointer is 0 then the
109        following values are the default:
110        <ul>
111          <li> <code>colub</code>: all columns have upper bound infinity
112          <li> <code>collb</code>: all columns have lower bound 0
113          <li> <code>rowub</code>: all rows have upper bound infinity
114          <li> <code>rowlb</code>: all rows have lower bound -infinity
115          <li> <code>obj</code>: all variables have 0 objective coefficient
116        </ul>
117    */
118  void loadProblem (  const ClpMatrixBase& matrix,
119                     const double* collb, const double* colub,   
120                     const double* obj,
121                     const double* rowlb, const double* rowub,
122                      const double * rowObjective=NULL);
123  void loadProblem (  const CoinPackedMatrix& matrix,
124                     const double* collb, const double* colub,   
125                     const double* obj,
126                     const double* rowlb, const double* rowub,
127                      const double * rowObjective=NULL);
128
129  /** Just like the other loadProblem() method except that the matrix is
130        given in a standard column major ordered format (without gaps). */
131  void loadProblem (  const int numcols, const int numrows,
132                     const CoinBigIndex* start, const int* index,
133                     const double* value,
134                     const double* collb, const double* colub,   
135                     const double* obj,
136                      const double* rowlb, const double* rowub,
137                      const double * rowObjective=NULL);
138  /// This one is for after presolve to save memory
139  void loadProblem (  const int numcols, const int numrows,
140                     const CoinBigIndex* start, const int* index,
141                      const double* value,const int * length,
142                     const double* collb, const double* colub,   
143                     const double* obj,
144                      const double* rowlb, const double* rowub,
145                      const double * rowObjective=NULL);
146  /// Read an mps file from the given filename
147  int readMps(const char *filename,
148              bool keepNames=false,
149              bool ignoreErrors = false);
150  /** Borrow model.  This is so we dont have to copy large amounts
151      of data around.  It assumes a derived class wants to overwrite
152      an empty model with a real one - while it does an algorithm.
153      This is same as ClpModel one, but sets scaling on etc. */
154  void borrowModel(ClpModel & otherModel);
155  void borrowModel(ClpSimplex & otherModel);
156   /// Pass in Event handler (cloned and deleted at end)
157   void passInEventHandler(const ClpEventHandler * eventHandler);
158  //@}
159
160  /**@name Functions most useful to user */
161  //@{
162  /** General solve algorithm which can do presolve.
163      See  ClpSolve.hpp for options
164   */
165  int initialSolve(ClpSolve & options);
166  /// Default initial solve
167  int initialSolve();
168  /// Dual initial solve
169  int initialDualSolve();
170  /// Primal initial solve
171  int initialPrimalSolve();
172  /** Dual algorithm - see ClpSimplexDual.hpp for method.
173      ifValuesPass==2 just does values pass and then stops */
174  int dual(int ifValuesPass=0);
175  /** Primal algorithm - see ClpSimplexPrimal.hpp for method.
176      ifValuesPass==2 just does values pass and then stops */
177  int primal(int ifValuesPass=0);
178  /** Solves nonlinear problem using SLP - may be used as crash
179      for other algorithms when number of iterations small.
180      Also exits if all problematical variables are changing
181      less than deltaTolerance
182  */
183  int nonlinearSLP(int numberPasses,double deltaTolerance);
184  /** Solves using barrier (assumes you have good cholesky factor code).
185      Does crossover to simplex if asked*/
186  int barrier(bool crossover=true);
187  /** Solves non-linear using reduced gradient.  Phase = 0 get feasible,
188      =1 use solution */
189  int reducedGradient(int phase=0);
190  /** Dual ranging.
191      This computes increase/decrease in cost for each given variable and corresponding
192      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
193      and numberColumns.. for artificials/slacks.
194      For non-basic variables the sequence number will be that of the non-basic variables.
195      The increase/decrease value is always >= 0.0
196
197      Up to user to provide correct length arrays.
198
199      Returns non-zero if infeasible unbounded etc
200  */
201  int dualRanging(int numberCheck,const int * which,
202                  double * costIncrease, int * sequenceIncrease,
203                  double * costDecrease, int * sequenceDecrease);
204  /** Primal ranging.
205      This computes increase/decrease in value for each given variable and corresponding
206      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
207      and numberColumns.. for artificials/slacks.
208      For basic variables the sequence number will be that of the basic variables.
209
210      Up to user to providen correct length arrays.
211
212      Returns non-zero if infeasible unbounded etc
213  */
214  int primalRanging(int numberCheck,const int * which,
215                  double * valueIncrease, int * sequenceIncrease,
216                  double * valueDecrease, int * sequenceDecrease);
217  /// Passes in factorization
218  void setFactorization( ClpFactorization & factorization);
219  /** Tightens primal bounds to make dual faster.  Unless
220      fixed, bounds are slightly looser than they could be.
221      This is to make dual go faster and is probably not needed
222      with a presolve.  Returns non-zero if problem infeasible.
223
224      Fudge for branch and bound - put bounds on columns of factor *
225      largest value (at continuous) - should improve stability
226      in branch and bound on infeasible branches (0.0 is off)
227  */
228  int tightenPrimalBounds(double factor=0.0);
229  /** Crash - at present just aimed at dual, returns
230      -2 if dual preferred and crash basis created
231      -1 if dual preferred and all slack basis preferred
232       0 if basis going in was not all slack
233       1 if primal preferred and all slack basis preferred
234       2 if primal preferred and crash basis created.
235       
236       if gap between bounds <="gap" variables can be flipped
237       ( If pivot -1 then can be made super basic!)
238
239       If "pivot" is
240       -1 No pivoting - always primal
241       0 No pivoting (so will just be choice of algorithm)
242       1 Simple pivoting e.g. gub
243       2 Mini iterations
244  */
245  int crash(double gap,int pivot);
246  /// Sets row pivot choice algorithm in dual
247  void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
248  /// Sets column pivot choice algorithm in primal
249  void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
250  /** For strong branching.  On input lower and upper are new bounds
251      while on output they are change in objective function values
252      (>1.0e50 infeasible).
253      Return code is 0 if nothing interesting, -1 if infeasible both
254      ways and +1 if infeasible one way (check values to see which one(s))
255      Solutions are filled in as well - even down, odd up - also
256      status and number of iterations
257  */
258  int strongBranching(int numberVariables,const int * variables,
259                      double * newLower, double * newUpper,
260                      double ** outputSolution,
261                      int * outputStatus, int * outputIterations,
262                      bool stopOnFirstInfeasible=true,
263                      bool alwaysFinish=false);
264  //@}
265
266  /**@name Needed for functionality of OsiSimplexInterface */
267  //@{
268  /** Pivot in a variable and out a variable.  Returns 0 if okay,
269      1 if inaccuracy forced re-factorization, -1 if would be singular.
270      Also updates primal/dual infeasibilities.
271      Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
272  */
273  int pivot();
274
275  /** Pivot in a variable and choose an outgoing one.  Assumes primal
276      feasible - will not go through a bound.  Returns step length in theta
277      Returns ray in ray_ (or NULL if no pivot)
278      Return codes as before but -1 means no acceptable pivot
279  */
280  int primalPivotResult();
281 
282  /** Pivot out a variable and choose an incoing one.  Assumes dual
283      feasible - will not go through a reduced cost. 
284      Returns step length in theta
285      Returns ray in ray_ (or NULL if no pivot)
286      Return codes as before but -1 means no acceptable pivot
287  */
288  int dualPivotResult();
289
290  /** Common bits of coding for dual and primal.  Return s0 if okay,
291      1 if bad matrix, 2 if very bad factorization
292  */
293  int startup(int ifValuesPass);
294  void finish();
295 
296  /** Factorizes and returns true if optimal.  Used by user */
297  bool statusOfProblem(bool initial=false);
298  //@}
299
300  /**@name most useful gets and sets */
301  //@{
302  /// If problem is primal feasible
303  inline bool primalFeasible() const
304         { return (numberPrimalInfeasibilities_==0);};
305  /// If problem is dual feasible
306  inline bool dualFeasible() const
307         { return (numberDualInfeasibilities_==0);};
308  /// factorization
309  inline ClpFactorization * factorization() const 
310          { return factorization_;};
311  /// Sparsity on or off
312  bool sparseFactorization() const;
313  void setSparseFactorization(bool value);
314  /// Factorization frequency
315  int factorizationFrequency() const;
316  void setFactorizationFrequency(int value);
317  /// Dual bound
318  inline double dualBound() const
319          { return dualBound_;};
320  void setDualBound(double value);
321  /// Infeasibility cost
322  inline double infeasibilityCost() const
323          { return infeasibilityCost_;};
324  void setInfeasibilityCost(double value);
325  /** Amount of print out:
326      0 - none
327      1 - just final
328      2 - just factorizations
329      3 - as 2 plus a bit more
330      4 - verbose
331      above that 8,16,32 etc just for selective debug
332  */
333  /** Perturbation:
334      50  - switch on perturbation
335      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
336      101 - we are perturbed
337      102 - don't try perturbing again
338      default is 100
339      others are for playing
340  */
341  inline int perturbation() const
342    { return perturbation_;};
343  void setPerturbation(int value);
344  /// Current (or last) algorithm
345  inline int algorithm() const 
346  {return algorithm_; } ;
347  /// Set algorithm
348  inline void setAlgorithm(int value)
349  {algorithm_=value; } ;
350  /// Sum of dual infeasibilities
351  inline double sumDualInfeasibilities() const 
352          { return sumDualInfeasibilities_;} ;
353  inline void setSumDualInfeasibilities(double value)
354          { sumDualInfeasibilities_=value;} ;
355  /// Sum of relaxed dual infeasibilities
356  inline double sumOfRelaxedDualInfeasibilities() const 
357          { return sumOfRelaxedDualInfeasibilities_;} ;
358  inline void setSumOfRelaxedDualInfeasibilities(double value)
359          { sumOfRelaxedDualInfeasibilities_=value;} ;
360  /// Number of dual infeasibilities
361  inline int numberDualInfeasibilities() const 
362          { return numberDualInfeasibilities_;} ;
363  inline void setNumberDualInfeasibilities(int value)
364          { numberDualInfeasibilities_=value;} ;
365  /// Sum of primal infeasibilities
366  inline double sumPrimalInfeasibilities() const 
367          { return sumPrimalInfeasibilities_;} ;
368  inline void setSumPrimalInfeasibilities(double value)
369          { sumPrimalInfeasibilities_=value;} ;
370  /// Sum of relaxed primal infeasibilities
371  inline double sumOfRelaxedPrimalInfeasibilities() const 
372          { return sumOfRelaxedPrimalInfeasibilities_;} ;
373  inline void setSumOfRelaxedPrimalInfeasibilities(double value)
374          { sumOfRelaxedPrimalInfeasibilities_=value;} ;
375  /// Number of primal infeasibilities
376  inline int numberPrimalInfeasibilities() const 
377          { return numberPrimalInfeasibilities_;} ;
378  inline void setNumberPrimalInfeasibilities(int value)
379          { numberPrimalInfeasibilities_=value;} ;
380  /** Save model to file, returns 0 if success.  This is designed for
381      use outside algorithms so does not save iterating arrays etc.
382  It does not save any messaging information.
383  Does not save scaling values.
384  It does not know about all types of virtual functions.
385  */
386  int saveModel(const char * fileName);
387  /** Restore model from file, returns 0 if success,
388      deletes current model */
389  int restoreModel(const char * fileName);
390 
391  /** Just check solution (for external use) - sets sum of
392      infeasibilities etc */
393  void checkSolution();
394  /// Useful row length arrays (0,1,2,3,4,5)
395  inline CoinIndexedVector * rowArray(int index) const
396  { return rowArray_[index];};
397  /// Useful column length arrays (0,1,2,3,4,5)
398  inline CoinIndexedVector * columnArray(int index) const
399  { return columnArray_[index];};
400  //@}
401
402  /******************** End of most useful part **************/
403  /**@name Functions less likely to be useful to casual user */
404  //@{
405  /** Given an existing factorization computes and checks
406      primal and dual solutions.  Uses input arrays for variables at
407      bounds.  Returns feasibility states */
408  int getSolution (  const double * rowActivities,
409                     const double * columnActivities);
410  /** Given an existing factorization computes and checks
411      primal and dual solutions.  Uses current problem arrays for
412      bounds.  Returns feasibility states */
413  int getSolution ();
414  /** Constructs a non linear cost from list of non-linearities (columns only)
415      First lower of each column is taken as real lower
416      Last lower is taken as real upper and cost ignored
417
418      Returns nonzero if bad data e.g. lowers not monotonic
419  */
420  int createPiecewiseLinearCosts(const int * starts,
421                   const double * lower, const double * gradient);
422  /** Return model - updates any scalars */
423  void returnModel(ClpSimplex & otherModel);
424  /** Factorizes using current basis. 
425      solveType - 1 iterating, 0 initial, -1 external
426      If 10 added then in primal values pass
427      Return codes are as from ClpFactorization unless initial factorization
428      when total number of singularities is returned
429  */
430  int internalFactorize(int solveType);
431  /// Save data
432  ClpDataSave saveData() ;
433  /// Restore data
434    void restoreData(ClpDataSave saved);
435  /// Clean up status
436  void cleanStatus();
437  /// Factorizes using current basis. For external use
438  int factorize();
439  /** Computes duals from scratch. If givenDjs then
440      allows for nonzero basic djs */
441  void computeDuals(double * givenDjs);
442  /// Computes primals from scratch
443  void computePrimals (  const double * rowActivities,
444                     const double * columnActivities);
445  /** Adds multiple of a column into an array */
446  void add(double * array,
447                   int column, double multiplier) const;
448  /**
449     Unpacks one column of the matrix into indexed array
450     Uses sequenceIn_
451     Also applies scaling if needed
452  */
453  void unpack(CoinIndexedVector * rowArray) const ;
454  /**
455     Unpacks one column of the matrix into indexed array
456     Slack if sequence>= numberColumns
457     Also applies scaling if needed
458  */
459  void unpack(CoinIndexedVector * rowArray,int sequence) const;
460  /**
461     Unpacks one column of the matrix into indexed array
462     ** as packed vector
463     Uses sequenceIn_
464     Also applies scaling if needed
465  */
466  void unpackPacked(CoinIndexedVector * rowArray) ;
467  /**
468     Unpacks one column of the matrix into indexed array
469     ** as packed vector
470     Slack if sequence>= numberColumns
471     Also applies scaling if needed
472  */
473  void unpackPacked(CoinIndexedVector * rowArray,int sequence);
474protected: 
475  /**
476      This does basis housekeeping and does values for in/out variables.
477      Can also decide to re-factorize
478  */
479  int housekeeping(double objectiveChange);
480  /** This sets largest infeasibility and most infeasible and sum
481      and number of infeasibilities (Primal) */
482  void checkPrimalSolution(const double * rowActivities=NULL,
483                           const double * columnActivies=NULL);
484  /** This sets largest infeasibility and most infeasible and sum
485      and number of infeasibilities (Dual) */
486  void checkDualSolution();
487public:
488  /** For advanced use.  When doing iterative solves things can get
489      nasty so on values pass if incoming solution has largest
490      infeasibility < incomingInfeasibility throw out variables
491      from basis until largest infeasibility < allowedInfeasibility
492      or incoming largest infeasibility.
493      If allowedInfeasibility>= incomingInfeasibility this is
494      always possible altough you may end up with an all slack basis.
495
496      Defaults are 1.0,10.0
497  */
498  void setValuesPassAction(float incomingInfeasibility,
499                           float allowedInfeasibility);
500  //@}
501  /**@name most useful gets and sets */
502  //@{
503  /// Worst column primal infeasibility
504  inline double columnPrimalInfeasibility() const 
505          { return columnPrimalInfeasibility_;} ;
506  /// Sequence of worst (-1 if feasible)
507  inline int columnPrimalSequence() const 
508          { return columnPrimalSequence_;} ;
509  /// Worst row primal infeasibility
510  inline double rowPrimalInfeasibility() const 
511          { return rowPrimalInfeasibility_;} ;
512  /// Sequence of worst (-1 if feasible)
513  inline int rowPrimalSequence() const 
514          { return rowPrimalSequence_;} ;
515  /** Worst column dual infeasibility (note - these may not be as meaningful
516      if the problem is primal infeasible */
517  inline double columnDualInfeasibility() const 
518          { return columnDualInfeasibility_;} ;
519  /// Sequence of worst (-1 if feasible)
520  inline int columnDualSequence() const 
521          { return columnDualSequence_;} ;
522  /// Worst row dual infeasibility
523  inline double rowDualInfeasibility() const 
524          { return rowDualInfeasibility_;} ;
525  /// Sequence of worst (-1 if feasible)
526  inline int rowDualSequence() const 
527          { return rowDualSequence_;} ;
528  /// Primal tolerance needed to make dual feasible (<largeTolerance)
529  inline double primalToleranceToGetOptimal() const 
530          { return primalToleranceToGetOptimal_;} ;
531  /// Remaining largest dual infeasibility
532  inline double remainingDualInfeasibility() const 
533          { return remainingDualInfeasibility_;} ;
534  /// Large bound value (for complementarity etc)
535  inline double largeValue() const 
536          { return largeValue_;} ;
537  void setLargeValue( double value) ;
538  /// Largest error on Ax-b
539  inline double largestPrimalError() const
540          { return largestPrimalError_;} ;
541  /// Largest error on basic duals
542  inline double largestDualError() const
543          { return largestDualError_;} ;
544  /// Largest difference between input primal solution and computed
545  inline double largestSolutionError() const
546          { return largestSolutionError_;} ;
547  /// Basic variables pivoting on which rows
548  inline int * pivotVariable() const
549          { return pivotVariable_;};
550  /// If automatic scaling on
551  inline bool automaticScaling() const
552  { return automaticScale_!=0;};
553  inline void setAutomaticScaling(bool onOff)
554  { automaticScale_ = onOff ? 1: 0;}; 
555  /// Current dual tolerance
556  inline double currentDualTolerance() const 
557          { return dualTolerance_;} ;
558  inline void setCurrentDualTolerance(double value)
559          { dualTolerance_ = value;} ;
560  /// Current primal tolerance
561  inline double currentPrimalTolerance() const 
562          { return primalTolerance_;} ;
563  inline void setCurrentPrimalTolerance(double value)
564          { primalTolerance_ = value;} ;
565  /// How many iterative refinements to do
566  inline int numberRefinements() const 
567          { return numberRefinements_;} ;
568  void setNumberRefinements( int value) ;
569  /// Alpha (pivot element) for use by classes e.g. steepestedge
570  inline double alpha() const { return alpha_;};
571  inline void setAlpha(double value) { alpha_ = value;};
572  /// Reduced cost of last incoming for use by classes e.g. steepestedge
573  inline double dualIn() const { return dualIn_;};
574  /// Pivot Row for use by classes e.g. steepestedge
575  inline int pivotRow() const{ return pivotRow_;};
576  inline void setPivotRow(int value) { pivotRow_=value;};
577  /// value of incoming variable (in Dual)
578  double valueIncomingDual() const;
579  //@}
580
581  protected:
582  /**@name protected methods */
583  //@{
584  /** May change basis and then returns number changed.
585      Computation of solutions may be overriden by given pi and solution
586  */
587  int gutsOfSolution ( double * givenDuals,
588                       const double * givenPrimals,
589                       bool valuesPass=false);
590  /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
591  void gutsOfDelete(int type);
592  /// Does most of copying
593  void gutsOfCopy(const ClpSimplex & rhs);
594  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
595      1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
596      8 bit does solution scaling in
597      16 bit does rowArray and columnArray indexed vectors
598      and makes row copy if wanted, also sets columnStart_ etc
599      Also creates scaling arrays if needed.  It does scaling if needed.
600      16 also moves solutions etc in to work arrays
601      On 16 returns false if problem "bad" i.e. matrix or bounds bad
602  */
603  bool createRim(int what,bool makeRowCopy=false);
604  /** releases above arrays and does solution scaling out.  May also
605      get rid of factorization data -
606      0 get rid of nothing, 1 get rid of arrays, 2 also factorization
607  */
608  void deleteRim(int getRidOfFactorizationData=2);
609  /// Sanity check on input rim data (after scaling) - returns true if okay
610  bool sanityCheck();
611  //@}
612  public:
613  /**@name public methods */
614  //@{
615  /** Return row or column sections - not as much needed as it
616      once was.  These just map into single arrays */
617  inline double * solutionRegion(int section) const
618  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
619  inline double * djRegion(int section) const
620  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
621  inline double * lowerRegion(int section) const
622  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
623  inline double * upperRegion(int section) const
624  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
625  inline double * costRegion(int section) const
626  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
627  /// Return region as single array
628  inline double * solutionRegion() const
629  { return solution_;};
630  inline double * djRegion() const
631  { return dj_;};
632  inline double * lowerRegion() const
633  { return lower_;};
634  inline double * upperRegion() const
635  { return upper_;};
636  inline double * costRegion() const
637  { return cost_;};
638  inline Status getStatus(int sequence) const
639  {return static_cast<Status> (status_[sequence]&7);};
640  inline void setStatus(int sequence, Status status)
641  {
642    unsigned char & st_byte = status_[sequence];
643    st_byte &= ~7;
644    st_byte |= status;
645  };
646  /** Normally the first factorization does sparse coding because
647      the factorization could be singular.  This allows initial dense
648      factorization when it is known to be safe
649  */
650  void setInitialDenseFactorization(bool onOff);
651  bool  initialDenseFactorization() const;
652  /** Return sequence In or Out */
653  inline int sequenceIn() const
654  {return sequenceIn_;};
655  inline int sequenceOut() const
656  {return sequenceOut_;};
657  /** Set sequenceIn or Out */
658  inline void  setSequenceIn(int sequence)
659  { sequenceIn_=sequence;};
660  inline void  setSequenceOut(int sequence)
661  { sequenceOut_=sequence;};
662  /** Return direction In or Out */
663  inline int directionIn() const
664  {return directionIn_;};
665  inline int directionOut() const
666  {return directionOut_;};
667  /** Set directionIn or Out */
668  inline void  setDirectionIn(int direction)
669  { directionIn_=direction;};
670  inline void  setDirectionOut(int direction)
671  { directionOut_=direction;};
672  /// Value of Out variable
673  inline double valueOut() const
674  { return valueOut_;};
675  /// Returns 1 if sequence indicates column
676  inline int isColumn(int sequence) const
677  { return sequence<numberColumns_ ? 1 : 0;};
678  /// Returns sequence number within section
679  inline int sequenceWithin(int sequence) const
680  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
681  /// Return row or column values
682  inline double solution(int sequence)
683  { return solution_[sequence];};
684  /// Return address of row or column values
685  inline double & solutionAddress(int sequence)
686  { return solution_[sequence];};
687  inline double reducedCost(int sequence)
688   { return dj_[sequence];};
689  inline double & reducedCostAddress(int sequence)
690   { return dj_[sequence];};
691  inline double lower(int sequence)
692  { return lower_[sequence];};
693  /// Return address of row or column lower bound
694  inline double & lowerAddress(int sequence)
695  { return lower_[sequence];};
696  inline double upper(int sequence)
697  { return upper_[sequence];};
698  /// Return address of row or column upper bound
699  inline double & upperAddress(int sequence)
700  { return upper_[sequence];};
701  inline double cost(int sequence)
702  { return cost_[sequence];};
703  /// Return address of row or column cost
704  inline double & costAddress(int sequence)
705  { return cost_[sequence];};
706  /// Return original lower bound
707  inline double originalLower(int iSequence) const
708  { if (iSequence<numberColumns_) return columnLower_[iSequence]; else
709    return rowLower_[iSequence-numberColumns_];};
710  /// Return original lower bound
711  inline double originalUpper(int iSequence) const
712  { if (iSequence<numberColumns_) return columnUpper_[iSequence]; else
713    return rowUpper_[iSequence-numberColumns_];};
714  /// Theta (pivot change)
715  inline double theta() const
716  { return theta_;};
717  /// Return pointer to details of costs
718  inline ClpNonLinearCost * nonLinearCost() const
719  { return nonLinearCost_;};
720  //@}
721  /**@name status methods */
722  //@{
723  inline void setFakeBound(int sequence, FakeBound fakeBound)
724  {
725    unsigned char & st_byte = status_[sequence];
726    st_byte &= ~24;
727    st_byte |= fakeBound<<3;
728  };
729  inline FakeBound getFakeBound(int sequence) const
730  {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
731  inline void setRowStatus(int sequence, Status status)
732  {
733    unsigned char & st_byte = status_[sequence+numberColumns_];
734    st_byte &= ~7;
735    st_byte |= status;
736  };
737  inline Status getRowStatus(int sequence) const
738  {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
739  inline void setColumnStatus(int sequence, Status status)
740  {
741    unsigned char & st_byte = status_[sequence];
742    st_byte &= ~7;
743    st_byte |= status;
744  };
745  inline Status getColumnStatus(int sequence) const
746  {return static_cast<Status> (status_[sequence]&7);};
747  inline void setPivoted( int sequence)
748  { status_[sequence] |= 32;};
749  inline void clearPivoted( int sequence)
750  { status_[sequence] &= ~32; };
751  inline bool pivoted(int sequence) const
752  {return (((status_[sequence]>>5)&1)!=0);};
753  /// To flag a variable (not inline to allow for column generation)
754  void setFlagged( int sequence);
755  inline void clearFlagged( int sequence)
756  {
757    status_[sequence] &= ~64;
758  };
759  inline bool flagged(int sequence) const
760  {return ((status_[sequence]&64)!=0);};
761  /// To say row active in primal pivot row choice
762  inline void setActive( int iRow)
763  {
764    status_[iRow] |= 128;
765  };
766  inline void clearActive( int iRow)
767  {
768    status_[iRow] &= ~128;
769  };
770  inline bool active(int iRow) const
771  {return ((status_[iRow]&128)!=0);};
772  /** Set up status array (can be used by OsiClp).
773      Also can be used to set up all slack basis */
774  void createStatus() ;
775  /** Sets up all slack basis and resets solution to
776      as it was after initial load or readMps */
777  void allSlackBasis(bool resetSolution=false);
778   
779  /// So we know when to be cautious
780  inline int lastBadIteration() const
781  {return lastBadIteration_;};
782  /// Progress flag - at present 0 bit says artificials out
783  inline int progressFlag() const
784  {return progressFlag_;};
785  /// Force re-factorization early
786  inline void forceFactorization(int value)
787  { forceFactorization_ = value;};
788  /// Raw objective value (so always minimize in primal)
789  inline double rawObjectiveValue() const
790  { return objectiveValue_;};
791  /** Number of extra rows.  These are ones which will be dynamically created
792      each iteration.  This is for GUB but may have other uses.
793  */
794  inline int numberExtraRows() const
795  { return numberExtraRows_;};
796  /** Maximum number of basic variables - can be more than number of rows if GUB
797  */
798  inline int maximumBasic() const
799  { return maximumBasic_;};
800  /** For advanced options
801      1 - Don't keep changing infeasibility weight
802      2 - Keep nonLinearCost round solves
803      4 - Force outgoing variables to exact bound (primal)
804      8 - Safe to use dense initial factorization
805      16 -Just use basic variables for operation if column generation
806      32 -Clean up with primal before strong branching
807      64 -Treat problem as feasible until last minute (i.e. minimize infeasibilities)
808  */
809  inline unsigned int specialOptions() const
810  { return specialOptions_;};
811  inline void setSpecialOptions(unsigned int value)
812  { specialOptions_=value;};
813  //@}
814
815////////////////// data //////////////////
816protected:
817
818  /**@name data.  Many arrays have a row part and a column part.
819   There is a single array with both - columns then rows and
820   then normally two arrays pointing to rows and columns.  The
821   single array is the owner of memory
822  */
823  //@{
824  /// Worst column primal infeasibility
825  double columnPrimalInfeasibility_;
826  /// Worst row primal infeasibility
827  double rowPrimalInfeasibility_;
828  /// Sequence of worst (-1 if feasible)
829  int columnPrimalSequence_;
830  /// Sequence of worst (-1 if feasible)
831  int rowPrimalSequence_;
832  /// Worst column dual infeasibility
833  double columnDualInfeasibility_;
834  /// Worst row dual infeasibility
835  double rowDualInfeasibility_;
836  /// Sequence of worst (-1 if feasible)
837  int columnDualSequence_;
838  /// Sequence of worst (-1 if feasible)
839  int rowDualSequence_;
840  /// Primal tolerance needed to make dual feasible (<largeTolerance)
841  double primalToleranceToGetOptimal_;
842  /// Remaining largest dual infeasibility
843  double remainingDualInfeasibility_;
844  /// Large bound value (for complementarity etc)
845  double largeValue_;
846  /// Largest error on Ax-b
847  double largestPrimalError_;
848  /// Largest error on basic duals
849  double largestDualError_;
850  /// Largest difference between input primal solution and computed
851  double largestSolutionError_;
852  /// Dual bound
853  double dualBound_;
854  /// Alpha (pivot element)
855  double alpha_;
856  /// Theta (pivot change)
857  double theta_;
858  /// Lower Bound on In variable
859  double lowerIn_;
860  /// Value of In variable
861  double valueIn_;
862  /// Upper Bound on In variable
863  double upperIn_;
864  /// Reduced cost of In variable
865  double dualIn_;
866  /// Lower Bound on Out variable
867  double lowerOut_;
868  /// Value of Out variable
869  double valueOut_;
870  /// Upper Bound on Out variable
871  double upperOut_;
872  /// Infeasibility (dual) or ? (primal) of Out variable
873  double dualOut_;
874  /// Current dual tolerance for algorithm
875  double dualTolerance_;
876  /// Current primal tolerance for algorithm
877  double primalTolerance_;
878  /// Sum of dual infeasibilities
879  double sumDualInfeasibilities_;
880  /// Sum of primal infeasibilities
881  double sumPrimalInfeasibilities_;
882  /// Weight assigned to being infeasible in primal
883  double infeasibilityCost_;
884  /// Sum of Dual infeasibilities using tolerance based on error in duals
885  double sumOfRelaxedDualInfeasibilities_;
886  /// Sum of Primal infeasibilities using tolerance based on error in primals
887  double sumOfRelaxedPrimalInfeasibilities_;
888  /// Working copy of lower bounds (Owner of arrays below)
889  double * lower_;
890  /// Row lower bounds - working copy
891  double * rowLowerWork_;
892  /// Column lower bounds - working copy
893  double * columnLowerWork_;
894  /// Working copy of upper bounds (Owner of arrays below)
895  double * upper_;
896  /// Row upper bounds - working copy
897  double * rowUpperWork_;
898  /// Column upper bounds - working copy
899  double * columnUpperWork_;
900  /// Working copy of objective (Owner of arrays below)
901  double * cost_;
902  /// Row objective - working copy
903  double * rowObjectiveWork_;
904  /// Column objective - working copy
905  double * objectiveWork_;
906  /// Useful row length arrays
907  CoinIndexedVector * rowArray_[6];
908  /// Useful column length arrays
909  CoinIndexedVector * columnArray_[6];
910  /// Sequence of In variable
911  int sequenceIn_;
912  /// Direction of In, 1 going up, -1 going down, 0 not a clude
913  int directionIn_;
914  /// Sequence of Out variable
915  int sequenceOut_;
916  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
917  int directionOut_;
918  /// Pivot Row
919  int pivotRow_;
920  /// Last good iteration (immediately after a re-factorization)
921  int lastGoodIteration_;
922  /// Working copy of reduced costs (Owner of arrays below)
923  double * dj_;
924  /// Reduced costs of slacks not same as duals (or - duals)
925  double * rowReducedCost_;
926  /// Possible scaled reduced costs
927  double * reducedCostWork_;
928  /// Working copy of primal solution (Owner of arrays below)
929  double * solution_;
930  /// Row activities - working copy
931  double * rowActivityWork_;
932  /// Column activities - working copy
933  double * columnActivityWork_;
934  /// Number of dual infeasibilities
935  int numberDualInfeasibilities_;
936  /// Number of dual infeasibilities (without free)
937  int numberDualInfeasibilitiesWithoutFree_;
938  /// Number of primal infeasibilities
939  int numberPrimalInfeasibilities_;
940  /// How many iterative refinements to do
941  int numberRefinements_;
942  /// dual row pivot choice
943  ClpDualRowPivot * dualRowPivot_;
944  /// primal column pivot choice
945  ClpPrimalColumnPivot * primalColumnPivot_;
946  /// Basic variables pivoting on which rows
947  int * pivotVariable_;
948  /// factorization
949  ClpFactorization * factorization_;
950  /// Saved version of solution
951  double * savedSolution_;
952  /// Number of times code has tentatively thought optimal
953  int numberTimesOptimal_;
954  /// If change has been made (first attempt at stopping looping)
955  int changeMade_;
956  /// Algorithm >0 == Primal, <0 == Dual
957  int algorithm_;
958  /** Now for some reliability aids
959      This forces re-factorization early */
960  int forceFactorization_;
961  /** Perturbation:
962      -50 to +50 - perturb by this power of ten (-6 sounds good)
963      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
964      101 - we are perturbed
965      102 - don't try perturbing again
966      default is 100
967  */
968  int perturbation_;
969  /// Saved status regions
970  unsigned char * saveStatus_;
971  /** Very wasteful way of dealing with infeasibilities in primal.
972      However it will allow non-linearities and use of dual
973      analysis.  If it doesn't work it can easily be replaced.
974  */
975  ClpNonLinearCost * nonLinearCost_;
976  /** For advanced options
977      1 - Don't keep changing infeasibility weight
978      2 - Keep nonLinearCost round solves
979      4 - Force outgoing variables to exact bound (primal)
980      8 - Safe to use dense initial factorization
981      16 -Just use basic variables for operation
982  */
983  unsigned int specialOptions_;
984  /// So we know when to be cautious
985  int lastBadIteration_;
986  /// So we know when to open up again
987  int lastFlaggedIteration_;
988  /// Can be used for count of fake bounds (dual) or fake costs (primal)
989  int numberFake_;
990  /// Progress flag - at present 0 bit says artificials out, 1 free in
991  int progressFlag_;
992  /// First free/super-basic variable (-1 if none)
993  int firstFree_;
994  /** Number of extra rows.  These are ones which will be dynamically created
995      each iteration.  This is for GUB but may have other uses.
996  */
997  int numberExtraRows_;
998  /** Maximum number of basic variables - can be more than number of rows if GUB
999  */
1000  int maximumBasic_;
1001  /** For advanced use.  When doing iterative solves things can get
1002      nasty so on values pass if incoming solution has largest
1003      infeasibility < incomingInfeasibility throw out variables
1004      from basis until largest infeasibility < allowedInfeasibility.
1005      if allowedInfeasibility>= incomingInfeasibility this is
1006      always possible altough you may end up with an all slack basis.
1007
1008      Defaults are 1.0,10.0
1009  */
1010  float incomingInfeasibility_;
1011  float allowedInfeasibility_;
1012  /// Automatic scaling of objective and rhs and bounds
1013  int automaticScale_;
1014  /// For dealing with all issues of cycling etc
1015  ClpSimplexProgress * progress_;
1016  //@}
1017};
1018//#############################################################################
1019/** A function that tests the methods in the ClpSimplex class. The
1020    only reason for it not to be a member method is that this way it doesn't
1021    have to be compiled into the library. And that's a gain, because the
1022    library should be compiled with optimization on, but this method should be
1023    compiled with debugging.
1024
1025    It also does some testing of ClpFactorization class
1026 */
1027void
1028ClpSimplexUnitTest(const std::string & mpsDir,
1029                   const std::string & netlibDir);
1030
1031
1032/// For saving extra information to see if looping.
1033class ClpSimplexProgress {
1034
1035public:
1036
1037
1038  /**@name Constructors and destructor and copy */
1039  //@{
1040  /// Default constructor
1041    ClpSimplexProgress (  );
1042
1043  /// Constructor from model
1044    ClpSimplexProgress ( ClpSimplex * model );
1045
1046  /// Copy constructor.
1047  ClpSimplexProgress(const ClpSimplexProgress &);
1048
1049  /// Assignment operator. This copies the data
1050    ClpSimplexProgress & operator=(const ClpSimplexProgress & rhs);
1051  /// Destructor
1052   ~ClpSimplexProgress (  );
1053  //@}
1054
1055  /**@name Check progress */
1056  //@{
1057  /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
1058      >=0 if give up and use as problem status
1059  */
1060    int looping (  );
1061  /// Start check at beginning of whileIterating
1062  void startCheck();
1063  /// Returns cycle length in whileIterating
1064  int cycle(int in, int out,int wayIn,int wayOut); 
1065
1066  /// Returns previous objective (if -1) - current if (0)
1067  double lastObjective(int back=1) const;
1068  /// Set real primal infeasibility and move back
1069  void setInfeasibility(double value);
1070  /// Returns real primal infeasibility (if -1) - current if (0)
1071  double lastInfeasibility(int back=1) const;
1072  /// Modify objective e.g. if dual infeasible in dual
1073  void modifyObjective(double value);
1074  /// Returns previous iteration number (if -1) - current if (0)
1075  int lastIterationNumber(int back=1) const;
1076  /// clears all iteration numbers (to switch off panic)
1077  void clearIterationNumbers();
1078  /// Odd state
1079  inline void newOddState()
1080  { oddState_= - oddState_-1;};
1081  inline void endOddState()
1082  { oddState_=abs(oddState_);};
1083  inline void clearOddState() 
1084  { oddState_=0;};
1085  inline int oddState() const
1086  { return oddState_;};
1087  /// number of bad times
1088  inline int badTimes() const
1089  { return numberBadTimes_;};
1090  inline void clearBadTimes()
1091  { numberBadTimes_=0;};
1092
1093  //@}
1094  /**@name Data  */
1095#define CLP_PROGRESS 5
1096  //@{
1097  /// Objective values
1098  double objective_[CLP_PROGRESS];
1099  /// Sum of infeasibilities for algorithm
1100  double infeasibility_[CLP_PROGRESS];
1101  /// Sum of real primal infeasibilities for primal
1102  double realInfeasibility_[CLP_PROGRESS];
1103#define CLP_CYCLE 12
1104  /// For cycle checking
1105  //double obj_[CLP_CYCLE];
1106  int in_[CLP_CYCLE];
1107  int out_[CLP_CYCLE];
1108  char way_[CLP_CYCLE];
1109  /// Pointer back to model so we can get information
1110  ClpSimplex * model_;
1111  /// Number of infeasibilities
1112  int numberInfeasibilities_[CLP_PROGRESS];
1113  /// Iteration number at which occurred
1114  int iterationNumber_[CLP_PROGRESS];
1115  /// Number of times checked (so won't stop too early)
1116  int numberTimes_;
1117  /// Number of times it looked like loop
1118  int numberBadTimes_;
1119  /// If things are in an odd state
1120  int oddState_;
1121  //@}
1122};
1123#endif
Note: See TracBrowser for help on using the repository browser.