source: trunk/Clp/examples/piece.cpp @ 800

Last change on this file since 800 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: 8.1 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/* This simple example takes a matrix read in by CoinMpsIo,
5   deletes every second column and solves the resulting problem */
6
7#include "ClpSimplex.hpp"
8#include "ClpNonLinearCost.hpp"
9#include "CoinMpsIO.hpp"
10#include <iomanip>
11
12int main (int argc, const char *argv[])
13{
14  int status;
15  CoinMpsIO m;
16  if (argc<2)
17    status=m.readMps("model1.mps","");
18  else
19    status=m.readMps(argv[1],"");
20
21  if (status) {
22    fprintf(stdout,"Bad readMps %s\n",argv[1]);
23    exit(1);
24  }
25
26  // Load up model1 - so we can use known good solution
27  ClpSimplex model1;
28  model1.loadProblem(*m.getMatrixByCol(),
29                    m.getColLower(),m.getColUpper(),
30                    m.getObjCoefficients(),
31                    m.getRowLower(),m.getRowUpper());
32  model1.dual();
33  // Get data arrays
34  const CoinPackedMatrix * matrix1 = m.getMatrixByCol();
35  const int * start1 = matrix1->getVectorStarts();
36  const int * length1 = matrix1->getVectorLengths();
37  const int * row1 = matrix1->getIndices();
38  const double * element1 = matrix1->getElements();
39
40  const double * columnLower1 = m.getColLower();
41  const double * columnUpper1 = m.getColUpper();
42  const double * rowLower1 = m.getRowLower();
43  const double * rowUpper1 = m.getRowUpper();
44  const double * objective1 = m.getObjCoefficients();
45
46  int numberColumns = m.getNumCols();
47  int numberRows = m.getNumRows();
48  int numberElements = m.getNumElements();
49
50  // Get new arrays
51  int numberColumns2 = (numberColumns+1);
52  int * start2 = new int[numberColumns2+1];
53  int * row2 = new int[numberElements];
54  double * element2 = new double[numberElements];
55  int * segstart = new int[numberColumns+1];
56  double * breakpt  = new double[2*numberColumns];
57  double * slope  = new double[2*numberColumns];
58
59  double * objective2 = new double[numberColumns2];
60  double * columnLower2 = new double[numberColumns2];
61  double * columnUpper2 = new double[numberColumns2];
62  double * rowLower2 = new double[numberRows];
63  double * rowUpper2 = new double[numberRows];
64
65  // We need to modify rhs
66  memcpy(rowLower2,rowLower1,numberRows*sizeof(double));
67  memcpy(rowUpper2,rowUpper1,numberRows*sizeof(double));
68  double objectiveOffset = 0.0;
69
70  // For new solution
71  double * newSolution = new double [numberColumns];
72  const double * oldSolution = model1.primalColumnSolution();
73
74  int iColumn;
75  for (iColumn=0;iColumn<numberColumns;iColumn++)
76    printf("%g ",oldSolution[iColumn]);
77  printf("\n");
78   
79  numberColumns2=0;
80  numberElements=0;
81  start2[0]=0;
82  int segptr=0;
83
84  segstart[0]=0;
85
86  // Now check for duplicates
87  for (iColumn=0;iColumn<numberColumns;iColumn++) {
88    // test if column identical to next column
89    bool ifcopy=1;
90    if (iColumn<numberColumns-1) {
91      int  joff = length1[iColumn];
92      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
93        if (row1[j] != row1[j+joff]){
94          ifcopy=0;
95          break;
96        }
97        if (element1[j] != element1[j+joff]){
98          ifcopy=0;
99          break;
100        }
101      }
102    } else {
103      ifcopy=0;
104    }
105    //if (iColumn>47||iColumn<45)
106    //ifcopy=0;
107    if (ifcopy) {
108      double lo1 = columnLower1[iColumn];
109      double up1 = columnUpper1[iColumn];
110      double obj1 = objective1[iColumn];
111      double sol1 = oldSolution[iColumn];
112      double lo2 = columnLower1[iColumn+1];
113      double up2 = columnUpper1[iColumn+1];
114      double obj2 = objective1[iColumn+1];
115      double sol2 = oldSolution[iColumn+1];
116      if(fabs(up1-lo2)>1.0e-8) {
117        // try other way
118        double temp;
119        temp = lo1;
120        lo1=lo2;
121        lo2=temp;
122        temp = up1;
123        up1=up2;
124        up2=temp;
125        temp = obj1;
126        obj1=obj2;
127        obj2=temp;
128        temp = sol1;
129        sol1=sol2;
130        sol2=temp;
131        assert(fabs(up1-lo2)<1.0e-8);
132      }
133      // subtract out from rhs
134      double fixed = up1;
135      // do offset
136      objectiveOffset += fixed*obj2;
137      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
138        int iRow = row1[j];
139        double value = element1[j];
140        if (rowLower2[iRow]>-1.0e30)
141          rowLower2[iRow] -= value*fixed;
142        if (rowUpper2[iRow]<1.0e30)
143          rowUpper2[iRow] -= value*fixed;
144      }
145      newSolution[numberColumns2]=fixed;
146      if (fabs(sol1-fixed)>1.0e-8)
147        newSolution[numberColumns2]=sol1;
148      if (fabs(sol2-fixed)>1.0e-8)
149        newSolution[numberColumns2]=sol2;
150      columnLower2[numberColumns2] = lo1;
151      columnUpper2[numberColumns2] = up2;
152      objective2[numberColumns2] = 0.0;
153      breakpt[segptr] = lo1;
154      slope[segptr++] = obj1;
155      breakpt[segptr] = lo2;
156      slope[segptr++] = obj2;
157      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
158        row2[numberElements]=row1[j];
159        element2[numberElements++] = element1[j];
160      }
161      start2[++numberColumns2]=numberElements;
162      breakpt[segptr] = up2;
163      slope[segptr++] = COIN_DBL_MAX;
164      segstart[numberColumns2] = segptr;
165      iColumn++; // skip next column
166    } else {
167      // normal column
168      columnLower2[numberColumns2] = columnLower1[iColumn];
169      columnUpper2[numberColumns2] = columnUpper1[iColumn];
170      objective2[numberColumns2] = objective1[iColumn];
171      breakpt[segptr] = columnLower1[iColumn];
172      slope[segptr++] = objective1[iColumn];
173      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
174        row2[numberElements]=row1[j];
175        element2[numberElements++] = element1[j];
176      }
177      newSolution[numberColumns2]=oldSolution[iColumn]; 
178      start2[++numberColumns2]=numberElements;
179      breakpt[segptr] = columnUpper1[iColumn];
180      slope[segptr++] = COIN_DBL_MAX;
181      segstart[numberColumns2] = segptr;
182    }
183  }
184
185  // print new number of columns, elements
186  printf("New number of columns  = %d\n",numberColumns2);
187  printf("New number of elements = %d\n",numberElements);
188  printf("Objective offset is %g\n",objectiveOffset);
189
190
191  ClpSimplex  model; 
192
193  // load up
194  model.loadProblem(numberColumns2,numberRows,
195                    start2,row2,element2,
196                    columnLower2,columnUpper2,
197                    objective2,
198                    rowLower2,rowUpper2);
199  model.scaling(0);
200  model.setDblParam(ClpObjOffset,-objectiveOffset);
201  // Create nonlinear objective
202  int returnCode= model.createPiecewiseLinearCosts(segstart, breakpt, slope);
203  assert (!returnCode);
204
205  // delete
206  delete [] segstart;
207  delete [] breakpt;
208  delete [] slope;
209  delete [] start2;
210  delete [] row2 ;
211  delete [] element2;
212
213  delete [] objective2;
214  delete [] columnLower2;
215  delete [] columnUpper2;
216  delete [] rowLower2;
217  delete [] rowUpper2;
218
219  // copy in solution - (should be optimal)
220  model.allSlackBasis();
221  memcpy(model.primalColumnSolution(),newSolution,numberColumns2*sizeof(double));
222  //memcpy(model.columnLower(),newSolution,numberColumns2*sizeof(double));
223  //memcpy(model.columnUpper(),newSolution,numberColumns2*sizeof(double));
224  delete [] newSolution;
225  //model.setLogLevel(63);
226
227  const double * solution = model.primalColumnSolution();
228  double * saveSol = new double[numberColumns2];
229  memcpy(saveSol,solution,numberColumns2*sizeof(double));
230  for (iColumn=0;iColumn<numberColumns2;iColumn++)
231    printf("%g ",solution[iColumn]);
232  printf("\n");
233  // solve
234  model.primal(1);
235  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
236    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
237      printf(" ** was %g ",saveSol[iColumn]);
238    printf("%g ",solution[iColumn]);
239  }
240  printf("\n");
241  model.primal(1);
242  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
243    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
244      printf(" ** was %g ",saveSol[iColumn]);
245    printf("%g ",solution[iColumn]);
246  }
247  printf("\n");
248  model.primal();
249  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
250    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
251      printf(" ** was %g ",saveSol[iColumn]);
252    printf("%g ",solution[iColumn]);
253  }
254  printf("\n");
255  model.allSlackBasis();
256  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
257    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
258      printf(" ** was %g ",saveSol[iColumn]);
259    printf("%g ",solution[iColumn]);
260  }
261  printf("\n");
262  model.setLogLevel(63);
263  model.primal();
264  for (iColumn=0;iColumn<numberColumns2;iColumn++) {
265    if (fabs(solution[iColumn]-saveSol[iColumn])>1.0e-3)
266      printf(" ** was %g ",saveSol[iColumn]);
267    printf("%g ",solution[iColumn]);
268  }
269  printf("\n");
270  return 0;
271}   
Note: See TracBrowser for help on using the repository browser.