Changeset 987 for branches


Ignore:
Timestamp:
Mar 8, 2007 4:59:45 PM (13 years ago)
Author:
forrest
Message:

trying to improve idiot

Location:
branches/devel/Clp/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Clp/src/Idiot.cpp

    r986 r987  
    1010#include "Idiot.hpp"
    1111#include "CoinTime.hpp"
     12#include "CoinSort.hpp"
    1213#include "CoinMessageHandler.hpp"
    1314#include "CoinHelperFunctions.hpp"
     
    6061  }
    6162}
     63// Deals with whenUsed and slacks
     64int
     65Idiot::cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd,
     66                      double * colsol, const double * lower, const double * upper,
     67                      const double * rowLower, const double * rowUpper,
     68                      const double * cost, const double * element, double fixTolerance,
     69                      double & objValue, double & infValue)
     70{
     71  int n=0;
     72  if ((strategy_&16384)==0) {
     73    for (int i=ordinaryStart;i<ordinaryEnd;i++) {
     74      if (colsol[i]>lower[i]+fixTolerance) {
     75        if (colsol[i]<upper[i]-fixTolerance) {
     76          n++;
     77        } else {
     78          colsol[i]=upper[i];
     79        }
     80        whenUsed_[i]=iteration;
     81      } else {
     82        colsol[i]=lower[i];
     83      }
     84    }
     85    return n;
     86  } else {
     87    printf("entering inf %g, obj %g\n",infValue,objValue);
     88    int nrows=model_->getNumRows();
     89    int ncols=model_->getNumCols();
     90    int * posSlack = whenUsed_+ncols;
     91    int * negSlack = posSlack+nrows;
     92    int * nextSlack = negSlack + nrows;
     93    double * rowsol = (double *) (nextSlack+ncols);
     94    memset(rowsol,0,nrows*sizeof(double));
     95#ifdef OSI_IDIOT
     96    const CoinPackedMatrix * matrix = model_->getMatrixByCol();
     97#else
     98    ClpMatrixBase * matrix = model_->clpMatrix();
     99#endif
     100    const int * row = matrix->getIndices();
     101    const CoinBigIndex * columnStart = matrix->getVectorStarts();
     102    const int * columnLength = matrix->getVectorLengths();
     103    //const double * element = matrix->getElements();
     104    int i;
     105    objValue=0.0;
     106    infValue=0.0;
     107    for ( i=0;i<ncols;i++) {
     108      if (nextSlack[i]==-1) {
     109        // not a slack
     110        if (colsol[i]>lower[i]+fixTolerance) {
     111          if (colsol[i]<upper[i]-fixTolerance) {
     112            n++;
     113            whenUsed_[i]=iteration;
     114          } else {
     115            colsol[i]=upper[i];
     116          }
     117          whenUsed_[i]=iteration;
     118        } else {
     119          colsol[i]=lower[i];
     120        }
     121        double value = colsol[i];
     122        if (value) {
     123          objValue += cost[i]*value;
     124          CoinBigIndex j;
     125          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     126            int iRow = row[j];
     127            rowsol[iRow] += value*element[j];
     128          }
     129        }
     130      }
     131    }
     132    for (i=0;i<nrows;i++) {
     133      double rowSave=rowsol[i];
     134      int iCol;
     135      iCol =posSlack[i];
     136      if (iCol>=0) {
     137        // slide all slack down
     138        double rowValue=rowsol[i];
     139        CoinBigIndex j=columnStart[iCol];
     140        rowSave += (colsol[iCol]-lower[iCol])*element[j];
     141        colsol[iCol]=lower[iCol];
     142        while (nextSlack[iCol]>=0) {
     143          iCol = nextSlack[iCol];
     144          j=columnStart[iCol];
     145          rowSave += (colsol[iCol]-lower[iCol])*element[j];
     146          colsol[iCol]=lower[iCol];
     147        }
     148        iCol =posSlack[i];
     149        while (rowValue<rowLower[i]&&iCol>=0) {
     150          // want to increase
     151          double distance = rowLower[i]-rowValue;
     152          double value = element[columnStart[iCol]];
     153          double thisCost = cost[iCol];
     154          if (distance<=value*(upper[iCol]-lower[iCol])) {
     155            // can get there
     156            double movement = distance/value;
     157            objValue += movement*thisCost;
     158            rowValue = rowLower[i];
     159            colsol[iCol] += movement;
     160          } else {
     161            // can't get there
     162            double movement = upper[iCol]-lower[iCol];
     163            objValue += movement*thisCost;
     164            rowValue += movement*value;
     165            colsol[iCol] = upper[iCol];
     166            iCol = nextSlack[iCol];
     167          }
     168        }
     169        if (iCol>=0) {
     170          // may want to carry on - because of cost?
     171          while (cost[iCol]<0&&rowValue<rowUpper[i]) {
     172            // want to increase
     173            double distance = rowUpper[i]-rowValue;
     174            double value = element[columnStart[iCol]];
     175            double thisCost = cost[iCol];
     176            if (distance<=value*(upper[iCol]-colsol[iCol])) {
     177              // can get there
     178              double movement = distance/value;
     179              objValue += movement*thisCost;
     180              rowValue = rowUpper[i];
     181              colsol[iCol] += movement;
     182              iCol=-1;
     183            } else {
     184              // can't get there
     185              double movement = upper[iCol]-colsol[iCol];
     186              objValue += movement*thisCost;
     187              rowValue += movement*value;
     188              colsol[iCol] = upper[iCol];
     189              iCol = nextSlack[iCol];
     190            }
     191          }
     192          if (iCol>=0&&colsol[iCol]>lower[iCol]+fixTolerance&&
     193              colsol[iCol]<upper[iCol]-fixTolerance) {
     194            whenUsed_[i]=iteration;
     195            n++;
     196          }
     197        }
     198        rowsol[i]=rowValue;
     199      }
     200      iCol =negSlack[i];
     201      if (iCol>=0) {
     202        // slide all slack down
     203        double rowValue=rowsol[i];
     204        CoinBigIndex j=columnStart[iCol];
     205        rowSave += (colsol[iCol]-lower[iCol])*element[j];
     206        colsol[iCol]=lower[iCol];
     207        while (nextSlack[iCol]>=0) {
     208          iCol = nextSlack[iCol];
     209          j=columnStart[iCol];
     210          rowSave += (colsol[iCol]-lower[iCol])*element[j];
     211          colsol[iCol]=lower[iCol];
     212        }
     213        iCol =negSlack[i];
     214        while (rowValue>rowUpper[i]&&iCol>=0) {
     215          // want to increase
     216          double distance = -(rowUpper[i]-rowValue);
     217          double value = -element[columnStart[iCol]];
     218          double thisCost = cost[iCol];
     219          if (distance<=value*(upper[iCol]-lower[iCol])) {
     220            // can get there
     221            double movement = distance/value;
     222            objValue += movement*thisCost;
     223            rowValue = rowUpper[i];
     224            colsol[iCol] += movement;
     225          } else {
     226            // can't get there
     227            double movement = upper[iCol]-lower[iCol];
     228            objValue += movement*thisCost;
     229            rowValue -= movement*value;
     230            colsol[iCol] = upper[iCol];
     231            iCol = nextSlack[iCol];
     232          }
     233        }
     234        if (iCol>=0) {
     235          // may want to carry on - because of cost?
     236          while (cost[iCol]<0&&rowValue>rowLower[i]) {
     237            // want to increase
     238            double distance = -(rowLower[i]-rowValue);
     239            double value = -element[columnStart[iCol]];
     240            double thisCost = cost[iCol];
     241            if (distance<=value*(upper[iCol]-colsol[iCol])) {
     242              // can get there
     243              double movement = distance/value;
     244              objValue += movement*thisCost;
     245              rowValue = rowLower[i];
     246              colsol[iCol] += movement;
     247              iCol=-1;
     248            } else {
     249              // can't get there
     250              double movement = upper[iCol]-colsol[iCol];
     251              objValue += movement*thisCost;
     252              rowValue -= movement*value;
     253              colsol[iCol] = upper[iCol];
     254              iCol = nextSlack[iCol];
     255            }
     256          }
     257          if (iCol>=0&&colsol[iCol]>lower[iCol]+fixTolerance&&
     258              colsol[iCol]<upper[iCol]-fixTolerance) {
     259            whenUsed_[i]=iteration;
     260            n++;
     261          }
     262        }
     263        rowsol[i]=rowValue;
     264      }
     265      infValue += CoinMax(CoinMax(0.0,rowLower[i]-rowsol[i]),rowsol[i]-rowUpper[i]);
     266      // just change
     267      rowsol[i] -= rowSave;
     268    }
     269    return n;
     270  }
     271}
    62272
    63 /* returns -1 or start of costed slacks */
     273/* returns -1 if none or start of costed slacks or -2 if
     274   there are costed slacks but it is messy */
    64275 static int countCostedSlacks(OsiSolverInterface * model)
    65276{
     
    250461  dj=new double[ncols];
    251462  delete [] whenUsed_;
    252   whenUsed=whenUsed_=new int[ncols];
     463  // See if any costed slacks
     464  int numberSlacks=0;
     465  for (i=0;i<ncols;i++) {
     466    if (columnLength[i]==1)
     467      numberSlacks++;
     468  }
     469  if (!numberSlacks||true) {
     470    whenUsed_=new int[ncols];
     471  } else {
     472    printf("%d slacks\n",numberSlacks);
     473    strategy_ |= 16384;
     474    int extra = (int) (nrows*sizeof(double)/sizeof(int));
     475    whenUsed_=new int[2*ncols+2*nrows+extra];
     476    int * posSlack = whenUsed_+ncols;
     477    int * negSlack = posSlack+nrows;
     478    int * nextSlack = negSlack + nrows;
     479    for (i=0;i<nrows;i++) {
     480      posSlack[i]=-1;
     481      negSlack[i]=-1;
     482    }
     483    for (i=0;i<ncols;i++)
     484      nextSlack[i]=-1;
     485    for (i=0;i<ncols;i++) {
     486      if (columnLength[i]==1) {
     487        CoinBigIndex j=columnStart[i];
     488        int iRow = row[j];
     489        if (element[j]>0.0) {
     490          if (posSlack[iRow]==-1) {
     491            posSlack[iRow]=i;
     492          } else {
     493            int iCol = posSlack[iRow];
     494            while (nextSlack[iCol]>=0)
     495              iCol = nextSlack[iCol];
     496            nextSlack[iCol]=i;
     497          }
     498        } else {
     499          if (negSlack[iRow]==-1) {
     500            negSlack[iRow]=i;
     501          } else {
     502            int iCol = negSlack[iRow];
     503            while (nextSlack[iCol]>=0)
     504              iCol = nextSlack[iCol];
     505            nextSlack[iCol]=i;
     506          }
     507        }
     508      }
     509    }
     510    // now sort
     511    for (i=0;i<nrows;i++) {
     512      int iCol;
     513      iCol = posSlack[i];
     514      if (iCol>=0) {
     515        CoinBigIndex j=columnStart[iCol];
     516        int iRow = row[j];
     517        assert (element[j]>0.0);
     518        assert (iRow==i);
     519        dj[0]= cost[iCol]/element[j];
     520        whenUsed_[0]=iCol;
     521        int n=1;
     522        while (nextSlack[iCol]>=0) {
     523          iCol = nextSlack[iCol];
     524          CoinBigIndex j=columnStart[iCol];
     525          int iRow = row[j];
     526          assert (element[j]>0.0);
     527          assert (iRow==i);
     528          dj[n]= cost[iCol]/element[j];
     529          whenUsed_[n++]=iCol;
     530        }
     531        for (j=0;j<n;j++) {
     532          int jCol = whenUsed_[j];
     533          nextSlack[jCol]=-2;
     534        }
     535        CoinSort_2(dj,dj+n,whenUsed_);
     536        // put back
     537        iCol = whenUsed_[0];
     538        posSlack[i]=iCol;
     539        for (j=1;j<n;j++) {
     540          int jCol = whenUsed_[j];
     541          nextSlack[iCol]=jCol;
     542          iCol=jCol;
     543        }
     544      }
     545      iCol = negSlack[i];
     546      if (iCol>=0) {
     547        CoinBigIndex j=columnStart[iCol];
     548        int iRow = row[j];
     549        assert (element[j]<0.0);
     550        assert (iRow==i);
     551        dj[0]= -cost[iCol]/element[j];
     552        whenUsed_[0]=iCol;
     553        int n=1;
     554        while (nextSlack[iCol]>=0) {
     555          iCol = nextSlack[iCol];
     556          CoinBigIndex j=columnStart[iCol];
     557          int iRow = row[j];
     558          assert (element[j]<0.0);
     559          assert (iRow==i);
     560          dj[n]= -cost[iCol]/element[j];
     561          whenUsed_[n++]=iCol;
     562        }
     563        for (j=0;j<n;j++) {
     564          int jCol = whenUsed_[j];
     565          nextSlack[jCol]=-2;
     566        }
     567        CoinSort_2(dj,dj+n,whenUsed_);
     568        // put back
     569        iCol = whenUsed_[0];
     570        negSlack[i]=iCol;
     571        for (j=1;j<n;j++) {
     572          int jCol = whenUsed_[j];
     573          nextSlack[iCol]=jCol;
     574          iCol=jCol;
     575        }
     576      }
     577    }
     578  }
     579  whenUsed=whenUsed_;
    253580  if (model_->getObjSense()==-1.0) {
    254581    maxmin=-1.0;
     
    388715                       0,mu,drop,
    389716                       maxmin,offset,strategy,djTol,djExit,djFlag);
    390   n=0;
    391   for (i=ordStart;i<ordEnd;i++) {
    392     if (colsol[i]>lower[i]+fixTolerance) {
    393       if (colsol[i]<upper[i]-fixTolerance) {
    394         n++;
    395       } else {
    396         colsol[i]=upper[i];
    397       }
    398       whenUsed[i]=iteration;
    399     } else {
    400       colsol[i]=lower[i];
    401     }
     717  // update whenUsed_
     718  n = cleanIteration(iteration, ordStart,ordEnd,
     719                     colsol,  lower,  upper,
     720                     rowlower, rowupper,
     721                     cost, elemXX, fixTolerance,lastResult.objval,lastResult.infeas);
     722  if ((strategy_&16384)!=0) {
     723    int * posSlack = whenUsed_+ncols;
     724    int * negSlack = posSlack+nrows;
     725    int * nextSlack = negSlack + nrows;
     726    double * rowsol2 = (double *) (nextSlack+ncols);
     727    for (i=0;i<nrows;i++)
     728      rowsol[i] += rowsol2[i];
    402729  }
    403730  if ((logLevel_&1)!=0) {
     
    435762      bestFeasible=lastResult.objval;
    436763    }
     764    if (lastResult.infeas+mu*lastResult.objval<bestWeighted) {
     765      bestWeighted=lastResult.objval+mu*lastResult.objval;
     766    }
    437767    if ((saveStrategy&4096)) strategy |=256;
    438768    if ((saveStrategy&4)!=0&&iteration>2) {
     
    449779                     maxIts,mu,drop,
    450780                     maxmin,offset,strategy,djTol,djExit,djFlag);
    451     n=0;
    452     for (i=ordStart;i<ordEnd;i++) {
    453       if (colsol[i]>lower[i]+fixTolerance) {
    454         if (colsol[i]<upper[i]-fixTolerance) {
    455           n++;
    456         } else {
    457           colsol[i]=upper[i];
    458         }
    459         whenUsed[i]=iteration;
    460       } else {
    461         colsol[i]=lower[i];
    462       }
     781    n = cleanIteration(iteration, ordStart,ordEnd,
     782                       colsol,  lower,  upper,
     783                       rowlower, rowupper,
     784                       cost, elemXX, fixTolerance,result.objval,result.infeas);
     785    if ((strategy_&16384)!=0) {
     786      int * posSlack = whenUsed_+ncols;
     787      int * negSlack = posSlack+nrows;
     788      int * nextSlack = negSlack + nrows;
     789      double * rowsol2 = (double *) (nextSlack+ncols);
     790      for (i=0;i<nrows;i++)
     791        rowsol[i] += rowsol2[i];
    463792    }
    464793    if ((logLevel_&1)!=0) {
     
    515844        if ((strategy_&1024)!=0) mu *= 1.0e-1;
    516845      }
    517     } else if (result.infeas<bestInfeas) {
    518       // Save best solution
    519       memcpy(saveSol,colsol,ncols*sizeof(double));
     846    } else {
     847      bool save=false;
     848      if (result.infeas<=smallInfeas) {
     849        if (result.objval<bestFeasible)
     850          save=true;
     851      } else if (result.infeas+mu*result.objval<bestWeighted) {
     852        save=true;
     853      }
     854      if (save) {
     855        // Save best solution
     856        memcpy(saveSol,colsol,ncols*sizeof(double));
     857      }
    520858    }
    521859    bestInfeas=CoinMin(bestInfeas,result.infeas);
     
    634972    if (iteration>=majorIterations_) {
    635973      // If not feasible and crash then dive dive dive
    636       if (mu_>1.0e-12&&result.infeas>1.0&&majorIterations_<40) {
    637         mu_=1.0e-30;
     974      if (mu>1.0e-12&&result.infeas>1.0&&majorIterations_<40) {
     975        mu=1.0e-30;
    638976        majorIterations_=iteration+1;
     977        stopMu=0.0;
    639978      } else {
    640979        if (logLevel>2)
     
    7051044  }
    7061045  muAtExit_=mu;
    707   n=0;
    708   for (i=ordStart;i<ordEnd;i++) {
    709     if (colsol[i]>lower[i]+fixTolerance) {
    710       n++;
    711       whenUsed[i]=iteration;
    712     } else {
    713       colsol[i]=lower[i];
    714     }
     1046  // not scaled
     1047  n = cleanIteration(iteration, ordStart,ordEnd,
     1048                     colsol,  lower,  upper,
     1049                     model_->rowLower(), model_->rowUpper(),
     1050                     cost, element, fixTolerance,lastResult.objval,lastResult.infeas);
     1051  if ((strategy_&16384)!=0) {
     1052    int * posSlack = whenUsed_+ncols;
     1053    int * negSlack = posSlack+nrows;
     1054    int * nextSlack = negSlack + nrows;
     1055    double * rowsol2 = (double *) (nextSlack+ncols);
     1056    for (i=0;i<nrows;i++)
     1057      rowsol[i] += rowsol2[i];
    7151058  }
    7161059  if ((logLevel&1)==0) {
  • branches/devel/Clp/src/Idiot.hpp

    r982 r987  
    205205                   double * solExtra, double * elemExtra, double * upperExtra,
    206206                   double * costExtra,double weight);
    207 
     207  // Deals with whenUsed and slacks
     208  int cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd,
     209                     double * colsol, const double * lower, const double * upper,
     210                     const double * rowLower, const double * rowUpper,
     211                     const double * cost, const double * element, double fixTolerance,double & objChange,
     212                     double & infChange);
    208213private:
    209214  /// Underlying model
     
    241246                  2048 - keep lambda across mu change
    242247                  4096 - return best solution (not last found)
    243                   8192 - always do a presolve in crossover */
     248                  8192 - always do a presolve in crossover
     249                 16384 - costed slacks found - so whenUsed_ longer */
    244250  int lightWeight_; // 0 - normal, 1 lightweight
    245251};
Note: See TracChangeset for help on using the changeset viewer.