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$
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 */
18 void printSolution(Cbc_Model *cbc_model);
19 Cbc_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 */
callBack(Cbc_Model * model,int messageNumber,int nDouble,const double * vDouble,int nInt,const int * vInt,int nString,char ** vString)24 static 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 */
getDefaultModel(int argc,const char * argv[])53 Cbc_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
printSolution(Cbc_Model * cbc_model)107 void 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
main(int argc,const char * argv[])219 int 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() */
428