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

Last change on this file since 1464 was 1464, checked in by stefan, 9 years ago

merge split branch into trunk; fix some examples

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