source: trunk/ClpSimplexPrimalQuadratic.cpp @ 162

Last change on this file since 162 was 162, checked in by forrest, 17 years ago

Piecewise linear

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.5 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include <math.h>
7
8#include "CoinHelperFunctions.hpp"
9#include "ClpSimplexPrimalQuadratic.hpp"
10#include "ClpFactorization.hpp"
11#include "ClpNonLinearCost.hpp"
12#include "CoinPackedMatrix.hpp"
13#include "CoinIndexedVector.hpp"
14#include "CoinWarmStartBasis.hpp"
15#include "ClpPrimalColumnPivot.hpp"
16#include "ClpMessage.hpp"
17#include <cfloat>
18#include <cassert>
19#include <string>
20#include <stdio.h>
21#include <iostream>
22
23// A sequential LP method
24int 
25ClpSimplexPrimalQuadratic::primalSLP(int numberPasses, double deltaTolerance)
26{
27  // Are we minimizing or maximizing
28  double whichWay=optimizationDirection();
29  // This is as a user would see
30
31  int numberColumns = this->numberColumns();
32  int numberRows = this->numberRows();
33  double * columnLower = this->columnLower();
34  double * columnUpper = this->columnUpper();
35  double * objective = this->objective();
36  double * solution = this->primalColumnSolution();
37 
38  // Save objective
39 
40  double * saveObjective = new double [numberColumns];
41  memcpy(saveObjective,objective,numberColumns*sizeof(double));
42
43  // Get list of non linear columns
44  CoinPackedMatrix * quadratic = quadraticObjective();
45  if (!quadratic) {
46    // no quadratic part
47    return primal(0);
48  }
49  int numberNonLinearColumns = 0;
50  int iColumn;
51  int * listNonLinearColumn = new int[numberColumns];
52  memset(listNonLinearColumn,0,numberColumns*sizeof(int));
53  const int * columnQuadratic = quadratic->getIndices();
54  const int * columnQuadraticStart = quadratic->getVectorStarts();
55  const int * columnQuadraticLength = quadratic->getVectorLengths();
56  const double * quadraticElement = quadratic->getElements();
57  for (iColumn=0;iColumn<numberColumns;iColumn++) {
58    int j;
59    for (j=columnQuadraticStart[iColumn];
60         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
61      int jColumn = columnQuadratic[j];
62      listNonLinearColumn[jColumn]=1;
63      listNonLinearColumn[iColumn]=1;
64    }
65  }
66  for (iColumn=0;iColumn<numberColumns;iColumn++) {
67    if(listNonLinearColumn[iColumn])
68      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
69  }
70 
71  if (!numberNonLinearColumns) {
72    delete [] listNonLinearColumn;
73    // no quadratic part
74    return primal(0);
75  }
76
77  // get feasible
78  if (numberPrimalInfeasibilities())
79    primal(1);
80  // still infeasible
81  if (numberPrimalInfeasibilities())
82    return 0;
83
84  int jNon;
85  int * last[3];
86 
87  double * trust = new double[numberNonLinearColumns];
88  double * trueLower = new double[numberNonLinearColumns];
89  double * trueUpper = new double[numberNonLinearColumns];
90  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
91    iColumn=listNonLinearColumn[jNon];
92    trust[jNon]=0.5;
93    trueLower[jNon]=columnLower[iColumn];
94    trueUpper[jNon]=columnUpper[iColumn];
95  }
96  int iPass;
97  double lastObjective=1.0e31;
98  double * saveSolution = new double [numberColumns];
99  double * savePi = new double [numberRows];
100  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
101  double targetDrop=1.0e31;
102  double objectiveOffset;
103  getDblParam(ClpObjOffset,objectiveOffset);
104  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
105  for (iPass=0;iPass<3;iPass++) {
106    last[iPass]=new int[numberNonLinearColumns];
107    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
108      last[iPass][jNon]=0;
109  }
110  double goodMove=false;
111  char * statusCheck = new char[numberColumns];
112  for (iPass=0;iPass<numberPasses;iPass++) {
113    // redo objective
114    double offset=0.0;
115    double objValue=-objectiveOffset;
116    memcpy(objective,saveObjective,numberColumns*sizeof(double));
117    for (iColumn=0;iColumn<numberColumns;iColumn++) 
118      objValue += objective[iColumn]*solution[iColumn];
119    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
120      iColumn=listNonLinearColumn[jNon];
121      if (getColumnStatus(iColumn)==basic) {
122        if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
123          statusCheck[iColumn]='l';
124        else if (solution[iColumn]>columnUpper[iColumn]-1.0e-8)
125          statusCheck[iColumn]='u';
126        else
127          statusCheck[iColumn]='B';
128      } else {
129        if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
130          statusCheck[iColumn]='L';
131        else
132          statusCheck[iColumn]='U';
133      }
134      double valueI = solution[iColumn];
135      int j;
136      for (j=columnQuadraticStart[iColumn];
137           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
138        int jColumn = columnQuadratic[j];
139        double valueJ = solution[jColumn];
140        double elementValue = quadraticElement[j];
141        objValue += 0.5*valueI*valueJ*elementValue;
142        offset += 0.5*valueI*valueJ*elementValue;
143        double gradientI = valueJ*elementValue;
144        double gradientJ = valueI*elementValue;
145        offset -= gradientI*valueI;
146        objective[iColumn] += gradientI;
147        offset -= gradientJ*valueJ;
148        objective[jColumn] += gradientJ;
149      }
150    }
151    printf("objective %g, objective offset %g\n",objValue,offset);
152    setDblParam(ClpObjOffset,objectiveOffset-offset);
153    objValue *= whichWay;
154    int * temp=last[2];
155    last[2]=last[1];
156    last[1]=last[0];
157    last[0]=temp;
158    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
159      iColumn=listNonLinearColumn[jNon];
160      double change = solution[iColumn]-saveSolution[iColumn];
161      if (change<-1.0e-5) {
162        if (fabs(change+trust[jNon])<1.0e-5) 
163          temp[jNon]=-1;
164        else
165          temp[jNon]=-2;
166      } else if(change>1.0e-5) {
167        if (fabs(change-trust[jNon])<1.0e-5) 
168          temp[jNon]=1;
169        else
170          temp[jNon]=2;
171      } else {
172        temp[jNon]=0;
173      }
174    } 
175    if (objValue<=lastObjective) 
176      goodMove=true;
177    else
178      goodMove=false;
179    double maxDelta=0.0;
180    double maxGap=0.0;
181    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
182      iColumn=listNonLinearColumn[jNon];
183      maxDelta = max(maxDelta,
184                     fabs(solution[iColumn]-saveSolution[iColumn]));
185      if (goodMove) {
186        if (last[0][jNon]*last[1][jNon]<0) {
187          // halve
188          trust[jNon] *= 0.5;
189        } else {
190          if (last[0][jNon]==last[1][jNon]&&
191              last[0][jNon]==last[2][jNon])
192            trust[jNon] *= 1.5; 
193        }
194      } else if (trust[jNon]>10.0*deltaTolerance) {
195        trust[jNon] *= 0.5;
196      }
197      maxGap = max(maxGap,trust[jNon]);
198    }
199    std::cout<<"largest gap is "<<maxGap<<std::endl;
200    if (goodMove) {
201      double drop = lastObjective-objValue;
202      std::cout<<"True drop was "<<drop<<std::endl;
203      std::cout<<"largest delta is "<<maxDelta<<std::endl;
204      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove) {
205        std::cout<<"Exiting"<<std::endl;
206        break;
207      }
208    } else {
209      //lastObjective += 1.0e-3; // to stop exit
210    }
211    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
212      iColumn=listNonLinearColumn[jNon];
213      columnLower[iColumn]=max(solution[iColumn]
214                               -trust[jNon],
215                               trueLower[jNon]);
216      columnUpper[iColumn]=min(solution[iColumn]
217                               +trust[jNon],
218                               trueUpper[jNon]);
219    }
220    if (goodMove) {
221      memcpy(saveSolution,solution,numberColumns*sizeof(double));
222      memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
223      memcpy(saveStatus,status_,numberRows+numberColumns);
224     
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      }
241      std::cout<<"Pass - "<<iPass
242               <<", target drop is "<<targetDrop
243               <<std::endl;
244      lastObjective = objValue;
245      if (targetDrop<1.0e-5&&goodMove&&iPass) {
246        printf("Exiting on target drop %g\n",targetDrop);
247        break;
248      }
249      {
250        double * r = this->dualColumnSolution();
251        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
252          iColumn=listNonLinearColumn[jNon];
253          printf("Trust %d %g - solution %d %g obj %g dj %g state %c\n",
254                 jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
255                 r[iColumn],statusCheck[iColumn]);
256        }
257      }
258      setLogLevel(63);
259      this->primal(1);
260    } else {
261      // bad pass - restore solution
262      printf("Backtracking\n");
263      memcpy(solution,saveSolution,numberColumns*sizeof(double));
264      memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
265      memcpy(status_,saveStatus,numberRows+numberColumns);
266      iPass--;
267    }
268  }
269  // restore solution
270  memcpy(solution,saveSolution,numberColumns*sizeof(double));
271  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
272    iColumn=listNonLinearColumn[jNon];
273    columnLower[iColumn]=max(solution[iColumn],
274                             trueLower[jNon]);
275    columnUpper[iColumn]=min(solution[iColumn],
276                             trueUpper[jNon]);
277  }
278  delete [] statusCheck;
279  delete [] savePi;
280  delete [] saveStatus;
281  // redo objective
282  double offset=0.0;
283  double objValue=-objectiveOffset;
284  memcpy(objective,saveObjective,numberColumns*sizeof(double));
285  for (iColumn=0;iColumn<numberColumns;iColumn++) 
286    objValue += objective[iColumn]*solution[iColumn];
287  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
288    iColumn=listNonLinearColumn[jNon];
289    double valueI = solution[iColumn];
290    int j;
291    for (j=columnQuadraticStart[iColumn];
292         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
293      int jColumn = columnQuadratic[j];
294      double valueJ = solution[jColumn];
295      double elementValue = quadraticElement[j];
296      objValue += 0.5*valueI*valueJ*elementValue;
297      offset += 0.5*valueI*valueJ*elementValue;
298      double gradientI = valueJ*elementValue;
299      double gradientJ = valueI*elementValue;
300      offset -= gradientI*valueI;
301      objective[iColumn] += gradientI;
302      offset -= gradientJ*valueJ;
303      objective[jColumn] += gradientJ;
304    }
305  }
306  printf("objective %g, objective offset %g\n",objValue,offset);
307  setDblParam(ClpObjOffset,objectiveOffset-offset);
308  this->primal(1);
309  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
310    iColumn=listNonLinearColumn[jNon];
311    columnLower[iColumn]= trueLower[jNon];
312    columnUpper[iColumn]= trueUpper[jNon];
313  }
314  delete [] saveSolution;
315  for (iPass=0;iPass<3;iPass++) 
316    delete [] last[iPass];
317  delete [] trust;
318  delete [] trueUpper;
319  delete [] trueLower;
320  delete [] saveObjective;
321  delete [] listNonLinearColumn;
322  return 0;
323}
324// Beale's method
325int 
326ClpSimplexPrimalQuadratic::primalBeale()
327{
328  // Wolfe's method looks easier - lets try that
329  // This is as a user would see
330
331  int numberColumns = this->numberColumns();
332  double * columnLower = this->columnLower();
333  double * columnUpper = this->columnUpper();
334  double * objective = this->objective();
335  double * solution = this->primalColumnSolution();
336  double * dj = this->dualColumnSolution();
337  double * pi = this->dualRowSolution();
338
339  int numberRows = this->numberRows();
340  double * rowLower = this->rowLower();
341  double * rowUpper = this->rowUpper();
342
343  // and elements
344  CoinPackedMatrix * matrix = this->matrix();
345  const int * row = matrix->getIndices();
346  const int * columnStart = matrix->getVectorStarts();
347  const double * element =  matrix->getElements();
348  const int * columnLength = matrix->getVectorLengths();
349
350  // Get list of non linear columns
351  CoinPackedMatrix * quadratic = quadraticObjective();
352  if (!quadratic||!quadratic->getNumElements()) {
353    // no quadratic part
354    return primal(1);
355  }
356
357  int iColumn;
358  const int * columnQuadratic = quadratic->getIndices();
359  const int * columnQuadraticStart = quadratic->getVectorStarts();
360  const int * columnQuadraticLength = quadratic->getVectorLengths();
361  const double * quadraticElement = quadratic->getElements();
362  // Get a feasible solution
363  if (numberPrimalInfeasibilities())
364    primal(1);
365  // still infeasible
366  if (numberPrimalInfeasibilities())
367    return 0;
368 
369  // Create larger problem
370  int newNumberRows = numberRows+numberColumns;
371  // See how many artificials we will need
372  // For now assume worst
373  int newNumberColumns = 3*numberColumns+ numberRows;
374  int numberElements = 2*matrix->getNumElements()
375    +2*quadratic->getNumElements()
376    + 2*numberColumns;
377  double * elements2 = new double[numberElements];
378  int * start2 = new int[newNumberColumns+1];
379  int * row2 = new int[numberElements];
380  double * objective2 = new double[newNumberColumns];
381  double * columnLower2 = new double[newNumberColumns];
382  double * columnUpper2 = new double[newNumberColumns];
383  double * rowLower2 = new double[newNumberRows];
384  double * rowUpper2 = new double[newNumberRows];
385  memset(rowLower2,0,newNumberRows*sizeof(double));
386  memcpy(rowLower2,rowLower,numberRows*sizeof(double));
387  memcpy(rowLower2+numberRows,objective,numberColumns*sizeof(double));
388  memset(rowUpper2,0,newNumberRows*sizeof(double));
389  memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
390  memcpy(rowUpper2+numberRows,objective,numberColumns*sizeof(double));
391  memset(objective2,0,newNumberColumns*sizeof(double));
392  // Get a row copy of quadratic objective in standard format
393  CoinPackedMatrix copyQ;
394  copyQ.reverseOrderedCopyOf(*quadratic);
395  const int * columnQ = copyQ.getIndices();
396  const CoinBigIndex * rowStartQ = copyQ.getVectorStarts();
397  const int * rowLengthQ = copyQ.getVectorLengths(); 
398  const double * elementByRowQ = copyQ.getElements();
399  newNumberColumns=0;
400  numberElements=0;
401  start2[0]=0;
402  // x
403  for (iColumn=0;iColumn<numberColumns;iColumn++) {
404    // Original matrix
405    columnLower2[iColumn]=columnLower[iColumn];
406    columnUpper2[iColumn]=columnUpper[iColumn];
407    int j;
408    for (j=columnStart[iColumn];
409         j<columnStart[iColumn]+columnLength[iColumn];
410         j++) {
411      elements2[numberElements]=element[j];
412      row2[numberElements++]=row[j];
413    }
414    // Quadratic and modify djs
415    for (j=columnQuadraticStart[iColumn];
416         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
417         j++) {
418      int jColumn = columnQuadratic[j];
419      double value = quadraticElement[j];
420      if (iColumn!=jColumn) 
421        value *= 0.5;
422      dj[iColumn] += solution[jColumn]*value;
423      elements2[numberElements]=-value;
424      row2[numberElements++]=jColumn+numberRows;
425    }
426    for (j=rowStartQ[iColumn];
427         j<rowStartQ[iColumn]+rowLengthQ[iColumn];
428         j++) {
429      int jColumn = columnQ[j];
430      double value = elementByRowQ[j];
431      if (iColumn!=jColumn) { 
432        value *= 0.5;
433        dj[iColumn] += solution[jColumn]*value;
434        elements2[numberElements]=-value;
435        row2[numberElements++]=jColumn+numberRows;
436      }
437    }
438    start2[iColumn+1]=numberElements;
439  }
440  newNumberColumns=numberColumns;
441  // pi
442  int iRow;
443  // Get a row copy in standard format
444  CoinPackedMatrix copy;
445  copy.reverseOrderedCopyOf(*(this->matrix()));
446  // get matrix data pointers
447  const int * column = copy.getIndices();
448  const CoinBigIndex * rowStart = copy.getVectorStarts();
449  const int * rowLength = copy.getVectorLengths(); 
450  const double * elementByRow = copy.getElements();
451  for (iRow=0;iRow<numberRows;iRow++) {
452    // should look at rows to get correct bounds
453    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
454    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
455    int j;
456    for (j=rowStart[iRow];
457         j<rowStart[iRow]+rowLength[iRow];
458         j++) {
459      elements2[numberElements]=elementByRow[j];
460      row2[numberElements++]=column[j]+numberRows;
461    }
462    newNumberColumns++;
463    start2[newNumberColumns]=numberElements;
464  }
465  // djs and artificials
466  double tolerance = dualTolerance();
467  for (iColumn=0;iColumn<numberColumns;iColumn++) {
468    // dj
469    columnLower2[newNumberColumns]=0.0;
470    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
471    elements2[numberElements]=1.0;
472    row2[numberElements++]=iColumn+numberRows;
473    newNumberColumns++;
474    start2[newNumberColumns]=numberElements;
475    // artificial (assuming no bounds)
476    if (getStatus(iColumn)==basic||solution[iColumn]>1.0e-7) {
477      if (fabs(dj[iColumn])>tolerance) {
478        columnLower2[newNumberColumns]=0.0;
479        columnUpper2[newNumberColumns]=COIN_DBL_MAX;
480        objective2[newNumberColumns]=1.0;
481        if (dj[iColumn]>0.0)
482          elements2[numberElements]=1.0;
483        else
484          elements2[numberElements]=-1.0;
485        row2[numberElements++]=iColumn+numberRows;
486        newNumberColumns++;
487        start2[newNumberColumns]=numberElements;
488      }
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
500  ClpSimplex model2;
501  model2.loadProblem(newNumberColumns,newNumberRows,
502                     start2,row2, elements2,
503                     columnLower2,columnUpper2,
504                     objective2,
505                     rowLower2,rowUpper2);
506  model2.allSlackBasis();
507  // Move solution across
508  double * solution2 = model2.primalColumnSolution();
509  memcpy(solution2,solution,numberColumns*sizeof(double));
510  memcpy(solution2+numberColumns,pi,numberRows*sizeof(double));
511  for (iRow=0;iRow<numberRows;iRow++) 
512    model2.setRowStatus(iRow,getRowStatus(iRow));
513  // djs and artificials
514  newNumberColumns = numberRows+numberColumns;
515  for (iColumn=0;iColumn<numberColumns;iColumn++) {
516    model2.setStatus(iColumn,getStatus(iColumn));
517    if (getStatus(iColumn)==basic) {
518      solution2[newNumberColumns++]=0.0;
519      if (fabs(dj[iColumn])>tolerance) {
520        solution2[newNumberColumns++]=fabs(dj[iColumn]);
521      }
522    } else if (dj[iColumn]<-tolerance) {
523      solution2[newNumberColumns++]=0.0;
524      solution2[newNumberColumns++]=-dj[iColumn];
525    } else {
526      solution2[newNumberColumns++]=dj[iColumn];
527    }
528  }
529  memset(model2.primalRowSolution(),0,newNumberRows*sizeof(double));
530  model2.times(1.0,model2.primalColumnSolution(),
531               model2.primalRowSolution());
532  // solve
533  model2.primal(1);
534  return 0;
535}
536 
537
Note: See TracBrowser for help on using the repository browser.