source: branches/pre/Samples/sprint.cpp @ 220

Last change on this file since 220 was 220, checked in by forrest, 16 years ago

Stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.5 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpSimplex.hpp"
5#include "CoinSort.hpp"
6#include <iomanip>
7
8int main (int argc, const char *argv[])
9{
10  ClpSimplex  model;
11  int status;
12  // Keep names
13  if (argc<2) {
14    status=model.readMps("small.mps",true);
15  } else {
16    status=model.readMps(argv[1],true);
17  }
18  if (status)
19    exit(10);
20  /*
21    This driver implements what I called Sprint.  Cplex calls it
22    "sifting" which is just as silly.  When I thought of this trivial idea
23    it reminded me of an LP code of the 60's called sprint which after
24    every factorization took a subset of the matrix into memory (all
25    64K words!) and then iterated very fast on that subset.  On the
26    problems of those days it did not work very well, but it worked very
27    well on aircrew scheduling problems where there were very large numbers
28    of columns all with the same flavor.
29  */
30
31  /* The idea works best if you can get feasible easily.  To make it
32     more general we can add in costed slacks */
33
34  int originalNumberColumns = model.numberColumns();
35  int numberRows = model.numberRows();
36
37  // We will need arrays to choose variables.  These are too big but ..
38  double * weight = new double [numberRows+originalNumberColumns];
39  int * sort = new int [numberRows+originalNumberColumns];
40  int numberSort=0;
41  // Say we are going to add slacks - if you can get a feasible
42  // solution then do that at the comment - Add in your own coding here
43  bool addSlacks = true;
44
45  if (addSlacks) {
46    // initial list will just be artificials
47    // first we will set all variables as close to zero as possible
48    int iColumn;
49    const double * columnLower = model.columnLower();
50    const double * columnUpper = model.columnUpper();
51    double * columnSolution = model.primalColumnSolution();
52
53    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
54      double value =0.0;
55      if (columnLower[iColumn]>0.0)
56        value = columnLower[iColumn];
57      else if (columnUpper[iColumn]<0.0)
58        value = columnUpper[iColumn];
59      columnSolution[iColumn]=value;
60    }
61    // now see what that does to row solution
62    double * rowSolution = model.primalRowSolution();
63    memset (rowSolution,0,numberRows*sizeof(double));
64    model.times(1.0,columnSolution,rowSolution);
65
66    int * addStarts = new int [numberRows+1];
67    int * addRow = new int[numberRows];
68    double * addElement = new double[numberRows];
69    const double * lower = model.rowLower();
70    const double * upper = model.rowUpper();
71    addStarts[0]=0;
72    int numberArtificials=0;
73    double * addCost = new double [numberRows];
74    const double penalty=1.0e8;
75    int iRow;
76    for (iRow=0;iRow<numberRows;iRow++) {
77      if (lower[iRow]>rowSolution[iRow]) {
78        addRow[numberArtificials]=iRow;
79        addElement[numberArtificials]=1.0;
80        addCost[numberArtificials]=penalty;
81        numberArtificials++;
82        addStarts[numberArtificials]=numberArtificials;
83      } else if (upper[iRow]<rowSolution[iRow]) {
84        addRow[numberArtificials]=iRow;
85        addElement[numberArtificials]=-1.0;
86        addCost[numberArtificials]=penalty;
87        numberArtificials++;
88        addStarts[numberArtificials]=numberArtificials;
89      }
90    }
91    model.addColumns(numberArtificials,NULL,NULL,addCost,
92                      addStarts,addRow,addElement);
93    delete [] addStarts;
94    delete [] addRow;
95    delete [] addElement;
96    delete [] addCost;
97    // Set up initial list
98    numberSort=numberArtificials;
99    int i;
100    for (i=0;i<numberSort;i++)
101      sort[i] = i+originalNumberColumns;
102  } else {
103    // Get initial list in some magical way
104    // Add in your own coding here
105    abort();
106  }
107
108  int numberColumns = model.numberColumns();
109  const double * columnLower = model.columnLower();
110  const double * columnUpper = model.columnUpper();
111  double * fullSolution = model.primalColumnSolution();
112
113  // Just do this number of passes
114  int maxPass=100;
115  int iPass;
116  double lastObjective=1.0e31;
117
118  // Just take this number of columns in small problem
119  int smallNumberColumns = min(3*numberRows,numberColumns);
120  // We will be using all rows
121  int * whichRows = new int [numberRows];
122  for (int iRow=0;iRow<numberRows;iRow++)
123    whichRows[iRow]=iRow;
124  double originalOffset;
125  model2->getDblParam(ClpObjOffset,originalOffset);
126
127  for (iPass=0;iPass<maxPass;iPass++) {
128    printf("Start of pass %d\n",iPass);
129    //printf("Bug until submodel new version\n");
130    CoinSort_2(sort,sort+numberSort,weight);
131    // Create small problem
132    ClpSimplex small(&model,numberRows,whichRows,numberSort,sort);
133    // now see what variables left out do to row solution
134    double * rowSolution = model.primalRowSolution();
135    memset (rowSolution,0,numberRows*sizeof(double));
136    int iRow,iColumn;
137    // zero out ones in small problem
138    for (iColumn=0;iColumn<numberSort;iColumn++) {
139      int kColumn = sort[iColumn];
140      fullSolution[kColumn]=0.0;
141    }
142    // Get objective offset
143    double offset=0.0;
144    const double * objective = model2->objective();
145    for (iColumn=0;iColumn<numberColumns_;iColumn++) 
146      offset += fullSolution[iColumn]*objective[iColumn];
147    small.setDblParam(ClpObjOffset,originalOffset-offset);
148    model.times(1.0,fullSolution,rowSolution);
149
150    double * lower = small.rowLower();
151    double * upper = small.rowUpper();
152    for (iRow=0;iRow<numberRows;iRow++) {
153      if (lower[iRow]>-1.0e50) 
154        lower[iRow] -= rowSolution[iRow];
155      if (upper[iRow]<1.0e50)
156        upper[iRow] -= rowSolution[iRow];
157    }
158    /* For some problems a useful variant is to presolve problem.
159       In this case you need to adjust smallNumberColumns to get
160       right size problem.  Also you can dispense with creating
161       small problem and fix variables in large problem and do presolve
162       on that. */
163    // Solve
164    small.primal();
165    // move solution back
166    const double * solution = small.primalColumnSolution();
167    for (iColumn=0;iColumn<numberSort;iColumn++) {
168      int kColumn = sort[iColumn];
169      model.setColumnStatus(kColumn,small.getColumnStatus(iColumn));
170      fullSolution[kColumn]=solution[iColumn];
171    }
172    for (iRow=0;iRow<numberRows;iRow++) 
173      model.setRowStatus(iRow,small.getRowStatus(iRow));
174    memcpy(model.primalRowSolution(),small.primalRowSolution(),
175           numberRows*sizeof(double));
176    if ((small.objectiveValue()>lastObjective-1.0e-7&&iPass>5)||
177        !small.numberIterations()||
178        iPass==maxPass-1) {
179
180      break; // finished
181    } else {
182      lastObjective = small.objectiveValue();
183      // get reduced cost for large problem
184      // this assumes minimization
185      memcpy(weight,model.objective(),numberColumns*sizeof(double));
186      model.transposeTimes(-1.0,small.dualRowSolution(),weight);
187      // now massage weight so all basic in plus good djs
188      for (iColumn=0;iColumn<numberColumns;iColumn++) {
189        double dj = weight[iColumn];
190        double value = fullSolution[iColumn];
191        if (model.getColumnStatus(iColumn)==ClpSimplex::basic)
192          dj = -1.0e50;
193        else if (dj<0.0&&value<columnUpper[iColumn])
194          dj = dj;
195        else if (dj>0.0&&value>columnLower[iColumn])
196          dj = -dj;
197        else if (columnUpper[iColumn]>columnLower[iColumn])
198          dj = fabs(dj);
199        else
200          dj = 1.0e50;
201        weight[iColumn] = dj;
202        sort[iColumn] = iColumn;
203      }
204      // sort
205      CoinSort_2(weight,weight+numberColumns,sort);
206      numberSort = smallNumberColumns;
207    }
208  }
209  if (addSlacks) {
210    int i;
211    int numberArtificials = numberColumns-originalNumberColumns;
212    for (i=0;i<numberArtificials;i++)
213      sort[i] = i + originalNumberColumns;
214    model.deleteColumns(numberArtificials,sort);
215  }
216  delete [] weight;
217  delete [] sort;
218  delete [] whichRows;
219  model.primal(1);
220  return 0;
221}   
Note: See TracBrowser for help on using the repository browser.