Changeset 364


Ignore:
Timestamp:
May 11, 2004 3:36:06 PM (14 years ago)
Author:
forrest
Message:

better PD but can't get centring right

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpPredictorCorrector.cpp

    r362 r364  
    9999  double * savePrimal=NULL;
    100100  int numberTotal = numberRows_+numberColumns_;
     101  // Extra regions for centering
     102  double * saveWeights = new double[numberTotal];
     103  double * saveZ = new double[numberTotal];
     104  double * saveW = new double[numberTotal];
    101105  while (problemStatus_<0) {
    102106    if (numberIterations_==xxxxxx) {
     
    309313    bool goodMove=false;
    310314    double bestNextGap=COIN_DBL_MAX;
     315    double nextGap=bestNextGap;
     316    int nextNumber=0;
    311317    worstDirectionAccuracy_=0.0;
    312318    while (!goodMove) {
     
    370376          worstDirectionAccuracy_=directionAccuracy;
    371377        }
    372         findStepLength(phase);
    373         int nextNumber; //number of complementarity pairs
    374         double nextGap=complementarityGap(nextNumber,1);
     378        findStepLength(phase,NULL);
     379        nextGap=complementarityGap(nextNumber,1);
    375380       
    376381        bestNextGap=max(nextGap,0.8*complementarityGap_);
     
    473478        }
    474479        if (goodMove) {
    475           findStepLength(phase);
     480          findStepLength(phase,NULL);
    476481          // just for debug
    477           int nextNumber;
    478           complementarityGap(nextNumber,1);
     482          nextGap = complementarityGap(nextNumber,1);
    479483          if (numberIterations_>=-77) {
    480484            goodMove=checkGoodMove(true,bestNextGap);
     
    507511          testValue=maximumRHSError_;
    508512        }
    509         findStepLength(2);
     513        findStepLength(2,NULL);
    510514        // just for debug
    511         int nextNumber;
    512         complementarityGap(nextNumber,2);
     515        nextGap=complementarityGap(nextNumber,2);
    513516      }
    514517    } /* endwhile */
     
    554557      }
    555558    }
     559#if 0
     560    // See if we can do centering steps
     561    memcpy(saveWeights,weights_,numberTotal*sizeof(double));
     562    memcpy(saveZ,deltaZ_,numberTotal*sizeof(double));
     563    memcpy(saveW,deltaW_,numberTotal*sizeof(double));
     564    double savePrimalStep = actualPrimalStep_;
     565    double saveDualStep = actualDualStep_;
     566    double saveMu = mu_;
     567    mu_=nextGap / (double) nextNumber;
     568    setupForSolve(3);
     569    findDirectionVector(3);
     570    memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
     571    memcpy(deltaW_,saveW,numberTotal*sizeof(double));
     572    findStepLength(3,saveWeights);
     573    double xGap = complementarityGap(nextNumber,3);
     574#if 0
     575    goodMove=checkGoodMove(true,bestNextGap);
     576#endif
     577    if (xGap>nextGap) {
     578      mu_=saveMu;
     579      actualPrimalStep_ = savePrimalStep;
     580      actualDualStep_ = saveDualStep;
     581      memcpy(weights_,saveWeights,numberTotal*sizeof(double));
     582      memcpy(deltaZ_,saveZ,numberTotal*sizeof(double));
     583      memcpy(deltaW_,saveW,numberTotal*sizeof(double));
     584    }
     585#endif
    556586    numberFixed=updateSolution();
    557587    numberFixedTotal+=numberFixed;
    558588  } /* endwhile */
     589  delete [] saveWeights;
     590  delete [] saveZ;
     591  delete [] saveW;
    559592  if (savePi) {
    560593    //std::cout<<"Restoring from iteration "<<saveIteration<<std::endl;
     
    595628//         1 corrector
    596629//         2 primal dual
    597 double ClpPredictorCorrector::findStepLength(const int phase)
     630double ClpPredictorCorrector::findStepLength(const int phase,const double * oldWeight)
    598631{
    599632  double directionNorm=0.0;
     
    870903    }
    871904    break;
    872   }
    873   printf("step - phase %d, norm %g, dual step %g, primal step %g\n",
    874          phase,directionNorm,maximumPrimalStep,maximumDualStep);
     905  case 3:
     906    //centering step
     907    for (int iColumn=0;iColumn<numberTotal;iColumn++) {
     908      if (!flagged(iColumn)) {
     909        // old deltas from previous step
     910        double z1=work1[iColumn];
     911        double w1=work2[iColumn];
     912        double work3Value=0.0;
     913        double work4Value=0.0;
     914        double directionElement=weights[iColumn];
     915        double oldElement = oldWeight[iColumn];
     916        double newElement = oldElement+directionElement;
     917        weights[iColumn]=newElement;
     918        double value=primal[iColumn];
     919        if (directionNorm<fabs(directionElement)) {
     920          directionNorm=fabs(directionElement);
     921        }
     922        if (lowerBound(iColumn)||
     923                        upperBound(iColumn)) {
     924          if (lowerBound(iColumn)) {
     925            double gap=lowerSlack[iColumn]+extra;
     926            double delta;
     927            double zDelta;
     928            if (!fakeLower(iColumn)) {
     929              zDelta = lower[iColumn]+lowerSlack[iColumn]
     930             -  value - directionElement;
     931            } else {
     932              zDelta = - directionElement;
     933            }
     934            delta = zDelta - oldElement;
     935            z1+=(work3[iColumn]-zVec[iColumn]*zDelta)/gap;
     936            work3Value=-delta;
     937            if (lowerSlack[iColumn]<maximumPrimalStep*delta) {
     938              maximumPrimalStep=lowerSlack[iColumn]/delta;
     939              chosenPrimalSequence=iColumn;
     940            }
     941            if (zVec[iColumn]>tolerance) {
     942              if (zVec[iColumn]<-z1*maximumDualStep) {
     943                maximumDualStep=-zVec[iColumn]/z1;
     944                chosenDualSequence=iColumn;
     945              }
     946            }
     947          }
     948          if (upperBound(iColumn)) {
     949            double gap=upperSlack[iColumn]+extra;
     950            double delta;
     951            double zDelta;
     952            if (!fakeUpper(iColumn)) {
     953              zDelta = -upper[iColumn]+upperSlack[iColumn]
     954                        + value + directionElement;
     955            } else {
     956              zDelta = directionElement;
     957            }
     958            delta = zDelta + oldElement;
     959            //double delta = -upper[iColumn]+upperSlack[iColumn]-
     960              //+ value + directionElement;
     961            w1+=(work4[iColumn]+wVec[iColumn]*zDelta)/gap;
     962            work4Value=-delta;
     963            if (upperSlack[iColumn]<maximumPrimalStep*delta) {
     964              maximumPrimalStep=upperSlack[iColumn]/delta;
     965              chosenPrimalSequence=iColumn;
     966            }
     967            if (wVec[iColumn]>tolerance) {
     968              if (wVec[iColumn]<-w1*maximumDualStep) {
     969                maximumDualStep=-wVec[iColumn]/w1;
     970                chosenDualSequence=iColumn;
     971              }
     972            }
     973          }
     974        } else {
     975          //free
     976          double gap=fabs(value);
     977          double multiplier=1.0/gap;
     978          if (gap<1.0) {
     979            multiplier=1,0;
     980          }
     981          z1+=multiplier*(work3[iColumn]-directionElement*zVec[iColumn]);
     982          w1+=multiplier*(work4[iColumn]+directionElement*wVec[iColumn]);
     983        }
     984        work1[iColumn]=z1;
     985        work2[iColumn]=w1;
     986        work3[iColumn]=work3Value;
     987        work4[iColumn]=work4Value;
     988      } else {
     989        work1[iColumn]=0.0;
     990        work2[iColumn]=0.0;
     991        work3[iColumn]=0.0;
     992        work4[iColumn]=0.0;
     993      }
     994    }
     995    break;
     996  }
     997  //printf("step - phase %d, norm %g, dual step %g, primal step %g\n",
     998  // phase,directionNorm,maximumPrimalStep,maximumDualStep);
    875999  actualPrimalStep_=stepLength_*maximumPrimalStep;
    8761000  if (phase>=0&&actualPrimalStep_>1.0) {
     
    9561080double ClpPredictorCorrector::findDirectionVector(const int phase)
    9571081{
     1082  if (phase==3) {
     1083    // centering step
     1084    double * work2 = deltaW_;
     1085    int numberTotal = numberRows_+numberColumns_;
     1086    double * tempRegion = new double [numberRows_];
     1087    double * newError = new double [numberRows_];
     1088    // For KKT separate out
     1089#ifndef KKT
     1090    int iColumn;
     1091    for (iColumn=0;iColumn<numberTotal;iColumn++)
     1092      weights_[iColumn] = work2[iColumn];
     1093    multiplyAdd(weights_+numberColumns_,numberRows_,-1.0,tempRegion,0.0);
     1094    matrix_->times(1.0,weights_,tempRegion);
     1095#else
     1096    // regions in will be work2 and newError
     1097    // regions out weights_ and tempRegion
     1098    memset(newError,0,numberRows_*sizeof(double));
     1099    double * region1Save=NULL;//for refinement
     1100#endif
     1101#ifndef KKT
     1102    double maximumRHS = maximumAbsElement(tempRegion,numberRows_);
     1103    double scale=1.0;
     1104    double unscale=1.0;
     1105    if (maximumRHS>1.0e-30) {
     1106      if (maximumRHS<=0.5) {
     1107        double factor=2.0;
     1108        while (maximumRHS<=0.5) {
     1109          maximumRHS*=factor;
     1110          scale*=factor;
     1111        } /* endwhile */
     1112      } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) {
     1113        double factor=0.5;
     1114        while (maximumRHS>=2.0) {
     1115          maximumRHS*=factor;
     1116          scale*=factor;
     1117        } /* endwhile */
     1118      }
     1119      unscale=diagonalScaleFactor_/scale;
     1120    } else {
     1121      //effectively zero
     1122      scale=0.0;
     1123      unscale=0.0;
     1124    }
     1125    //printf("--putting scales to 1.0\n");
     1126    //scale=1.0;
     1127    //unscale=1.0;
     1128    multiplyAdd(NULL,numberRows_,0.0,tempRegion,scale);
     1129    cholesky_->solve(tempRegion);
     1130    multiplyAdd(NULL,numberRows_,0.0,tempRegion,unscale);
     1131    //CoinZeroN(newError,numberRows_);
     1132    multiplyAdd(tempRegion,numberRows_,-1.0,weights_+numberColumns_,0.0);
     1133    CoinZeroN(weights_,numberColumns_);
     1134    matrix_->transposeTimes(1.0,tempRegion,weights_);
     1135    //if flagged then entries zero so can do
     1136    for (iColumn=0;iColumn<numberTotal;iColumn++)
     1137      weights_[iColumn] = weights_[iColumn]*diagonal_[iColumn]
     1138        -work2[iColumn];
     1139#else
     1140    solveSystem(weights_, tempRegion,
     1141                work2,newError,region1Save,regionSave,lastError>1.0e-5);
     1142#endif
     1143    multiplyAdd(weights_+numberColumns_,numberRows_,-1.0,newError,0.0);
     1144    matrix_->times(1.0,weights_,newError);
     1145    delete [] newError;
     1146    delete [] tempRegion;
     1147    return 0.0;
     1148  }
    9581149  double projectionTolerance=projectionTolerance_;
    9591150  //temporary
     
    16221813    numberComplementarityPairs=1;
    16231814  }
    1624   printf("gap %g - smallest %g, largest %g, pairs %d\n",
    1625         gap,smallestGap,largestGap,numberComplementarityPairs);
     1815  //printf("gap %g - smallest %g, largest %g, pairs %d\n",
     1816  // gap,smallestGap,largestGap,numberComplementarityPairs);
    16261817  return gap;
    16271818}
     
    16441835
    16451836  int iColumn;
    1646   printf("phase %d in setupForSolve, mu %g\n",phase,mu_);
     1837  //printf("phase %d in setupForSolve, mu %g\n",phase,mu_);
    16471838  switch (phase) {
    16481839  case 0:
     
    17391930    break;
    17401931  case 3:
    1741     for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
    1742       if (!flagged(iColumn)) {
    1743         double value= 0.0;
    1744         if (lowerBound(iColumn)) {
    1745           value-=work3[iColumn]/(lowerSlack[iColumn]+extra);
    1746         }
    1747         if (upperBound(iColumn)) {
    1748           value+=work4[iColumn]/(upperSlack[iColumn]+extra);
    1749         }
     1932    {
     1933      double minBeta = 0.1*mu_;
     1934      double maxBeta = 10.0*mu_;
     1935      double * work1 = deltaZ_;
     1936      double * weights = weights_;
     1937      for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     1938        work3[iColumn]=0.0;
     1939        work4[iColumn]=0.0;
     1940        if (!flagged(iColumn)) {
     1941          double value= 0.0;
     1942          if (lowerBound(iColumn)) {
     1943            double change;
     1944            if (!fakeLower(iColumn)) {
     1945              change =primal[iColumn]+weights[iColumn]-lowerSlack[iColumn]-lower[iColumn];
     1946            } else {
     1947              change =weights[iColumn];
     1948            }
     1949            double dualValue=zVec[iColumn]+actualDualStep_*work1[iColumn];
     1950            double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change;
     1951            double gapProduct=dualValue*primalValue;
     1952            if (gapProduct<minBeta) {
     1953              value-= (minBeta-gapProduct)/(lowerSlack[iColumn]+extra);
     1954            } else if (gapProduct>maxBeta) {
     1955              value-= max(maxBeta-gapProduct,-maxBeta)/(lowerSlack[iColumn]+extra);
     1956            }
     1957            value=-value;
     1958            work3[iColumn]=value;
     1959          } 
     1960          if (upperBound(iColumn)) {
     1961            double change;
     1962            if (!fakeUpper(iColumn)) {
     1963              change =upper[iColumn]-primal[iColumn]-weights[iColumn]-upperSlack[iColumn];
     1964            } else {
     1965              change =-weights[iColumn];
     1966            }
     1967            double dualValue=wVec[iColumn]+actualDualStep_*work2[iColumn];
     1968            double primalValue=upperSlack[iColumn]+actualPrimalStep_*change;
     1969            double gapProduct=dualValue*primalValue;
     1970            if (gapProduct<minBeta) {
     1971              value+= (minBeta-gapProduct)/(upperSlack[iColumn]+extra);
     1972            } else if (gapProduct>maxBeta) {
     1973              value+= max(maxBeta-gapProduct,-maxBeta)/(upperSlack[iColumn]+extra);
     1974            }
     1975            value=-value;
     1976            work4[iColumn]=value-work3[iColumn];
     1977          }
    17501978#ifndef KKT
    1751         work2[iColumn]=diagonal_[iColumn]*(dual[iColumn]+value);
     1979          work2[iColumn]=diagonal_[iColumn]*value;
    17521980#else
    1753         work2[iColumn]=dual[iColumn]+value;
     1981          work2[iColumn]=value;
    17541982#endif
    1755       } else {
    1756         work2[iColumn]=0.0;
    1757       }
    1758     }
     1983        } else {
     1984          work2[iColumn]=0.0;
     1985        }
     1986      }
     1987    }
    17591988    break;
    17601989  } /* endswitch */
  • trunk/include/ClpPredictorCorrector.hpp

    r328 r364  
    4444  //         1 corrector
    4545  //         2 primal dual
    46   double findStepLength(const int phase);
     46  double findStepLength(const int phase,const double * oldWeight);
    4747  /// findDirectionVector.
    4848  double findDirectionVector(const int phase);
Note: See TracChangeset for help on using the changeset viewer.