Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpPEPrimalColumnSteepest.cpp

    r2150 r2385  
    2929// Default Constructor
    3030//-------------------------------------------------------------------
    31 ClpPEPrimalColumnSteepest::ClpPEPrimalColumnSteepest (double psi, int mode)
    32 : ClpPrimalColumnSteepest(mode), modelPE_(NULL), psi_(psi),
    33 iCurrent_(0), iInterval_(100),
    34   coDegenCompatibles_(0), coConsecutiveCompatibles_(0),
    35 updateCompatibles_(true)
     31ClpPEPrimalColumnSteepest::ClpPEPrimalColumnSteepest(double psi, int mode)
     32  : ClpPrimalColumnSteepest(mode)
     33  , modelPE_(NULL)
     34  , psi_(psi)
     35  , iCurrent_(0)
     36  , iInterval_(100)
     37  , coDegenCompatibles_(0)
     38  , coConsecutiveCompatibles_(0)
     39  , updateCompatibles_(true)
    3640{
    3741}
     
    4044// Copy constructor
    4145//-------------------------------------------------------------------
    42 ClpPEPrimalColumnSteepest::ClpPEPrimalColumnSteepest (const ClpPEPrimalColumnSteepest &source)
    43 : ClpPrimalColumnSteepest(source)
    44 {
    45         modelPE_    = NULL;
    46         psi_        = source.psi_;
    47         iCurrent_  = source.iCurrent_;
    48         iInterval_ = source.iInterval_;
    49         updateCompatibles_ = source.updateCompatibles_;
    50         coDegenCompatibles_ = source.coDegenCompatibles_;
    51         coConsecutiveCompatibles_ = source.coConsecutiveCompatibles_;
     46ClpPEPrimalColumnSteepest::ClpPEPrimalColumnSteepest(const ClpPEPrimalColumnSteepest &source)
     47  : ClpPrimalColumnSteepest(source)
     48{
     49  modelPE_ = NULL;
     50  psi_ = source.psi_;
     51  iCurrent_ = source.iCurrent_;
     52  iInterval_ = source.iInterval_;
     53  updateCompatibles_ = source.updateCompatibles_;
     54  coDegenCompatibles_ = source.coDegenCompatibles_;
     55  coConsecutiveCompatibles_ = source.coConsecutiveCompatibles_;
    5256}
    5357
     
    5559// Destructor
    5660//-------------------------------------------------------------------
    57 ClpPEPrimalColumnSteepest::~ClpPEPrimalColumnSteepest ()
     61ClpPEPrimalColumnSteepest::~ClpPEPrimalColumnSteepest()
    5862{
    5963  delete modelPE_;
     
    6468//-------------------------------------------------------------------
    6569ClpPEPrimalColumnSteepest &
    66 ClpPEPrimalColumnSteepest::operator=(const ClpPEPrimalColumnSteepest& rhs)
    67 {
    68         if (this != &rhs) {
    69                 ClpPrimalColumnSteepest::operator=(rhs);
    70                 delete modelPE_;
    71                 modelPE_=NULL;
    72         }
    73         return *this;
     70ClpPEPrimalColumnSteepest::operator=(const ClpPEPrimalColumnSteepest &rhs)
     71{
     72  if (this != &rhs) {
     73    ClpPrimalColumnSteepest::operator=(rhs);
     74    delete modelPE_;
     75    modelPE_ = NULL;
     76  }
     77  return *this;
    7478}
    7579
     
    7781// Clone
    7882//-------------------------------------------------------------------
    79 ClpPrimalColumnPivot * ClpPEPrimalColumnSteepest::clone(bool CopyData) const
    80 {
    81         if (CopyData) {
    82                 return new ClpPEPrimalColumnSteepest(*this);
    83         } else {
    84                 return new ClpPEPrimalColumnSteepest(psi_);
    85         }
    86 }
    87 
     83ClpPrimalColumnPivot *ClpPEPrimalColumnSteepest::clone(bool CopyData) const
     84{
     85  if (CopyData) {
     86    return new ClpPEPrimalColumnSteepest(*this);
     87  } else {
     88    return new ClpPEPrimalColumnSteepest(psi_);
     89  }
     90}
    8891
    8992// These have to match ClpPackedMatrix version
    9093#define TRY_NORM 1.0e-4
    9194#define ADD_ONE 1.0
    92 
    9395
    9496// Returns pivot column, -1 if none
     
    9698that is just +-weight where a feasibility changed.  It also has
    9799reduced cost from last iteration in pivot row*/
    98 int     ClpPEPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
    99         CoinIndexedVector * spareRow1,
    100         CoinIndexedVector * spareRow2,
    101         CoinIndexedVector * spareColumn1,
    102         CoinIndexedVector * spareColumn2)
    103 {
    104 // check that the model exists and has the proper form
    105         assert(model_);
    106         assert (!model_->nonLinearCost()->lookBothWays() && model_->algorithm() != 2);
    107 
    108 
    109 #if PE_DEBUG >=1
    110         std::cout << "@@@@@@@@@@ClpPEPrimalColumnSteepest::pivotColumn\n";
    111 #endif
    112 
    113 
    114         int number = 0;
    115         int * index;
    116 
    117         double tolerance = model_->currentDualTolerance();
    118 // we can't really trust infeasibilities if there is dual error
    119 // this coding has to mimic coding in checkDualSolution
    120         double error = CoinMin(1.0e-2, model_->largestDualError());
    121 // allow tolerance at least slightly bigger than standard
    122         tolerance = tolerance +  error;
    123         int pivotRow = model_->pivotRow();
    124         int anyUpdates;
    125         double * infeas = infeasible_->denseVector();
    126 
    127 // Local copy of mode so can decide what to do
    128         int switchType;
    129         if (mode_ == 4) switchType = 5 - numberSwitched_;
    130         else if (mode_ >= 10) switchType = 3;
    131         else switchType = mode_;
    132 
    133         /* switchType -
     100int ClpPEPrimalColumnSteepest::pivotColumn(CoinIndexedVector *updates,
     101  CoinIndexedVector *spareRow1,
     102  CoinIndexedVector *spareRow2,
     103  CoinIndexedVector *spareColumn1,
     104  CoinIndexedVector *spareColumn2)
     105{
     106  // check that the model exists and has the proper form
     107  assert(model_);
     108  assert(!model_->nonLinearCost()->lookBothWays() && model_->algorithm() != 2);
     109
     110#if PE_DEBUG >= 1
     111  std::cout << "@@@@@@@@@@ClpPEPrimalColumnSteepest::pivotColumn\n";
     112#endif
     113
     114  int number = 0;
     115  int *index;
     116
     117  double tolerance = model_->currentDualTolerance();
     118  // we can't really trust infeasibilities if there is dual error
     119  // this coding has to mimic coding in checkDualSolution
     120  double error = CoinMin(1.0e-2, model_->largestDualError());
     121  // allow tolerance at least slightly bigger than standard
     122  tolerance = tolerance + error;
     123  int pivotRow = model_->pivotRow();
     124  int anyUpdates;
     125  double *infeas = infeasible_->denseVector();
     126
     127  // Local copy of mode so can decide what to do
     128  int switchType;
     129  if (mode_ == 4)
     130    switchType = 5 - numberSwitched_;
     131  else if (mode_ >= 10)
     132    switchType = 3;
     133  else
     134    switchType = mode_;
     135
     136  /* switchType -
    134137        0 - all exact devex
    135138        1 - all steepest
     
    141144        */
    142145
    143         // Definition in ClpGubMatrix, mode is set to 4
    144         model_->clpMatrix()->dualExpanded(model_, updates, NULL, 4);
    145 
    146         if (updates->getNumElements() > 1) {
    147         // would have to have two goes for devex, three for steepest
    148                 anyUpdates = 2;
    149         }
    150         else if (updates->getNumElements()) {
    151                 if (updates->getIndices()[0] == pivotRow && fabs(updates->denseVector()[0]) > 1.0e-6) {
    152         // reasonable size
    153                         anyUpdates = 1;
    154                 }
    155                 else {
    156         // too small
    157                         anyUpdates = 2;
    158                 }
    159         }
    160         else if (pivotSequence_ >= 0) {
    161         // just after re-factorization
    162                 anyUpdates = -1;
    163         }
    164         else {
    165         // sub flip - nothing to do
    166                 anyUpdates = 0;
    167         }
    168         int sequenceOut = model_->sequenceOut();
    169 
    170         if (anyUpdates == 1) {
    171                 if (switchType < 4) {
    172         // exact etc when can use dj
    173                         djsAndSteepest(updates, spareRow2,
    174                                 spareColumn1, spareColumn2);
    175                 } else {
    176         // devex etc when can use dj
    177                         djsAndDevex(updates, spareRow2,
    178                                 spareColumn1, spareColumn2);
    179                 }
    180         } else if (anyUpdates == -1) {
    181                 if (switchType < 4) {
    182         // exact etc when djs okay
    183                         justSteepest(updates, spareRow2,
    184                                 spareColumn1, spareColumn2);
    185                 } else {
    186         // devex etc when djs okay
    187                         justDevex(updates, spareRow2,
    188                                 spareColumn1, spareColumn2);
    189                 }
    190         } else if (anyUpdates == 2) {
    191                 if (switchType < 4) {
    192         // exact etc when have to use pivot
    193                         djsAndSteepest2(updates, spareRow2,
    194                                 spareColumn1, spareColumn2);
    195                 } else {
    196         // devex etc when have to use pivot
    197                         djsAndDevex2(updates, spareRow2,
    198                                 spareColumn1, spareColumn2);
    199                 }
    200         }
    201 
    202         // make sure outgoing from last iteration okay
    203         if (sequenceOut >= 0) {
    204                 ClpSimplex::Status status = model_->getStatus(sequenceOut);
    205                 double value = model_->reducedCost(sequenceOut);
    206 
    207                 switch(status) {
    208 
    209                         case ClpSimplex::basic:
    210                         case ClpSimplex::isFixed:
    211                         break;
    212                         case ClpSimplex::isFree:
    213                         case ClpSimplex::superBasic:
    214                         if (fabs(value) > FREE_ACCEPT * tolerance) {
    215         // we are going to bias towards free (but only if reasonable)
    216                                 value *= FREE_BIAS;
    217         // store square in list
    218         if (infeas[sequenceOut]) infeas[sequenceOut] = value * value; // already there
    219         else infeasible_->quickAdd(sequenceOut, value * value);
    220         } else {
    221                 infeasible_->zero(sequenceOut);
    222         }
    223         break;
    224         case ClpSimplex::atUpperBound:
    225         if (value > tolerance) {
    226         // store square in list
    227         if (infeas[sequenceOut])        infeas[sequenceOut] = value * value; // already there
    228         else    infeasible_->quickAdd(sequenceOut, value * value);
    229         } else {
    230                 infeasible_->zero(sequenceOut);
    231         }
    232         break;
    233         case ClpSimplex::atLowerBound:
    234         if (value < -tolerance) {
    235         // store square in list
    236         if (infeas[sequenceOut])        infeas[sequenceOut] = value * value; // already there
    237         else    infeasible_->quickAdd(sequenceOut, value * value);
    238         } else {
    239                 infeasible_->zero(sequenceOut);
    240         }
    241         }
    242         }
    243 
    244         /* Determine whether the set of compatible variables should be updated
     146  // Definition in ClpGubMatrix, mode is set to 4
     147  model_->clpMatrix()->dualExpanded(model_, updates, NULL, 4);
     148
     149  if (updates->getNumElements() > 1) {
     150    // would have to have two goes for devex, three for steepest
     151    anyUpdates = 2;
     152  } else if (updates->getNumElements()) {
     153    if (updates->getIndices()[0] == pivotRow && fabs(updates->denseVector()[0]) > 1.0e-6) {
     154      // reasonable size
     155      anyUpdates = 1;
     156    } else {
     157      // too small
     158      anyUpdates = 2;
     159    }
     160  } else if (pivotSequence_ >= 0) {
     161    // just after re-factorization
     162    anyUpdates = -1;
     163  } else {
     164    // sub flip - nothing to do
     165    anyUpdates = 0;
     166  }
     167  int sequenceOut = model_->sequenceOut();
     168
     169  if (anyUpdates == 1) {
     170    if (switchType < 4) {
     171      // exact etc when can use dj
     172      djsAndSteepest(updates, spareRow2,
     173        spareColumn1, spareColumn2);
     174    } else {
     175      // devex etc when can use dj
     176      djsAndDevex(updates, spareRow2,
     177        spareColumn1, spareColumn2);
     178    }
     179  } else if (anyUpdates == -1) {
     180    if (switchType < 4) {
     181      // exact etc when djs okay
     182      justSteepest(updates, spareRow2,
     183        spareColumn1, spareColumn2);
     184    } else {
     185      // devex etc when djs okay
     186      justDevex(updates, spareRow2,
     187        spareColumn1, spareColumn2);
     188    }
     189  } else if (anyUpdates == 2) {
     190    if (switchType < 4) {
     191      // exact etc when have to use pivot
     192      djsAndSteepest2(updates, spareRow2,
     193        spareColumn1, spareColumn2);
     194    } else {
     195      // devex etc when have to use pivot
     196      djsAndDevex2(updates, spareRow2,
     197        spareColumn1, spareColumn2);
     198    }
     199  }
     200
     201  // make sure outgoing from last iteration okay
     202  if (sequenceOut >= 0) {
     203    ClpSimplex::Status status = model_->getStatus(sequenceOut);
     204    double value = model_->reducedCost(sequenceOut);
     205
     206    switch (status) {
     207
     208    case ClpSimplex::basic:
     209    case ClpSimplex::isFixed:
     210      break;
     211    case ClpSimplex::isFree:
     212    case ClpSimplex::superBasic:
     213      if (fabs(value) > FREE_ACCEPT * tolerance) {
     214        // we are going to bias towards free (but only if reasonable)
     215        value *= FREE_BIAS;
     216        // store square in list
     217        if (infeas[sequenceOut])
     218          infeas[sequenceOut] = value * value; // already there
     219        else
     220          infeasible_->quickAdd(sequenceOut, value * value);
     221      } else {
     222        infeasible_->zero(sequenceOut);
     223      }
     224      break;
     225    case ClpSimplex::atUpperBound:
     226      if (value > tolerance) {
     227        // store square in list
     228        if (infeas[sequenceOut])
     229          infeas[sequenceOut] = value * value; // already there
     230        else
     231          infeasible_->quickAdd(sequenceOut, value * value);
     232      } else {
     233        infeasible_->zero(sequenceOut);
     234      }
     235      break;
     236    case ClpSimplex::atLowerBound:
     237      if (value < -tolerance) {
     238        // store square in list
     239        if (infeas[sequenceOut])
     240          infeas[sequenceOut] = value * value; // already there
     241        else
     242          infeasible_->quickAdd(sequenceOut, value * value);
     243      } else {
     244        infeasible_->zero(sequenceOut);
     245      }
     246    }
     247  }
     248
     249  /* Determine whether the set of compatible variables should be updated
    245250        */
    246         // store the number of degenerate pivots on compatible variables and the
    247         // overall number of degenerate pivots
    248         //double progress = fabs(modelPE_->lastObjectiveValue() - model_->objectiveValue());
    249         //bool isLastDegenerate = progress <= 1.0e-12*fabs(model_->objectiveValue()) ? true:false;
    250         bool isLastDegenerate = fabs(model_->theta()) < 1.01e-7;
    251         bool isLastDegenerate2;
    252         // could fine tune a bit e.g. use 0.0 rather than tolerance
    253         if (model_->directionOut()<0) {
    254           // going out to lower bound
    255           isLastDegenerate2 = (model_->valueOut()-model_->lowerOut()<model_->primalTolerance());
    256         } else {
    257           // going out to upper bound
    258           isLastDegenerate2 = (model_->upperOut()-model_->valueOut()<model_->primalTolerance());
    259         }
     251  // store the number of degenerate pivots on compatible variables and the
     252  // overall number of degenerate pivots
     253  //double progress = fabs(modelPE_->lastObjectiveValue() - model_->objectiveValue());
     254  //bool isLastDegenerate = progress <= 1.0e-12*fabs(model_->objectiveValue()) ? true:false;
     255  bool isLastDegenerate = fabs(model_->theta()) < 1.01e-7;
     256  bool isLastDegenerate2;
     257  // could fine tune a bit e.g. use 0.0 rather than tolerance
     258  if (model_->directionOut() < 0) {
     259    // going out to lower bound
     260    isLastDegenerate2 = (model_->valueOut() - model_->lowerOut() < model_->primalTolerance());
     261  } else {
     262    // going out to upper bound
     263    isLastDegenerate2 = (model_->upperOut() - model_->valueOut() < model_->primalTolerance());
     264  }
    260265#if 0
    261266        if ( isLastDegenerate2 &&!isLastDegenerate)
     
    268273                 (model_->directionOut()<0) ? "lower" : "upper");
    269274#else
    270         isLastDegenerate=isLastDegenerate2;
     275  isLastDegenerate = isLastDegenerate2;
    271276#endif
    272         if ( isLastDegenerate ) {
    273                 modelPE_->addDegeneratePivot();
    274                 modelPE_->addDegeneratePivotConsecutive();
    275 
    276                 if ( modelPE_->isLastPivotCompatible() ) {
    277                         modelPE_->addDegenerateCompatiblePivot();
    278                 }
    279         }
    280         else {
    281                 modelPE_->resetDegeneratePivotsConsecutive();
    282         }
    283 
    284         // the compatible variables are updated when the number of degenerate pivots
    285         // on compatible variables is more than 20% the overall number of degenerate pivots
    286         if ( modelPE_->isLastPivotCompatible() ) {
    287                 coConsecutiveCompatibles_++;
    288                 if (isLastDegenerate) {
    289                         coDegenCompatibles_++;
    290                         if ( coConsecutiveCompatibles_ >= 10 && 5*coDegenCompatibles_*model_->numberIterations() > modelPE_->coDegeneratePivots()*coConsecutiveCompatibles_ )  {
    291                                 updateCompatibles_ = true;
    292                         }
    293                 }
    294         }
    295 
    296         if (modelPE_->doStatistics()) {
    297         /* For results comparison.
     277  if (isLastDegenerate) {
     278    modelPE_->addDegeneratePivot();
     279    modelPE_->addDegeneratePivotConsecutive();
     280
     281    if (modelPE_->isLastPivotCompatible()) {
     282      modelPE_->addDegenerateCompatiblePivot();
     283    }
     284  } else {
     285    modelPE_->resetDegeneratePivotsConsecutive();
     286  }
     287
     288  // the compatible variables are updated when the number of degenerate pivots
     289  // on compatible variables is more than 20% the overall number of degenerate pivots
     290  if (modelPE_->isLastPivotCompatible()) {
     291    coConsecutiveCompatibles_++;
     292    if (isLastDegenerate) {
     293      coDegenCompatibles_++;
     294      if (coConsecutiveCompatibles_ >= 10 && 5 * coDegenCompatibles_ * model_->numberIterations() > modelPE_->coDegeneratePivots() * coConsecutiveCompatibles_) {
     295        updateCompatibles_ = true;
     296      }
     297    }
     298  }
     299
     300  if (modelPE_->doStatistics()) {
     301    /* For results comparison.
    298302        count the number of degenerate variables if psi==1
    299303        add the time spent in counting to the time limit
    300304        */
    301         modelPE_->startTimer();
    302         if (psi_ >= 1 && iCurrent_ >= 100) {
    303                 modelPE_->updateDualDegenerates();
    304                 modelPE_->updateDualDegeneratesAvg(100);
    305                 //model_->setMaximumSeconds(36000.0+modelPE_->timeCompatibility()-CoinCpuTime());
    306                 iCurrent_ = 0;
    307         }
    308         modelPE_->stopTimer();
    309         }
    310 
    311         /* Update the set of compatible variables if necessary
     305    modelPE_->startTimer();
     306    if (psi_ >= 1 && iCurrent_ >= 100) {
     307      modelPE_->updateDualDegenerates();
     308      modelPE_->updateDualDegeneratesAvg(100);
     309      //model_->setMaximumSeconds(36000.0+modelPE_->timeCompatibility()-CoinCpuTime());
     310      iCurrent_ = 0;
     311    }
     312    modelPE_->stopTimer();
     313  }
     314
     315  /* Update the set of compatible variables if necessary
    312316        */
    313         if (modelPE_->doStatistics())
    314         modelPE_->startTimer();   
    315         double psiTmp = psi_;
    316         if ( (psi_ < 1.0) && (iCurrent_ >= iInterval_) && (updateCompatibles_||iCurrent_ >= 1000) )
    317         {
    318                 // the compatible variables are never updated if the last pivot is non degenerate
    319                 // this could be counterproductive
    320                 if (isLastDegenerate ) {
    321                         modelPE_->updatePrimalDegenerates();
    322                         for (int i=0;i<4;i++)
    323                           assert(!model_->rowArray(i)->getNumElements());
    324                         modelPE_->identifyCompatibleCols(model_->numberRows()+model_->numberColumns(), NULL, spareRow2,spareRow1);
    325                         for (int i=0;i<4;i++)
    326                           assert(!model_->rowArray(i)->getNumElements());
    327 
    328         if (modelPE_->doStatistics()) {
    329         // update the average number of degenerates and compatibles for statistics
    330                         modelPE_->updatePrimalDegeneratesAvg(iCurrent_);
    331                         modelPE_->updateCompatibleColsAvg(iCurrent_);
    332                         if (modelPE_->doStatistics()>3) {
    333                           char generalPrint[100];
    334                          
    335                           sprintf(generalPrint,"coDegen = %d; coComp = %d; iCurrent_ = %d; compatibleColumns = %d",
    336                                   coDegenCompatibles_,coConsecutiveCompatibles_,
    337                                   iCurrent_,modelPE_->coCompatibleCols());
    338                           model_->messageHandler()->message(CLP_GENERAL,
    339                                                             *model_->messagesPointer())
    340                             << generalPrint << CoinMessageEol;
    341                         }
    342         }
    343         // switch off ? - maybe just until drops back
    344         if (modelPE_->coCompatibleCols()*1.0>model_->numberColumns()) {
    345           printf("switching off pe\n");
    346           psi_ = 1.0;
    347         }
    348 
    349                         // the checking frequency is modified to reflect what appears to be needed
    350                         if (iCurrent_ == iInterval_) iInterval_ = std::max(50, iInterval_-50);
    351                         else                                                              iInterval_ = std::min(300, iInterval_+50);
    352 
    353 
    354                         // reset all the indicators that are used to decide whether the compatible
    355                         // variables must be updated
    356                         iCurrent_ = 0;
    357                         updateCompatibles_ = false;
    358                         coConsecutiveCompatibles_ = 0;
    359                         coDegenCompatibles_ = 0;
    360                 }
    361                 else iInterval_++;
    362         }
    363                 /* otherwise, update the value of the priority factor depending on the number of
     317  if (modelPE_->doStatistics())
     318    modelPE_->startTimer();
     319  double psiTmp = psi_;
     320  if ((psi_ < 1.0) && (iCurrent_ >= iInterval_) && (updateCompatibles_ || iCurrent_ >= 1000)) {
     321    // the compatible variables are never updated if the last pivot is non degenerate
     322    // this could be counterproductive
     323    if (isLastDegenerate) {
     324      modelPE_->updatePrimalDegenerates();
     325      for (int i = 0; i < 4; i++)
     326        assert(!model_->rowArray(i)->getNumElements());
     327      modelPE_->identifyCompatibleCols(model_->numberRows() + model_->numberColumns(), NULL, spareRow2, spareRow1);
     328      for (int i = 0; i < 4; i++)
     329        assert(!model_->rowArray(i)->getNumElements());
     330
     331      if (modelPE_->doStatistics()) {
     332        // update the average number of degenerates and compatibles for statistics
     333        modelPE_->updatePrimalDegeneratesAvg(iCurrent_);
     334        modelPE_->updateCompatibleColsAvg(iCurrent_);
     335        if (modelPE_->doStatistics() > 3) {
     336          char generalPrint[100];
     337
     338          sprintf(generalPrint, "coDegen = %d; coComp = %d; iCurrent_ = %d; compatibleColumns = %d",
     339            coDegenCompatibles_, coConsecutiveCompatibles_,
     340            iCurrent_, modelPE_->coCompatibleCols());
     341          model_->messageHandler()->message(CLP_GENERAL,
     342            *model_->messagesPointer())
     343            << generalPrint << CoinMessageEol;
     344        }
     345      }
     346      // switch off ? - maybe just until drops back
     347      if (modelPE_->coCompatibleCols() * 1.0 > model_->numberColumns()) {
     348        printf("switching off pe\n");
     349        psi_ = 1.0;
     350      }
     351
     352      // the checking frequency is modified to reflect what appears to be needed
     353      if (iCurrent_ == iInterval_)
     354        iInterval_ = std::max(50, iInterval_ - 50);
     355      else
     356        iInterval_ = std::min(300, iInterval_ + 50);
     357
     358      // reset all the indicators that are used to decide whether the compatible
     359      // variables must be updated
     360      iCurrent_ = 0;
     361      updateCompatibles_ = false;
     362      coConsecutiveCompatibles_ = 0;
     363      coDegenCompatibles_ = 0;
     364    } else
     365      iInterval_++;
     366  }
     367  /* otherwise, update the value of the priority factor depending on the number of
    364368         * consecutive degenerate pivots
    365369         */
    366         else {
    367                 // the idea is that when a lot of consecutive pivots are degenerate, there is
    368                 // a potentially high added value in putting a very high priority on compatible
    369                 // variables
    370                 if (modelPE_->coDegeneratePivotsConsecutive() >= 10) {
    371                         psiTmp = 0.01;
    372                         psiTmp = 0.25*psi_;
    373 
    374                         #if PE_DEBUG >= 1
    375                         std::cout << "Lower the priority factor to " << psiTmp << std::endl;
    376                         std::cout << "Consecutive degenerate pivots=" << modelPE_->coDegeneratePivotsConsecutive() << std::endl;
    377                         #endif
    378                 }
    379         }
    380         iCurrent_++;
    381         if (modelPE_->doStatistics())
    382         modelPE_->stopTimer();
    383 
    384 
    385         /* Update of the dual solution and compatible variables is finished,
     370  else {
     371    // the idea is that when a lot of consecutive pivots are degenerate, there is
     372    // a potentially high added value in putting a very high priority on compatible
     373    // variables
     374    if (modelPE_->coDegeneratePivotsConsecutive() >= 10) {
     375      psiTmp = 0.01;
     376      psiTmp = 0.25 * psi_;
     377
     378#if PE_DEBUG >= 1
     379      std::cout << "Lower the priority factor to " << psiTmp << std::endl;
     380      std::cout << "Consecutive degenerate pivots=" << modelPE_->coDegeneratePivotsConsecutive() << std::endl;
     381#endif
     382    }
     383  }
     384  iCurrent_++;
     385  if (modelPE_->doStatistics())
     386    modelPE_->stopTimer();
     387
     388  /* Update of the dual solution and compatible variables is finished,
    386389        ** now do the pricing */
    387390
    388         // See what sort of pricing
    389         int numberWanted = 10;
    390         number = infeasible_->getNumElements();
    391         int numberColumns = model_->numberColumns();
    392         int numberRows = model_->numberRows();
    393         // ratio is done on number of rows here
    394         double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
    395         if(switchType == 4) {
    396         // Still in devex mode
    397         // Go to steepest if lot of iterations?
    398                 if (ratio < 5.0) {
    399                         numberWanted = CoinMax(2000, number / 10);
    400                         numberWanted = CoinMax(numberWanted, numberColumns / 20);
    401                 } else if (ratio < 7.0) {
    402                         numberWanted = CoinMax(2000, number / 5);
    403                         numberWanted = CoinMax(numberWanted, numberColumns / 10);
    404                 } else {
    405         // we can zero out
    406                         updates->clear();
    407                         spareColumn1->clear();
    408                         switchType = 3;
    409         // initialize
    410                         pivotSequence_ = -1;
    411                         pivotRow = -1;
    412                         numberSwitched_++;
    413         // Make sure will re-do
    414                         delete [] weights_;
    415                         weights_ = NULL;
    416                         saveWeights(model_, 4);
    417                         COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
    418                         updates->clear();
    419                 }
    420                 if (model_->numberIterations() % 1000 == 0)
    421                         COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
    422         }
    423         if (switchType < 4) {
    424                 if (switchType < 2 ) {
    425                         numberWanted = number + 1;
    426                 } else if (switchType == 2) {
    427                         numberWanted = CoinMax(2000, number / 8);
    428                 } else {
    429                         if (ratio < 1.0) {
    430                                 numberWanted = CoinMax(2000, number / 20);
    431                         } else if (ratio < 5.0) {
    432                                 numberWanted = CoinMax(2000, number / 10);
    433                                 numberWanted = CoinMax(numberWanted, numberColumns / 40);
    434                         } else if (ratio < 10.0) {
    435                                 numberWanted = CoinMax(2000, number / 8);
    436                                 numberWanted = CoinMax(numberWanted, numberColumns / 20);
    437                         } else {
    438                                 ratio = number * (ratio / 80.0);
    439                                 if (ratio > number) {
    440                                         numberWanted = number + 1;
    441                                 } else {
    442                                         numberWanted = CoinMax(2000, static_cast<int> (ratio));
    443                                         numberWanted = CoinMax(numberWanted, numberColumns / 10);
    444                                 }
    445                         }
    446                 }
    447         }
    448 
    449         // initialize the best reduced cost values
    450         double bestDj = 1.0e-30;
    451         int bestSequence = -1;
    452         double bestDjComp = 1.0e-30;
    453         int bestSequenceComp = -1;
    454 
    455         int i, iSequence;
    456         index = infeasible_->getIndices();
    457         number = infeasible_->getNumElements();
    458         // Re-sort infeasible every 100 pivots
    459         // Not a good idea
    460         if (0 && model_->factorization()->pivots() > 0 &&
    461                 (model_->factorization()->pivots() % 100) == 0) {
    462                 int nLook = model_->numberRows() + numberColumns;
    463         number = 0;
    464         for (i = 0; i < nLook; i++) {
    465                 if (infeas[i]) {
    466                         if (fabs(infeas[i]) > COIN_INDEXED_TINY_ELEMENT)
    467                                 index[number++] = i;
    468                         else
    469                                 infeas[i] = 0.0;
    470                 }
    471         }
    472         infeasible_->setNumElements(number);
    473         }
    474         if(model_->numberIterations() < model_->lastBadIteration() + 200 &&
    475                 model_->factorization()->pivots() > 10) {
    476         // we can't really trust infeasibilities if there is dual error
    477                 double checkTolerance = 1.0e-8;
    478         if (model_->largestDualError() > checkTolerance)
    479                 tolerance *= model_->largestDualError() / checkTolerance;
    480         // But cap
    481         tolerance = CoinMin(1000.0, tolerance);
    482         }
    483 
    484         // stop last one coming immediately
    485         double saveOutInfeasibility = 0.0;
    486         if (sequenceOut >= 0) {
    487                 saveOutInfeasibility = infeas[sequenceOut];
    488                 infeas[sequenceOut] = 0.0;
    489         }
    490         if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
    491                 tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
    492         tolerance *= tolerance; // as we are using squares
    493 
    494         // only check the compatible variables when the bidimensional factor is less than 1
    495         // and the ratio of compatible variables is larger than 0.01
    496         bool checkCompatibles = true;
    497         double ratioCompatibles = static_cast<double> ( modelPE_->coCompatibleCols())/
    498           static_cast<double> ((model_->numberRows()+model_->numberColumns()));
    499         double ratioCompatibles2 = static_cast<double> ( modelPE_->coCompatibleCols())/
    500           static_cast<double> (model_->numberColumns());
    501 
    502         if (psi_ >= 1.0 || ratioCompatibles < 0.01)     checkCompatibles = false;
    503         if (ratioCompatibles2>0.5) checkCompatibles=false;
    504         // proceed to a partial pricing in two phases checking the weighted reduced
    505         // costs of the variables until numberWanted variables have been checked,
    506         // the number of variables explored in the first phase is randomly generated
    507         int iPass;
    508         int start[4];
    509         start[1] = number;
    510         start[2] = 0;
    511         double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
    512         start[0] = static_cast<int> (dstart);
    513         start[3] = start[0];
    514         for (iPass = 0; iPass < 2; iPass++) {
    515                 int end = start[2*iPass+1];
    516 
    517                 for (i = start[2*iPass]; i < end; i++) {
    518                         iSequence = index[i];
    519                         double value = infeas[iSequence];
    520                         double weight = weights_[iSequence];
    521                         double weightedDj = weight * bestDj;
    522                         double largestWeightedDj =  std::max(psi_*weightedDj, weight * bestDjComp);
    523                         if (value > tolerance) {
    524                                 if (value > largestWeightedDj ) {
    525                                         if (model_->flagged(iSequence)) {
    526         // just to make sure we don't exit before got something
    527                                                 numberWanted++;
    528                                         }
    529                                         else if (checkCompatibles && modelPE_->isCompatibleCol(iSequence) ) {
    530         // the condition on value is sufficient if the variable is compatible
    531                                                 bestDjComp = value / weight;
    532                                                 bestSequenceComp = iSequence;
    533                                         }
    534                                         else if (value >  weightedDj) {
    535         // the usual condition is checked if the variable is not compatible
    536                                                 bestDj = value / weight;
    537                                                 bestSequence = iSequence;
    538                                         }
    539                                 }
    540                                 numberWanted--;
    541                         }
    542                         if (!numberWanted)
    543                                 break;
    544                 }
    545                 if (!numberWanted)
    546                         break;
    547         }
    548         if (bestSequence>=0&&model_->getStatus(bestSequence)==
    549             ClpSimplex::isFree) {
    550           printf("Free in %d compat %c dj %g\n",
    551                  bestSequence,modelPE_->isCompatibleCol(bestSequence) ? 'y' : 'n',bestDj);
    552           bestDjComp=0.0;
    553         }
    554         // choose a compatible or an incompatible row depending on the
    555         // largest values and on the factor of preference psi_
    556         if ( bestSequenceComp >= 0 && bestDjComp >= psiTmp * bestDj)  {
     391  // See what sort of pricing
     392  int numberWanted = 10;
     393  number = infeasible_->getNumElements();
     394  int numberColumns = model_->numberColumns();
     395  int numberRows = model_->numberRows();
     396  // ratio is done on number of rows here
     397  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
     398  if (switchType == 4) {
     399    // Still in devex mode
     400    // Go to steepest if lot of iterations?
     401    if (ratio < 5.0) {
     402      numberWanted = CoinMax(2000, number / 10);
     403      numberWanted = CoinMax(numberWanted, numberColumns / 20);
     404    } else if (ratio < 7.0) {
     405      numberWanted = CoinMax(2000, number / 5);
     406      numberWanted = CoinMax(numberWanted, numberColumns / 10);
     407    } else {
     408      // we can zero out
     409      updates->clear();
     410      spareColumn1->clear();
     411      switchType = 3;
     412      // initialize
     413      pivotSequence_ = -1;
     414      pivotRow = -1;
     415      numberSwitched_++;
     416      // Make sure will re-do
     417      delete[] weights_;
     418      weights_ = NULL;
     419      saveWeights(model_, 4);
     420      COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
     421      updates->clear();
     422    }
     423    if (model_->numberIterations() % 1000 == 0)
     424      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
     425  }
     426  if (switchType < 4) {
     427    if (switchType < 2) {
     428      numberWanted = number + 1;
     429    } else if (switchType == 2) {
     430      numberWanted = CoinMax(2000, number / 8);
     431    } else {
     432      if (ratio < 1.0) {
     433        numberWanted = CoinMax(2000, number / 20);
     434      } else if (ratio < 5.0) {
     435        numberWanted = CoinMax(2000, number / 10);
     436        numberWanted = CoinMax(numberWanted, numberColumns / 40);
     437      } else if (ratio < 10.0) {
     438        numberWanted = CoinMax(2000, number / 8);
     439        numberWanted = CoinMax(numberWanted, numberColumns / 20);
     440      } else {
     441        ratio = number * (ratio / 80.0);
     442        if (ratio > number) {
     443          numberWanted = number + 1;
     444        } else {
     445          numberWanted = CoinMax(2000, static_cast< int >(ratio));
     446          numberWanted = CoinMax(numberWanted, numberColumns / 10);
     447        }
     448      }
     449    }
     450  }
     451
     452  // initialize the best reduced cost values
     453  double bestDj = 1.0e-30;
     454  int bestSequence = -1;
     455  double bestDjComp = 1.0e-30;
     456  int bestSequenceComp = -1;
     457
     458  int i, iSequence;
     459  index = infeasible_->getIndices();
     460  number = infeasible_->getNumElements();
     461  // Re-sort infeasible every 100 pivots
     462  // Not a good idea
     463  if (0 && model_->factorization()->pivots() > 0 && (model_->factorization()->pivots() % 100) == 0) {
     464    int nLook = model_->numberRows() + numberColumns;
     465    number = 0;
     466    for (i = 0; i < nLook; i++) {
     467      if (infeas[i]) {
     468        if (fabs(infeas[i]) > COIN_INDEXED_TINY_ELEMENT)
     469          index[number++] = i;
     470        else
     471          infeas[i] = 0.0;
     472      }
     473    }
     474    infeasible_->setNumElements(number);
     475  }
     476  if (model_->numberIterations() < model_->lastBadIteration() + 200 && model_->factorization()->pivots() > 10) {
     477    // we can't really trust infeasibilities if there is dual error
     478    double checkTolerance = 1.0e-8;
     479    if (model_->largestDualError() > checkTolerance)
     480      tolerance *= model_->largestDualError() / checkTolerance;
     481    // But cap
     482    tolerance = CoinMin(1000.0, tolerance);
     483  }
     484
     485  // stop last one coming immediately
     486  double saveOutInfeasibility = 0.0;
     487  if (sequenceOut >= 0) {
     488    saveOutInfeasibility = infeas[sequenceOut];
     489    infeas[sequenceOut] = 0.0;
     490  }
     491  if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
     492    tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
     493  tolerance *= tolerance; // as we are using squares
     494
     495  // only check the compatible variables when the bidimensional factor is less than 1
     496  // and the ratio of compatible variables is larger than 0.01
     497  bool checkCompatibles = true;
     498  double ratioCompatibles = static_cast< double >(modelPE_->coCompatibleCols()) / static_cast< double >((model_->numberRows() + model_->numberColumns()));
     499  double ratioCompatibles2 = static_cast< double >(modelPE_->coCompatibleCols()) / static_cast< double >(model_->numberColumns());
     500
     501  if (psi_ >= 1.0 || ratioCompatibles < 0.01)
     502    checkCompatibles = false;
     503  if (ratioCompatibles2 > 0.5)
     504    checkCompatibles = false;
     505  // proceed to a partial pricing in two phases checking the weighted reduced
     506  // costs of the variables until numberWanted variables have been checked,
     507  // the number of variables explored in the first phase is randomly generated
     508  int iPass;
     509  int start[4];
     510  start[1] = number;
     511  start[2] = 0;
     512  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
     513  start[0] = static_cast< int >(dstart);
     514  start[3] = start[0];
     515  for (iPass = 0; iPass < 2; iPass++) {
     516    int end = start[2 * iPass + 1];
     517
     518    for (i = start[2 * iPass]; i < end; i++) {
     519      iSequence = index[i];
     520      double value = infeas[iSequence];
     521      double weight = weights_[iSequence];
     522      double weightedDj = weight * bestDj;
     523      double largestWeightedDj = std::max(psi_ * weightedDj, weight * bestDjComp);
     524      if (value > tolerance) {
     525        if (value > largestWeightedDj) {
     526          if (model_->flagged(iSequence)) {
     527            // just to make sure we don't exit before got something
     528            numberWanted++;
     529          } else if (checkCompatibles && modelPE_->isCompatibleCol(iSequence)) {
     530            // the condition on value is sufficient if the variable is compatible
     531            bestDjComp = value / weight;
     532            bestSequenceComp = iSequence;
     533          } else if (value > weightedDj) {
     534            // the usual condition is checked if the variable is not compatible
     535            bestDj = value / weight;
     536            bestSequence = iSequence;
     537          }
     538        }
     539        numberWanted--;
     540      }
     541      if (!numberWanted)
     542        break;
     543    }
     544    if (!numberWanted)
     545      break;
     546  }
     547  if (bestSequence >= 0 && model_->getStatus(bestSequence) == ClpSimplex::isFree) {
     548    printf("Free in %d compat %c dj %g\n",
     549      bestSequence, modelPE_->isCompatibleCol(bestSequence) ? 'y' : 'n', bestDj);
     550    bestDjComp = 0.0;
     551  }
     552  // choose a compatible or an incompatible row depending on the
     553  // largest values and on the factor of preference psi_
     554  if (bestSequenceComp >= 0 && bestDjComp >= psiTmp * bestDj) {
    557555
    558556#if 0
     
    564562                 bestSequenceComp,bestDjComp,bestDj);
    565563#endif
    566 // record the number of pivots done on compatible variables
    567 // that would not have been done without positive edge
    568         if (modelPE_->doStatistics()) {
    569                 if (bestDjComp < bestDj) modelPE_->addPriorityPivot();
    570         }
    571                 bestSequence = bestSequenceComp;
    572         }
    573         if ( bestSequence >= 0 && psi_ < 1.0 && modelPE_->isCompatibleCol(bestSequence) )  {
    574                 modelPE_->isLastPivotCompatible(true);
    575                 modelPE_->addCompatiblePivot();
    576         }
    577         else
    578                 modelPE_->isLastPivotCompatible(false);
    579 
    580         model_->clpMatrix()->setSavedBestSequence(bestSequence);
    581         if (bestSequence >= 0)
    582                 model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
    583         if (sequenceOut >= 0) {
    584                 infeas[sequenceOut] = saveOutInfeasibility;
    585         }
    586 
    587         modelPE_->updateLastObjectiveValue();
    588         return bestSequence;
     564    // record the number of pivots done on compatible variables
     565    // that would not have been done without positive edge
     566    if (modelPE_->doStatistics()) {
     567      if (bestDjComp < bestDj)
     568        modelPE_->addPriorityPivot();
     569    }
     570    bestSequence = bestSequenceComp;
     571  }
     572  if (bestSequence >= 0 && psi_ < 1.0 && modelPE_->isCompatibleCol(bestSequence)) {
     573    modelPE_->isLastPivotCompatible(true);
     574    modelPE_->addCompatiblePivot();
     575  } else
     576    modelPE_->isLastPivotCompatible(false);
     577
     578  model_->clpMatrix()->setSavedBestSequence(bestSequence);
     579  if (bestSequence >= 0)
     580    model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
     581  if (sequenceOut >= 0) {
     582    infeas[sequenceOut] = saveOutInfeasibility;
     583  }
     584
     585  modelPE_->updateLastObjectiveValue();
     586  return bestSequence;
    589587}
    590588/* Save weights - this may initialize weights as well
    591589   This is as parent but may initialize ClpPESimplex
    592590*/
    593 void
    594 ClpPEPrimalColumnSteepest::saveWeights(ClpSimplex * model, int mode)
     591void ClpPEPrimalColumnSteepest::saveWeights(ClpSimplex *model, int mode)
    595592{
    596593  // See if we need to initialize ClpPESimplex
    597   if (!modelPE_||model!=modelPE_->clpModel()||
    598       !modelPE_->checkSize()) {
     594  if (!modelPE_ || model != modelPE_->clpModel() || !modelPE_->checkSize()) {
    599595    delete modelPE_;
    600596    modelPE_ = new ClpPESimplex(model);
    601597  }
    602   ClpPrimalColumnSteepest::saveWeights(model,mode);
     598  ClpPrimalColumnSteepest::saveWeights(model, mode);
    603599}
    604600// Updates weights - as ordinary but checks for zero moves
    605 void
    606 ClpPEPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
     601void ClpPEPrimalColumnSteepest::updateWeights(CoinIndexedVector *input)
    607602{
    608603  //if (modelPE_->isLastPivotCompatible()) {
     
    613608}
    614609
    615 
     610/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     611*/
Note: See TracChangeset for help on using the changeset viewer.