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   sepa_convexproj.c
26   	 * @ingroup DEFPLUGINS_SEPA
27   	 * @brief  convexproj separator
28   	 * @author Felipe Serrano
29   	 *
30   	 * @todo should separator only be run when SCIPallColsInLP is true?
31   	 * @todo check if it makes sense to implement the copy callback
32   	 * @todo add SCIPisStopped(scip) to the condition of time consuming loops
33   	 */
34   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
35   	
36   	#include <assert.h>
37   	#include <string.h>
38   	
39   	#include "blockmemshell/memory.h"
40   	#include "scip/scip_expr.h"
41   	#include "scip/scip_nlpi.h"
42   	#include "scip/expr_varidx.h"
43   	#include "scip/expr_pow.h"
44   	#include "scip/expr_sum.h"
45   	#include "scip/pub_message.h"
46   	#include "scip/pub_misc.h"
47   	#include "scip/pub_nlp.h"
48   	#include "scip/pub_sepa.h"
49   	#include "scip/pub_var.h"
50   	#include "scip/scip_cut.h"
51   	#include "scip/scip_general.h"
52   	#include "scip/scip_lp.h"
53   	#include "scip/scip_mem.h"
54   	#include "scip/scip_message.h"
55   	#include "scip/scip_nlp.h"
56   	#include "scip/scip_numerics.h"
57   	#include "scip/scip_param.h"
58   	#include "scip/scip_prob.h"
59   	#include "scip/scip_sepa.h"
60   	#include "scip/scip_sol.h"
61   	#include "scip/scip_solvingstats.h"
62   	#include "scip/scip_timing.h"
63   	#include "scip/scip_tree.h"
64   	#include "scip/sepa_convexproj.h"
65   	#include <string.h>
66   	
67   	
68   	#define SEPA_NAME              "convexproj"
69   	#define SEPA_DESC              "separate at projection of point onto convex region"
70   	#define SEPA_PRIORITY                 0
71   	#define SEPA_FREQ                    -1
72   	#define SEPA_MAXBOUNDDIST           1.0
73   	#define SEPA_USESSUBSCIP          FALSE      /**< does the separator use a secondary SCIP instance? */
74   	#define SEPA_DELAY                 TRUE      /**< should separation method be delayed, if other separators found cuts? */
75   	
76   	#define DEFAULT_MAXDEPTH             -1      /**< maximum depth at which the separator is applied; -1 means no limit */
77   	#define DEFAULT_NLPITERLIM          250      /**< default NLP iteration limit */
78   	
79   	#define VIOLATIONFAC                100      /**< points regarded violated if max violation > VIOLATIONFAC*SCIPfeastol() */
80   	
81   	/*
82   	 * Data structures
83   	 */
84   	
85   	/** side that makes an nlrow convex */
86   	enum ConvexSide
87   	{
88   	   LHS = 0,                                  /**< left hand side */
89   	   RHS = 1                                   /**< right hand side */
90   	};
91   	typedef enum ConvexSide CONVEXSIDE;
92   	
93   	/** separator data
94   	 * it keeps the nlpi which represents the projection problem (see sepa_convexproj.h); it also keeps the convex nlrows
95   	 * and the side which actually makes them convex; when separating, we use the nlpi to compute the projection and then
96   	 * the convex nlrows to compute the actual gradient cuts */
97   	struct SCIP_SepaData
98   	{
99   	   SCIP_NLPI*            nlpi;               /**< nlpi used to create the nlpi problem */
100  	   SCIP_NLPIPROBLEM*     nlpiprob;           /**< nlpi problem representing the convex NLP relaxation */
101  	   SCIP_VAR**            nlpivars;           /**< array containing all variables of the nlpi */
102  	   SCIP_HASHMAP*         var2nlpiidx;        /**< mapping between variables and nlpi indices */
103  	   int                   nlpinvars;          /**< total number of nlpi variables */
104  	
105  	   SCIP_Bool             skipsepa;           /**< should separator be skipped? */
106  	
107  	   SCIP_NLROW**          nlrows;             /**< convex nlrows */
108  	   CONVEXSIDE*           convexsides;        /**< which sides make the nlrows convex */
109  	   SCIP_Real*            constraintviolation;/**< array storing the violation of constraint by current solution; 0.0 if it is not violated */
110  	   int                   nnlrows;            /**< total number of nlrows */
111  	   int                   nlrowssize;         /**< memory allocated for nlrows, convexsides and nlrowsidx */
112  	
113  	   /* parameter */
114  	   int                   nlpiterlimit;       /**< iteration limit of NLP solver; 0 for no limit */
115  	   int                   maxdepth;           /**< maximal depth at which the separator is applied */
116  	
117  	   int                   ncuts;              /**< number of cuts generated */
118  	};
119  	
120  	
121  	/*
122  	 * Local methods
123  	 */
124  	
125  	/** clears the sepadata data */
126  	static
127  	SCIP_RETCODE sepadataClear(
128  	   SCIP*                 scip,               /**< SCIP data structure */
129  	   SCIP_SEPADATA*        sepadata            /**< separator data */
130  	   )
131  	{
132  	   assert(sepadata != NULL);
133  	
134  	   /* nlrowssize gets allocated first and then its decided whether to create the nlpiprob */
135  	   if( sepadata->nlrowssize > 0 )
136  	   {
137  	      SCIPfreeBlockMemoryArray(scip, &sepadata->constraintviolation, sepadata->nlrowssize);
138  	      SCIPfreeBlockMemoryArray(scip, &sepadata->convexsides, sepadata->nlrowssize);
139  	      SCIPfreeBlockMemoryArray(scip, &sepadata->nlrows, sepadata->nlrowssize);
140  	      sepadata->nlrowssize = 0;
141  	   }
142  	
143  	   if( sepadata->nlpiprob != NULL )
144  	   {
145  	      assert(sepadata->nlpi != NULL);
146  	
147  	      SCIPfreeBlockMemoryArray(scip, &sepadata->nlpivars, sepadata->nlpinvars);
148  	
149  	      SCIPhashmapFree(&sepadata->var2nlpiidx);
150  	      SCIP_CALL( SCIPfreeNlpiProblem(scip, sepadata->nlpi, &sepadata->nlpiprob) );
151  	
152  	      sepadata->nlpinvars = 0;
153  	      sepadata->nnlrows = 0;
154  	   }
155  	   assert(sepadata->nlpinvars == 0);
156  	   assert(sepadata->nnlrows == 0);
157  	   assert(sepadata->nlrowssize == 0);
158  	
159  	   sepadata->skipsepa = FALSE;
160  	
161  	   return SCIP_OKAY;
162  	}
163  	
164  	/** computes gradient cut (linearization) of nlrow at projection */
165  	static
166  	SCIP_RETCODE generateCut(
167  	   SCIP*                 scip,               /**< SCIP data structure */
168  	   SCIP_SEPA*            sepa,               /**< the cut separator itself */
169  	   SCIP_SOL*             projection,         /**< point where we compute gradient cut */
170  	   SCIP_NLROW*           nlrow,              /**< constraint for which we generate gradient cut */
171  	   CONVEXSIDE            convexside,         /**< which side makes the nlrow convex */
172  	   SCIP_Real             activity,           /**< activity of constraint at projection */
173  	   SCIP_EXPRITER*        exprit,             /**< expression iterator that can be used */
174  	   SCIP_ROW**            row                 /**< storage for cut */
175  	   )
176  	{
177  	   char rowname[SCIP_MAXSTRLEN];
178  	   SCIP_SEPADATA* sepadata;
179  	   SCIP_Real gradx0; /* <grad f(x_0), x_0> */
180  	   SCIP_EXPR* expr;
181  	   int i;
182  	
183  	   assert(scip != NULL);
184  	   assert(sepa != NULL);
185  	   assert(nlrow != NULL);
186  	   assert(row != NULL);
187  	
188  	   sepadata = SCIPsepaGetData(sepa);
189  	
190  	   assert(sepadata != NULL);
191  	
192  	   gradx0 = 0.0;
193  	
194  	   /* an nlrow has a linear part and expression; ideally one would just build the gradient but we
195  	    * do not know if the different parts share variables or not, so we can't just build the gradient; for this reason
196  	    * we create the row right away and compute the gradients of each part independently and add them to the row; the
197  	    * row takes care to add coeffs corresponding to the same variable when they appear in different parts of the nlrow
198  	    * NOTE: a gradient cut is globally valid whenever the constraint from which it is deduced is globally valid; since
199  	    *       we build the convex relaxation using only globally valid constraints, the cuts are globally valid
200  	    */
201  	   (void) SCIPsnprintf(rowname, SCIP_MAXSTRLEN, "proj_cut_%s_%u", SCIPnlrowGetName(nlrow), ++(sepadata->ncuts));
202  	   SCIP_CALL( SCIPcreateEmptyRowSepa(scip, row, sepa, rowname, -SCIPinfinity(scip), SCIPinfinity(scip), TRUE, FALSE ,
203  	            TRUE) );
204  	
205  	   SCIP_CALL( SCIPcacheRowExtensions(scip, *row) );
206  	
207  	   /* linear part */
208  	   for( i = 0; i < SCIPnlrowGetNLinearVars(nlrow); i++ )
209  	   {
210  	      gradx0 += SCIPgetSolVal(scip, projection, SCIPnlrowGetLinearVars(nlrow)[i]) * SCIPnlrowGetLinearCoefs(nlrow)[i];
211  	      SCIP_CALL( SCIPaddVarToRow(scip, *row, SCIPnlrowGetLinearVars(nlrow)[i], SCIPnlrowGetLinearCoefs(nlrow)[i]) );
212  	   }
213  	
214  	   expr = SCIPnlrowGetExpr(nlrow);
215  	   assert(expr != NULL);
216  	
217  	   SCIP_CALL( SCIPevalExprGradient(scip, expr, projection, 0L) );
218  	
219  	   SCIP_CALL( SCIPexpriterInit(exprit, expr, SCIP_EXPRITER_DFS, FALSE) );
220  	   for( ; !SCIPexpriterIsEnd(exprit); expr = SCIPexpriterGetNext(exprit) )  /*lint !e441*/ /*lint !e440*/
221  	   {
222  	      SCIP_Real grad;
223  	      SCIP_VAR* var;
224  	
225  	      if( !SCIPisExprVar(scip, expr) )
226  	         continue;
227  	
228  	      grad = SCIPexprGetDerivative(expr);
229  	      var = SCIPgetVarExprVar(expr);
230  	      assert(var != NULL);
231  	
232  	      gradx0 += grad * SCIPgetSolVal(scip, projection, var);
233  	      SCIP_CALL( SCIPaddVarToRow(scip, *row, var, grad) );
234  	   }
235  	
236  	   SCIP_CALL( SCIPflushRowExtensions(scip, *row) );
237  	
238  	   SCIPdebugPrintf("gradient: ");
239  	   SCIPdebug( SCIP_CALL( SCIPprintRow(scip, *row, NULL) ) );
240  	   SCIPdebugPrintf("gradient dot x_0: %g\n", gradx0);
241  	
242  	   /* gradient cut is f(x_0) - <grad f(x_0), x_0> + <grad f(x_0), x> <= rhs or >= lhs */
243  	   if( convexside == RHS )
244  	   {
245  	      assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)));
246  	      SCIP_CALL( SCIPchgRowRhs(scip, *row, SCIPnlrowGetRhs(nlrow) - activity + gradx0) );
247  	   }
248  	   else
249  	   {
250  	      assert(convexside == LHS);
251  	      assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow)));
252  	      SCIP_CALL( SCIPchgRowLhs(scip, *row, SCIPnlrowGetLhs(nlrow) - activity + gradx0) );
253  	   }
254  	
255  	   SCIPdebugPrintf("gradient cut: ");
256  	   SCIPdebug( SCIP_CALL( SCIPprintRow(scip, *row, NULL) ) );
257  	
258  	   return SCIP_OKAY;
259  	}
260  	
261  	/** set quadratic part of objective function: \f$ \sum_i x_i^2 \f$
262  	 *
263  	 * the objective function is \f$ ||x - x_0||^2 \f$,
264  	 * where \f$ x_0 \f$ is the point to separate; the only part that changes is the term \f$ -2 \langle x_0, x \rangle \f$
265  	 * which is linear and is set every time we want to separate a point, see separateCuts()
266  	 */
267  	static
268  	SCIP_RETCODE setQuadraticObj(
269  	   SCIP*                 scip,               /**< SCIP data structure */
270  	   SCIP_SEPADATA*        sepadata            /**< the cut separator data */
271  	   )
272  	{
273  	   SCIP_EXPR* exprsum;
274  	   SCIP_EXPR** exprspow;
275  	   int i;
276  	
277  	   assert(scip != NULL);
278  	   assert(sepadata != NULL);
279  	   assert(sepadata->nlpi != NULL);
280  	   assert(sepadata->nlpiprob != NULL);
281  	   assert(sepadata->var2nlpiidx != NULL);
282  	   assert(sepadata->nlpinvars > 0);
283  	
284  	   SCIP_CALL( SCIPallocBufferArray(scip, &exprspow, sepadata->nlpinvars) );
285  	   for( i = 0; i < sepadata->nlpinvars; i++ )
286  	   {
287  	      SCIP_VAR* var;
288  	      SCIP_EXPR* varexpr;
289  	
290  	      var = sepadata->nlpivars[i];
291  	      assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) );
292  	
293  	      SCIP_CALL( SCIPcreateExprVaridx(scip, &varexpr, SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void*)var), NULL, NULL) );
294  	      SCIP_CALL( SCIPcreateExprPow(scip, &exprspow[i], varexpr, 2.0, NULL, NULL) );
295  	      SCIP_CALL( SCIPreleaseExpr(scip, &varexpr) );
296  	   }
297  	
298  	   SCIP_CALL( SCIPcreateExprSum(scip, &exprsum, sepadata->nlpinvars, exprspow, NULL, 0.0, NULL, NULL) );
299  	
300  	   /* set quadratic part of objective function */
301  	   SCIP_CALL( SCIPsetNlpiObjective(scip, sepadata->nlpi, sepadata->nlpiprob, 0, NULL, NULL, exprsum, 0.0) );
302  	
303  	   /* free memory */
304  	   SCIP_CALL( SCIPreleaseExpr(scip, &exprsum) );
305  	   for( i = sepadata->nlpinvars-1; i >= 0; --i )
306  	   {
307  	      SCIP_CALL( SCIPreleaseExpr(scip, &exprspow[i]) );
308  	   }
309  	   SCIPfreeBufferArray(scip, &exprspow);
310  	
311  	   return SCIP_OKAY;
312  	}
313  	
314  	/** projects sol onto convex relaxation (stored in sepadata) and tries to generate gradient cuts at the projection
315  	 *
316  	 * it generates cuts only for the constraints that were violated by the LP solution and are now active or still
317  	 * violated (in case we don't solve to optimality).
318  	 * @todo: store a feasible solution if one is found to use as warmstart
319  	 */
320  	static
321  	SCIP_RETCODE separateCuts(
322  	   SCIP*                 scip,               /**< SCIP data structure */
323  	   SCIP_SEPA*            sepa,               /**< the cut separator itself */
324  	   SCIP_SOL*             sol,                /**< solution that should be separated */
325  	   SCIP_RESULT*          result              /**< pointer to store the result of the separation call */
326  	   )
327  	{
328  	   SCIP_SEPADATA* sepadata;
329  	   SCIP_SOL*      projection;
330  	   SCIP_Real*     linvals;
331  	   SCIP_Real*     nlpisol;
332  	   int            nlpinvars;
333  	   int            i;
334  	   int*           lininds;
335  	   SCIP_Bool      nlpunstable;
336  	   SCIP_EXPRITER* exprit;
337  	
338  	   nlpunstable = FALSE;
339  	
340  	   assert(sepa != NULL);
341  	
342  	   sepadata = SCIPsepaGetData(sepa);
343  	
344  	   assert(result != NULL);
345  	   assert(sepadata != NULL);
346  	   assert(sepadata->nnlrows > 0);
347  	   assert(sepadata->nlpi != NULL);
348  	   assert(sepadata->nlpinvars > 0);
349  	   assert(sepadata->nlrows != NULL);
350  	   assert(sepadata->nlpiprob != NULL);
351  	   assert(sepadata->var2nlpiidx != NULL);
352  	   assert(sepadata->convexsides != NULL);
353  	   assert(sepadata->constraintviolation != NULL);
354  	
355  	   nlpinvars = sepadata->nlpinvars;
356  	   /* set linear part of objective function: \norm(x - x^0)^2 = \norm(x)^2 - \sum 2 * x_i * x^0_i + const
357  	    * we ignore the constant; x0 is `sol`
358  	    */
359  	   SCIP_CALL( SCIPallocBufferArray(scip, &linvals, nlpinvars) );
360  	   SCIP_CALL( SCIPallocBufferArray(scip, &lininds, nlpinvars) );
361  	   for( i = 0; i < nlpinvars; i++ )
362  	   {
363  	      SCIP_VAR* var;
364  	
365  	      var = sepadata->nlpivars[i];
366  	      assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) );
367  	
368  	      lininds[i] = SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void*)var);
369  	      linvals[i] = - 2.0 * SCIPgetSolVal(scip, sol, var);
370  	
371  	      /* if coefficient is too large, don't separate */
372  	      if( SCIPisHugeValue(scip, REALABS(linvals[i])) )
373  	      {
374  	         SCIPdebugMsg(scip, "Don't separate points too close to infinity\n");
375  	         goto CLEANUP;
376  	      }
377  	   }
378  	
379  	   /* set linear part of objective function */
380  	   SCIP_CALL( SCIPchgNlpiLinearCoefs(scip, sepadata->nlpi, sepadata->nlpiprob, -1, nlpinvars, lininds, linvals) );
381  	
382  	   /* compute the projection onto the convex NLP relaxation */
383  	   SCIP_CALL( SCIPsolveNlpi(scip, sepadata->nlpi, sepadata->nlpiprob,
384  	      .iterlimit = sepadata->nlpiterlimit > 0 ? sepadata->nlpiterlimit : INT_MAX,
385  	      .feastol = SCIPfeastol(scip) / 10.0, /* use tighter tolerances for the NLP solver */
386  	      .opttol = MAX(SCIPfeastol(scip), SCIPdualfeastol(scip))) );  /*lint !e666*/
387  	   SCIPdebugMsg(scip, "NLP solstat = %d\n", SCIPgetNlpiSolstat(scip, sepadata->nlpi, sepadata->nlpiprob));
388  	
389  	   /* if solution is feasible, add cuts */
390  	   switch( SCIPgetNlpiSolstat(scip, sepadata->nlpi, sepadata->nlpiprob) )
391  	   {
392  	      case SCIP_NLPSOLSTAT_GLOBOPT:
393  	      case SCIP_NLPSOLSTAT_LOCOPT:
394  	         /* @todo: if solution is optimal, we might as well add the cut <x - P(x_0), x_0 - P(x_0)> <= 0
395  	          * even though this cut is implied by all the gradient cuts of the rows active at the projection,
396  	          * we do not add them all (only the gradient cuts of constraints that violated the LP solution */
397  	      case SCIP_NLPSOLSTAT_FEASIBLE:
398  	         /* generate cuts for violated constraints (at sol) that are active or still violated at the projection, since
399  	          * a suboptimal solution or numerical issues could give a solution of the projection problem where constraints
400  	          * are not active; if the solution of the projection problem is in the interior of the region, we do nothing
401  	          */
402  	
403  	         /* get solution: build SCIP_SOL out of nlpi sol */
404  	         SCIP_CALL( SCIPgetNlpiSolution(scip, sepadata->nlpi, sepadata->nlpiprob, &nlpisol, NULL, NULL, NULL, NULL) );
405  	         assert(nlpisol != NULL);
406  	
407  	         SCIP_CALL( SCIPcreateSol(scip, &projection, NULL) );
408  	         for( i = 0; i < nlpinvars; i++ )
409  	         {
410  	            SCIP_VAR* var;
411  	
412  	            var = sepadata->nlpivars[i];
413  	            assert(SCIPhashmapExists(sepadata->var2nlpiidx, (void*)var) );
414  	
415  	            SCIP_CALL( SCIPsetSolVal(scip, projection, var,
416  	                     nlpisol[SCIPhashmapGetImageInt(sepadata->var2nlpiidx, (void *)var)]) );
417  	         }
418  	         SCIPdebug( SCIPprintSol(scip, projection, NULL, TRUE) );
419  	
420  	         /** @todo this could just be created inside generateCut and the extra argument removed */
421  	         SCIP_CALL( SCIPcreateExpriter(scip, &exprit) );
422  	
423  	         /* check for active or violated constraints */
424  	         for( i = 0; i < sepadata->nnlrows; ++i )
425  	         {
426  	            SCIP_NLROW* nlrow;
427  	            CONVEXSIDE convexside;
428  	            SCIP_Real activity;
429  	
430  	            /* ignore constraints that are not violated by `sol` */
431  	            if( SCIPisFeasZero(scip, sepadata->constraintviolation[i]) )
432  	               continue;
433  	
434  	            convexside = sepadata->convexsides[i];
435  	            nlrow = sepadata->nlrows[i];
436  	            assert(nlrow != NULL);
437  	
438  	            /* check for currently active constraints at projected point */
439  	            SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, projection, &activity) );
440  	
441  	            SCIPdebugMsg(scip, "NlRow activity at nlpi solution: %g <= %g <= %g\n", SCIPnlrowGetLhs(nlrow), activity,
442  	                  SCIPnlrowGetRhs(nlrow) );
443  	
444  	            /* if nlrow is active or violates the projection, build gradient cut at projection */
445  	            if( (convexside == RHS && SCIPisFeasGE(scip, activity, SCIPnlrowGetRhs(nlrow)))
446  	               || (convexside == LHS && SCIPisFeasLE(scip, activity, SCIPnlrowGetLhs(nlrow))) )
447  	            {
448  	               SCIP_ROW* row;
449  	
450  	               SCIP_CALL( generateCut(scip, sepa, projection, nlrow, convexside, activity, exprit,
451  	                        &row) );
452  	
453  	               SCIPdebugMsg(scip, "active or violated nlrow: (sols vio: %e)\n", sepadata->constraintviolation[i]);
454  	               SCIPdebug( SCIP_CALL( SCIPprintNlRow(scip, nlrow, NULL) ) );
455  	               SCIPdebugMsg(scip, "cut with efficacy %g generated\n", SCIPgetCutEfficacy(scip, sol, row));
456  	               SCIPdebug( SCIPprintRow(scip, row, NULL) );
457  	
458  	               /* add cut if it is efficacious for the point we want to separate (sol) */
459  	               if( SCIPisCutEfficacious(scip, sol, row) )
460  	               {
461  	                  SCIP_Bool infeasible;
462  	
463  	                  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
464  	
465  	                  if( infeasible )
466  	                  {
467  	                     *result = SCIP_CUTOFF;
468  	                     SCIP_CALL( SCIPreleaseRow(scip, &row) );
469  	                     break;
470  	                  }
471  	                  else
472  	                  {
473  	                     *result = SCIP_SEPARATED;
474  	                  }
475  	               }
476  	
477  	               /* release the row */
478  	               SCIP_CALL( SCIPreleaseRow(scip, &row) );
479  	            }
480  	         }
481  	
482  	         SCIPfreeExpriter(&exprit);
483  	
484  	#ifdef SCIP_DEBUG
485  	         {
486  	            SCIP_Real distance;
487  	
488  	            /* compute distance between LP sol and its projection (only makes sense when it is optimal) */
489  	            distance = 0.0;
490  	            for( i = 0; i < SCIPgetNNLPVars(scip); ++i )
491  	            {
492  	               SCIP_VAR* var;
493  	
494  	               var = SCIPgetNLPVars(scip)[i];
495  	               assert(var != NULL);
496  	
497  	               /* assert NLP solution is within the bounds of the variable (only make sense when sol is optimal) */
498  	               if( !SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)) )
499  	                  assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(var), SCIPvarGetNLPSol(var)));
500  	               if( !SCIPisInfinity(scip, SCIPvarGetUbLocal(var)) )
501  	                  assert(SCIPisFeasLE(scip, SCIPvarGetNLPSol(var), SCIPvarGetUbLocal(var)));
502  	
503  	               /*SCIPdebugMsg(scip, "NLP sol (LP sol): %s = %f (%g)\n", SCIPvarGetName(var),
504  	                *     SCIPvarGetNLPSol(var), SCIPgetSolVal(scip, sol, var));
505  	                */
506  	
507  	               distance += SQR( SCIPvarGetNLPSol(var) - SCIPgetSolVal(scip, sol, var) );
508  	            }
509  	
510  	            SCIPdebugMsg(scip, "NLP objval: %e, distance: %e\n", SCIPgetNLPObjval(scip), distance);
511  	         }
512  	#endif
513  	
514  	         /* free solution */
515  	         SCIP_CALL( SCIPfreeSol(scip, &projection) );
516  	         break;
517  	
518  	      case SCIP_NLPSOLSTAT_GLOBINFEASIBLE:
519  	      case SCIP_NLPSOLSTAT_LOCINFEASIBLE:
520  	         /* fallthrough;
521  	          * @todo: write what it means to be locinfeasible and why it can't be used to cutoff the node */
522  	      case SCIP_NLPSOLSTAT_UNKNOWN:
523  	         /* unknown... assume numerical issues */
524  	         nlpunstable = TRUE;
525  	         break;
526  	
527  	      case SCIP_NLPSOLSTAT_UNBOUNDED:
528  	      default:
529  	         SCIPerrorMessage("Projection NLP is not unbounded by construction, should not get here!\n");
530  	         SCIPABORT();
531  	         nlpunstable = TRUE;
532  	   }
533  	
534  	   /* if nlp is detected to be unstable, don't try to separate again */
535  	   if( nlpunstable )
536  	   {
537  	      /* @todo: maybe change objective function to \sum [(x_i - x_i^*)/max(|x_i^*|, 1)]^2
538  	       * or some other scaling when unstable and try again.
539  	       *       maybe free it here */
540  	      sepadata->skipsepa = TRUE;
541  	   }
542  	
543  	   /* reset objective */
544  	   BMSclearMemoryArray(linvals, nlpinvars);
545  	   SCIP_CALL( SCIPchgNlpiLinearCoefs(scip, sepadata->nlpi, sepadata->nlpiprob, -1, nlpinvars, lininds, linvals) );
546  	
547  	CLEANUP:
548  	   /* free memory */
549  	   SCIPfreeBufferArray(scip, &lininds);
550  	   SCIPfreeBufferArray(scip, &linvals);
551  	
552  	   return SCIP_OKAY;
553  	}
554  	
555  	/** computes the violation and maximum violation of the convex nlrows stored in sepadata wrt sol */
556  	static
557  	SCIP_RETCODE computeMaxViolation(
558  	   SCIP*                 scip,               /**< SCIP data structure */
559  	   SCIP_SEPADATA*        sepadata,           /**< separator data */
560  	   SCIP_SOL*             sol,                /**< solution that should be separated */
561  	   SCIP_Real*            maxviolation        /**< buffer to store maximum violation */
562  	   )
563  	{
564  	   SCIP_NLROW*    nlrow;
565  	   int            i;
566  	
567  	   assert(sepadata != NULL);
568  	   assert(sepadata->nnlrows > 0);
569  	   assert(sepadata->nlrows != NULL);
570  	   assert(sepadata->convexsides != NULL);
571  	   assert(sepadata->constraintviolation != NULL);
572  	
573  	   *maxviolation = 0.0;
574  	   for( i = 0; i < sepadata->nnlrows; i++ )
575  	   {
576  	      SCIP_Real activity;
577  	      SCIP_Real violation;
578  	
579  	      nlrow = sepadata->nlrows[i];
580  	
581  	      /* get activity of nlrow */
582  	      SCIP_CALL( SCIPgetNlRowSolActivity(scip, nlrow, sol, &activity) );
583  	
584  	      /* violation = max{activity - rhs, 0.0} when convex and max{lhs - activity, 0.0} when concave */
585  	      if( sepadata->convexsides[i] == RHS )
586  	      {
587  	         assert(SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONVEX);
588  	         assert(!SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)));
589  	
590  	         violation = activity - SCIPnlrowGetRhs(nlrow);
591  	         sepadata->constraintviolation[i] = MAX(violation, 0.0);
592  	      }
593  	      if( sepadata->convexsides[i] == LHS )
594  	      {
595  	         assert(SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONCAVE);
596  	         assert(!SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow)));
597  	
598  	         violation = SCIPnlrowGetLhs(nlrow) - activity;
599  	         sepadata->constraintviolation[i] = MAX(violation, 0.0);
600  	      }
601  	
602  	      /* compute maximum */
603  	      if( *maxviolation < sepadata->constraintviolation[i] )
604  	         *maxviolation = sepadata->constraintviolation[i];
605  	   }
606  	
607  	   SCIPdebugMsg(scip, "Maximum violation %g\n", *maxviolation);
608  	
609  	   return SCIP_OKAY;
610  	}
611  	
612  	
613  	/** stores, from the constraints represented by nlrows, the nonlinear convex ones in sepadata */
614  	static
615  	SCIP_RETCODE storeNonlinearConvexNlrows(
616  	   SCIP*                 scip,               /**< SCIP data structure */
617  	   SCIP_SEPADATA*        sepadata,           /**< separator data */
618  	   SCIP_NLROW**          nlrows,             /**< nlrows from which to store convex ones */
619  	   int                   nnlrows             /**< number of nlrows */
620  	   )
621  	{
622  	   int i;
623  	
624  	   assert(scip != NULL);
625  	   assert(sepadata != NULL);
626  	
627  	   SCIPdebugMsg(scip, "storing convex nlrows\n");
628  	
629  	   sepadata->nlrowssize = nnlrows;
630  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->nlrows), nnlrows) );
631  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->convexsides), nnlrows) );
632  	   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(sepadata->constraintviolation), nnlrows) );
633  	
634  	   /* count the number of nonlinear convex rows and store them */
635  	   sepadata->nnlrows = 0;
636  	   for( i = 0; i < nnlrows; ++i )
637  	   {
638  	      SCIP_NLROW* nlrow;
639  	
640  	      nlrow = nlrows[i];
641  	      assert(nlrow != NULL);
642  	
643  	      /* linear case */
644  	      if( SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_LINEAR || SCIPnlrowGetExpr(nlrow) == NULL )
645  	         continue;
646  	
647  	      /* nonlinear case */
648  	      if( !SCIPisInfinity(scip, SCIPnlrowGetRhs(nlrow)) && SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONVEX )
649  	      {
650  	         sepadata->convexsides[sepadata->nnlrows] = RHS;
651  	         sepadata->nlrows[sepadata->nnlrows] = nlrow;
652  	         ++(sepadata->nnlrows);
653  	      }
654  	      else if( !SCIPisInfinity(scip, -SCIPnlrowGetLhs(nlrow)) && SCIPnlrowGetCurvature(nlrow) == SCIP_EXPRCURV_CONCAVE )
655  	      {
656  	         sepadata->convexsides[sepadata->nnlrows] = LHS;
657  	         sepadata->nlrows[sepadata->nnlrows] = nlrow;
658  	         ++(sepadata->nnlrows);
659  	      }
660  	   }
661  	
662  	   return SCIP_OKAY;
663  	}
664  	
665  	
666  	/*
667  	 * Callback methods of separator
668  	 */
669  	
670  	
671  	/** destructor of separator to free user data (called when SCIP is exiting) */
672  	static
673  	SCIP_DECL_SEPAFREE(sepaFreeConvexproj)
674  	{  /*lint --e{715}*/
675  	   SCIP_SEPADATA* sepadata;
676  	
677  	   assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
678  	
679  	   /* free separator data */
680  	   sepadata = SCIPsepaGetData(sepa);
681  	   assert(sepadata != NULL);
682  	
683  	   SCIP_CALL( sepadataClear(scip, sepadata) );
684  	
685  	   SCIPfreeBlockMemory(scip, &sepadata);
686  	
687  	   SCIPsepaSetData(sepa, NULL);
688  	
689  	   return SCIP_OKAY;
690  	}
691  	
692  	/** solving process deinitialization method of separator (called before branch and bound process data is freed) */
693  	static
694  	SCIP_DECL_SEPAEXITSOL(sepaExitsolConvexproj)
695  	{  /*lint --e{715}*/
696  	   SCIP_SEPADATA* sepadata;
697  	
698  	   assert(sepa != NULL);
699  	
700  	   sepadata = SCIPsepaGetData(sepa);
701  	
702  	   assert(sepadata != NULL);
703  	
704  	   SCIP_CALL( sepadataClear(scip, sepadata) );
705  	
706  	   return SCIP_OKAY;
707  	}
708  	
709  	
710  	/** LP solution separation method of separator */
711  	static
712  	SCIP_DECL_SEPAEXECLP(sepaExeclpConvexproj)
713  	{  /*lint --e{715}*/
714  	   SCIP_Real maxviolation;
715  	   SCIP_SOL* lpsol;
716  	   SCIP_SEPADATA* sepadata;
717  	
718  	   *result = SCIP_DIDNOTRUN;
719  	
720  	   sepadata = SCIPsepaGetData(sepa);
721  	   assert(sepadata != NULL);
722  	
723  	   /* do not run if there is no interesting convex relaxation (with at least one nonlinear convex constraint),
724  	    * or if we have found it to be numerically unstable
725  	    * @todo: should it be with at least 2 nonlinear convex constraints?
726  	    */
727  	   if( sepadata->skipsepa )
728  	   {
729  	      SCIPdebugMsg(scip, "not running because convex relaxation is uninteresting or numerically unstable\n");
730  	      return SCIP_OKAY;
731  	   }
732  	
733  	   /* the separator needs an NLP solver */
734  	   if( SCIPgetNNlpis(scip) == 0 )
735  	      return SCIP_OKAY;
736  	
737  	   /* only call separator up to a maximum depth */
738  	   if( sepadata->maxdepth >= 0 && depth > sepadata->maxdepth )
739  	      return SCIP_OKAY;
740  	
741  	   /* only call separator, if we are not close to terminating */
742  	   if( SCIPisStopped(scip) )
743  	      return SCIP_OKAY;
744  	
745  	   /* do not run if SCIP does not have constructed an NLP */
746  	   if( !SCIPisNLPConstructed(scip) )
747  	   {
748  	      SCIPdebugMsg(scip, "NLP not constructed, skipping convex projection separator\n");
749  	      return SCIP_OKAY;
750  	   }
751  	
752  	   /* recompute convex NLP relaxation if the variable set changed and we are still at the root node */
753  	   if( sepadata->nlpiprob != NULL && SCIPgetNVars(scip) != sepadata->nlpinvars  && SCIPgetDepth(scip) == 0 )
754  	   {
755  	      SCIP_CALL( sepadataClear(scip, sepadata) );
756  	      assert(sepadata->nlpiprob == NULL);
757  	   }
758  	
759  	   /* create or update convex NLP relaxation */
760  	   if( sepadata->nlpiprob == NULL )
761  	   {
762  	      /* store convex nonlinear constraints */
763  	      SCIP_CALL( storeNonlinearConvexNlrows(scip, sepadata, SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip)) );
764  	
765  	      /* check that convex NLP relaxation is interesting (more than one nonlinear constraint) */
766  	      if( sepadata->nnlrows < 1 )
767  	      {
768  	         SCIPdebugMsg(scip, "convex relaxation uninteresting, don't run\n");
769  	         sepadata->skipsepa = TRUE;
770  	         return SCIP_OKAY;
771  	      }
772  	
773  	      sepadata->nlpinvars = SCIPgetNVars(scip);
774  	      sepadata->nlpi = SCIPgetNlpis(scip)[0];
775  	      assert(sepadata->nlpi != NULL);
776  	
777  	      SCIP_CALL( SCIPhashmapCreate(&sepadata->var2nlpiidx, SCIPblkmem(scip), sepadata->nlpinvars) );
778  	      SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &sepadata->nlpivars, SCIPgetVars(scip), sepadata->nlpinvars) ); /*lint !e666*/
779  	
780  	      SCIP_CALL( SCIPcreateNlpiProblemFromNlRows(scip, sepadata->nlpi, &sepadata->nlpiprob, "convexproj-nlp", SCIPgetNLPNlRows(scip), SCIPgetNNLPNlRows(scip),
781  	            sepadata->var2nlpiidx, NULL, NULL, SCIPgetCutoffbound(scip), FALSE, TRUE) );
782  	
783  	      /* add rows of the LP
784  	       * we do not sue the depth argument of the callback because we want to build a globally valid initia lrelaxation
785  	       */
786  	      if( SCIPgetDepth(scip) == 0 )
787  	      {
788  	         SCIP_CALL( SCIPaddNlpiProblemRows(scip, sepadata->nlpi, sepadata->nlpiprob, sepadata->var2nlpiidx,
789  	                  SCIPgetLPRows(scip), SCIPgetNLPRows(scip)) );
790  	      }
791  	
792  	      /* set quadratic part of objective function */
793  	      SCIP_CALL( setQuadraticObj(scip, sepadata) );
794  	   }
795  	   else
796  	   {
797  	      SCIP_CALL( SCIPupdateNlpiProblem(scip, sepadata->nlpi, sepadata->nlpiprob, sepadata->var2nlpiidx,
798  	            sepadata->nlpivars, sepadata->nlpinvars, SCIPgetCutoffbound(scip)) );
799  	   }
800  	
801  	   /* assert that the lp solution satisfies the cutoff bound; if this fails then we shouldn't have a cutoff bound in the
802  	    * nlpi, since then the projection could be in the interior of the actual convex relaxation */
803  	   assert(SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ||
804  	         SCIPisLE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)));
805  	
806  	   /* get current sol: LP or pseudo solution if LP sol is not available */
807  	   SCIP_CALL( SCIPcreateCurrentSol(scip, &lpsol, NULL) );
808  	
809  	   /* do not run if current solution's violation is small */
810  	   SCIP_CALL( computeMaxViolation(scip, sepadata, lpsol, &maxviolation) );
811  	   if( maxviolation < VIOLATIONFAC * SCIPfeastol(scip) )
812  	   {
813  	      SCIPdebugMsg(scip, "solution doesn't violate constraints enough, do not separate\n");
814  	      SCIP_CALL( SCIPfreeSol(scip, &lpsol) );
815  	      return SCIP_OKAY;
816  	   }
817  	
818  	   /* run the separator */
819  	   *result = SCIP_DIDNOTFIND;
820  	
821  	   /* separateCuts computes the projection and then gradient cuts on each constraint that was originally violated */
822  	   SCIP_CALL( separateCuts(scip, sepa, lpsol, result) );
823  	
824  	   /* free memory */
825  	   SCIP_CALL( SCIPfreeSol(scip, &lpsol) );
826  	
827  	   return SCIP_OKAY;
828  	}
829  	
830  	
831  	/*
832  	 * separator specific interface methods
833  	 */
834  	
835  	/** creates the convexproj separator and includes it in SCIP */
836  	SCIP_RETCODE SCIPincludeSepaConvexproj(
837  	   SCIP*                 scip                /**< SCIP data structure */
838  	   )
839  	{
840  	   SCIP_SEPADATA* sepadata;
841  	   SCIP_SEPA* sepa;
842  	
843  	   /* create convexproj separator data */
844  	   SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
845  	
846  	   /* this sets all data in sepadata to 0 */
847  	   BMSclearMemory(sepadata);
848  	
849  	   /* include separator */
850  	   SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
851  	         SEPA_USESSUBSCIP, SEPA_DELAY,
852  	         sepaExeclpConvexproj, NULL,
853  	         sepadata) );
854  	   assert(sepa != NULL);
855  	
856  	   /* set non fundamental callbacks via setter functions */
857  	   SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeConvexproj) );
858  	   SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolConvexproj) );
859  	
860  	   /* add convexproj separator parameters */
861  	   SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxdepth",
862  	         "maximal depth at which the separator is applied (-1: unlimited)",
863  	         &sepadata->maxdepth, FALSE, DEFAULT_MAXDEPTH, -1, INT_MAX, NULL, NULL) );
864  	
865  	   SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nlpiterlimit",
866  	         "iteration limit of NLP solver; 0 for no limit",
867  	         &sepadata->nlpiterlimit, TRUE, DEFAULT_NLPITERLIM, 0, INT_MAX, NULL, NULL) );
868  	
869  	   return SCIP_OKAY;
870  	}
871