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