Changeset 54


Ignore:
Timestamp:
Nov 5, 2002 2:01:40 PM (17 years ago)
Author:
forrest
Message:

First attempt at perturbation

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpNonLinearCost.cpp

    r50 r54  
    3333  largestInfeasibility_(0.0),
    3434  sumInfeasibilities_(0.0),
    35   convex_(true)
     35  convex_(true),
     36  infeasible_(NULL)
    3637{
    3738
     
    7576  lower_ = new double [put];
    7677  cost_ = new double [put];
     78  infeasible_ = new unsigned int[(put+31)>>5];
     79  memset(infeasible_,0,((put+31)>>5)*sizeof(unsigned int));
    7780
    7881  put=0;
     
    8386
    8487    lower_[put] = -DBL_MAX;
     88    setInfeasible(put,true);
     89
    8590    cost_[put++] = cost[iSequence]-infeasibilityCost;
    8691    whichRange_[iSequence]=put;
     
    9196    if (upper[iSequence]<1.0e20) {
    9297      lower_[put] = DBL_MAX;
     98      setInfeasible(put-1,true);
    9399      cost_[put++] = 1.0e50;
    94100    }
     
    138144  lower_ = new double [put];
    139145  cost_ = new double [put];
     146  infeasible_ = new unsigned int[(put+31)>>5];
     147  memset(infeasible_,0,((put+31)>>5)*sizeof(unsigned int));
    140148
    141149  // now fill in
     
    146154  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    147155    lower_[put] = -DBL_MAX;
     156    setInfeasible(put,true);
    148157    cost_[put++] = cost[iSequence]-infeasibilityCost;
    149158    whichRange_[iSequence]=put;
     
    177186    if (upper[iSequence]<1.0e20) {
    178187      lower_[put] = DBL_MAX;
     188      setInfeasible(put-1,true);
    179189      cost_[put++] = 1.0e50;
    180190    }
     
    201211  largestInfeasibility_(0.0),
    202212  sumInfeasibilities_(0.0),
    203   convex_(true)
     213  convex_(true),
     214  infeasible_(NULL)
    204215
    205216  if (numberRows_) {
     
    222233    sumInfeasibilities_ = rhs.sumInfeasibilities_;
    223234    convex_ = rhs.convex_;
     235    infeasible_ = new unsigned int[(numberEntries+31)>>5];
     236    memcpy(infeasible_,rhs.infeasible_,
     237           ((numberEntries+31)>>5)*sizeof(unsigned int));
    224238  }
    225239}
     
    235249  delete [] lower_;
    236250  delete [] cost_;
    237 
     251  delete [] infeasible_;
    238252}
    239253
     
    252266    delete [] lower_;
    253267    delete []cost_;
     268    delete [] infeasible_;
    254269    start_ = NULL;
    255270    whichRange_ = NULL;
    256271    lower_ = NULL;
    257272    cost_= NULL;
     273    infeasible_=NULL;
    258274    if (numberRows_) {
    259275      int numberTotal = numberRows_+numberColumns_;
     
    269285      cost_ = new double [numberEntries];
    270286      memcpy(cost_,rhs.cost_,numberEntries*sizeof(double));
     287      infeasible_ = new unsigned int[(numberEntries+31)>>5];
     288      memcpy(infeasible_,rhs.infeasible_,
     289             ((numberEntries+31)>>5)*sizeof(unsigned int));
    271290    }
    272291    model_ = rhs.model_;
     
    302321  // nonbasic should be at a valid bound
    303322  for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++) {
    304     double lowerValue=lower[iSequence];
    305     double upperValue=upper[iSequence];
     323    double lowerValue;
     324    double upperValue;
    306325    double value=solution[iSequence];
    307326    int iRange;
    308     // correct costs
    309    
    310     if (lowerValue>-1.0e20) {
    311       iRange=start_[iSequence];
    312       cost_[iRange] = cost[iSequence]-infeasibilityCost;
    313     }
    314     if (upperValue<1.0e20) {
    315       iRange=start_[iSequence+1]-2;
    316       cost_[iRange] = cost[iSequence]+infeasibilityCost;
    317     }
    318327    // get correct place
    319328    int start = start_[iSequence];
    320329    int end = start_[iSequence+1]-1;
     330    // correct costs for this infeasibility weight
     331    if (infeasible(start))
     332      cost_[start] = cost_[start+1]-infeasibilityCost;
     333    if (infeasible(end-1))
     334      cost_[end-1] = cost_[end-2]+infeasibilityCost;
    321335    for (iRange=start; iRange<end;iRange++) {
    322336      if (value<lower_[iRange+1]+primalTolerance) {
     
    329343    }
    330344    assert(iRange<end);
     345    lowerValue = lower_[iRange];
     346    upperValue = lower_[iRange+1];
    331347    switch(model_->getStatus(iSequence)) {
    332348     
     
    335351      // iRange is in correct place
    336352      // slot in here
    337       if (value<lower[iSequence]-primalTolerance) {
    338         value = lower[iSequence]-value;
    339         sumInfeasibilities_ += value;
    340         largestInfeasibility_ = max(largestInfeasibility_,value);
    341         changeCost_ -= lower[iSequence]*
    342           (cost_[iRange]-cost[iSequence]);
    343         numberInfeasibilities_++;
    344       } else if (value>upper[iSequence]+primalTolerance) {
    345         value = value-upper[iSequence];
    346         sumInfeasibilities_ += value;
    347         largestInfeasibility_ = max(largestInfeasibility_,value);
    348         changeCost_ -= upper[iSequence]*
    349           (cost_[iRange]-cost[iSequence]);
    350         numberInfeasibilities_++;
     353      if (infeasible(iRange)) {
     354        if (lower_[iRange]<-1.0e50) {
     355          //cost_[iRange] = cost_[iRange+1]-infeasibilityCost;
     356          // possibly below
     357          lowerValue = lower_[iRange+1];
     358          if (value<lowerValue-primalTolerance) {
     359            value = lowerValue-value;
     360            sumInfeasibilities_ += value;
     361            largestInfeasibility_ = max(largestInfeasibility_,value);
     362            changeCost_ -= lowerValue*
     363              (cost_[iRange]-cost[iSequence]);
     364            numberInfeasibilities_++;
     365          }
     366        } else {
     367          //cost_[iRange] = cost_[iRange-1]+infeasibilityCost;
     368          // possibly above
     369          upperValue = lower_[iRange];
     370          if (value>upperValue+primalTolerance) {
     371            value = value-upperValue;
     372            sumInfeasibilities_ += value;
     373            largestInfeasibility_ = max(largestInfeasibility_,value);
     374            changeCost_ -= upperValue*
     375              (cost_[iRange]-cost[iSequence]);
     376            numberInfeasibilities_++;
     377          }
     378        }
    351379      }
    352380      lower[iSequence] = lower_[iRange];
     
    366394        }
    367395      } else {
    368         solution[iSequence] = upperValue;
     396        if (fabs(value-upperValue)<=fabs(value-lowerValue)) {
     397          solution[iSequence] = upperValue;
     398        } else {
     399          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
     400          solution[iSequence] = lowerValue;
     401        }
    369402      }
    370403      break;
     
    377410        }
    378411      } else {
    379         solution[iSequence] = lowerValue;
     412        if (fabs(value-lowerValue)<=fabs(value-upperValue)) {
     413          solution[iSequence] = lowerValue;
     414        } else {
     415          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
     416          solution[iSequence] = upperValue;
     417        }
    380418      }
    381419      break;
  • trunk/ClpSimplexDual.cpp

    r52 r54  
    618618          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
    619619#endif
    620 
    621620        if (whatNext==1) {
    622621          problemStatus_ =-2; // refactorize
     
    870869  pivotRow_=dualRowPivot_->pivotRow();
    871870  if (pivotRow_>=0) {
    872     int iPivot=pivotVariable_[pivotRow_];
    873     sequenceOut_ = iPivot;
    874     if (iPivot>=numberColumns_) {
    875       // slack
    876       iPivot-=numberColumns_;
    877       valueOut_=rowActivityWork_[iPivot];
    878       lowerOut_=rowLowerWork_[iPivot];
    879       upperOut_=rowUpperWork_[iPivot];
    880     } else {
    881       // column
    882       valueOut_=columnActivityWork_[iPivot];
    883       lowerOut_=columnLowerWork_[iPivot];
    884       upperOut_=columnUpperWork_[iPivot];
    885     }
     871    sequenceOut_ = pivotVariable_[pivotRow_];
     872    valueOut_ = solution_[sequenceOut_];
     873    lowerOut_ = lower_[sequenceOut_];
     874    upperOut_ = upper_[sequenceOut_];
    886875    // if we have problems we could try other way and hope we get a
    887876    // zero pivot?
     
    11781167        break;
    11791168      case atUpperBound:
    1180 #ifdef CHECK_DJ
    1181         // For Debug so we can find out where problem is
    1182         perturbation_ = iSequence+addSequence;
    1183 #endif
    11841169        assert (oldValue<=dualTolerance_*1.0001);
    11851170        if (value>newTolerance) {
     
    11911176        break;
    11921177      case atLowerBound:
    1193 #ifdef CHECK_DJ
    1194         // For Debug so we can find out where problem is
    1195         perturbation_ = iSequence+addSequence;
    1196 #endif
    11971178        assert (oldValue>=-dualTolerance_*1.0001);
    11981179        if (value<-newTolerance) {
     
    13681349                if (-alpha>=acceptablePivot) {
    13691350                  upperTheta = (oldValue-newTolerance)/alpha;
     1351                  // recompute value and make sure works
     1352                  value = oldValue-upperTheta*alpha;
     1353                  if (value<0)
     1354                    upperTheta *= 1.0 +1.0e-11; // must be large
    13701355                }
    13711356              }
     
    13751360                if (alpha>=acceptablePivot) {
    13761361                  upperTheta = (oldValue+newTolerance)/alpha;
     1362                  // recompute value and make sure works
     1363                  value = oldValue-upperTheta*alpha;
     1364                  if (value>0)
     1365                    upperTheta *= 1.0 +1.0e-11; // must be large
    13771366                }
    13781367              }
     
    18611850        if (numberChangedBounds<=0) {
    18621851          //looks optimal - do we need to reset tolerance
     1852          if (perturbation_==101) {
     1853            perturbation_=102; // stop any perturbations
     1854            cleanDuals=true;
     1855            createRim(4);
     1856          }
    18631857          if (lastCleaned<numberIterations_&&numberTimesOptimal_<4) {
    18641858            doOriginalTolerance=2;
     
    18841878            problemStatus_ = checkUnbounded(rowArray_[0],rowArray_[1],
    18851879                                            changeCost);
     1880            if (problemStatus_==2&&perturbation_==101) {
     1881              perturbation_=102; // stop any perturbations
     1882              cleanDuals=true;
     1883              createRim(4);
     1884              problemStatus_=-1;
     1885            }
    18861886            if (problemStatus_==2) {
    18871887              // it is unbounded - restore solution
     
    19121912          <<CoinMessageEol;
    19131913        numberChangedBounds=changeBounds(false,NULL,changeCost);
     1914        if (perturbation_==101) {
     1915          perturbation_=102; // stop any perturbations
     1916          cleanDuals=true;
     1917          numberChangedBounds=1;
     1918          createRim(4);
     1919        }
    19141920        if (numberChangedBounds<=0||dualBound_>1.0e20||
    19151921            (largestPrimalError_>1.0&&dualBound_>1.0e17)) {
     
    21452151  // maximum fraction of cost to perturb
    21462152  double maximumFraction = 1.0e-4;
    2147   if (perturbation_==100) {
     2153  if (perturbation_>=50) {
    21482154    perturbation = 1.0e-4;
    21492155    for (iRow=0;iRow<numberRows_;iRow++) {
    2150       double value = fabs(rowActivityWork_[iRow]*rowObjectiveWork_[iRow]);
     2156      double value = fabs(rowObjectiveWork_[iRow]);
    21512157      perturbation = max(perturbation,value);
    21522158    }
    21532159    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    21542160      double value =
    2155         fabs(columnActivityWork_[iColumn]*objectiveWork_[iColumn]);
     2161        fabs(objectiveWork_[iColumn]);
    21562162      perturbation = max(perturbation,value);
    21572163    }
     
    22002206  }
    22012207  // say perturbed
    2202   perturbation_=102;
     2208  perturbation_=101;
    22032209
    22042210}
  • trunk/ClpSimplexPrimal.cpp

    r52 r54  
    5353      3) Degeneracy.  Gill et al helps but may not be enough.  We
    5454      may need more.  Also it can improve speed a lot if we perturb
    55       the costs significantly. 
     55      the rhs and bounds significantly. 
    5656
    5757  References:
     
    221221   
    222222   
     223   
    223224    // If user asked for perturbation - do it
    224225   
    225226    if (perturbation_<100)
    226227      perturb();
    227    
     228
    228229    // for primal we will change bounds using infeasibilityCost_
    229230    if (nonLinearCost_==NULL) {
     
    668669  }
    669670  int tentativeStatus = problemStatus_;
    670 
    671671  if (problemStatus_>-3||problemStatus_==-4) {
    672672    // factorize
     
    698698  // at this stage status is -3 or -5 if looks unbounded
    699699  // get primal and dual solutions
    700   // put back original bounds and then check
    701   createRim(7);
     700  // put back original costs and then check
     701  createRim(4);
    702702  gutsOfSolution(rowActivityWork_, columnActivityWork_);
    703703  // Check if looping
     
    736736  }
    737737  if (dualFeasible()||problemStatus_==-4||(type==3&&problemStatus_!=-5)) {
     738   
    738739    if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
    739740      //may need infeasiblity cost changed
     
    745746      // do twice to make sure Primal solution has settled
    746747      // put non-basics to bounds in case tolerance moved
    747       // put back original bounds
    748       createRim(7);
     748      // put back original costs
     749      createRim(4);
    749750      nonLinearCost_->checkInfeasibilities(true);
    750751      gutsOfSolution(rowActivityWork_, columnActivityWork_);
    751752
    752753      infeasibilityCost_=1.0e100;
    753       // put back original bounds
    754       createRim(7);
     754      // put back original costs
     755      createRim(4);
    755756      nonLinearCost_->checkInfeasibilities(true);
    756757      nonLinearCost_=NULL;
     
    762763      nonLinearCost_=nonLinear;
    763764      infeasibilityCost_=saveWeight;
     765      if ((infeasibilityCost_>=1.0e18||
     766          numberDualInfeasibilities_==0)&&perturbation_==101) {
     767        unPerturb(); // stop any further perturbation
     768        numberDualInfeasibilities_=1; // carry on
     769      }
    764770      if (infeasibilityCost_>=1.0e20||
    765771          numberDualInfeasibilities_==0) {
     
    770776        // and get feasible duals
    771777        infeasibilityCost_=0.0;
    772         createRim(7);
     778        createRim(4);
    773779        nonLinearCost_->checkInfeasibilities(true);
    774780        gutsOfSolution(rowActivityWork_, columnActivityWork_);
     
    783789          <<infeasibilityCost_
    784790          <<CoinMessageEol;
    785         // put back original bounds and then check
    786         createRim(7);
     791        // put back original costs and then check
     792        createRim(4);
    787793        nonLinearCost_->checkInfeasibilities(true);
    788794        gutsOfSolution(rowActivityWork_, columnActivityWork_);
     
    794800    } else {
    795801      // may be optimal
    796       if ( lastCleaned!=numberIterations_) {
     802      if (perturbation_==101) {
     803        unPerturb(); // stop any further perturbation
     804        lastCleaned=-1; // carry on
     805      }
     806      if ( lastCleaned!=numberIterations_||unflag()) {
    797807        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
    798808          <<primalTolerance_
     
    810820          primalTolerance_=dblParam_[ClpPrimalTolerance];
    811821         
    812           // put back original bounds and then check
    813           createRim(7);
     822          // put back original costs and then check
     823          createRim(4);
    814824          nonLinearCost_->checkInfeasibilities(true);
    815825          gutsOfSolution(rowActivityWork_, columnActivityWork_);
     
    830840    if (problemStatus_==-5) {
    831841      if (nonLinearCost_->numberInfeasibilities()) {
     842        if (infeasibilityCost_>1.0e18&&perturbation_==101) {
     843          // back off weight
     844          infeasibilityCost_ = 1.0e13;
     845          unPerturb(); // stop any further perturbation
     846        }
    832847        //we need infeasiblity cost changed
    833848        if (infeasibilityCost_<1.0e20) {
     
    837852            <<infeasibilityCost_
    838853            <<CoinMessageEol;
    839           // put back original bounds and then check
    840           createRim(7);
     854          // put back original costs and then check
     855          createRim(4);
    841856          gutsOfSolution(rowActivityWork_, columnActivityWork_);
    842857          problemStatus_=-1; //continue
     
    14171432  if (perturbation_>100)
    14181433    return; //perturbed already
    1419   abort();
     1434  int i;
     1435  // primal perturbation
     1436  double perturbation=1.0e-20;
     1437  // maximum fraction of rhs/bounds to perturb
     1438  double maximumFraction = 1.0e-4;
     1439  if (perturbation_>=50) {
     1440    perturbation = 1.0e-4;
     1441    for (i=0;i<numberColumns_+numberRows_;i++) {
     1442      if (upper_[i]>lower_[i]+primalTolerance_) {
     1443        double lowerValue, upperValue;
     1444        if (lower_[i]>-1.0e20)
     1445          lowerValue = fabs(lower_[i]);
     1446        else
     1447          lowerValue=0.0;
     1448        if (upper_[i]<1.0e20)
     1449          upperValue = fabs(upper_[i]);
     1450        else
     1451          upperValue=0.0;
     1452        double value = max(lowerValue,upperValue);
     1453        perturbation = max(perturbation,value);
     1454      }
     1455    }
     1456    perturbation *= 1.0e-8;
     1457  } else if (perturbation_<100) {
     1458    perturbation = pow(10.0,perturbation_);
     1459    // user is in charge
     1460    maximumFraction = 1.0e100;
     1461  }
     1462  // modify bounds
     1463  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
     1464    <<perturbation
     1465    <<CoinMessageEol;
     1466  for (i=0;i<numberColumns_+numberRows_;i++) {
     1467    double lowerValue=lower_[i], upperValue=upper_[i];
     1468    if (upperValue>lowerValue+primalTolerance_) {
     1469      double value = drand48()*perturbation;
     1470      if (lowerValue>-1.0e20)
     1471        lowerValue -= value * (max(1.0,1.0e-3*fabs(lowerValue)));
     1472      if (upperValue<1.0e20)
     1473        upperValue += value * (max(1.0,1.0e-3*fabs(upperValue)));
     1474    }
     1475    lower_[i]=lowerValue;
     1476    upper_[i]=upperValue;
     1477  }
     1478  // say perturbed
     1479  perturbation_=101;
     1480}
     1481// un perturb
     1482void
     1483ClpSimplexPrimal::unPerturb()
     1484{
     1485  // put back original bounds and costs
     1486  createRim(7);
     1487  // unflag
     1488  unflag();
     1489  // get a valid nonlinear cost function
     1490  delete nonLinearCost_;
     1491  nonLinearCost_= new ClpNonLinearCost(this);
     1492  perturbation_ = 102; // stop any further perturbation
     1493  // move non basic variables to new bounds
     1494  nonLinearCost_->checkInfeasibilities(true);
     1495  gutsOfSolution(rowActivityWork_, columnActivityWork_);
     1496}
     1497// Unflag all variables and return number unflagged
     1498int
     1499ClpSimplexPrimal::unflag()
     1500{
     1501  int i;
     1502  int number = numberRows_+numberColumns_;
     1503  int numberFlagged=0;
     1504  for (i=0;i<number;i++) {
     1505    if (flagged(i)) {
     1506      clearFlagged(i);
     1507      numberFlagged++;
     1508    }
     1509  }
     1510  return numberFlagged;
    14201511}
    14211512// Do not change infeasibility cost and always say optimal
  • trunk/include/ClpNonLinearCost.hpp

    r50 r54  
    143143  {changeCost_ = value;};
    144144  //@}
     145  ///@name Private functions to deal with infeasible regions
     146  inline bool infeasible(int i) const {
     147    return ((infeasible_[i>>5]>>(i&31))&1)!=0;
     148  }
     149  inline void setInfeasible(int i,bool trueFalse) {
     150    unsigned int & value = infeasible_[i>>5];
     151    int bit = i&31;
     152    if (trueFalse)
     153      value |= (1<<bit);
     154    else
     155      value &= ~(1<<bit);
     156  }
     157  //@}
    145158   
    146159private:
     
    175188  /// If all non-linear costs convex
    176189  bool convex_;
     190  // Array to say which regions are infeasible
     191  unsigned int * infeasible_;
    177192  //@}
    178193};
  • trunk/include/ClpSimplexPrimal.hpp

    r50 r54  
    187187  /// Perturbs problem (method depends on perturbation())
    188188  void perturb();
     189  /// Take off effect of perturbation
     190  void unPerturb();
     191  /// Unflag all variables and return number unflagged
     192  int unflag();
    189193  /// Sets sequenceIn_ to next superBasic (input by first..) and updates
    190194  void nextSuperBasic(int & firstSuperBasic);
Note: See TracChangeset for help on using the changeset viewer.