source: trunk/Samples/useVolume.cpp @ 225

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

This should break everything

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.4 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpSimplex.hpp"
5#include "ClpFactorization.hpp"
6#include "VolVolume.hpp"
7
8//#############################################################################
9
10class lpHook : public VOL_user_hooks {
11private:
12   lpHook(const lpHook&);
13   lpHook& operator=(const lpHook&);
14private:
15   /// Pointer to dense vector of structural variable upper bounds
16   double  *colupper_;
17   /// Pointer to dense vector of structural variable lower bounds
18   double  *collower_;
19   /// Pointer to dense vector of objective coefficients
20   double  *objcoeffs_;
21   /// Pointer to dense vector of right hand sides
22   double  *rhs_;
23   /// Pointer to dense vector of senses
24   char    *sense_;
25
26   /// The problem matrix in a row ordered form
27   CoinPackedMatrix rowMatrix_;
28   /// The problem matrix in a column ordered form
29   CoinPackedMatrix colMatrix_;
30
31public:
32   lpHook(const double* clb, const double* cub, const double* obj,
33          const double* rhs, const char* sense, const CoinPackedMatrix& mat);
34   virtual ~lpHook();
35   
36public:
37   // for all hooks: return value of -1 means that volume should quit
38   /** compute reduced costs   
39       @param u (IN) the dual variables
40       @param rc (OUT) the reduced cost with respect to the dual values
41   */
42   virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc);
43
44   /** Solve the subproblem for the subgradient step.
45       @param dual (IN) the dual variables
46       @param rc (IN) the reduced cost with respect to the dual values
47       @param lcost (OUT) the lagrangean cost with respect to the dual values
48       @param x (OUT) the primal result of solving the subproblem
49       @param v (OUT) b-Ax for the relaxed constraints
50       @param pcost (OUT) the primal objective value of <code>x</code>
51   */
52   virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
53                                double& lcost, VOL_dvector& x, VOL_dvector& v,
54                                double& pcost);
55   /** Starting from the primal vector x, run a heuristic to produce
56       an integer solution 
57       @param x (IN) the primal vector
58       @param heur_val (OUT) the value of the integer solution (return
59       <code>DBL_MAX</code> here if no feas sol was found
60   */
61   virtual int heuristics(const VOL_problem& p, 
62                          const VOL_dvector& x, double& heur_val) {
63      return 0;
64   }
65};
66 
67//#############################################################################
68
69lpHook::lpHook(const double* clb, const double* cub, const double* obj,
70               const double* rhs, const char* sense,
71               const CoinPackedMatrix& mat)
72{
73   const int colnum = mat.getNumCols();
74   const int rownum = mat.getNumRows();
75
76   colupper_ = new double[colnum];
77   collower_ = new double[colnum];
78   objcoeffs_ = new double[colnum];
79   rhs_ = new double[rownum];
80   sense_ = new char[rownum];
81
82   std::copy(clb, clb + colnum, collower_);
83   std::copy(cub, cub + colnum, colupper_);
84   std::copy(obj, obj + colnum, objcoeffs_);
85   std::copy(rhs, rhs + rownum, rhs_);
86   std::copy(sense, sense + rownum, sense_);
87
88   if (mat.isColOrdered()) {
89      colMatrix_.copyOf(mat);
90      rowMatrix_.reverseOrderedCopyOf(mat);
91   } else {
92      rowMatrix_.copyOf(mat);
93      colMatrix_.reverseOrderedCopyOf(mat);
94   }
95}
96
97//-----------------------------------------------------------------------------
98
99lpHook::~lpHook()
100{
101   delete[] colupper_;
102   delete[] collower_;
103   delete[] objcoeffs_;
104   delete[] rhs_;
105   delete[] sense_;
106}
107
108//#############################################################################
109
110int
111lpHook::compute_rc(const VOL_dvector& u, VOL_dvector& rc)
112{
113   rowMatrix_.transposeTimes(u.v, rc.v);
114   const int psize = rowMatrix_.getNumCols();
115
116   for (int i = 0; i < psize; ++i)
117      rc[i] = objcoeffs_[i] - rc[i];
118   return 0;
119}
120
121//-----------------------------------------------------------------------------
122
123int
124lpHook::solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
125                         double& lcost, VOL_dvector& x, VOL_dvector& v,
126                         double& pcost)
127{
128   int i;
129   const int psize = x.size();
130   const int dsize = v.size();
131
132   // compute the lagrangean solution corresponding to the reduced costs
133   for (i = 0; i < psize; ++i)
134      x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];
135
136   // compute the lagrangean value (rhs*dual + primal*rc)
137   lcost = 0;
138   for (i = 0; i < dsize; ++i)
139      lcost += rhs_[i] * dual[i];
140   for (i = 0; i < psize; ++i)
141      lcost += x[i] * rc[i];
142
143   // compute the rhs - lhs
144   colMatrix_.times(x.v, v.v);
145   for (i = 0; i < dsize; ++i)
146      v[i] = rhs_[i] - v[i];
147
148   // compute the lagrangean primal objective
149   pcost = 0;
150   for (i = 0; i < psize; ++i)
151      pcost += x[i] * objcoeffs_[i];
152
153   return 0;
154}
155
156//#############################################################################
157
158int main (int argc, const char *argv[])
159{
160  ClpSimplex  model;
161  int status;
162
163  if (argc<2)
164    status=model.readMps("../../Mps/Sample/p0033.mps",true);
165  else
166    status=model.readMps(argv[1],true);
167  /*
168    This driver uses volume algorithm
169    then does dual - after adjusting costs
170    then solves real problem
171  */
172
173  // do volume for a bit
174  VOL_problem volprob;
175  const CoinPackedMatrix* mat = model.matrix();
176  const int psize = mat->getNumCols();
177  const int dsize = mat->getNumRows();
178  char * sense = new char[dsize];
179  double * rhs = new double[dsize];
180  const double * rowLower = model.rowLower();
181  const double * rowUpper = model.rowUpper();
182  // Set the lb/ub on the duals
183  volprob.dsize = dsize;
184  volprob.psize = psize;
185  volprob.dual_lb.allocate(dsize);
186  volprob.dual_ub.allocate(dsize);
187  volprob.dsol.allocate(dsize);
188  int i;
189  for (i = 0; i < dsize; ++i) {
190    if (rowUpper[i]==rowLower[i]) {
191    // 'E':
192      volprob.dual_lb[i] = -1.0e31;
193      volprob.dual_ub[i] = 1.0e31;
194      rhs[i]=rowUpper[i];
195      sense[i]='E';
196    } else if (rowLower[i]<-0.99e10&&rowUpper[i]<0.99e10) {
197      // 'L':
198      volprob.dual_lb[i] = -1.0e31;
199      volprob.dual_ub[i] = 0.0;
200      rhs[i]=rowUpper[i];
201      sense[i]='L';
202    } else if (rowLower[i]>-0.99e10&&rowUpper[i]>0.99e10) {
203      // 'G':
204      volprob.dual_lb[i] = 0.0;
205      volprob.dual_ub[i] = 1.0e31;
206      rhs[i]=rowLower[i];
207      sense[i]='G';
208    } else {
209      printf("Volume Algorithm can't work if there is a non ELG row\n");
210      abort();
211    }
212  }
213  // Can't use read_param as private
214  // anyway I want automatic use - so maybe this is problem
215#if 0
216  FILE* infile = fopen("parameters", "r");
217  if (!infile) {
218    printf("Failure to open parameter file\n");
219  } else {
220    volprob.read_params("parameters");
221  }
222#endif
223#if 0
224  // should save and restore bounds
225  model.tightenPrimalBounds();
226#else
227  double * colUpper = model.columnUpper();
228  for (i=0;i<psize;i++)
229    colUpper[i]=1.0;
230#endif
231  lpHook myHook(model.getColLower(), model.getColUpper(),
232                model.getObjCoefficients(),
233                rhs, sense, *mat);
234  // move duals
235  double * pi = model.dualRowSolution();
236  memcpy(volprob.dsol.v,pi,dsize*sizeof(double));
237  volprob.solve(myHook,  false /* not warmstart */);
238  // For now stop as not doing any good
239  exit(77);
240  // create objectives
241  int numberRows = model.numberRows();
242  int numberColumns = model.numberColumns();
243  memcpy(pi,volprob.dsol.v,numberRows*sizeof(double));
244#define MODIFYCOSTS
245#ifdef MODIFYCOSTS
246  double * saveObj = new double[numberColumns];
247  memcpy(saveObj,model.objective(),numberColumns*sizeof(double));
248  memcpy(model.dualColumnSolution(),model.objective(),
249         numberColumns*sizeof(double));
250  model.clpMatrix()->transposeTimes(-1.0,pi,model.dualColumnSolution());
251  memcpy(model.objective(),model.dualColumnSolution(),
252         numberColumns*sizeof(double));
253  const double * rowsol = model.primalRowSolution();
254  //const double * rowLower = model.rowLower();
255  //const double * rowUpper = model.rowUpper();
256  double offset=0.0;
257  for (i=0;i<numberRows;i++) {
258    offset += pi[i]*rowsol[i];
259  }
260  double value2;
261  model.getDblParam(ClpObjOffset,value2);
262  printf("Offset %g %g\n",offset,value2);
263  model.setRowObjective(pi);
264  // zero out pi
265  memset(pi,0,numberRows*sizeof(double));
266#endif
267  // Could put some in basis - only partially tested
268  model.allSlackBasis(); 
269  model.factorization()->maximumPivots(1000);
270  //model.setLogLevel(63);
271  // solve
272  model.dual(1);
273  //model.primal(1);
274#ifdef MODIFYCOSTS
275  memcpy(model.objective(),saveObj,numberColumns*sizeof(double));
276  // zero out pi
277  memset(pi,0,numberRows*sizeof(double));
278  model.setRowObjective(pi);
279  delete [] saveObj;
280  model.primal();
281#endif
282 
283  return 0;
284}   
Note: See TracBrowser for help on using the repository browser.