source: branches/sandbox/Cbc/examples/cbc_driverC_sos.c @ 1273

Last change on this file since 1273 was 713, checked in by jpfasano, 12 years ago

Changes to get Cbc/examples/cbc_driverC_sos.c to build and run

File size: 15.2 KB
Line 
1/* Copyright (C) 2004-2007 EPRI
2Corporation and others.  All Rights Reserved. */
3
4/* This example shows the use of the "C" interface for CBC. */
5
6#include "Cbc_C_Interface.h"
7#include <stdio.h>
8#include <stdlib.h>
9#include <string.h>
10#include <assert.h>
11
12// prototypes
13void printSolution(Cbc_Model *cbc_model);
14Cbc_Model * getDefaultModel(int argc, const char *argv[]);
15
16
17/* Call back function - just says whenever it gets Clp0005 or Coin0002 */
18// TODO: It seems that Cbc gives callbacks but not Coin
19static void callBack(Cbc_Model * model, int messageNumber,
20                     int nDouble, const double * vDouble,
21                     int nInt, const int * vInt,
22                     int nString, char ** vString) 
23{
24  const char prefix[] = "cbc_driverC_sos.cpp::callBack(): ";
25  const int  VERBOSE = 0;
26  if (VERBOSE>0) printf("%s begin\n",prefix);
27  if (VERBOSE>1) printf("%s messageNumber %i\n",prefix,messageNumber);
28
29  if (messageNumber==1000002) {
30    /* Coin0002 */
31    assert (nString==1&&nInt==3);
32    printf("Name of problem is %s\n",vString[0]);
33    printf("row %d col %d el %d\n",vInt[0],vInt[1],vInt[2]);
34  } else if (messageNumber==5) {
35    /* Clp0005 */
36    int i;
37    assert (nInt==4&&nDouble==3); /* they may not all print */
38    for (i=0;i<3;i++)
39      printf("%d %g\n",vInt[i],vDouble[i]);
40  }
41
42  if (VERBOSE>0) printf("%s return\n",prefix);
43}
44
45/**
46* Get default model
47*/
48Cbc_Model * getDefaultModel(int argc, const char *argv[])
49{
50  const char prefix[] = "cbc_driverC_sos.cpp::getDefaultModel(): ";
51  const int  VERBOSE = 0;
52  Cbc_Model  *model;
53  int status; 
54
55  if (VERBOSE>0) printf("%s begin\n",prefix);
56  model = Cbc_newModel(); 
57
58  /** Amount of print out:
59  0 - none
60  1 - just final
61  2 - just factorizations
62  3 - as 2 plus a bit more
63  4 - verbose
64  above that 8,16,32 etc just for selective debug
65  */
66  Cbc_setLogLevel(model, 1);
67  if (VERBOSE>0) printf("%s Log Level %i\n", prefix, Cbc_logLevel(model));
68  // register callback
69  Cbc_registerCallBack(model,callBack);
70  // Keep names when reading an mps file
71  if (argc<2)
72    // status=Cbc_readMps(model,"../../../../../Data/Sample/ltw.mps");
73    status=Cbc_readMps(model,"../../../../../Data/Sample/p0033.mps");
74  else
75    status=Cbc_readMps(model,argv[1]);
76
77  if (status) {
78    fprintf(stderr,"Bad readMps %s\n",argv[1]);
79    fprintf(stdout,"Bad readMps %s\n",argv[1]);
80    exit(1);
81  }
82  Cbc_setOptimizationDirection(model, 1);
83  if (1) Cbc_setIntegerTolerance(model,  1.0e-5);
84
85  /// Solve initial LP relaxation
86  Cbc_initialSolve(model);
87
88  if (VERBOSE>0) printf("%s return\n",prefix);
89  return model;
90} //  getDefaultModel()
91
92void printSolution(Cbc_Model *cbc_model) {
93  //
94  //  Now to print out solution.  The methods used return modifiable
95  //  arrays while the alternative names return const pointers -
96  //  which is of course much more virtuous.
97  // 
98  //  This version just does non-zero column and row values.
99  //
100  //
101
102  // If we have not kept names (parameter to readMps) this will be 0
103  //   assert(Cbc_lengthNames(cbc_model));
104  {
105    int  name_length = 256;
106    char model_name[256];
107    Cbc_problemName(cbc_model, name_length, model_name);
108    printf("Model Name = '%s'\n", model_name);
109  }
110  printf("Iteration Count = %i\n",Cbc_getIterationCount(cbc_model));
111  printf("Iteration Limit = %i\n",Cbc_maximumIterations(cbc_model));
112  printf("Is Abandoned = %i\n",Cbc_isAbandoned(cbc_model));
113  printf("Is Proven Optimal = %i\n",Cbc_isProvenOptimal(cbc_model));
114  printf("Is Proven Infeasible = %i\n",Cbc_isProvenPrimalInfeasible(cbc_model));
115  printf("Is Proven Dual Infeasible = %i\n",Cbc_isProvenDualInfeasible(cbc_model));
116  printf("Is Proven Unbounded = %i\n",(Cbc_infeasibilityRay(cbc_model) == NULL) ? 0 : 1);
117  printf("Is Primal Objective Limit Reached = %i\n",Cbc_isPrimalObjectiveLimitReached(cbc_model));
118  printf("Is Dual Objective Limit Reached = %i\n",Cbc_isDualObjectiveLimitReached(cbc_model));
119  printf("Is Iteration Limit Reached = %i\n",Cbc_isIterationLimitReached(cbc_model));
120  printf("Objective Sense = %g\n",Cbc_getObjSense(cbc_model));  // (1 - minimize, -1 - maximize, 0 - ignore)
121  printf("Primal Feasible = %i\n",Cbc_primalFeasible(cbc_model));
122  printf("Dual Feasible = %i\n",Cbc_dualFeasible(cbc_model));
123  printf("Dual Bound = %g\n",Cbc_dualBound(cbc_model));
124  printf("Infeasibility Cost = %g\n",Cbc_infeasibilityCost(cbc_model));
125  printf("Sum Dual Infeasibilities = %g\n",Cbc_sumDualInfeasibilities(cbc_model));
126  printf("Number Dual Infeasibilities = %i\n",Cbc_numberDualInfeasibilities(cbc_model));
127  printf("Sum Primal Infeasibilities = %g\n",Cbc_sumPrimalInfeasibilities(cbc_model));
128  printf("Number Primal Infeasibilities = %i\n",Cbc_numberPrimalInfeasibilities(cbc_model));
129  printf("Objective Value = %g\n",Cbc_objectiveValue(cbc_model)); 
130  printf("Optimization Direction = %g\n", Cbc_optimizationDirection(cbc_model));
131  printf("  (1 - minimize, -1 - maximize, 0 - ignore)\n");
132  printf("--------------------------------------\n");
133
134  //  Rows
135  {
136    int numberRows = Cbc_numberRows(cbc_model);
137    int iRow;
138
139    const double * rowPrimal = Cbc_getRowActivity(cbc_model);
140    const double * rowDual   = Cbc_getRowPrice(cbc_model);
141    const double * rowLower  = Cbc_getRowLower(cbc_model);
142    const double * rowUpper  = Cbc_getRowUpper(cbc_model);
143
144    assert(rowPrimal != NULL);
145    assert(rowDual   != NULL);
146    assert(rowLower  != NULL);
147    assert(rowUpper  != NULL);
148
149    printf("                       Primal          Dual         Lower         Upper\n");
150    for (iRow=0;iRow<numberRows;iRow++) {
151      double value;
152      value = rowDual[iRow];
153      if (value>1.0e-8||value<-1.0e-8) {
154        char name[20];
155        sprintf(name," Row%-4i",iRow);
156        printf("%6d %8s",iRow,name);
157        printf(" %13g",rowPrimal[iRow]);
158        printf(" %13g",rowDual[iRow]);
159        printf(" %13g",rowLower[iRow]);
160        printf(" %13g",rowUpper[iRow]);
161        printf("\n");
162      }
163    }
164  }
165  printf("--------------------------------------\n");
166  // Columns
167  {
168    int numberColumns = Cbc_numberColumns(cbc_model);
169    int iColumn;
170
171    const double * columnPrimal    = Cbc_getColSolution(cbc_model);
172    const double * columnDual      = Cbc_getReducedCost(cbc_model);
173    const double * columnLower     = Cbc_getColLower(cbc_model);
174    const double * columnUpper     = Cbc_getColUpper(cbc_model);
175    const double * columnObjective = Cbc_getObjCoefficients(cbc_model);
176
177    assert(columnPrimal    != NULL);
178    assert(columnDual      != NULL);
179    assert(columnLower     != NULL);
180    assert(columnUpper     != NULL);
181    assert(columnObjective != NULL);
182
183    printf("                       Primal          Dual         Lower         Upper          Cost\n");
184    for (iColumn=0;iColumn<numberColumns;iColumn++) {
185      double value;
186      value = columnPrimal[iColumn];
187      if (value>1.0e-8||value<-1.0e-8) {
188        char name[20];
189        sprintf(name," Col%-4i",iColumn);
190        //      Cbc_columnName(cbc_model,iColumn,name);
191        printf("%6d %8s",iColumn,name);
192        printf(" %13g",columnPrimal[iColumn]);
193        printf(" %13g",columnDual[iColumn]);
194        printf(" %13g",columnLower[iColumn]);
195        printf(" %13g",columnUpper[iColumn]);
196        printf(" %13g",columnObjective[iColumn]);
197        printf("\n");
198      }
199    }
200  }
201  printf("--------------------------------------\n");
202} //  printSolution()
203
204int main (int argc, const char *argv[])
205{
206  const char prefix[] = "cbc_driverC_sos.cpp:main(): ";
207  const int  VERBOSE = 0;
208  Cbc_Model * model, * model2;
209  double time1;
210  char modelName[80];
211  int numberIntegers;
212  int * integerVariable;
213
214  if (VERBOSE>0) printf("%s begin\n",prefix);
215  if (VERBOSE>0) printf("%s Version %g\n",prefix,Cbc_getVersion());
216
217  // set model using the local routine for reading an MPS file
218  model = getDefaultModel(argc, argv);
219  model2 = NULL;  // used to keep model around
220
221  // This clause ought to set the initial basis, but does not yet work.
222  {
223    int i;
224    int row_length = Cbc_getNumRows(model);
225    int col_length = Cbc_getNumCols(model);
226    int rim_length = row_length + col_length;
227    int elem_length = Cbc_getNumElements(model);
228
229    int * cbc_rowStatus    = NULL; //(int *) malloc(row_length*sizeof(int));
230    int * cbc_columnStatus = NULL; //(int *) malloc(col_length*sizeof(int));
231    //    Cbc_getBasisStatus(model, cbc_columnStatus, cbc_rowStatus);
232
233    if (0) {
234      fprintf(stdout,"%s row_length = %i\n", prefix, row_length);
235      fprintf(stdout,"%s col_length = %i\n", prefix, col_length);
236      fprintf(stdout,"%s rim_length = %i\n", prefix, rim_length);
237      fprintf(stdout,"%s elem_length = %i\n", prefix, elem_length);
238      fflush(stdout);
239    }
240    // print solution status variables
241    if (0) {
242      if (cbc_rowStatus) {
243        for (i = 0; i < row_length; i++) {
244          printf("%s cbc_rowStatus[%i] = %i\n", prefix, i, cbc_rowStatus[i]);
245          fflush(stdout);
246        }
247      } else {
248        fprintf(stdout,"%s cbc_rowStatus = %p\n", prefix, cbc_rowStatus);
249        fflush(stdout);
250      }
251      if (cbc_columnStatus) {
252        for (i = 0; i < row_length; i++) {
253          fprintf(stdout,"%s cbc_rowStatus[%i] = %i\n", prefix, i, cbc_columnStatus[i]);
254          fflush(stdout);
255        }
256      } else {
257        fprintf(stdout,"%s cbc_columnStatus = %p\n", prefix, cbc_columnStatus);
258        fflush(stdout);
259      }
260    }
261    //    free( cbc_rowStatus );
262    //    free( cbc_columnStatus );
263  }
264
265  // Save model as a clone (does not work as of 2004).
266  model2 = Cbc_clone(model);
267
268  // Store which variables are integer as defined in the MPS file.
269  // They will be set to Continuous before adding the SOS constraints.
270
271  {
272    int i;
273    int numberColumns=Cbc_getNumCols(model);
274    numberIntegers = 0;
275    integerVariable = malloc(numberColumns*sizeof(int[1]));
276    for ( i=0;i<numberColumns;i++) {
277      if (Cbc_isInteger(model,i)) {
278        integerVariable[numberIntegers++]=i;
279        if (VERBOSE>3) printf("%s integerVariable[%i] = %i\n",prefix,
280          numberIntegers-1, integerVariable[numberIntegers-1]);
281      }
282    }
283  }
284  Cbc_problemName(model,80,modelName);
285
286  // Solve the MPS version of the problem
287  if (1) { 
288    if (VERBOSE>1) {
289      printf("%s Solve MPS version of the problem\n",prefix);
290      printf("%s Optimization Direction = %g (1 - minimize, -1 - maximize, 0 - ignore)\n",
291        prefix, Cbc_optimizationDirection(model)); 
292    }
293    //    Cbc_setLogLevel(model, VERBOSE); // 4 is verbose
294    time1 = Cbc_cpuTime(model) ;
295    Cbc_branchAndBound(model);
296    if (VERBOSE>1) printf("Model %s has %d rows and %d columns\n",
297      modelName,Cbc_numberRows(model),Cbc_numberColumns(model));
298
299    if (VERBOSE>1) printf("%s Solving model %s took %g seconds, %i nodes with objective %g\n",
300      prefix, modelName, (Cbc_cpuTime(model)-time1), 
301      Cbc_getNodeCount(model), Cbc_getObjValue(model));
302    if (VERBOSE>0) (!Cbc_status(model)) ? printf(" Finished\n") : printf(" Not finished\n");
303    if (VERBOSE>2) printSolution(model);
304  }
305  {
306    //
307    // SOS specification data
308    // specify numberSets, numPoints, and whichRanges explicitly
309    // NOTE: These need to be commented according to the MPS file.
310    //   Example_1_1: Cbc0004I MPS reads 1081 only columns out of 1085
311    //                Then the 4th range goes out of bounds
312    //   Example_2: Cbc0006I The LP relaxation is infeasible or too expensive
313    //   Mod_RT_1: Cbc0006I The LP relaxation is infeasible or too expensive
314    //
315    const int numberSets = 
316      1; // dummy
317    //    4; // Example_1_1
318    //    2;  // Example_2
319    //    2;  // Mod_RT_1
320    //    2;  // SITH
321    const int numPoints = 
322      1; // dummy
323    //    257; // Example_1_1
324    //    256;  // Example_2
325    //    257;  // Mod_RT_1
326    //    257;  // SITH
327
328    //const int whichRanges[numberSets][2]
329    const int whichRanges[1][2] = { // counting from zero?
330      {0,0} // dummy
331      //    {56,312}, {313, 569}, {572, 828}, {829, 1085} // Example_1_1
332      //    {48, 303}, {304, 559} // Example_2
333      //    {45, 301}, {302, 558} // Mod_RT_1
334      //    {45, 301}, {302, 558} // SITH
335    };
336    // the rest is determined parametrically
337    int *len = malloc(numberSets* sizeof(int[1]));
338    //int **which = new int * [numberSets];
339    int **which = malloc(numberSets* sizeof(int[1]));
340    int setNum, pointNum;
341    double *weights;
342    int i, j;
343    for (setNum = 0; setNum < numberSets; setNum++) {
344      len[setNum] = whichRanges[setNum][1] - whichRanges[setNum][0] + 1;
345      if (len[setNum] != numPoints) {
346        printf("%s ERROR: len[%i] (%i) != numPoints (%i)\n",
347          prefix, setNum, len[setNum], numPoints);
348        return 1;
349      }
350      which[setNum] = malloc(numPoints*sizeof(int[1]));
351      for (j = 0; j < len[setNum]; j++)
352        //      which[setNum][j] = whichRanges[setNum][0] + j - 6; // hack for missing columns in example 1_1
353        which[setNum][j] = whichRanges[setNum][0] + j; // Example_2
354    }
355    weights = malloc(numPoints*sizeof(double[1])); 
356    for (pointNum = 0; pointNum < numPoints; pointNum++) weights[pointNum] = pointNum+1;
357
358    // Now use SOS2
359    // NOTE: Only enable this if known good SOS Specification (above)
360    if (1) { 
361      if (VERBOSE>1) printf("%s Use SOS2\n",prefix);
362
363      // Restore model
364      Cbc_deleteModel(model);
365      if (0) {
366        model = Cbc_clone(model2);
367      } else {
368        model = getDefaultModel(argc, argv); 
369      }
370      //    Cbc_setLogLevel(model, 4); // 4 is verbose
371      if (VERBOSE>1) printf("%s Model %s has %d rows and %d columns\n",
372        prefix, modelName,Cbc_numberRows(model),Cbc_numberColumns(model));
373
374      // make SOS cuts
375      for (i=0;i<numberIntegers;i++) {
376        int iColumn = integerVariable[i];
377        // Stop being integer
378        Cbc_setContinuous(model, iColumn);
379      }
380      // add cut (0 - use dense, 1 - use sparse (TODO: test this))
381      if (0) for (i=0;i<numberSets;i++) {
382        printf(
383          "%s calling Cbc_addSOS(), len[%i] = %i, which[%i][0] = %i, which[%i][%i] = %i,\n  weights[0] = %g, weights[%i] = %g\n",
384          prefix, i, len[i], i, which[i][0], i, len[i]-1, which[i][len[i]-1], weights[0], len[i]-1, weights[len[i]-1]);
385        //      Cbc_addSOS_Sparse(model,len[i],which[i],weights,i,2);
386      } else {
387        int numObjects = numberSets; // cannot pass const int
388        Cbc_addSOS_Dense(model, numObjects, len, (const int**)which, (const double*)weights, 2);
389      }
390    }
391
392    Cbc_setOptimizationDirection(model, 1); // 1 minimize
393    if (VERBOSE>1) {
394      printf("%s Solve MPS version of the problem\n",prefix);
395      printf("%s Optimization Direction = %g (1 - minimize, -1 - maximize, 0 - ignore)\n",
396        prefix, Cbc_optimizationDirection(model)); 
397    }
398    if (VERBOSE>1) printf("%s calling Cbc_scaling()\n",prefix);
399    Cbc_scaling(model,1);
400    time1 = Cbc_cpuTime(model) ;
401    if (VERBOSE>1) printf("%s calling Cbc_initialSolve()\n",prefix);
402    Cbc_initialSolve(model);
403    if (VERBOSE>3) Cbc_printModel(model,prefix);
404    if (VERBOSE>1) printf("%s calling Cbc_branchAndBound()\n",prefix);
405    Cbc_branchAndBound(model);
406    if (VERBOSE>1) printf("%s %s took %g seconds, %i nodes with objective %g\n",
407      prefix, modelName, (Cbc_cpuTime(model)-time1), 
408      Cbc_getNodeCount(model), Cbc_getObjValue(model));
409    if (VERBOSE>0) (!Cbc_status(model)) ? printf(" Finished\n") : printf(" Not finished\n");
410    if (VERBOSE>2) printSolution(model);
411  }
412
413  if (VERBOSE>1) printf("%s Log Level %i\n", prefix, Cbc_logLevel(model));
414  if (VERBOSE>0) printf("%s return 0\n",prefix);
415  return 0;
416} //  main()
Note: See TracBrowser for help on using the repository browser.