Changeset 156 for trunk


Ignore:
Timestamp:
Apr 4, 2003 1:48:05 PM (17 years ago)
Author:
forrest
Message:

Start of piecewiae

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpNonLinearCost.cpp

    r137 r156  
    100100
    101101}
    102 
    103102ClpNonLinearCost::ClpNonLinearCost(ClpSimplex * model,const int * starts,
    104103                   const double * lowerNon, const double * costNon)
     
    116115  offset_ = new int [numberTotal];
    117116  memset(offset_,0,numberTotal*sizeof(int));
     117 
     118  double whichWay = model_->optimizationDirection();
     119  printf("Direction %g\n",whichWay);
    118120
    119121  numberInfeasibilities_=0;
     
    124126
    125127  int iSequence;
    126   double * upper = model_->upperRegion();
    127   double * lower = model_->lowerRegion();
    128   double * cost = model_->costRegion();
    129   // First see how much space we need - we know column part
    130   // but not infeasibilities - see how much extra
    131   // may be over-estimate
    132   int put=starts[numberColumns_]-numberColumns_;
    133 
    134   for (iSequence=0;iSequence<numberTotal;iSequence++) {
    135     if (lower[iSequence]>-1.0e20)
    136       if (upper[iSequence]<1.0e20)
    137         put++;
    138     put += 3;
    139   }
    140 
     128  assert (!model_->rowObjective());
     129  double * cost = model_->objective();
     130
     131  // First see how much space we need
     132  // Also set up feasible bounds
     133  int put=starts[numberColumns_];
     134
     135  double * columnUpper = model_->columnUpper();
     136  double * columnLower = model_->columnLower();
     137  for (iSequence=0;iSequence<numberColumns_;iSequence++) {
     138    if (columnLower[iSequence]>-1.0e20)
     139      put++;
     140    if (columnUpper[iSequence]<1.0e20)
     141      put++;
     142  }
     143
     144  double * rowUpper = model_->rowUpper();
     145  double * rowLower = model_->rowLower();
     146  for (iSequence=0;iSequence<numberRows_;iSequence++) {
     147    if (rowLower[iSequence]>-1.0e20)
     148      put++;
     149    if (rowUpper[iSequence]<1.0e20)
     150      put++;
     151    put +=2;
     152  }
     153  printf("put %d\n",put);
    141154  lower_ = new double [put];
    142155  cost_ = new double [put];
     
    148161
    149162  start_[0]=0;
    150 
    151163  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    152164    lower_[put] = -COIN_DBL_MAX;
    153165    setInfeasible(put,true);
    154     cost_[put++] = cost[iSequence]-infeasibilityCost;
    155     whichRange_[iSequence]=put;
     166    whichRange_[iSequence]=put+1;
    156167    double thisCost;
     168    double lowerValue;
     169    double upperValue;
    157170    if (iSequence>=numberColumns_) {
    158171      // rows
    159       lower_[put] = lower[iSequence];
    160       cost_[put++] = cost[iSequence];
    161       thisCost = cost[iSequence];
     172      lowerValue = rowLower[iSequence-numberColumns_];
     173      upperValue = rowUpper[iSequence-numberColumns_];
     174      if (lowerValue>-1.0e30) {
     175        cost_[put++] = -infeasibilityCost;
     176        lower_[put] = lowerValue;
     177      }
     178      cost_[put++] = 0.0;
     179      thisCost = 0.0;
    162180    } else {
    163181      // columns - move costs and see if convex
     182      lowerValue = columnLower[iSequence];
     183      upperValue = columnUpper[iSequence];
     184      if (lowerValue>-1.0e30) {
     185        cost_[put++] = whichWay*cost[iSequence]-infeasibilityCost;
     186        lower_[put] = lowerValue;
     187      }
    164188      int iIndex = starts[iSequence];
    165189      int end = starts[iSequence+1];
    166       assert (fabs(lower[iSequence]-lowerNon[iIndex])<1.0e-8);
     190      assert (fabs(columnLower[iSequence]-lowerNon[iIndex])<1.0e-8);
    167191      thisCost = -COIN_DBL_MAX;
    168192      for (;iIndex<end;iIndex++) {
    169         if (lowerNon[iIndex]<upper[iSequence]-1.0e-8) {
     193        if (lowerNon[iIndex]<columnUpper[iSequence]-1.0e-8) {
    170194          lower_[put] = lowerNon[iIndex];
    171           cost_[put++] = costNon[iIndex];
     195          cost_[put++] = whichWay*costNon[iIndex];
    172196          // check convexity
    173           if (costNon[iIndex]<thisCost-1.0e-12)
     197          if (whichWay*costNon[iIndex]<thisCost-1.0e-12)
    174198            convex_ = false;
    175           thisCost = costNon[iIndex];
     199          thisCost = whichWay*costNon[iIndex];
    176200        } else {
    177201          break;
     
    179203      }
    180204    }
    181     lower_[put] = upper[iSequence];
     205    lower_[put] = upperValue;
    182206    cost_[put++] = thisCost+infeasibilityCost;
    183     if (upper[iSequence]<1.0e20) {
     207    if (upperValue<1.0e20) {
    184208      lower_[put] = COIN_DBL_MAX;
    185209      setInfeasible(put-1,true);
     
    188212    start_[iSequence+1]=put;
    189213  }
     214  printf("put %d\n",put);
    190215  // can't handle non-convex at present
    191216  assert(convex_);
  • trunk/ClpSimplex.cpp

    r153 r156  
    41614161  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
    41624162}
    4163 // Pass in an existing non linear cost
    4164 void
    4165 ClpSimplex::passInNonLinearCost(ClpNonLinearCost * cost)
     4163/* Constructs a non linear cost from list of non-linearities (columns only)
     4164   First lower of each column is taken as real lower
     4165   Last lower is taken as real upper and cost ignored
     4166   
     4167   Returns nonzero if bad data e.g. lowers not monotonic
     4168*/
     4169int
     4170ClpSimplex::createPiecewiseLinearCosts(const int * starts,
     4171                                       const double * lower, const double * gradient)
    41664172{
    41674173  delete nonLinearCost_;
    4168   nonLinearCost_ = cost;
    4169 }
     4174  // Set up feasible bounds and check monotonicity
     4175  int iColumn;
     4176  int returnCode=0;
     4177
     4178  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     4179    int iIndex = starts[iColumn];
     4180    int end = starts[iColumn+1]-1;
     4181    columnLower_[iColumn] = lower[iIndex];
     4182    columnUpper_[iColumn] = lower[end];
     4183    double value = columnLower_[iColumn];
     4184    iIndex++;
     4185    for (;iIndex<end;iIndex++) {
     4186      if (lower[iIndex]<value)
     4187        returnCode++; // not monotonic
     4188      value=lower[iIndex];
     4189    }
     4190  }
     4191  nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);
     4192  return returnCode;
     4193}
  • trunk/ClpSimplexPrimal.cpp

    r145 r156  
    430430          return;
    431431        }
     432#if 1
     433        internalFactorize(0);
     434#else
     435
    432436        // no - restore previous basis
    433437        assert (type==1);
     
    440444        type = 2;
    441445        assert (internalFactorize(1)==0);
     446#endif
    442447        changeMade_++; // say change made
    443448      }
     
    635640    }
    636641  }
    637   if (type==0||type==1) {
    638     if (!type) {
     642  if (type==0||type==1||(type==3&&!numberIterations_)) {
     643    if (!type||(type==3&&!numberIterations_)) {
    639644      // create save arrays
    640645      delete [] saveStatus_;
     
    12561261    specialOptions_ &= ~1;
    12571262}
     1263bool
     1264ClpSimplexPrimal::alwaysOptimal() const
     1265{
     1266  return (specialOptions_&1)!=0;
     1267}
    12581268/*
    12591269  Reasons to come out (normal mode/user mode):
  • trunk/ClpSimplexPrimalQuadratic.cpp

    r152 r156  
    227227      // bad pass - restore solution
    228228      memcpy(solution,saveSolution,numberColumns*sizeof(double));
     229      iPass--;
    229230    }
    230231  }
     
    285286ClpSimplexPrimalQuadratic::primalBeale()
    286287{
    287   abort();
     288  // Wolfe's method looks easier - lets try that
     289  // This is as a user would see
     290
     291  int numberColumns = this->numberColumns();
     292  double * columnLower = this->columnLower();
     293  double * columnUpper = this->columnUpper();
     294  double * objective = this->objective();
     295  double * solution = this->primalColumnSolution();
     296  double * dj = this->dualColumnSolution();
     297  double * pi = this->dualRowSolution();
     298
     299  int numberRows = this->numberRows();
     300  double * rowLower = this->rowLower();
     301  double * rowUpper = this->rowUpper();
     302
     303  // and elements
     304  CoinPackedMatrix * matrix = this->matrix();
     305  const int * row = matrix->getIndices();
     306  const int * columnStart = matrix->getVectorStarts();
     307  const double * element =  matrix->getElements();
     308  const int * columnLength = matrix->getVectorLengths();
     309
     310  // Get list of non linear columns
     311  CoinPackedMatrix * quadratic = quadraticObjective();
     312  if (!quadratic||!quadratic->getNumElements()) {
     313    // no quadratic part
     314    return primal(1);
     315  }
     316
     317  int iColumn;
     318  const int * columnQuadratic = quadratic->getIndices();
     319  const int * columnQuadraticStart = quadratic->getVectorStarts();
     320  const int * columnQuadraticLength = quadratic->getVectorLengths();
     321  const double * quadraticElement = quadratic->getElements();
     322  // Get a feasible solution
     323  if (numberPrimalInfeasibilities())
     324    primal(1);
     325  // still infeasible
     326  if (numberPrimalInfeasibilities())
     327    return 0;
     328 
     329  // Create larger problem
     330  int newNumberRows = numberRows+numberColumns;
     331  // See how many artificials we will need
     332  // For now assume worst
     333  int newNumberColumns = 3*numberColumns+ numberRows;
     334  int numberElements = 2*matrix->getNumElements()
     335    +quadratic->getNumElements()
     336    + 2*numberColumns;
     337  double * elements2 = new double[numberElements];
     338  int * start2 = new int[newNumberColumns+1];
     339  int * row2 = new int[numberElements];
     340  double * objective2 = new double[newNumberColumns];
     341  double * columnLower2 = new double[newNumberColumns];
     342  double * columnUpper2 = new double[newNumberColumns];
     343  double * rowLower2 = new double[newNumberRows];
     344  double * rowUpper2 = new double[newNumberRows];
     345  memset(rowLower2,0,newNumberRows*sizeof(double));
     346  memcpy(rowLower2,rowLower,numberRows*sizeof(double));
     347  memcpy(rowLower2+numberRows,objective,numberColumns*sizeof(double));
     348  memset(rowUpper2,0,newNumberRows*sizeof(double));
     349  memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
     350  memcpy(rowUpper2+numberRows,objective,numberColumns*sizeof(double));
     351  memset(objective2,0,newNumberColumns*sizeof(double));
     352  // Get a row copy of quadratic objective in standard format
     353  CoinPackedMatrix copyQ;
     354  copyQ.reverseOrderedCopyOf(*quadratic);
     355  const int * columnQ = copyQ.getIndices();
     356  const CoinBigIndex * rowStartQ = copyQ.getVectorStarts();
     357  const int * rowLengthQ = copyQ.getVectorLengths();
     358  const double * elementByRowQ = copyQ.getElements();
     359  newNumberColumns=0;
     360  numberElements=0;
     361  start2[0]=0;
     362  // x
     363  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     364    // Original matrix
     365    columnLower2[iColumn]=columnLower[iColumn];
     366    columnUpper2[iColumn]=columnUpper[iColumn];
     367    int j;
     368    for (j=columnStart[iColumn];
     369         j<columnStart[iColumn]+columnLength[iColumn];
     370         j++) {
     371      elements2[numberElements]=element[j];
     372      row2[numberElements++]=row[j];
     373    }
     374    // Quadratic and modify djs
     375    for (j=columnQuadraticStart[iColumn];
     376         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
     377         j++) {
     378      int jColumn = columnQuadratic[j];
     379      double value = quadraticElement[j];
     380      if (iColumn!=jColumn)
     381        value *= 0.5;
     382      dj[iColumn] += solution[jColumn]*value;
     383      elements2[numberElements]=-value;
     384      row2[numberElements++]=jColumn+numberRows;
     385    }
     386    for (j=rowStartQ[iColumn];
     387         j<rowStartQ[iColumn]+rowLengthQ[iColumn];
     388         j++) {
     389      int jColumn = columnQ[j];
     390      double value = elementByRowQ[j];
     391      if (iColumn!=jColumn) {
     392        value *= 0.5;
     393        dj[iColumn] += solution[jColumn]*value;
     394        elements2[numberElements]=-value;
     395        row2[numberElements++]=jColumn+numberRows;
     396      }
     397    }
     398    start2[iColumn+1]=numberElements;
     399  }
     400  newNumberColumns=numberColumns;
     401  // pi
     402  int iRow;
     403  // Get a row copy in standard format
     404  CoinPackedMatrix copy;
     405  copy.reverseOrderedCopyOf(*(this->matrix()));
     406  // get matrix data pointers
     407  const int * column = copy.getIndices();
     408  const CoinBigIndex * rowStart = copy.getVectorStarts();
     409  const int * rowLength = copy.getVectorLengths();
     410  const double * elementByRow = copy.getElements();
     411  for (iRow=0;iRow<numberRows;iRow++) {
     412    // should look at rows to get correct bounds
     413    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
     414    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     415    int j;
     416    for (j=rowStart[iRow];
     417         j<rowStart[iRow]+rowLength[iRow];
     418         j++) {
     419      elements2[numberElements]=elementByRow[j];
     420      row2[numberElements++]=column[j]+numberRows;
     421    }
     422    newNumberColumns++;
     423    start2[newNumberColumns]=numberElements;
     424  }
     425  // djs and artificials
     426  double tolerance = dualTolerance();
     427  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     428    // dj
     429    columnLower2[newNumberColumns]=0.0;
     430    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     431    elements2[numberElements]=1.0;
     432    row2[numberElements++]=iColumn+numberRows;
     433    newNumberColumns++;
     434    start2[newNumberColumns]=numberElements;
     435    // artificial (assuming no bounds)
     436    if (getStatus(iColumn)==basic||solution[iColumn]>1.0e-7) {
     437      if (fabs(dj[iColumn])>tolerance) {
     438        columnLower2[newNumberColumns]=0.0;
     439        columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     440        objective2[newNumberColumns]=1.0;
     441        if (dj[iColumn]>0.0)
     442          elements2[numberElements]=1.0;
     443        else
     444          elements2[numberElements]=-1.0;
     445        row2[numberElements++]=iColumn+numberRows;
     446        newNumberColumns++;
     447        start2[newNumberColumns]=numberElements;
     448      }
     449    } else if (dj[iColumn]<-tolerance) {
     450      columnLower2[newNumberColumns]=0.0;
     451      columnUpper2[newNumberColumns]=COIN_DBL_MAX;
     452      objective2[newNumberColumns]=1.0;
     453      elements2[numberElements]=-1.0;
     454      row2[numberElements++]=iColumn+numberRows;
     455      newNumberColumns++;
     456      start2[newNumberColumns]=numberElements;
     457    }
     458  }
     459  // Create model
     460  ClpSimplex model2;
     461  model2.loadProblem(newNumberColumns,newNumberRows,
     462                     start2,row2, elements2,
     463                     columnLower2,columnUpper2,
     464                     objective2,
     465                     rowLower2,rowUpper2);
     466  model2.allSlackBasis();
     467  // Move solution across
     468  double * solution2 = model2.primalColumnSolution();
     469  memcpy(solution2,solution,numberColumns*sizeof(double));
     470  memcpy(solution2+numberColumns,pi,numberRows*sizeof(double));
     471  for (iRow=0;iRow<numberRows;iRow++)
     472    model2.setRowStatus(iRow,getRowStatus(iRow));
     473  // djs and artificials
     474  newNumberColumns = numberRows+numberColumns;
     475  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     476    model2.setStatus(iColumn,getStatus(iColumn));
     477    if (getStatus(iColumn)==basic) {
     478      solution2[newNumberColumns++]=0.0;
     479      if (fabs(dj[iColumn])>tolerance) {
     480        solution2[newNumberColumns++]=fabs(dj[iColumn]);
     481      }
     482    } else if (dj[iColumn]<-tolerance) {
     483      solution2[newNumberColumns++]=0.0;
     484      solution2[newNumberColumns++]=-dj[iColumn];
     485    } else {
     486      solution2[newNumberColumns++]=dj[iColumn];
     487    }
     488  }
     489  memset(model2.primalRowSolution(),0,newNumberRows*sizeof(double));
     490  model2.times(1.0,model2.primalColumnSolution(),
     491               model2.primalRowSolution());
     492  // solve
     493  model2.primal(1);
    288494  return 0;
    289495}
  • trunk/include/ClpSimplex.hpp

    r153 r156  
    319319      bounds.  Returns feasibility states */
    320320  int getSolution ();
    321   /// Pass in an existing non linear cost
    322   void passInNonLinearCost(ClpNonLinearCost * cost);
     321  /** Constructs a non linear cost from list of non-linearities (columns only)
     322      First lower of each column is taken as real lower
     323      Last lower is taken as real upper and cost ignored
     324
     325      Returns nonzero if bad data e.g. lowers not monotonic
     326  */
     327  int createPiecewiseLinearCosts(const int * starts,
     328                   const double * lower, const double * gradient);
    323329  /** Return model - updates any scalars */
    324330  void returnModel(ClpSimplex & otherModel);
  • trunk/include/ClpSimplexPrimal.hpp

    r112 r156  
    118118  /// Do not change infeasibility cost and always say optimal
    119119  void alwaysOptimal(bool onOff);
     120  bool alwaysOptimal() const;
    120121  //@}
    121122
Note: See TracChangeset for help on using the changeset viewer.