probdata_rpa.c

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*  Copyright (c) 2002-2025 Zuse Institute Berlin (ZIB)                      */
/*                                                                           */
/*  Licensed under the Apache License, Version 2.0 (the "License");          */
/*  you may not use this file except in compliance with the License.         */
/*  You may obtain a copy of the License at                                  */
/*                                                                           */
/*      http://www.apache.org/licenses/LICENSE-2.0                           */
/*                                                                           */
/*  Unless required by applicable law or agreed to in writing, software      */
/*  distributed under the License is distributed on an "AS IS" BASIS,        */
/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
/*  See the License for the specific language governing permissions and      */
/*  limitations under the License.                                           */
/*                                                                           */
/*  You should have received a copy of the Apache-2.0 license                */
/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
/**@file   probdata_rpa.c
 * @brief  Problem data for ringpacking problem
 * @author Benjamin Mueller
 *
 * This file handles the main problem data used in that project. For more details see \ref RINGPACKING_PROBLEMDATA page.
 *
 * @page RINGPACKING_PROBLEMDATA Main problem data
 *
 * The problem data is accessible in all plugins. The function SCIPgetProbData() returns the pointer to that
 * structure. We use this data structure to store all the information of the ringpacking problem. Since this structure is
 * not visible in the other plugins, we implemented setter and getter functions to access most data.
 *
 * The function SCIPprobdataCreate(), which is called in the \ref reader_bpa.c "reader plugin" after the input file was
 * parsed, initializes the problem data structure. Afterwards, the problem is setup in SCIPprobdataSetupProblem. For this,
 * it enumerates all dominating circular patterns, selects a set of initial rectangular patterns and creates the
 * corresponding variables and constraints. Note that the pattern constraints have to have the
 * <code>modifiable</code>-flag set to TRUE. This is necessary to tell the solver that these constraints are not
 * completed yet. This means, during the search new variables/patterns might be added. The solver needs this information
 * because certain reductions are not allowed.
 *
 * A list of all interface methods can be found in probdata_binpacking.h.
 */
 
/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
 
#include "scip/scip.h"
#include "scip/scipdefplugins.h"
#include "probdata_rpa.h"
#include "pricer_rpa.h"
 
#include <string.h>
#include <math.h>
 
/* properties of the ringpacking statistics table */
#define TABLE_NAME_RPA                       "ringpacking"
#define TABLE_DESC_RPA                       "ringpacking statistics"
#define TABLE_POSITION_RPA                   12500                  /**< the position of the statistics table */
#define TABLE_EARLIEST_STAGE_RPA             SCIP_STAGE_TRANSFORMED /**< output of the statistics table is only printed from this stage onwards */
 
#ifndef M_PI
#define M_PI           3.141592653589793238462643
#endif
 
/** @brief Problem data which is accessible in all places
 *
 * This problem data is used to store the input of the ringpacking, all variables which are created, and all
 * constraints.
 */
struct SCIP_ProbData
{
   int*                  demands;            /**< array of demands */
   SCIP_Real*            rints;              /**< internal radii of each ring */
   SCIP_Real*            rexts;              /**< external radii of each ring */
   int                   ntypes;             /**< number of different types */
 
   SCIP_Real             width;              /**< height of each rectangle */
   SCIP_Real             height;             /**< width of each rectangle */
 
   SCIP_CONS**           patternconss;       /**< pattern constraints for each type */
 
   /* circular pattern data */
   SCIP_PATTERN**        cpatterns;          /**< array containing all circular patterns */
   SCIP_VAR**            cvars;              /**< variables corresponding to circular patterns */
   int                   ncpatterns;         /**< total number of circular patterns */
   int                   cpatternsize;       /**< size of cpatterns and cvars array */
 
   /* rectangular pattern data */
   SCIP_PATTERN**        rpatterns;          /**< array containing all rectangular patterns */
   SCIP_VAR**            rvars;              /**< variables corresponding to rectangular patterns */
   int                   nrpatterns;         /**< total number of rectangular patterns */
   int                   rpatternsize;       /**< size of rpatterns and rvars array */
 
   /* variables for statistics */
   int                   ncppatternsunknownbeg;/**< number of unknown circular patterns after enumeration step */
   SCIP_Real             enumtime;           /**< time spend for enumerating circular patterns */
   SCIP_Bool             isdualinvalid;      /**< whether the following reported dual bounds are valid */
   SCIP_Real             dualbound;          /**< valid dual bound for RCPP instance */
 
   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generator */
 
 
   /* parameters */
   SCIP_Real             nlptilimsoft;       /**< soft time limit for verification NLP */
   SCIP_Real             heurtilimsoft;      /**< soft time limit for verification heuristic */
   SCIP_Real             totaltilimsoft;     /**< soft time limit for enumerating circular patterns */
   SCIP_Longint          nlpnodelimsoft;     /**< soft node limit for verification NLP */
   int                   heuriterlimsoft;    /**< soft iteration limit for verification heuristic */
};
 
 
/**@name Local methods
 *
 * @{
 */
 
/** auxiliary function to create problem data;
 *
 * @note captures patterns and corresponding variables
 */
static
SCIP_RETCODE probdataCreate(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA**       probdata,           /**< pointer to problem data */
   SCIP_CONS**           patternconss,       /**< pattern constraints */
   SCIP_PATTERN**        cpatterns,          /**< circular patterns */
   SCIP_VAR**            cvars,              /**< variables corresponding to circular patterns */
   int                   ncpatterns,         /**< total number of circular patterns */
   SCIP_PATTERN**        rpatterns,          /**< rectangular patterns */
   SCIP_VAR**            rvars,              /**< variables corresponding to rectangular patterns */
   int                   nrpatterns,         /**< total number of rectangular patterns */
   int*                  demands,            /**< array containing the demands */
   SCIP_Real*            rints,              /**< interal radii of each ring */
   SCIP_Real*            rexts,              /**< external radii of each ring */
   int                   ntypes,             /**< number of different types */
   SCIP_Real             width,              /**< width of each rectangle */
   SCIP_Real             height              /**< height of each rectangle */
   )
{
   int i;
 
   assert(probdata != NULL);
   assert(demands != NULL);
   assert(rints != NULL);
   assert(rexts != NULL);
   assert(ntypes > 0);
   assert(height > 0.0);
   assert(width >= height);
 
   /* allocate memory */
   SCIP_CALL( SCIPallocBlockMemory(scip, probdata) );
   BMSclearMemory(*probdata);
 
   if( ncpatterns > 0 )
   {
      SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->cvars, cvars, ncpatterns) );
      SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->cpatterns, cpatterns, ncpatterns) );
      (*probdata)->ncpatterns = ncpatterns;
      (*probdata)->cpatternsize = ncpatterns;
 
      /* capture circular patterns */
      for( i = 0; i < ncpatterns; ++i )
         SCIPpatternCapture(cpatterns[i]);
   }
 
   if( nrpatterns > 0 )
   {
      SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->rvars, rvars, nrpatterns) );
      SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->rpatterns, rpatterns, nrpatterns) );
      (*probdata)->nrpatterns = nrpatterns;
      (*probdata)->rpatternsize = nrpatterns;
 
      /* capture rectangular patterns */
      for( i = 0; i < nrpatterns; ++i )
         SCIPpatternCapture(rpatterns[i]);
   }
 
   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->demands, demands, ntypes) );
   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->rints, rints, ntypes) );
   SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->rexts, rexts, ntypes) );
 
   /* copy pattern constraints if available, otherwise allocate enough memory */
   if( patternconss != NULL )
   {
      SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*probdata)->patternconss, patternconss, ntypes) );
   }
   else
   {
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*probdata)->patternconss, ntypes) );
      BMSclearMemoryArray((*probdata)->patternconss, ntypes);
   }
 
   /* create random number generator */
   SCIP_CALL( SCIPcreateRandom(scip, &(*probdata)->randnumgen, 0, TRUE) );
 
   (*probdata)->ntypes = ntypes;
   (*probdata)->width = width;
   (*probdata)->height = height;
 
   return SCIP_OKAY;
}
 
/** auxiliary function to free problem data */
static
SCIP_RETCODE probdataFree(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA**       probdata            /**< pointer to release the probdata */
   )
{
   int i;
 
   assert(scip != NULL);
   assert(probdata != NULL);
   assert(*probdata != NULL);
 
   /* free random number generator */
   if( (*probdata)->randnumgen != NULL )
   {
      SCIPfreeRandom(scip, &(*probdata)->randnumgen);
   }
 
   /* release pattern constraints */
   if( (*probdata)->patternconss != NULL )
   {
      for( i = 0; i < SCIPprobdataGetNTypes(*probdata); ++i )
      {
         if( (*probdata)->patternconss[i] != NULL )
         {
            SCIP_CALL( SCIPreleaseCons(scip, &(*probdata)->patternconss[i]));
         }
      }
      SCIPfreeBlockMemoryArray(scip, &(*probdata)->patternconss, (*probdata)->ntypes);
   }
 
   /* release circular patterns */
   for( i = 0; i < (*probdata)->ncpatterns; ++i )
   {
      assert((*probdata)->cpatterns[i] != NULL);
 
      if( (*probdata)->cvars[i] != NULL )
      {
         SCIP_CALL( SCIPreleaseVar(scip, &(*probdata)->cvars[i]) );
      }
 
      SCIPpatternRelease(scip, &(*probdata)->cpatterns[i]);
   }
 
   SCIPfreeBlockMemoryArrayNull(scip, &(*probdata)->cpatterns, (*probdata)->cpatternsize);
   SCIPfreeBlockMemoryArrayNull(scip, &(*probdata)->cvars, (*probdata)->cpatternsize);
 
   /* release rectangular patterns */
   for( i = 0; i < (*probdata)->nrpatterns; ++i )
   {
      assert((*probdata)->rpatterns[i] != NULL);
      assert((*probdata)->rvars[i] != NULL);
 
      if( (*probdata)->rvars[i] != NULL )
      {
         SCIP_CALL( SCIPreleaseVar(scip, &(*probdata)->rvars[i]) );
      }
 
      SCIPpatternRelease(scip, &(*probdata)->rpatterns[i]);
   }
 
   SCIPfreeBlockMemoryArrayNull(scip, &(*probdata)->rpatterns, (*probdata)->rpatternsize);
   SCIPfreeBlockMemoryArrayNull(scip, &(*probdata)->rvars, (*probdata)->rpatternsize);
 
   SCIPfreeBlockMemoryArray(scip, &(*probdata)->rexts, (*probdata)->ntypes);
   SCIPfreeBlockMemoryArray(scip, &(*probdata)->rints, (*probdata)->ntypes);
   SCIPfreeBlockMemoryArray(scip, &(*probdata)->demands, (*probdata)->ntypes);
 
   SCIPfreeBlockMemory(scip, probdata);
   SCIP_CALL( SCIPsetProbData(scip, NULL) );
 
   return SCIP_OKAY;
}
 
/** counts the number of circular patterns with a given packable status */
static
int getNCPatterns(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PACKABLE         status              /**< packable status */
   )
{
   int count = 0;
   int p;
 
   assert(probdata != NULL);
 
   for( p = 0; p < probdata->ncpatterns; ++p )
   {
      if( SCIPpatternGetPackableStatus(probdata->cpatterns[p]) == status )
         ++count;
   }
 
   return count;
}
 
/** ensures a minimum size of the pattern and variable arrays */
static
SCIP_RETCODE ensureSize(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PATTERNTYPE      type,               /**< pattern type */
   int                   size                /**< required size */
   )
{
   int newsize;
 
   assert(probdata != NULL);
   assert(size > 0);
 
   if( type == SCIP_PATTERNTYPE_CIRCULAR && size > probdata->cpatternsize )
   {
      newsize = MAX(size, 2 * probdata->cpatternsize);
 
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &probdata->cpatterns, probdata->cpatternsize, newsize) );
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &probdata->cvars, probdata->cpatternsize, newsize) );
      probdata->cpatternsize = newsize;
   }
   else if( type == SCIP_PATTERNTYPE_RECTANGULAR && size > probdata->rpatternsize )
   {
      newsize = MAX(size, 2 * probdata->rpatternsize);
 
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &probdata->rpatterns, probdata->rpatternsize, newsize) );
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &probdata->rvars, probdata->rpatternsize, newsize) );
      probdata->rpatternsize = newsize;
   }
 
   return SCIP_OKAY;
}
 
/** create variables for all existing circular and rectangular patterns */
static
SCIP_RETCODE createPatternVars(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   SCIP_VAR* var;
   SCIP_PATTERN* pattern;
   char name[SCIP_MAXSTRLEN];
   int k;
 
   assert(probdata != NULL);
   assert(probdata->ncpatterns > 0);
   assert(probdata->nrpatterns > 0);
 
   /* create variables for circular patterns */
   for( k = 0; k < probdata->ncpatterns; ++k )
   {
      SCIP_Real ub;
      int type;
      int i;
 
      pattern = probdata->cpatterns[k];
      assert(pattern != NULL);
 
      type = SCIPpatternGetCircleType(pattern);
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "c%d", type);
      ub = (SCIP_Real)SCIPprobdataGetDemands(probdata)[type];
 
      /* create variable name */
      for( i = 0; i < SCIPpatternGetNElemens(pattern); ++i )
      {
         char strtmp[SCIP_MAXSTRLEN];
         int elemtype = SCIPpatternGetElementType(pattern, i);
         (void) SCIPsnprintf(strtmp, SCIP_MAXSTRLEN, "_%d", elemtype);
         (void) strcat(name, strtmp);
      }
 
      /* create variable */
      SCIP_CALL( SCIPcreateVarBasic(scip, &var, name, 0.0, ub, 0.0, SCIP_VARTYPE_INTEGER) );
      SCIP_CALL( SCIPaddVar(scip, var) );
 
      /* store variables in problem data */
      probdata->cvars[k] = var;
   }
 
   /* create variables for rectangular patterns */
   for( k = 0; k < probdata->nrpatterns; ++k )
   {
      int i;
 
      pattern = probdata->rpatterns[k];
      assert(pattern != NULL);
 
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "r");
 
      /* create variable name */
      for( i = 0; i < SCIPpatternGetNElemens(pattern); ++i )
      {
         char strtmp[SCIP_MAXSTRLEN];
         int elemtype = SCIPpatternGetElementType(pattern, i);
         (void) SCIPsnprintf(strtmp, SCIP_MAXSTRLEN, "_%d", elemtype);
         (void) strcat(name, strtmp);
      }
 
      /* create variable */
      SCIP_CALL( SCIPcreateVarBasic(scip, &var, name, 0.0, SCIPinfinity(scip), 1.0, SCIP_VARTYPE_INTEGER) );
      SCIP_CALL( SCIPaddVar(scip, var) );
 
      /* store variables in problem data */
      probdata->rvars[k] = var;
   }
 
   return SCIP_OKAY;
}
 
/** upper bound on the number of circles of a single type that fit into a circular pattern of a given type */
static
int maxCircles(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   int                   type,               /**< type of the circular pattern */
   int                   elemtype            /**< type of element to pack */
   )
{
   SCIP_Real _rint;
   SCIP_Real rext;
   SCIP_Real rintscaled;
   int demand;
   int n;
 
   assert(type >= 0 && type < SCIPprobdataGetNTypes(probdata));
   assert(elemtype >= 0 && elemtype < SCIPprobdataGetNTypes(probdata));
 
   _rint = SCIPprobdataGetRints(probdata)[type];
   rext = SCIPprobdataGetRexts(probdata)[elemtype];
   demand = SCIPprobdataGetDemands(probdata)[elemtype];
 
   /* volume-bsaed bound */
   n = MIN(demand, (int) SCIPceil(scip, SQR(_rint) / SQR(rext)));
 
   if( n <= 1 )
      return 1;
 
   /* use proven bounds on the density */
   rintscaled = _rint / rext;
   assert(rintscaled >= 1.0);
 
   if( SCIPisLT(scip, rintscaled, 2.0) )
      return MIN(1, n);
   else if( SCIPisLT(scip, rintscaled, 2.1547005383792515) )
      return MIN(2, n);
   else if( SCIPisLT(scip, rintscaled, 2.414213562373095) )
      return MIN(3, n);
   else if( SCIPisLT(scip, rintscaled, 2.7013016167040798) )
      return MIN(4, n);
   else if( SCIPisLT(scip, rintscaled, 3.0) )
      return MIN(5, n);
   else if( SCIPisLT(scip, rintscaled, 3.3047648709624866) )
      return MIN(7, n); /* note that here is a jump and 7 is correct */
   else if( SCIPisLT(scip, rintscaled, 3.613125929752753) )
      return MIN(8, n);
 
   return n;
}
 
/** helper function to compare two patterns; returns
 *
 *   -1 if p dominates q
 *   +1 if q dominates p
 *    0 otherwise
 */
static
int isPatternDominating(
   SCIP_PATTERN*         p,                  /**< pattern */
   SCIP_PATTERN*         q,                  /**< pattern */
   int*                  count,              /**< array for counting elements of patterns */
   int                   ntypes              /**< total number of types */
   )
{
   SCIP_Bool pdomq;
   SCIP_Bool qdomp;
   int i;
 
   /* patterns can only dominate each other if they have the same type */
   if(SCIPpatternGetCircleType(p) != SCIPpatternGetCircleType(q) )
      return 0;
 
   /* reset count array */
   BMSclearMemoryArray(count, ntypes);
 
   /* increase array entry for each element in p */
   for( i = 0; i < SCIPpatternGetNElemens(p); ++i )
   {
      int t = SCIPpatternGetElementType(p, i);
      count[t] += 1;
   }
 
   /* decrease array entry for each element in q */
   for( i = 0; i < SCIPpatternGetNElemens(q); ++i )
   {
      int t = SCIPpatternGetElementType(q, i);
      count[t] -= 1;
   }
 
   pdomq = TRUE;
   qdomp = TRUE;
 
   for( i = 0; i < ntypes && (pdomq || qdomp); ++i )
   {
      if( count[i] < 0 )
         pdomq = FALSE;
      else if( count[i] > 0 )
         qdomp = FALSE;
   }
 
   if( pdomq && (SCIPpatternGetPackableStatus(p) == SCIP_PACKABLE_YES || SCIPpatternGetPackableStatus(q) == SCIP_PACKABLE_UNKNOWN) )
      return -1;
   else if( qdomp && (SCIPpatternGetPackableStatus(q) == SCIP_PACKABLE_YES || SCIPpatternGetPackableStatus(p) == SCIP_PACKABLE_UNKNOWN) )
      return 1;
   return 0;
}
 
/** filter dominated patterns */
static
SCIP_RETCODE filterPatterns(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   SCIP_PATTERN** cpatterns;
   SCIP_Bool* deleted;
   int* count;
   int ncpatterns;
   int i;
 
   SCIP_CALL( SCIPallocBufferArray(scip, &count, SCIPprobdataGetNTypes(probdata)) );
   SCIP_CALL( SCIPallocBufferArray(scip, &cpatterns, probdata->ncpatterns) );
   SCIP_CALL( SCIPallocBufferArray(scip, &deleted, probdata->ncpatterns) );
   BMSclearMemoryArray(deleted, probdata->ncpatterns);
 
   for( i = 0; i < probdata->ncpatterns - 1; ++i )
   {
      SCIP_PATTERN* p = probdata->cpatterns[i];
      int j;
 
      if( deleted[i] )
         continue;
 
      for( j = i + 1; j < probdata->ncpatterns; ++j )
      {
         SCIP_PATTERN* q = probdata->cpatterns[j];
         int res;
 
         if( deleted[j] )
            continue;
 
         res = isPatternDominating(p, q, count, SCIPprobdataGetNTypes(probdata));
 
         /* p dominates q */
         if( res == -1 )
            deleted[j] = TRUE;
         else if( res == 1 ) /* q dominates p */
            deleted[i] = TRUE;
      }
   }
 
   /* remove filtered patterns */
   ncpatterns = 0;
   for( i = 0; i < probdata->ncpatterns; ++i )
   {
      if( deleted[i] )
      {
         SCIPpatternRelease(scip, &probdata->cpatterns[i]);
      }
      else
      {
         cpatterns[ncpatterns] = probdata->cpatterns[i];
         ++ncpatterns;
      }
   }
   assert(ncpatterns > 0);
 
   BMScopyMemoryArray(probdata->cpatterns, cpatterns, ncpatterns);
   probdata->ncpatterns = ncpatterns;
 
   /* free memory */
   SCIPfreeBufferArray(scip, &deleted);
   SCIPfreeBufferArray(scip, &cpatterns);
   SCIPfreeBufferArray(scip, &count);
 
   return SCIP_OKAY;
}
 
/** enumerates all circular patterns for a given type */
static
SCIP_RETCODE enumeratePatterns(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PATTERN*         pattern,            /**< pattern (passed for performance reasons) */
   int*                  ms,                 /**< maximum number of elements for each type (passed for performance reasons) */
   int*                  nselected,          /**< number of selected elements for each type (passed for performance reasons) */
   SCIP_Real             nlptilim,           /**< time limit for each NLP verification */
   SCIP_Real             heurtilim,          /**< time limit for each call of the heuristics */
   SCIP_Longint          nlpnodelim,         /**< node limit for each NLP verification */
   int                   heuriterlim,        /**< iteration limit for each call of the heuristics */
   SCIP_Real*            timeleft            /**< pointer to update the remaining time for the enumeration */
   )
{
   SCIP_Real* rexts;
   SCIP_Real* _rints;
   SCIP_Real maxvolume;
   SCIP_Real volume;
   int ntypes;
   int type;
   int lasttype;
 
   assert(ms != NULL);
   assert(pattern != NULL);
   assert(timeleft != NULL);
 
   type = SCIPpatternGetCircleType(pattern);
   assert(type >= 0 && type < SCIPprobdataGetNTypes(probdata));
 
   /* get problem data */
   rexts = SCIPprobdataGetRexts(probdata);
   _rints = SCIPprobdataGetRints(probdata);
   ntypes = SCIPprobdataGetNTypes(probdata);
   lasttype = ntypes -1;
   volume = 0.0;
   maxvolume = SQR(_rints[SCIPpatternGetCircleType(pattern)]) * M_PI; /*lint !e666*/
 
   /* main loop */
   while( TRUE )
   {
      SCIP_Real timelim;
      int t = lasttype;
 
      /* reset packable status */
      SCIPpatternSetPackableStatus(pattern, SCIP_PACKABLE_UNKNOWN);
 
      SCIPdebugMsg(scip, "volume = %g <= %g\n", volume, maxvolume);
 
      {
         int j;
         SCIPdebugMsg(scip, "verify c%d", type);
 
         for( j = 0; j < SCIPpatternGetNElemens(pattern); ++j )
            SCIPdebugMsgPrint(scip, "_%d", SCIPpatternGetElementType(pattern, j));
         SCIPdebugMsgPrint(scip, "\n");
      }
 
      /* check volume */
      if( SCIPisLE(scip, volume, maxvolume) )
      {
         /*
          * try to verify with heuristic
          */
 
         /* compute time limit */
         timelim = MIN(heurtilim, *timeleft);
 
         /* verify pattern */
         *timeleft += SCIPgetTotalTime(scip);
         SCIP_CALL( SCIPverifyCircularPatternHeuristic(scip, probdata, pattern, timelim, heuriterlim) );
         *timeleft -= SCIPgetTotalTime(scip);
 
         /*
          * try to verify with NLP
          */
         if( SCIPpatternGetPackableStatus(pattern) == SCIP_PACKABLE_UNKNOWN )
         {
            /* compute time limit */
            timelim = MIN(*timeleft, nlptilim);
 
            /* verify pattern */
            *timeleft += SCIPgetTotalTime(scip);
            SCIP_CALL( SCIPverifyCircularPatternNLP(scip, probdata, pattern, timelim, nlpnodelim) );
            *timeleft -= SCIPgetTotalTime(scip);
         }
 
         /* pattern is not packable -> don't add more elements */
         if( SCIPpatternGetPackableStatus(pattern) == SCIP_PACKABLE_NO )
         {
            SCIPpatternRemoveLastElements(pattern, nselected[t]);
            volume -= SQR(rexts[t]) * M_PI * nselected[t];
            nselected[t] = 0;
            --t;
         }
         /* otherwise add the pattern (and hope for filtering) */
         else
         {
            SCIP_CALL( SCIPprobdataAddVar(scip, probdata, pattern, NULL) );
         }
      }
 
      /* update selection */
      while( t > type && (nselected[t] == ms[t] || SCIPisGT(scip, volume, maxvolume)) )
      {
         SCIPpatternRemoveLastElements(pattern, nselected[t]);
         volume -= SQR(rexts[t]) * M_PI * nselected[t];
         nselected[t] = 0;
         t--;
      }
 
      /* check termination criterion */
      if( t == type )
         break;
 
      /* add element of type i to the pattern */
      assert(nselected[t] < ms[t]);
      ++(nselected[t]);
      volume += SQR(rexts[t]) * M_PI;
      SCIP_CALL( SCIPpatternAddElement(pattern, t, SCIP_INVALID, SCIP_INVALID) );
   }
 
   assert(SCIPpatternGetNElemens(pattern) == 0);
   assert(SCIPisZero(scip, volume));
 
   return SCIP_OKAY;
}
 
/** auxiliary function to setup the master problem */
static
SCIP_RETCODE setupProblem(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   char name[SCIP_MAXSTRLEN];
   SCIP_Real* rexts;
   SCIP_Real* rints;
   int* demands;
   SCIP_Real dualbound;
   SCIP_Real minrext;
   SCIP_Real volume;
   int ntypes;
   int p;
   int t;
 
   assert(probdata != NULL);
   assert(SCIPprobdataGetNTypes(probdata) > 0);
 
   /* set objective sense; tell SCIP that the objective will be always integral */
   SCIP_CALL( SCIPsetObjsense(scip, SCIP_OBJSENSE_MINIMIZE) );
   SCIP_CALL( SCIPsetObjIntegral(scip) );
 
   /* get problem data */
   ntypes = SCIPprobdataGetNTypes(probdata);
   rexts = SCIPprobdataGetRexts(probdata);
   rints = SCIPprobdataGetRints(probdata);
   demands = SCIPprobdataGetDemands(probdata);
 
   /* compute all non-dominated circular patterns */
   probdata->enumtime -= SCIPgetTotalTime(scip);
   SCIP_CALL( SCIPprobdataEnumeratePatterns(scip, probdata, probdata->nlptilimsoft, probdata->heurtilimsoft,
      probdata->totaltilimsoft, probdata->nlpnodelimsoft, probdata->heuriterlimsoft) );
   probdata->enumtime += SCIPgetTotalTime(scip);
   probdata->ncppatternsunknownbeg = getNCPatterns(scip, probdata, SCIP_PACKABLE_UNKNOWN);
 
   SCIPinfoMessage(scip, NULL, "+++++++++++++ starting with |CP|=%d\n", probdata->ncpatterns);
 
   /* create initial rectangular patterns */
   for( t = 0; t < ntypes; ++t )
   {
      SCIP_PATTERN* pattern;
 
      /* create a pattern containing a single circle of type t; set position of the circle to the left-bottom */
      SCIP_CALL( SCIPpatternCreateRectangular(scip, &pattern) );
      SCIP_CALL( SCIPpatternAddElement(pattern, t, rexts[t], rexts[t]) );
      SCIPpatternSetPackableStatus(pattern, SCIP_PACKABLE_YES);
 
      /* add and release pattern */
      SCIP_CALL( SCIPprobdataAddVar(scip, probdata, pattern, NULL) );
      SCIPpatternRelease(scip, &pattern);
   }
 
   /* create variables for all existing patterns */
   SCIP_CALL( createPatternVars(scip, probdata) );
 
   /* create demand constraints */
   for( t = 0; t < ntypes; ++t )
   {
      SCIP_CONS* cons;
 
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "demand_%d", t);
      SCIP_CALL( SCIPcreateConsBasicLinear(scip, &cons, name, 0, NULL, NULL, (SCIP_Real)demands[t], SCIPinfinity(scip) ) );
 
      for( p = 0; p < probdata->ncpatterns; ++p )
      {
         SCIP_PATTERN* pattern;
         SCIP_VAR* var;
 
         pattern = probdata->cpatterns[p];
         assert(pattern != NULL);
         assert(SCIPpatternGetPatternType(pattern) == SCIP_PATTERNTYPE_CIRCULAR);
 
         var = probdata->cvars[p];
         assert(var != NULL);
 
         /* add coefficient to the pattern if the pattern is of type t */
         if(SCIPpatternGetCircleType(pattern) == t )
         {
            SCIP_CALL( SCIPaddCoefLinear(scip, cons, var, 1.0) );
         }
      }
 
      /* add and release constraint */
      SCIP_CALL( SCIPaddCons(scip, cons) );
      SCIP_CALL( SCIPreleaseCons(scip, &cons) );
   }
 
   /* create pattern constraints */
   for( t = 0; t < ntypes; ++t )
   {
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "patterncons_%d", t);
      SCIP_CALL( SCIPcreateConsBasicLinear(scip, &probdata->patternconss[t], name, 0, NULL, NULL, 0.0,
            SCIPinfinity(scip) ) );
 
      /* declare constraint modifiable for adding variables during pricing */
      SCIP_CALL( SCIPsetConsModifiable(scip, probdata->patternconss[t], TRUE) );
      SCIP_CALL( SCIPaddCons(scip, probdata->patternconss[t]) );
   }
 
   /* add coefficients for circular patterns */
   for( p = 0; p < probdata->ncpatterns; ++p )
   {
      SCIP_PATTERN* pattern = probdata->cpatterns[p];
      SCIP_VAR* var = probdata->cvars[p];
      int type;
 
      assert(pattern != NULL);
      assert(SCIPpatternGetPatternType(pattern) == SCIP_PATTERNTYPE_CIRCULAR);
      assert(var != NULL);
 
      type = SCIPpatternGetCircleType(pattern);
      assert(type >= 0 && type < ntypes);
 
      /* - z_C */
      SCIP_CALL( SCIPaddCoefLinear(scip, probdata->patternconss[type], var, -1.0) );
 
      for( t = 0; t < ntypes; ++t )
      {
         int nelems = SCIPpatternCountElements(pattern, t);
 
         if( nelems > 0 )
         {
            /* + P_t z_C */
            SCIP_CALL( SCIPaddCoefLinear(scip, probdata->patternconss[t], var, (SCIP_Real)nelems) );
         }
      }
   }
 
   /* add coefficients for rectangular patterns */
   for( p = 0; p < probdata->nrpatterns; ++p )
   {
      SCIP_PATTERN* pattern = probdata->rpatterns[p];
      SCIP_VAR* var = probdata->rvars[p];
 
      assert(pattern != NULL);
      assert(SCIPpatternGetPatternType(pattern) == SCIP_PATTERNTYPE_RECTANGULAR);
      assert(var != NULL);
 
      for( t = 0; t < ntypes; ++t )
      {
         int nelems = SCIPpatternCountElements(pattern, t);
 
         if( nelems > 0 )
         {
            /* + P_t z_P */
            SCIP_CALL( SCIPaddCoefLinear(scip, probdata->patternconss[t], var, (SCIP_Real)nelems) );
         }
      }
   }
 
   /* compute an initial dual bound by considering the volume of all rings */
   minrext = rexts[ntypes-1];
   volume = 0.0;
   for( t = 0; t < ntypes; ++t )
   {
      SCIP_Real vol;
 
      /* consider ring as circle if there is no ring with a smaller radius than than inner one */
      if( SCIPisFeasLT(scip, rints[t], minrext) )
         vol = M_PI * SQR(rexts[t]);
      else
         vol = M_PI * (SQR(rexts[t]) - SQR(rints[t]));
 
      volume += vol * demands[t];
   }
 
   /* update initial dual bound */
   dualbound = SCIPfeasCeil(scip, volume / (SCIPprobdataGetWidth(probdata) * SCIPprobdataGetHeight(probdata)));
   SCIP_CALL( SCIPupdateLocalDualbound(scip, dualbound) );
   SCIPinfoMessage(scip, NULL, "+++++++++++++ volume-based bound = ceil(%g / %g) = %g\n", volume,
      SCIPprobdataGetWidth(probdata) * SCIPprobdataGetHeight(probdata), dualbound);
   SCIPprobdataUpdateDualbound(scip, probdata, dualbound);
 
   return SCIP_OKAY;
}
 
/** output method of statistics table to output file stream 'file' */
static
SCIP_DECL_TABLEOUTPUT(tableOutputRpa)
{ /*lint --e{715}*/
   SCIP_PROBDATA* probdata;
   SCIP_Real* rexts;
   SCIP_Real* rints;
   SCIP_Real dualbound;
   SCIP_Real maxrint;
   SCIP_Real minrext;
   int* demands;
   int ntypes;
   int nrings;
   int t;
 
   probdata = SCIPgetProbData(scip);
   assert(probdata != NULL);
 
   ntypes = SCIPprobdataGetNTypes(probdata);
   demands = SCIPprobdataGetDemands(probdata);
   rexts = SCIPprobdataGetRexts(probdata);
   rints = SCIPprobdataGetRints(probdata);
   nrings = 0;
   maxrint = 0.0;
   minrext = SCIPinfinity(scip);
 
   /* use global dual bound if it is still valid */
   if( !probdata->isdualinvalid )
   {
      assert(SCIPisGE(scip, SCIPgetDualbound(scip), probdata->dualbound));
      dualbound = SCIPgetDualbound(scip);
   }
   else
      dualbound = probdata->dualbound;
 
   /* count the number of rings */
   for( t = 0; t < ntypes; ++t )
   {
      nrings += demands[t];
      maxrint = MAX(maxrint, rints[t]);
      minrext = MIN(minrext, rexts[t]);
   }
 
   SCIPinfoMessage(scip, file, "Ringpacking        : %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
      "dual", "ntypes", "nrings", "width", "height", "CP", "CP_unk", "CP_unk_end" ,"CP_infeas", "RP", "enumtime", "radiiratio");
 
   SCIPinfoMessage(scip, file, "  %-17s:", "");
   SCIPinfoMessage(scip, file, " %10.2f", dualbound);
   SCIPinfoMessage(scip, file, " %10d", ntypes);
   SCIPinfoMessage(scip, file, " %10d", nrings);
   SCIPinfoMessage(scip, file, " %10.2f", SCIPprobdataGetWidth(probdata));
   SCIPinfoMessage(scip, file, " %10.2f", SCIPprobdataGetHeight(probdata));
   SCIPinfoMessage(scip, file, " %10d", probdata->ncpatterns);
   SCIPinfoMessage(scip, file, " %10d", probdata->ncppatternsunknownbeg);
   SCIPinfoMessage(scip, file, " %10d", getNCPatterns(scip, probdata, SCIP_PACKABLE_UNKNOWN));
   SCIPinfoMessage(scip, file, " %10d", getNCPatterns(scip, probdata, SCIP_PACKABLE_NO));
   SCIPinfoMessage(scip, file, " %10d", probdata->nrpatterns);
   SCIPinfoMessage(scip, file, " %10.2f", probdata->enumtime);
   SCIPinfoMessage(scip, file, " %10.1f", maxrint / minrext);
   SCIPinfoMessage(scip, file, "\n");
 
   return SCIP_OKAY;
}
 
/** auxiliary function to update the best known candidate */
static
void updateBestCandidate(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Real*            xs,                 /**< x-coordinates of packed elements */
   SCIP_Real*            ys,                 /**< y-coordinates of packed elements */
   SCIP_Real*            rexts,              /**< radii of packed elements */
   SCIP_Real             rext,               /**< radii of element that should be packed */
   SCIP_Real             rbounding,          /**< inner radius of bounding circle (ignored for rectangular patterns) */
   SCIP_Real             wbounding,          /**< width of bounding rectangular (ignored for circular patterns) */
   SCIP_Real             hbounding,          /**< height of bounding rectangular (ignored for circular patterns) */
   SCIP_Real             rmax,               /**< maximum radius of elements in the pattern */
   SCIP_PATTERNTYPE      patterntype,        /**< pattern type */
   SCIP_Bool*            ispacked,           /**< array indicating which elements are already packed */
   int*                  elements,           /**< the order of the elements in the pattern */
   int                   nelements,          /**< the total number of elements */
   SCIP_Real*            bestx,              /**< buffer to update best x-coordinate */
   SCIP_Real*            besty,              /**< buffer to update best y-coordinate */
   SCIP_Real             x,                  /**< x-coordinate of a candidate point */
   SCIP_Real             y,                  /**< y-coordinate of a candidate point */
   int                   ncalls              /**< total number of calls of the packing heuristic */
   )
{
   SCIP_Real threshold;
   SCIP_Bool isoverthreshold;
   int i;
 
   /* candidate is not valid -> skip */
   if( x == SCIP_INVALID || y == SCIP_INVALID ) /*lint !e777*/
      return;
 
   /* check whether there is an intersection with the boundary */
   if( patterntype == SCIP_PATTERNTYPE_CIRCULAR )
   {
      if( SCIPisGT(scip, x*x + y*y, SQR(rbounding - rext)) )
         return;
   }
   else
   {
      if( SCIPisLT(scip, x, rext) || SCIPisGT(scip, x, wbounding - rext)
         || SCIPisLT(scip, y, rext) || SCIPisGT(scip, y, hbounding - rext) )
         return;
   }
 
   /* check whether circle intersects other circles */
   for( i = 0; i < nelements; ++i )
   {
      SCIP_Real dist;
 
      /* only consider packed elements */
      if( !ispacked[i] )
         continue;
 
      dist = SQR(x - xs[i]) + SQR(y - ys[i]);
 
      /* check if the distance between mid points is smaller than the sum of the radii */
      if( SCIPisLT(scip, dist, SQR(rext + rexts[elements[i]])) )
         return;
   }
 
   threshold = (patterntype == SCIP_PATTERNTYPE_RECTANGULAR ? wbounding - 2.0*rmax - rext : rbounding - 2.0*rmax - rext);
   isoverthreshold = (ncalls % 2) == 1 && SCIPisGT(scip, *bestx, threshold) && SCIPisGT(scip, x, threshold);
 
   /* check whether the candidate is better than the best known candidate */
   if( *bestx == SCIP_INVALID || *besty == SCIP_INVALID /*lint !e777*/
      || ((!isoverthreshold || SCIPisEQ(scip, y, *besty)) && SCIPisLT(scip, x, *bestx)) /*lint !e777*/
      || ((isoverthreshold || SCIPisEQ(scip, x, *bestx)) && SCIPisLT(scip, y, *besty)) ) /*lint !e777*/
   {
      *bestx = x;
      *besty = y;
   }
}
 
/** auxiliary function for computing a candidate position between a circle and the outer ring */
static
void computePosRingCircle(
   SCIP*                 scip,               /**< SCIP data structure */
   int*                  elements,           /**< types of elements that have been packed */
   int                   nelements,          /**< the total number of elements */
   SCIP_Real*            rexts,              /**< external radii */
   SCIP_Real*            xs,                 /**< x-coordinate of circle */
   SCIP_Real*            ys,                 /**< y-coordinate of circle */
   int                   pos,                /**< position of element in the elements array */
   SCIP_Bool*            ispacked,           /**< array indicating whether an element has been packed already */
   SCIP_Real             rmax,               /**< maximum radius of elements in the pattern */
   SCIP_Real             rbound,             /**< radius of bounding circle */
   SCIP_Real*            bestx,              /**< pointer to store the best x-coordinate */
   SCIP_Real*            besty,              /**< pointer to store the best y-coordinate */
   int                   ncalls              /**< total number of calls of the packing heuristic */
   )
{
   int i;
 
   /* consider already packed patterns */
   for( i = 0; i < nelements; ++i )
   {
      SCIP_Real alpha, a, b, c, h, u, v, n1, n2;
 
      /* only consider packed elements */
      if( !ispacked[i] )
         continue;
 
      c = sqrt(xs[i]*xs[i] + ys[i]*ys[i]);
 
      /* inner ring is too far away from boundary or both rings can not fit */
      if( !SCIPisGE(scip, c + rexts[elements[i]] + 2.0*rexts[elements[pos]], rbound)
         || SCIPisGT(scip, rexts[elements[pos]] + rexts[elements[i]], rbound) )
         continue;
 
      a = rexts[elements[pos]] + rexts[elements[i]];
      b = rbound - rexts[elements[pos]];
 
      /* if a ring is in the center than there are infinitely many solutions; take an arbitrary point */
      if( SCIPisZero(scip, c) )
      {
         updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], rbound, -1.0, -1.0, rmax, SCIP_PATTERNTYPE_CIRCULAR,
            ispacked, elements, nelements, bestx, besty, -rbound + rexts[elements[pos]], 0.0, ncalls);
         updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], rbound, -1.0, -1.0, rmax, SCIP_PATTERNTYPE_CIRCULAR,
            ispacked, elements, nelements, bestx, besty, +rbound - rexts[elements[pos]], 0.0, ncalls);
      }
      else
      {
         assert(c != 0.0);
         alpha = (c*c - b*b + a*a) / (2*c);
 
         if( a*a >= alpha*alpha )
         {
            h = sqrt(MAX(a*a - alpha*alpha, 0.0));
            u = (c - alpha) * xs[i] / c;
            v = (c - alpha) * ys[i] / c;
 
            n1 = SCIPisZero(scip, v) ? 0.0 : h * (v / sqrt(v*v + u*u));
            n2 = SCIPisZero(scip, u) ? 0.0 : h * (-u / sqrt(v*v + u*u));
 
            updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], rbound, -1.0, -1.0, rmax, SCIP_PATTERNTYPE_CIRCULAR,
               ispacked, elements, nelements, bestx, besty, u + n1, v + n2, ncalls);
            updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], rbound, -1.0, -1.0, rmax, SCIP_PATTERNTYPE_CIRCULAR,
               ispacked, elements, nelements, bestx, besty, u - n1, v - n2, ncalls);
         }
      }
   }
}
 
/** auxiliary function for computing trivial candidate positions */
static
void computePosTrivial(
   SCIP*                 scip,               /**< SCIP data structure */
   int*                  elements,           /**< types of elements that have been packed */
   int                   nelements,          /**< the total number of elements */
   SCIP_Real*            rexts,              /**< external radii */
   SCIP_Real*            xs,                 /**< x-coordinate of circle */
   SCIP_Real*            ys,                 /**< y-coordinate of circle */
   int                   pos,                /**< position of element in the elements array */
   SCIP_Bool*            ispacked,           /**< array indicating whether an element has been packed already */
   SCIP_Real             rmax,               /**< maximum radius of elements in the pattern */
   SCIP_Real             rbound,             /**< radius of bounding circle */
   SCIP_Real             width,              /**< width of the rectangle */
   SCIP_Real             height,             /**< height of the rectangle */
   SCIP_PATTERNTYPE      patterntype,        /**< the pattern type (rectangular or circular) */
   SCIP_Real*            bestx,              /**< pointer to store the best x-coordinate */
   SCIP_Real*            besty,              /**< pointer to store the best y-coordinate */
   int                   ncalls              /**< total number of calls of the packing heuristic */
   )
{
   SCIP_Real rext = rexts[elements[pos]];
   int i;
 
   if( patterntype == SCIP_PATTERNTYPE_CIRCULAR )
   {
      SCIP_Real xcands[4] = {-rbound + rext, +rbound - rext, 0.0, 0.0};
      SCIP_Real ycands[4] = {0.0, 0.0, -rbound + rext, +rbound - rext};
 
      for( i = 0; i < 4; ++i )
         updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], rbound, width, height, rmax, patterntype,
            ispacked, elements, nelements, bestx, besty, xcands[i], ycands[i], ncalls);
   }
   else
   {
      SCIP_Real xcands[4] = {rext, width - rext, rext, width - rext};
      SCIP_Real ycands[4] = {rext, rext, height - rext, height - rext};
 
      for( i = 0; i < 4; ++i )
         updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], rbound, width, height, rmax, patterntype,
            ispacked, elements, nelements, bestx, besty, xcands[i], ycands[i], ncalls);
   }
}
 
/** auxiliary function for computing a candidate position between a circle and the rectangle */
static
void computePosRectangleCircle(
   SCIP*                 scip,               /**< SCIP data structure */
   int*                  elements,           /**< types of elements that have been packed */
   int                   nelements,          /**< the total number of elements */
   SCIP_Real*            rexts,              /**< external radii */
   SCIP_Real*            xs,                 /**< x-coordinate of circle */
   SCIP_Real*            ys,                 /**< y-coordinate of circle */
   int                   pos,                /**< position of element in the elements array */
   SCIP_Bool*            ispacked,           /**< array indicating whether an element has been packed already */
   SCIP_Real             rmax,               /**< maximum radius of elements in the pattern */
   SCIP_Real             width,              /**< width of the rectangle */
   SCIP_Real             height,             /**< height of the rectangle */
   SCIP_Real*            bestx,              /**< pointer to store the best x-coordinate */
   SCIP_Real*            besty,              /**< pointer to store the best y-coordinate */
   int                   ncalls              /**< total number of calls of the packing heuristic */
   )
{
   SCIP_Real rext;
   int i;
 
   rext = rexts[elements[pos]];
 
   for( i = 0; i < nelements; ++i )
   {
      SCIP_Real xfix[2] = {rext, width - rext};
      SCIP_Real yfix[2] = {rext, height - rext};
      SCIP_Real Ri;
      int k;
 
      if( !ispacked[i] )
         continue;
 
      Ri = rexts[elements[i]];
 
      /* fix x */
      for( k = 0; k < 2; ++k )
      {
         SCIP_Real alpha = SQR(rext + Ri) - SQR(xfix[k] - xs[i]);
 
         if( alpha < 0.0 )
            continue;
 
         updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], -1.0, width, height, rmax,
            SCIP_PATTERNTYPE_RECTANGULAR, ispacked, elements, nelements, bestx, besty, xfix[k], ys[i] + sqrt(alpha), ncalls);
 
         updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], -1.0, width, height, rmax,
            SCIP_PATTERNTYPE_RECTANGULAR, ispacked, elements, nelements, bestx, besty, xfix[k], ys[i] - sqrt(alpha), ncalls);
      }
 
      /* fix y */
      for( k = 0; k < 2; ++k )
      {
         SCIP_Real alpha = SQR(rext + Ri) - SQR(yfix[k] - ys[i]);
 
         if( alpha < 0.0 )
            continue;
 
         updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], -1.0, width, height, rmax,
            SCIP_PATTERNTYPE_RECTANGULAR, ispacked, elements, nelements, bestx, besty, xs[i] + sqrt(alpha), yfix[k], ncalls);
 
         updateBestCandidate(scip, xs, ys, rexts, rexts[elements[pos]], -1.0, width, height, rmax,
            SCIP_PATTERNTYPE_RECTANGULAR, ispacked, elements, nelements, bestx, besty, xs[i] - sqrt(alpha), yfix[k], ncalls);
      }
   }
}
 
/** auxiliary function for computing a candidate position between two circles */
static
void computePosCircleCircle(
   SCIP*                 scip,               /**< SCIP data structure */
   int*                  elements,           /**< types of elements that have been packed */
   int                   nelements,          /**< the total number of elements */
   SCIP_Real*            rexts,              /**< external radii */
   SCIP_Real*            xs,                 /**< x-coordinate of circle */
   SCIP_Real*            ys,                 /**< y-coordinate of circle */
   int                   pos,                /**< position of element in the elements array */
   SCIP_Bool*            ispacked,           /**< array indicating whether an element has been packed already */
   SCIP_Real             rmax,               /**< maximum radius of elements in the pattern */
   SCIP_Real             rbound,             /**< radius of bounding circle */
   SCIP_Real             width,              /**< width of the rectangle */
   SCIP_Real             height,             /**< height of the rectangle */
   SCIP_PATTERNTYPE      patterntype,        /**< the pattern type (rectangular or circular) */
   SCIP_Real*            bestx,              /**< pointer to store the best x-coordinate */
   SCIP_Real*            besty,              /**< pointer to store the best y-coordinate */
   int                   ncalls              /**< total number of calls of the packing heuristic */
   )
{
   SCIP_Real rext;
   int i;
 
   rext = rexts[elements[pos]];
 
   /* consider all pairs of already packed circles */
   for( i = 0; i < nelements - 1; ++i )
   {
      SCIP_Real alpha, a, b, h, u, v, n1, n2;
      SCIP_Real Ri;
      int j;
 
      if( !ispacked[i] )
         continue;
 
      Ri = rexts[elements[i]];
 
      for( j = i + 1; j < nelements; ++j )
      {
         SCIP_Real Rj;
         SCIP_Real dist;
 
         if( !ispacked[j] )
            continue;
 
         Rj = rexts[elements[j]];
         dist = sqrt(SQR(xs[i] - xs[j]) + SQR(ys[i] - ys[j]));
 
         /* circles are too far away */
         if( SCIPisGE(scip, dist, Ri + Rj + 2.0 * rext) )
            continue;
 
         a = Ri + rext;
         b = Rj + rext;
         assert(dist != 0.0);
         alpha = (dist*dist - b*b + a*a) / (2.0*dist);
         assert(a*a >= alpha*alpha);
         h = sqrt(MAX(a*a - alpha*alpha, 0.0));
         u = xs[i] + (alpha / dist) * (xs[j] - xs[i]);
         v = ys[i] + (alpha / dist) * (ys[j] - ys[i]);
         n1 = h * ((ys[j] - ys[i]) / dist);
         n2 = h * ((xs[i] - xs[j]) / dist);
         assert(n1*n1 + n2*n2 > 0.0);
 
         updateBestCandidate(scip, xs, ys, rexts, rext, rbound, width, height, rmax, patterntype, ispacked, elements,
            nelements, bestx, besty, u + n1, v + n2, ncalls);
         updateBestCandidate(scip, xs, ys, rexts, rext, rbound, width, height, rmax, patterntype, ispacked, elements,
            nelements, bestx, besty, u - n1, v - n2, ncalls);
      }
   }
}
 
/** array to compute the score of each element */
static
void computeScores(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   int*                  elements,           /**< type of each element */
   int                   nelements,          /**< total number of elements */
   SCIP_Real*            scores,             /**< array to store the score of each element */
   int                   iter                /**< iteration round */
   )
{
   SCIP_Real* rexts;
   int i;
 
   rexts = SCIPprobdataGetRexts(probdata);
   assert(rexts != NULL);
 
   for( i = 0; i < nelements; ++i )
   {
      SCIP_Real rext = rexts[elements[i]];
      /* use largest elements first */
      if( iter == 0 )
         scores[i] = rext;
 
      /* use smallest elements first */
      else if( iter == 1 )
         scores[i] = -rext;
 
      /* use [0,1] * radius */
      else if( iter <= 10 )
         scores[i] = SCIPrandomGetReal(probdata->randnumgen, 0.0, 1.0) * rext;
 
      /* use [-1,0] * radius */
      else if( iter <= 20 )
         scores[i] = SCIPrandomGetReal(probdata->randnumgen, -1.0, 0.0) * rext;
 
      /* use a random order */
      else
         scores[i] = SCIPrandomGetReal(probdata->randnumgen, 0.0, 1.0);
   }
}
 
/**@} */
 
/**@name Callback methods of problem data
 *
 * @{
 */
 
/** frees user data of original problem (called when the original problem is freed) */
static
SCIP_DECL_PROBDELORIG(probdelorigRingpacking)
{
   SCIPdebugMessage("free original problem data\n");
   SCIP_CALL( probdataFree(scip, probdata) );
 
   return SCIP_OKAY;
}
 
/** frees user data of transformed problem (called when the transformed problem is freed) */
static
SCIP_DECL_PROBDELTRANS(probdeltransRingpacking)
{
   SCIPdebugMessage("free transformed problem data\n");
   SCIP_CALL( probdataFree(scip, probdata) );
 
   return SCIP_OKAY;
}
 
/** creates user data of transformed problem by transforming the original user problem data
 *  (called after problem was transformed) */
static
SCIP_DECL_PROBTRANS(probtransRingpacking)
{  /*lint --e{715}*/
   /* create transformed problem data */
   SCIP_CALL( probdataCreate(scip, targetdata, sourcedata->patternconss, sourcedata->cpatterns, sourcedata->cvars,
         sourcedata->ncpatterns, sourcedata->rpatterns, sourcedata->rvars, sourcedata->nrpatterns,
         sourcedata->demands, sourcedata->rints, sourcedata->rexts, sourcedata->ntypes,
         sourcedata->width, sourcedata->height) );
 
   /* transform pattern constraints */
   SCIP_CALL( SCIPtransformConss(scip, (*targetdata)->ntypes, (*targetdata)->patternconss,
         (*targetdata)->patternconss) );
 
   /* transform all variables */
   SCIP_CALL( SCIPtransformVars(scip, (*targetdata)->ncpatterns, (*targetdata)->cvars, (*targetdata)->cvars) );
   SCIP_CALL( SCIPtransformVars(scip, (*targetdata)->nrpatterns, (*targetdata)->rvars, (*targetdata)->rvars) );
 
   /* copy statistics to transformed problem data */
   (*targetdata)->ncppatternsunknownbeg = sourcedata->ncppatternsunknownbeg;
   (*targetdata)->enumtime = sourcedata->enumtime;
   (*targetdata)->dualbound = sourcedata->dualbound;
   (*targetdata)->isdualinvalid = sourcedata->isdualinvalid;
 
   return SCIP_OKAY;
}
 
/**@} */
 
/**@name Interface methods
 *
 * @{
 */
 
/** sets up the problem data */
SCIP_RETCODE SCIPprobdataCreate(
   SCIP*                 scip,               /**< SCIP data structure */
   const char*           probname,           /**< problem name */
   int*                  demands,            /**< array containing the demands */
   SCIP_Real*            rints,              /**< internal radii of each ring */
   SCIP_Real*            rexts,              /**< external radii of each ring (assumed to be sorted) */
   int                   ntypes,             /**< number of different types */
   SCIP_Real             width,              /**< width of each rectangle */
   SCIP_Real             height              /**< height of each rectangle */
   )
{
   SCIP_PROBDATA* probdata;
 
#ifndef NDEBUG
   {
      int t;
 
      for( t = 0; t < ntypes -1; ++t )
         assert(rexts[t] >= rexts[t+1]);
   }
#endif
 
   /* create SCIP problem */
   SCIP_CALL( SCIPcreateProbBasic(scip, probname) );
 
   /* create and set problem data */
   SCIP_CALL( probdataCreate(scip, &probdata, NULL, NULL, NULL, 0, NULL, NULL, 0, demands, rints, rexts, ntypes, width,
         height) );
   SCIP_CALL( SCIPsetProbData(scip, probdata) );
 
   SCIP_CALL( SCIPsetProbDelorig(scip, probdelorigRingpacking) );
   SCIP_CALL( SCIPsetProbTrans(scip, probtransRingpacking) );
   SCIP_CALL( SCIPsetProbDeltrans(scip, probdeltransRingpacking) );
 
   /* activate pricer */
   SCIP_CALL( SCIPpricerRpaActivate(scip) );
 
   /* add table output */
   assert(SCIPfindTable(scip, TABLE_NAME_RPA) == NULL);
   SCIP_CALL( SCIPincludeTable(scip, TABLE_NAME_RPA, TABLE_DESC_RPA, TRUE,
         NULL, NULL, NULL, NULL, NULL, NULL, tableOutputRpa, NULL,
         NULL, TABLE_POSITION_RPA, TABLE_EARLIEST_STAGE_RPA) );
 
   return SCIP_OKAY;
}
 
/** enumerates circular patterns and creates restricted master problem */
SCIP_RETCODE SCIPprobdataSetupProblem(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_PROBDATA* probdata = SCIPgetProbData(scip);
   assert(probdata != NULL);
 
   /* collect parameters for verification */
   SCIP_CALL( SCIPgetRealParam(scip, "ringpacking/verification/nlptilimsoft", &probdata->nlptilimsoft) );
   SCIP_CALL( SCIPgetRealParam(scip, "ringpacking/verification/heurtilimsoft", &probdata->heurtilimsoft) );
   SCIP_CALL( SCIPgetLongintParam(scip, "ringpacking/verification/nlpnodelimsoft", &probdata->nlpnodelimsoft) );
   SCIP_CALL( SCIPgetIntParam(scip, "ringpacking/verification/heuriterlimsoft", &probdata->heuriterlimsoft) );
   SCIP_CALL( SCIPgetRealParam(scip, "ringpacking/verification/totaltilimsoft", &probdata->totaltilimsoft) );
 
   SCIP_CALL( setupProblem(scip, probdata) );
 
   return SCIP_OKAY;
}
 
/** enumerate all non-dominated circular patterns */
SCIP_RETCODE SCIPprobdataEnumeratePatterns(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_Real             nlptilim,           /**< time limit for each NLP verification */
   SCIP_Real             heurtilim,          /**< time limit for each call of the heuristics */
   SCIP_Real             totaltilim,         /**< total time limit for enumeration */
   SCIP_Longint          nlpnodelim,         /**< node limit for each NLP verification */
   int                   heuriterlim         /**< iteration limit for each call of the heuristics */
   )
{
   SCIP_PATTERN* pattern;
   int* ms;
   int* nselected;
   SCIP_Real timeleft;
   int ntypes;
   int t;
 
   assert(probdata != NULL);
   ntypes = SCIPprobdataGetNTypes(probdata);
   assert(ntypes > 0);
 
   /* create data that is used for the whole enumeration algorithm */
   SCIP_CALL( SCIPpatternCreateCircular(scip, &pattern, 0) );
   SCIP_CALL( SCIPallocBufferArray(scip, &ms, ntypes) );
   SCIP_CALL( SCIPallocBufferArray(scip, &nselected, ntypes) );
   BMSclearMemoryArray(nselected, ntypes);
   BMSclearMemoryArray(ms, ntypes);
 
   /* compute time limit */
   SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timeleft) );
   timeleft = MAX(0.0, MIN(timeleft - SCIPgetTotalTime(scip), totaltilim)); /*lint !e666*/
 
   /* find all circlular patterns of each type separately */
   for( t = 0; t < ntypes; ++t )
   {
      int k;
 
      for( k = t+1; k < ntypes; ++k )
         ms[k] = maxCircles(scip, probdata, t, k);
 
      SCIPpatternSetType(pattern, t);
      SCIP_CALL( enumeratePatterns(scip, probdata, pattern, ms, nselected, nlptilim, heurtilim, nlpnodelim,
         heuriterlim, &timeleft) );
   }
 
   /* release memory */
   SCIPfreeBufferArray(scip, &nselected);
   SCIPfreeBufferArray(scip, &ms);
   SCIPpatternRelease(scip, &pattern);
 
   /* filter circular patterns */
   SCIP_CALL( filterPatterns(scip, probdata) );
 
   return SCIP_OKAY;
}
 
/** returns number of different types */
int SCIPprobdataGetNTypes(
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   return probdata->ntypes;
}
 
/** returns all external radii */
SCIP_Real* SCIPprobdataGetRexts(
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   return probdata->rexts;
}
 
/** returns all internal radii */
SCIP_Real* SCIPprobdataGetRints(
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   return probdata->rints;
}
 
/** returns all demands */
int* SCIPprobdataGetDemands(
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   return probdata->demands;
}
 
/** returns the width of each rectangle */
SCIP_Real SCIPprobdataGetWidth(
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   return probdata->width;
}
 
 
/** returns the height of each rectangle */
SCIP_Real SCIPprobdataGetHeight(
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   return probdata->height;
}
 
/** returns all information about circular patterns */
void SCIPprobdataGetCInfos(
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PATTERN***       cpatterns,          /**< pointer to store the circular patterns (might be NULL) */
   SCIP_VAR***           cvars,              /**< pointer to store the variables corresponding circular patterns (might be NULL) */
   int*                  ncpatterns          /**< pointer to store the number of circular patterns (might be NULL) */
   )
{
   assert(probdata != NULL);
 
   if( cpatterns != NULL )
      *cpatterns = probdata->cpatterns;
   if( cvars != NULL )
      *cvars= probdata->cvars;
   if( ncpatterns != NULL )
      *ncpatterns = probdata->ncpatterns;
}
 
/** returns all information about rectangular patterns */
void SCIPprobdataGetRInfos(
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PATTERN***       rpatterns,          /**< pointer to store the rectangular patterns (might be NULL) */
   SCIP_VAR***           rvars,              /**< pointer to store the variables corresponding rectangular patterns (might be NULL) */
   int*                  nrpatterns          /**< pointer to store the number of rectangular patterns (might be NULL) */
   )
{
   assert(probdata != NULL);
 
   if( rpatterns != NULL )
      *rpatterns = probdata->rpatterns;
   if( rvars != NULL )
      *rvars= probdata->rvars;
   if( nrpatterns != NULL )
      *nrpatterns = probdata->nrpatterns;
}
 
/** returns array of set pattern constraints */
SCIP_CONS** SCIPprobdataGetPatternConss(
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   return probdata->patternconss;
}
 
/** adds given variable to the problem data */
SCIP_RETCODE SCIPprobdataAddVar(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PATTERN*         pattern,            /**< pattern */
   SCIP_VAR*             var                 /**< variables to add */
   )
{
   SCIP_PATTERN* copy;
 
   assert(probdata != NULL);
   assert(pattern != NULL);
   assert(SCIPpatternGetPackableStatus(pattern) != SCIP_PACKABLE_NO);
 
   /* copy pattern */
   SCIP_CALL( SCIPpatternCopy(scip, pattern, &copy) );
   SCIPcheckPattern(scip, probdata, copy);
 
   if( SCIPpatternGetPatternType(pattern) == SCIP_PATTERNTYPE_CIRCULAR )
   {
      SCIP_CALL( ensureSize(scip, probdata, SCIP_PATTERNTYPE_CIRCULAR, probdata->ncpatterns + 1) );
      probdata->cpatterns[probdata->ncpatterns] = copy;
      probdata->cvars[probdata->ncpatterns] = var;
      ++(probdata->ncpatterns);
   }
   else
   {
      SCIP_CALL( ensureSize(scip, probdata, SCIP_PATTERNTYPE_RECTANGULAR, probdata->nrpatterns + 1) );
      probdata->rpatterns[probdata->nrpatterns] = copy;
      probdata->rvars[probdata->nrpatterns] = var;
      ++(probdata->nrpatterns);
   }
 
   /* capture variable and pattern */
   if( var != NULL )
   {
      SCIP_CALL( SCIPcaptureVar(scip, var) );
   }
 
   return SCIP_OKAY;
}
 
/** updates the dual bound */
void SCIPprobdataUpdateDualbound(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_Real             dualbound           /**< new dual bound */
   )
{
   assert(probdata != NULL);
 
   if( !probdata->isdualinvalid && SCIPisFeasLT(scip, probdata->dualbound, dualbound) )
   {
      SCIPinfoMessage(scip, NULL, "+++++++++++++ update dual bound to %g\n", dualbound);
      probdata->dualbound = dualbound;
   }
}
 
/** marks that further reported dual bounds are not valid */
void SCIPprobdataInvalidateDualbound(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   if( !probdata->isdualinvalid )
   {
      SCIPinfoMessage(scip, NULL, "+++++++++++++ invalidate dual bound\n");
      probdata->isdualinvalid = TRUE;
   }
}
 
/** returns whether dual bound is marked to be invalid */
SCIP_Bool SCIPprobdataIsDualboundInvalid(
   SCIP_PROBDATA*        probdata            /**< problem data */
   )
{
   assert(probdata != NULL);
 
   return probdata->isdualinvalid;
}
 
/** Tries to pack a list of elements into a specified boundary circle by using a simple left-first bottom-second
 *  heuristic. Returns the number of elements that could be stored and indicated which ones these are in the buffer
 *  parameter ispacked. This auxiliary method can be used both to find such a packing or to verify a certain pattern.
 */
void SCIPpackCirclesGreedy(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_Real*            rexts,              /**< outer radii of elements (in original order of probdata) */
   SCIP_Real*            xs,                 /**< buffer to store the resulting x-coordinates */
   SCIP_Real*            ys,                 /**< buffer to store the resulting y-coordinates */
   SCIP_Real             rbounding,          /**< inner radius of bounding circle (ignored for rectangular patterns) */
   SCIP_Real             width,              /**< width of the rectangle */
   SCIP_Real             height,             /**< height of the rectangle */
   SCIP_Bool*            ispacked,           /**< buffer to store which elements could be packed */
   int*                  elements,           /**< the order of the elements in the pattern */
   int                   nelements,          /**< number of elements in the pattern */
   SCIP_PATTERNTYPE      patterntype,        /**< the pattern type (rectangular or circular) */
   int*                  npacked,            /**< pointer to store the number of packed elements */
   int                   ncalls              /**< total number of calls of the packing heuristic */
   )
{
   SCIP_Real rmax;
   SCIP_Bool added;
   int i;
 
   assert(rexts != NULL);
   assert(xs != NULL);
   assert(ys != NULL);
   assert(ispacked != NULL);
   assert(elements != NULL);
   assert(nelements > 0);
   assert(npacked != NULL);
 
   /* no element packed so far */
   BMSclearMemoryArray(ispacked, nelements);
 
   /* place first element at left-most position */
   if( patterntype == SCIP_PATTERNTYPE_CIRCULAR )
   {
      assert(rexts[elements[0]] <= rbounding);
      xs[0] = rexts[elements[0]] - rbounding;
      ys[0] = 0.0;
   }
   else
   {
      assert(2.0 * rexts[elements[0]] <= width);
      assert(2.0 * rexts[elements[0]] <= height);
      xs[0] = rexts[elements[0]];
      ys[0] = rexts[elements[0]];
   }
 
   /* initialize results */
   (*npacked) = 1;
   ispacked[0] = TRUE;
   added = TRUE;
 
   /* find max radius */
   rmax = rexts[elements[0]];
   for( i = 1; i < nelements; ++i )
   {
      if( rexts[elements[i]] > rmax )
         rmax = rexts[elements[i]];
   }
 
   /* iterate over all elements and try to pack them */
   while( added )
   {
      added = FALSE;
 
      for( i = 1; i < nelements; ++i )
      {
         SCIP_Real bestx = SCIP_INVALID;
         SCIP_Real besty = SCIP_INVALID;
 
         /* skip packed elements */
         if( ispacked[i] )
            continue;
 
         /* use trivial candidates */
         computePosTrivial(scip, elements, nelements, rexts, xs, ys, i, ispacked, rmax, rbounding, width, height,
            patterntype, &bestx, &besty, ncalls);
 
         /* consider circles intersection a previous circle and the boundary ring */
         if( patterntype == SCIP_PATTERNTYPE_CIRCULAR )
            computePosRingCircle(scip, elements, nelements, rexts, xs, ys, i, ispacked, rmax, rbounding, &bestx,
               &besty, ncalls);
         else
            computePosRectangleCircle(scip, elements, nelements, rexts, xs, ys, i, ispacked, rmax, width, height, &bestx,
               &besty, ncalls);
 
         /* consider circles that have been packed already */
         computePosCircleCircle(scip, elements, nelements, rexts, xs, ys, i, ispacked, rmax, rbounding, width, height,
            patterntype, &bestx, &besty, ncalls);
 
         /* pack circle if a possible position has been found */
         if( bestx != SCIP_INVALID && besty != SCIP_INVALID ) /*lint !e777*/
         {
            assert(!ispacked[i]);
            ispacked[i] = TRUE;
            xs[i] = bestx;
            ys[i] = besty;
            ++(*npacked);
            added = TRUE;
         }
      }
   }
 
   return;
}
 
/** verifies a circular pattern heuristically */
SCIP_RETCODE SCIPverifyCircularPatternHeuristic(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PATTERN*         pattern,            /**< pattern */
   SCIP_Real             timelim,            /**< time limit */
   int                   iterlim             /**< iteration limit */
   )
{
   SCIP_Real* rexts;
   SCIP_Real* rints;
   SCIP_Real* scores;
   SCIP_Real* xs;
   SCIP_Real* ys;
   SCIP_Bool* ispacked;
   int* elements;
   int* pos;
   SCIP_Real timestart;
   int nelements;
   int niters;
   int type;
   int i;
 
   assert(probdata != NULL);
   assert(pattern != NULL);
   assert(iterlim > 0);
   assert(SCIPpatternGetPatternType(pattern) == SCIP_PATTERNTYPE_CIRCULAR);
   assert(SCIPpatternGetPackableStatus(pattern) == SCIP_PACKABLE_UNKNOWN);
   assert(SCIPpatternGetCircleType(pattern) < SCIPprobdataGetNTypes(probdata));
 
   /* check whether there is any time left */
   if( timelim <= 0.0 )
      return SCIP_OKAY;
 
   rexts = SCIPprobdataGetRexts(probdata);
   rints = SCIPprobdataGetRints(probdata);
   nelements = SCIPpatternGetNElemens(pattern);
   type = SCIPpatternGetCircleType(pattern);
   assert(type >= 0 && type < SCIPprobdataGetNTypes(probdata));
 
   /* pattern is empty -> set status to packable */
   if( SCIPpatternGetNElemens(pattern) == 0 )
   {
      SCIPpatternSetPackableStatus(pattern, SCIP_PACKABLE_YES);
      SCIPcheckPattern(scip, probdata, pattern);
      return SCIP_OKAY;
   }
 
   /* pattern contains only one element -> compare radii */
   if( SCIPpatternGetNElemens(pattern) == 1 )
   {
      int elemtype;
 
      elemtype = SCIPpatternGetElementType(pattern, 0);
      assert(elemtype >= 0 && elemtype < SCIPprobdataGetNTypes(probdata));
 
      /* check whether element fits into the circular pattern */
      if( SCIPisGE(scip, rints[type], rexts[elemtype]) )
      {
         SCIPpatternSetElementPos(pattern, 0, rexts[elemtype]-rints[type], 0.0);
         SCIPpatternSetPackableStatus(pattern, SCIP_PACKABLE_YES);
      }
      else
         SCIPpatternSetPackableStatus(pattern, SCIP_PACKABLE_NO);
 
      SCIPcheckPattern(scip, probdata, pattern);
      return SCIP_OKAY;
   }
 
   timestart = SCIPgetTotalTime(scip);
   niters = 0;
 
   /* store elements in a separate array; remember positions of elements in the pattern */
   SCIP_CALL( SCIPallocBufferArray(scip, &pos, nelements) );
   SCIP_CALL( SCIPallocBufferArray(scip, &scores, nelements) );
   SCIP_CALL( SCIPallocBufferArray(scip, &ispacked, nelements) );
   SCIP_CALL( SCIPallocBufferArray(scip, &elements, nelements) );
   SCIP_CALL( SCIPallocBufferArray(scip, &xs, nelements) );
   SCIP_CALL( SCIPallocBufferArray(scip, &ys, nelements) );
   for( i = 0; i < nelements; ++i )
   {
      elements[i] = SCIPpatternGetElementType(pattern, i);
      ispacked[i] = FALSE;
      pos[i] = i;
   }
 
   /* main loop for calling heuristic verification */
   while( SCIPpatternGetPackableStatus(pattern) == SCIP_PACKABLE_UNKNOWN
      && niters < iterlim
      && SCIPgetTotalTime(scip) - timestart <= timelim )
   {
      int npacked;
 
      /* compute scores depending on iteration counter */
      computeScores(scip, probdata, elements, nelements, scores, niters);
 
      /* sort elements in non-increasing order */
      SCIPsortDownRealIntInt(scores, elements, pos, nelements);
 
      /* call heuristic */
      SCIPpackCirclesGreedy(scip, rexts, xs, ys, rints[type], SCIPprobdataGetWidth(probdata),
         SCIPprobdataGetHeight(probdata), ispacked, elements, nelements, SCIP_PATTERNTYPE_CIRCULAR, &npacked, niters);
 
      /* check whether all elements could have been packed */
      if( npacked == nelements )
      {
         for( i = 0; i < nelements; ++i )
         {
            assert(elements[i] == SCIPpatternGetElementType(pattern, pos[i]));
            SCIPpatternSetElementPos(pattern, pos[i], xs[i], ys[i]);
         }
         SCIPpatternSetPackableStatus(pattern, SCIP_PACKABLE_YES);
 
         SCIPdebugMsg(scip, "heuristic verified pattern after %d iterations\n", niters + 1);
      }
 
      ++niters;
   }
 
   SCIPcheckPattern(scip, probdata, pattern);
 
   /* free memory */
   SCIPfreeBufferArray(scip, &ys);
   SCIPfreeBufferArray(scip, &xs);
   SCIPfreeBufferArray(scip, &elements);
   SCIPfreeBufferArray(scip, &ispacked);
   SCIPfreeBufferArray(scip, &scores);
   SCIPfreeBufferArray(scip, &pos);
 
   return SCIP_OKAY;
}
 
/** verifies a circular pattern via a verification NLP */
SCIP_RETCODE SCIPverifyCircularPatternNLP(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PATTERN*         pattern,            /**< pattern */
   SCIP_Real             timelim,            /**< time limit */
   SCIP_Longint          nodelim             /**< node limit */
   )
{
   SCIP* subscip;
   SCIP_CONS* cons;
   SCIP_VAR** xvars;
   SCIP_VAR** yvars;
   SCIP_VAR* quadvars1[6];
   SCIP_VAR* quadvars2[6];
   SCIP_Real quadcoefs[6];
   SCIP_Real* rexts;
   SCIP_Real* rints;
   char name[SCIP_MAXSTRLEN];
   int nelems;
   int type;
   int k;
 
   assert(probdata != NULL);
   assert(pattern != NULL);
   assert(SCIPpatternGetPatternType(pattern) == SCIP_PATTERNTYPE_CIRCULAR);
   assert(SCIPpatternGetPackableStatus(pattern) == SCIP_PACKABLE_UNKNOWN);
 
   /* check whether there is any time left */
   if( timelim <= 0.0 )
      return SCIP_OKAY;
 
   rexts = SCIPprobdataGetRexts(probdata);
   rints = SCIPprobdataGetRints(probdata);
   type = SCIPpatternGetCircleType(pattern);
   nelems = SCIPpatternGetNElemens(pattern);
 
   /* set up the sub-SCIP */
   SCIP_CALL( SCIPcreate(&subscip) );
   SCIP_CALL( SCIPcreateProbBasic(subscip, "verify") );
   SCIP_CALL( SCIPincludeDefaultPlugins(subscip) );
 
   /* allocate memory for (x,y) variables */
   SCIP_CALL( SCIPallocBufferArray(scip, &xvars, nelems) );
   SCIP_CALL( SCIPallocBufferArray(scip, &yvars, nelems) );
 
   /* set feasibility emphasis settings */
   SCIP_CALL( SCIPsetEmphasis(subscip, SCIP_PARAMEMPHASIS_FEASIBILITY, TRUE) );
 
   /* set working limit */
   SCIP_CALL( SCIPsetIntParam(subscip, "limits/solutions", 1) );
   SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", timelim) );
   SCIP_CALL( SCIPsetLongintParam(subscip, "limits/nodes", nodelim) );
 
#ifndef SCIP_DEBUG
   SCIPsetMessagehdlrQuiet(subscip, TRUE);
#endif
 
   /* create (x,y) variables */
   for( k = 0; k < nelems; ++k )
   {
      int elemtype;
 
      elemtype = SCIPpatternGetElementType(pattern, k);
      assert(elemtype >= 0 && elemtype < SCIPprobdataGetNTypes(probdata));
 
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "x_%d", k);
      SCIP_CALL( SCIPcreateVarBasic(subscip, &xvars[k], name, rexts[elemtype] - rints[type],
            rints[type] - rexts[elemtype], 0.0, SCIP_VARTYPE_CONTINUOUS) );
      SCIP_CALL( SCIPaddVar(subscip, xvars[k]) );
 
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "y_%d", k);
      SCIP_CALL( SCIPcreateVarBasic(subscip, &yvars[k], name, rexts[elemtype] - rints[type],
            rints[type] - rexts[elemtype], 1.0, SCIP_VARTYPE_CONTINUOUS) );
      SCIP_CALL( SCIPaddVar(subscip, yvars[k]) );
   }
 
   /* create non-overlapping constraints */
   for( k = 0; k < nelems; ++k )
   {
      int elemtype1;
      int l;
 
      elemtype1 = SCIPpatternGetElementType(pattern, k);
      assert(elemtype1 >= 0 && elemtype1 < SCIPprobdataGetNTypes(probdata));
 
      for( l = k + 1; l < nelems; ++l )
      {
         int elemtype2;
 
         elemtype2 = SCIPpatternGetElementType(pattern, l);
         assert(elemtype2 >= 0 && elemtype2 < SCIPprobdataGetNTypes(probdata));
 
         quadvars1[0] = xvars[k]; quadvars2[0] = xvars[k]; quadcoefs[0] =  1.0;
         quadvars1[1] = xvars[k]; quadvars2[1] = xvars[l]; quadcoefs[1] = -2.0;
         quadvars1[2] = xvars[l]; quadvars2[2] = xvars[l]; quadcoefs[2] =  1.0;
         quadvars1[3] = yvars[k]; quadvars2[3] = yvars[k]; quadcoefs[3] =  1.0;
         quadvars1[4] = yvars[k]; quadvars2[4] = yvars[l]; quadcoefs[4] = -2.0;
         quadvars1[5] = yvars[l]; quadvars2[5] = yvars[l]; quadcoefs[5] =  1.0;
 
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "over_%d_%d", k, l);
         SCIP_CALL( SCIPcreateConsQuadraticNonlinear(subscip, &cons, name, 0, NULL, NULL, 6, quadvars1, quadvars2,
               quadcoefs, SQR(rexts[elemtype1] + rexts[elemtype2]), SCIPinfinity(subscip),
               TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
 
         SCIP_CALL( SCIPaddCons(subscip, cons) );
         SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
      }
   }
 
   /* create non-overlapping constraints with outer ring */
   for( k = 0; k < nelems; ++k )
   {
      int elemtype;
 
      elemtype = SCIPpatternGetElementType(pattern, k);
      assert(elemtype >= 0 && elemtype < SCIPprobdataGetNTypes(probdata));
 
      quadvars1[0] = xvars[k]; quadvars2[0] = xvars[k]; quadcoefs[0] = 1.0;
      quadvars1[1] = yvars[k]; quadvars2[1] = yvars[k]; quadcoefs[1] = 1.0;
 
      (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "bound_%d", k);
      SCIP_CALL( SCIPcreateConsQuadraticNonlinear(subscip, &cons, name, 0, NULL, NULL, 2, quadvars1, quadvars2, quadcoefs,
            0.0, SQR(rints[type] - rexts[elemtype]),
            TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE) );
 
      SCIP_CALL( SCIPaddCons(subscip, cons) );
      SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
   }
 
   /* sort circles in x direction if they have the same type */
   for( k = 0; k < nelems - 1; ++k )
   {
      int elemtype1;
      int l;
 
      elemtype1 = SCIPpatternGetElementType(pattern, k);
      assert(elemtype1 >= 0 && elemtype1 < SCIPprobdataGetNTypes(probdata));
 
      for( l = k + 1; l < nelems; ++l )
      {
         int elemtype2;
 
         elemtype2 = SCIPpatternGetElementType(pattern, k+1);
         assert(elemtype2 >= 0 && elemtype2 < SCIPprobdataGetNTypes(probdata));
 
         if( elemtype1 != elemtype2 )
            continue;
 
         (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "sortx_%d_%d", k, l);
         SCIP_CALL( SCIPcreateConsBasicLinear(subscip, &cons, name, 0, NULL, NULL, -SCIPinfinity(subscip), 0.0) );
         SCIP_CALL( SCIPaddCoefLinear(subscip, cons, xvars[k], 1.0) );
         SCIP_CALL( SCIPaddCoefLinear(subscip, cons, xvars[l], -1.0) );
 
         SCIP_CALL( SCIPaddCons(subscip, cons) );
         SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
      }
   }
 
   /* solve verification NLP */
   SCIPdebugMsg(scip, "--------------------- SOLVE VERIFICATION NLP -------------------\n");
   SCIP_CALL( SCIPsolve(subscip) );
   SCIPdebugMsg(scip, "----------------------------------------------------------------\n");
 
   SCIPdebugMsg(scip, "result of verification NLP: nsols=%d solstat=%d\n",
      SCIPgetNSols(subscip), SCIPgetStatus(subscip));
 
   /* check whether a solution could be found or whether the problem is proven to be infeasible */
   if( SCIPgetNSols(subscip) > 0 )
   {
      SCIPpatternSetPackableStatus(pattern, SCIP_PACKABLE_YES);
 
      for( k = 0; k < nelems; ++k )
      {
         SCIP_Real solx = SCIPgetSolVal(subscip, SCIPgetBestSol(subscip), xvars[k]);
         SCIP_Real soly = SCIPgetSolVal(subscip, SCIPgetBestSol(subscip), yvars[k]);
 
         SCIPpatternSetElementPos(pattern, k, solx, soly);
      }
 
      SCIPcheckPattern(scip, probdata, pattern);
   }
   else if( SCIPgetStatus(subscip) == SCIP_STATUS_INFEASIBLE )
      SCIPpatternSetPackableStatus(pattern, SCIP_PACKABLE_NO);
 
   /* free all variables */
   for( k = 0; k < nelems; ++k )
   {
      SCIP_CALL( SCIPreleaseVar(subscip, &yvars[k]) );
      SCIP_CALL( SCIPreleaseVar(subscip, &xvars[k]) );
   }
 
   /* free memory */
   SCIPfreeBufferArray(scip, &yvars);
   SCIPfreeBufferArray(scip, &xvars);
   SCIP_CALL( SCIPfree(&subscip) );
 
   return SCIP_OKAY;
}
 
/** check a pattern for consistency */
void SCIPcheckPattern(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_PROBDATA*        probdata,           /**< problem data */
   SCIP_PATTERN*         pattern             /**< pattern */
   )
{  /*lint --e{715}*/
#ifndef NDEBUG
   SCIP_Real* rexts;
   SCIP_Real* rints;
   SCIP_Real width;
   SCIP_Real height;
   int i;
 
   assert(probdata != NULL);
   assert(pattern != NULL);
 
   rexts = SCIPprobdataGetRexts(probdata);
   rints = SCIPprobdataGetRints(probdata);
   width = SCIPprobdataGetWidth(probdata);
   height = SCIPprobdataGetHeight(probdata);
 
   /* check types */
   for( i = 0; i < SCIPpatternGetNElemens(pattern); ++i )
   {
      int type = SCIPpatternGetElementType(pattern, i);
 
      assert(type >= 0);
      assert(type < SCIPprobdataGetNTypes(probdata));
   }
 
   /* check positions iff packable */
   if( SCIPpatternGetPackableStatus(pattern) != SCIP_PACKABLE_YES )
      return;
 
   for( i = 0; i < SCIPpatternGetNElemens(pattern); ++i )
   {
      SCIP_Real xi = SCIPpatternGetElementPosX(pattern, i);
      SCIP_Real yi = SCIPpatternGetElementPosY(pattern, i);
      int typei = SCIPpatternGetElementType(pattern, i);
      int j;
 
      /* check distance between circles */
      for( j = i + 1; j < SCIPpatternGetNElemens(pattern); ++j )
      {
         SCIP_Real xj = SCIPpatternGetElementPosX(pattern, j);
         SCIP_Real yj = SCIPpatternGetElementPosY(pattern, j);
         int typej = SCIPpatternGetElementType(pattern, j);
 
         assert(SCIPisFeasGE(scip, sqrt(SQR(xi - xj) + SQR(yi - yj)), rexts[typei] + rexts[typej]));
      }
 
      /* check distance to boundary */
      if( SCIPpatternGetPatternType(pattern) == SCIP_PATTERNTYPE_CIRCULAR )
      {
         SCIP_Real distance = sqrt(SQR(xi) + SQR(yi));
         int patterntype = SCIPpatternGetCircleType(pattern);
 
         assert(patterntype >= 0);
         assert(patterntype < SCIPprobdataGetNTypes(probdata));
         assert(SCIPisFeasLE(scip, distance, rints[patterntype] - rexts[typei]));
      }
      else
      {
         assert(SCIPisFeasGE(scip, xi, rexts[typei]));
         assert(SCIPisFeasLE(scip, xi, width - rexts[typei]));
         assert(SCIPisFeasGE(scip, yi, rexts[typei]));
         assert(SCIPisFeasLE(scip, yi, height - rexts[typei]));
      }
   }
#endif
}
 
/**@} */