source: trunk/Samples/piece.cpp @ 159

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

Still not very good

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 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
7void sortSegments(int, int*, double*, double*);
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  numberColumns2=0;
77  numberElements=0;
78  start2[0]=0;
79  int segptr=0;
80  int numseg=0;
81
82  // treat first column separately
83  iColumn=0;
84  segstart[0]=0;
85  columnLower2[numberColumns2] = columnLower1[iColumn];
86  columnUpper2[numberColumns2] = columnUpper1[iColumn];
87  objective2[numberColumns2] = objective1[iColumn];
88  breakpt[segptr] = columnLower1[iColumn];
89  slope[segptr++] = objective1[iColumn];
90  double maxUpper = columnUpper1[iColumn];
91  numseg=1;
92  for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
93    row2[numberElements]=row1[j];
94    element2[numberElements++] = element1[j];
95  }
96  newSolution[0]=oldSolution[0];
97  start2[++numberColumns2]=numberElements;
98
99  // Now check for duplicates
100  for (iColumn=1;iColumn<numberColumns;iColumn++) {
101    maxUpper = max(maxUpper, columnUpper1[iColumn]);
102    // test if column identical to previous column
103    bool ifcopy=1;
104    int  joff = length1[iColumn-1];
105    for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
106      if (row1[j] != row1[j-joff]){
107        ifcopy=0;
108        break;
109      }
110      if (element1[j] != element1[j-joff]){
111        ifcopy=0;
112        break;
113      }
114    }
115    if (ifcopy) {
116      // subtract out from rhs
117      double fixed = columnLower1[iColumn];
118      if(fabs(fixed-columnUpper1[iColumn-1])>1.0e-8) {
119        // try other way
120        fixed = columnUpper1[iColumn];
121        assert(fabs(fixed-columnLower1[iColumn-1])<1.0e-8);
122      }
123      // do offset
124      objectiveOffset += fixed*objective1[iColumn];
125      for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
126        int iRow = row1[j];
127        double value = element1[j];
128        if (rowLower2[iRow]>-1.0e30)
129          rowLower2[iRow] -= value*fixed;
130        if (rowUpper2[iRow]<1.0e30)
131          rowUpper2[iRow] -= value*fixed;
132      }
133      if (fabs(oldSolution[iColumn]-fixed)>1.0e-8)
134        newSolution[numberColumns2-1]=oldSolution[iColumn]; 
135      breakpt[segptr] = columnLower1[iColumn];
136      slope[segptr++] = objective1[iColumn];
137      numseg++;
138      continue;
139    }
140    breakpt[segptr] = maxUpper;
141    slope[segptr++] = COIN_DBL_MAX;
142    segstart[numberColumns2] = segptr;
143    // new column found
144    columnLower2[numberColumns2] = columnLower1[iColumn];
145    columnUpper2[numberColumns2] = columnUpper1[iColumn];
146    objective2[numberColumns2] = objective1[iColumn];
147    breakpt[segptr] = columnLower1[iColumn];
148    slope[segptr++] = objective1[iColumn];
149    maxUpper = columnUpper1[iColumn];
150    for (int j=start1[iColumn];j<start1[iColumn]+length1[iColumn];j++) {
151      row2[numberElements]=row1[j];
152      element2[numberElements++] = element1[j];
153    }
154    newSolution[numberColumns2]=oldSolution[iColumn]; 
155    start2[++numberColumns2]=numberElements;
156  }
157  breakpt[segptr] = maxUpper;
158  slope[segptr++] = COIN_DBL_MAX;
159  segstart[numberColumns2] = segptr;
160
161  // print new number of columns, elements
162  printf("New number of columns  = %d\n",numberColumns2);
163  printf("New number of elements = %d\n",numberElements);
164  printf("Objective offset is %g\n",objectiveOffset);
165
166  /*  for (int k=0; k<20; k++)
167    printf("%d  %e %e\n",segstart[k],breakpt[k],slope[k]);
168  */
169  sortSegments(numberColumns2, segstart, breakpt, slope);
170
171  printf("After sorting\n");
172  /*  for (int k=0; k<20; k++)
173    printf("%d  %e %e\n",segstart[k],breakpt[k],slope[k]);
174  */
175
176  ClpSimplex  model; 
177
178  // load up
179  model.loadProblem(numberColumns2,numberRows,
180                    start2,row2,element2,
181                    columnLower2,columnUpper2,
182                    objective2,
183                    rowLower2,rowUpper2);
184  model.scaling(0);
185  // Create nonlinear objective
186  int returnCode= model.createPiecewiseLinearCosts(segstart, breakpt, slope);
187  assert (!returnCode);
188
189  // delete
190  delete [] segstart;
191  delete [] breakpt;
192  delete [] slope;
193  delete [] start2;
194  delete [] row2 ;
195  delete [] element2;
196
197  delete [] objective2;
198  delete [] columnLower2;
199  delete [] columnUpper2;
200  delete [] rowLower2;
201  delete [] rowUpper2;
202
203  // copy in solution - (should be optimal)
204  model.allSlackBasis();
205  memcpy(model.primalColumnSolution(),newSolution,numberColumns2*sizeof(double));
206  //memcpy(model.columnLower(),newSolution,numberColumns2*sizeof(double));
207  //memcpy(model.columnUpper(),newSolution,numberColumns2*sizeof(double));
208  delete [] newSolution;
209
210  const double * solution = model.primalColumnSolution();
211  for (iColumn=0;iColumn<numberColumns2;iColumn++)
212    printf("%g ",solution[iColumn]);
213  printf("\n");
214  // solve
215  model.primal(1);
216  for (iColumn=0;iColumn<numberColumns2;iColumn++)
217    printf("%g ",solution[iColumn]);
218  printf("\n");
219  model.primal(1);
220  for (iColumn=0;iColumn<numberColumns2;iColumn++)
221    printf("%g ",solution[iColumn]);
222  printf("\n");
223  model.primal();
224  for (iColumn=0;iColumn<numberColumns2;iColumn++)
225    printf("%g ",solution[iColumn]);
226  printf("\n");
227  model.allSlackBasis();
228  for (iColumn=0;iColumn<numberColumns2;iColumn++)
229    printf("%g ",solution[iColumn]);
230  printf("\n");
231  model.primal();
232  for (iColumn=0;iColumn<numberColumns2;iColumn++)
233    printf("%g ",solution[iColumn]);
234  printf("\n");
235  return 0;
236}   
237
238// Quick and dirty sort for 2 segments only
239void sortSegments(int n_, int* ptr_, double* breakpt_, double* slope_){
240  for (int i=0; i<n_; i++){
241    if ( (ptr_[i+1] - ptr_[i]) <= 2)
242      continue;
243    if ( (ptr_[i+1] - ptr_[i]) > 3){
244      printf("More than 2 real segments not implemented %d\n", i);
245      exit(2);
246    }
247    if (breakpt_[ptr_[i]] < breakpt_[ptr_[i]+1])
248      continue;
249    double bsave = breakpt_[ptr_[i]];
250    double ssave = slope_[ptr_[i]];
251    breakpt_[ptr_[i]] = breakpt_[ptr_[i]+1];
252    slope_[ptr_[i]] = slope_[ptr_[i]+1];
253    breakpt_[ptr_[i]+1] = bsave;
254    slope_[ptr_[i]+1] = ssave;
255  }
256  return;
257}
Note: See TracBrowser for help on using the repository browser.