source: trunk/Clp/examples/piece.cpp

Last change on this file was 2278, checked in by forrest, 15 months ago

COIN_BIG_INDEX 2 changes

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