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

Last change on this file since 1559 was 1559, checked in by stefan, 10 years ago

merge split branch into trunk

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