Ignore:
Timestamp:
Jun 26, 2007 3:00:48 AM (13 years ago)
Author:
forrest
Message:

moving branches/devel to trunk

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk

    • Property svn:externals
      •  

        old new  
        11MSVisualStudio   https://projects.coin-or.org/svn/MSVisualStudio/trunk/ExternalsDirs/Clp
        22BuildTools    https://projects.coin-or.org/svn/BuildTools/trunk
        3 Data/Netlib   https://projects.coin-or.org/svn/Data/trunk/Netlib
        4 Data/Sample   https://projects.coin-or.org/svn/Data/trunk/Sample
        5 CoinUtils     https://projects.coin-or.org/svn/CoinUtils/stable/2.0/CoinUtils
         3ThirdParty/Blas https://projects.coin-or.org/svn/BuildTools/ThirdParty/Blas/stable/1.0
         4ThirdParty/Lapack https://projects.coin-or.org/svn/BuildTools/ThirdParty/Lapack/stable/1.0
         5Data/Netlib   https://projects.coin-or.org/svn/Data/stable/1.0/Netlib
         6Data/Sample   https://projects.coin-or.org/svn/Data/stable/1.0/Sample
         7CoinUtils     https://projects.coin-or.org/svn/CoinUtils/trunk/CoinUtils
  • trunk/Clp/src/ClpSimplexNonlinear.cpp

    r800 r1034  
    1111#include "ClpNonLinearCost.hpp"
    1212#include "ClpLinearObjective.hpp"
     13#include "ClpConstraint.hpp"
    1314#include "ClpQuadraticObjective.hpp"
    1415#include "CoinPackedMatrix.hpp"
     
    589590  //if (goToDual)
    590591  //problemStatus_=10; // try dual
     592  // Allow matrices to be sorted etc
     593  int fake=-999; // signal sort
     594  matrix_->correctSequence(this,fake,fake);
    591595}
    592596/*
     
    27102714  double offset=0.0;
    27112715  double objValue=0.0;
     2716  int exitPass = 2*numberPasses+10;
    27122717  for (iPass=0;iPass<numberPasses;iPass++) {
     2718    exitPass--;
    27132719    // redo objective
    27142720    offset=0.0;
    27152721    objValue=-objectiveOffset;
     2722    // make sure x updated
     2723    trueObjective->newXValues();
    27162724    double theta=-1.0;
    27172725    double maxTheta=COIN_DBL_MAX;
     
    27962804    double theta2 = trueObjective->stepLength(this,saveSolution,changeRegion,maxTheta2,
    27972805                                             objValue,predictedObj,thetaObj);
    2798    
     2806    int lastMoveStatus = goodMove;
    27992807    if (goodMove>=0) {
    28002808      theta = CoinMin(theta2,maxTheta);
     
    28192827#define LOCAL
    28202828#ifdef LOCAL
     2829        bool absolutelyOptimal=false;
    28212830        int saveScaling = scalingFlag();
    28222831        scaling(0);
     
    28842893            directionVector(longArray,spare,rowArray,0,
    28852894                            normFlagged,normUnflagged,numberNonBasic);
    2886             if (normFlagged+normUnflagged<1.0e-8)
     2895            if (normFlagged+normUnflagged<1.0e-8) {
     2896              absolutelyOptimal=true;
    28872897              break;  //looks optimal
     2898            }
    28882899            double djNorm=0.0;
    28892900            int iIndex;
     
    28982909            if (!numberNonBasic) {
    28992910              // looks optimal
     2911              absolutelyOptimal=true;
    29002912              break;
    29012913            }
     
    30053017        memset(rowActivity_,0,numberRows_*sizeof(double));
    30063018        times(1.0,solution,rowActivity_);
    3007         if (theta>0.99999&&theta2<1.9) {
     3019        if (theta>0.99999&&theta2<1.9&&!absolutelyOptimal) {
    30083020          // If big changes then tighten
    30093021          /*  thetaObj is objvalue + a*2*2 +b*2
     
    30243036    memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
    30253037           numberColumns*sizeof(double));
     3038    //printf("offset comp %g orig %g\n",offset,objectiveOffset);
    30263039    // could do faster
    30273040    trueObjective->stepLength(this,solution,changeRegion,0.0,
    30283041                                             objValue,predictedObj,thetaObj);
     3042    printf("offset comp %g orig %g - obj from stepLength %g\n",offset,objectiveOffset,objValue);
    30293043    setDblParam(ClpObjOffset,objectiveOffset+offset);
    30303044    int * temp=last[2];
     
    30523066    double maxDelta=0.0;
    30533067    if (goodMove>=0) {
    3054       if (objValue<=lastObjective+1.0e-15*fabs(lastObjective))
     3068      if (objValue-lastObjective<=1.0e-15*fabs(lastObjective))
    30553069        goodMove=1;
    30563070      else
     
    30953109        trust[jNon] *=0.0001;
    30963110    }
     3111    if(lastMoveStatus==-1&&goodMove==-1)
     3112      goodMove=1; // to force solve
    30973113    if (goodMove>0) {
    30983114      double drop = lastObjective-objValue;
     
    31673183    }
    31683184#endif
    3169     if (goodMove) {
     3185    if (goodMove>0) {
    31703186      memcpy(saveSolution,solution,numberColumns*sizeof(double));
    31713187      memcpy(saveRowSolution,rowActivity_,numberRows*sizeof(double));
     
    32573273      memcpy(status_,saveStatus,numberRows+numberColumns);
    32583274      iPass--;
     3275      assert (exitPass>0);
    32593276      goodMove=-1;
    32603277    }
     
    33073324  return 0;
    33083325}
     3326/* Primal algorithm for nonlinear constraints
     3327   Using a semi-trust region approach as for pooling problem
     3328   This is in because I have it lying around
     3329   
     3330*/
     3331int
     3332ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint ** constraints,
     3333                               int numberPasses, double deltaTolerance)
     3334{
     3335  if (!numberConstraints) {
     3336    // no nonlinear constraints - may be nonlinear objective
     3337    return primalSLP(numberPasses,deltaTolerance);
     3338  }
     3339  // check all matrix for odd rows is empty
     3340  int iConstraint;
     3341  int numberBad=0;
     3342  CoinPackedMatrix * columnCopy = matrix();
     3343  // Get a row copy in standard format
     3344  CoinPackedMatrix copy;
     3345  copy.reverseOrderedCopyOf(*columnCopy);
     3346  // get matrix data pointers
     3347  //const int * column = copy.getIndices();
     3348  //const CoinBigIndex * rowStart = copy.getVectorStarts();
     3349  const int * rowLength = copy.getVectorLengths();
     3350  //const double * elementByRow = copy.getElements();
     3351  int numberArtificials=0;
     3352  // see how many extra we need
     3353  CoinBigIndex numberExtra=0;
     3354  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3355    ClpConstraint * constraint = constraints[iConstraint];
     3356    int iRow = constraint->rowNumber();
     3357    assert (iRow>=0);
     3358    int n = constraint->numberCoefficients() -rowLength[iRow];
     3359    numberExtra += n;
     3360    if (iRow>=numberRows_)
     3361      numberBad++;
     3362    else if (rowLength[iRow]&&n)
     3363      numberBad++;
     3364    if (rowLower_[iRow]>-1.0e20)
     3365      numberArtificials++;
     3366    if (rowUpper_[iRow]<1.0e20)
     3367      numberArtificials++;
     3368  }
     3369  if (numberBad)
     3370    return numberBad;
     3371  ClpObjective * trueObjective = NULL;
     3372  if (objective_->type()>=2) {
     3373    // Replace objective
     3374    trueObjective = objective_;
     3375    objective_=new ClpLinearObjective(NULL,numberColumns_);
     3376  }
     3377  ClpSimplex newModel(*this);
     3378  // we can put back true objective
     3379  if (trueObjective) {
     3380    // Replace objective
     3381    delete objective_;
     3382    objective_=trueObjective;
     3383  }
     3384  int numberColumns2 = numberColumns_;
     3385  double penalty=1.0e9;
     3386  if (numberArtificials) {
     3387    numberColumns2 += numberArtificials;
     3388    int * addStarts = new int [numberArtificials+1];
     3389    int * addRow = new int[numberArtificials];
     3390    double * addElement = new double[numberArtificials];
     3391    addStarts[0]=0;
     3392    double * addCost = new double [numberArtificials];
     3393    numberArtificials=0;
     3394    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3395      ClpConstraint * constraint = constraints[iConstraint];
     3396      int iRow = constraint->rowNumber();
     3397      if (rowLower_[iRow]>-1.0e20) {
     3398        addRow[numberArtificials]=iRow;
     3399        addElement[numberArtificials]=1.0;
     3400        addCost[numberArtificials]=penalty;
     3401        numberArtificials++;
     3402        addStarts[numberArtificials]=numberArtificials;
     3403      }
     3404      if (rowUpper_[iRow]<1.0e20) {
     3405        addRow[numberArtificials]=iRow;
     3406        addElement[numberArtificials]=-1.0;
     3407        addCost[numberArtificials]=penalty;
     3408        numberArtificials++;
     3409        addStarts[numberArtificials]=numberArtificials;
     3410      }
     3411    }
     3412    newModel.addColumns(numberArtificials,NULL,NULL,addCost,
     3413                       addStarts,addRow,addElement);
     3414    delete [] addStarts;
     3415    delete [] addRow;
     3416    delete [] addElement;
     3417    delete [] addCost;
     3418  }
     3419  // find nonlinear columns
     3420  int * listNonLinearColumn = new int [numberColumns_];
     3421  char * mark = new char[numberColumns_];
     3422  memset(mark,0,numberColumns_);
     3423  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3424    ClpConstraint * constraint = constraints[iConstraint];
     3425    constraint->markNonlinear(mark);
     3426  }
     3427  if (trueObjective)
     3428    trueObjective->markNonlinear(mark);
     3429  int iColumn;
     3430  int numberNonLinearColumns=0;
     3431  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     3432    if (mark[iColumn])
     3433      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
     3434  }
     3435  // build row copy of matrix with space for nonzeros
     3436  // Get column copy
     3437  columnCopy = newModel.matrix();
     3438  copy.reverseOrderedCopyOf(*columnCopy);
     3439  // get matrix data pointers
     3440  const int * column = copy.getIndices();
     3441  const CoinBigIndex * rowStart = copy.getVectorStarts();
     3442  rowLength = copy.getVectorLengths();
     3443  const double * elementByRow = copy.getElements();
     3444  numberExtra +=copy.getNumElements();
     3445  CoinBigIndex * newStarts = new CoinBigIndex [numberRows_+1];
     3446  int * newColumn = new int[numberExtra];
     3447  double * newElement = new double[numberExtra];
     3448  newStarts[0]=0;
     3449  int * backRow = new int [numberRows_];
     3450  int iRow;
     3451  for (iRow=0;iRow<numberRows_;iRow++)
     3452    backRow[iRow]=-1;
     3453  for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3454    ClpConstraint * constraint = constraints[iConstraint];
     3455    iRow = constraint->rowNumber();
     3456    backRow[iRow]=iConstraint;
     3457  }
     3458  numberExtra=0;
     3459  for (iRow=0;iRow<numberRows_;iRow++) {
     3460    if (backRow[iRow]<0) {
     3461      // copy normal
     3462      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
     3463           j++) {
     3464        newColumn[numberExtra]=column[j];
     3465        newElement[numberExtra++] = elementByRow[j];
     3466      }
     3467    } else {
     3468      ClpConstraint * constraint = constraints[backRow[iRow]];
     3469      assert(iRow == constraint->rowNumber());
     3470      int numberArtificials=0;
     3471      if (rowLower_[iRow]>-1.0e20)
     3472        numberArtificials++;
     3473      if (rowUpper_[iRow]<1.0e20)
     3474        numberArtificials++;
     3475      if (numberArtificials==rowLength[iRow]) {
     3476        // all possible
     3477        memset(mark,0,numberColumns_);
     3478        constraint->markNonzero(mark);
     3479        for (int k=0;k<numberColumns_;k++) {
     3480          if (mark[k]) {
     3481            newColumn[numberExtra]=k;
     3482            newElement[numberExtra++] = 1.0;
     3483          }
     3484        }
     3485        // copy artificials
     3486        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
     3487             j++) {
     3488          newColumn[numberExtra]=column[j];
     3489          newElement[numberExtra++] = elementByRow[j];
     3490        }
     3491      } else {
     3492        // there already
     3493        // copy
     3494        for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];
     3495             j++) {
     3496          newColumn[numberExtra]=column[j];
     3497          assert (elementByRow[j]);
     3498          newElement[numberExtra++] = elementByRow[j];
     3499        }
     3500      }
     3501    }
     3502    newStarts[iRow+1]=numberExtra;
     3503  }
     3504  delete [] backRow;
     3505  CoinPackedMatrix saveMatrix(false,numberColumns2,numberRows_,
     3506                              numberExtra,newElement,newColumn,newStarts,NULL,0.0,0.0);
     3507  delete [] newStarts;
     3508  delete [] newColumn;
     3509  delete [] newElement;
     3510  delete [] mark;
     3511  // get feasible
     3512  newModel.primal(1);
     3513  // still infeasible
     3514  if (newModel.numberPrimalInfeasibilities()) {
     3515    delete [] listNonLinearColumn;
     3516    return 0;
     3517  }
     3518  int numberRows = newModel.numberRows();
     3519  double * columnLower = newModel.columnLower();
     3520  double * columnUpper = newModel.columnUpper();
     3521  double * objective = newModel.objective();
     3522  double * rowLower = newModel.rowLower();
     3523  double * rowUpper = newModel.rowUpper();
     3524  double * solution = newModel.primalColumnSolution();
     3525  int jNon;
     3526  int * last[3];
     3527 
     3528  double * trust = new double[numberNonLinearColumns];
     3529  double * trueLower = new double[numberNonLinearColumns];
     3530  double * trueUpper = new double[numberNonLinearColumns];
     3531  double objectiveOffset;
     3532  double objectiveOffset2;
     3533  getDblParam(ClpObjOffset,objectiveOffset);
     3534  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3535    iColumn=listNonLinearColumn[jNon];
     3536    double upper = columnUpper[iColumn];
     3537    double lower = columnLower[iColumn];
     3538    if (solution[iColumn]<trueLower[jNon])
     3539      solution[iColumn]=trueLower[jNon];
     3540    else if (solution[iColumn]>trueUpper[jNon])
     3541      solution[iColumn]=trueUpper[jNon];
     3542#if 0
     3543    double large = CoinMax(1000.0,10.0*fabs(solution[iColumn]));
     3544    if (upper>1.0e10)
     3545      upper = solution[iColumn]+large;
     3546    if (lower<-1.0e10)
     3547      lower = solution[iColumn]-large;
     3548#else
     3549    upper = solution[iColumn]+0.5;
     3550    lower = solution[iColumn]-0.5;
     3551#endif
     3552    //columnUpper[iColumn]=upper;
     3553    trust[jNon]=0.05*(1.0+upper-lower);
     3554    trueLower[jNon]=columnLower[iColumn];
     3555    //trueUpper[jNon]=upper;
     3556    trueUpper[jNon]=columnUpper[iColumn];
     3557  }
     3558  bool tryFix=false;
     3559  int iPass;
     3560  double lastObjective=1.0e31;
     3561  double * saveSolution = new double [numberColumns2+numberRows];
     3562  char * saveStatus = new char [numberColumns2+numberRows];
     3563  double targetDrop=1.0e31;
     3564  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
     3565  for (iPass=0;iPass<3;iPass++) {
     3566    last[iPass]=new int[numberNonLinearColumns];
     3567    for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     3568      last[iPass][jNon]=0;
     3569  }
     3570  int numberZeroPasses=0;
     3571  bool zeroTargetDrop=false;
     3572  double * gradient = new double [numberColumns_];
     3573#define SMALL_FIX 0.0
     3574  for (iPass=0;iPass<numberPasses;iPass++) {
     3575    objectiveOffset2 = objectiveOffset;
     3576    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3577      iColumn=listNonLinearColumn[jNon];
     3578      if (solution[iColumn]<trueLower[jNon])
     3579        solution[iColumn]=trueLower[jNon];
     3580      else if (solution[iColumn]>trueUpper[jNon])
     3581        solution[iColumn]=trueUpper[jNon];
     3582      columnLower[iColumn]=CoinMax(solution[iColumn]
     3583                                   -trust[jNon],
     3584                                   trueLower[jNon]);
     3585      if (!trueLower[jNon]&&columnLower[iColumn]<SMALL_FIX)
     3586        columnLower[iColumn]=SMALL_FIX;
     3587      columnUpper[iColumn]=CoinMin(solution[iColumn]
     3588                                   +trust[jNon],
     3589                                   trueUpper[jNon]);
     3590      if (!trueLower[jNon])
     3591        columnUpper[iColumn] = CoinMax(columnUpper[iColumn],
     3592                                       columnLower[iColumn]+SMALL_FIX);
     3593      if (!trueLower[jNon]&&tryFix&&
     3594          columnLower[iColumn]==SMALL_FIX&&
     3595          columnUpper[iColumn]<3.0*SMALL_FIX) {
     3596        columnLower[iColumn]=0.0;
     3597        solution[iColumn]=0.0;
     3598        columnUpper[iColumn]=0.0;
     3599        printf("fixing %d\n",iColumn);
     3600      }
     3601    }
     3602    // redo matrix
     3603    double offset;
     3604    CoinPackedMatrix newMatrix(saveMatrix);
     3605    // get matrix data pointers
     3606    column = newMatrix.getIndices();
     3607    rowStart = newMatrix.getVectorStarts();
     3608    rowLength = newMatrix.getVectorLengths();
     3609    // make sure x updated
     3610    if (numberConstraints)
     3611      constraints[0]->newXValues();
     3612    else
     3613      trueObjective->newXValues();
     3614    double * changeableElement = newMatrix.getMutableElements();
     3615    if (trueObjective) {
     3616      memcpy(objective,trueObjective->gradient(this,solution,offset,true,2),
     3617             numberColumns_*sizeof(double));
     3618    } else {
     3619      memcpy(objective,objective_->gradient(this,solution,offset,true,2),
     3620             numberColumns_*sizeof(double));
     3621    }
     3622    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3623      ClpConstraint * constraint = constraints[iConstraint];
     3624      int iRow = constraint->rowNumber();
     3625      double functionValue;
     3626#ifndef NDEBUG
     3627      int numberErrors =
     3628#endif
     3629      constraint->gradient(&newModel,solution,gradient,functionValue,offset);
     3630      assert (!numberErrors);
     3631      double dualValue = newModel.dualRowSolution()[iRow];
     3632      int numberCoefficients = constraint->numberCoefficients();
     3633      for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+numberCoefficients;j++) {
     3634        int iColumn = column[j];
     3635        changeableElement[j] = gradient[iColumn];
     3636        //objective[iColumn] -= dualValue*gradient[iColumn];
     3637        gradient[iColumn]=0.0;
     3638      }
     3639      for (int k=0;k<numberColumns_;k++)
     3640        assert (!gradient[k]);
     3641      if (rowLower_[iRow]>-1.0e20)
     3642        rowLower[iRow] = rowLower_[iRow] - offset;
     3643      if (rowUpper_[iRow]<1.0e20)
     3644        rowUpper[iRow] = rowUpper_[iRow] - offset;
     3645    }
     3646    // Replace matrix
     3647    // Get a column copy in standard format
     3648    CoinPackedMatrix * columnCopy = new CoinPackedMatrix();;
     3649    columnCopy->reverseOrderedCopyOf(newMatrix);
     3650    newModel.replaceMatrix(columnCopy,true);
     3651    // solve
     3652    newModel.primal(1);
     3653    if (newModel.status()==1) {
     3654      // Infeasible!
     3655      newModel.allSlackBasis();
     3656      newModel.primal();
     3657      newModel.writeMps("infeas.mps");
     3658      assert(!newModel.status());
     3659    }
     3660    double sumInfeas=0.0;
     3661    int numberInfeas=0;
     3662    for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn++) {
     3663      if (solution[iColumn]>1.0e-8) {
     3664        numberInfeas++;
     3665        sumInfeas += solution[iColumn];
     3666      }
     3667    }
     3668    printf("%d artificial infeasibilities - summing to %g\n",
     3669           numberInfeas,sumInfeas);
     3670    double infValue=0.0;
     3671    double objValue=0.0;
     3672    // make sure x updated
     3673    if (numberConstraints)
     3674      constraints[0]->newXValues();
     3675    else
     3676      trueObjective->newXValues();
     3677    if (trueObjective) {
     3678      objValue=trueObjective->objectiveValue(this,solution);
     3679      printf("objective offset %g\n",offset);
     3680      objectiveOffset2 =objectiveOffset+offset; // ? sign
     3681      newModel.setObjectiveOffset(objectiveOffset2);
     3682    } else {
     3683      objValue=objective_->objectiveValue(this,solution);
     3684    }
     3685    for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) {
     3686      ClpConstraint * constraint = constraints[iConstraint];
     3687      int iRow = constraint->rowNumber();
     3688      double functionValue=constraint->functionValue(this,solution);
     3689      double dualValue = newModel.dualRowSolution()[iRow];
     3690      if (numberConstraints<50)
     3691        printf("For row %d current value is %g (activity %g) , dual is %g - offset %g\n",iRow,functionValue,
     3692               newModel.primalRowSolution()[iRow],
     3693               dualValue,offset);
     3694      if (functionValue<rowLower_[iRow]-1.0e-5) {
     3695        double under = rowLower_[iRow]-functionValue;
     3696        infValue += under;
     3697      } else if (functionValue>rowUpper_[iRow]+1.0e-5) {
     3698        double over = functionValue-rowUpper_[iRow];
     3699        infValue += over;
     3700      }
     3701    }
     3702    if (infValue)
     3703      printf("Sum infeasibilities %g ",infValue);
     3704    if (infValue<0.1)
     3705      infValue=0.0;
     3706    infValue *= penalty;
     3707    if (infValue)
     3708      printf("Infeasible obj %g ",infValue);
     3709    if (objectiveOffset2)
     3710      printf("offset2 %g\n",objectiveOffset2);
     3711    objValue -= objectiveOffset2;
     3712    printf("True objective %g\n",objValue);
     3713    objValue += infValue;
     3714    if (iPass) {
     3715      double drop = lastObjective-objValue;
     3716      std::cout<<"True drop was "<<drop<<std::endl;
     3717      if (drop<-0.05*fabs(objValue)-1.0e-4) {
     3718        // pretty bad - go back and halve
     3719        memcpy(solution,saveSolution,numberColumns2*sizeof(double));
     3720        memcpy(newModel.primalRowSolution(),saveSolution+numberColumns2,
     3721               numberRows*sizeof(double));
     3722        memcpy(newModel.statusArray(),saveStatus,
     3723               numberColumns2+numberRows);
     3724        for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     3725          if (trust[jNon]>0.1)
     3726            trust[jNon] *= 0.5;
     3727          else
     3728            trust[jNon] *= 0.9;
     3729       
     3730        printf("** Halving trust\n");
     3731        objValue=lastObjective;
     3732        continue;
     3733      } else if ((iPass%25)==-1) {
     3734        for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     3735          trust[jNon] *= 2.0;
     3736        printf("** Doubling trust\n");
     3737      }
     3738      int * temp=last[2];
     3739      last[2]=last[1];
     3740      last[1]=last[0];
     3741      last[0]=temp;
     3742      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3743        iColumn=listNonLinearColumn[jNon];
     3744        double change = solution[iColumn]-saveSolution[iColumn];
     3745        if (change<-1.0e-5) {
     3746          if (fabs(change+trust[jNon])<1.0e-5)
     3747            temp[jNon]=-1;
     3748          else
     3749            temp[jNon]=-2;
     3750        } else if(change>1.0e-5) {
     3751          if (fabs(change-trust[jNon])<1.0e-5)
     3752            temp[jNon]=1;
     3753          else
     3754            temp[jNon]=2;
     3755        } else {
     3756          temp[jNon]=0;
     3757        }
     3758      }
     3759      double maxDelta=0.0;
     3760      double smallestTrust=1.0e31;
     3761      double smallestNonLinearGap=1.0e31;
     3762      double smallestGap=1.0e31;
     3763      bool increasing=false;
     3764      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     3765        double gap = columnUpper[iColumn]-columnLower[iColumn];
     3766        assert (gap>=0.0);
     3767        if (gap)
     3768          smallestGap = CoinMin(smallestGap,gap);
     3769      }
     3770      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3771        iColumn=listNonLinearColumn[jNon];
     3772        double gap = columnUpper[iColumn]-columnLower[iColumn];
     3773        assert (gap>=0.0);
     3774        if (gap) {
     3775          smallestNonLinearGap = CoinMin(smallestNonLinearGap,gap);
     3776          if (gap<1.0e-7&&iPass==1) {
     3777            printf("Small gap %d %d %g %g %g\n",
     3778                   jNon,iColumn,columnLower[iColumn],columnUpper[iColumn],
     3779                   gap);
     3780            //trueUpper[jNon]=trueLower[jNon];
     3781            //columnUpper[iColumn]=columnLower[iColumn];
     3782          }
     3783        }
     3784        maxDelta = CoinMax(maxDelta,
     3785                       fabs(solution[iColumn]-saveSolution[iColumn]));
     3786        if (last[0][jNon]*last[1][jNon]<0) {
     3787          // halve
     3788          if (trust[jNon]>1.0)
     3789            trust[jNon] *= 0.5;
     3790          else
     3791            trust[jNon] *= 0.7;
     3792        } else {
     3793          // ? only increase if +=1 ?
     3794          if (last[0][jNon]==last[1][jNon]&&
     3795              last[0][jNon]==last[2][jNon]&&
     3796              last[0][jNon]) {
     3797            trust[jNon] *= 1.8;
     3798            increasing=true;
     3799          }
     3800        }
     3801        smallestTrust = CoinMin(smallestTrust,trust[jNon]);
     3802      }
     3803      std::cout<<"largest delta is "<<maxDelta
     3804               <<", smallest trust is "<<smallestTrust
     3805               <<", smallest gap is "<<smallestGap
     3806               <<", smallest nonlinear gap is "<<smallestNonLinearGap
     3807               <<std::endl;
     3808      if (maxDelta<deltaTolerance&&!increasing&&iPass>100) {
     3809        numberZeroPasses++;
     3810        if (numberZeroPasses==3) {
     3811          if (tryFix) {
     3812            std::cout<<"Exiting"<<std::endl;
     3813            break;
     3814          } else {
     3815            tryFix=true;
     3816            for (jNon=0;jNon<numberNonLinearColumns;jNon++)
     3817              trust[jNon] = CoinMax(trust[jNon],1.0e-1);
     3818            numberZeroPasses=0;
     3819          }
     3820        }
     3821      } else {
     3822        numberZeroPasses=0;
     3823      }
     3824    }
     3825    memcpy(saveSolution,solution,numberColumns2*sizeof(double));
     3826    memcpy(saveSolution+numberColumns2,newModel.primalRowSolution(),
     3827           numberRows*sizeof(double));
     3828    memcpy(saveStatus,newModel.statusArray(),
     3829           numberColumns2+numberRows);
     3830   
     3831    targetDrop=0.0;
     3832    if (iPass) {
     3833      // get reduced costs
     3834      const double * pi = newModel.dualRowSolution();
     3835      newModel.matrix()->transposeTimes(pi,
     3836                                        newModel.dualColumnSolution());
     3837      double * r = newModel.dualColumnSolution();
     3838      for (iColumn=0;iColumn<numberColumns_;iColumn++)
     3839        r[iColumn] = objective[iColumn]-r[iColumn];
     3840      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3841        iColumn=listNonLinearColumn[jNon];
     3842        double dj = r[iColumn];
     3843        if (dj<-1.0e-6) {
     3844          double drop = -dj*(columnUpper[iColumn]-solution[iColumn]);
     3845          //double upper = CoinMin(trueUpper[jNon],solution[iColumn]+0.1);
     3846          //double drop2 = -dj*(upper-solution[iColumn]);
     3847#if 0
     3848          if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2&&iPass>100))
     3849            printf("Big drop %d %g %g %g %g T %g\n",
     3850                   iColumn,columnLower[iColumn],solution[iColumn],
     3851                   columnUpper[iColumn],dj,trueUpper[jNon]);
     3852#endif
     3853          targetDrop += drop;
     3854          if (dj<-1.0e-1&&trust[jNon]<1.0e-3
     3855              &&trueUpper[jNon]-solution[iColumn]>1.0e-2) {
     3856            trust[jNon] *= 1.5;
     3857            //printf("Increasing trust on %d to %g\n",
     3858            //     iColumn,trust[jNon]);
     3859          }
     3860        } else if (dj>1.0e-6) {
     3861          double drop = -dj*(columnLower[iColumn]-solution[iColumn]);
     3862          //double lower = CoinMax(trueLower[jNon],solution[iColumn]-0.1);
     3863          //double drop2 = -dj*(lower-solution[iColumn]);
     3864#if 0
     3865          if (drop>1.0e8||drop2>100.0*drop||(drop>1.0e-2))
     3866            printf("Big drop %d %g %g %g %g T %g\n",
     3867                   iColumn,columnLower[iColumn],solution[iColumn],
     3868                   columnUpper[iColumn],dj,trueLower[jNon]);
     3869#endif
     3870          targetDrop += drop;
     3871          if (dj>1.0e-1&&trust[jNon]<1.0e-3
     3872              &&solution[iColumn]-trueLower[jNon]>1.0e-2) {
     3873            trust[jNon] *= 1.5;
     3874            printf("Increasing trust on %d to %g\n",
     3875                   iColumn,trust[jNon]);
     3876          }
     3877        }
     3878      }
     3879    }
     3880    std::cout<<"Pass - "<<iPass
     3881             <<", target drop is "<<targetDrop
     3882             <<std::endl;
     3883    if (iPass>1&&targetDrop<1.0e-5&&zeroTargetDrop)
     3884      break;
     3885    if (iPass>1&&targetDrop<1.0e-5)
     3886      zeroTargetDrop = true;
     3887    else
     3888      zeroTargetDrop = false;
     3889    //if (iPass==5)
     3890    //newModel.setLogLevel(63);
     3891    lastObjective = objValue;
     3892    // take out when ClpPackedMatrix changed
     3893    //newModel.scaling(false);
     3894#if 0
     3895    CoinMpsIO writer;
     3896    writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX,
     3897                      newModel.getColLower(), newModel.getColUpper(),
     3898                      newModel.getObjCoefficients(),
     3899                      (const char*) 0 /*integrality*/,
     3900                      newModel.getRowLower(), newModel.getRowUpper(),
     3901                      NULL,NULL);
     3902    writer.writeMps("xx.mps");
     3903#endif
     3904  }
     3905  delete [] saveSolution;
     3906  delete [] saveStatus;
     3907  for (iPass=0;iPass<3;iPass++)
     3908    delete [] last[iPass];
     3909  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     3910    iColumn=listNonLinearColumn[jNon];
     3911    columnLower[iColumn]=trueLower[jNon];
     3912    columnUpper[iColumn]=trueUpper[jNon];
     3913  }
     3914  delete [] trust;
     3915  delete [] trueUpper;
     3916  delete [] trueLower;
     3917  // Simplest way to get true row activity ?
     3918  double * rowActivity = newModel.primalRowSolution();
     3919  for (iRow=0;iRow<numberRows;iRow++) {
     3920    double difference;
     3921    if (fabs(rowLower_[iRow])<fabs(rowUpper_[iRow]))
     3922      difference = rowLower_[iRow]-rowLower[iRow];
     3923    else
     3924      difference = rowUpper_[iRow]-rowUpper[iRow];
     3925    rowLower[iRow]=rowLower_[iRow];
     3926    rowUpper[iRow]=rowUpper_[iRow];
     3927    if (difference) {
     3928      if (numberRows<50)
     3929        printf("For row %d activity changes from %g to %g\n",
     3930               iRow,rowActivity[iRow],rowActivity[iRow]+difference);
     3931      rowActivity[iRow]+= difference;
     3932    }
     3933  }
     3934  delete [] listNonLinearColumn;
     3935  delete [] gradient;
     3936  printf("solution still in newModel - do objective etc!\n");
     3937  objectiveValue_=newModel.objectiveValue();
     3938  numberIterations_ =newModel.numberIterations();
     3939  problemStatus_ =newModel.problemStatus();
     3940  secondaryStatus_ =newModel.secondaryStatus();
     3941  memcpy(columnActivity_,newModel.primalColumnSolution(),numberColumns_*sizeof(double));
     3942  // should do status region
     3943  CoinZeroN(rowActivity_,numberRows_);
     3944  matrix_->times(1.0,columnActivity_,rowActivity_);
     3945  return 0;
     3946}
Note: See TracChangeset for help on using the changeset viewer.