source: trunk/Cbc/examples/cbc_driverC_sos.c @ 1468

Last change on this file since 1468 was 1468, checked in by stefan, 11 years ago

do not require CbcConfig?.h in example to decide whether sample or miplib3 is present - do this in makefile

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