Ignore:
Timestamp:
Jan 29, 2010 9:25:07 AM (10 years ago)
Author:
forrest
Message:

moving sandbox stuff to trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpSimplexDual.hpp

    r1402 r1502  
    33// Corporation and others.  All Rights Reserved.
    44
    5 /* 
     5/*
    66   Authors
    7    
     7
    88   John Forrest
    99
     
    1616/** This solves LPs using the dual simplex method
    1717
    18     It inherits from ClpSimplex.  It has no data of its own and 
    19     is never created - only cast from a ClpSimplex object at algorithm time. 
     18    It inherits from ClpSimplex.  It has no data of its own and
     19    is never created - only cast from a ClpSimplex object at algorithm time.
    2020
    2121*/
     
    2525public:
    2626
    27   /**@name Description of algorithm */
    28   //@{
    29   /** Dual algorithm
    30 
    31       Method
    32 
    33      It tries to be a single phase approach with a weight of 1.0 being
    34      given to getting optimal and a weight of updatedDualBound_ being
    35      given to getting dual feasible.  In this version I have used the
    36      idea that this weight can be thought of as a fake bound.  If the
    37      distance between the lower and upper bounds on a variable is less
    38      than the feasibility weight then we are always better off flipping
    39      to other bound to make dual feasible.  If the distance is greater
    40      then we make up a fake bound updatedDualBound_ away from one bound.
    41      If we end up optimal or primal infeasible, we check to see if
    42      bounds okay.  If so we have finished, if not we increase updatedDualBound_
    43      and continue (after checking if unbounded). I am undecided about
    44      free variables - there is coding but I am not sure about it.  At
    45      present I put them in basis anyway.
    46 
    47      The code is designed to take advantage of sparsity so arrays are
    48      seldom zeroed out from scratch or gone over in their entirety.
    49      The only exception is a full scan to find outgoing variable for
    50      Dantzig row choice.  For steepest edge we keep an updated list
    51      of infeasibilities (actually squares). 
    52      On easy problems we don't need full scan - just
    53      pick first reasonable.
    54 
    55      One problem is how to tackle degeneracy and accuracy.  At present
    56      I am using the modification of costs which I put in OSL and some
    57      of what I think is the dual analog of Gill et al.
    58      I am still not sure of the exact details.
    59 
    60      The flow of dual is three while loops as follows:
    61 
    62      while (not finished) {
    63 
    64        while (not clean solution) {
    65 
    66           Factorize and/or clean up solution by flipping variables so
    67           dual feasible.  If looks finished check fake dual bounds.
    68           Repeat until status is iterating (-1) or finished (0,1,2)
     27    /**@name Description of algorithm */
     28    //@{
     29    /** Dual algorithm
     30
     31        Method
     32
     33       It tries to be a single phase approach with a weight of 1.0 being
     34       given to getting optimal and a weight of updatedDualBound_ being
     35       given to getting dual feasible.  In this version I have used the
     36       idea that this weight can be thought of as a fake bound.  If the
     37       distance between the lower and upper bounds on a variable is less
     38       than the feasibility weight then we are always better off flipping
     39       to other bound to make dual feasible.  If the distance is greater
     40       then we make up a fake bound updatedDualBound_ away from one bound.
     41       If we end up optimal or primal infeasible, we check to see if
     42       bounds okay.  If so we have finished, if not we increase updatedDualBound_
     43       and continue (after checking if unbounded). I am undecided about
     44       free variables - there is coding but I am not sure about it.  At
     45       present I put them in basis anyway.
     46
     47       The code is designed to take advantage of sparsity so arrays are
     48       seldom zeroed out from scratch or gone over in their entirety.
     49       The only exception is a full scan to find outgoing variable for
     50       Dantzig row choice.  For steepest edge we keep an updated list
     51       of infeasibilities (actually squares).
     52       On easy problems we don't need full scan - just
     53       pick first reasonable.
     54
     55       One problem is how to tackle degeneracy and accuracy.  At present
     56       I am using the modification of costs which I put in OSL and some
     57       of what I think is the dual analog of Gill et al.
     58       I am still not sure of the exact details.
     59
     60       The flow of dual is three while loops as follows:
     61
     62       while (not finished) {
     63
     64         while (not clean solution) {
     65
     66            Factorize and/or clean up solution by flipping variables so
     67        dual feasible.  If looks finished check fake dual bounds.
     68        Repeat until status is iterating (-1) or finished (0,1,2)
     69
     70         }
     71
     72         while (status==-1) {
     73
     74           Iterate until no pivot in or out or time to re-factorize.
     75
     76           Flow is:
     77
     78           choose pivot row (outgoing variable).  if none then
     79       we are primal feasible so looks as if done but we need to
     80       break and check bounds etc.
     81
     82       Get pivot row in tableau
     83
     84           Choose incoming column.  If we don't find one then we look
     85       primal infeasible so break and check bounds etc.  (Also the
     86       pivot tolerance is larger after any iterations so that may be
     87       reason)
     88
     89           If we do find incoming column, we may have to adjust costs to
     90       keep going forwards (anti-degeneracy).  Check pivot will be stable
     91       and if unstable throw away iteration and break to re-factorize.
     92       If minor error re-factorize after iteration.
     93
     94       Update everything (this may involve flipping variables to stay
     95       dual feasible.
     96
     97         }
    6998
    7099       }
    71100
    72        while (status==-1) {
    73 
    74          Iterate until no pivot in or out or time to re-factorize.
    75 
    76          Flow is:
    77 
    78          choose pivot row (outgoing variable).  if none then
    79          we are primal feasible so looks as if done but we need to
    80          break and check bounds etc.
    81 
    82          Get pivot row in tableau
    83 
    84          Choose incoming column.  If we don't find one then we look
    85          primal infeasible so break and check bounds etc.  (Also the
    86          pivot tolerance is larger after any iterations so that may be
    87          reason)
    88 
    89          If we do find incoming column, we may have to adjust costs to
    90          keep going forwards (anti-degeneracy).  Check pivot will be stable
    91          and if unstable throw away iteration and break to re-factorize.
    92          If minor error re-factorize after iteration.
    93 
    94          Update everything (this may involve flipping variables to stay
    95          dual feasible.
    96 
    97        }
    98 
    99      }
    100 
    101      TODO's (or maybe not)
    102 
    103      At present we never check we are going forwards.  I overdid that in
    104      OSL so will try and make a last resort.
    105 
    106      Needs partial scan pivot out option.
    107 
    108      May need other anti-degeneracy measures, especially if we try and use
    109      loose tolerances as a way to solve in fewer iterations.
    110 
    111      I like idea of dynamic scaling.  This gives opportunity to decouple
    112      different implications of scaling for accuracy, iteration count and
    113      feasibility tolerance.
    114 
    115      for use of exotic parameter startFinishoptions see Clpsimplex.hpp
    116   */
    117 
    118   int dual(int ifValuesPass,int startFinishOptions=0);
    119   /** For strong branching.  On input lower and upper are new bounds
    120       while on output they are change in objective function values
    121       (>1.0e50 infeasible).
    122       Return code is 0 if nothing interesting, -1 if infeasible both
    123       ways and +1 if infeasible one way (check values to see which one(s))
    124       Solutions are filled in as well - even down, odd up - also
    125       status and number of iterations
    126   */
    127   int strongBranching(int numberVariables,const int * variables,
    128                       double * newLower, double * newUpper,
    129                       double ** outputSolution,
    130                       int * outputStatus, int * outputIterations,
    131                       bool stopOnFirstInfeasible=true,
    132                       bool alwaysFinish=false,
    133                       int startFinishOptions=0);
    134   /// This does first part of StrongBranching
    135   ClpFactorization * setupForStrongBranching(char * arrays, int numberRows,
    136                                              int numberColumns,bool solveLp=false);
    137   /// This cleans up after strong branching
    138   void cleanupAfterStrongBranching(ClpFactorization * factorization);
    139   //@}
    140 
    141   /**@name Functions used in dual */
    142   //@{
    143   /** This has the flow between re-factorizations
    144       Broken out for clarity and will be used by strong branching
    145 
    146       Reasons to come out:
    147       -1 iterations etc
    148       -2 inaccuracy
    149       -3 slight inaccuracy (and done iterations)
    150       +0 looks optimal (might be unbounded - but we will investigate)
    151       +1 looks infeasible
    152       +3 max iterations
    153 
    154       If givenPi not NULL then in values pass
    155    */
    156   int whileIterating(double * & givenPi,int ifValuesPass);
    157   /** The duals are updated by the given arrays.
    158       Returns number of infeasibilities.
    159       After rowArray and columnArray will just have those which
    160       have been flipped.
    161       Variables may be flipped between bounds to stay dual feasible.
    162       The output vector has movement of primal
    163       solution (row length array) */
    164   int updateDualsInDual(CoinIndexedVector * rowArray,
    165                   CoinIndexedVector * columnArray,
    166                   CoinIndexedVector * outputArray,
    167                   double theta,
    168                   double & objectiveChange,
    169                         bool fullRecompute);
    170   /** The duals are updated by the given arrays.
    171       This is in values pass - so no changes to primal is made
    172   */
    173   void updateDualsInValuesPass(CoinIndexedVector * rowArray,
    174                   CoinIndexedVector * columnArray,
    175                   double theta);
    176   /** While updateDualsInDual sees what effect is of flip
    177       this does actual flipping.
    178   */
    179   void flipBounds(CoinIndexedVector * rowArray,
    180                   CoinIndexedVector * columnArray);
    181   /**
    182       Row array has row part of pivot row
    183       Column array has column part.
    184       This chooses pivot column.
    185       Spare arrays are used to save pivots which will go infeasible
    186       We will check for basic so spare array will never overflow.
    187       If necessary will modify costs
    188       For speed, we may need to go to a bucket approach when many
    189       variables are being flipped.
    190       Returns best possible pivot value
    191   */
    192   double dualColumn(CoinIndexedVector * rowArray,
    193                   CoinIndexedVector * columnArray,
    194                   CoinIndexedVector * spareArray,
    195                   CoinIndexedVector * spareArray2,
    196                   double accpetablePivot,
    197                   CoinBigIndex * dubiousWeights);
    198   /// Does first bit of dualColumn
    199   int dualColumn0(const CoinIndexedVector * rowArray,
    200                   const CoinIndexedVector * columnArray,
    201                   CoinIndexedVector * spareArray,
    202                   double acceptablePivot,
    203                   double & upperReturn, double &bestReturn,double & badFree);
    204   /**
    205       Row array has row part of pivot row
    206       Column array has column part.
    207       This sees what is best thing to do in dual values pass
    208       if sequenceIn==sequenceOut can change dual on chosen row and leave variable in basis
    209   */
    210   void checkPossibleValuesMove(CoinIndexedVector * rowArray,
    211                                CoinIndexedVector * columnArray,
    212                               double acceptablePivot);
    213   /**
    214       Row array has row part of pivot row
    215       Column array has column part.
    216       This sees what is best thing to do in branch and bound cleanup
    217       If sequenceIn_ < 0 then can't do anything
    218   */
    219   void checkPossibleCleanup(CoinIndexedVector * rowArray,
    220                                CoinIndexedVector * columnArray,
    221                               double acceptablePivot);
    222   /**
    223       This sees if we can move duals in dual values pass.
    224       This is done before any pivoting
    225   */
    226   void doEasyOnesInValuesPass(double * givenReducedCosts);
    227   /**
    228       Chooses dual pivot row
    229       Would be faster with separate region to scan
    230       and will have this (with square of infeasibility) when steepest
    231       For easy problems we can just choose one of the first rows we look at
    232      
    233       If alreadyChosen >=0 then in values pass and that row has been
    234       selected
    235   */
    236   void dualRow(int alreadyChosen);
    237   /** Checks if any fake bounds active - if so returns number and modifies
    238       updatedDualBound_ and everything.
    239       Free variables will be left as free
    240       Returns number of bounds changed if >=0
    241       Returns -1 if not initialize and no effect
    242       Fills in changeVector which can be used to see if unbounded
    243       and cost of change vector
    244       If 2 sets to original (just changed)
    245   */
    246   int changeBounds(int initialize,CoinIndexedVector * outputArray,
    247                    double & changeCost);
    248   /** As changeBounds but just changes new bounds for a single variable.
    249       Returns true if change */
    250   bool changeBound( int iSequence);
    251   /// Restores bound to original bound
    252   void originalBound(int iSequence);
    253   /** Checks if tentative optimal actually means unbounded in dual
    254       Returns -3 if not, 2 if is unbounded */
    255   int checkUnbounded(CoinIndexedVector * ray,CoinIndexedVector * spare,
    256                      double changeCost);
    257   /**  Refactorizes if necessary
    258        Checks if finished.  Updates status.
    259        lastCleaned refers to iteration at which some objective/feasibility
    260        cleaning too place.
    261 
    262        type - 0 initial so set up save arrays etc
    263             - 1 normal -if good update save
    264             - 2 restoring from saved
    265   */
    266   void statusOfProblemInDual(int & lastCleaned, int type,
    267                              double * givenDjs, ClpDataSave & saveData,
    268                              int ifValuesPass);
    269   /** Perturbs problem (method depends on perturbation())
    270       returns nonzero if should go to dual */
    271   int perturb();
    272   /** Fast iterations.  Misses out a lot of initialization.
    273       Normally stops on maximum iterations, first re-factorization
    274       or tentative optimum.  If looks interesting then continues as
    275       normal.  Returns 0 if finished properly, 1 otherwise.
    276   */
    277   int fastDual(bool alwaysFinish=false);
    278   /** Checks number of variables at fake bounds.  This is used by fastDual
    279       so can exit gracefully before end */
    280   int numberAtFakeBound();
    281 
    282   /** Pivot in a variable and choose an outgoing one.  Assumes dual
    283       feasible - will not go through a reduced cost.  Returns step length in theta
    284       Returns ray in ray_ (or NULL if no pivot)
    285       Return codes as before but -1 means no acceptable pivot
    286   */
    287   int pivotResult();
    288   /** Get next free , -1 if none */
    289   int nextSuperBasic();
    290   /** Startup part of dual (may be extended to other algorithms)
    291       returns 0 if good, 1 if bad */
    292   int startupSolve(int ifValuesPass,double * saveDuals,int startFinishOptions);
    293   void finishSolve(int startFinishOptions);
    294   void gutsOfDual(int ifValuesPass,double * & saveDuals,int initialStatus,
    295                   ClpDataSave & saveData);
    296   //int dual2(int ifValuesPass,int startFinishOptions=0);
    297   void resetFakeBounds(int type);
    298  
    299   //@}
     101       TODO's (or maybe not)
     102
     103       At present we never check we are going forwards.  I overdid that in
     104       OSL so will try and make a last resort.
     105
     106       Needs partial scan pivot out option.
     107
     108       May need other anti-degeneracy measures, especially if we try and use
     109       loose tolerances as a way to solve in fewer iterations.
     110
     111       I like idea of dynamic scaling.  This gives opportunity to decouple
     112       different implications of scaling for accuracy, iteration count and
     113       feasibility tolerance.
     114
     115       for use of exotic parameter startFinishoptions see Clpsimplex.hpp
     116    */
     117
     118    int dual(int ifValuesPass, int startFinishOptions = 0);
     119    /** For strong branching.  On input lower and upper are new bounds
     120        while on output they are change in objective function values
     121        (>1.0e50 infeasible).
     122        Return code is 0 if nothing interesting, -1 if infeasible both
     123        ways and +1 if infeasible one way (check values to see which one(s))
     124        Solutions are filled in as well - even down, odd up - also
     125        status and number of iterations
     126    */
     127    int strongBranching(int numberVariables, const int * variables,
     128                        double * newLower, double * newUpper,
     129                        double ** outputSolution,
     130                        int * outputStatus, int * outputIterations,
     131                        bool stopOnFirstInfeasible = true,
     132                        bool alwaysFinish = false,
     133                        int startFinishOptions = 0);
     134    /// This does first part of StrongBranching
     135    ClpFactorization * setupForStrongBranching(char * arrays, int numberRows,
     136            int numberColumns, bool solveLp = false);
     137    /// This cleans up after strong branching
     138    void cleanupAfterStrongBranching(ClpFactorization * factorization);
     139    //@}
     140
     141    /**@name Functions used in dual */
     142    //@{
     143    /** This has the flow between re-factorizations
     144        Broken out for clarity and will be used by strong branching
     145
     146        Reasons to come out:
     147        -1 iterations etc
     148        -2 inaccuracy
     149        -3 slight inaccuracy (and done iterations)
     150        +0 looks optimal (might be unbounded - but we will investigate)
     151        +1 looks infeasible
     152        +3 max iterations
     153
     154        If givenPi not NULL then in values pass
     155     */
     156    int whileIterating(double * & givenPi, int ifValuesPass);
     157    /** The duals are updated by the given arrays.
     158        Returns number of infeasibilities.
     159        After rowArray and columnArray will just have those which
     160        have been flipped.
     161        Variables may be flipped between bounds to stay dual feasible.
     162        The output vector has movement of primal
     163        solution (row length array) */
     164    int updateDualsInDual(CoinIndexedVector * rowArray,
     165                          CoinIndexedVector * columnArray,
     166                          CoinIndexedVector * outputArray,
     167                          double theta,
     168                          double & objectiveChange,
     169                          bool fullRecompute);
     170    /** The duals are updated by the given arrays.
     171        This is in values pass - so no changes to primal is made
     172    */
     173    void updateDualsInValuesPass(CoinIndexedVector * rowArray,
     174                                 CoinIndexedVector * columnArray,
     175                                 double theta);
     176    /** While updateDualsInDual sees what effect is of flip
     177        this does actual flipping.
     178    */
     179    void flipBounds(CoinIndexedVector * rowArray,
     180                    CoinIndexedVector * columnArray);
     181    /**
     182        Row array has row part of pivot row
     183        Column array has column part.
     184        This chooses pivot column.
     185        Spare arrays are used to save pivots which will go infeasible
     186        We will check for basic so spare array will never overflow.
     187        If necessary will modify costs
     188        For speed, we may need to go to a bucket approach when many
     189        variables are being flipped.
     190        Returns best possible pivot value
     191    */
     192    double dualColumn(CoinIndexedVector * rowArray,
     193                      CoinIndexedVector * columnArray,
     194                      CoinIndexedVector * spareArray,
     195                      CoinIndexedVector * spareArray2,
     196                      double accpetablePivot,
     197                      CoinBigIndex * dubiousWeights);
     198    /// Does first bit of dualColumn
     199    int dualColumn0(const CoinIndexedVector * rowArray,
     200                    const CoinIndexedVector * columnArray,
     201                    CoinIndexedVector * spareArray,
     202                    double acceptablePivot,
     203                    double & upperReturn, double &bestReturn, double & badFree);
     204    /**
     205        Row array has row part of pivot row
     206        Column array has column part.
     207        This sees what is best thing to do in dual values pass
     208        if sequenceIn==sequenceOut can change dual on chosen row and leave variable in basis
     209    */
     210    void checkPossibleValuesMove(CoinIndexedVector * rowArray,
     211                                 CoinIndexedVector * columnArray,
     212                                 double acceptablePivot);
     213    /**
     214        Row array has row part of pivot row
     215        Column array has column part.
     216        This sees what is best thing to do in branch and bound cleanup
     217        If sequenceIn_ < 0 then can't do anything
     218    */
     219    void checkPossibleCleanup(CoinIndexedVector * rowArray,
     220                              CoinIndexedVector * columnArray,
     221                              double acceptablePivot);
     222    /**
     223        This sees if we can move duals in dual values pass.
     224        This is done before any pivoting
     225    */
     226    void doEasyOnesInValuesPass(double * givenReducedCosts);
     227    /**
     228        Chooses dual pivot row
     229        Would be faster with separate region to scan
     230        and will have this (with square of infeasibility) when steepest
     231        For easy problems we can just choose one of the first rows we look at
     232
     233        If alreadyChosen >=0 then in values pass and that row has been
     234        selected
     235    */
     236    void dualRow(int alreadyChosen);
     237    /** Checks if any fake bounds active - if so returns number and modifies
     238        updatedDualBound_ and everything.
     239        Free variables will be left as free
     240        Returns number of bounds changed if >=0
     241        Returns -1 if not initialize and no effect
     242        Fills in changeVector which can be used to see if unbounded
     243        and cost of change vector
     244        If 2 sets to original (just changed)
     245    */
     246    int changeBounds(int initialize, CoinIndexedVector * outputArray,
     247                     double & changeCost);
     248    /** As changeBounds but just changes new bounds for a single variable.
     249        Returns true if change */
     250    bool changeBound( int iSequence);
     251    /// Restores bound to original bound
     252    void originalBound(int iSequence);
     253    /** Checks if tentative optimal actually means unbounded in dual
     254        Returns -3 if not, 2 if is unbounded */
     255    int checkUnbounded(CoinIndexedVector * ray, CoinIndexedVector * spare,
     256                       double changeCost);
     257    /**  Refactorizes if necessary
     258         Checks if finished.  Updates status.
     259         lastCleaned refers to iteration at which some objective/feasibility
     260         cleaning too place.
     261
     262         type - 0 initial so set up save arrays etc
     263              - 1 normal -if good update save
     264          - 2 restoring from saved
     265    */
     266    void statusOfProblemInDual(int & lastCleaned, int type,
     267                               double * givenDjs, ClpDataSave & saveData,
     268                               int ifValuesPass);
     269    /** Perturbs problem (method depends on perturbation())
     270        returns nonzero if should go to dual */
     271    int perturb();
     272    /** Fast iterations.  Misses out a lot of initialization.
     273        Normally stops on maximum iterations, first re-factorization
     274        or tentative optimum.  If looks interesting then continues as
     275        normal.  Returns 0 if finished properly, 1 otherwise.
     276    */
     277    int fastDual(bool alwaysFinish = false);
     278    /** Checks number of variables at fake bounds.  This is used by fastDual
     279        so can exit gracefully before end */
     280    int numberAtFakeBound();
     281
     282    /** Pivot in a variable and choose an outgoing one.  Assumes dual
     283        feasible - will not go through a reduced cost.  Returns step length in theta
     284        Returns ray in ray_ (or NULL if no pivot)
     285        Return codes as before but -1 means no acceptable pivot
     286    */
     287    int pivotResult();
     288    /** Get next free , -1 if none */
     289    int nextSuperBasic();
     290    /** Startup part of dual (may be extended to other algorithms)
     291        returns 0 if good, 1 if bad */
     292    int startupSolve(int ifValuesPass, double * saveDuals, int startFinishOptions);
     293    void finishSolve(int startFinishOptions);
     294    void gutsOfDual(int ifValuesPass, double * & saveDuals, int initialStatus,
     295                    ClpDataSave & saveData);
     296    //int dual2(int ifValuesPass,int startFinishOptions=0);
     297    void resetFakeBounds(int type);
     298
     299    //@}
    300300};
    301301#endif
Note: See TracChangeset for help on using the changeset viewer.