Changeset 1884 for trunk


Ignore:
Timestamp:
Nov 5, 2012 12:38:48 PM (7 years ago)
Author:
forrest
Message:

minor update to improve accuarcy in parametrics

Location:
trunk/Clp/src
Files:
2 edited

Legend:

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

    r1878 r1884  
    28762876  // allow for unscaled - even if not needed
    28772877  int lengthArrays = 4*numberTotal+(3*numberTotal+2)*ratio+2*numberRows_+1;
     2878  int unscaledChangesOffset=lengthArrays; // need extra copy of change
     2879  lengthArrays += numberTotal;
     2880 
    28782881  /*
    28792882    Save information and modify
     
    29102913  paramData.upperGap = upperGap;
    29112914  paramData.upperCoefficient = upperCoefficient;
     2915  paramData.unscaledChangesOffset = unscaledChangesOffset-numberTotal;
     2916  paramData.firstIteration=true;
    29122917  // Find theta when bounds will cross over and create arrays
    29132918  memset(lowerChange, 0, numberTotal * sizeof(double));
     
    29232928    memcpy(upperChange+numberColumns_,
    29242929           upperChangeRhs,numberRows_*sizeof(double));
     2930  // clean
     2931  for (int iRow = 0; iRow < numberRows_; iRow++) {
     2932    double lower = rowLower_[iRow];
     2933    double upper = rowUpper_[iRow];
     2934    if (lower<-1.0e30)
     2935      lowerChange[numberColumns_+iRow]=0.0;
     2936    if (upper>1.0e30)
     2937      upperChange[numberColumns_+iRow]=0.0;
     2938  }
     2939  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2940    double lower = columnLower_[iColumn];
     2941    double upper = columnUpper_[iColumn];
     2942    if (lower<-1.0e30)
     2943      lowerChange[iColumn]=0.0;
     2944    if (upper>1.0e30)
     2945      upperChange[iColumn]=0.0;
     2946  }
     2947  // save unscaled version of changes
     2948  memcpy(saveLower+unscaledChangesOffset,lowerChange,numberTotal*sizeof(double));
     2949  memcpy(saveUpper+unscaledChangesOffset,upperChange,numberTotal*sizeof(double));
    29252950  int nLowerChange=0;
    29262951  int nUpperChange=0;
     
    29632988           rowUpper_,numberRows_*sizeof(double));
    29642989  }
    2965   double maxTheta = 1.0e50;
    2966   double largestChange=0.0;
    2967   for (int iRow = 0; iRow < numberRows_; iRow++) {
    2968     double lower = rowLower_[iRow];
    2969     double upper = rowUpper_[iRow];
    2970     if (lower<-1.0e30)
    2971       lowerChange[numberColumns_+iRow]=0.0;
    2972     double chgLower = lowerChange[numberColumns_+iRow];
    2973     largestChange=CoinMax(largestChange,fabs(chgLower));
    2974     if (upper>1.0e30)
    2975       upperChange[numberColumns_+iRow]=0.0;
    2976     double chgUpper = upperChange[numberColumns_+iRow];
    2977     largestChange=CoinMax(largestChange,fabs(chgUpper));
    2978     if (lower > -1.0e30 && upper < 1.0e30) {
    2979       if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
    2980         maxTheta = (upper - lower) / (chgLower - chgUpper);
    2981       }
    2982     }
    2983     lower+=startingTheta*chgLower;
    2984     upper+=startingTheta*chgUpper;
    2985 #ifndef CLP_USER_DRIVEN
    2986     if (lower > upper) {
    2987       maxTheta = -1.0;
    2988       break;
    2989     }
    2990 #endif
    2991     rowLower_[iRow]=lower;
    2992     rowUpper_[iRow]=upper;
    2993   }
    2994   for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2995     double lower = columnLower_[iColumn];
    2996     double upper = columnUpper_[iColumn];
    2997     if (lower<-1.0e30)
    2998       lowerChange[iColumn]=0.0;
    2999     double chgLower = lowerChange[iColumn];
    3000     largestChange=CoinMax(largestChange,fabs(chgLower));
    3001     if (upper>1.0e30)
    3002       upperChange[iColumn]=0.0;
    3003     double chgUpper = upperChange[iColumn];
    3004     largestChange=CoinMax(largestChange,fabs(chgUpper));
    3005     if (lower > -1.0e30 && upper < 1.0e30) {
    3006       if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
    3007         maxTheta = (upper - lower) / (chgLower - chgUpper);
    3008       }
    3009     }
    3010     lower+=startingTheta*chgLower;
    3011     upper+=startingTheta*chgUpper;
    3012 #ifndef CLP_USER_DRIVEN
    3013     if (lower > upper) {
    3014       maxTheta = -1.0;
    3015       break;
    3016     }
    3017 #endif
    3018     columnLower_[iColumn]=lower;
    3019     columnUpper_[iColumn]=upper;
    3020   }
    3021   if (maxTheta == 1.0e50)
    3022     maxTheta = COIN_DBL_MAX;
    30232990  int returnCode=0;
    3024   /*
    3025     If user is in charge then they are sophisticated enough
    3026     to know what they are doing.  They must be doing something clever
    3027     with event handlers!
    3028    */
    3029 #ifndef CLP_USER_DRIVEN
    3030   if (maxTheta < 0.0) {
    3031     // bad ranges or initial
    3032     returnCode = -1;
    3033   }
    3034   if (maxTheta < endingTheta) {
    3035     char line[100];
    3036     sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n",
    3037             endingTheta,maxTheta);
    3038     handler_->message(CLP_GENERAL,messages_)
    3039       << line << CoinMessageEol;
    3040     endingTheta = maxTheta;
    3041   }
    3042   if (endingTheta < startingTheta) {
    3043     // bad initial
    3044     returnCode = -2;
    3045   }
    3046   paramData.maxTheta=CoinMin(maxTheta,endingTheta);
    3047 #else
    3048   // accept user's version
     2991  paramData.startingTheta=startingTheta;
     2992  paramData.endingTheta=endingTheta;
    30492993  paramData.maxTheta=endingTheta;
    3050 #endif
    3051   /* given largest change element choose acceptable end
    3052      be safe and make sure difference < 0.1*tolerance */
    3053   double acceptableDifference=0.1*primalTolerance_/
    3054     CoinMax(largestChange,1.0);
    3055   paramData.acceptableMaxTheta=maxTheta-acceptableDifference;
     2994  computeRhsEtc(paramData);
    30562995  bool swapped=false;
    30572996  // Dantzig
     
    36473586          useTheta += theta_;
    36483587          double change = useTheta - lastTheta;
     3588          if (paramData.firstIteration) {
     3589            // redo rhs etc to make as accurate as possible
     3590            paramData.firstIteration=false;
     3591            if (change>1.0e-14) {
     3592              startingTheta=useTheta;
     3593              lastTheta=startingTheta;
     3594              change=0.0;
     3595              // restore rhs
     3596              const double * saveLower = paramData.lowerChange+2*numberTotal;
     3597              memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
     3598              memcpy(rowLower_,saveLower+numberColumns_,numberRows_*sizeof(double));
     3599              const double * saveUpper = paramData.upperChange+2*numberTotal;
     3600              memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
     3601              memcpy(rowUpper_,saveUpper+numberColumns_,numberRows_*sizeof(double));
     3602              paramData.startingTheta=startingTheta;
     3603              computeRhsEtc(paramData);
     3604              redoInternalArrays();
     3605              // Update solution
     3606              rowArray_[4]->clear();
     3607              for (int i=0;i<numberTotal;i++) {
     3608                if (getStatus(i)==atLowerBound||getStatus(i)==isFixed)
     3609                  solution_[i]=lower_[i];
     3610                else if (getStatus(i)==atUpperBound)
     3611                  solution_[i]=upper_[i];
     3612              }
     3613              gutsOfSolution(NULL,NULL);
     3614            }
     3615          }
    36493616          if (change>1.0e-14) {
    36503617            int n;
     
    44294396     startingTheta = lastTheta+theta_;
    44304397     return returnCode;
     4398}
     4399// Compute new rowLower_ etc (return negative if infeasible - otherwise largest change)
     4400double
     4401ClpSimplexOther::computeRhsEtc( parametricsData & paramData)
     4402{
     4403  double maxTheta = COIN_DBL_MAX;
     4404  double largestChange=0.0;
     4405  double startingTheta = paramData.startingTheta;
     4406  const double * lowerChange = paramData.lowerChange+
     4407    paramData.unscaledChangesOffset;
     4408  const double * upperChange = paramData.upperChange+
     4409    paramData.unscaledChangesOffset;
     4410  for (int iRow = 0; iRow < numberRows_; iRow++) {
     4411    double lower = rowLower_[iRow];
     4412    double upper = rowUpper_[iRow];
     4413    double chgLower = lowerChange[numberColumns_+iRow];
     4414    largestChange=CoinMax(largestChange,fabs(chgLower));
     4415    double chgUpper = upperChange[numberColumns_+iRow];
     4416    largestChange=CoinMax(largestChange,fabs(chgUpper));
     4417    if (lower > -1.0e30 && upper < 1.0e30) {
     4418      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
     4419        maxTheta = (upper - lower) / (chgLower - chgUpper);
     4420      }
     4421    }
     4422    lower+=startingTheta*chgLower;
     4423    upper+=startingTheta*chgUpper;
     4424#ifndef CLP_USER_DRIVEN
     4425    if (lower > upper) {
     4426      maxTheta = -1.0;
     4427      break;
     4428    }
     4429#endif
     4430    rowLower_[iRow]=lower;
     4431    rowUpper_[iRow]=upper;
     4432  }
     4433  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     4434    double lower = columnLower_[iColumn];
     4435    double upper = columnUpper_[iColumn];
     4436    double chgLower = lowerChange[iColumn];
     4437    largestChange=CoinMax(largestChange,fabs(chgLower));
     4438    double chgUpper = upperChange[iColumn];
     4439    largestChange=CoinMax(largestChange,fabs(chgUpper));
     4440    if (lower > -1.0e30 && upper < 1.0e30) {
     4441      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
     4442        maxTheta = (upper - lower) / (chgLower - chgUpper);
     4443      }
     4444    }
     4445    lower+=startingTheta*chgLower;
     4446    upper+=startingTheta*chgUpper;
     4447#ifndef CLP_USER_DRIVEN
     4448    if (lower > upper) {
     4449      maxTheta = -1.0;
     4450      break;
     4451    }
     4452#endif
     4453    columnLower_[iColumn]=lower;
     4454    columnUpper_[iColumn]=upper;
     4455  }
     4456#ifndef CLP_USER_DRIVEN
     4457  paramData.maxTheta=maxTheta;
     4458  if (maxTheta<0)
     4459    largestChange=-1.0; // signal infeasible
     4460#else
     4461  // maxTheta already set
     4462  /* given largest change element choose acceptable end
     4463     be safe and make sure difference < 0.1*tolerance */
     4464  double acceptableDifference=0.1*primalTolerance_/
     4465    CoinMax(largestChange,1.0);
     4466  paramData.acceptableMaxTheta=maxTheta-acceptableDifference;
     4467#endif
     4468  return largestChange;
     4469}
     4470// Redo lower_ from rowLower_ etc
     4471void
     4472ClpSimplexOther::redoInternalArrays()
     4473{
     4474  double * lowerSave = lower_;
     4475  double * upperSave = upper_;
     4476  memcpy(lowerSave,columnLower_,numberColumns_*sizeof(double));
     4477  memcpy(lowerSave+numberColumns_,rowLower_,numberRows_*sizeof(double));
     4478  memcpy(upperSave,columnUpper_,numberColumns_*sizeof(double));
     4479  memcpy(upperSave+numberColumns_,rowUpper_,numberRows_*sizeof(double));
     4480  if (rowScale_) {
     4481    // scale arrays
     4482    for (int i=0;i<numberColumns_;i++) {
     4483      double multiplier = inverseColumnScale_[i];
     4484      if (lowerSave[i]>-1.0e20)
     4485        lowerSave[i] *= multiplier;
     4486      if (upperSave[i]<1.0e20)
     4487        upperSave[i] *= multiplier;
     4488    }
     4489    lowerSave += numberColumns_;
     4490    upperSave += numberColumns_;
     4491    for (int i=0;i<numberRows_;i++) {
     4492      double multiplier = rowScale_[i];
     4493      if (lowerSave[i]>-1.0e20)
     4494        lowerSave[i] *= multiplier;
     4495      if (upperSave[i]<1.0e20)
     4496        upperSave[i] *= multiplier;
     4497    }
     4498  }
    44314499}
    44324500#if 0
  • trunk/Clp/src/ClpSimplexOther.hpp

    r1860 r1884  
    122122    double * upperGap;
    123123    double * upperCoefficient;
     124    int unscaledChangesOffset;
     125    bool firstIteration; // so can update rhs for accuracy
    124126  } parametricsData;
    125127
     
    173175     void originalBound(int iSequence, double theta, const double * changeLower,
    174176                     const double * changeUpper);
     177     /// Compute new rowLower_ etc (return negative if infeasible - otherwise largest change)
     178     double computeRhsEtc(parametricsData & paramData);
     179     /// Redo lower_ from rowLower_ etc
     180     void redoInternalArrays();
    175181     /**
    176182         Row array has row part of pivot row
Note: See TracChangeset for help on using the changeset viewer.