source: trunk/ClpSimplexPrimalQuadratic.cpp @ 156

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

Start of piecewiae

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 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  double * columnLower = this->columnLower();
33  double * columnUpper = this->columnUpper();
34  double * objective = this->objective();
35  double * solution = this->primalColumnSolution();
36 
37  // Save objective
38 
39  double * saveObjective = new double [numberColumns];
40  memcpy(saveObjective,objective,numberColumns*sizeof(double));
41
42  // Get list of non linear columns
43  CoinPackedMatrix * quadratic = quadraticObjective();
44  if (!quadratic) {
45    // no quadratic part
46    return primal(0);
47  }
48  int numberNonLinearColumns = 0;
49  int iColumn;
50  int * listNonLinearColumn = new int[numberColumns];
51  memset(listNonLinearColumn,0,numberColumns*sizeof(int));
52  const int * columnQuadratic = quadratic->getIndices();
53  const int * columnQuadraticStart = quadratic->getVectorStarts();
54  const int * columnQuadraticLength = quadratic->getVectorLengths();
55  const double * quadraticElement = quadratic->getElements();
56  for (iColumn=0;iColumn<numberColumns;iColumn++) {
57    int j;
58    for (j=columnQuadraticStart[iColumn];
59         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
60      int jColumn = columnQuadratic[j];
61      listNonLinearColumn[jColumn]=1;
62      listNonLinearColumn[iColumn]=1;
63    }
64  }
65  for (iColumn=0;iColumn<numberColumns;iColumn++) {
66    if(listNonLinearColumn[iColumn])
67      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
68  }
69 
70  if (!numberNonLinearColumns) {
71    delete [] listNonLinearColumn;
72    // no quadratic part
73    return primal(0);
74  }
75
76  // get feasible
77  if (numberPrimalInfeasibilities())
78    primal(1);
79  // still infeasible
80  if (numberPrimalInfeasibilities())
81    return 0;
82
83  int jNon;
84  int * last[3];
85 
86  double * trust = new double[numberNonLinearColumns];
87  double * trueLower = new double[numberNonLinearColumns];
88  double * trueUpper = new double[numberNonLinearColumns];
89  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
90    iColumn=listNonLinearColumn[jNon];
91    trust[jNon]=0.5;
92    trueLower[jNon]=columnLower[iColumn];
93    trueUpper[jNon]=columnUpper[iColumn];
94  }
95  int iPass;
96  double lastObjective=1.0e31;
97  double * saveSolution = new double [numberColumns];
98  double targetDrop=1.0e31;
99  double objectiveOffset;
100  getDblParam(ClpObjOffset,objectiveOffset);
101  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
102  for (iPass=0;iPass<3;iPass++) {
103    last[iPass]=new int[numberNonLinearColumns];
104    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
105      last[iPass][jNon]=0;
106  }
107  double goodMove=false;
108  for (iPass=0;iPass<numberPasses;iPass++) {
109    // redo objective
110    double offset=0.0;
111    double objValue=-objectiveOffset;
112    memcpy(objective,saveObjective,numberColumns*sizeof(double));
113    for (iColumn=0;iColumn<numberColumns;iColumn++) 
114      objValue += objective[iColumn]*solution[iColumn];
115    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
116      iColumn=listNonLinearColumn[jNon];
117      double valueI = solution[iColumn];
118      int j;
119      for (j=columnQuadraticStart[iColumn];
120           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
121        int jColumn = columnQuadratic[j];
122        double valueJ = solution[jColumn];
123        double elementValue = quadraticElement[j];
124        objValue += 0.5*valueI*valueJ*elementValue;
125        offset += 0.5*valueI*valueJ*elementValue;
126        double gradientI = valueJ*elementValue;
127        double gradientJ = valueI*elementValue;
128        offset -= gradientI*valueI;
129        objective[iColumn] += gradientI;
130        offset -= gradientJ*valueJ;
131        objective[jColumn] += gradientJ;
132      }
133    }
134    printf("objective %g, objective offset %g\n",objValue,offset);
135    setDblParam(ClpObjOffset,objectiveOffset-offset);
136    objValue *= whichWay;
137    int * temp=last[2];
138    last[2]=last[1];
139    last[1]=last[0];
140    last[0]=temp;
141    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
142      iColumn=listNonLinearColumn[jNon];
143      double change = solution[iColumn]-saveSolution[iColumn];
144      if (change<-1.0e-5) {
145        if (fabs(change+trust[jNon])<1.0e-5) 
146          temp[jNon]=-1;
147        else
148          temp[jNon]=-2;
149      } else if(change>1.0e-5) {
150        if (fabs(change-trust[jNon])<1.0e-5) 
151          temp[jNon]=1;
152        else
153          temp[jNon]=2;
154      } else {
155        temp[jNon]=0;
156      }
157    } 
158    if (objValue<=lastObjective) 
159      goodMove=true;
160    else
161      goodMove=false;
162    double maxDelta=0.0;
163    double maxGap=0.0;
164    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
165      iColumn=listNonLinearColumn[jNon];
166      maxDelta = max(maxDelta,
167                     fabs(solution[iColumn]-saveSolution[iColumn]));
168      if (goodMove) {
169        if (last[0][jNon]*last[1][jNon]<0) {
170          // halve
171          trust[jNon] *= 0.5;
172        } else {
173          if (last[0][jNon]==last[1][jNon]&&
174              last[0][jNon]==last[2][jNon])
175            trust[jNon] *= 1.5; 
176        }
177      } else if (trust[jNon]>10.0*deltaTolerance) {
178        trust[jNon] *= 0.5;
179      }
180      maxGap = max(maxGap,trust[jNon]);
181    }
182    std::cout<<"largest gap is "<<maxGap<<std::endl;
183    if (goodMove) {
184      double drop = lastObjective-objValue;
185      std::cout<<"True drop was "<<drop<<std::endl;
186      std::cout<<"largest delta is "<<maxDelta<<std::endl;
187      if (maxDelta<deltaTolerance&&drop<1.0e-4) {
188        std::cout<<"Exiting"<<std::endl;
189        break;
190      }
191    } else {
192      lastObjective += 1.0e-3; // to stop exit
193    }
194    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
195      iColumn=listNonLinearColumn[jNon];
196      columnLower[iColumn]=max(solution[iColumn]
197                               -trust[jNon],
198                               trueLower[jNon]);
199      columnUpper[iColumn]=min(solution[iColumn]
200                               +trust[jNon],
201                               trueUpper[jNon]);
202    }
203    if (goodMove) {
204      memcpy(saveSolution,solution,numberColumns*sizeof(double));
205     
206      targetDrop=0.0;
207      if (iPass) {
208        // get reduced costs
209        this->matrix()->transposeTimes(this->dualRowSolution(),
210                                       this->dualColumnSolution());
211        double * r = this->dualColumnSolution();
212        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
213          iColumn=listNonLinearColumn[jNon];
214          double dj = objective[iColumn]-r[iColumn];
215          if (dj<0.0) 
216            targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
217          else
218            targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
219        }
220      }
221      std::cout<<"Pass - "<<iPass
222               <<", target drop is "<<targetDrop
223               <<std::endl;
224      lastObjective = objValue;
225      this->primal(1);
226    } else {
227      // bad pass - restore solution
228      memcpy(solution,saveSolution,numberColumns*sizeof(double));
229      iPass--;
230    }
231  }
232  // restore solution
233  memcpy(solution,saveSolution,numberColumns*sizeof(double));
234  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
235    iColumn=listNonLinearColumn[jNon];
236    columnLower[iColumn]=max(solution[iColumn],
237                             trueLower[jNon]);
238    columnUpper[iColumn]=min(solution[iColumn],
239                             trueUpper[jNon]);
240  }
241  // redo objective
242  double offset=0.0;
243  double objValue=-objectiveOffset;
244  memcpy(objective,saveObjective,numberColumns*sizeof(double));
245  for (iColumn=0;iColumn<numberColumns;iColumn++) 
246    objValue += objective[iColumn]*solution[iColumn];
247  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
248    iColumn=listNonLinearColumn[jNon];
249    double valueI = solution[iColumn];
250    int j;
251    for (j=columnQuadraticStart[iColumn];
252         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
253      int jColumn = columnQuadratic[j];
254      double valueJ = solution[jColumn];
255      double elementValue = quadraticElement[j];
256      objValue += 0.5*valueI*valueJ*elementValue;
257      offset += 0.5*valueI*valueJ*elementValue;
258      double gradientI = valueJ*elementValue;
259      double gradientJ = valueI*elementValue;
260      offset -= gradientI*valueI;
261      objective[iColumn] += gradientI;
262      offset -= gradientJ*valueJ;
263      objective[jColumn] += gradientJ;
264    }
265  }
266  printf("objective %g, objective offset %g\n",objValue,offset);
267  setDblParam(ClpObjOffset,objectiveOffset-offset);
268  this->primal(1);
269  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
270    iColumn=listNonLinearColumn[jNon];
271    columnLower[iColumn]= trueLower[jNon];
272    columnUpper[iColumn]= trueUpper[jNon];
273  }
274  delete [] saveSolution;
275  for (iPass=0;iPass<3;iPass++) 
276    delete [] last[iPass];
277  delete [] trust;
278  delete [] trueUpper;
279  delete [] trueLower;
280  delete [] saveObjective;
281  delete [] listNonLinearColumn;
282  return 0;
283}
284// Beale's method
285int 
286ClpSimplexPrimalQuadratic::primalBeale()
287{
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);
494  return 0;
495}
496 
497
Note: See TracBrowser for help on using the repository browser.