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

Last change on this file since 1898 was 1898, checked in by stefan, 6 years ago

fixup examples

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