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   heur_multistart.c
26   	 * @ingroup DEFPLUGINS_HEUR
27   	 * @brief  multistart heuristic for convex and nonconvex MINLPs
28   	 * @author Benjamin Mueller
29   	 */
30   	
31   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32   	
33   	#include "blockmemshell/memory.h"
34   	#include "scip/scip_expr.h"
35   	#include "scip/pub_expr.h"
36   	#include "scip/heur_multistart.h"
37   	#include "scip/heur_subnlp.h"
38   	#include "scip/pub_heur.h"
39   	#include "scip/pub_message.h"
40   	#include "scip/pub_misc.h"
41   	#include "scip/pub_misc_sort.h"
42   	#include "scip/pub_nlp.h"
43   	#include "scip/pub_var.h"
44   	#include "scip/scip_general.h"
45   	#include "scip/scip_heur.h"
46   	#include "scip/scip_mem.h"
47   	#include "scip/scip_message.h"
48   	#include "scip/scip_nlp.h"
49   	#include "scip/scip_nlpi.h"
50   	#include "scip/scip_numerics.h"
51   	#include "scip/scip_param.h"
52   	#include "scip/scip_prob.h"
53   	#include "scip/scip_randnumgen.h"
54   	#include "scip/scip_sol.h"
55   	#include "scip/scip_timing.h"
56   	#include <string.h>
57   	
58   	
59   	#define HEUR_NAME             "multistart"
60   	#define HEUR_DESC             "multistart heuristic for convex and nonconvex MINLPs"
61   	#define HEUR_DISPCHAR         SCIP_HEURDISPCHAR_LNS
62   	#define HEUR_PRIORITY         -2100000
63   	#define HEUR_FREQ             0
64   	#define HEUR_FREQOFS          0
65   	#define HEUR_MAXDEPTH         -1
66   	#define HEUR_TIMING           SCIP_HEURTIMING_AFTERNODE
67   	#define HEUR_USESSUBSCIP      TRUE           /**< does the heuristic use a secondary SCIP instance? */
68   	
69   	#define DEFAULT_RANDSEED      131            /**< initial random seed */
70   	#define DEFAULT_NRNDPOINTS    100            /**< default number of generated random points per call */
71   	#define DEFAULT_MAXBOUNDSIZE  2e+4           /**< default maximum variable domain size for unbounded variables */
72   	#define DEFAULT_MAXITER       300            /**< default number of iterations to reduce the violation of a point */
73   	#define DEFAULT_MINIMPRFAC    0.05           /**< default minimum required improving factor to proceed in improvement of a point */
74   	#define DEFAULT_MINIMPRITER   10             /**< default number of iteration when checking the minimum improvement */
75   	#define DEFAULT_MAXRELDIST    0.15           /**< default maximum distance between two points in the same cluster */
76   	#define DEFAULT_GRADLIMIT     5e+6           /**< default limit for gradient computations for all improvePoint() calls */
77   	#define DEFAULT_MAXNCLUSTER   3              /**< default maximum number of considered clusters per heuristic call */
78   	#define DEFAULT_ONLYNLPS      TRUE           /**< should the heuristic run only on continuous problems? */
79   	
80   	#define MINFEAS               -1e+4          /**< minimum feasibility for a point; used for filtering and improving
81   	                                              *   feasibility */
82   	#define MINIMPRFAC            0.95           /**< improvement factor used to discard randomly generated points with a
83   	                                              *   too large objective value */
84   	#define GRADCOSTFAC_LINEAR    1.0            /**< gradient cost factor for the number of linear variables */
85   	#define GRADCOSTFAC_NONLINEAR 3.0            /**< gradient cost factor for the number of nodes in nonlinear expression */
86   	
87   	/*
88   	 * Data structures
89   	 */
90   	
91   	/** primal heuristic data */
92   	struct SCIP_HeurData
93   	{
94   	   int                   nrndpoints;         /**< number of random points generated per execution call */
95   	   SCIP_Real             maxboundsize;       /**< maximum variable domain size for unbounded variables */
96   	   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generator */
97   	   SCIP_HEUR*            heursubnlp;         /**< sub-NLP heuristic */
98   	
99   	   int                   maxiter;            /**< number of iterations to reduce the maximum violation of a point */
100  	   SCIP_Real             minimprfac;         /**< minimum required improving factor to proceed in the improvement of a single point */
101  	   int                   minimpriter;        /**< number of iteration when checking the minimum improvement */
102  	
103  	   SCIP_Real             maxreldist;         /**< maximum distance between two points in the same cluster */
104  	   SCIP_Real             gradlimit;          /**< limit for gradient computations for all improvePoint() calls (0 for no limit) */
105  	   int                   maxncluster;        /**< maximum number of considered clusters per heuristic call */
106  	   SCIP_Bool             onlynlps;           /**< should the heuristic run only on continuous problems? */
107  	};
108  	
109  	
110  	/*
111  	 * Local methods
112  	 */
113  	
114  	
115  	/** returns an unique index of a variable in the range of 0,..,SCIPgetNVars(scip)-1 */
116  	#ifndef NDEBUG
117  	static
118  	int getVarIndex(
119  	   SCIP_HASHMAP*         varindex,           /**< maps variables to indicies between 0,..,SCIPgetNVars(scip)-1 */
120  	   SCIP_VAR*             var                 /**< variable */
121  	   )
122  	{
123  	   assert(varindex != NULL);
124  	   assert(var != NULL);
125  	   assert(SCIPhashmapExists(varindex, (void*)var));
126  	
127  	   return SCIPhashmapGetImageInt(varindex, (void*)var);
128  	}
129  	#else
130  	#define getVarIndex(varindex,var) (SCIPhashmapGetImageInt((varindex), (void*)(var)))
131  	#endif
132  	
133  	/** samples and stores random points; stores points which have a better objective value than the current incumbent
134  	 *  solution
135  	 */
136  	static
137  	SCIP_RETCODE sampleRandomPoints(
138  	   SCIP*                 scip,               /**< SCIP data structure */
139  	   SCIP_SOL**            rndpoints,          /**< array to store all random points */
140  	   int                   nmaxrndpoints,      /**< maximum number of random points to compute */
141  	   SCIP_Real             maxboundsize,       /**< maximum variable domain size for unbounded variables */
142  	   SCIP_RANDNUMGEN*      randnumgen,         /**< random number generator */
143  	   SCIP_Real             bestobj,            /**< objective value in the transformed space of the current incumbent */
144  	   int*                  nstored             /**< pointer to store the number of randomly generated points */
145  	   )
146  	{
147  	   SCIP_VAR** vars;
148  	   SCIP_SOL* sol;
149  	   SCIP_Real val;
150  	   SCIP_Real lb;
151  	   SCIP_Real ub;
152  	   int nvars;
153  	   int niter;
154  	   int i;
155  	
156  	   assert(scip != NULL);
157  	   assert(rndpoints != NULL);
158  	   assert(nmaxrndpoints > 0);
159  	   assert(maxboundsize > 0.0);
160  	   assert(randnumgen != NULL);
161  	   assert(nstored != NULL);
162  	
163  	   vars = SCIPgetVars(scip);
164  	   nvars = SCIPgetNVars(scip);
165  	   *nstored = 0;
166  	
167  	   SCIP_CALL( SCIPcreateSol(scip, &sol, NULL) );
168  	
169  	   for( niter = 0; niter < 3 * nmaxrndpoints && *nstored < nmaxrndpoints; ++niter )
170  	   {
171  	      /* reset solution, in case the old one had infinite objective, which can give difficulties in updating the obj value */
172  	      SCIP_CALL( SCIPclearSol(scip, sol) );
173  	
174  	      for( i = 0; i < nvars; ++i )
175  	      {
176  	         lb = MIN(SCIPvarGetLbLocal(vars[i]), SCIPvarGetUbLocal(vars[i])); /*lint !e666*/
177  	         ub = MAX(SCIPvarGetLbLocal(vars[i]), SCIPvarGetUbLocal(vars[i])); /*lint !e666*/
178  	
179  	         if( SCIPisFeasEQ(scip, lb, ub) )
180  	            val = (lb + ub) / 2.0;
181  	         /* use a smaller domain for unbounded variables */
182  	         else if( !SCIPisInfinity(scip, -lb) && !SCIPisInfinity(scip, ub) )
183  	            val = SCIPrandomGetReal(randnumgen, lb, ub);
184  	         else if( !SCIPisInfinity(scip, -lb) )
185  	            val = lb + SCIPrandomGetReal(randnumgen, 0.0, maxboundsize);
186  	         else if( !SCIPisInfinity(scip, ub) )
187  	            val = ub - SCIPrandomGetReal(randnumgen, 0.0, maxboundsize);
188  	         else
189  	         {
190  	            assert(SCIPisInfinity(scip, -lb) && SCIPisInfinity(scip, ub));
191  	            val = SCIPrandomGetReal(randnumgen, -0.5*maxboundsize, 0.5*maxboundsize);
192  	         }
193  	         assert(SCIPisFeasGE(scip, val, lb) && SCIPisFeasLE(scip, val, ub));
194  	
195  	         /* set solution value; round the sampled point for integer variables */
196  	         if( SCIPvarGetType(vars[i]) < SCIP_VARTYPE_CONTINUOUS )
197  	            val = SCIPfeasRound(scip, val);
198  	         SCIP_CALL( SCIPsetSolVal(scip, sol, vars[i], val) );
199  	      }
200  	
201  	      /* add solution if it is good enough */
202  	      if( SCIPisLE(scip, SCIPgetSolTransObj(scip, sol), bestobj) )
203  	      {
204  	         SCIP_CALL( SCIPcreateSolCopy(scip, &rndpoints[*nstored], sol) );
205  	         ++(*nstored);
206  	      }
207  	   }
208  	   assert(*nstored <= nmaxrndpoints);
209  	   SCIPdebugMsg(scip, "found %d randomly generated points\n", *nstored);
210  	
211  	   SCIP_CALL( SCIPfreeSol(scip, &sol) );
212  	
213  	   return SCIP_OKAY;
214  	}
215  	
216  	/** computes the minimum feasibility of a given point; a negative value means that there is an infeasibility */
217  	static
218  	SCIP_RETCODE getMinFeas(
219  	   SCIP*                 scip,               /**< SCIP data structure */
220  	   SCIP_NLROW**          nlrows,             /**< array containing all nlrows */
221  	   int                   nnlrows,            /**< total number of nlrows */
222  	   SCIP_SOL*             sol,                /**< solution */
223  	   SCIP_Real*            minfeas             /**< buffer to store the minimum feasibility */
224  	   )
225  	{
226  	   SCIP_Real tmp;
227  	   int i;
228  	
229  	   assert(scip != NULL);
230  	   assert(sol != NULL);
231  	   assert(minfeas != NULL);
232  	   assert(nlrows != NULL);
233  	   assert(nnlrows > 0);
234  	
235  	   *minfeas = SCIPinfinity(scip);
236  	
237  	   for( i = 0; i < nnlrows; ++i )
238  	   {
239  	      assert(nlrows[i] != NULL);
240  	
241  	      SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrows[i], sol, &tmp) );
242  	      *minfeas = MIN(*minfeas, tmp);
243  	   }
244  	
245  	   return SCIP_OKAY;
246  	}
247  	
248  	/** computes the gradient for a given point and nonlinear row */
249  	static
250  	SCIP_RETCODE computeGradient(
251  	   SCIP*                 scip,               /**< SCIP data structure */
252  	   SCIP_NLROW*           nlrow,              /**< nonlinear row */
253  	   SCIP_SOL*             sol,                /**< solution to compute the gradient for */
254  	   SCIP_HASHMAP*         varindex,           /**< maps variables to indicies between 0,..,SCIPgetNVars(scip)-1 uniquely */
255  	   SCIP_EXPRITER*        exprit,             /**< expression iterator that can be used */
256  	   SCIP_Real*            grad,               /**< buffer to store the gradient; grad[varindex(i)] corresponds to SCIPgetVars(scip)[i] */
257  	   SCIP_Real*            norm                /**< buffer to store ||grad||^2  */
258  	   )
259  	{
260  	   SCIP_EXPR* expr;
261  	   SCIP_VAR* var;
262  	   int i;
263  	
264  	   assert(scip != NULL);
265  	   assert(nlrow != NULL);
266  	   assert(varindex != NULL);
267  	   assert(sol != NULL);
268  	   assert(norm != NULL);
269  	
270  	   BMSclearMemoryArray(grad, SCIPgetNVars(scip));
271  	   *norm = 0.0;
272  	
273  	   /* linear part */
274  	   for( i = 0; i < SCIPnlrowGetNLinearVars(nlrow); i++ )
275  	   {
276  	      var = SCIPnlrowGetLinearVars(nlrow)[i];
277  	      assert(var != NULL);
278  	      assert(getVarIndex(varindex, var) >= 0 && getVarIndex(varindex, var) < SCIPgetNVars(scip));
279  	
280  	      grad[getVarIndex(varindex, var)] += SCIPnlrowGetLinearCoefs(nlrow)[i];
281  	   }
282  	
283  	   /* expression part */
284  	   expr = SCIPnlrowGetExpr(nlrow);
285  	
286  	   if( expr != NULL )
287  	   {
288  	      assert(exprit != NULL);
289  	
290  	      SCIP_CALL( SCIPevalExprGradient(scip, expr, sol, 0L) );
291  	
292  	      /* TODO: change this when nlrows store the vars */
293  	      SCIP_CALL( SCIPexpriterInit(exprit, expr, SCIP_EXPRITER_DFS, FALSE) );
294  	      for( ; !SCIPexpriterIsEnd(exprit); expr = SCIPexpriterGetNext(exprit) )  /*lint !e441*/ /*lint !e440*/
295  	      {
296  	         if( !SCIPisExprVar(scip, expr) )
297  	            continue;
298  	
299  	         var = SCIPgetVarExprVar(expr);
300  	         assert(var != NULL);
301  	         assert(getVarIndex(varindex, var) >= 0 && getVarIndex(varindex, var) < SCIPgetNVars(scip));
302  	
303  	         grad[getVarIndex(varindex, var)] += SCIPexprGetDerivative(expr);
304  	      }
305  	   }
306  	
307  	   /* compute ||grad||^2 */
308  	   for( i = 0; i < SCIPgetNVars(scip); ++i )
309  	      *norm += SQR(grad[i]);
310  	
311  	   return SCIP_OKAY;
312  	}
313  	
314  	/** use consensus vectors to improve feasibility for a given starting point */
315  	static
316  	SCIP_RETCODE improvePoint(
317  	   SCIP*                 scip,               /**< SCIP data structure */
318  	   SCIP_NLROW**          nlrows,             /**< array containing all nlrows */
319  	   int                   nnlrows,            /**< total number of nlrows */
320  	   SCIP_HASHMAP*         varindex,           /**< maps variables to indicies between 0,..,SCIPgetNVars(scip)-1 */
321  	   SCIP_SOL*             point,              /**< random generated point */
322  	   int                   maxiter,            /**< maximum number of iterations */
323  	   SCIP_Real             minimprfac,         /**< minimum required improving factor to proceed in the improvement of a single point */
324  	   int                   minimpriter,        /**< number of iteration when checking the minimum improvement */
325  	   SCIP_Real*            minfeas,            /**< pointer to store the minimum feasibility */
326  	   SCIP_Real*            nlrowgradcosts,     /**< estimated costs for each gradient computation */
327  	   SCIP_Real*            gradcosts           /**< pointer to store the estimated gradient costs */
328  	   )
329  	{
330  	   SCIP_VAR** vars;
331  	   SCIP_EXPRITER* exprit;
332  	   SCIP_Real* grad;
333  	   SCIP_Real* updatevec;
334  	   SCIP_Real lastminfeas;
335  	   int nvars;
336  	   int r;
337  	   int i;
338  	
339  	   assert(varindex != NULL);
340  	   assert(point != NULL);
341  	   assert(maxiter > 0);
342  	   assert(minfeas != NULL);
343  	   assert(nlrows != NULL);
344  	   assert(nnlrows > 0);
345  	   assert(nlrowgradcosts != NULL);
346  	   assert(gradcosts != NULL);
347  	
348  	   *gradcosts = 0.0;
349  	
350  	   SCIP_CALL( getMinFeas(scip, nlrows, nnlrows, point, minfeas) );
351  	#ifdef SCIP_DEBUG_IMPROVEPOINT
352  	   printf("start minfeas = %e\n", *minfeas);
353  	#endif
354  	
355  	   /* stop since start point is feasible */
356  	   if( !SCIPisFeasLT(scip, *minfeas, 0.0) )
357  	   {
358  	#ifdef SCIP_DEBUG_IMPROVEPOINT
359  	      printf("start point is feasible");
360  	#endif
361  	      return SCIP_OKAY;
362  	   }
363  	
364  	   lastminfeas = *minfeas;
365  	   vars = SCIPgetVars(scip);
366  	   nvars = SCIPgetNVars(scip);
367  	
368  	   SCIP_CALL( SCIPallocBufferArray(scip, &grad, nvars) );
369  	   SCIP_CALL( SCIPallocBufferArray(scip, &updatevec, nvars) );
370  	   SCIP_CALL( SCIPcreateExpriter(scip, &exprit) );
371  	
372  	   /* main loop */
373  	   for( r = 0; r < maxiter && SCIPisFeasLT(scip, *minfeas, 0.0); ++r )
374  	   {
375  	      SCIP_Real feasibility;
376  	      SCIP_Real activity;
377  	      SCIP_Real nlrownorm;
378  	      SCIP_Real scale;
379  	      int nviolnlrows;
380  	
381  	      BMSclearMemoryArray(updatevec, nvars);
382  	      nviolnlrows = 0;
383  	
384  	      for( i = 0; i < nnlrows; ++i )
385  	      {
386  	         int j;
387  	
388  	         SCIP_CALL( SCIPgetNlRowSolFeasibility(scip, nlrows[i], point, &feasibility) );
389  	
390  	         /* do not consider non-violated constraints */
391  	         if( SCIPisFeasGE(scip, feasibility, 0.0) )
392  	            continue;
393  	
394  	         /* increase number of violated nlrows */
395  	         ++nviolnlrows;
396  	
397  	         SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrows[i], point, &activity) );
398  	         SCIP_CALL( computeGradient(scip, nlrows[i], point, varindex, exprit, grad, &nlrownorm) );
399  	
400  	         /* update estimated costs for computing gradients */
401  	         *gradcosts += nlrowgradcosts[i];
402  	
403  	         /* stop if the gradient disappears at the current point */
404  	         if( SCIPisZero(scip, nlrownorm) )
405  	         {
406  	#ifdef SCIP_DEBUG_IMPROVEPOINT
407  	            printf("gradient vanished at current point -> stop\n");
408  	#endif
409  	            goto TERMINATE;
410  	         }
411  	
412  	         /* compute -g(x_k) / ||grad(g)(x_k)||^2 for a constraint g(x_k) <= 0 */
413  	         scale = -feasibility / nlrownorm;
414  	         if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrows[i])) && SCIPisGT(scip, activity, SCIPnlrowGetRhs(nlrows[i])) )
415  	            scale *= -1.0;
416  	
417  	         /* skip nonliner row if the scaler is too small or too large */
418  	         if( SCIPisEQ(scip, scale, 0.0) || SCIPisHugeValue(scip, REALABS(scale)) )
419  	            continue;
420  	
421  	         for( j = 0; j < nvars; ++j )
422  	            updatevec[j] += scale * grad[j];
423  	      }
424  	
425  	      /* if there are no violated rows, stop since start point is feasible */
426  	      if( nviolnlrows == 0 )
427  	      {
428  	         assert(updatevec[i] == 0.0);
429  	         return SCIP_OKAY;
430  	      }
431  	
432  	      for( i = 0; i < nvars; ++i )
433  	      {
434  	         /* adjust point */
435  	         updatevec[i] = SCIPgetSolVal(scip, point, vars[i]) + updatevec[i] / nviolnlrows;
436  	         updatevec[i] = MIN(updatevec[i], SCIPvarGetUbLocal(vars[i])); /*lint !e666*/
437  	         updatevec[i] = MAX(updatevec[i], SCIPvarGetLbLocal(vars[i])); /*lint !e666*/
438  	
439  	         SCIP_CALL( SCIPsetSolVal(scip, point, vars[i], updatevec[i]) );
440  	      }
441  	
442  	      /* update feasibility */
443  	      SCIP_CALL( getMinFeas(scip, nlrows, nnlrows, point, minfeas) );
444  	
445  	      /* check stopping criterion */
446  	      if( r % minimpriter == 0 && r > 0 )
447  	      {
448  	         if( *minfeas <= MINFEAS
449  	            || (*minfeas-lastminfeas) / MAX(REALABS(*minfeas), REALABS(lastminfeas)) < minimprfac ) /*lint !e666*/
450  	            break;
451  	         lastminfeas = *minfeas;
452  	      }
453  	   }
454  	
455  	TERMINATE:
456  	#ifdef SCIP_DEBUG_IMPROVEPOINT
457  	   printf("niter=%d minfeas=%e\n", r, *minfeas);
458  	#endif
459  	
460  	   SCIPfreeExpriter(&exprit);
461  	
462  	   SCIPfreeBufferArray(scip, &updatevec);
463  	   SCIPfreeBufferArray(scip, &grad);
464  	
465  	   return SCIP_OKAY;
466  	}
467  	
468  	/** sorts points w.r.t their feasibilities; points with a feasibility which is too small (w.r.t. the geometric mean of
469  	 *  all feasibilities) will be filtered out
470  	 */
471  	static
472  	SCIP_RETCODE filterPoints(
473  	   SCIP*                 scip,               /**< SCIP data structure */
474  	   SCIP_SOL**            points,             /**< array containing improved points */
475  	   SCIP_Real*            feasibilities,      /**< array containing feasibility for each point (sorted) */
476  	   int                   npoints,            /**< total number of points */
477  	   int*                  nusefulpoints       /**< pointer to store the total number of useful points */
478  	   )
479  	{
480  	   SCIP_Real minfeas;
481  	   SCIP_Real meanfeas;
482  	   int i;
483  	
484  	   assert(points != NULL);
485  	   assert(feasibilities != NULL);
486  	   assert(npoints > 0);
487  	   assert(nusefulpoints != NULL);
488  	
489  	   /* sort points w.r.t their feasibilities; non-negative feasibility correspond to feasible points for the NLP */
490  	   SCIPsortDownRealPtr(feasibilities, (void**)points, npoints);
491  	   minfeas = feasibilities[npoints - 1];
492  	
493  	   /* check if all points are feasible */
494  	   if( SCIPisFeasGE(scip, minfeas, 0.0) )
495  	   {
496  	      *nusefulpoints = npoints;
497  	      return SCIP_OKAY;
498  	   }
499  	
500  	   *nusefulpoints = 0;
501  	
502  	   /* compute shifted geometric mean of feasibilities (shift value = 1 - minfeas) */
503  	   meanfeas = 1.0;
504  	   for( i = 0; i < npoints; ++i )
505  	   {
506  	      assert(feasibilities[i] - minfeas + 1.0 > 0.0);
507  	      meanfeas *= pow(feasibilities[i] - minfeas + 1.0, 1.0 / npoints);
508  	   }
509  	   meanfeas += minfeas - 1.0;
510  	   SCIPdebugMsg(scip, "meanfeas = %e\n", meanfeas);
511  	
512  	   /* keep all points with which have a feasibility not much below the geometric mean of infeasibilities */
513  	   for( i = 0; i < npoints; ++i )
514  	   {
515  	      if( SCIPisFeasLT(scip, feasibilities[i], 0.0)
516  	         && (feasibilities[i] <= 1.05 * meanfeas || SCIPisLE(scip, feasibilities[i], MINFEAS)) )
517  	         break;
518  	
519  	      ++(*nusefulpoints);
520  	   }
521  	
522  	   return SCIP_OKAY;
523  	}
524  	
525  	/** returns the relative distance between two points; considers a smaller bounded domain for unbounded variables */
526  	static
527  	SCIP_Real getRelDistance(
528  	   SCIP*                 scip,               /**< SCIP data structure */
529  	   SCIP_SOL*             x,                  /**< first point */
530  	   SCIP_SOL*             y,                  /**< second point */
531  	   SCIP_Real             maxboundsize        /**< maximum variable domain size for unbounded variables */
532  	   )
533  	{
534  	   SCIP_VAR** vars;
535  	   int nvars;
536  	   SCIP_Real distance;
537  	   SCIP_Real solx;
538  	   SCIP_Real soly;
539  	   SCIP_Real lb;
540  	   SCIP_Real ub;
541  	   int i;
542  	
543  	   assert(x != NULL);
544  	   assert(y != NULL);
545  	
546  	   vars = SCIPgetVars(scip);
547  	   nvars = SCIPgetNVars(scip);
548  	   distance = 0.0;
549  	
550  	   if( nvars == 0 )
551  	      return 0.0;
552  	
553  	   for( i = 0; i < nvars; ++i )
554  	   {
555  	      lb = SCIPvarGetLbLocal(vars[i]);
556  	      ub = SCIPvarGetUbLocal(vars[i]);
557  	      solx = SCIPgetSolVal(scip, x, vars[i]);
558  	      soly = SCIPgetSolVal(scip, y, vars[i]);
559  	
560  	      /* adjust lower and upper bounds for unbounded variables*/
561  	      if( SCIPisInfinity(scip, -lb) && SCIPisInfinity(scip, ub) )
562  	      {
563  	         lb = -maxboundsize / 2.0;
564  	         ub = +maxboundsize / 2.0;
565  	      }
566  	      else if( SCIPisInfinity(scip, -lb) )
567  	      {
568  	         lb = ub - maxboundsize;
569  	      }
570  	      else if( SCIPisInfinity(scip, ub) )
571  	      {
572  	         ub = lb + maxboundsize;
573  	      }
574  	
575  	      /* project solution values to the variable domain */
576  	      solx = MIN(MAX(solx, lb), ub);
577  	      soly = MIN(MAX(soly, lb), ub);
578  	
579  	      distance += REALABS(solx - soly) / MAX(1.0, ub - lb);
580  	   }
581  	
582  	   return distance / nvars;
583  	}
584  	
585  	/** cluster useful points with a greedy algorithm */
586  	static
587  	SCIP_RETCODE clusterPointsGreedy(
588  	   SCIP*                 scip,               /**< SCIP data structure */
589  	   SCIP_SOL**            points,             /**< array containing improved points */
590  	   int                   npoints,            /**< total number of points */
591  	   int*                  clusteridx,         /**< array to store for each point the index of the cluster */
592  	   int*                  ncluster,           /**< pointer to store the total number of cluster */
593  	   SCIP_Real             maxboundsize,       /**< maximum variable domain size for unbounded variables */
594  	   SCIP_Real             maxreldist,         /**< maximum relative distance between any two points of the same cluster */
595  	   int                   maxncluster         /**< maximum number of clusters to compute */
596  	   )
597  	{
598  	   int i;
599  	
600  	   assert(points != NULL);
601  	   assert(npoints > 0);
602  	   assert(clusteridx != NULL);
603  	   assert(ncluster != NULL);
604  	   assert(maxreldist >= 0.0);
605  	   assert(maxncluster >= 0);
606  	
607  	   /* initialize cluster indices */
608  	   for( i = 0; i < npoints; ++i )
609  	      clusteridx[i] = INT_MAX;
610  	
611  	   *ncluster = 0;
612  	
613  	   for( i = 0; i < npoints && (*ncluster < maxncluster); ++i )
614  	   {
615  	      int j;
616  	
617  	      /* point is already assigned to a cluster */
618  	      if( clusteridx[i] != INT_MAX )
619  	         continue;
620  	
621  	      /* create a new cluster for i */
622  	      clusteridx[i] = *ncluster;
623  	
624  	      for( j = i + 1; j < npoints; ++j )
625  	      {
626  	         if( clusteridx[j] == INT_MAX && getRelDistance(scip, points[i], points[j], maxboundsize) <= maxreldist )
627  	            clusteridx[j] = *ncluster;
628  	      }
629  	
630  	      ++(*ncluster);
631  	   }
632  	
633  	#ifndef NDEBUG
634  	   for( i = 0; i < npoints; ++i )
635  	   {
636  	      assert(clusteridx[i] >= 0);
637  	      assert(clusteridx[i] < *ncluster || clusteridx[i] == INT_MAX);
638  	   }
639  	#endif
640  	
641  	   return SCIP_OKAY;
642  	}
643  	
644  	/** calls the sub-NLP heuristic for a given cluster */
645  	static
646  	SCIP_RETCODE solveNLP(
647  	   SCIP*                 scip,               /**< SCIP data structure */
648  	   SCIP_HEUR*            heur,               /**< multi-start heuristic */
649  	   SCIP_HEUR*            nlpheur,            /**< pointer to NLP local search heuristics */
650  	   SCIP_SOL**            points,             /**< array containing improved points */
651  	   int                   npoints,            /**< total number of points */
652  	   SCIP_Bool*            success             /**< pointer to store if we could find a solution */
653  	   )
654  	{
655  	   SCIP_VAR** vars;
656  	   SCIP_SOL* refpoint;
657  	   SCIP_RESULT nlpresult;
658  	   SCIP_Real val;
659  	   int nbinvars;
660  	   int nintvars;
661  	   int nvars;
662  	   int i;
663  	
664  	   assert(points != NULL);
665  	   assert(npoints > 0);
666  	
667  	   SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, &nbinvars, &nintvars, NULL, NULL) );
668  	   *success = FALSE;
669  	
670  	   SCIP_CALL( SCIPcreateSol(scip, &refpoint, heur) );
671  	
672  	   /* compute reference point */
673  	   for( i = 0; i < nvars; ++i )
674  	   {
675  	      int p;
676  	
677  	      val = 0.0;
678  	
679  	      for( p = 0; p < npoints; ++p )
680  	      {
681  	         assert(points[p] != NULL);
682  	         val += SCIPgetSolVal(scip, points[p], vars[i]);
683  	      }
684  	
685  	      SCIP_CALL( SCIPsetSolVal(scip, refpoint, vars[i], val / npoints) );
686  	   }
687  	
688  	   /* round point for sub-NLP heuristic */
689  	   SCIP_CALL( SCIProundSol(scip, refpoint, success) );
690  	   SCIPdebugMsg(scip, "rounding of refpoint successfully? %u\n", *success);
691  	
692  	   /* round variables manually if the locks did not allow us to round them */
693  	   if( !(*success) )
694  	   {
695  	      for( i = 0; i < nbinvars + nintvars; ++i )
696  	      {
697  	         val = SCIPgetSolVal(scip, refpoint, vars[i]);
698  	
699  	         if( !SCIPisFeasIntegral(scip, val) )
700  	         {
701  	            assert(SCIPisFeasIntegral(scip, SCIPvarGetLbLocal(vars[i])));
702  	            assert(SCIPisFeasIntegral(scip, SCIPvarGetUbLocal(vars[i])));
703  	
704  	            /* round and adjust value */
705  	            val = SCIPround(scip, val);
706  	            val = MIN(val, SCIPvarGetUbLocal(vars[i])); /*lint !e666*/
707  	            val = MAX(val, SCIPvarGetLbLocal(vars[i])); /*lint !e666*/
708  	            assert(SCIPisFeasIntegral(scip, val));
709  	
710  	            SCIP_CALL( SCIPsetSolVal(scip, refpoint, vars[i], val) );
711  	         }
712  	      }
713  	   }
714  	
715  	   /* call sub-NLP heuristic */
716  	   SCIP_CALL( SCIPapplyHeurSubNlp(scip, nlpheur, &nlpresult, refpoint, NULL) );
717  	   SCIP_CALL( SCIPfreeSol(scip, &refpoint) );
718  	
719  	   /* let sub-NLP heuristic decide whether the solution is feasible or not */
720  	   *success = nlpresult == SCIP_FOUNDSOL;
721  	
722  	   return SCIP_OKAY;
723  	}
724  	
725  	/** recursive helper function to count the number of nodes in a sub-expr */
726  	static
727  	int getExprSize(
728  	   SCIP_EXPR*            expr                /**< expression */
729  	   )
730  	{
731  	   int sum;
732  	   int i;
733  	
734  	   assert(expr != NULL);
735  	
736  	   sum = 0;
737  	   for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
738  	   {
739  	      SCIP_EXPR* child = SCIPexprGetChildren(expr)[i];
740  	      sum += getExprSize(child);
741  	   }
742  	   return 1 + sum;
743  	}
744  	
745  	/** main function of the multi-start heuristic (see @ref heur_multistart.h for more details); it consists of the
746  	 *  following four steps:
747  	 *
748  	 *  1. sampling points in the current domain; for unbounded variables we use a bounded box
749  	 *
750  	 *  2. reduce infeasibility by using a gradient descent method
751  	 *
752  	 *  3. cluster points; filter points with a too large infeasibility
753  	 *
754  	 *  4. compute start point for each cluster and use it in the sub-NLP heuristic (@ref heur_subnlp.h)
755  	 */
756  	static
757  	SCIP_RETCODE applyHeur(
758  	   SCIP*                 scip,               /**< SCIP data structure */
759  	   SCIP_HEUR*            heur,               /**< heuristic */
760  	   SCIP_HEURDATA*        heurdata,           /**< heuristic data */
761  	   SCIP_RESULT*          result              /**< pointer to store the result */
762  	   )
763  	{
764  	   SCIP_NLROW** nlrows;
765  	   SCIP_SOL** points;
766  	   SCIP_HASHMAP* varindex;
767  	   SCIP_Real* feasibilities;
768  	   SCIP_Real* nlrowgradcosts;
769  	   int* clusteridx;
770  	   SCIP_Real gradlimit;
771  	   SCIP_Real bestobj;
772  	   int nusefulpoints;
773  	   int nrndpoints;
774  	   int ncluster;
775  	   int nnlrows;
776  	   int npoints;
777  	   int start;
778  	   int i;
779  	
780  	   assert(scip != NULL);
781  	   assert(heur != NULL);
782  	   assert(result != NULL);
783  	   assert(heurdata != NULL);
784  	
785  	   SCIPdebugMsg(scip, "call applyHeur()\n");
786  	
787  	   nlrows = SCIPgetNLPNlRows(scip);
788  	   nnlrows = SCIPgetNNLPNlRows(scip);
789  	   bestobj = SCIPgetNSols(scip) > 0 ? MINIMPRFAC * SCIPgetSolTransObj(scip, SCIPgetBestSol(scip)) : SCIPinfinity(scip);
790  	
791  	   SCIP_CALL( SCIPallocBufferArray(scip, &points, heurdata->nrndpoints) );
792  	   SCIP_CALL( SCIPallocBufferArray(scip, &nlrowgradcosts, nnlrows) );
793  	   SCIP_CALL( SCIPallocBufferArray(scip, &feasibilities, heurdata->nrndpoints) );
794  	   SCIP_CALL( SCIPallocBufferArray(scip, &clusteridx, heurdata->nrndpoints) );
795  	   SCIP_CALL( SCIPhashmapCreate(&varindex, SCIPblkmem(scip), SCIPgetNVars(scip)) );
796  	
797  	   /* create an unique mapping of all variables to 0,..,SCIPgetNVars(scip)-1 */
798  	   for( i = 0; i < SCIPgetNVars(scip); ++i )
799  	   {
800  	      SCIP_CALL( SCIPhashmapInsertInt(varindex, (void*)SCIPgetVars(scip)[i], i) );
801  	   }
802  	
803  	   /* compute estimated costs of computing a gradient for each nlrow */
804  	   for( i = 0; i < nnlrows; ++i )
805  	   {
806  	      nlrowgradcosts[i] = GRADCOSTFAC_LINEAR * SCIPnlrowGetNLinearVars(nlrows[i]);
807  	      if( SCIPnlrowGetExpr(nlrows[i]) != NULL )
808  	         nlrowgradcosts[i] += GRADCOSTFAC_NONLINEAR * getExprSize(SCIPnlrowGetExpr(nlrows[i]));
809  	   }
810  	
811  	   /*
812  	    * 1. sampling points in the current domain; for unbounded variables we use a bounded box
813  	    */
814  	   SCIP_CALL( sampleRandomPoints(scip, points, heurdata->nrndpoints, heurdata->maxboundsize, heurdata->randnumgen,
815  	         bestobj, &nrndpoints) );
816  	   assert(nrndpoints >= 0);
817  	
818  	   if( nrndpoints == 0 )
819  	      goto TERMINATE;
820  	
821  	   /*
822  	    * 2. improve points via consensus vectors
823  	    */
824  	   gradlimit = heurdata->gradlimit == 0.0 ? SCIPinfinity(scip) : heurdata->gradlimit;
825  	   for( npoints = 0; npoints < nrndpoints && gradlimit >= 0 && !SCIPisStopped(scip); ++npoints )
826  	   {
827  	      SCIP_Real gradcosts;
828  	
829  	      SCIP_CALL( improvePoint(scip, nlrows, nnlrows, varindex, points[npoints],
830  	            heurdata->maxiter, heurdata->minimprfac, heurdata->minimpriter, &feasibilities[npoints], nlrowgradcosts,
831  	            &gradcosts) );
832  	
833  	      gradlimit -= gradcosts;
834  	      SCIPdebugMsg(scip, "improve point %d / %d gradlimit = %g\n", npoints, nrndpoints, gradlimit);
835  	   }
836  	   assert(npoints >= 0 && npoints <= nrndpoints);
837  	
838  	   if( npoints == 0 )
839  	      goto TERMINATE;
840  	
841  	   /*
842  	    * 3. filter and cluster points
843  	    */
844  	   SCIP_CALL( filterPoints(scip, points, feasibilities, npoints, &nusefulpoints) );
845  	   assert(nusefulpoints >= 0);
846  	   SCIPdebugMsg(scip, "nusefulpoints = %d\n", nusefulpoints);
847  	
848  	   if( nusefulpoints == 0 )
849  	      goto TERMINATE;
850  	
851  	   SCIP_CALL( clusterPointsGreedy(scip, points, nusefulpoints, clusteridx, &ncluster, heurdata->maxboundsize,
852  	         heurdata->maxreldist, heurdata->maxncluster) );
853  	   assert(ncluster >= 0 && ncluster <= heurdata->maxncluster);
854  	   SCIPdebugMsg(scip, "ncluster = %d\n", ncluster);
855  	
856  	   SCIPsortIntPtr(clusteridx, (void**)points, nusefulpoints);
857  	
858  	   /*
859  	    * 4. compute start point for each cluster and use it in the sub-NLP heuristic (@ref heur_subnlp.h)
860  	    */
861  	   start = 0;
862  	   while( start < nusefulpoints && clusteridx[start] != INT_MAX && !SCIPisStopped(scip) )
863  	   {
864  	      SCIP_Bool success;
865  	      int end;
866  	
867  	      end = start;
868  	      while( end < nusefulpoints && clusteridx[start] == clusteridx[end] )
869  	         ++end;
870  	
871  	      assert(end - start > 0);
872  	
873  	      /* call sub-NLP heuristic */
874  	      SCIP_CALL( solveNLP(scip, heur, heurdata->heursubnlp, &points[start], end - start, &success) );
875  	      SCIPdebugMsg(scip, "solveNLP result = %u\n", success);
876  	
877  	      if( success )
878  	         *result = SCIP_FOUNDSOL;
879  	
880  	      /* go to the next cluster */
881  	      start = end;
882  	   }
883  	
884  	TERMINATE:
885  	   /* free memory */
886  	   for( i = nrndpoints - 1; i >= 0 ; --i )
887  	   {
888  	      assert(points[i] != NULL);
889  	      SCIP_CALL( SCIPfreeSol(scip, &points[i]) );
890  	   }
891  	
892  	   SCIPhashmapFree(&varindex);
893  	   SCIPfreeBufferArray(scip, &clusteridx);
894  	   SCIPfreeBufferArray(scip, &feasibilities);
895  	   SCIPfreeBufferArray(scip, &nlrowgradcosts);
896  	   SCIPfreeBufferArray(scip, &points);
897  	
898  	   return SCIP_OKAY;
899  	}
900  	
901  	/*
902  	 * Callback methods of primal heuristic
903  	 */
904  	
905  	/** copy method for primal heuristic plugins (called when SCIP copies plugins) */
906  	static
907  	SCIP_DECL_HEURCOPY(heurCopyMultistart)
908  	{  /*lint --e{715}*/
909  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
910  	
911  	   /* call inclusion method of primal heuristic */
912  	   SCIP_CALL( SCIPincludeHeurMultistart(scip) );
913  	
914  	   return SCIP_OKAY;
915  	}
916  	
917  	/** destructor of primal heuristic to free user data (called when SCIP is exiting) */
918  	static
919  	SCIP_DECL_HEURFREE(heurFreeMultistart)
920  	{  /*lint --e{715}*/
921  	   SCIP_HEURDATA* heurdata;
922  	
923  	   /* free heuristic data */
924  	   heurdata = SCIPheurGetData(heur);
925  	
926  	   SCIPfreeBlockMemory(scip, &heurdata);
927  	   SCIPheurSetData(heur, NULL);
928  	
929  	   return SCIP_OKAY;
930  	}
931  	
932  	/** initialization method of primal heuristic (called after problem was transformed) */
933  	static
934  	SCIP_DECL_HEURINIT(heurInitMultistart)
935  	{  /*lint --e{715}*/
936  	   SCIP_HEURDATA* heurdata;
937  	
938  	   assert( heur != NULL );
939  	
940  	   heurdata = SCIPheurGetData(heur);
941  	   assert(heurdata != NULL);
942  	
943  	   SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen,
944  	         DEFAULT_RANDSEED, TRUE) );
945  	
946  	   /* try to find sub-NLP heuristic */
947  	   heurdata->heursubnlp = SCIPfindHeur(scip, "subnlp");
948  	
949  	   return SCIP_OKAY;
950  	}
951  	
952  	/** deinitialization method of primal heuristic (called before transformed problem is freed) */
953  	static
954  	SCIP_DECL_HEUREXIT(heurExitMultistart)
955  	{  /*lint --e{715}*/
956  	   SCIP_HEURDATA* heurdata;
957  	
958  	   assert( heur != NULL );
959  	
960  	   heurdata = SCIPheurGetData(heur);
961  	   assert(heurdata != NULL);
962  	   assert(heurdata->randnumgen != NULL);
963  	
964  	   SCIPfreeRandom(scip, &heurdata->randnumgen);
965  	
966  	   return SCIP_OKAY;
967  	}
968  	
969  	/** execution method of primal heuristic */
970  	static
971  	SCIP_DECL_HEUREXEC(heurExecMultistart)
972  	{  /*lint --e{715}*/
973  	   SCIP_HEURDATA* heurdata;
974  	
975  	   assert( heur != NULL );
976  	
977  	   heurdata = SCIPheurGetData(heur);
978  	   assert(heurdata != NULL);
979  	
980  	   *result = SCIP_DIDNOTRUN;
981  	
982  	   /* check cases for which the heuristic is not applicable */
983  	   if( !SCIPisNLPConstructed(scip) || heurdata->heursubnlp == NULL || SCIPgetNNlpis(scip) <= 0 )
984  	      return SCIP_OKAY;
985  	
986  	   /* check whether the heuristic should be applied for a problem containing integer variables */
987  	   if( heurdata->onlynlps && (SCIPgetNBinVars(scip) > 0 || SCIPgetNIntVars(scip) > 0) )
988  	      return SCIP_OKAY;
989  	
990  	   *result = SCIP_DIDNOTFIND;
991  	
992  	   SCIP_CALL( applyHeur(scip, heur, heurdata, result) );
993  	
994  	   return SCIP_OKAY;
995  	}
996  	
997  	/*
998  	 * primal heuristic specific interface methods
999  	 */
1000 	
1001 	/** creates the multistart primal heuristic and includes it in SCIP */
1002 	SCIP_RETCODE SCIPincludeHeurMultistart(
1003 	   SCIP*                 scip                /**< SCIP data structure */
1004 	   )
1005 	{
1006 	   SCIP_HEURDATA* heurdata;
1007 	   SCIP_HEUR* heur;
1008 	
1009 	   /* create multistart primal heuristic data */
1010 	   SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) );
1011 	   BMSclearMemory(heurdata);
1012 	
1013 	   /* include primal heuristic */
1014 	   SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
1015 	         HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
1016 	         HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecMultistart, heurdata) );
1017 	
1018 	   assert(heur != NULL);
1019 	
1020 	   /* set non fundamental callbacks via setter functions */
1021 	   SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyMultistart) );
1022 	   SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeMultistart) );
1023 	   SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitMultistart) );
1024 	   SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitMultistart) );
1025 	
1026 	   /* add multistart primal heuristic parameters */
1027 	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/nrndpoints",
1028 	         "number of random points generated per execution call",
1029 	         &heurdata->nrndpoints, FALSE, DEFAULT_NRNDPOINTS, 0, INT_MAX, NULL, NULL) );
1030 	
1031 	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxboundsize",
1032 	         "maximum variable domain size for unbounded variables",
1033 	         &heurdata->maxboundsize, FALSE, DEFAULT_MAXBOUNDSIZE, 0.0, SCIPinfinity(scip), NULL, NULL) );
1034 	
1035 	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxiter",
1036 	         "number of iterations to reduce the maximum violation of a point",
1037 	         &heurdata->maxiter, FALSE, DEFAULT_MAXITER, 0, INT_MAX, NULL, NULL) );
1038 	
1039 	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/minimprfac",
1040 	         "minimum required improving factor to proceed in improvement of a single point",
1041 	         &heurdata->minimprfac, FALSE, DEFAULT_MINIMPRFAC, -SCIPinfinity(scip), SCIPinfinity(scip), NULL, NULL) );
1042 	
1043 	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/minimpriter",
1044 	         "number of iteration when checking the minimum improvement",
1045 	         &heurdata->minimpriter, FALSE, DEFAULT_MINIMPRITER, 1, INT_MAX, NULL, NULL) );
1046 	
1047 	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxreldist",
1048 	         "maximum distance between two points in the same cluster",
1049 	         &heurdata->maxreldist, FALSE, DEFAULT_MAXRELDIST, 0.0, SCIPinfinity(scip), NULL, NULL) );
1050 	
1051 	   SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/gradlimit",
1052 	         "limit for gradient computations for all improvePoint() calls (0 for no limit)",
1053 	         &heurdata->gradlimit, FALSE, DEFAULT_GRADLIMIT, 0.0, SCIPinfinity(scip), NULL, NULL) );
1054 	
1055 	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/maxncluster",
1056 	         "maximum number of considered clusters per heuristic call",
1057 	         &heurdata->maxncluster, FALSE, DEFAULT_MAXNCLUSTER, 0, INT_MAX, NULL, NULL) );
1058 	
1059 	   SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/onlynlps",
1060 	         "should the heuristic run only on continuous problems?",
1061 	         &heurdata->onlynlps, FALSE, DEFAULT_ONLYNLPS, NULL, NULL) );
1062 	
1063 	   return SCIP_OKAY;
1064 	}
1065