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-2019 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, 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 one = 1.0;
339  int i, j;
340 
341  /* if both width and height are fixed, then we don't optimize the area anymore, but the number of circles */
342  minarea = (fixwidth == SCIP_INVALID || fixheight == SCIP_INVALID);
343 
344  /* create empty problem */
345  SCIP_CALL( SCIPcreateProbBasic(scip, "Packing circles") );
346 
347  /* change to maximization if optimizing number of circles instead of rectangle area */
348  if( !minarea )
349  {
351  }
352 
353  /* Create variables, setup variable bounds, and setup objective function:
354  * We add variables x[i], y[i] for the circle midpoints and, if optimizing the number of circles in the rectangle,
355  * a binary variable b[i] to indicate whether circle i should be in the rectangle or not.
356  * Additionally, we add variables for the rectangle area, width, and height.
357  *
358  * Since the rectangle lower-left corner is assumed to be at (0,0), we can set a lower bound
359  * r[i] for both variables x[i] and y[i].
360  *
361  * In the objective function, we have 1.0*a, if minimizing the area of the rectangle,
362  * otherwise the area does not contribute to the objective, but every b[i] will be present instead.
363  *
364  * Further below we fix the width and height of the rectangle, if fixwidth and fixheight are valid.
365  * If both are valid, then we can also fix the area of the rectangle.
366  */
370  for( i = 0; i < ncircles; ++i )
371  {
372  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x_%d", i);
373  SCIP_CALL( SCIPcreateVarBasic(scip, &x[i], name, r[i], SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
374  SCIP_CALL( SCIPaddVar(scip, x[i]) );
375 
376  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "y_%d", i);
377  SCIP_CALL( SCIPcreateVarBasic(scip, &y[i], name, r[i], SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
378  SCIP_CALL( SCIPaddVar(scip, y[i]) );
379 
380  if( !minarea )
381  {
382  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "b_%d", i);
383  SCIP_CALL( SCIPcreateVarBasic(scip, &b[i], name, 0.0, 1.0, 1.0, SCIP_VARTYPE_BINARY) );
384  SCIP_CALL( SCIPaddVar(scip, b[i]) );
385  }
386  }
387  SCIP_CALL( SCIPcreateVarBasic(scip, &a, "a", 0.0, SCIPinfinity(scip), minarea ? 1.0 : 0.0, SCIP_VARTYPE_CONTINUOUS) );
388  SCIP_CALL( SCIPaddVar(scip, a) );
389 
390  SCIP_CALL( SCIPcreateVarBasic(scip, &w, "w", 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
391  SCIP_CALL( SCIPaddVar(scip, w) );
392 
393  SCIP_CALL( SCIPcreateVarBasic(scip, &h, "h", 0.0, SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
394  SCIP_CALL( SCIPaddVar(scip, h) );
395 
396  /* fix width if a valid value for fixwidth is given */
397  if( fixwidth != SCIP_INVALID )
398  {
399  SCIP_Bool infeas;
400  SCIP_Bool fixed;
401 
402  SCIP_CALL( SCIPfixVar(scip, w, fixwidth, &infeas, &fixed) );
403 
404  assert(!infeas);
405  assert(fixed);
406  }
407 
408  /* fix height if a valid value for fixheight is given */
409  if( fixheight != SCIP_INVALID )
410  {
411  SCIP_Bool infeas;
412  SCIP_Bool fixed;
413 
414  SCIP_CALL( SCIPfixVar(scip, h, fixheight, &infeas, &fixed) );
415 
416  assert(!infeas);
417  assert(fixed);
418  }
419 
420  /* fix area if both width and height are fixed */
421  if( !minarea )
422  {
423  SCIP_Bool infeas;
424  SCIP_Bool fixed;
425 
426  SCIP_CALL( SCIPfixVar(scip, a, fixwidth * fixheight, &infeas, &fixed) );
427 
428  assert(!infeas);
429  assert(fixed);
430  }
431 
432  /* boundary constraints on circle coordinates
433  *
434  * If minimizing the area of the rectangle, then the coordinates of every circle must
435  * satisfy the boundary conditions r_i <= x_i <= w - r_i and r_i <= y_i <= h - r_i.
436  * The lower bounds r_i are already required by the variable bounds (see above).
437  * For the upper bounds, we add the corresponding linear constraints.
438  *
439  * If the area of the rectangle is fixed, then it would be sufficient to place only
440  * circles into the rectangle that have been decided to be put into the rectangle.
441  * We could model this as a big-M constraint on x_i and y_i, but on the other hand,
442  * we can also require that the circle coordinates always satisfy the boundary conditions,
443  * even if the circle is not placed into the rectangle (b_i=0).
444  * As the non-overlapping constraints do not apply for circles that are not placed into
445  * the rectangle, satisfying these boundary conditions is always possible, unless the
446  * circle itself is too big to be placed into the rectangle. In this case, though,
447  * we can decide a-priori that the circle is not placed into the rectangle, i.e., fix b_i to 0.
448  */
449  for( i = 0; i < ncircles; ++i )
450  {
451  if( !minarea && SCIPisLT(scip, MIN(fixwidth, fixheight), 2*r[i]) )
452  {
453  SCIP_Bool infeas;
454  SCIP_Bool fixed;
455 
456  SCIP_CALL( SCIPfixVar(scip, b[i], 0.0, &infeas, &fixed) );
457 
458  assert(!infeas);
459  assert(fixed);
460 
461  continue;
462  }
463 
464  /* linear constraint: x_i + r_i <= w --> r_i <= w - x_i */
465  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "boundaryright_%d", i, i);
466  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 0, NULL, NULL, r[i], SCIPinfinity(scip)) );
467  SCIP_CALL( SCIPaddCoefLinear(scip, cons, w, 1.0) );
468  SCIP_CALL( SCIPaddCoefLinear(scip, cons, x[i], -1.0) );
469  SCIP_CALL( SCIPaddCons(scip, cons) );
470  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
471 
472  /* linear constraint: y_i + r_i <= h --> r_i <= h - y_i */
473  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "boundarytop_%d", i, i);
474  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 0, NULL, NULL, r[i], SCIPinfinity(scip)) );
475  SCIP_CALL( SCIPaddCoefLinear(scip, cons, h, 1.0) );
476  SCIP_CALL( SCIPaddCoefLinear(scip, cons, y[i], -1.0) );
477  SCIP_CALL( SCIPaddCons(scip, cons) );
478  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
479  }
480 
481  /* constraint that defines the area of the rectangle
482  *
483  * We could add the quadratic constraint w * h - a = 0.
484  * But if we are minimizing a, then we can relax this constraint to w * h - a <= 0.
485  * If the size the rectangle is fixed, then w, h, and a have been fixed above.
486  * We could actually omit this constraint, but here SCIP presolve will take care of removing it.
487  */
488  SCIP_CALL( SCIPcreateConsBasicQuadratic(scip, &cons, "area", 0, NULL, NULL, 1, &w, &h, &one, -SCIPinfinity(scip), 0.0) );
489  SCIP_CALL( SCIPaddLinearVarQuadratic(scip, cons, a, -1.0) );
490  SCIP_CALL( SCIPaddCons(scip, cons) );
491  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
492 
493  /* quadratic constraints that require that circles are not overlapping.
494  * For circles i and j, i<j, we require that the euclidean distance of the circle middle points
495  * is at least the sum of the circle radii, i.e., || (x_i,y_i) - (x_j,y_j) || >= r_i + r_j.
496  * Equivalently, (x_i - x_j)^2 + (y_i - y_j)^2 >= (r_i + r_j)^2, which can be expanded to
497  * 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
498  *
499  * When not minimizing the area of the circles, however, then this constraint only needs
500  * to hold if both circles are placed into the rectangle, that is if b_i=1 and b_j=1.
501  * We can achieve this by relaxing the right-hand-side to 0 or a negative value if b_i + b_j <= 1:
502  * (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
503  * 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
504  */
505  for( i = 0; i < ncircles; ++i )
506  {
507  for( j = i + 1; j < ncircles; ++j )
508  {
509  /* create empty quadratic constraint with right-hand-side +/- (r_i - r_j)^2 */
510  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "nooverlap_%d,%d", i, j);
511  SCIP_CALL( SCIPcreateConsBasicQuadratic(scip, &cons, name, 0, NULL, NULL, 0, NULL, NULL, NULL, (minarea ? 1.0 : -1.0) * SQR(r[i] + r[j]), SCIPinfinity(scip)) );
512 
513  SCIP_CALL( SCIPaddSquareCoefQuadratic(scip, cons, x[i], 1.0) ); /* x_i^2 */
514  SCIP_CALL( SCIPaddSquareCoefQuadratic(scip, cons, x[j], 1.0) ); /* x_j^2 */
515  SCIP_CALL( SCIPaddBilinTermQuadratic(scip, cons, x[i], x[j], -2.0) ); /* - 2 x_i x_j */
516 
517  SCIP_CALL( SCIPaddSquareCoefQuadratic(scip, cons, y[i], 1.0) ); /* y_i^2 */
518  SCIP_CALL( SCIPaddSquareCoefQuadratic(scip, cons, y[j], 1.0) ); /* y_j^2 */
519  SCIP_CALL( SCIPaddBilinTermQuadratic(scip, cons, y[i], y[j], -2.0) ); /* - 2 y_i y_j */
520 
521  if( !minarea )
522  {
523  /* add -(r_i+r_j)^2 (b_i + b_j) */
524  SCIP_CALL( SCIPaddLinearVarQuadratic(scip, cons, b[i], -SQR(r[i] + r[j])) );
525  SCIP_CALL( SCIPaddLinearVarQuadratic(scip, cons, b[j], -SQR(r[i] + r[j])) );
526  }
527 
528  SCIP_CALL( SCIPaddCons(scip, cons) );
529  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
530  }
531  }
532 
533  return SCIP_OKAY;
534 }
535 
536 /** run circle packing example
537  *
538  * Sets up SCIP and the SCIP problem, solves the problem, and shows the solution.
539  */
541  SCIP_Real fixwidth, /**< a given fixed width for the rectangle, or SCIP_INVALID if not fixed */
542  SCIP_Real fixheight, /**< a given fixed height for the rectangle, or SCIP_INVALID if not fixed */
543  SCIP_Bool dognuplot, /**< whether to draw best solution with gnuplot */
544  SCIP_Bool domatplotlib /**< whether to draw best solution with python/matplotlib */
545  )
546 {
547  SCIP* scip;
548  int i;
549 
550  SCIP_CALL( SCIPcreate(&scip) );
553 
554  SCIPinfoMessage(scip, NULL, "\n");
555  SCIPinfoMessage(scip, NULL, "***************************\n");
556  SCIPinfoMessage(scip, NULL, "* Running Packing Circles *\n");
557  SCIPinfoMessage(scip, NULL, "***************************\n");
558  SCIPinfoMessage(scip, NULL, "\n");
559  SCIPinfoMessage(scip, NULL, "%d circles given with radii", ncircles);
560  for( i = 0; i < ncircles; ++i )
561  SCIPinfoMessage(scip, NULL, " %.2f", r[i]);
562  SCIPinfoMessage(scip, NULL, "\n\n");
563 
564  SCIP_CALL( setupProblem(scip, fixwidth, fixheight) );
565 
566  SCIPinfoMessage(scip, NULL, "Original problem:\n");
567  SCIP_CALL( SCIPprintOrigProblem(scip, NULL, "cip", FALSE) );
568 
569  /* By default, SCIP tries to close the gap between primal and dual bound completely.
570  * This can take very long for this example, so we increase the gap tolerance to have
571  * SCIP stop when the distance between primal and dual bound is already below 1e-4.
572  */
573  SCIP_CALL( SCIPsetRealParam(scip, "limits/gap", 1e-4) );
574 
575  SCIPinfoMessage(scip, NULL, "\nSolving...\n");
576  SCIP_CALL( SCIPsolve(scip) );
577 
578  if( SCIPgetNSols(scip) > 0 )
579  {
580  SCIPinfoMessage(scip, NULL, "\nSolution:\n");
581  SCIP_CALL( SCIPprintSol(scip, SCIPgetBestSol(scip), NULL, FALSE) );
582 
583  if( dognuplot )
585 
586  if( domatplotlib )
588  }
589 
590  /* release variables */
591  SCIP_CALL( SCIPreleaseVar(scip, &a) );
592  SCIP_CALL( SCIPreleaseVar(scip, &w) );
593  SCIP_CALL( SCIPreleaseVar(scip, &h) );
594  for( i = 0; i < ncircles; ++i )
595  {
596  SCIP_CALL( SCIPreleaseVar(scip, &x[i]) );
597  SCIP_CALL( SCIPreleaseVar(scip, &y[i]) );
598  if( !minarea )
599  {
600  SCIP_CALL( SCIPreleaseVar(scip, &b[i]) );
601  }
602  }
603 
604  /* free memory arrays */
605  SCIPfreeMemoryArray(scip, &b);
606  SCIPfreeMemoryArray(scip, &y);
607  SCIPfreeMemoryArray(scip, &x);
608 
609  SCIP_CALL( SCIPfree(&scip) );
610 
611  return SCIP_OKAY;
612 }
613 
614 /** main method starting SCIP */
615 int main(
616  int argc, /**< number of arguments from the shell */
617  char** argv /**< array of shell arguments */
618  )
619 { /*lint --e{715}*/
620  SCIP_RETCODE retcode;
621  SCIP_Bool dognuplot = FALSE;
622  SCIP_Bool domatplotlib = FALSE;
623  SCIP_Real fixwidth = SCIP_INVALID;
624  SCIP_Real fixheight = SCIP_INVALID;
625  char* endptr;
626  int i;
627 
628  ncircles = 0;
629  rsize = 5;
631 
632  for( i = 1; i < argc; ++i )
633  {
634  if( strcmp(argv[i], "--help") == 0 )
635  {
636  printf("usage: %s [--help] [-w <width>] [-h <height>]", argv[0]);
637 #if _POSIX_C_SOURCE >= 2
638  printf(" [-g] [-m]");
639 #endif
640  puts(" { <radius> }");
641  puts(" --help shows this help and exits");
642  puts(" -w <width> fix rectangle width to given value");
643  puts(" -h <height> fix rectangle height to given value");
644 #if _POSIX_C_SOURCE >= 2
645  puts(" -g show final solution with gnuplot");
646  puts(" -m show final solution with matplotlib");
647 #endif
648  puts("If no radii are given, then 5 circles with radii 0.25, 0.25, 0.4, 0.7, and 0.1 used.");
649  puts("If both width and height are fixed, then the number of circles that fit into the rectangle is maximized.");
650 
651  return EXIT_SUCCESS;
652  }
653 
654 #if _POSIX_C_SOURCE >= 2
655  if( strcmp(argv[i], "-g") == 0 )
656  {
657  dognuplot = TRUE;
658  continue;
659  }
660 
661  if( strcmp(argv[i], "-m") == 0 )
662  {
663  domatplotlib = TRUE;
664  continue;
665  }
666 #endif
667 
668  if( strcmp(argv[i], "-w") == 0 )
669  {
670  if( i == argc-1 )
671  {
672  fprintf(stderr, "ERROR: Missing argument for -w.\n");
673  return EXIT_FAILURE;
674  }
675 
676  fixwidth = strtod(argv[i+1], &endptr);
677  if( *endptr != '\0' )
678  {
679  fprintf(stderr, "ERROR: Could not parse argument %s into a number.\n", argv[i+1]);
680  return EXIT_FAILURE;
681  }
682 
683  ++i;
684  continue;
685  }
686 
687  if( strcmp(argv[i], "-h") == 0 )
688  {
689  if( i == argc-1 )
690  {
691  fprintf(stderr, "ERROR: Missing argument for -h.\n");
692  return EXIT_FAILURE;
693  }
694 
695  fixheight = strtod(argv[i+1], &endptr);
696  if( *endptr != '\0' )
697  {
698  fprintf(stderr, "ERROR: Could not parse argument %s into a number.\n", argv[i+1]);
699  return EXIT_FAILURE;
700  }
701 
702  ++i;
703  continue;
704  }
705 
706  /* see whether the argument can be parsed as a positive real number */
707  if( rsize <= ncircles )
708  {
709  rsize += 5;
711  }
712  r[ncircles] = strtod(argv[i], &endptr);
713  if( *endptr == '\0' && endptr != argv[i] && r[ncircles] > 0.0 )
714  {
715  ++ncircles;
716  continue;
717  }
718 
719  fprintf(stderr, "ERROR: Unknown option %s.\n", argv[i]);
720  return EXIT_FAILURE;
721  }
722 
723  if( ncircles == 0 )
724  {
725  assert(rsize >= 5);
726  r[0] = 0.25;
727  r[1] = 0.25;
728  r[2] = 0.4;
729  r[3] = 0.7;
730  r[4] = 0.1;
731  ncircles = 5;
732  }
733 
734  /* run the actual circle packing example (setting up SCIP, solving the problem, showing the solution) */
735  retcode = runPacking(fixwidth, fixheight, dognuplot, domatplotlib);
736 
737  /* evaluate return code of the SCIP process */
738  if( retcode != SCIP_OKAY )
739  {
740  /* write error back trace */
741  SCIPprintError(retcode);
742  return EXIT_FAILURE;
743  }
744 
745  return EXIT_SUCCESS;
746 }
static SCIP_RETCODE includeEventHdlrDispsol(SCIP *scip)
SCIP_Real SCIPround(SCIP *scip, SCIP_Real val)
void SCIPprintError(SCIP_RETCODE retcode)
Definition: scip_general.c:210
static void visualizeSolutionMatplotlib(SCIP *scip, SCIP_SOL *sol)
Definition: circlepacking.c:65
SCIP_Bool minarea
Definition: circlepacking.c:60
#define NULL
Definition: def.h:253
SCIP_RETCODE SCIPcreateConsBasicQuadratic(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)
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:69
static SCIP_DECL_EVENTINIT(eventInitDispsol)
#define SCIP_MAXSTRLEN
Definition: def.h:274
SCIP_RETCODE SCIPsetObjsense(SCIP *scip, SCIP_OBJSENSE objsense)
Definition: scip_prob.c:1241
SCIP_RETCODE SCIPaddSquareCoefQuadratic(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real coef)
#define SQR(x)
Definition: def.h:205
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:53
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1352
SCIP_RETCODE SCIPaddCoefLinear(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real val)
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:184
SCIPInterval cos(const SCIPInterval &x)
#define FALSE
Definition: def.h:73
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 TRUE
Definition: def.h:72
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
static SCIP_DECL_EVENTEXEC(eventExecDispsol)
SCIP_RETCODE SCIPsolve(SCIP *scip)
Definition: scip_solve.c:2535
#define BMSallocMemoryArray(ptr, num)
Definition: memory.h:113
SCIP_Real SCIPceil(SCIP *scip, SCIP_Real val)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:123
SCIP_VAR ** x
Definition: circlepacking.c:54
int SCIPgetNSols(SCIP *scip)
Definition: scip_sol.c:2205
SCIP_RETCODE SCIPsetRealParam(SCIP *scip, const char *name, SCIP_Real value)
Definition: scip_param.c:612
#define SCIP_ALLOC_ABORT(x)
Definition: def.h:355
SCIP_VAR * w
Definition: circlepacking.c:58
SCIP_RETCODE SCIPcreate(SCIP **scip)
Definition: scip_general.c:282
static SCIP_DECL_EVENTEXIT(eventExitDispsol)
#define SCIPerrorMessage
Definition: pub_message.h:45
SCIP_RETCODE SCIPcatchEvent(SCIP *scip, SCIP_EVENTTYPE eventtype, SCIP_EVENTHDLR *eventhdlr, SCIP_EVENTDATA *eventdata, int *filterpos)
Definition: scip_event.c:276
static SCIP_Real phi(SCIP *scip, SCIP_Real val, SCIP_Real lb, SCIP_Real ub)
Definition: sepa_eccuts.c:747
int rsize
Definition: circlepacking.c:51
SCIP_RETCODE SCIPprintOrigProblem(SCIP *scip, FILE *file, const char *extension, SCIP_Bool genericnames)
#define SCIP_CALL(x)
Definition: def.h:365
SCIP_RETCODE SCIPaddBilinTermQuadratic(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var1, SCIP_VAR *var2, SCIP_Real coef)
SCIP_EVENTTYPE SCIPeventGetType(SCIP_EVENT *event)
Definition: event.c:995
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)
SCIP_RETCODE SCIPgetIntParam(SCIP *scip, const char *name, int *value)
Definition: scip_param.c:259
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:111
SCIP_Real SCIPinfinity(SCIP *scip)
#define SCIP_Bool
Definition: def.h:70
SCIP_RETCODE SCIPincludeDefaultPlugins(SCIP *scip)
SCIPInterval sin(const SCIPInterval &x)
static void visualizeSolutionGnuplot(SCIP *scip, SCIP_SOL *sol)
static SCIP_RETCODE setupProblem(SCIP *scip, SCIP_Real fixwidth, SCIP_Real fixheight)
#define MIN(x, y)
Definition: def.h:223
SCIP_RETCODE SCIPfixVar(SCIP *scip, SCIP_VAR *var, SCIP_Real fixedval, SCIP_Bool *infeasible, SCIP_Bool *fixed)
Definition: scip_var.c:8180
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_prob.c:1667
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:94
SCIP_RETCODE SCIPaddLinearVarQuadratic(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real coef)
SCIP_Real * r
Definition: circlepacking.c:50
SCIP_VAR ** b
Definition: circlepacking.c:56
#define SCIP_EVENTTYPE_BESTSOLFOUND
Definition: type_event.h:88
#define MAX(x, y)
Definition: def.h:222
SCIP_RETCODE SCIPsetEventhdlrInit(SCIP *scip, SCIP_EVENTHDLR *eventhdlr, SCIP_DECL_EVENTINIT((*eventinit)))
Definition: scip_event.c:154
SCIP_RETCODE SCIPdropEvent(SCIP *scip, SCIP_EVENTTYPE eventtype, SCIP_EVENTHDLR *eventhdlr, SCIP_EVENTDATA *eventdata, int filterpos)
Definition: scip_event.c:310
int ncircles
Definition: circlepacking.c:47
SCIP_VAR * a
Definition: circlepacking.c:57
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10263
SCIP_RETCODE SCIPprintSol(SCIP *scip, SCIP_SOL *sol, FILE *file, SCIP_Bool printzeros)
Definition: scip_sol.c:1766
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1251
#define SCIP_Real
Definition: def.h:164
SCIP_VAR ** y
Definition: circlepacking.c:55
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
int main(int argc, char **argv)
#define SCIP_INVALID
Definition: def.h:184
#define BMSreallocMemoryArray(ptr, num)
Definition: memory.h:117
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2765
SCIP_RETCODE SCIPfree(SCIP **scip)
Definition: scip_general.c:314
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:198
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1109
default SCIP plugins
SCIP_SOL * SCIPgetBestSol(SCIP *scip)
Definition: scip_sol.c:2304
SCIP_SOL * SCIPeventGetSol(SCIP_EVENT *event)
Definition: event.c:1259
SCIP_RETCODE SCIPcreateProbBasic(SCIP *scip, const char *name)
Definition: scip_prob.c:166
SCIP callable library.
SCIP_RETCODE SCIPsetEventhdlrExit(SCIP *scip, SCIP_EVENTHDLR *eventhdlr, SCIP_DECL_EVENTEXIT((*eventexit)))
Definition: scip_event.c:168
SCIP_Real SCIPgetSolOrigObj(SCIP *scip, SCIP_SOL *sol)
Definition: scip_sol.c:1435