Scippy

SCIP

Solving Constraint Integer Programs

circlepacking.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file circlepacking.c
17  * @brief Packing circles in a rectangle of minimal size.
18  * @author Jose Salmeron
19  * @author Stefan Vigerske
20  *
21  * This example shows how to setup quadratic constraints in SCIP when using SCIP as callable library.
22  * The example implements a model for the computation of a smallest rectangle that contains a number of
23  * given circles in the plane or the computation of the maximal number of circles that can be placed
24  * into a given rectangle.
25  *
26  * Suppose that n circles with radii \f$r_i\f$ are given.
27  * The task is to find coordinates \f$(x_i, y_i)\f$ for the circle midpoints and a rectangle of
28  * width \f$W \geq 0\f$ and height \f$H \geq 0\f$, such that
29  * - every circle is placed within the rectangle (\f$r_i \leq x_i \leq W-r_i\f$, \f$r_i \leq y_i \leq H-r_i\f$),
30  * - circles are not overlapping \f$\left((x_i-x_j)^2 + (y_i-y_j)^2 \geq (r_i + r_j)^2\right)\f$, and
31  * - the area of the rectangle is minimal.
32  *
33  * Alternatively, one may fix the size of the rectangle and maximize the number of circles that
34  * can be fit into the rectangle without being overlapping.
35  */
36 
37 #include <stdio.h>
38 #include <string.h>
39 #include <math.h>
40 
41 #include "scip/scip.h"
42 #include "scip/scipdefplugins.h"
43 
44 /* Model parameters */
45 
46 /** Number of possible circles **/
47 int ncircles = 0;
48 
49 /** Radii **/
51 int rsize = 0;
52 
53 /* Variables */
54 SCIP_VAR** x; /**< x coordinates */
55 SCIP_VAR** y; /**< y coordinates */
56 SCIP_VAR** b; /**< whether circle is placed into rectangle */
57 SCIP_VAR* a; /**< area of rectangle */
58 SCIP_VAR* w; /**< width of rectangle */
59 SCIP_VAR* h; /**< height of rectangle */
60 SCIP_Bool minarea; /**< whether we minimize the area (TRUE) or maximize the number of circles in the rectangle (FALSE) */
61 
62 
63 /** plots solution by use of Python/Matplotlib */
64 static
66  SCIP* scip, /**< SCIP data structure */
67  SCIP_SOL* sol /**< solution to plot */
68  )
69 {
70 #if _POSIX_C_SOURCE < 2
71  SCIPinfoMessage(scip, NULL, "No POSIX version 2. Try http://distrowatch.com/.");
72 #else
73  FILE* stream;
74  int i;
75 
76  stream = popen("python", "w");
77  if( stream == NULL )
78  {
79  SCIPerrorMessage("Could not open pipe to python.\n");
80  return;
81  }
82 
83  fputs("import numpy as np\n"
84  "import matplotlib\n"
85  "import matplotlib.pyplot as plt\n",
86  stream);
87 
88  fputs("fig, ax = plt.subplots()\n"
89  "patches = []\n",
90  stream);
91 
92  for( i = 0; i < ncircles; ++i )
93  {
94  /* The binary variable b[i] indicates which circles should be included in the current solution.
95  * Here we want to skip circles that are not included, that is b[i] is zero (or close to zero due to tolerances).
96  */
97  if( !minarea && SCIPgetSolVal(scip, sol, b[i]) < 0.5 )
98  continue;
99 
100  fprintf(stream, "patches.append(matplotlib.patches.Circle((%g, %g), %g))\n",
101  SCIPgetSolVal(scip, sol, x[i]),
102  SCIPgetSolVal(scip, sol, y[i]),
103  r[i]);
104  }
105 
106  fputs("colors = 100*np.random.rand(len(patches))\n"
107  "p = matplotlib.collections.PatchCollection(patches, alpha=0.4)\n"
108  "p.set_array(np.array(colors))\n"
109  "ax.add_collection(p)\n",
110  stream);
111 
112  fprintf(stream, "plt.xlim(xmax=%g)\n", SCIPgetSolVal(scip, sol, w));
113  fprintf(stream, "plt.ylim(ymax=%g)\n", SCIPgetSolVal(scip, sol, h));
114  if( minarea )
115  fprintf(stream, "plt.title('Area = %.4f')\n", SCIPgetSolVal(scip, sol, a));
116  else
117  fprintf(stream, "plt.title('Number of circles = %d')\n", (int)SCIPround(scip, SCIPgetSolOrigObj(scip, sol)));
118  fputs("plt.gca().set_aspect(1)\n", stream);
119  fputs("plt.show()\n", stream);
120 
121  pclose(stream);
122 #endif
123 }
124 
125 /** plots solution by use of gnuplot */
126 static
128  SCIP* scip, /**< SCIP data structure */
129  SCIP_SOL* sol /**< solution to plot */
130  )
131 {
132 #if _POSIX_C_SOURCE < 2
133  SCIPinfoMessage(scip, NULL, "No POSIX version 2. Try http://distrowatch.com/.");
134 #else
135  SCIP_Real wval;
136  SCIP_Real hval;
137  FILE* stream;
138  int i;
139 
140  /* -p (persist) to keep the plot open after gnuplot terminates */
141  stream = popen("gnuplot -p", "w");
142  if( stream == NULL )
143  {
144  SCIPerrorMessage("Could not open pipe to gnuplot.\n");
145  return;
146  }
147 
148  fputs("unset xtics\n"
149  "unset ytics\n"
150  "unset border\n"
151  "set size ratio 1\n",
152  stream);
153 
154  wval = SCIPgetSolVal(scip, sol, w);
155  hval = SCIPgetSolVal(scip, sol, h);
156 
157  fprintf(stream, "set xrange [0:%.2f]\n", MAX(wval, hval));
158  fprintf(stream, "set yrange [0:%.2f]\n", MAX(wval, hval));
159  fprintf(stream, "set object rectangle from 0,0 to %.2f,%.2f\n", wval, hval);
160  if( minarea )
161  fprintf(stream, "set xlabel \"Area = %.4f\"\n", SCIPgetSolVal(scip, sol, a));
162  else
163  fprintf(stream, "set xlabel \"Number of circles = %d\"\n", (int)SCIPround(scip, SCIPgetSolOrigObj(scip, sol)));
164 
165  fputs("plot '-' with circles notitle\n", stream);
166  for( i = 0; i < ncircles; ++i )
167  {
168  /* skip circles that are not included in current solution, that is b[i] is close to zero */
169  if( !minarea && SCIPgetSolVal(scip, sol, b[i]) < 0.5 )
170  continue;
171 
172  fprintf(stream, "%g %g %g\n",
173  SCIPgetSolVal(scip, sol, x[i]),
174  SCIPgetSolVal(scip, sol, y[i]),
175  r[i]);
176  }
177  fputs("e\n", stream);
178 
179  pclose(stream);
180 #endif
181 }
182 
183 /** plots solution by use of ascii graphics */
184 static
186  SCIP* scip, /**< SCIP data structure */
187  SCIP_SOL* sol /**< solution to plot */
188  )
189 {
190  SCIP_Real wval;
191  SCIP_Real hval;
192  SCIP_Real xval;
193  SCIP_Real yval;
194  SCIP_Real radius;
195  SCIP_Real scalex;
196  SCIP_Real scaley;
197  int dispwidth;
198  int width;
199  int height;
200  char* picture;
201  int i;
202 
203  wval = SCIPgetSolVal(scip, sol, w);
204  hval = SCIPgetSolVal(scip, sol, h);
205 
206  /* scale so that picture is about as wide as SCIP B&B log */
207  SCIP_CALL( SCIPgetIntParam(scip, "display/width", &dispwidth) );
208  scalex = (dispwidth-3) / wval;
209  scaley = scalex / 2.0; /* this gives almost round looking circles on my (SV) terminal */
210 
211  width = SCIPceil(scip, scalex*wval)+3; /* +2 for left and right border and +1 for \n */
212  height = SCIPceil(scip, scaley*hval)+2; /* +2 for top and bottom border */
213 
214  SCIP_CALL( SCIPallocBufferArray(scip, &picture, width * height + 1) );
215 
216  /* initialize with white space and termination */
217  memset(picture, ' ', width * height);
218  picture[width*height] = '\0';
219 
220  /* add border and linebreaks */
221  memset(picture, '*', width-1); /* top border */
222  memset(picture + (height-1) * width, '*', width-1); /* bottom border */
223  for( i = 0; i < height; ++i )
224  {
225  picture[i*width] = '*'; /* left border */
226  picture[i*width + width-2] = '*'; /* right border */
227  picture[i*width + width-1] = '\n';
228  }
229 
230  /* add circles */
231  for( i = 0; i < ncircles; ++i )
232  {
233  SCIP_Real phi;
234  int xcoord;
235  int ycoord;
236 
237  /* skip circles that are not included in current solution, that is b[i] close to zero */
238  if( !minarea && SCIPgetSolVal(scip, sol, b[i]) < 0.5 )
239  continue;
240 
241  xval = SCIPgetSolVal(scip, sol, x[i]);
242  yval = SCIPgetSolVal(scip, sol, y[i]);
243  radius = r[i];
244 
245  for( phi = 0.0; phi < 6.283 /* 2*pi */; phi += 0.01 )
246  {
247  xcoord = SCIPround(scip, scalex * (xval + radius * cos(phi))) + 1; /* +1 for border */
248  ycoord = SCIPround(scip, scaley * (yval + radius * sin(phi))) + 1; /* +1 for border */
249 
250  /* feasible solutions should be within box
251  * due to rounding, they can be on the border
252  */
253  assert(xcoord >= 0);
254  assert(ycoord >= 0);
255  assert(xcoord < width);
256  assert(ycoord < height);
257 
258  picture[ycoord * width + xcoord] = 'a' + i;
259  }
260  }
261 
262  /* print objective value inside top border */
263  i = SCIPsnprintf(picture + width/2 - 8, width/2 + 8,
264  minarea ? " Area = %g " : " #Circles = %.0f ", SCIPgetSolOrigObj(scip, sol));
265  picture[width/2-8+i] = '*';
266 
267  /* show plot */
268  SCIPinfoMessage(scip, NULL, "%s", picture);
269 
270  SCIPfreeBufferArray(scip, &picture);
271 
272  return SCIP_OKAY;
273 }
274 
275 /** initialization method of event handler (called after problem was transformed) */
276 static
277 SCIP_DECL_EVENTINIT(eventInitDispsol)
278 {
280 
281  return SCIP_OKAY;
282 }
283 
284 /** deinitialization method of event handler (called before transformed problem is freed) */
285 static
286 SCIP_DECL_EVENTEXIT(eventExitDispsol)
287 {
289 
290  return SCIP_OKAY;
291 }
292 
293 /** execution method of event handler */
294 static
295 SCIP_DECL_EVENTEXEC(eventExecDispsol)
296 { /*lint --e{715}*/
297  SCIP_SOL* sol;
298 
300 
301  sol = SCIPeventGetSol(event);
302  assert(sol != NULL);
303 
305 
306  return SCIP_OKAY;
307 }
308 
309 /** creates event handler for dispsol event */
310 static
312  SCIP* scip /**< SCIP data structure */
313  )
314 {
315  SCIP_EVENTHDLR* eventhdlr = NULL;
316 
317  /* include event handler into SCIP */
318  SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &eventhdlr, "dispsol", "displays new solutions",
319  eventExecDispsol, NULL) );
320  assert(eventhdlr != NULL);
321 
322  /* set non fundamental callbacks via setter functions */
323  SCIP_CALL( SCIPsetEventhdlrInit(scip, eventhdlr, eventInitDispsol) );
324  SCIP_CALL( SCIPsetEventhdlrExit(scip, eventhdlr, eventExitDispsol) );
325 
326  return SCIP_OKAY;
327 }
328 
329 /** create problem in given SCIP and add all variables and constraints to model the circle packing problem */
331  SCIP* scip, /**< SCIP data structure */
332  SCIP_Real fixwidth, /**< a given fixed width for the rectangle, or SCIP_INVALID if not fixed */
333  SCIP_Real fixheight /**< a given fixed height for the rectangle, or SCIP_INVALID if not fixed */
334  )
335 {
336  char name[SCIP_MAXSTRLEN];
337  SCIP_CONS* cons;
338  SCIP_Real minusone = -1.0;
339  SCIP_Real one = 1.0;
340  int i, j;
341 
342  /* if both width and height are fixed, then we don't optimize the area anymore, but the number of circles */
343  minarea = (fixwidth == SCIP_INVALID || fixheight == SCIP_INVALID);
344 
345  /* create empty problem */
346  SCIP_CALL( SCIPcreateProbBasic(scip, "Packing circles") );
347 
348  /* change to maximization if optimizing number of circles instead of rectangle area */
349  if( !minarea )
350  {
352  }
353 
354  /* Create variables, setup variable bounds, and setup objective function:
355  * We add variables x[i], y[i] for the circle midpoints and, if optimizing the number of circles in the rectangle,
356  * a binary variable b[i] to indicate whether circle i should be in the rectangle or not.
357  * Additionally, we add variables for the rectangle area, width, and height.
358  *
359  * Since the rectangle lower-left corner is assumed to be at (0,0), we can set a lower bound
360  * r[i] for both variables x[i] and y[i].
361  *
362  * In the objective function, we have 1.0*a, if minimizing the area of the rectangle,
363  * otherwise the area does not contribute to the objective, but every b[i] will be present instead.
364  *
365  * Further below we fix the width and height of the rectangle, if fixwidth and fixheight are valid.
366  * If both are valid, then we can also fix the area of the rectangle.
367  */
371  for( i = 0; i < ncircles; ++i )
372  {
373  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x_%d", i);
374  SCIP_CALL( SCIPcreateVarBasic(scip, &x[i], name, r[i], SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
375  SCIP_CALL( SCIPaddVar(scip, x[i]) );
376 
377  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "y_%d", i);
378  SCIP_CALL( SCIPcreateVarBasic(scip, &y[i], name, r[i], SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
379  SCIP_CALL( SCIPaddVar(scip, y[i]) );
380 
381  if( !minarea )
382  {
383  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "b_%d", i);
384  SCIP_CALL( SCIPcreateVarBasic(scip, &b[i], name, 0.0, 1.0, 1.0, SCIP_VARTYPE_BINARY) );
385  SCIP_CALL( SCIPaddVar(scip, b[i]) );
386  }
387  }
388  SCIP_CALL( SCIPcreateVarBasic(scip, &a, "a", 0.0, SCIPinfinity(scip), minarea ? 1.0 : 0.0, SCIP_VARTYPE_CONTINUOUS) );
389  SCIP_CALL( SCIPaddVar(scip, a) );
390 
391  SCIP_CALL( SCIPcreateVarBasic(scip, &w, "w", 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
392  SCIP_CALL( SCIPaddVar(scip, w) );
393 
394  SCIP_CALL( SCIPcreateVarBasic(scip, &h, "h", 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
395  SCIP_CALL( SCIPaddVar(scip, h) );
396 
397  /* fix width if a valid value for fixwidth is given */
398  if( fixwidth != SCIP_INVALID )
399  {
400  SCIP_Bool infeas;
401  SCIP_Bool fixed;
402 
403  SCIP_CALL( SCIPfixVar(scip, w, fixwidth, &infeas, &fixed) );
404 
405  assert(!infeas);
406  assert(fixed);
407  }
408 
409  /* fix height if a valid value for fixheight is given */
410  if( fixheight != SCIP_INVALID )
411  {
412  SCIP_Bool infeas;
413  SCIP_Bool fixed;
414 
415  SCIP_CALL( SCIPfixVar(scip, h, fixheight, &infeas, &fixed) );
416 
417  assert(!infeas);
418  assert(fixed);
419  }
420 
421  /* fix area if both width and height are fixed */
422  if( !minarea )
423  {
424  SCIP_Bool infeas;
425  SCIP_Bool fixed;
426 
427  SCIP_CALL( SCIPfixVar(scip, a, fixwidth * fixheight, &infeas, &fixed) );
428 
429  assert(!infeas);
430  assert(fixed);
431  }
432 
433  /* boundary constraints on circle coordinates
434  *
435  * If minimizing the area of the rectangle, then the coordinates of every circle must
436  * satisfy the boundary conditions r_i <= x_i <= w - r_i and r_i <= y_i <= h - r_i.
437  * The lower bounds r_i are already required by the variable bounds (see above).
438  * For the upper bounds, we add the corresponding linear constraints.
439  *
440  * If the area of the rectangle is fixed, then it would be sufficient to place only
441  * circles into the rectangle that have been decided to be put into the rectangle.
442  * We could model this as a big-M constraint on x_i and y_i, but on the other hand,
443  * we can also require that the circle coordinates always satisfy the boundary conditions,
444  * even if the circle is not placed into the rectangle (b_i=0).
445  * As the non-overlapping constraints do not apply for circles that are not placed into
446  * the rectangle, satisfying these boundary conditions is always possible, unless the
447  * circle itself is too big to be placed into the rectangle. In this case, though,
448  * we can decide a-priori that the circle is not placed into the rectangle, i.e., fix b_i to 0.
449  */
450  for( i = 0; i < ncircles; ++i )
451  {
452  if( !minarea && SCIPisLT(scip, MIN(fixwidth, fixheight), 2*r[i]) )
453  {
454  SCIP_Bool infeas;
455  SCIP_Bool fixed;
456 
457  SCIP_CALL( SCIPfixVar(scip, b[i], 0.0, &infeas, &fixed) );
458 
459  assert(!infeas);
460  assert(fixed);
461 
462  continue;
463  }
464 
465  /* linear constraint: x_i + r_i <= w --> r_i <= w - x_i */
466  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "boundaryright_%d", i, i);
467  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 0, NULL, NULL, r[i], SCIPinfinity(scip)) );
468  SCIP_CALL( SCIPaddCoefLinear(scip, cons, w, 1.0) );
469  SCIP_CALL( SCIPaddCoefLinear(scip, cons, x[i], -1.0) );
470  SCIP_CALL( SCIPaddCons(scip, cons) );
471  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
472 
473  /* linear constraint: y_i + r_i <= h --> r_i <= h - y_i */
474  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "boundarytop_%d", i, i);
475  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 0, NULL, NULL, r[i], SCIPinfinity(scip)) );
476  SCIP_CALL( SCIPaddCoefLinear(scip, cons, h, 1.0) );
477  SCIP_CALL( SCIPaddCoefLinear(scip, cons, y[i], -1.0) );
478  SCIP_CALL( SCIPaddCons(scip, cons) );
479  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
480  }
481 
482  /* constraint that defines the area of the rectangle
483  *
484  * We could add the quadratic constraint w * h - a = 0.
485  * But if we are minimizing a, then we can relax this constraint to w * h - a <= 0.
486  * If the size the rectangle is fixed, then w, h, and a have been fixed above.
487  * We could actually omit this constraint, but here SCIP presolve will take care of removing it.
488  */
489  SCIP_CALL( SCIPcreateConsQuadraticNonlinear(scip, &cons, "area", 1, &a, &minusone, 1, &w, &h, &one, -SCIPinfinity(scip),
490  0.0, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
491  SCIP_CALL( SCIPaddCons(scip, cons) );
492  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
493 
494  /* quadratic constraints that require that circles are not overlapping.
495  * For circles i and j, i<j, we require that the euclidean distance of the circle middle points
496  * is at least the sum of the circle radii, i.e., || (x_i,y_i) - (x_j,y_j) || >= r_i + r_j.
497  * Equivalently, (x_i - x_j)^2 + (y_i - y_j)^2 >= (r_i + r_j)^2, which can be expanded to
498  * x_i^2 - 2 x_i x_j + x_j^2 + y_i^2 - 2 y_i y_j + y_j^2 >= (r_i + r_j)^2
499  *
500  * When not minimizing the area of the circles, however, then this constraint only needs
501  * to hold if both circles are placed into the rectangle, that is if b_i=1 and b_j=1.
502  * We can achieve this by relaxing the right-hand-side to 0 or a negative value if b_i + b_j <= 1:
503  * (x_i - x_j)^2 + (y_i - y_j)^2 >= (r_i + r_j)^2 * (b_i+b_j-1), which can be expanded to
504  * x_i^2 - 2 x_i x_j + x_j^2 + y_i^2 - 2 y_i y_j + y_j^2 - (r_i+r_j)^2 b_i - (r_i+r_j)^2 b_j >= -(r_i + r_j)^2
505  */
506  for( i = 0; i < ncircles; ++i )
507  {
508  for( j = i + 1; j < ncircles; ++j )
509  {
510  SCIP_VAR* quadvars1[6] = {x[i], x[j], x[i], y[i], y[j], y[i]};
511  SCIP_VAR* quadvars2[6] = {x[i], x[j], x[j], y[i], y[j], y[j]};
512  SCIP_Real quadcoefs[6] = {1.0, 1.0, -2.0, 1.0, 1.0, -2.0};
513  SCIP_VAR* linvars[2] = {NULL, NULL};
514  SCIP_Real lincoefs[2] = {0.0, 0.0};
515  int nlinvars = 0;
516 
517  /* create empty quadratic constraint with right-hand-side +/- (r_i - r_j)^2 */
518  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nooverlap_%d,%d", i, j);
519 
520  if( !minarea )
521  {
522  linvars[0] = b[i];
523  lincoefs[0] = -SQR(r[i] + r[j]);
524  linvars[1] = b[j];
525  lincoefs[1] = -SQR(r[i] + r[j]);
526  nlinvars = 2;
527  }
528 
529  SCIP_CALL( SCIPcreateConsQuadraticNonlinear(scip, &cons, name, nlinvars, linvars, lincoefs, 6, quadvars1,
530  quadvars2, quadcoefs, (minarea ? 1.0 : -1.0) * SQR(r[i] + r[j]), SCIPinfinity(scip),
531  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
532 
533  SCIP_CALL( SCIPaddCons(scip, cons) );
534  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
535  }
536  }
537 
538  return SCIP_OKAY;
539 }
540 
541 /** run circle packing example
542  *
543  * Sets up SCIP and the SCIP problem, solves the problem, and shows the solution.
544  */
546  SCIP_Real fixwidth, /**< a given fixed width for the rectangle, or SCIP_INVALID if not fixed */
547  SCIP_Real fixheight, /**< a given fixed height for the rectangle, or SCIP_INVALID if not fixed */
548  SCIP_Bool dognuplot, /**< whether to draw best solution with gnuplot */
549  SCIP_Bool domatplotlib /**< whether to draw best solution with python/matplotlib */
550  )
551 {
552  SCIP* scip;
553  int i;
554 
555  SCIP_CALL( SCIPcreate(&scip) );
558 
559  SCIPinfoMessage(scip, NULL, "\n");
560  SCIPinfoMessage(scip, NULL, "***************************\n");
561  SCIPinfoMessage(scip, NULL, "* Running Packing Circles *\n");
562  SCIPinfoMessage(scip, NULL, "***************************\n");
563  SCIPinfoMessage(scip, NULL, "\n");
564  SCIPinfoMessage(scip, NULL, "%d circles given with radii", ncircles);
565  for( i = 0; i < ncircles; ++i )
566  SCIPinfoMessage(scip, NULL, " %.2f", r[i]);
567  SCIPinfoMessage(scip, NULL, "\n\n");
568 
569  SCIP_CALL( setupProblem(scip, fixwidth, fixheight) );
570 
571  SCIPinfoMessage(scip, NULL, "Original problem:\n");
572  SCIP_CALL( SCIPprintOrigProblem(scip, NULL, "cip", FALSE) );
573 
574  /* By default, SCIP tries to close the gap between primal and dual bound completely.
575  * This can take very long for this example, so we increase the gap tolerance to have
576  * SCIP stop when the distance between primal and dual bound is already below 1e-4.
577  */
578  SCIP_CALL( SCIPsetRealParam(scip, "limits/gap", 1e-4) );
579 
580  SCIPinfoMessage(scip, NULL, "\nSolving...\n");
581  SCIP_CALL( SCIPsolve(scip) );
582 
583  if( SCIPgetNSols(scip) > 0 )
584  {
585  SCIPinfoMessage(scip, NULL, "\nSolution:\n");
586  SCIP_CALL( SCIPprintSol(scip, SCIPgetBestSol(scip), NULL, FALSE) );
587 
588  if( dognuplot )
590 
591  if( domatplotlib )
593  }
594 
595  /* release variables */
596  SCIP_CALL( SCIPreleaseVar(scip, &a) );
597  SCIP_CALL( SCIPreleaseVar(scip, &w) );
598  SCIP_CALL( SCIPreleaseVar(scip, &h) );
599  for( i = 0; i < ncircles; ++i )
600  {
601  SCIP_CALL( SCIPreleaseVar(scip, &x[i]) );
602  SCIP_CALL( SCIPreleaseVar(scip, &y[i]) );
603  if( !minarea )
604  {
605  SCIP_CALL( SCIPreleaseVar(scip, &b[i]) );
606  }
607  }
608 
609  /* free memory arrays */
610  SCIPfreeMemoryArray(scip, &b);
611  SCIPfreeMemoryArray(scip, &y);
612  SCIPfreeMemoryArray(scip, &x);
613 
614  SCIP_CALL( SCIPfree(&scip) );
615 
616  return SCIP_OKAY;
617 }
618 
619 /** main method starting SCIP */
620 int main(
621  int argc, /**< number of arguments from the shell */
622  char** argv /**< array of shell arguments */
623  )
624 { /*lint --e{715}*/
625  SCIP_RETCODE retcode;
626  SCIP_Bool dognuplot = FALSE;
627  SCIP_Bool domatplotlib = FALSE;
628  SCIP_Real fixwidth = SCIP_INVALID;
629  SCIP_Real fixheight = SCIP_INVALID;
630  char* endptr;
631  int i;
632 
633  ncircles = 0;
634  rsize = 5;
636 
637  for( i = 1; i < argc; ++i )
638  {
639  if( strcmp(argv[i], "--help") == 0 )
640  {
641  printf("usage: %s [--help] [-w <width>] [-h <height>]", argv[0]);
642 #if _POSIX_C_SOURCE >= 2
643  printf(" [-g] [-m]");
644 #endif
645  puts(" { <radius> }");
646  puts(" --help shows this help and exits");
647  puts(" -w <width> fix rectangle width to given value");
648  puts(" -h <height> fix rectangle height to given value");
649 #if _POSIX_C_SOURCE >= 2
650  puts(" -g show final solution with gnuplot");
651  puts(" -m show final solution with matplotlib");
652 #endif
653  puts("If no radii are given, then 5 circles with radii 0.25, 0.25, 0.4, 0.7, and 0.1 used.");
654  puts("If both width and height are fixed, then the number of circles that fit into the rectangle is maximized.");
655 
656  return EXIT_SUCCESS;
657  }
658 
659 #if _POSIX_C_SOURCE >= 2
660  if( strcmp(argv[i], "-g") == 0 )
661  {
662  dognuplot = TRUE;
663  continue;
664  }
665 
666  if( strcmp(argv[i], "-m") == 0 )
667  {
668  domatplotlib = TRUE;
669  continue;
670  }
671 #endif
672 
673  if( strcmp(argv[i], "-w") == 0 )
674  {
675  if( i == argc-1 )
676  {
677  fprintf(stderr, "ERROR: Missing argument for -w.\n");
678  return EXIT_FAILURE;
679  }
680 
681  fixwidth = strtod(argv[i+1], &endptr);
682  if( *endptr != '\0' )
683  {
684  fprintf(stderr, "ERROR: Could not parse argument %s into a number.\n", argv[i+1]);
685  return EXIT_FAILURE;
686  }
687 
688  ++i;
689  continue;
690  }
691 
692  if( strcmp(argv[i], "-h") == 0 )
693  {
694  if( i == argc-1 )
695  {
696  fprintf(stderr, "ERROR: Missing argument for -h.\n");
697  return EXIT_FAILURE;
698  }
699 
700  fixheight = strtod(argv[i+1], &endptr);
701  if( *endptr != '\0' )
702  {
703  fprintf(stderr, "ERROR: Could not parse argument %s into a number.\n", argv[i+1]);
704  return EXIT_FAILURE;
705  }
706 
707  ++i;
708  continue;
709  }
710 
711  /* see whether the argument can be parsed as a positive real number */
712  if( rsize <= ncircles )
713  {
714  rsize += 5;
716  }
717  r[ncircles] = strtod(argv[i], &endptr);
718  if( *endptr == '\0' && endptr != argv[i] && r[ncircles] > 0.0 )
719  {
720  ++ncircles;
721  continue;
722  }
723 
724  fprintf(stderr, "ERROR: Unknown option %s.\n", argv[i]);
725  return EXIT_FAILURE;
726  }
727 
728  if( ncircles == 0 )
729  {
730  assert(rsize >= 5);
731  r[0] = 0.25;
732  r[1] = 0.25;
733  r[2] = 0.4;
734  r[3] = 0.7;
735  r[4] = 0.1;
736  ncircles = 5;
737  }
738 
739  /* run the actual circle packing example (setting up SCIP, solving the problem, showing the solution) */
740  retcode = runPacking(fixwidth, fixheight, dognuplot, domatplotlib);
741 
742  /* evaluate return code of the SCIP process */
743  if( retcode != SCIP_OKAY )
744  {
745  /* write error back trace */
746  SCIPprintError(retcode);
747  return EXIT_FAILURE;
748  }
749 
750  return EXIT_SUCCESS;
751 }
static SCIP_RETCODE includeEventHdlrDispsol(SCIP *scip)
static void visualizeSolutionMatplotlib(SCIP *scip, SCIP_SOL *sol)
Definition: circlepacking.c:65
SCIP_Bool minarea
Definition: circlepacking.c:60
SCIP_RETCODE SCIPcreateConsBasicLinear(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, SCIP_VAR **vars, SCIP_Real *vals, SCIP_Real lhs, SCIP_Real rhs)
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
static SCIP_DECL_EVENTINIT(eventInitDispsol)
#define SCIP_MAXSTRLEN
Definition: def.h:293
SCIP_RETCODE SCIPsetEventhdlrExit(SCIP *scip, SCIP_EVENTHDLR *eventhdlr, SCIP_DECL_EVENTEXIT((*eventexit)))
Definition: scip_event.c:169
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
SCIP_RETCODE SCIPincludeEventhdlrBasic(SCIP *scip, SCIP_EVENTHDLR **eventhdlrptr, const char *name, const char *desc, SCIP_DECL_EVENTEXEC((*eventexec)), SCIP_EVENTHDLRDATA *eventhdlrdata)
Definition: scip_event.c:95
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1245
#define FALSE
Definition: def.h:87
SCIP_Real SCIPinfinity(SCIP *scip)
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10755
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
static SCIP_DECL_EVENTEXEC(eventExecDispsol)
SCIP_RETCODE SCIPcreateVarBasic(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype)
Definition: scip_var.c:185
#define BMSallocMemoryArray(ptr, num)
Definition: memory.h:116
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
SCIP_RETCODE SCIPcreate(SCIP **scip)
Definition: scip_general.c:283
SCIP_RETCODE SCIPsetRealParam(SCIP *scip, const char *name, SCIP_Real value)
Definition: scip_param.c:594
SCIP_VAR ** x
Definition: circlepacking.c:54
SCIP_RETCODE SCIPaddCoefLinear(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real val)
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:199
SCIP_RETCODE SCIPcreateProbBasic(SCIP *scip, const char *name)
Definition: scip_prob.c:170
#define SCIP_ALLOC_ABORT(x)
Definition: def.h:374
SCIP_VAR * w
Definition: circlepacking.c:58
SCIP_RETCODE SCIPsetObjsense(SCIP *scip, SCIP_OBJSENSE objsense)
Definition: scip_prob.c:1240
static SCIP_DECL_EVENTEXIT(eventExitDispsol)
SCIP_RETCODE SCIPsolve(SCIP *scip)
Definition: scip_solve.c:2613
#define SCIPerrorMessage
Definition: pub_message.h:55
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2768
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPsetEventhdlrInit(SCIP *scip, SCIP_EVENTHDLR *eventhdlr, SCIP_DECL_EVENTINIT((*eventinit)))
Definition: scip_event.c:155
SCIP_RETCODE SCIPcreateConsQuadraticNonlinear(SCIP *scip, SCIP_CONS **cons, const char *name, int nlinvars, SCIP_VAR **linvars, SCIP_Real *lincoefs, int nquadterms, SCIP_VAR **quadvars1, SCIP_VAR **quadvars2, SCIP_Real *quadcoefs, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool initial, SCIP_Bool separate, SCIP_Bool enforce, SCIP_Bool check, SCIP_Bool propagate, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool dynamic, SCIP_Bool removable)
static SCIP_Real phi(SCIP *scip, SCIP_Real val, SCIP_Real lb, SCIP_Real ub)
Definition: sepa_eccuts.c:837
int rsize
Definition: circlepacking.c:51
SCIP_RETCODE SCIPprintOrigProblem(SCIP *scip, FILE *file, const char *extension, SCIP_Bool genericnames)
#define NULL
Definition: lpi_spx1.cpp:155
SCIP_RETCODE SCIPgetIntParam(SCIP *scip, const char *name, int *value)
Definition: scip_param.c:260
#define SCIP_CALL(x)
Definition: def.h:384
SCIP_VAR * h
Definition: circlepacking.c:59
static SCIP_RETCODE runPacking(SCIP_Real fixwidth, SCIP_Real fixheight, SCIP_Bool dognuplot, SCIP_Bool domatplotlib)
static SCIP_RETCODE visualizeSolutionAscii(SCIP *scip, SCIP_SOL *sol)
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
#define SCIP_Bool
Definition: def.h:84
SCIP_RETCODE SCIPincludeDefaultPlugins(SCIP *scip)
SCIP_RETCODE SCIPcatchEvent(SCIP *scip, SCIP_EVENTTYPE eventtype, SCIP_EVENTHDLR *eventhdlr, SCIP_EVENTDATA *eventdata, int *filterpos)
Definition: scip_event.c:277
SCIP_EVENTTYPE SCIPeventGetType(SCIP_EVENT *event)
Definition: event.c:1021
static void visualizeSolutionGnuplot(SCIP *scip, SCIP_SOL *sol)
static SCIP_RETCODE setupProblem(SCIP *scip, SCIP_Real fixwidth, SCIP_Real fixheight)
#define MAX(x, y)
Definition: tclique_def.h:83
SCIP_RETCODE SCIPdropEvent(SCIP *scip, SCIP_EVENTTYPE eventtype, SCIP_EVENTHDLR *eventhdlr, SCIP_EVENTDATA *eventdata, int filterpos)
Definition: scip_event.c:311
int SCIPgetNSols(SCIP *scip)
Definition: scip_sol.c:2205
SCIP_RETCODE SCIPfixVar(SCIP *scip, SCIP_VAR *var, SCIP_Real fixedval, SCIP_Bool *infeasible, SCIP_Bool *fixed)
Definition: scip_var.c:8273
SCIP_Real SCIPgetSolOrigObj(SCIP *scip, SCIP_SOL *sol)
Definition: scip_sol.c:1435
SCIP_Real * r
Definition: circlepacking.c:50
SCIP_VAR ** b
Definition: circlepacking.c:56
#define SCIP_EVENTTYPE_BESTSOLFOUND
Definition: type_event.h:96
SCIP_SOL * SCIPgetBestSol(SCIP *scip)
Definition: scip_sol.c:2304
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_prob.c:1666
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1110
int ncircles
Definition: circlepacking.c:47
SCIP_VAR * a
Definition: circlepacking.c:57
#define SCIP_Real
Definition: def.h:177
SCIP_VAR ** y
Definition: circlepacking.c:55
int main(int argc, char **argv)
#define SCIP_INVALID
Definition: def.h:197
#define BMSreallocMemoryArray(ptr, num)
Definition: memory.h:120
SCIP_Real SCIPceil(SCIP *scip, SCIP_Real val)
void SCIPprintError(SCIP_RETCODE retcode)
Definition: scip_general.c:211
SCIP_Real SCIPround(SCIP *scip, SCIP_Real val)
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1352
default SCIP plugins
SCIP_SOL * SCIPeventGetSol(SCIP_EVENT *event)
Definition: event.c:1328
SCIP callable library.
SCIP_RETCODE SCIPfree(SCIP **scip)
Definition: scip_general.c:315
SCIP_RETCODE SCIPprintSol(SCIP *scip, SCIP_SOL *sol, FILE *file, SCIP_Bool printzeros)
Definition: scip_sol.c:1766