Changeset 163 for trunk


Ignore:
Timestamp:
Apr 16, 2003 11:20:06 AM (17 years ago)
Author:
forrest
Message:

slow progress

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpSimplexPrimalQuadratic.cpp

    r162 r163  
    1313#include "CoinIndexedVector.hpp"
    1414#include "CoinWarmStartBasis.hpp"
     15#include "CoinMpsIO.hpp"
    1516#include "ClpPrimalColumnPivot.hpp"
    1617#include "ClpMessage.hpp"
     
    108109      last[iPass][jNon]=0;
    109110  }
    110   double goodMove=false;
     111  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
     112  int goodMove=-2;
    111113  char * statusCheck = new char[numberColumns];
    112114  for (iPass=0;iPass<numberPasses;iPass++) {
     
    114116    double offset=0.0;
    115117    double objValue=-objectiveOffset;
     118    double lambda=-1.0;
     119    if (goodMove>=0) {
     120      // get best interpolation
     121      double coeff0=-objectiveOffset,coeff1=0.0,coeff2=0.0;
     122      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     123        coeff0 += saveObjective[iColumn]*solution[iColumn];
     124        coeff1 += saveObjective[iColumn]*(saveSolution[iColumn]-solution[iColumn]);
     125      }
     126      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     127        iColumn=listNonLinearColumn[jNon];
     128        double valueI = solution[iColumn];
     129        double valueISave = saveSolution[iColumn];
     130        int j;
     131        for (j=columnQuadraticStart[iColumn];
     132             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     133          int jColumn = columnQuadratic[j];
     134          double valueJ = solution[jColumn];
     135          double valueJSave = saveSolution[jColumn];
     136          double elementValue = 0.5*quadraticElement[j];
     137          coeff0 += valueI*valueJ*elementValue;
     138          coeff1 += (valueI*valueJSave+valueISave*valueJ-2.0*valueI*valueJ)*elementValue;
     139          coeff2 += (valueISave*valueJSave+valueI*valueJ-valueISave*valueJ-valueI*valueJSave)*elementValue;
     140        }
     141      }
     142      double lambdaValue;
     143      if (fabs(coeff2)<1.0e-9) {
     144        if (coeff1+coeff2>=0.0)
     145          lambda = 0.0;
     146        else
     147          lambda = 1.0;
     148      } else {
     149        lambda = -(0.5*coeff1)/coeff2;
     150        if (lambda>1.0||lambda<0.0) {
     151          if (coeff1+coeff2>=0.0)
     152            lambda = 0.0;
     153          else
     154            lambda = 1.0;
     155        }
     156      }
     157      lambdaValue = lambda*lambda*coeff2+lambda*coeff1+coeff0;
     158      printf("coeffs %g %g %g - lastobj %g\n",coeff0,coeff1,coeff2,lastObjective);
     159      printf("obj at saved %g, obj now %g zero deriv at %g - value %g\n",
     160             coeff0+coeff1+coeff2,coeff0,lambda,lambdaValue);
     161      if (lambda>0.0&&lambda<=1.0) {
     162        // update solution
     163        for (iColumn=0;iColumn<numberColumns;iColumn++)
     164          solution[iColumn] = lambda * saveSolution[iColumn]
     165            + (1.0-lambda) * solution[iColumn];
     166        if (lambda>0.999) {
     167          memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
     168          memcpy(status_,saveStatus,numberRows+numberColumns);
     169        }
     170        if (lambda>0.99999&&fabs(coeff1+coeff2)>1.0e-2) {
     171          // tighten all
     172          goodMove=-1;
     173        }
     174      }
     175    }
    116176    memcpy(objective,saveObjective,numberColumns*sizeof(double));
    117177    for (iColumn=0;iColumn<numberColumns;iColumn++)
     
    173233      }
    174234    }
    175     if (objValue<=lastObjective)
    176       goodMove=true;
    177     else
    178       goodMove=false;
     235    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
    179236    double maxDelta=0.0;
     237    if (goodMove>=0) {
     238      if (objValue<=lastObjective)
     239        goodMove=1;
     240      else
     241        goodMove=0;
     242    } else {
     243      maxDelta=1.0e10;
     244    }
    180245    double maxGap=0.0;
    181246    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     
    183248      maxDelta = max(maxDelta,
    184249                     fabs(solution[iColumn]-saveSolution[iColumn]));
    185       if (goodMove) {
     250      if (goodMove>0) {
    186251        if (last[0][jNon]*last[1][jNon]<0) {
    187252          // halve
     
    192257            trust[jNon] *= 1.5;
    193258        }
    194       } else if (trust[jNon]>10.0*deltaTolerance) {
     259      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
    195260        trust[jNon] *= 0.5;
    196261      }
     
    198263    }
    199264    std::cout<<"largest gap is "<<maxGap<<std::endl;
    200     if (goodMove) {
     265    if (goodMove>0) {
    201266      double drop = lastObjective-objValue;
    202267      std::cout<<"True drop was "<<drop<<std::endl;
    203268      std::cout<<"largest delta is "<<maxDelta<<std::endl;
    204       if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove) {
     269      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&lambda<0.99999) {
    205270        std::cout<<"Exiting"<<std::endl;
    206271        break;
    207272      }
    208     } else {
    209       //lastObjective += 1.0e-3; // to stop exit
    210     }
     273    }
     274    if (!iPass)
     275      goodMove=1;
     276    targetDrop=0.0;
     277    double * r = this->dualColumnSolution();
    211278    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    212279      iColumn=listNonLinearColumn[jNon];
     
    218285                               trueUpper[jNon]);
    219286    }
     287    if (iPass) {
     288      // get reduced costs
     289      this->matrix()->transposeTimes(savePi,
     290                                     this->dualColumnSolution());
     291      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     292        iColumn=listNonLinearColumn[jNon];
     293        double dj = objective[iColumn]-r[iColumn];
     294        r[iColumn]=dj;
     295        if (dj<0.0)
     296          targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
     297        else
     298          targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
     299      }
     300    } else {
     301      memset(r,0,numberColumns*sizeof(double));
     302    }
     303    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
     304      iColumn=listNonLinearColumn[jNon];
     305      if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
     306        columnLower[iColumn]=max(solution[iColumn],
     307                                 trueLower[jNon]);
     308        columnUpper[iColumn]=min(solution[iColumn]
     309                                 +trust[jNon],
     310                                 trueUpper[jNon]);
     311      } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
     312        columnLower[iColumn]=max(solution[iColumn]
     313                                 -trust[jNon],
     314                                 trueLower[jNon]);
     315        columnUpper[iColumn]=min(solution[iColumn],
     316                                 trueUpper[jNon]);
     317      } else {
     318        columnLower[iColumn]=max(solution[iColumn]
     319                                 -trust[jNon],
     320                                 trueLower[jNon]);
     321        columnUpper[iColumn]=min(solution[iColumn]
     322                                 +trust[jNon],
     323                                 trueUpper[jNon]);
     324      }
     325    }
    220326    if (goodMove) {
    221327      memcpy(saveSolution,solution,numberColumns*sizeof(double));
     
    223329      memcpy(saveStatus,status_,numberRows+numberColumns);
    224330     
    225       targetDrop=0.0;
    226       if (iPass) {
    227         // get reduced costs
    228         this->matrix()->transposeTimes(savePi,
    229                                      this->dualColumnSolution());
    230         double * r = this->dualColumnSolution();
    231         for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    232           iColumn=listNonLinearColumn[jNon];
    233           double dj = objective[iColumn]-r[iColumn];
    234           r[iColumn]=dj;
    235           if (dj<0.0)
    236             targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
    237           else
    238             targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
    239         }
    240       }
    241331      std::cout<<"Pass - "<<iPass
    242332               <<", target drop is "<<targetDrop
     
    251341        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    252342          iColumn=listNonLinearColumn[jNon];
    253           printf("Trust %d %g - solution %d %g obj %g dj %g state %c\n",
     343          printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
    254344                 jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
    255                  r[iColumn],statusCheck[iColumn]);
     345                 r[iColumn],statusCheck[iColumn],columnLower[iColumn],
     346                 columnUpper[iColumn]);
    256347        }
    257348      }
    258349      setLogLevel(63);
     350      this->scaling(false);
    259351      this->primal(1);
     352      goodMove=1;
    260353    } else {
    261354      // bad pass - restore solution
     
    265358      memcpy(status_,saveStatus,numberRows+numberColumns);
    266359      iPass--;
     360      goodMove=-1;
    267361    }
    268362  }
     
    307401  setDblParam(ClpObjOffset,objectiveOffset-offset);
    308402  this->primal(1);
     403  // redo values
     404  setDblParam(ClpObjOffset,objectiveOffset);
     405  objectiveValue_ += optimizationDirection_*offset;
    309406  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
    310407    iColumn=listNonLinearColumn[jNon];
     
    312409    columnUpper[iColumn]= trueUpper[jNon];
    313410  }
     411  memcpy(objective,saveObjective,numberColumns*sizeof(double));
    314412  delete [] saveSolution;
    315413  for (iPass=0;iPass<3;iPass++)
     
    375473    +2*quadratic->getNumElements()
    376474    + 2*numberColumns;
     475  // type array
     476  // >=0 points to dj or x (whichever this isn't)
     477  // -2 row variable
     478  // -1 artificial (see previous)
     479  int * type = new int[newNumberColumns];
    377480  double * elements2 = new double[numberElements];
    378481  int * start2 = new int[newNumberColumns+1];
     
    451554  for (iRow=0;iRow<numberRows;iRow++) {
    452555    // should look at rows to get correct bounds
    453     columnLower2[newNumberColumns]=-COIN_DBL_MAX;
    454     columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     556    if (rowLower[iRow]==rowUpper[iRow]) {
     557      columnLower2[newNumberColumns]=-COIN_DBL_MAX;
     558      columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     559    } else if (rowLower[iRow]<-1.0e20) {
     560      assert(rowUpper[iRow]<1.0e20);
     561      columnLower2[newNumberColumns]=-COIN_DBL_MAX;
     562      columnUpper2[newNumberColumns]=0.0;
     563    } else if (rowUpper[iRow]>1.0e20) {
     564      columnLower2[newNumberColumns]=0.0;
     565      columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     566    } else {
     567      // can't do ranges just now
     568      abort();
     569    }
    455570    int j;
    456571    for (j=rowStart[iRow];
     
    460575      row2[numberElements++]=column[j]+numberRows;
    461576    }
     577    type[newNumberColumns]=-2;
    462578    newNumberColumns++;
    463579    start2[newNumberColumns]=numberElements;
     
    467583  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    468584    // dj
    469     columnLower2[newNumberColumns]=0.0;
     585    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
    470586    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
    471587    elements2[numberElements]=1.0;
    472588    row2[numberElements++]=iColumn+numberRows;
     589    type[newNumberColumns]=iColumn;
     590    type[iColumn]=newNumberColumns;
    473591    newNumberColumns++;
    474592    start2[newNumberColumns]=numberElements;
    475593    // artificial (assuming no bounds)
    476     if (getStatus(iColumn)==basic||solution[iColumn]>1.0e-7) {
     594    if (getStatus(iColumn)==basic||
     595        (solution[iColumn]>columnLower[iColumn]+1.0e-7&&
     596         solution[iColumn]<columnUpper[iColumn]-1.0e-7)) {
     597      columnUpper2[newNumberColumns-1]=0.0; // fix for now
     598      columnLower2[newNumberColumns-1]=0.0; // fix for now
    477599      if (fabs(dj[iColumn])>tolerance) {
    478600        columnLower2[newNumberColumns]=0.0;
     
    487609        start2[newNumberColumns]=numberElements;
    488610      }
    489     } else if (dj[iColumn]<-tolerance) {
    490       columnLower2[newNumberColumns]=0.0;
    491       columnUpper2[newNumberColumns]=COIN_DBL_MAX;
    492       objective2[newNumberColumns]=1.0;
    493       elements2[numberElements]=-1.0;
    494       row2[numberElements++]=iColumn+numberRows;
    495       newNumberColumns++;
    496       start2[newNumberColumns]=numberElements;
    497     }
    498   }
    499   // Create model
     611    } else {
     612      if (solution[iColumn]<columnLower[iColumn]+1.0e-7) {
     613        columnUpper2[iColumn]=columnLower[iColumn]; // fix for now
     614        if (dj[iColumn]<-tolerance) {
     615          columnLower2[newNumberColumns]=0.0;
     616          columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     617          objective2[newNumberColumns]=1.0;
     618          elements2[numberElements]=-1.0;
     619          row2[numberElements++]=iColumn+numberRows;
     620          newNumberColumns++;
     621          start2[newNumberColumns]=numberElements;
     622        }
     623      } else {
     624        columnLower2[iColumn]=columnUpper[iColumn]; // fix for now
     625        if (dj[iColumn]>tolerance) {
     626          columnLower2[newNumberColumns]=0.0;
     627          columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     628          objective2[newNumberColumns]=1.0;
     629          elements2[numberElements]=1.0;
     630          row2[numberElements++]=iColumn+numberRows;
     631          newNumberColumns++;
     632          start2[newNumberColumns]=numberElements;
     633        }
     634      }
     635    }
     636    // temp
     637    //columnLower2[iColumn]=solution[iColumn];
     638    //columnUpper2[iColumn]=solution[iColumn];
     639  }
     640  // Create model
    500641  ClpSimplex model2;
    501642  model2.loadProblem(newNumberColumns,newNumberRows,
     
    515656  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    516657    model2.setStatus(iColumn,getStatus(iColumn));
    517     if (getStatus(iColumn)==basic) {
     658    if (getStatus(iColumn)==basic||
     659        (solution[iColumn]>columnLower[iColumn]+1.0e-7&&
     660         solution[iColumn]<columnUpper[iColumn]-1.0e-7)) {
    518661      solution2[newNumberColumns++]=0.0;
    519       if (fabs(dj[iColumn])>tolerance) {
     662      if (fabs(dj[iColumn])>tolerance)
    520663        solution2[newNumberColumns++]=fabs(dj[iColumn]);
    521       }
    522664    } else if (dj[iColumn]<-tolerance) {
    523665      solution2[newNumberColumns++]=0.0;
     
    531673               model2.primalRowSolution());
    532674  // solve
     675#if 0
     676  CoinMpsIO writer;
     677  writer.setMpsData(*model2.matrix(), COIN_DBL_MAX,
     678                    model2.getColLower(), model2.getColUpper(),
     679                    model2.getObjCoefficients(),
     680                    (const char*) 0 /*integrality*/,
     681                    model2.getRowLower(), model2.getRowUpper(),
     682                    NULL,NULL);
     683  writer.writeMps("xx.mps");
     684#endif 
     685  model2.scaling(false);
    533686  model2.primal(1);
     687  // If objective 0.0 then we can drop artificials
     688  // See if any s sub j have wrong sign and/or use djs from infeasibility objective
     689  double objectiveOffset;
     690  getDblParam(ClpObjOffset,objectiveOffset);
     691  double objValue = -objectiveOffset;
     692  for (iColumn=0;iColumn<numberColumns;iColumn++)
     693    objValue += objective[iColumn]*solution2[iColumn];
     694  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     695    double valueI = solution2[iColumn];
     696    if (fabs(valueI)>1.0e-5)
     697      assert(solution2[type[iColumn]]<1.0e-7);
     698    int j;
     699    for (j=columnQuadraticStart[iColumn];
     700         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     701      int jColumn = columnQuadratic[j];
     702      double valueJ = solution2[jColumn];
     703      double elementValue = quadraticElement[j];
     704      objValue += 0.5*valueI*valueJ*elementValue;
     705    }
     706  }
     707  printf("Objective value %g\n",objValue);
     708  for (iColumn=0;iColumn<newNumberColumns;iColumn++)
     709    printf("%d %g\n",iColumn,solution2[iColumn]);
     710  delete [] type;
    534711  return 0;
    535712}
Note: See TracChangeset for help on using the changeset viewer.