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

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

Trying

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.2 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
125  for (iPass=0;iPass<maxPass;iPass++) {
126    printf("Start of pass %d\n",iPass);
127    //printf("Bug until submodel new version\n");
128    CoinSort_2(sort,sort+numberSort,weight);
129    // Create small problem
130    ClpSimplex small(&model,numberRows,whichRows,numberSort,sort);
131    // now see what variables left out do to row solution
132    double * rowSolution = model.primalRowSolution();
133    memset (rowSolution,0,numberRows*sizeof(double));
134    int iRow,iColumn;
135    // zero out ones in small problem
136    for (iColumn=0;iColumn<numberSort;iColumn++) {
137      int kColumn = sort[iColumn];
138      fullSolution[kColumn]=0.0;
139    }
140    model.times(1.0,fullSolution,rowSolution);
141
142    double * lower = small.rowLower();
143    double * upper = small.rowUpper();
144    for (iRow=0;iRow<numberRows;iRow++) {
145      if (lower[iRow]>-1.0e50) 
146        lower[iRow] -= rowSolution[iRow];
147      if (upper[iRow]<1.0e50)
148        upper[iRow] -= rowSolution[iRow];
149    }
150    /* For some problems a useful variant is to presolve problem.
151       In this case you need to adjust smallNumberColumns to get
152       right size problem.  Also you can dispense with creating
153       small problem and fix variables in large problem and do presolve
154       on that. */
155    // Solve
156    small.primal();
157    // move solution back
158    const double * solution = small.primalColumnSolution();
159    for (iColumn=0;iColumn<numberSort;iColumn++) {
160      int kColumn = sort[iColumn];
161      model.setColumnStatus(kColumn,small.getColumnStatus(iColumn));
162      fullSolution[kColumn]=solution[iColumn];
163    }
164    for (iRow=0;iRow<numberRows;iRow++) 
165      model.setRowStatus(iRow,small.getRowStatus(iRow));
166    memcpy(model.primalRowSolution(),small.primalRowSolution(),
167           numberRows*sizeof(double));
168    if ((small.objectiveValue()>lastObjective-1.0e-7&&iPass>5)||
169        !small.numberIterations()||
170        iPass==maxPass-1) {
171
172      break; // finished
173    } else {
174      lastObjective = small.objectiveValue();
175      // get reduced cost for large problem
176      // this assumes minimization
177      memcpy(weight,model.objective(),numberColumns*sizeof(double));
178      model.transposeTimes(-1.0,small.dualRowSolution(),weight);
179      // now massage weight so all basic in plus good djs
180      for (iColumn=0;iColumn<numberColumns;iColumn++) {
181        double dj = weight[iColumn];
182        double value = fullSolution[iColumn];
183        if (model.getColumnStatus(iColumn)==ClpSimplex::basic)
184          dj = -1.0e50;
185        else if (dj<0.0&&value<columnUpper[iColumn])
186          dj = dj;
187        else if (dj>0.0&&value>columnLower[iColumn])
188          dj = -dj;
189        else if (columnUpper[iColumn]>columnLower[iColumn])
190          dj = fabs(dj);
191        else
192          dj = 1.0e50;
193        weight[iColumn] = dj;
194        sort[iColumn] = iColumn;
195      }
196      // sort
197      CoinSort_2(weight,weight+numberColumns,sort);
198      numberSort = smallNumberColumns;
199    }
200  }
201  if (addSlacks) {
202    int i;
203    int numberArtificials = numberColumns-originalNumberColumns;
204    for (i=0;i<numberArtificials;i++)
205      sort[i] = i + originalNumberColumns;
206    model.deleteColumns(numberArtificials,sort);
207  }
208  delete [] weight;
209  delete [] sort;
210  delete [] whichRows;
211  model.primal(1);
212  return 0;
213}   
Note: See TracBrowser for help on using the repository browser.