Scippy

SCIP

Solving Constraint Integer Programs

sepa_lagromory.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (c) 2002-2024 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_lagromory.c
26  * @ingroup DEFPLUGINS_SEPA
27  * @brief Lagromory separator
28  * @author Suresh Bolusani
29  * @author Mark Ruben Turner
30  * @author Mathieu Besançon
31  *
32  * This separator is based on the following article that discusses Lagromory separation using the relax-and-cut
33  * framework. Multiple enhancements have been implemented on top of the basic algorithm described in the article.
34  *
35  * Fischetti M. and Salvagnin D. (2011).@n
36  * A relax-and-cut framework for Gomory mixed-integer cuts.@n
37  * Mathematical Programming Computation, 3, 79-102.
38  *
39  * Consider the following linear relaxation at a node:
40  *
41  * \f[
42  * \begin{array}{rrl}
43  * \min & c^T x &\\
44  * & x & \in P,
45  * \end{array}
46  * \f]
47  *
48  * where \f$P\f$ is the feasible region of the relaxation. Let the following be the cuts generated so far in the current
49  * separation round.
50  *
51  * \f[
52  * {\alpha^i}^T x \leq \alpha^i_0, i = 1, 2, \hdots, M
53  * \f]
54  *
55  * Then, the following is the Lagrangian dual problem considered in the relax-and-cut framework used in the separator.
56  *
57  * \f[
58  * z_D := \max\limits_{u \geq 0} \left\{L(u) := \min \left\{c^T x + \sum\limits_{i = 1}^{M} \left(u_i
59  * \left({\alpha^i}^T x - \alpha^i_0\right) \right) \mid x \in P\right\} \right\},
60  * \f]
61  *
62  * where \f$u\f$ are the Lagrangian multipliers (referred to as \a dualvector in this separator) used for penalizing the
63  * violation of the generated cuts, and \f$z_D\f$ is the optimal objective value (which is approximated via \a ubparam in this separator).
64  * Then, the following are the steps of the relax-and-cut algorithm implemented in this separator.
65  *
66  * \begin{itemize}
67  * \item Generate an initial pool of cuts to build the initial Lagrangian dual problem.
68  * \item Select initial values for Lagrangian multipliers \f$u^0\f$ (e.g., all zeroes vector).
69  * \item In the outer main loop \f$i\f$ of the algorithm:
70  * \begin{enumerate}
71  * \item Solve the Lagrangian dual problem until certain termination criterion is met. This results in an inner
72  * subgradient loop, whose iteration \f$j\f$ is described below.
73  * \begin{enumerate}
74  * \item Fix \f$u^j\f$, and solve the LP corresponding to the Lagrangian dual with fixed multipliers.
75  * Gather its optimal simplex tableau and optimal objective value (i.e., the Lagrangian value)
76  * \f$L(u^j)\f$.
77  * \item Update \f$u^j\f$ to \f$u^{j+1}\f$ as follows.
78  * \f[
79  * u^{j+1} = \left(u^j + \lambda_j s^k\right)_+,
80  * \f]
81  * where \f$\lambda_j\f$ is the step length:
82  * \f[
83  * \lambda_j = \frac{\mu_j (UB - L(u^j))}{\|s^j\|^2_2},
84  * \f]
85  * where \f$mu_j\f$ is a factor (i.e., \a muparam) such that \f$0 < \mu_j \leq 2\f$, UB is \p ubparam,
86  * and \f$s^j\f$ is the subgradient vector defined as:
87  * \f[
88  * s^j_k = \left({\alpha^k}^T x - \alpha^k_0\right), k = 1, 2, \hdots, M.
89  * \f]
90  * The factor \f$mu_j\f$ is updated as below.
91  * \f[
92  * mu_j = \begin{cases}
93  * \p mubacktrackfactor * mu_j & \text{if } L(u^j) < bestLB - \delta\\
94  * \begin{cases}
95  * \p muslab1factor * mu_j & \text{if } bestLB - avgLB < \p deltaslab1ub * delta\\
96  * \p muslab2factor * mu_j & \text{if } \p deltaslab1ub * \delta \leq bestLB - avgLB < \p deltaslab2ub * delta\\
97  * \p muslab3factor * mu_j & \text{otherwise}
98  * \end{cases} & \text{otherwise},
99  * \end{cases}
100  * \f]
101  * where \f$bestLB\f$ and \f$avgLB\f$ are best and average Lagrangian values found so far, and
102  * \f$\delta = UB - bestLB\f$.
103  * \item Stabilize \f$u^{j+1}\f$ by projecting onto a norm ball followed by taking a convex combination
104  * with a core vector of Lagrangian multipliers.
105  * \item Generate GMI cuts based on the optimal simplex tableau.
106  * \item Relax the newly generated cuts by penalizing and adding them to the objective function.
107  * \item Go to the next iteration \f$j+1\f$.
108  * \end{enumerate}
109  * \item Gather all the generated cuts and build an LP by adding all these cuts to the node relaxation.
110  * \item Solve this LP to obtain its optimal primal and dual solutions.
111  * \item If this primal solution is MIP primal feasible, then add this solution to the solution pool, add all
112  * the generated cuts to the cutpool or sepastore as needed, and exit the separator.
113  * \item Otherwise, update the Lagrangian multipliers based on this optimal dual solution, and go to the next
114  * iteration \f$i+1\f$.
115  * \end{enumerate}
116  * \end{itemize}
117  *
118  * @todo store all LP sols in a data structure, and send them to fix-and-propagate at the end.
119  *
120  * @todo test heuristics such as feasibility pump with multiple input solutions.
121  *
122  * @todo find dual degenerate problems and test the separator on these problems.
123  *
124  * @todo identify instance classes where these cuts work better.
125  *
126  * @todo add termination criteria based on failed efforts.
127  *
128  * @todo for warm starting, if there are additional rows/cols, set their basis status to non-basic and then set WS info.
129  *
130  * @todo filter cuts using multiple explored LP solutions.
131  *
132  * @todo for bases on optimal face only, aggregate to get a new basis and separate it.
133  *
134  * @todo generate other separators in addition to GMI cuts (0-1/2)
135  *
136  * @todo: convert iters from int to SCIP_Longint
137  */
138 
139 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
140 
141 #include "scip/heur_trysol.h"
142 #include "scip/sepa_lagromory.h"
143 
144 #define SEPA_NAME "lagromory"
145 #define SEPA_DESC "separator for Lagromory cuts for MIP relaxations"
146 #define SEPA_PRIORITY -8000
147 #define SEPA_FREQ -1
148 #define SEPA_MAXBOUNDDIST 1.0
149 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */
150 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */
151 
152 /* generic parameters */
153 #define DEFAULT_AWAY 0.01 /**< minimal integrality violation of a basis variable to try separation */
154 #define DEFAULT_DELAYEDCUTS FALSE /**< should cuts be added to the delayed cut pool? */
155 #define DEFAULT_SEPARATEROWS TRUE /**< separate rows with integral slack? */
156 #define DEFAULT_SORTCUTOFFSOL TRUE /**< sort fractional integer columns based on fractionality? */
157 #define DEFAULT_SIDETYPEBASIS TRUE /**< choose side types of row (lhs/rhs) based on basis information? */
158 #define DEFAULT_DYNAMICCUTS TRUE /**< should generated cuts be removed from the LP if they are no longer tight? */
159 #define DEFAULT_MAKEINTEGRAL FALSE /**< try to scale all cuts to integral coefficients? */
160 #define DEFAULT_FORCECUTS FALSE /**< force cuts to be added to the LP? */
161 #define DEFAULT_ALLOWLOCAL FALSE /**< should locally valid cuts be generated? */
162 
163 /* parameters related to the separator's termination check */
164 #define DEFAULT_MAXROUNDSROOT 1 /**< maximal number of separation rounds in the root node (-1: unlimited) */
165 #define DEFAULT_MAXROUNDS 1 /**< maximal number of separation rounds per node (-1: unlimited) */
166 #define DEFAULT_DUALDEGENERACYRATETHRESHOLD 0.5 /**< minimum dual degeneracy rate for separator execution */
167 #define DEFAULT_VARCONSRATIOTHRESHOLD 1.0 /**< minimum variable-constraint ratio on optimal face for separator execution */
168 #define DEFAULT_MINRESTART 1 /**< minimum restart round for separator execution (0: from beginning of the
169  instance solving, >= n with n >= 1: from restart round n) */
170 #define DEFAULT_PERLPMAXCUTSROOT 50 /**< maximal number of cuts separated per Lagromory LP in the root node */
171 #define DEFAULT_PERLPMAXCUTS 10 /**< maximal number of cuts separated per Lagromory LP in the non-root node */
172 #define DEFAULT_PERROUNDLPITERLIMITFACTOR -1.0 /**< factor w.r.t. root node LP iterations for maximal separating LP iterations
173  per separation round (negative for no limit) */
174 #define DEFAULT_ROOTLPITERLIMITFACTOR -1.0 /**< factor w.r.t. root node LP iterations for maximal separating LP iterations
175  in the root node (negative for no limit) */
176 #define DEFAULT_TOTALLPITERLIMITFACTOR -1.0 /**< factor w.r.t. root node LP iterations for maximal separating LP iterations
177  in the tree (negative for no limit) */
178 #define DEFAULT_PERROUNDMAXLPITERS 50000 /**< maximal number of separating LP iterations per separation round (-1: unlimited) */
179 #define DEFAULT_PERROUNDCUTSFACTORROOT 1.0 /**< factor w.r.t. number of integer columns for number of cuts separated per
180  separation round in root node */
181 #define DEFAULT_PERROUNDCUTSFACTOR 0.5 /**< factor w.r.t. number of integer columns for number of cuts separated per
182  separation round at a non-root node */
183 #define DEFAULT_TOTALCUTSFACTOR 50.0 /**< factor w.r.t. number of integer columns for total number of cuts separated */
184 #define DEFAULT_MAXMAINITERS 4 /**< maximal number of main loop iterations of the relax-and-cut algorithm */
185 #define DEFAULT_MAXSUBGRADIENTITERS 6 /**< maximal number of subgradient loop iterations of the relax-and-cut algorithm */
187 /* parameters related to the relax-and-cut algorithm */
188 #define DEFAULT_MUPARAMCONST TRUE /**< is the mu parameter (factor for step length) constant? */
189 #define DEFAULT_MUPARAMINIT 0.01 /**< initial value of the mu parameter (factor for step length) */
190 #define DEFAULT_MUPARAMLB 0.0 /**< lower bound for the mu parameter (factor for step length) */
191 #define DEFAULT_MUPARAMUB 2.0 /**< upper bound for the mu parameter (factor for step length) */
192 #define DEFAULT_MUBACKTRACKFACTOR 0.5 /**< factor of mu while backtracking the mu parameter (factor for step length) -
193  see updateMuSteplengthParam() */
194 #define DEFAULT_MUSLAB1FACTOR 10.0 /**< factor of mu parameter (factor for step length) for larger increment - see
195  updateMuSteplengthParam() */
196 #define DEFAULT_MUSLAB2FACTOR 2.0 /**< factor of mu parameter (factor for step length) for smaller increment - see
197  updateMuSteplengthParam() */
198 #define DEFAULT_MUSLAB3FACTOR 0.5 /**< factor of mu parameter (factor for step length) for reduction - see
199  updateMuSteplengthParam() */
200 #define DEFAULT_DELTASLAB1UB 0.001 /**< factor of delta deciding larger increment of mu parameter (factor for step
201  length) - see updateMuSteplengthParam() */
202 #define DEFAULT_DELTASLAB2UB 0.01 /**< factor of delta deciding smaller increment of mu parameter (factor for step
203  length) - see updateMuSteplengthParam() */
204 #define DEFAULT_UBPARAMPOSFACTOR 2.0 /**< factor for positive upper bound used as an estimate for the optimal
205  Lagrangian dual value */
206 #define DEFAULT_UBPARAMNEGFACTOR 0.5 /**< factor for negative upper bound used as an estimate for the optimal
207  Lagrangian dual value */
208 #define DEFAULT_MAXLAGRANGIANVALSFORAVG 2 /**< maximal number of iterations for rolling average of Lagrangian value */
209 #define DEFAULT_MAXCONSECITERSFORMUUPDATE 10 /**< consecutive number of iterations used to determine if mu needs to be backtracked */
210 #define DEFAULT_PERROOTLPITERFACTOR 0.2 /**< factor w.r.t. root node LP iterations for iteration limit of each separating
211  LP (negative for no limit) */
212 #define DEFAULT_PERLPITERFACTOR 0.1 /**< factor w.r.t. non-root node LP iterations for iteration limit of each
213  separating LP (negative for no limit) */
214 #define DEFAULT_CUTGENFREQ 1 /**< frequency of subgradient iterations for generating cuts */
215 #define DEFAULT_CUTADDFREQ 1 /**< frequency of subgradient iterations for adding cuts to objective function */
216 #define DEFAULT_CUTSFILTERFACTOR 1.0 /**< fraction of generated cuts per explored basis to accept from separator */
217 #define DEFAULT_OPTIMALFACEPRIORITY 2 /**< priority of the optimal face for separator execution (0: low priority, 1:
218  medium priority, 2: high priority) */
219 #define DEFAULT_AGGREGATECUTS TRUE /**< aggregate all generated cuts using the Lagrangian multipliers? */
220 /* parameters for stabilization of the Lagrangian multipliers */
221 #define DEFAULT_PROJECTIONTYPE 2 /**< the ball into which the Lagrangian multipliers are projected for
222  stabilization (0: no projection, 1: L1-norm ball projection, 2: L2-norm ball
223  projection, 3: L_inf-norm ball projection) */
224 #define DEFAULT_STABILITYCENTERTYPE 1 /**< type of stability center for taking weighted average of Lagrangian multipliers for
225  stabilization (0: no weighted stabilization, 1: best Lagrangian multipliers) */
226 #define DEFAULT_RADIUSINIT 0.5 /**< initial radius of the ball used in stabilization of Lagrangian multipliers */
227 #define DEFAULT_RADIUSMAX 20.0 /**< maximum radius of the ball used in stabilization of Lagrangian multipliers */
228 #define DEFAULT_RADIUSMIN 1e-6 /**< minimum radius of the ball used in stabilization of Lagrangian multipliers */
229 #define DEFAULT_CONST 2.0 /**< a constant for stablity center based stabilization of Lagrangian multipliers */
230 #define DEFAULT_RADIUSUPDATEWEIGHT 0.98 /**< multiplier to evaluate cut violation score used for updating ball radius */
232 /* macros that are used directly */
233 #define RANDSEED 42 /**< random seed */
234 #define MAKECONTINTEGRAL FALSE /**< convert continuous variable to integral variables in SCIPmakeRowIntegral()? */
235 #define POSTPROCESS TRUE /**< apply postprocessing after MIR calculation? - see SCIPcalcMIR() */
236 #define BOUNDSWITCH 0.9999 /**< threshold for bound switching - see SCIPcalcMIR() */
237 #define USEVBDS TRUE /**< use variable bounds? - see SCIPcalcMIR() */
238 #define FIXINTEGRALRHS FALSE /**< try to generate an integral rhs? - see SCIPcalcMIR() */
239 #define MAXAGGRLEN(ncols) (0.1*(ncols)+1000) /**< maximal length of base inequality */
240 
241 /*
242  * Data structures
243  */
244 
245 /** separator data */
246 struct SCIP_SepaData
247 {
248  /* generic variables */
249  SCIP_Real away; /**< minimal integrality violation of a basis variable to try separation */
250  SCIP_Bool delayedcuts; /**< should cuts be added to the delayed cut pool? */
251  SCIP_Bool separaterows; /**< separate rows with integral slack? */
252  SCIP_Bool sortcutoffsol; /**< sort fractional integer columns based on fractionality? */
253  SCIP_Bool sidetypebasis; /**< choose side types of row (lhs/rhs) based on basis information? */
254  SCIP_Bool dynamiccuts; /**< should generated cuts be removed from the LP if they are no longer tight? */
255  SCIP_Bool makeintegral; /**< try to scale all cuts to integral coefficients? */
256  SCIP_Bool forcecuts; /**< force cuts to be added to the LP? */
257  SCIP_Bool allowlocal; /**< should locally valid cuts be generated? */
258  SCIP_RANDNUMGEN* randnumgen; /**< random number generator */
259  SCIP_HEUR* heurtrysol; /**< a pointer to the trysol heuristic, if available */
260 
261  /* used to define separating LPs */
262  SCIP_LPI* lpiwithsoftcuts; /**< pointer to the lpi interface of Lagrangian dual with fixed multipliers */
263  SCIP_LPI* lpiwithhardcuts; /**< pointer to the lpi interface of LP with all generated cuts */
264  int nrowsinhardcutslp; /**< nrows of \a lpiwithhardcuts */
265  int nrunsinsoftcutslp; /**< number of branch-and-bound runs on current instance */
266 
267  /* used for termination checks */
268  SCIP_Longint ncalls; /**< number of calls to the separator */
269  int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */
270  int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */
271  SCIP_Real dualdegeneracyratethreshold; /**< minimum dual degeneracy rate for separator execution */
272  SCIP_Real varconsratiothreshold; /**< minimum variable-constraint ratio on optimal face for separator execution */
273  int minrestart; /**< minimum restart round for separator execution (0: from beginning of the instance
274  solving, >= n with n >= 1: from restart round n) */
275  int nmaxcutsperlproot; /**< maximal number of cuts separated per Lagromory LP in the root node */
276  int nmaxcutsperlp; /**< maximal number of cuts separated per Lagromory LP in the non-root node */
277  SCIP_Real perroundlpiterlimitfactor; /**< factor w.r.t. root node LP iterations for maximal separating LP iterations per
278  separation round (negative for no limit) */
279  SCIP_Real rootlpiterlimitfactor; /**< factor w.r.t. root node LP iterations for maximal separating LP iterations in the
280  root node (negative for no limit) */
281  SCIP_Real totallpiterlimitfactor; /**< factor w.r.t. root node LP iterations for maximal separating LP iterations in the
282  tree (negative for no limit) */
283  int perroundnmaxlpiters; /**< maximal number of separating LP iterations per separation round (-1: unlimited) */
284  SCIP_Real perroundcutsfactorroot; /**< factor w.r.t. number of integer columns for number of cuts separated per
285  separation round in root node */
286  SCIP_Real perroundcutsfactor; /**< factor w.r.t. number of integer columns for number of cuts separated per
287  separation round at a non-root node */
288  SCIP_Real totalcutsfactor; /**< factor w.r.t. number of integer columns for total number of cuts separated */
289  int nmaxmainiters; /**< maximal number of main loop iterations of the relax-and-cut algorithm */
290  int nmaxsubgradientiters; /**< maximal number of subgradient loop iterations of the relax-and-cut algorithm */
291  int nmaxperroundlpiters; /**< maximal number of separating LP iterations per separation round */
292  int nmaxrootlpiters; /**< maximal number of separating LP iterations in the root node */
293  int nrootlpiters; /**< number of separating LP iterations in the root node */
294  int nmaxtotallpiters; /**< maximal number of separating LP iterations in the tree */
295  int ntotallpiters; /**< number of separating LP iterations in the tree */
296  int nmaxperroundcutsroot; /**< maximal number of cuts separated per separation round in root node */
297  int nmaxperroundcuts; /**< maximal number of cuts separated per separation round */
298  int nmaxtotalcuts; /**< maximal number of cuts separated in the tree */
299  int ntotalcuts; /**< number of cuts separated in the tree */
300 
301  /* used for the relax-and-cut algorithm */
302  SCIP_Bool muparamconst; /**< is the mu parameter (factor for step length) constant? */
303  SCIP_Real muparaminit; /**< initial value of the mu parameter (factor for step length) */
304  SCIP_Real muparamlb; /**< lower bound for the mu parameter (factor for step length) */
305  SCIP_Real muparamub; /**< upper bound for the mu parameter (factor for step length) */
306  SCIP_Real mubacktrackfactor; /**< factor of mu while backtracking the mu parameter (factor for step length) - see
307  updateMuSteplengthParam() */
308  SCIP_Real muslab1factor; /**< factor of mu parameter (factor for step length) for larger increment - see
309  updateMuSteplengthParam() */
310  SCIP_Real muslab2factor; /**< factor of mu parameter (factor for step length) for smaller increment - see
311  updateMuSteplengthParam() */
312  SCIP_Real muslab3factor; /**< factor of mu parameter (factor for step length) for reduction - see updateMuSteplengthParam() */
313  SCIP_Real deltaslab1ub; /**< factor of delta deciding larger increment of mu parameter (factor for step
314  length) - see updateMuSteplengthParam() */
315  SCIP_Real deltaslab2ub; /**< factor of delta deciding smaller increment of mu parameter (factor for step
316  length) - see updateMuSteplengthParam() */
317  SCIP_Real ubparamposfactor; /**< factor for positive upper bound used as an estimate for the optimal Lagrangian
318  dual value */
319  SCIP_Real ubparamnegfactor; /**< factor for negative upper bound used as an estimate for the optimal Lagrangian
320  dual value */
321  int nmaxlagrangianvalsforavg; /**< maximal number of iterations for rolling average of Lagrangian value */
322  int nmaxconsecitersformuupdate; /**< consecutive number of iterations used to determine if mu needs to be backtracked */
323  SCIP_Real perrootlpiterfactor; /**< factor w.r.t. root node LP iterations for iteration limit of each separating LP
324  (negative for no limit) */
325  SCIP_Real perlpiterfactor; /**< factor w.r.t. non-root node LP iterations for iteration limit of each separating
326  LP (negative for no limit) */
327  int cutgenfreq; /**< frequency of subgradient iterations for generating cuts */
328  int cutaddfreq; /**< frequency of subgradient iterations for adding cuts to objective function */
329  SCIP_Real cutsfilterfactor; /**< fraction of generated cuts per explored basis to accept from separator */
330  int optimalfacepriority; /**< priority of the optimal face for separator execution (0: low priority, 1: medium
331  priority, 2: high priority) */
332  SCIP_Bool aggregatecuts; /**< aggregate all generated cuts using the Lagrangian multipliers? */
333 
334  /* for stabilization of Lagrangian multipliers */
335  int projectiontype; /**< the ball into which the Lagrangian multipliers are projected for stabilization
336  (0: no projection, 1: L1-norm ball projection, 2: L2-norm ball projection, 3:
337  L_inf-norm ball projection) */
338  int stabilitycentertype; /**< type of stability center for taking weighted average of Lagrangian multipliers for
339  stabilization (0: no weighted stabilization, 1: best Lagrangian multipliers) */
340  SCIP_Real radiusinit; /**< initial radius of the ball used in stabilization of Lagrangian multipliers */
341  SCIP_Real radiusmax; /**< maximum radius of the ball used in stabilization of Lagrangian multipliers */
342  SCIP_Real radiusmin; /**< minimum radius of the ball used in stabilization of Lagrangian multipliers */
343  SCIP_Real constant; /**< a constant for stablity center based stabilization of Lagrangian multipliers */
344  SCIP_Real radiusupdateweight; /**< multiplier to evaluate cut violation score used for updating ball radius */
345 };
346 
347 
348 /*
349  * Local methods
350  */
351 
352 /** start the diving mode for solving LPs corresponding to the Lagrangian dual with fixed multipliers in the subgradient
353  * loop of the separator, and update some sepadata values */
354 static
356  SCIP* scip, /**< SCIP data structure */
357  SCIP_SEPADATA* sepadata /**< separator data structure */
358  )
359 {
360  SCIP_ROW** rows;
361  SCIP_COL** cols;
362  int nruns;
363  int runnum;
364  int nrows;
365  int ncols;
366  unsigned int nintcols;
367  SCIP_Longint nrootlpiters;
368 
369  assert(scip != NULL);
370  assert(sepadata != NULL);
371 
372  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
373  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
374  runnum = SCIPgetNRuns(scip); /* current run number of SCIP (starts from 1) indicating restarts */
375  nruns = sepadata->nrunsinsoftcutslp; /* previous run number of SCIP in which the diving LP was created */
376  nintcols = 0;
377 
378  /* start diving mode so that the diving LP can be used for all subgradient iterations */
379  SCIP_CALL( SCIPstartDive(scip) );
380 
381  /* store LPI pointer to be able to use LPI functions directly (e.g., setting time limit) */
382  SCIP_CALL( SCIPgetLPI(scip, &(sepadata->lpiwithsoftcuts)) );
383 
384  /* if called for the first time in a restart (including the initial run), set certain sepadata values */
385  if( nruns != runnum )
386  {
387  /* get number of LP iterations of root node's first LP solving */
388  nrootlpiters = SCIPgetNRootFirstLPIterations(scip);
389 
390  /* calculate maximum number of LP iterations allowed for all separation calls in the root node */
391  if( (sepadata->rootlpiterlimitfactor >= 0.0) && !SCIPisInfinity(scip, sepadata->rootlpiterlimitfactor) )
392  {
393  sepadata->nmaxrootlpiters = (int)(sepadata->rootlpiterlimitfactor * nrootlpiters);
394  }
395  else
396  {
397  sepadata->nmaxrootlpiters = -1; /* no finite limit */
398  }
399 
400  /* calculate maximum number of LP iterations allowed for all separation calls in the entire tree */
401  if( (sepadata->totallpiterlimitfactor >= 0.0) && !SCIPisInfinity(scip, sepadata->totallpiterlimitfactor) )
402  {
403  sepadata->nmaxtotallpiters = (int)(sepadata->totallpiterlimitfactor * nrootlpiters);
404  }
405  else
406  {
407  sepadata->nmaxtotallpiters = -1; /* no finite limit */
408  }
409 
410  /* calculate maximum number of LP iterations allowed per separation call */
411  if( (sepadata->perroundlpiterlimitfactor >= 0.0) && !SCIPisInfinity(scip, sepadata->perroundlpiterlimitfactor) )
412  {
413  sepadata->nmaxperroundlpiters = (int)(sepadata->perroundlpiterlimitfactor * nrootlpiters);
414  }
415  else
416  {
417  sepadata->nmaxperroundlpiters = -1; /* no finite limit */
418  }
419 
420  /* update maximum number of LP iterations allowed per separation call using absolute limits */
421  if( sepadata->perroundnmaxlpiters > 0 )
422  {
423  sepadata->nmaxperroundlpiters = ((sepadata->nmaxperroundlpiters >= 0) ? MIN(sepadata->nmaxperroundlpiters,
424  sepadata->perroundnmaxlpiters) : sepadata->perroundnmaxlpiters);
425  }
426 
427  /* set maximum number of cuts allowed to generate per round in root and non-root nodes as well as the total tree */
428  for( int i = 0; i < ncols; ++i )
429  {
430  nintcols += SCIPcolIsIntegral(cols[i]);
431  }
432  sepadata->nmaxperroundcutsroot = (int)(sepadata->perroundcutsfactorroot * nintcols);
433  sepadata->nmaxperroundcuts = (int)(sepadata->perroundcutsfactor * nintcols);
434 
435  if( sepadata->ncalls == 0 )
436  {
437  sepadata->nmaxtotalcuts = (int)(sepadata->totalcutsfactor * nintcols);
438  sepadata->ntotalcuts = 0;
439  }
440 
441  /* update the run number of solving to represent the restart number in which the above limits were set */
442  sepadata->nrunsinsoftcutslp = runnum;
443  }
444 
445  return SCIP_OKAY;
446 }
447 
448 /** end the diving mode that was used for solving LPs corresponding to the Lagrangian dual with fixed multipliers */
449 static
451  SCIP* scip, /**< SCIP data structure */
452  SCIP_SEPADATA* sepadata /**< separator data structure */
453  )
454 {
455  assert(scip != NULL);
456  assert(sepadata != NULL);
457 
458  SCIP_CALL( SCIPendDive(scip) );
459 
460  return SCIP_OKAY;
461 }
462 
463 /** set up LP interface to solve LPs in the (outer) main loop of the relax-and-cut algorithm; these LPs are built by
464  * adding all the generated cuts to the node relaxation */
465 /* @todo add lpi iters to global statistics */
466 static
468  SCIP* scip, /**< SCIP data structure */
469  SCIP_SEPADATA* sepadata, /**< separator data structure */
470  SCIP_ROW** cuts, /**< generated cuts to be added to the LP */
471  int ncuts /**< number of generated cuts to be added to the LP */
472  )
473 {
474  SCIP_ROW** rows;
475  SCIP_COL** cols;
476  SCIP_Real* collb;
477  SCIP_Real* colub;
478  SCIP_Real* colobj;
479  SCIP_Real* rowlhs;
480  SCIP_Real* rowrhs;
481  SCIP_Real rowconst;
482  SCIP_Real* rowvals;
483  SCIP_Real* rowval;
484  SCIP_COL** rowcols;
485  SCIP_LPI* wslpi;
486  SCIP_LPISTATE* lpistate;
487  BMS_BLKMEM* blkmem;
488  SCIP_Real pinf;
489  SCIP_Real ninf;
490  int nrows;
491  int ncols;
492  int nrownonz;
493  int collppos;
494  int* rowcolinds;
495  int* rowbegs;
496 
497  assert(scip != NULL);
498  assert(sepadata != NULL);
499 
500  SCIP_LPI** lpi = &(sepadata->lpiwithhardcuts);
501  blkmem = SCIPblkmem(scip);
502 
503  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
504  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
505 
506  /* if this function is called for the first time in this separation round, then create an LPI and add cols & rows */
507  if( ncuts == 0 )
508  {
509  if( *lpi != NULL )
510  {
511  SCIP_CALL( SCIPlpiFree(lpi) );
512  *lpi = NULL;
513  }
514  assert(*lpi == NULL);
515 
516  /* create an LPI with appropriate objective sense */
518  {
519  SCIP_CALL( SCIPlpiCreate(lpi, SCIPgetMessagehdlr(scip), "node LP with generated cuts", SCIP_OBJSEN_MAXIMIZE) );
520  }
521  else
522  {
523  SCIP_CALL( SCIPlpiCreate(lpi, SCIPgetMessagehdlr(scip), "node LP with generated cuts", SCIP_OBJSEN_MINIMIZE) );
524  }
525 
526  /* add cols to the LP interface */
527  SCIP_CALL( SCIPallocBufferArray(scip, &colobj, ncols) );
528  SCIP_CALL( SCIPallocBufferArray(scip, &collb, ncols) );
529  SCIP_CALL( SCIPallocBufferArray(scip, &colub, ncols) );
530  /* gather required column information */
531  for( int i = 0; i < ncols; ++i )
532  {
533  colobj[i] = SCIPcolGetObj(cols[i]);
534  collb[i] = SCIPcolGetLb(cols[i]);
535  colub[i] = SCIPcolGetUb(cols[i]);
536  }
537  /* add cols */
538  SCIP_CALL( SCIPlpiAddCols(*lpi, ncols, colobj, collb, colub, NULL, 0, NULL, NULL, NULL) );
539  SCIPfreeBufferArray(scip, &colub);
540  SCIPfreeBufferArray(scip, &collb);
541  SCIPfreeBufferArray(scip, &colobj);
542 
543  /* add rows to the LP interface */
544  /* find number of nonzeroes in rows */
545  nrownonz = 0;
546  for( int i = 0; i < nrows; ++i )
547  {
548  assert(!(SCIPisInfinity(scip, -SCIProwGetLhs(rows[i])) && SCIPisInfinity(scip, SCIProwGetRhs(rows[i]))));
549  nrownonz += SCIProwGetNLPNonz(rows[i]);
550  }
551  SCIP_CALL( SCIPallocBufferArray(scip, &rowcolinds, nrownonz) );
552  SCIP_CALL( SCIPallocBufferArray(scip, &rowvals, nrownonz) );
553  SCIP_CALL( SCIPallocBufferArray(scip, &rowbegs, nrows + 1) );
554  SCIP_CALL( SCIPallocBufferArray(scip, &rowlhs, nrows) );
555  SCIP_CALL( SCIPallocBufferArray(scip, &rowrhs, nrows) );
556  /* gather required row information */
557  rowbegs[0] = 0;
558  pinf = SCIPlpiInfinity(*lpi);
559  ninf = -SCIPlpiInfinity(*lpi);
560  for( int i = 0; i < nrows; ++i )
561  {
562  nrownonz = SCIProwGetNLPNonz(rows[i]);
563  assert(nrownonz <= ncols);
564  rowval = SCIProwGetVals(rows[i]);
565  rowcols = SCIProwGetCols(rows[i]);
566 
567  rowbegs[i + 1] = rowbegs[i] + nrownonz;
568  rowconst = SCIProwGetConstant(rows[i]);
569  rowlhs[i] = SCIPisInfinity(scip, -SCIProwGetLhs(rows[i])) ? ninf : SCIProwGetLhs(rows[i]) - rowconst;
570  rowrhs[i] = SCIPisInfinity(scip, SCIProwGetRhs(rows[i])) ? pinf : SCIProwGetRhs(rows[i]) - rowconst;
571 
572  for( int j = 0; j < nrownonz; ++j )
573  {
574  collppos = SCIPcolGetLPPos(rowcols[j]);
575  assert(collppos >= 0);
576  assert(collppos <= ncols);
577 
578  rowcolinds[rowbegs[i] + j] = collppos;
579  rowvals[rowbegs[i] + j] = rowval[j];
580  }
581  }
582  /* add rows */
583  SCIP_CALL( SCIPlpiAddRows(*lpi, nrows, rowlhs, rowrhs, NULL, rowbegs[nrows], rowbegs, rowcolinds, rowvals) );
584 
585  /* get warm starting info */
586  SCIP_CALL( SCIPgetLPI(scip, &wslpi) );
587  SCIP_CALL( SCIPlpiGetState(wslpi, blkmem, &lpistate) );
588  }
589  /* if there are any cuts, then add the cuts that were not added earlier to the LPI */
590  else
591  {
592  assert(nrows + ncuts >= sepadata->nrowsinhardcutslp);
593 
594  /* get warm starting info */
595  wslpi = *lpi;
596  SCIP_CALL( SCIPlpiGetState(wslpi, blkmem, &lpistate) );
597 
598  /* find number of nonzeros in cuts and allocate memory */
599  nrownonz = 0;
600  pinf = SCIPlpiInfinity(*lpi);
601  ninf = -SCIPlpiInfinity(*lpi);
602  for( int i = sepadata->nrowsinhardcutslp - nrows; i < ncuts; ++i )
603  {
604  assert(!(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) && SCIPisInfinity(scip, SCIProwGetRhs(cuts[i]))));
605  assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])));
606  assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])));
607  nrownonz += SCIProwGetNNonz(cuts[i]);
608  }
609  SCIP_CALL( SCIPallocBufferArray(scip, &rowcolinds, nrownonz) );
610  SCIP_CALL( SCIPallocBufferArray(scip, &rowvals, nrownonz) );
611  SCIP_CALL( SCIPallocBufferArray(scip, &rowbegs, (ncuts - sepadata->nrowsinhardcutslp + nrows) + 1) );
612  SCIP_CALL( SCIPallocBufferArray(scip, &rowlhs, (ncuts - sepadata->nrowsinhardcutslp + nrows)) );
613  SCIP_CALL( SCIPallocBufferArray(scip, &rowrhs, (ncuts - sepadata->nrowsinhardcutslp + nrows)) );
614 
615  /* gather required cut information */
616  rowbegs[0] = 0;
617  for( int i = sepadata->nrowsinhardcutslp - nrows; i < ncuts; ++i )
618  {
619  nrownonz = SCIProwGetNNonz(cuts[i]);
620  assert(nrownonz <= ncols);
621  rowval = SCIProwGetVals(cuts[i]);
622  rowcols = SCIProwGetCols(cuts[i]);
623 
624  rowbegs[i - sepadata->nrowsinhardcutslp + nrows + 1] = rowbegs[i - sepadata->nrowsinhardcutslp + nrows] +
625  nrownonz;
626  rowconst = SCIProwGetConstant(cuts[i]);
627  rowlhs[i - sepadata->nrowsinhardcutslp + nrows] = SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) ? ninf :
628  SCIProwGetLhs(cuts[i]) - rowconst;
629  rowrhs[i - sepadata->nrowsinhardcutslp + nrows] = SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])) ? pinf :
630  SCIProwGetRhs(cuts[i]) - rowconst;
631 
632  for( int j = 0; j < nrownonz; ++j )
633  {
634  collppos = SCIPcolGetLPPos(rowcols[j]);
635  assert(collppos >= 0);
636  assert(collppos <= ncols);
637 
638  rowcolinds[rowbegs[i - sepadata->nrowsinhardcutslp + nrows] + j] = collppos;
639  rowvals[rowbegs[i - sepadata->nrowsinhardcutslp + nrows] + j] = rowval[j];
640  }
641  }
642 
643  /* add cuts */
644  SCIP_CALL( SCIPlpiAddRows(*lpi, (ncuts - sepadata->nrowsinhardcutslp + nrows), rowlhs, rowrhs, NULL,
645  rowbegs[(ncuts - sepadata->nrowsinhardcutslp + nrows)], rowbegs, rowcolinds, rowvals) );
646  }
647 
648  /* set warm starting basis */
649  SCIP_CALL( SCIPlpiSetState(*lpi, blkmem, lpistate) );
650 
651  /* reset remaining sepadata values */
652  sepadata->nrowsinhardcutslp = nrows + ncuts;
653 
654  /* free memory */
655  SCIP_CALL( SCIPlpiFreeState(*lpi, blkmem, &lpistate) );
656  SCIPfreeBufferArray(scip, &rowrhs);
657  SCIPfreeBufferArray(scip, &rowlhs);
658  SCIPfreeBufferArray(scip, &rowbegs);
659  SCIPfreeBufferArray(scip, &rowvals);
660  SCIPfreeBufferArray(scip, &rowcolinds);
661 
662  return SCIP_OKAY;
663 }
664 
665 /** free separator data */
666 static
668  SCIP* scip, /**< SCIP data structure */
669  SCIP_SEPADATA** sepadata /**< separator data structure */
670  )
671 {
672  assert(scip != NULL);
673  assert(sepadata != NULL);
674  assert(*sepadata != NULL);
675 
676  if( (*sepadata)->lpiwithhardcuts != NULL )
677  {
678  SCIP_CALL( SCIPlpiFree(&((*sepadata)->lpiwithhardcuts)) );
679  }
680 
681  (*sepadata)->nrowsinhardcutslp = 0;
682  (*sepadata)->nrunsinsoftcutslp = 0;
683  (*sepadata)->ncalls = 0;
684  (*sepadata)->nmaxperroundlpiters = 0;
685  (*sepadata)->nmaxrootlpiters = 0;
686  (*sepadata)->nrootlpiters = 0;
687  (*sepadata)->nmaxtotallpiters = 0;
688  (*sepadata)->ntotallpiters = 0;
689  (*sepadata)->nmaxperroundcutsroot = 0;
690  (*sepadata)->nmaxperroundcuts = 0;
691  (*sepadata)->nmaxtotalcuts = 0;
692  (*sepadata)->ntotalcuts = 0;
693 
694  SCIPfreeBlockMemory(scip, sepadata);
695 
696  return SCIP_OKAY;
697 }
698 
699 /** update mu parameter which is used as a factor in the step length calculation; refer to the top of the file for a
700  * description of the formula.
701  */
702 /* @todo some adaptive strategy like constant after certain changes? */
703 static
705  SCIP* scip, /**< SCIP data structure */
706  SCIP_SEPADATA* sepadata, /**< separator data structure */
707  int subgradientiternum, /**< subgradient iteration number */
708  SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */
709  SCIP_Real* lagrangianvals, /**< vector of Lagrangian values found so far */
710  SCIP_Real bestlagrangianval, /**< best Lagrangian value found so far */
711  SCIP_Real avglagrangianval, /**< rolling average of the Lagrangian values found so far */
712  SCIP_Real* muparam, /**< mu parameter to be updated */
713  SCIP_Bool* backtrack /**< whether mu parameter has been backtracked */
714  )
715 {
716  SCIP_Real delta;
717  SCIP_Real deltaslab1ub;
718  SCIP_Real deltaslab2ub;
719  SCIP_Real muslab1factor;
720  SCIP_Real muslab2factor;
721  SCIP_Real muslab3factor;
722  int maxiters;
723  int i;
725  *backtrack = FALSE;
726 
727  /* update the mu parameter only if it is not set to be a constant value */
728  if( !sepadata->muparamconst )
729  {
730  delta = ubparam - bestlagrangianval;
731  deltaslab1ub = MIN(sepadata->deltaslab1ub, sepadata->deltaslab2ub);
732  deltaslab2ub = MAX(sepadata->deltaslab1ub, sepadata->deltaslab2ub);
733  /* ensure that the ordering of different user input parameters is as expected */
734  if( SCIPisPositive(scip, sepadata->muslab1factor - sepadata->muslab2factor) )
735  {
736  if( SCIPisPositive(scip, sepadata->muslab2factor - sepadata->muslab3factor) )
737  {
738  muslab1factor = sepadata->muslab1factor;
739  muslab2factor = sepadata->muslab2factor;
740  muslab3factor = sepadata->muslab3factor;
741  }
742  else
743  {
744  if( SCIPisPositive(scip, sepadata->muslab1factor - sepadata->muslab3factor) )
745  {
746  muslab1factor = sepadata->muslab1factor;
747  muslab2factor = sepadata->muslab3factor;
748  muslab3factor = sepadata->muslab2factor;
749  }
750  else
751  {
752  muslab1factor = sepadata->muslab3factor;
753  muslab2factor = sepadata->muslab1factor;
754  muslab3factor = sepadata->muslab2factor;
755  }
756  }
757  }
758  else
759  {
760  if( SCIPisPositive(scip, sepadata->muslab1factor - sepadata->muslab3factor) )
761  {
762  muslab1factor = sepadata->muslab2factor;
763  muslab2factor = sepadata->muslab1factor;
764  muslab3factor = sepadata->muslab3factor;
765  }
766  else
767  {
768  if( SCIPisPositive(scip, sepadata->muslab2factor - sepadata->muslab3factor) )
769  {
770  muslab1factor = sepadata->muslab2factor;
771  muslab2factor = sepadata->muslab3factor;
772  muslab3factor = sepadata->muslab1factor;
773  }
774  else
775  {
776  muslab1factor = sepadata->muslab3factor;
777  muslab2factor = sepadata->muslab2factor;
778  muslab3factor = sepadata->muslab1factor;
779  }
780  }
781  }
782 
783  maxiters = MIN(sepadata->nmaxconsecitersformuupdate, sepadata->nmaxlagrangianvalsforavg);
784  i = -1;
785 
786  /* if certain number of iterations are done, then check for a possibility of backtracking and apply accordingly */
787  if( subgradientiternum >= maxiters )
788  {
789  for( i = subgradientiternum - maxiters; i < subgradientiternum; i++ )
790  {
791  if( SCIPisGE(scip, lagrangianvals[i], bestlagrangianval - delta) )
792  break;
793  }
794 
795  if( i == subgradientiternum )
796  {
797  *muparam *= sepadata->mubacktrackfactor;
798  *backtrack = TRUE;
799  }
800  }
801 
802  /* update mu parameter based on the different between best and average Lagrangian values */
803  if( (subgradientiternum < maxiters) || (i >= 0 && i < subgradientiternum) )
804  {
805  if( bestlagrangianval - avglagrangianval < deltaslab1ub * delta )
806  *muparam *= muslab1factor;
807  else if( bestlagrangianval - avglagrangianval < deltaslab2ub * delta )
808  *muparam *= muslab2factor;
809  else
810  *muparam *= muslab3factor;
811  }
812 
813  /* reset the mu parameter to within its bounds */
814  *muparam = MAX(*muparam, sepadata->muparamlb);
815  *muparam = MIN(*muparam, sepadata->muparamub);
816  }
817 
818  return SCIP_OKAY;
819 }
820 
821 /** update subgradient, i.e., residuals of generated cuts */
822 /* @note: assumed that \f$i^{th}\f$ cut is of the form \f${\alpha^i}^T x \leq {\alpha^i_0}\f$ */
823 static
824 void updateSubgradient(
825  SCIP* scip, /**< SCIP data structure */
826  SCIP_SOL* sol, /**< LP solution used in updating subgradient vector */
827  SCIP_ROW** cuts, /**< cuts generated so far */
828  int ncuts, /**< number of cuts generated so far */
829  SCIP_Real* subgradient, /**< vector of subgradients to be updated */
830  SCIP_Real* dualvector, /**< Lagrangian multipliers */
831  SCIP_Bool* subgradientzero, /**< whether the subgradient vector is all zero */
832  int* ncutviols, /**< number of violations of generated cuts */
833  SCIP_Real* maxcutviol, /**< maximum violation of generated cuts */
834  int* nnzsubgradientdualprod, /**< number of nonzero products of subgradient vector and Lagrangian multipliers (i.e.,
835  number of complementarity slackness violations) */
836  SCIP_Real* maxnzsubgradientdualprod /**< maximum value of nonzero products of subgradient vector and Lagrangian multipliers
837  (i.e., maximum value of complementarity slackness violations) */
838  )
839 {
840  int nzerosubgradient;
841  SCIP_Real term;
842 
843  assert(subgradientzero != NULL);
844  assert(ncutviols != NULL);
845  assert(maxcutviol != NULL);
846  assert(nnzsubgradientdualprod != NULL);
847  assert(maxnzsubgradientdualprod != NULL);
848 
849  *ncutviols = 0;
850  *maxcutviol = 0.0;
851  *nnzsubgradientdualprod = 0;
852  *maxnzsubgradientdualprod = 0.0;
853  nzerosubgradient = 0;
854  *subgradientzero = FALSE;
855 
856  /* for each cut, calculate the residual along with various violation metrics */
857  for( int i = 0; i < ncuts; i++ )
858  {
859  assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])));
860  assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])));
861  subgradient[i] = SCIPgetRowSolActivity(scip, cuts[i], sol) + SCIProwGetConstant(cuts[i]) - SCIProwGetRhs(cuts[i]);
862  if( SCIPisFeasZero(scip, subgradient[i]) )
863  {
864  subgradient[i] = 0.0;
865  nzerosubgradient++;
866  }
867  else
868  {
869  /* check for cut violation */
870  if( SCIPisFeasPositive(scip, subgradient[i]) )
871  {
872  (*ncutviols)++;
873  *maxcutviol = MAX(*maxcutviol, subgradient[i]);
874  }
875 
876  /* check for violation of complementarity slackness associated with the cut */
877  if( !SCIPisZero(scip, subgradient[i] * dualvector[i]) )
878  {
879  (*nnzsubgradientdualprod)++;
880  term = REALABS(subgradient[i] * dualvector[i]);
881  *maxnzsubgradientdualprod = MAX(*maxnzsubgradientdualprod, term);
882  }
883  }
884  }
885 
886  /* indicator for all zero subgradient vector */
887  if( nzerosubgradient == ncuts )
888  {
889  *subgradientzero = TRUE;
890  }
891 }
892 
893 /** update Lagrangian value, i.e., optimal value of the Lagrangian dual with fixed multipliers */
894 static
896  SCIP* scip, /**< SCIP data structure */
897  SCIP_Real objval, /**< objective value of the Lagrangian dual with fixed multipliers */
898  SCIP_Real* dualvector, /**< Lagrangian multipliers */
899  SCIP_ROW** cuts, /**< cuts generated so far */
900  int ncuts, /**< number of cuts generated so far */
901  SCIP_Real* lagrangianval /**< Lagrangian value to be updated */
902  )
903 {
904  *lagrangianval = objval;
905 
906  for( int i = 0; i < ncuts; i++ )
907  {
908  assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])));
909  assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])));
910  *lagrangianval += dualvector[i] * (SCIProwGetConstant(cuts[i]) - SCIProwGetRhs(cuts[i]));
911  }
912 }
913 
914 /** update step length based on various input arguments; refer to the top of the file for an expression */
915 static
916 void updateStepLength(
917  SCIP* scip, /**< SCIP data structure */
918  SCIP_Real muparam, /**< mu parameter used as a factor for step length */
919  SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */
920  SCIP_Real lagrangianval, /**< Lagrangian value */
921  SCIP_Real* subgradient, /**< subgradient vector */
922  int ncuts, /**< number of cuts generated so far */
923  SCIP_Real* steplength /**< step length to be updated */
924  )
925 {
926  SCIP_Real normsquared = 0.0;
927 
928  for( int i = 0; i < ncuts; i++ )
929  {
930  normsquared += SQR(subgradient[i]);
931  }
932 
933  if( !SCIPisFeasZero(scip, normsquared) )
934  {
935  *steplength = (muparam * (ubparam - lagrangianval))/(normsquared); /*lint !e795*/
936  }
937 }
938 
939 /** update the ball radius (based on various violation metrics) that is used for stabilization of Lagrangian multipliers */
940 static
941 void updateBallRadius(
942  SCIP* scip, /**< SCIP data structure */
943  SCIP_SEPADATA* sepadata, /**< separator data structure */
944  SCIP_Real maxviolscore, /**< weighted average of maximum value of generated cut violations and maximum value of
945  complementarity slackness violations, in the current iteration */
946  SCIP_Real maxviolscoreold, /**< weighted average of maximum value of generated cut violations and maximum value of
947  complementarity slackness violations, in the previous iteration */
948  SCIP_Real nviolscore, /**< weighted average of number of generated cut violations and number of complementarity
949  slackness violations, in the current iteration */
950  SCIP_Real nviolscoreold, /**< weighted average of number of generated cut violations and number of complementarity
951  slackness violations, in the previous iteration */
952  int nlpiters, /**< number of LP iterations taken for solving the Lagrangian dual with fixed multipliers in
953  current iteration */
954  SCIP_Real* ballradius /**< norm ball radius to be updated */
955  )
956 {
957  SCIP_Bool maxviolscoreimproved;
958  SCIP_Bool nviolscoreimproved;
959 
960  assert(ballradius != NULL);
962  maxviolscoreimproved = !SCIPisNegative(scip, maxviolscoreold - maxviolscore);
963  nviolscoreimproved = !SCIPisNegative(scip, nviolscoreold - nviolscore);
964 
965  if( maxviolscoreimproved && nviolscoreimproved )
966  {
967  /* both the maximum violation and number of violations scores have become better, so, increase the radius */
968  if( sepadata->optimalfacepriority <= 1 )
969  {
970  *ballradius *= 2.0;
971  *ballradius = MIN(*ballradius, sepadata->radiusmax);
972  }
973  else
974  {
975  *ballradius *= 1.5;
976  *ballradius = MIN(*ballradius, sepadata->radiusmax/2.0);
977  }
978  }
979  else if( !maxviolscoreimproved && !nviolscoreimproved )
980  {
981  /* both the maximum violation and number of violations scores have become worse, so, decrease the radius */
982  *ballradius *= 0.5;
983  *ballradius = MAX(*ballradius, sepadata->radiusmin);
984  }
985  else if( nlpiters == 0 )
986  {
987  /* only one among the maximum violation and number of violations scores has become better, and the LP basis did
988  * not change (i.e., nlpters = 0), so, increase the radius slightly */
989  if( sepadata->optimalfacepriority <= 1 )
990  {
991  *ballradius *= 1.5;
992  *ballradius = MIN(*ballradius, sepadata->radiusmax);
993  }
994  else
995  {
996  *ballradius *= 1.2;
997  *ballradius = MIN(*ballradius, sepadata->radiusmax/2.0);
998  }
999  }
1000 }
1001 
1002 /** projection of Lagrangian multipliers onto L1-norm ball. This algorithm is based on the following article.
1003  *
1004  * Condat L. (2016).@n
1005  * Fast projection onto the simplex and the \f$l_1\f$ ball.@n
1006  * Mathematical Programming, 158, 1-2, 575-585.
1007  *
1008  */
1009 static
1011  SCIP* scip, /**< SCIP data structure */
1012  SCIP_Real* dualvector, /**< Lagrangian multipliers to be projected onto L1-norm vall */
1013  int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1014  SCIP_Real radius /**< radius of the L1-norm ball */
1015  )
1016 {
1017  SCIP_Real* temp1vals;
1018  SCIP_Real* temp2vals;
1019  SCIP_Real pivotparam;
1020  SCIP_Real val;
1021  SCIP_Real term;
1022  SCIP_Bool temp1changed;
1023  int ntemp1removed;
1024  int* temp1inds;
1025  int* temp2inds;
1026  int temp1len;
1027  int temp2len;
1028 
1029  assert(!SCIPisNegative(scip, radius));
1030  val = REALABS(dualvector[0]);
1031  /* calculate the L1-norm of the Lagrangian multipliers */
1032  for( int i = 1; i < dualvectorlen; i++ )
1033  {
1034  val += REALABS(dualvector[i]);
1035  }
1036 
1037  /* if the vector of Lagrangian multipliers lies outside the L1-norm ball, then do the projection */
1038  if( SCIPisGT(scip, val, radius) )
1039  {
1040  SCIP_CALL( SCIPallocCleanBufferArray(scip, &temp1vals, dualvectorlen) );
1041  SCIP_CALL( SCIPallocCleanBufferArray(scip, &temp2vals, dualvectorlen) );
1042  SCIP_CALL( SCIPallocBufferArray(scip, &temp1inds, dualvectorlen) );
1043  SCIP_CALL( SCIPallocBufferArray(scip, &temp2inds, dualvectorlen) );
1044  for( int i = 0; i < dualvectorlen; i++ )
1045  {
1046  temp1inds[i] = -1;
1047  temp2inds[i] = -1;
1048  }
1049  temp2len = 0;
1050 
1051  temp1vals[0] = REALABS(dualvector[0]);
1052  temp1inds[0] = 0;
1053  temp1len = 1;
1054  pivotparam = REALABS(dualvector[0]) - radius;
1055 
1056  for( int i = 1; i < dualvectorlen; i++ )
1057  {
1058  if( SCIPisGT(scip, REALABS(dualvector[i]), pivotparam) )
1059  {
1060  pivotparam += ((REALABS(dualvector[i]) - pivotparam) / (temp1len + 1));
1061  if( SCIPisGT(scip, pivotparam, REALABS(dualvector[i]) - radius) )
1062  {
1063  temp1vals[temp1len] = REALABS(dualvector[i]);
1064  temp1inds[temp1len] = i;
1065  temp1len++;
1066  }
1067  else
1068  {
1069  for( int j = 0; j < temp1len; j++ )
1070  {
1071  temp2vals[temp2len + j] = temp1vals[j];
1072  temp2inds[temp2len + j] = temp1inds[j];
1073  }
1074  temp2len += temp1len;
1075  temp1vals[0] = REALABS(dualvector[i]);
1076  temp1inds[0] = i;
1077  temp1len = 1;
1078  pivotparam = REALABS(dualvector[i]) - radius;
1079  }
1080  }
1081  }
1082 
1083  for( int i = 0; i < temp2len; i++ )
1084  {
1085  if( SCIPisGT(scip, temp2vals[i], pivotparam) )
1086  {
1087  temp1vals[temp1len] = temp2vals[i];
1088  temp1inds[temp1len] = temp2inds[i];
1089  temp1len++;
1090  pivotparam += ((temp2vals[i] - pivotparam) / temp1len);
1091  }
1092  }
1093 
1094  temp1changed = TRUE;
1095  ntemp1removed = 0;
1096  while( temp1changed )
1097  {
1098  temp1changed = FALSE;
1099 
1100  for( int i = 0; i < temp1len; i++ )
1101  {
1102  /* @note: the third condition (temp1len - ntemp1removed > 0) is true as long as the first condition
1103  * (temp1inds[i] >= 0) is true.
1104  */
1105  if( (temp1inds[i] >= 0) && SCIPisLE(scip, temp1vals[i], pivotparam) )
1106  {
1107  temp1inds[i] = -1;
1108  temp1changed = TRUE;
1109  ntemp1removed++;
1110  assert(temp1len - ntemp1removed > 0);
1111  /* coverity[divide_by_zero] */
1112  pivotparam += ((pivotparam - temp1vals[i]) / (temp1len - ntemp1removed));
1113  }
1114  }
1115  }
1116 
1117  for( int i = 0; i < dualvectorlen; i++ )
1118  {
1119  term = REALABS(dualvector[i]);
1120  val = MAX(term - pivotparam, 0.0);
1121 
1122  if( SCIPisPositive(scip, dualvector[i]) )
1123  {
1124  dualvector[i] = val;
1125  }
1126  else if( SCIPisNegative(scip, dualvector[i]) )
1127  {
1128  dualvector[i] = -val;
1129  }
1130  }
1131 
1132  /* free memory */
1133  for( int i = 0; i < dualvectorlen; i++ )
1134  {
1135  temp2vals[i] = 0.0;
1136  temp1vals[i] = 0.0;
1137  }
1138  SCIPfreeBufferArray(scip, &temp2inds);
1139  SCIPfreeBufferArray(scip, &temp1inds);
1140  SCIPfreeCleanBufferArray(scip, &temp2vals);
1141  SCIPfreeCleanBufferArray(scip, &temp1vals);
1142  }
1143 
1144  return SCIP_OKAY;
1145 }
1146 
1147 /** projection of Lagrangian multipliers onto L2-norm ball */
1148 static
1149 void l2BallProjection(
1150  SCIP* scip, /**< SCIP data structure */
1151  SCIP_Real* dualvector, /**< Lagrangian multipliers to be projected onto L2-norm vall */
1152  int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1153  SCIP_Real radius /**< radius of the L2-norm ball */
1154  )
1155 {
1156  SCIP_Real l2norm;
1157  SCIP_Real factor;
1158 
1159  assert(!SCIPisNegative(scip, radius));
1160 
1161  l2norm = 0.0;
1162  /* calculate the L2-norm of the Lagrangian multipliers */
1163  for( int i = 0; i < dualvectorlen; i++ )
1164  {
1165  l2norm += SQR(dualvector[i]);
1166  }
1167  l2norm = sqrt(l2norm);
1168  factor = radius/(1.0 + l2norm);
1170  /* if the vector of Lagrangian multipliers is outside the L2-norm ball, then do the projection */
1171  if( SCIPisGT(scip, l2norm, radius) && SCIPisLT(scip, factor, 1.0) )
1172  {
1173  for( int i = 0; i < dualvectorlen; i++ )
1174  {
1175  dualvector[i] *= factor;
1176  }
1177  }
1178 }
1179 
1180 /** projection of Lagrangian multipliers onto L_infinity-norm ball */
1181 static
1182 void linfBallProjection(
1183  SCIP* scip, /**< SCIP data structure */
1184  SCIP_Real* dualvector, /**< Lagrangian multipliers to be projected onto L_infinity-norm vall */
1185  int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1186  SCIP_Real radius /**< radius of the L_infinity-norm ball */
1187  )
1188 {
1189  assert(!SCIPisNegative(scip, radius));
1190 
1191  /* if the vector of Lagrangian multipliers is outside the L_infinity-norm ball, then do the projection */
1192  for( int i = 0; i < dualvectorlen; i++ )
1193  {
1194  if( SCIPisLT(scip, dualvector[i], -radius) )
1195  {
1196  dualvector[i] = -radius;
1197  }
1198  else if( SCIPisGT(scip, dualvector[i], radius) )
1199  {
1200  dualvector[i] = radius;
1201  }
1202  }
1203 }
1204 
1205 /** weighted Lagrangian multipliers based on a given vector as stability center */
1206 /* @todo calculate weight outside this function and pass it (so that this function becomes generic and independent of
1207  * the terminology related to best Lagrangian multipliers)
1208  */
1209 static
1211  SCIP* scip, /**< SCIP data structure */
1212  SCIP_SEPADATA* sepadata, /**< separator data structure */
1213  SCIP_Real* dualvector, /**< Lagrangian multipliers */
1214  int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1215  SCIP_Real* stabilitycenter, /**< stability center (i.e., core vector of Lagrangian multipliers) */
1216  int stabilitycenterlen, /**< length of the stability center */
1217  int nbestdualupdates, /**< number of best Lagrangian values found so far */
1218  int totaliternum /**< total number of iterations of the relax-and-cut algorithm performed so far */
1219  )
1220 {
1221  SCIP_Real constant;
1222  SCIP_Real weight;
1223  SCIP_Real alpha;
1224 
1225  constant = MAX(2.0, sepadata->constant);
1226  /* weight factor from the literature on Dantzig-Wolfe decomposition stabilization schemes */
1227  weight = MIN(constant, (totaliternum + 1 + nbestdualupdates) / 2.0);
1228  alpha = 1.0 / weight;
1229 
1230  assert(dualvectorlen >= stabilitycenterlen);
1231 
1232  /* weighted Lagrangian multipliers */
1233  for( int i = 0; i < stabilitycenterlen; i++ )
1234  {
1235  dualvector[i] = alpha * dualvector[i] + (1 - alpha) * stabilitycenter[i];
1236  }
1237  for( int i = stabilitycenterlen; i < dualvectorlen; i++ )
1238  {
1239  dualvector[i] = alpha * dualvector[i];
1240  }
1241 
1242  return SCIP_OKAY;
1243 }
1244 
1245 /** stabilize Lagrangian multipliers */
1246 static
1248  SCIP* scip, /**< SCIP data structure */
1249  SCIP_SEPADATA* sepadata, /**< separator data structure */
1250  SCIP_Real* dualvector, /**< Lagrangian multipliers */
1251  int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1252  SCIP_Real* bestdualvector, /**< best Lagrangian multipliers found so far */
1253  int bestdualvectorlen, /**< length of the best Lagrangian multipliers vector */
1254  int nbestdualupdates, /**< number of best Lagrangian values found so far */
1255  int subgradientiternum, /**< iteration number of the subgradient algorithm */
1256  int totaliternum, /**< total number of iterations of the relax-and-cut algorithm performed so far */
1257  SCIP_Real maxviolscore, /**< weighted average of maximum value of generated cut violations and maximum value of
1258  complementarity slackness violations, in the current iteration */
1259  SCIP_Real maxviolscoreold, /**< weighted average of maximum value of generated cut violations and maximum value of
1260  complementarity slackness violations, in the previous iteration */
1261  SCIP_Real nviolscore, /**< weighted average of number of generated cut violations and number of complementarity
1262  slackness violations, in the current iteration */
1263  SCIP_Real nviolscoreold, /**< weighted average of number of generated cut violations and number of complementarity
1264  slackness violations, in the previous iteration */
1265  int nlpiters, /**< number of LP iterations taken for solving the Lagrangian dual with fixed multipliers in
1266  current iteration */
1267  SCIP_Real* ballradius /**< norm ball radius */
1268  )
1269 {
1270  if( sepadata->projectiontype > 0 )
1271  {
1272  if( subgradientiternum >= 1 )
1273  {
1274  /* update the ball radius */
1275  updateBallRadius(scip, sepadata, maxviolscore, maxviolscoreold, nviolscore, nviolscoreold, nlpiters,
1276  ballradius);
1277  }
1278 
1279  if( sepadata->projectiontype == 1 )
1280  {
1281  /* projection of Lagrangian multipliers onto L1-norm ball */
1282  SCIP_CALL( l1BallProjection(scip, dualvector, dualvectorlen, *ballradius) );
1283  }
1284  else if( sepadata->projectiontype == 2 )
1285  {
1286  /* projection of Lagrangian multipliers onto L2-norm ball */
1287  l2BallProjection(scip, dualvector, dualvectorlen, *ballradius);
1288  }
1289  else if( sepadata->projectiontype == 3 )
1290  {
1291  /* projection of Lagrangian multipliers onto L_inf-norm ball */
1292  linfBallProjection(scip, dualvector, dualvectorlen, *ballradius);
1293  }
1294  }
1295 
1296  if( sepadata->stabilitycentertype == 1 )
1297  {
1298  /* weighted Lagrangian multipliers based on best Langrangian multipliers as stability center */
1299  SCIP_CALL( weightedDualVector(scip, sepadata, dualvector, dualvectorlen, bestdualvector,
1300  bestdualvectorlen, nbestdualupdates, totaliternum) );
1301  }
1302 
1303  return SCIP_OKAY;
1304 }
1305 
1306 /** update Lagrangian multipliers */
1307 static
1309  SCIP* scip, /**< SCIP data structure */
1310  SCIP_SEPADATA* sepadata, /**< separator data structure */
1311  SCIP_Real* dualvector1, /**< Lagrangian multipliers vector to be updated */
1312  SCIP_Real* dualvector2, /**< Lagrangian multipliers vector used for backtracking */
1313  int dualvector2len, /**< length of the Lagrangian multipliers vector used for backtracking */
1314  int ndualvector2updates,/**< number of best Lagrangian values found so far */
1315  int subgradientiternum, /**< iteration number of the subgradient algorithm */
1316  int totaliternum, /**< total number of iterations of the relax-and-cut algorithm performed so far */
1317  SCIP_Real steplength, /**< step length used for updating Lagrangian multipliers */
1318  SCIP_Real* subgradient, /**< subgradient vector */
1319  int ncuts, /**< number of generated cuts so far */
1320  SCIP_Bool backtrack, /**< whether the Lagrangian multipliers need to be backtracked */
1321  SCIP_Real maxviolscore, /**< weighted average of maximum value of generated cut violations and maximum value of
1322  complementarity slackness violations, in the current iteration */
1323  SCIP_Real maxviolscoreold, /**< weighted average of maximum value of generated cut violations and maximum value of
1324  complementarity slackness violations, in the previous iteration */
1325  SCIP_Real nviolscore, /**< weighted average of number of generated cut violations and number of complementarity
1326  slackness violations, in the current iteration */
1327  SCIP_Real nviolscoreold, /**< weighted average of number of generated cut violations and number of complementarity
1328  slackness violations, in the previous iteration */
1329  int nlpiters, /**< number of LP iterations taken for solving the Lagrangian dual with fixed multipliers in
1330  current iteration */
1331  SCIP_Bool* dualvecsdiffer, /**< whether the updated Lagrangian multipliers differ from the old one */
1332  SCIP_Real* ballradius /**< norm ball radius */
1333  )
1334 {
1335  SCIP_Real* dualvector1copy;
1336 
1337  assert(dualvector2len <= ncuts);
1338  assert(dualvecsdiffer != NULL);
1339 
1340  *dualvecsdiffer = FALSE;
1341  /* @todo do allocation and free operations outside only once instead of every time this function is called? */
1342  /* copy of the Lagrangian multipliers to be used to check if the updated vector is different than this */
1343  SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualvector1copy, ncuts) );
1344  for( int i = 0; i < ncuts; i++ )
1345  {
1346  dualvector1copy[i] = dualvector1[i];
1347  }
1348 
1349  /* if backtracking was not identified at the time of the mu parameter update, then update the Lagrangian multipliers
1350  * based on the given subgradient vector
1351  */
1352  if( !backtrack )
1353  {
1354  assert((subgradient != NULL) || (ncuts == 0));
1355  assert(subgradientiternum >= 0);
1356 
1357  /* update Lagrangian multipliers */
1358  for( int i = 0; i < ncuts; i++ )
1359  {
1360  dualvector1[i] += steplength * subgradient[i];
1361  }
1362 
1363  /* projection onto non-negative orthant */
1364  for( int i = 0; i < ncuts; i++ )
1365  {
1366  dualvector1[i] = MAX(dualvector1[i], 0.0);
1367  }
1368 
1369  /* stabilization of Lagrangian multipliers */
1370  SCIP_CALL( stabilizeDualVector(scip, sepadata, dualvector1, ncuts, dualvector2, dualvector2len,
1371  ndualvector2updates, subgradientiternum, totaliternum, maxviolscore, maxviolscoreold, nviolscore,
1372  nviolscoreold, nlpiters, ballradius) );
1373 
1374  /* projection onto non-negative orthant again in case stabilization changed some components negative*/
1375  for( int i = 0; i < ncuts; i++ )
1376  {
1377  dualvector1[i] = MAX(dualvector1[i], 0.0);
1378  }
1379  }
1380  /* if backtracking was identified at the time of the mu parameter update, then backtrack the Lagrangian multipliers
1381  * based on the given backtracking multipliers
1382  */
1383  else
1384  {
1385  for( int i = 0; i < dualvector2len; i++ )
1386  {
1387  dualvector1[i] = dualvector2[i];
1388  }
1389 
1390  for( int i = dualvector2len; i < ncuts; i++ )
1391  {
1392  dualvector1[i] = 0.0;
1393  }
1394  }
1395 
1396  /* identify if the vector of Lagrangian multipliers is indeed different after updating */
1397  for( int i = 0; i < ncuts; i++ )
1398  {
1399  if( !SCIPisEQ(scip, dualvector1[i], dualvector1copy[i]) )
1400  {
1401  *dualvecsdiffer = TRUE;
1402  break;
1403  }
1404  }
1405 
1406  /* free memory */
1407  for( int i = 0; i < ncuts; i++ )
1408  {
1409  dualvector1copy[i] = 0.0;
1410  }
1411  SCIPfreeCleanBufferArray(scip, &dualvector1copy);
1412 
1413  return SCIP_OKAY;
1414 }
1415 
1416 /** check different termination criteria */
1417 /* @note: the criterion based on objvecsdiffer assumes deterministic solving process (i.e., we would get same LP solution
1418  * for "Lagrangian dual with fixed Lagrangian multipliers" when the objective vector remains the same across iterations).
1419  */
1420 /* @todo nlpssolved criterion? */
1421 static
1423  SCIP_SEPADATA* sepadata, /**< separator data structure */
1424  int nnewaddedsoftcuts, /**< number of cuts that were recently penalized and added to the Lagrangian dual's
1425  objective function */
1426  int nyettoaddsoftcuts, /**< number of cuts that are yet to be penalized and added to the Lagrangian dual's
1427  objective function */
1428  SCIP_Bool objvecsdiffer, /**< whether the Lagrangian dual's objective function has changed */
1429  int ngeneratedcurrroundcuts, /**< number of cuts generated in the current separation round */
1430  int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
1431  int ncurrroundlpiters, /**< number of separating LP iterations in the current separation round */
1432  int depth, /**< depth of the current node */
1433  SCIP_Bool* terminate /**< whether to terminate the subgradient algorithm loop */
1434  )
1435 {
1436  *terminate = FALSE;
1437 
1438  /* check if no new cuts were added to the Lagrangian dual, no cuts are remaining to be added, and the objective
1439  * function of the Lagrangian dual with fixed multipliers had not changed from the previous iteration
1440  */
1441  if( (nnewaddedsoftcuts == 0) && (nyettoaddsoftcuts == 0) && !objvecsdiffer )
1442  *terminate = TRUE;
1443 
1444  /* check if allowed number of cuts in this separation round have already been generated */
1445  if( ngeneratedcurrroundcuts >= nmaxgeneratedperroundcuts )
1446  *terminate = TRUE;
1447 
1448  /* check if allowed number of cuts in the tree have already been generated */
1449  if( sepadata->ntotalcuts >= sepadata->nmaxtotalcuts )
1450  *terminate = TRUE;
1451 
1452  /* check if allowed number of simplex iterations in this separation round have already been used up */
1453  if( (sepadata->nmaxperroundlpiters >= 0) && (ncurrroundlpiters >= sepadata->nmaxperroundlpiters) )
1454  *terminate = TRUE;
1455 
1456  /* check if allowed number of simplex iterations in the root node have already been used up */
1457  if( (depth == 0) && (sepadata->nmaxrootlpiters >= 0) && (sepadata->nrootlpiters >= sepadata->nmaxrootlpiters) )
1458  *terminate = TRUE;
1459 
1460  /* check if allowed number of simplex iterations in the tree have already been used up */
1461  if( (sepadata->nmaxtotallpiters >= 0) && (sepadata->ntotallpiters >= sepadata->nmaxtotallpiters) )
1462  *terminate = TRUE;
1463 
1464  return SCIP_OKAY;
1465 }
1466 
1467 /** solve the LP corresponding to the Lagrangian dual with fixed Lagrangian multipliers */
1468 static
1470  SCIP* scip, /**< SCIP data structure */
1471  SCIP_SEPADATA* sepadata, /**< separator data structure */
1472  int depth, /**< depth of the current node in the tree */
1473  SCIP_Real origobjoffset, /**< objective offset in the current node's relaxation */
1474  SCIP_Bool* solfound, /**< whether an LP optimal solution has been found */
1475  SCIP_SOL* sol, /**< data structure to store LP optimal solution, if found */
1476  SCIP_Real* solvals, /**< values of the LP optimal solution, if found */
1477  SCIP_Real* objval, /**< optimal objective value of the LP optimal solution, if found */
1478  int* ncurrroundlpiters /**< number of LP iterations taken for solving Lagrangian dual problems with fixed multipliers
1479  in the current separator round */
1480  )
1481 {
1482  SCIP_Real timelimit;
1483  SCIP_COL** cols;
1484  SCIP_COL* col;
1485  SCIP_VAR* var;
1486  SCIP_LPI* lpi;
1487  SCIP_Bool lperror;
1488  SCIP_Bool cutoff;
1490  SCIP_Longint ntotallpiters;
1491  SCIP_Longint nlpiters;
1492  int ncols;
1493  int iterlimit;
1494 
1495  assert(solfound != NULL);
1496  assert(sol != NULL);
1497  assert(solvals != NULL);
1498  assert(ncurrroundlpiters != NULL);
1499  assert(objval != NULL);
1500 
1501  *solfound = FALSE;
1502  lperror = FALSE;
1503  cutoff = FALSE;
1504  iterlimit = -1;
1505  lpi = sepadata->lpiwithsoftcuts;
1506 
1507  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
1508 
1509  /* set time limit */
1510  SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
1511  if( !SCIPisInfinity(scip, timelimit) )
1512  {
1513  timelimit -= SCIPgetSolvingTime(scip);
1514  if( timelimit <= 0.0 )
1515  {
1516  SCIPdebugMsg(scip, "skip Lagromory cut generation since no time left\n");
1517  goto TERMINATE;
1518  }
1519  /* @note: the following direct LPI call is being used because of the lack of an equivalent function call in
1520  * scip_lp.c (lpSetRealpar exists in lp.c though)
1521  */
1522  SCIP_CALL( SCIPlpiSetRealpar(lpi, SCIP_LPPAR_LPTILIM, timelimit) );
1523  }
1524 
1525  /* find iteration limit */
1526  if( (depth == 0) &&
1527  (sepadata->perrootlpiterfactor >= 0.0 && !SCIPisInfinity(scip, sepadata->perrootlpiterfactor)) )
1528  {
1529  iterlimit = (int)(sepadata->perrootlpiterfactor * SCIPgetNRootFirstLPIterations(scip));
1530  }
1531  else if( (depth > 0) &&
1532  (sepadata->perlpiterfactor >= 0.0 && !SCIPisInfinity(scip, sepadata->perlpiterfactor)) )
1533  {
1534  iterlimit = (int)(sepadata->perlpiterfactor * SCIPgetNNodeInitLPIterations(scip));
1535  }
1536  if( sepadata->nmaxperroundlpiters >= 0 )
1537  {
1538  if( sepadata->nmaxperroundlpiters - *ncurrroundlpiters >= 0 )
1539  {
1540  if( iterlimit >= 0 )
1541  {
1542  iterlimit = MIN(iterlimit, sepadata->nmaxperroundlpiters - *ncurrroundlpiters);
1543  }
1544  else
1545  {
1546  iterlimit = sepadata->nmaxperroundlpiters - *ncurrroundlpiters;
1547  }
1548  }
1549  else
1550  {
1551  iterlimit = 0;
1552  }
1553  }
1554  /* @todo impose a finite iteration limit only when the dualvector changes from zero to non-zero for the first time because
1555  * many simplex pivots are performed in this case even with warm starting (compared to the case when the
1556  * dualvector changes from non-zero to non-zero).
1557  */
1558 
1559  /* solve the LP with an iteration limit and get number of simplex iterations taken */
1560  ntotallpiters = SCIPgetNLPIterations(scip);
1561 
1562  SCIP_CALL( SCIPsolveDiveLP(scip, iterlimit, &lperror, &cutoff) );
1563 
1564  nlpiters = SCIPgetNLPIterations(scip) - ntotallpiters;
1565 
1566  /* get the solution and objective value if optimal */
1567  stat = SCIPgetLPSolstat(scip);
1568  /* @todo is there any way to accept terminations due to iterlimit and timelimit as well? It is not possible
1569  * currently because primal sol is not saved in these cases.
1570  */
1571  /* @note: ideally, only primal feasibility is sufficient. But, there is no such option with SCIPgetLPSolstat. */
1572  if( stat == SCIP_LPSOLSTAT_OPTIMAL )
1573  {
1574  if( SCIPisLPSolBasic(scip) )
1575  {
1576  *solfound = TRUE;
1577 
1578  /* update sol */
1579  for( int i = 0; i < ncols; ++i )
1580  {
1581  col = cols[i];
1582  assert(col != NULL);
1583 
1584  var = SCIPcolGetVar(col);
1585 
1586  solvals[i] = SCIPcolGetPrimsol(col);
1587  SCIP_CALL( SCIPsetSolVal(scip, sol, var, solvals[i]) );
1588  }
1589 
1590  *objval = SCIPgetLPObjval(scip);
1591  *objval = *objval + origobjoffset;
1592  }
1593  }
1594 
1595  /* update some statistics */
1596  if( depth == 0 )
1597  {
1598  sepadata->nrootlpiters += (int)nlpiters;
1599  }
1600  sepadata->ntotallpiters += (int)nlpiters;
1601  *ncurrroundlpiters += (int)nlpiters;
1602 
1603 TERMINATE:
1604  return SCIP_OKAY;
1605 }
1606 
1607 /** solve the LP corresponding to the node relaxation upon adding all the generated cuts */
1608 static
1610  SCIP* scip, /**< SCIP data structure */
1611  SCIP_SEPADATA* sepadata, /**< separator data structure */
1612  SCIP_Bool* solfound, /**< whether an LP optimal solution has been found */
1613  SCIP_SOL* sol, /**< data structure to store LP optimal solution, if found */
1614  SCIP_Real* solvals /**< values of the LP optimal solution, if found */
1615  )
1616 {
1617  SCIP_Real timelimit;
1618  SCIP_COL** cols;
1619  SCIP_COL* col;
1620  SCIP_VAR* var;
1621  int ncols;
1622 
1623  assert(solfound != NULL);
1624  assert(sol != NULL);
1625  assert(solvals != NULL);
1626 
1627  *solfound = FALSE;
1628 
1629  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
1630 
1631  /* set time limit */
1632  SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
1633  if( !SCIPisInfinity(scip, timelimit) )
1634  {
1635  timelimit -= SCIPgetSolvingTime(scip);
1636  if( timelimit <= 0.0 )
1637  {
1638  SCIPdebugMsg(scip, "skip Lagromory cut generation since no time left\n");
1639  goto TERMINATE;
1640  }
1641  SCIP_CALL( SCIPlpiSetRealpar(sepadata->lpiwithhardcuts, SCIP_LPPAR_LPTILIM, timelimit) );
1642  }
1643 
1644  /* solve the LP */
1645  SCIP_CALL( SCIPlpiSolvePrimal(sepadata->lpiwithhardcuts) );
1646 
1647  /* get the solution if primal feasible */
1648  if( SCIPlpiIsPrimalFeasible(sepadata->lpiwithhardcuts) )
1649  {
1650  *solfound = TRUE;
1651  SCIP_CALL( SCIPlpiGetSol(sepadata->lpiwithhardcuts, NULL, solvals, NULL, NULL, NULL) );
1652 
1653  /* update sol */
1654  for( int i = 0; i < ncols; ++i )
1655  {
1656  col = cols[i];
1657  assert(col != NULL);
1658 
1659  var = SCIPcolGetVar(col);
1660 
1661  SCIP_CALL( SCIPsetSolVal(scip, sol, var, solvals[i]) );
1662  }
1663  }
1664 
1665 TERMINATE:
1666  return SCIP_OKAY;
1667 }
1668 
1669 /** construct a cut based on the input cut coefficients, sides, etc */
1670 static
1672  SCIP* scip, /**< SCIP data structure */
1673  SCIP_SEPA* sepa, /**< pointer to the separator */
1674  SCIP_SEPADATA* sepadata, /**< separator data structure */
1675  int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */
1676  int subgradientiternum, /**< iteration number of the subgradient algorithm */
1677  int cutnnz, /**< number of nonzeros in cut */
1678  int* cutinds, /**< column indices in cut */
1679  SCIP_Real* cutcoefs, /**< cut cofficients */
1680  SCIP_Real cutefficacy, /**< cut efficacy */
1681  SCIP_Real cutrhs, /**< RHS of cut */
1682  SCIP_Bool cutislocal, /**< whether cut is local */
1683  int cutrank, /**< rank of cut */
1684  SCIP_ROW** generatedcuts, /**< array of generated cuts */
1685  SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut
1686  generations */
1687  int ngeneratedcurrroundcuts, /**< number of cuts generated until the previous basis in the current separation round */
1688  int* ngeneratednewcuts, /**< number of new cuts generated using the current basis */
1689  SCIP_Bool* cutoff /**< should the current node be cutoff? */
1690  )
1692  SCIP_COL** cols;
1693  SCIP_Real minactivity;
1694  SCIP_Real maxactivity;
1695  SCIP_Real lhs;
1696  SCIP_Real rhs;
1697 
1698  assert(scip != NULL);
1699  assert(mainiternum >= 0);
1700  assert(ngeneratednewcuts != NULL);
1701  assert(cutoff != NULL);
1702 
1703  *cutoff = FALSE;
1704 
1705  if( cutnnz == 0 && SCIPisFeasNegative(scip, cutrhs) ) /*lint !e644*/
1706  {
1707  SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", cutrhs);
1708  *cutoff = TRUE;
1709  return SCIP_OKAY;
1710  }
1711 
1712  /* only take efficient cuts */
1713  if( SCIPisEfficacious(scip, cutefficacy) )
1714  {
1715  SCIP_ROW* cut;
1716  SCIP_VAR* var;
1717  char cutname[SCIP_MAXSTRLEN];
1718  int v;
1719 
1720  /* construct cut name */
1721  if( subgradientiternum >= 0 )
1722  {
1723  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_%d" "_%d" "_%d", SCIPsepaGetName(sepa),
1724  sepadata->ncalls, mainiternum, subgradientiternum, ngeneratedcurrroundcuts + *ngeneratednewcuts);
1725  }
1726  else
1727  {
1728  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_%d" "_%d", SCIPsepaGetName(sepa),
1729  sepadata->ncalls, mainiternum, ngeneratedcurrroundcuts + *ngeneratednewcuts);
1730  }
1731 
1732  /* create empty cut */
1733  SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE,
1734  sepadata->dynamiccuts) );
1735 
1736  /* set cut rank */
1737  SCIProwChgRank(cut, cutrank); /*lint !e644*/
1738 
1739  /* cache the row extension and only flush them if the cut gets added */
1740  SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
1741 
1742  cols = SCIPgetLPCols(scip);
1743 
1744  /* collect all non-zero coefficients */
1745  for( v = 0; v < cutnnz; ++v )
1746  {
1747  var = SCIPcolGetVar(cols[cutinds[v]]);
1748  SCIP_CALL( SCIPaddVarToRow(scip, cut, var, cutcoefs[v]) );
1749  }
1750 
1751  /* flush all changes before adding the cut */
1752  SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
1753 
1754  if( SCIProwGetNNonz(cut) == 0 )
1755  {
1756  assert( SCIPisFeasNegative(scip, cutrhs) );
1757  SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", cutrhs);
1758  *cutoff = TRUE;
1759  return SCIP_OKAY;
1760  }
1761  else
1762  {
1763  /* gathering lhs and rhs again in case the separator is extended later to add other cuts/constraints that may
1764  * have non-inf lhs or inf rhs */
1765  lhs = SCIProwGetLhs(cut);
1766  rhs = SCIProwGetRhs(cut);
1767  assert(SCIPisInfinity(scip, -lhs));
1768  assert(!SCIPisInfinity(scip, rhs));
1769 
1770  SCIPdebugMsg(scip, " -> %s cut <%s>: rhs=%f, eff=%f\n", "lagromory", cutname, cutrhs, cutefficacy);
1771 
1772  /* check if the cut leads to infeasibility (i.e., *cutoff = TRUE) */
1773  {
1774  /* modifiable cuts cannot be declared infeasible, since we don't know all coefficients */
1775  if( SCIProwIsModifiable(cut) )
1776  *cutoff = FALSE;
1777 
1778  /* check for activity infeasibility */
1779  minactivity = SCIPgetRowMinActivity(scip, cut);
1780  maxactivity = SCIPgetRowMaxActivity(scip, cut);
1781 
1782  if( (!SCIPisInfinity(scip, rhs) && SCIPisFeasPositive(scip, minactivity - rhs)) ||
1783  (!SCIPisInfinity(scip, -lhs) && SCIPisFeasNegative(scip, maxactivity - lhs)) )
1784  {
1785  SCIPdebugMsg(scip, "cut <%s> is infeasible (sides=[%g,%g], act=[%g,%g])\n",
1786  SCIProwGetName(cut), lhs, rhs, minactivity, maxactivity);
1787  *cutoff = TRUE;
1788  }
1789  }
1790 
1791  /* store the newly generated cut in an array and update some statistics */
1792  generatedcuts[ngeneratedcurrroundcuts + *ngeneratednewcuts] = cut;
1793  generatedcutefficacies[ngeneratedcurrroundcuts + *ngeneratednewcuts] = cutefficacy;
1794  ++(*ngeneratednewcuts);
1795  }
1796  }
1797 
1798  return SCIP_OKAY;
1799 }
1800 
1801 /** aggregated generated cuts based on the best Lagrangian multipliers */
1802 static
1804  SCIP* scip, /**< SCIP data structure */
1805  SCIP_SEPA* sepa, /**< pointer to the separator */
1806  SCIP_SEPADATA* sepadata, /**< separator data structure */
1807  SCIP_ROW** generatedcuts, /**< cuts generated in the current separation round */
1808  SCIP_Real* bestdualvector, /**< best Lagrangian multipliers vector */
1809  int bestdualvectorlen, /**< length of the best Lagrangian multipliers vector */
1810  SCIP_ROW** aggrcuts, /**< aggregated cuts generated so far in the current separation round */
1811  int* naggrcuts, /**< number of aggregated cuts generated so far in the current separation round */
1812  SCIP_Bool* cutoff /**< should the current node be cutoff? */
1813  )
1814 {
1815  SCIP_Real* cutvals; /**< cut cofficients */
1816  SCIP_COL** cutcols;
1817  SCIP_COL** cols;
1818  SCIP_VAR* var;
1819  SCIP_ROW* cut;
1820  SCIP_ROW* aggrcut;
1821  SCIP_Real minactivity;
1822  SCIP_Real maxactivity;
1823  SCIP_Real cutlhs;
1824  SCIP_Real cutrhs;
1825  SCIP_Real cutconst;
1826  SCIP_Real aggrcutlhs;
1827  SCIP_Real aggrcutrhs;
1828  SCIP_Real aggrcutconst;
1829  SCIP_Real* aggrcutvals;
1830  SCIP_Real* aggrcutcoefs;
1831  SCIP_Real multiplier;
1832  SCIP_Bool aggrcutislocal;
1833  SCIP_Bool aggrindicator;
1834  SCIP_Real* tmpcutvals;
1835  SCIP_Real QUAD(quadterm);
1836  SCIP_Real QUAD(tmpcutconst);
1837  SCIP_Real QUAD(tmpcutrhs);
1838  SCIP_Real QUAD(quadprod);
1839  char aggrcutname[SCIP_MAXSTRLEN];
1840  int cutnnz; /**< number of nonzeros in cut */
1841  int aggrcutnnz;
1842  int* aggrcutinds; /**< column indices in cut */
1843  int aggrcutrank; /**< rank of cut */
1844  int cutrank;
1845  int ncols;
1846  int collppos;
1847  int nlocalcuts;
1848 
1849  assert(scip != NULL);
1850  assert(naggrcuts != NULL);
1851  assert(cutoff != NULL);
1852 
1853  *cutoff = FALSE;
1854  aggrcutlhs = -SCIPinfinity(scip);
1855  aggrcutnnz = 0;
1856  aggrcutrank = -1;
1857  nlocalcuts = 0;
1858  aggrindicator = FALSE;
1859 
1860  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
1861 
1862  /* allocate memory */
1863  SCIP_CALL( SCIPallocCleanBufferArray(scip, &tmpcutvals, QUAD_ARRAY_SIZE(ncols)) );
1864  QUAD_ASSIGN(tmpcutconst, 0.0);
1865  QUAD_ASSIGN(tmpcutrhs, 0.0);
1866  SCIP_CALL( SCIPallocCleanBufferArray(scip, &aggrcutvals, ncols) );
1867  SCIP_CALL( SCIPallocBufferArray(scip, &aggrcutcoefs, ncols) );
1868  SCIP_CALL( SCIPallocBufferArray(scip, &aggrcutinds, ncols) );
1869 
1870  /* aggregate cuts based on the input Lagrangian multipliers */
1871  for( int i = 0; i < bestdualvectorlen; i++ )
1872  {
1873  multiplier = bestdualvector[i];
1874  if( SCIPisGE(scip, multiplier, 1e-4) )
1875  {
1876  cut = generatedcuts[i];
1877  cutnnz = SCIProwGetNNonz(cut);
1878  cutlhs = SCIProwGetLhs(cut);
1879  cutrhs = SCIProwGetRhs(cut);
1880  assert(SCIPisInfinity(scip, -cutlhs));
1881  assert(!SCIPisInfinity(scip, cutrhs));
1882  cutvals = SCIProwGetVals(cut);
1883  cutcols = SCIProwGetCols(cut);
1884  cutconst = SCIProwGetConstant(cut);
1885 
1886  for( int j = 0; j < cutnnz; j++ )
1887  {
1888  collppos = SCIPcolGetLPPos(cutcols[j]);
1889  assert(collppos >= 0);
1890  assert(collppos <= ncols);
1891 
1892  QUAD_ARRAY_LOAD(quadterm, tmpcutvals, collppos);
1893  SCIPquadprecProdDD(quadprod, multiplier, cutvals[j]);
1894  SCIPquadprecSumQQ(quadterm, quadterm, quadprod);
1895  QUAD_ARRAY_STORE(tmpcutvals, collppos, quadterm);
1896  }
1897 
1898  SCIPquadprecProdDD(quadprod, multiplier, cutconst);
1899  SCIPquadprecSumQQ(tmpcutconst, tmpcutconst, quadprod);
1900  SCIPquadprecProdDD(quadprod, multiplier, cutrhs);
1901  SCIPquadprecSumQQ(tmpcutrhs, tmpcutrhs, quadprod);
1902 
1903  cutrank = SCIProwGetRank(cut);
1904  aggrcutrank = MAX(aggrcutrank, cutrank);
1905  nlocalcuts += (SCIProwIsLocal(cut) ? 1 : 0);
1906  aggrindicator = TRUE;
1907  }
1908  }
1909  /* if at least one cut is local, then the aggregated cut is local */
1910  aggrcutislocal = (nlocalcuts > 0 ? TRUE : FALSE);
1911 
1912  /* if the aggregation was successful, then create an empty row and build a cut */
1913  if( aggrindicator )
1914  {
1915  aggrcutconst = QUAD_TO_DBL(tmpcutconst);
1916  aggrcutrhs = QUAD_TO_DBL(tmpcutrhs);
1917 
1918  /* build sparse representation of the aggregated cut */
1919  for( int i = 0; i < ncols; i++ )
1920  {
1921  QUAD_ARRAY_LOAD(quadterm, tmpcutvals, i);
1922  aggrcutvals[i] = QUAD_TO_DBL(quadterm);
1923  if( !SCIPisZero(scip, aggrcutvals[i]) )
1924  {
1925  aggrcutcoefs[aggrcutnnz] = aggrcutvals[i];
1926  aggrcutinds[aggrcutnnz] = i;
1927  aggrcutnnz++;
1928  }
1929  }
1930 
1931  if( aggrcutnnz == 0 && SCIPisFeasNegative(scip, aggrcutrhs) ) /*lint !e644*/
1932  {
1933  SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", aggrcutrhs);
1934  *cutoff = TRUE;
1935  goto TERMINATE;
1936  }
1937 
1938  /* construct aggregated cut name */
1939  (void) SCIPsnprintf(aggrcutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_aggregated", SCIPsepaGetName(sepa),
1940  sepadata->ncalls);
1941 
1942  /* create empty cut */
1943  SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &aggrcut, sepa, aggrcutname, aggrcutlhs - aggrcutconst, aggrcutrhs -
1944  aggrcutconst, aggrcutislocal, FALSE, sepadata->dynamiccuts) );
1945 
1946  /* set cut rank */
1947  SCIProwChgRank(aggrcut, aggrcutrank); /*lint !e644*/
1948 
1949  /* cache the row extension and only flush them if the cut gets added */
1950  SCIP_CALL( SCIPcacheRowExtensions(scip, aggrcut) );
1951 
1952  /* collect all non-zero coefficients */
1953  for( int i = 0; i < aggrcutnnz; i++ )
1954  {
1955  var = SCIPcolGetVar(cols[aggrcutinds[i]]);
1956  SCIP_CALL( SCIPaddVarToRow(scip, aggrcut, var, aggrcutcoefs[i]) );
1957  }
1958 
1959  /* flush all changes before adding the cut */
1960  SCIP_CALL( SCIPflushRowExtensions(scip, aggrcut) );
1961 
1962  if( SCIProwGetNNonz(aggrcut) == 0 )
1963  {
1964  assert( SCIPisFeasNegative(scip, aggrcutrhs) );
1965  SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", aggrcutrhs);
1966  *cutoff = TRUE;
1967  goto TERMINATE;
1968  }
1969  else
1970  {
1971  /* gathering lhs and rhs again in case the separator is extended later to add other cuts/constraints that may
1972  * have non-inf lhs or inf rhs */
1973  cutlhs = SCIProwGetLhs(aggrcut);
1974  cutrhs = SCIProwGetRhs(aggrcut);
1975  assert(SCIPisInfinity(scip, -cutlhs));
1976  assert(!SCIPisInfinity(scip, cutrhs));
1977 
1978  SCIPdebugMsg(scip, " -> %s cut <%s>: rhs=%f\n", "lagromory", aggrcutname, aggrcutrhs);
1979 
1980  /* check if the cut leads to infeasibility (i.e., *cutoff = TRUE) */
1981  {
1982  /* modifiable cuts cannot be declared infeasible, since we don't know all coefficients */
1983  if( SCIProwIsModifiable(aggrcut) )
1984  *cutoff = FALSE;
1985 
1986  /* check for activity infeasibility */
1987  minactivity = SCIPgetRowMinActivity(scip, aggrcut);
1988  maxactivity = SCIPgetRowMaxActivity(scip, aggrcut);
1989 
1990  if( (!SCIPisInfinity(scip, cutrhs) && SCIPisFeasPositive(scip, minactivity - cutrhs)) ||
1991  (!SCIPisInfinity(scip, -cutlhs) && SCIPisFeasNegative(scip, maxactivity - cutlhs)) )
1992  {
1993  SCIPdebugMsg(scip, "cut <%s> is infeasible (sides=[%g,%g], act=[%g,%g])\n",
1994  SCIProwGetName(aggrcut), cutlhs, cutrhs, minactivity, maxactivity);
1995  *cutoff = TRUE;
1996  }
1997  }
1998 
1999  /* add the aggregated cut to a separate data structure */
2000  aggrcuts[*naggrcuts] = aggrcut;
2001  (*naggrcuts)++;
2002  }
2003 
2004  QUAD_ASSIGN(quadterm, 0.0);
2005  for( int i = 0; i < ncols; i++ )
2006  {
2007  aggrcutvals[i] = 0.0;
2008  QUAD_ARRAY_STORE(tmpcutvals, i, quadterm);
2009  }
2010  }
2011 
2012 TERMINATE:
2013  /* free memory */
2014  SCIPfreeBufferArray(scip, &aggrcutinds);
2015  SCIPfreeBufferArray(scip, &aggrcutcoefs);
2016  SCIPfreeCleanBufferArray(scip, &aggrcutvals);
2017  SCIPfreeCleanBufferArray(scip, &tmpcutvals);
2018 
2019  return SCIP_OKAY;
2020 }
2021 
2022 /** main method: LP solution separation method of separator */
2023 static
2025  SCIP* scip, /**< SCIP data structure */
2026  SCIP_SEPA* sepa, /**< pointer to the separator */
2027  SCIP_SEPADATA* sepadata, /**< separator data structure */
2028  int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */
2029  int subgradientiternum, /**< iteration number of the subgradient algorithm */
2030  SCIP_SOL* sol, /**< LP solution to be used for cut generation */
2031  SCIP_Real* solvals, /**< values of the LP solution to be used for cut generation */
2032  int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
2033  SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */
2034  SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */
2035  SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut
2036  generations */
2037  int ngeneratedcurrroundcuts, /**< number of cuts generated until the previous basis in the current separation round */
2038  int* ngeneratednewcuts, /**< number of new cuts generated using the current basis */
2039  int depth, /**< depth of the current node in the tree */
2040  SCIP_Bool* cutoff /**< should the current node be cutoff? */
2041  )
2042 {
2043  SCIP_Real minfrac;
2044  SCIP_Real maxfrac;
2045  SCIP_Real frac;
2046  SCIP_Real* basisfrac;
2047  SCIP_Real* binvrow;
2048  SCIP_AGGRROW* aggrrow;
2049  SCIP_Bool success;
2050  SCIP_Real* cutcoefs;
2051  SCIP_Real cutrhs;
2052  SCIP_Real cutefficacy;
2053  SCIP_Bool cutislocal;
2054  SCIP_ROW** rows;
2055  SCIP_COL** cols;
2056  SCIP_VAR* var;
2057  SCIP_ROW* row;
2058  int cutrank;
2059  int cutnnz;
2060  int* cutinds;
2061  int* basisind;
2062  int* inds;
2063  int nrows;
2064  int ncols;
2065  int* basisperm;
2066  int k;
2067  int c;
2068  int ninds;
2069  int nmaxcutsperlp;
2070 
2071  assert(ngeneratednewcuts != NULL);
2072 
2073  minfrac = sepadata->away;
2074  maxfrac = 1.0 - sepadata->away;
2075  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
2076  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2077  *ngeneratednewcuts = 0;
2078  nmaxcutsperlp = ((depth == 0) ? sepadata->nmaxcutsperlproot : sepadata->nmaxcutsperlp);
2079 
2080  /* allocate memory */
2081  SCIP_CALL( SCIPallocBufferArray(scip, &basisperm, nrows) );
2082  SCIP_CALL( SCIPallocBufferArray(scip, &basisfrac, nrows) );
2083  SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) );
2084  SCIP_CALL( SCIPallocBufferArray(scip, &inds, nrows) );
2085  SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
2086  SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, ncols) );
2087  SCIP_CALL( SCIPallocBufferArray(scip, &cutinds, ncols) );
2088  SCIP_CALL( SCIPaggrRowCreate(scip, &aggrrow) );
2089 
2090  SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
2091 
2092  /* check if the rows of the simplex tableau are suitable for cut generation and build an array of fractions */
2093  for( int i = 0; i < nrows; ++i )
2094  {
2095  frac = 0.0;
2096 
2097  c = basisind[i];
2098 
2099  basisperm[i] = i;
2100 
2101  /* if the simplex tableau row corresponds to an LP column */
2102  if( c >= 0 )
2103  {
2104  assert(c < ncols);
2105 
2106  var = SCIPcolGetVar(cols[c]);
2107  /* if the column is non-continuous one */
2109  {
2110  frac = SCIPfeasFrac(scip, solvals[c]);
2111  frac = MIN(frac, 1.0 - frac);
2112  }
2113  }
2114  /* if the simplex tableau row corresponds to an LP row and separation on rows is allowed */
2115  else if( sepadata->separaterows )
2116  {
2117  assert(0 <= -c-1);
2118  assert(-c-1 < nrows);
2119 
2120  row = rows[-c-1];
2121  /* if the row is suitable for cut generation */
2122  if( SCIProwIsIntegral(row) && !SCIProwIsModifiable(row) )
2123  {
2124  frac = SCIPfeasFrac(scip, SCIPgetRowActivity(scip, row));
2125  frac = MIN(frac, 1.0 - frac);
2126  }
2127  }
2128 
2129  if( frac >= minfrac )
2130  {
2131  /* slightly change fractionality to have random order for equal fractions */
2132  basisfrac[i] = frac + SCIPrandomGetReal(sepadata->randnumgen, -1e-6, 1e-6);
2133  }
2134  else
2135  {
2136  basisfrac[i] = 0.0;
2137  }
2138  }
2139 
2140  /* if there is a need to sort the fractionalities */
2141  if( sepadata->sortcutoffsol )
2142  {
2143  /* sort basis indices by fractionality */
2144  SCIPsortDownRealInt(basisfrac, basisperm, nrows);
2145  }
2146 
2147  /* for all basic columns belonging to integer variables, try to generate a GMI cut */
2148  for( int i = 0; i < nrows && !SCIPisStopped(scip) && !*cutoff; ++i )
2149  {
2150  if( (ngeneratedcurrroundcuts + *ngeneratednewcuts >= nmaxgeneratedperroundcuts) ||
2151  (sepadata->ntotalcuts + *ngeneratednewcuts >= sepadata->nmaxtotalcuts) ||
2152  (*ngeneratednewcuts >= nmaxcutsperlp) )
2153  break;
2154 
2155  ninds = -1;
2156  cutefficacy = 0.0;
2157 
2158  /* either break the loop or proceed to the next iteration if the fractionality is zero */
2159  if( basisfrac[i] == 0.0 )
2160  {
2161  if( sepadata->sortcutoffsol )
2162  break;
2163  else
2164  continue;
2165  }
2166 
2167  k = basisperm[i];
2168 
2169  /* get the row of B^-1 for this basic integer variable with fractional solution value and call aggregate function */
2170  SCIP_CALL( SCIPgetLPBInvRow(scip, k, binvrow, inds, &ninds) );
2171 
2172  SCIP_CALL( SCIPaggrRowSumRows(scip, aggrrow, binvrow, inds, ninds, sepadata->sidetypebasis, allowlocal, 2,
2173  (int) MAXAGGRLEN(ncols), &success) );
2174 
2175  if( !success )
2176  continue;
2177 
2178  /* @todo Currently we are using the SCIPcalcMIR() function to compute the coefficients of the Gomory
2179  * cut. Alternatively, we could use the direct version (see thesis of Achterberg formula (8.4)) which
2180  * leads to cut a of the form \sum a_i x_i \geq 1. Rumor has it that these cuts are better.
2181  */
2182 
2183  /* try to create GMI cut out of the aggregation row */
2185  NULL, minfrac, maxfrac, 1.0, aggrrow, cutcoefs, &cutrhs, cutinds, &cutnnz, &cutefficacy, &cutrank,
2186  &cutislocal, &success) );
2187 
2188  if( success )
2189  {
2190  assert(allowlocal || !cutislocal); /*lint !e644*/
2191 
2192  SCIP_CALL( constructCutRow(scip, sepa, sepadata, mainiternum, subgradientiternum, cutnnz, cutinds, cutcoefs,
2193  cutefficacy, cutrhs, cutislocal, cutrank, generatedcurrroundcuts, generatedcutefficacies,
2194  ngeneratedcurrroundcuts, ngeneratednewcuts, cutoff));
2195  }
2196  }
2197 
2198  /* free temporary memory */
2199  SCIPfreeBufferArray(scip, &cutinds);
2200  SCIPfreeBufferArray(scip, &cutcoefs);
2201  SCIPfreeBufferArray(scip, &basisind);
2202  SCIPfreeBufferArray(scip, &inds);
2203  SCIPfreeBufferArray(scip, &binvrow);
2204  SCIPfreeBufferArray(scip, &basisfrac);
2205  SCIPfreeBufferArray(scip, &basisperm);
2206  SCIPaggrRowFree(scip, &aggrrow);
2207 
2208  return SCIP_OKAY;
2209 }
2210 
2211 /** update objective vector w.r.t. the fixed Lagrangian multipliers */
2212 static
2214  SCIP* scip, /**< SCIP data structure */
2215  SCIP_Real* dualvector, /**< Lagrangian multipliers vector */
2216  SCIP_ROW** cuts, /**< cuts generated so far in the current separation round */
2217  int ncuts, /**< number of cuts generated so far in the current separation round */
2218  SCIP_Real* origobjcoefs, /**< original objective function coefficients of the node linear relaxation */
2219  SCIP_Bool* objvecsdiffer /**< whether the updated objective function coefficients differ from the old ones */
2220  )
2221 {
2222  SCIP_Real* objvals;
2223  SCIP_Real* prod;
2224  SCIP_Real* oldobjvals;
2225  SCIP_Real* cutvals;
2226  SCIP_COL** cutcols;
2227  SCIP_COL** cols;
2228  SCIP_VAR* var;
2229  int cutnnz;
2230  int collppos;
2231  int ncols;
2232 
2233  assert(objvecsdiffer != NULL);
2234  assert(ncuts > 0);
2235 
2236  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2237 
2238  /* allocate memory */
2239  SCIP_CALL( SCIPallocBufferArray(scip, &objvals, ncols) );
2240  SCIP_CALL( SCIPallocBufferArray(scip, &oldobjvals, ncols) );
2241  SCIP_CALL( SCIPallocCleanBufferArray(scip, &prod, ncols) );
2242  *objvecsdiffer = FALSE;
2243 
2244  /* find the product of Lagrangian multipliers and cut coefficients */
2245  for( int i = 0; i < ncuts; i++ )
2246  {
2247  if( !SCIPisZero(scip, dualvector[i]) )
2248  {
2249  assert(!(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) && SCIPisInfinity(scip, SCIProwGetRhs(cuts[i]))));
2250  assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])));
2251  assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])));
2252 
2253  cutnnz = SCIProwGetNNonz(cuts[i]);
2254  assert(cutnnz <= ncols);
2255  cutvals = SCIProwGetVals(cuts[i]);
2256  cutcols = SCIProwGetCols(cuts[i]);
2257 
2258  for( int j = 0; j < cutnnz; ++j )
2259  {
2260  collppos = SCIPcolGetLPPos(cutcols[j]);
2261  assert(collppos >= 0);
2262  assert(collppos <= ncols);
2263 
2264  prod[collppos] += dualvector[i] * cutvals[j];
2265  }
2266  }
2267  }
2268 
2269  /* change objective coefficients */
2270  for( int i = 0; i < ncols; i++ )
2271  {
2272  var = SCIPcolGetVar(cols[i]);
2273  oldobjvals[i] = SCIPgetVarObjDive(scip, var);
2274  objvals[i] = origobjcoefs[i] + prod[i];
2275 
2276  SCIP_CALL( SCIPchgVarObjDive(scip, var, objvals[i]) );
2277 
2278  /* identify if the updated objective vector is indeed different from the previous one */
2279  if( !(*objvecsdiffer) && !SCIPisEQ(scip, oldobjvals[i], objvals[i]) )
2280  *objvecsdiffer = TRUE;
2281  }
2282 
2283  for( int i = 0; i < ncols; i++)
2284  {
2285  prod[i] = 0.0;
2286  }
2287 
2288  /* free memory */
2289  SCIPfreeCleanBufferArray(scip, &prod);
2290  SCIPfreeBufferArray(scip, &oldobjvals);
2291  SCIPfreeBufferArray(scip, &objvals);
2292 
2293  return SCIP_OKAY;
2294 }
2295 
2296 /** add GMI cuts to the objective function of the Lagrangian dual problem by introducing new Lagrangian multipliers */
2297 static
2299  SCIP_Real* dualvector, /**< Lagrangian multipliers vector */
2300  int ngeneratedcuts, /**< number of cuts generated so far in the current separation round */
2301  int* naddedcuts, /**< number of cuts added so far in the current separation round to the Lagrangian dual problem
2302  upon penalization */
2303  int* nnewaddedsoftcuts /**< number of cuts added newly to the Lagrangian dual problem upon penalization */
2304  )
2305 {
2306  assert(*naddedcuts <= ngeneratedcuts);
2307 
2308  /* set the initial penalty of the newly penalized cuts as zero */
2309  for( int i = *naddedcuts; i < ngeneratedcuts; i++ )
2310  dualvector[i] = 0.0;
2311 
2312  *nnewaddedsoftcuts = ngeneratedcuts - *naddedcuts;
2313  *naddedcuts = ngeneratedcuts;
2314 
2315  return SCIP_OKAY;
2316 }
2317 
2318 /** solve the Lagrangian dual problem */
2319 static
2321  SCIP* scip, /**< SCIP data structure */
2322  SCIP_SEPA* sepa, /**< pointer to the separator */
2323  SCIP_SEPADATA* sepadata, /**< separator data structure */
2324  SCIP_SOL* sol, /**< data structure to store an LP solution upon solving a Lagrangian dual problem with
2325  fixed Lagrangian multipliers */
2326  SCIP_Real* solvals, /**< values of the LP solution obtained upon solving a Lagrangian dual problem with fixed
2327  Lagrangian multipliers */
2328  int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */
2329  SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */
2330  int depth, /**< depth of the current node in the tree */
2331  SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */
2332  int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
2333  SCIP_Real* origobjcoefs, /**< original objective function coefficients of the node linear relaxation */
2334  SCIP_Real origobjoffset, /**< original objective function offset of the node linear relaxation */
2335  SCIP_Real* dualvector, /**< Lagrangian multipliers vector */
2336  int* nsoftcuts, /**< number of generated cuts that were penalized and added to the Lagrangian dual problem */
2337  SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */
2338  SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut
2339  generations */
2340  int* ngeneratedcutsperiter, /**< number of cuts generated per subgradient iteration in the current separation round */
2341  int* ngeneratedcurrroundcuts, /**< number of cuts generated so far in the current separation round */
2342  int* ncurrroundlpiters, /**< number of LP iterations taken for solving Lagrangian dual problems with fixed
2343  multipliers in the current separator round */
2344  SCIP_Bool* cutoff, /**< should the current node be cutoff? */
2345  SCIP_Real* bestlagrangianval, /**< best Lagrangian value found so far */
2346  SCIP_Real* bestdualvector, /**< Lagrangian multipliers corresponding to the best Lagrangian value found so far */
2347  int* bestdualvectorlen, /**< length of the Lagrangian multipliers corresponding to the best Lagrangian value
2348  found so far */
2349  int* nbestdualupdates, /**< number of best Lagrangian values found so far */
2350  int* totaliternum /**< total number of iterations of the relax-and-cut algorithm performed so far */
2351  )
2352 {
2353  SCIP_Real* subgradient;
2354  SCIP_Real muparam;
2355  SCIP_Real steplength;
2356  SCIP_Real objval;
2357  SCIP_Real lagrangianval;
2358  SCIP_Real* lagrangianvals;
2359  SCIP_Real avglagrangianval;
2360  SCIP_Real maxsoftcutviol;
2361  SCIP_Real maxnzsubgradientdualprod;
2362  SCIP_Real maxviolscore;
2363  SCIP_Real maxviolscoreold;
2364  SCIP_Real nviolscore;
2365  SCIP_Real nviolscoreold;
2366  SCIP_Real scoreweight;
2367  SCIP_Real ballradius;
2368  SCIP_Bool solfound;
2369  SCIP_Bool backtrack;
2370  SCIP_Bool terminate;
2371  SCIP_Bool subgradientzero;
2372  SCIP_Bool objvecsdiffer;
2373  SCIP_Bool dualvecsdiffer;
2374  SCIP_Bool solvelp;
2375  int ncurrroundlpiterslast;
2376  int nlpiters;
2377  int ngeneratednewcuts;
2378  int nnewaddedsoftcuts;
2379  int nsoftcutviols;
2380  int nnzsubgradientdualprod;
2381 
2382  SCIP_CALL( SCIPallocBufferArray(scip, &lagrangianvals, sepadata->nmaxsubgradientiters) );
2383  SCIP_CALL( SCIPallocCleanBufferArray(scip, &subgradient, nmaxgeneratedperroundcuts) );
2384 
2385  muparam = sepadata->muparaminit;
2386  steplength = 0.0;
2387  objval = 0.0;
2388  avglagrangianval = 0.0;
2389  maxsoftcutviol = 0.0;
2390  maxnzsubgradientdualprod = 0.0;
2391  maxviolscore = 0.0;
2392  nviolscore = 0.0;
2393  scoreweight = 1.0;
2394  ballradius = sepadata->radiusinit;
2395  ngeneratednewcuts = 0;
2396  nsoftcutviols = 0;
2397  nnzsubgradientdualprod = 0;
2398  terminate = FALSE;
2399  subgradientzero = FALSE;
2400  objvecsdiffer = FALSE;
2401  dualvecsdiffer = FALSE;
2402  solvelp = TRUE;
2403 
2404  /* update objective vector based on input Lagrangian multipliers */
2405  if( *nsoftcuts > 0 )
2406  {
2407  SCIP_CALL( updateObjectiveVector(scip, dualvector, generatedcurrroundcuts, *nsoftcuts, origobjcoefs, &objvecsdiffer) );
2408  }
2409 
2410  /* termination check */
2411  SCIP_CALL( checkLagrangianDualTermination(sepadata, -1, -1, FALSE, *ngeneratedcurrroundcuts,
2412  nmaxgeneratedperroundcuts, *ncurrroundlpiters, depth, &terminate) );
2413 
2414  /* the subgradient algorithm loop */
2415  for( int i = 0; i < sepadata->nmaxsubgradientiters && !SCIPisStopped(scip) && !terminate; i++ )
2416  {
2417  solfound = FALSE;
2418  subgradientzero = FALSE;
2419  objvecsdiffer = FALSE;
2420  dualvecsdiffer = FALSE;
2421  nnewaddedsoftcuts = 0;
2422  scoreweight *= sepadata->radiusupdateweight;
2423 
2424  ncurrroundlpiterslast = *ncurrroundlpiters;
2425  if( solvelp )
2426  {
2427  /* solve Lagrangian dual for fixed Lagrangian multipliers */
2428  SCIP_CALL( solveLagromoryLP(scip, sepadata, depth, origobjoffset, &solfound, sol, solvals, &objval,
2429  ncurrroundlpiters) );
2430  }
2431  nlpiters = *ncurrroundlpiters - ncurrroundlpiterslast;
2432 
2433  /* if an optimal solution has been found, then generate cuts and do other operations */
2434  if( solfound )
2435  {
2436  /* generate GMI cuts if a new basis solution is found */
2437  if( (nlpiters >= 1) && (i % sepadata->cutgenfreq == 0) )
2438  {
2439  ngeneratednewcuts = 0;
2440  SCIP_CALL( generateGMICuts(scip, sepa, sepadata, mainiternum, i, sol, solvals,
2441  nmaxgeneratedperroundcuts, allowlocal, generatedcurrroundcuts, generatedcutefficacies,
2442  *ngeneratedcurrroundcuts, &ngeneratednewcuts, depth, cutoff));
2443  sepadata->ntotalcuts += ngeneratednewcuts;
2444  *ngeneratedcurrroundcuts += ngeneratednewcuts;
2445  ngeneratedcutsperiter[mainiternum * sepadata->nmaxsubgradientiters + i + 1] = ngeneratednewcuts;
2446  }
2447 
2448  /* update subgradient, i.e., find the residuals of the penalized cuts, and determine various violations */
2449  updateSubgradient(scip, sol, generatedcurrroundcuts, *nsoftcuts, subgradient, dualvector, &subgradientzero,
2450  &nsoftcutviols, &maxsoftcutviol, &nnzsubgradientdualprod, &maxnzsubgradientdualprod);
2451 
2452  /* calculate Lagrangian value for the fixed Lagrangian multipliers, and update best and avg values */
2453  updateLagrangianValue(scip, objval, dualvector, generatedcurrroundcuts, *nsoftcuts, &lagrangianval);
2454  if( SCIPisPositive(scip, lagrangianval - *bestlagrangianval) )
2455  {
2456  *bestlagrangianval = lagrangianval;
2457  for( int j = 0; j < *nsoftcuts; j++ )
2458  {
2459  bestdualvector[j] = dualvector[j];
2460  }
2461  *bestdualvectorlen = *nsoftcuts;
2462  (*nbestdualupdates)++;
2463  }
2464  lagrangianvals[i] = lagrangianval;
2465  if( i < sepadata->nmaxlagrangianvalsforavg )
2466  {
2467  avglagrangianval = (avglagrangianval * i + lagrangianval)/(i+1);
2468  }
2469  else
2470  {
2471  avglagrangianval = (avglagrangianval * sepadata->nmaxlagrangianvalsforavg -
2472  lagrangianvals[i - sepadata->nmaxlagrangianvalsforavg] +
2473  lagrangianval)/(sepadata->nmaxlagrangianvalsforavg);
2474  }
2475 
2476  /* if the subgradient vector is non-zero, then update the mu parameter and the Lagrangian multipliers */
2477  if( !subgradientzero )
2478  {
2479  /* update mu param */
2480  SCIP_CALL( updateMuSteplengthParam(scip, sepadata, i, ubparam, lagrangianvals, *bestlagrangianval, avglagrangianval,
2481  &muparam, &backtrack) );
2482 
2483  /* update step length */
2484  updateStepLength(scip, muparam, ubparam, lagrangianval, subgradient, *nsoftcuts, &steplength);
2485 
2486  /* update scores to determine how to update the stabilization ball radius */
2487  maxviolscoreold = maxviolscore;
2488  nviolscoreold = nviolscore;
2489  maxviolscore = (1.0 - scoreweight) * maxsoftcutviol + scoreweight * maxnzsubgradientdualprod;
2490  nviolscore = (1.0 - scoreweight) * nsoftcutviols + scoreweight * nnzsubgradientdualprod;
2491 
2492  /* update Lagrangian multipliers */
2493  SCIP_CALL( updateDualVector(scip, sepadata, dualvector, bestdualvector, *bestdualvectorlen,
2494  *nbestdualupdates, i, *totaliternum, steplength, subgradient, *nsoftcuts, backtrack, maxviolscore,
2495  maxviolscoreold, nviolscore, nviolscoreold, nlpiters, &dualvecsdiffer, &ballradius) );
2496 
2497  /* update objective vector based on updated Lagrangian multipliers */
2498  if( dualvecsdiffer )
2499  {
2500  SCIP_CALL( updateObjectiveVector(scip, dualvector, generatedcurrroundcuts, *nsoftcuts, origobjcoefs, &objvecsdiffer) );
2501  }
2502  }
2503  /* if the subgradient vector if zero, then simply mark that the Lagrangian multipliers and the objective
2504  * function of the Lagrangian dual did not change */
2505  else
2506  {
2507  dualvecsdiffer = FALSE;
2508  objvecsdiffer = FALSE;
2509  }
2510 
2511 
2512  /* add generated GMI cuts to the objective function of the Lagrangian dual problem by introducing new
2513  * Lagrangian multipliers */
2514  if( (i % sepadata->cutaddfreq == 0) || (!dualvecsdiffer && !objvecsdiffer &&
2515  (*ngeneratedcurrroundcuts - *nsoftcuts > 0)) )
2516  {
2517  SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) );
2518  }
2519  }
2520  else
2521  {
2522  /* add any remaining generated GMI cuts to the objective function of the Lagrangian dual problem by introducing
2523  * new Lagrangian multipliers */
2524  if( (*ngeneratedcurrroundcuts - *nsoftcuts) > 0 )
2525  {
2526  SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) );
2527  }
2528 
2529  solvelp = FALSE;
2530  }
2531 
2532  /* termination check */
2533  SCIP_CALL( checkLagrangianDualTermination(sepadata, nnewaddedsoftcuts, *ngeneratedcurrroundcuts - *nsoftcuts,
2534  objvecsdiffer, *ngeneratedcurrroundcuts, nmaxgeneratedperroundcuts, *ncurrroundlpiters, depth,
2535  &terminate) );
2536 
2537  (*totaliternum)++;
2538  }
2539 
2540  /* add any remaining generated GMI cuts to the objective function of the Lagrangian dual problem by introducing new
2541  * Lagrangian multipliers */
2542  if( (*ngeneratedcurrroundcuts - *nsoftcuts) > 0 )
2543  {
2544  SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) );
2545  }
2546 
2547  /* free memory */
2548  for( int i = 0; i < nmaxgeneratedperroundcuts; i++)
2549  {
2550  subgradient[i] = 0.0;
2551  }
2552  SCIPfreeCleanBufferArray(scip, &subgradient);
2553  SCIPfreeBufferArray(scip, &lagrangianvals);
2554 
2555  return SCIP_OKAY;
2556 }
2557 
2558 /** generates initial cut pool before solving the Lagrangian dual */
2559 static
2561  SCIP* scip, /**< SCIP data structure */
2562  SCIP_SEPA* sepa, /**< separator */
2563  SCIP_SEPADATA* sepadata, /**< separator data structure */
2564  int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */
2565  SCIP_SOL* sol, /**< LP solution to be used for cut generation */
2566  SCIP_Real* solvals, /**< values of the LP solution to be used for cut generation */
2567  int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
2568  SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */
2569  SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */
2570  SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut
2571  generations */
2572  int* ngeneratedcutsperiter, /**< number of cuts generated per subgradient iteration in the current separation round */
2573  int* ngeneratedcurrroundcuts, /**< number of cuts generated so far in the current separation round */
2574  int depth, /**< depth of the current node in the tree */
2575  SCIP_Bool* cutoff /**< should the current node be cutoff? */
2576  )
2577 {
2578  int ngeneratednewcuts;
2579 
2580  ngeneratednewcuts = 0;
2581 
2582  /* generate initial set of cuts */
2583  SCIP_CALL( generateGMICuts(scip, sepa, sepadata, mainiternum, -1, sol, solvals, nmaxgeneratedperroundcuts,
2584  allowlocal, generatedcurrroundcuts, generatedcutefficacies, *ngeneratedcurrroundcuts, &ngeneratednewcuts,
2585  depth, cutoff) );
2586 
2587  /* update certain statistics */
2588  sepadata->ntotalcuts += ngeneratednewcuts;
2589  *ngeneratedcurrroundcuts += ngeneratednewcuts;
2590  ngeneratedcutsperiter[sepadata->nmaxsubgradientiters * mainiternum] = ngeneratednewcuts;
2591 
2592  return SCIP_OKAY;
2593 }
2594 
2595 /** add cuts to SCIP */
2596 static
2598  SCIP* scip, /**< SCIP data structure */
2599  SCIP_SEPADATA* sepadata, /**< separator data structure */
2600  SCIP_ROW** cuts, /**< cuts generated so far in the current separation round */
2601  int ncuts, /**< number of cuts generated so far in the current separation round */
2602  SCIP_Longint maxdnom, /**< maximum denominator in the rational representation of cuts */
2603  SCIP_Real maxscale, /**< maximal scale factor to scale the cuts to integral values */
2604  int* naddedcuts, /**< number of cuts added to either global cutpool or sepastore */
2605  SCIP_Bool* cutoff /**< should the current node be cutoff? */
2606  )
2607 {
2608  SCIP_ROW* cut;
2609  SCIP_Bool madeintegral;
2610 
2611  assert(scip != NULL);
2612  assert(sepadata != NULL);
2613  assert(naddedcuts != NULL);
2614  assert(cutoff != NULL);
2615 
2616  *cutoff = FALSE;
2617  madeintegral = FALSE;
2618 
2619  for( int i = 0; i < ncuts && !*cutoff; i++ )
2620  {
2621  cut = cuts[i];
2622 
2623  if( SCIProwGetNNonz(cut) == 1 )
2624  {
2625  /* Add the bound change as cut to avoid that the LP gets modified. This would mean that the LP is not flushed
2626  * and the method SCIPgetLPBInvRow() fails; SCIP internally will apply this bound change automatically. */
2627  SCIP_CALL( SCIPaddRow(scip, cut, TRUE, cutoff) );
2628  ++(*naddedcuts);
2629  }
2630  else
2631  {
2632  if( sepadata->makeintegral && SCIPgetRowNumIntCols(scip, cut) == SCIProwGetNNonz(cut) )
2633  {
2634  /* try to scale the cut to integral values */
2635  SCIP_CALL( SCIPmakeRowIntegral(scip, cut, -SCIPepsilon(scip), SCIPsumepsilon(scip),
2636  maxdnom, maxscale, MAKECONTINTEGRAL, &madeintegral) );
2637 
2638  /* if RHS = plus infinity (due to scaling), the cut is useless, so we are not adding it */
2639  if( madeintegral && SCIPisInfinity(scip, SCIProwGetRhs(cut)) )
2640  return SCIP_OKAY;
2641  }
2642 
2643  if( SCIPisCutNew(scip, cut) )
2644  {
2645  /* add global cuts which are not implicit bound changes to the cut pool */
2646  if( !SCIProwIsLocal(cut) )
2647  {
2648  if( sepadata->delayedcuts )
2649  {
2650  SCIP_CALL( SCIPaddDelayedPoolCut(scip, cut) );
2651  }
2652  else
2653  {
2654  SCIP_CALL( SCIPaddPoolCut(scip, cut) );
2655  }
2656  }
2657  else
2658  {
2659  /* local cuts we add to the sepastore */
2660  SCIP_CALL( SCIPaddRow(scip, cut, sepadata->forcecuts, cutoff) );
2661  }
2662  ++(*naddedcuts);
2663  }
2664  }
2665  }
2666 
2667  return SCIP_OKAY;
2668 }
2669 
2670 /** check different termination criteria */
2671 /* @todo nlpssolved criterion? */
2672 static
2674  SCIP_SEPADATA* sepadata, /**< separator data structure */
2675  SCIP_Bool cutoff, /**< should the current node be cutoff? */
2676  SCIP_Bool dualvecsdiffer, /**< whether the updated Lagrangian multipliers differ from the old one */
2677  int ngeneratedcurrroundcuts, /**< number of cuts generated in the current separation round */
2678  int nsoftcuts, /**< number of generated cuts that were penalized and added to the Lagrangian dual problem */
2679  int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
2680  int ncurrroundlpiters, /**< number of LP iterations taken for solving Lagrangian dual problems with fixed
2681  multipliers in the current separator round */
2682  int depth, /**< depth of the current node in the tree */
2683  SCIP_Bool* terminate /**< whether to terminate the relax-and-cut algorithm */
2684  )
2685 {
2686  *terminate = FALSE;
2687 
2688  /* check if the node has been identified to be cutoff */
2689  if( cutoff )
2690  *terminate = TRUE;
2691 
2692  /* check if the Lagrangian multipliers do not differ from the previous iteration and no new cuts exist for penalizing */
2693  if( !dualvecsdiffer && (ngeneratedcurrroundcuts == nsoftcuts) )
2694  *terminate = TRUE;
2695 
2696  /* check if allowed number of cuts in this separation round have already been generated */
2697  if( ngeneratedcurrroundcuts >= nmaxgeneratedperroundcuts )
2698  *terminate = TRUE;
2699 
2700  /* check if allowed number of cuts in the tree have already been generated */
2701  if( sepadata->ntotalcuts >= sepadata->nmaxtotalcuts )
2702  *terminate = TRUE;
2703 
2704  /* check if allowed number of simplex iterations in this separation round have already been used up */
2705  if( (sepadata->nmaxperroundlpiters >= 0) && (ncurrroundlpiters >= sepadata->nmaxperroundlpiters) )
2706  *terminate = TRUE;
2707 
2708  /* check if allowed number of simplex iterations in the root node have already been used up */
2709  if( (depth == 0) && (sepadata->nmaxrootlpiters >= 0) && (sepadata->nrootlpiters >= sepadata->nmaxrootlpiters) )
2710  *terminate = TRUE;
2711 
2712  /* check if allowed number of simplex iterations in the tree have already been used up */
2713  if( (sepadata->nmaxtotallpiters >= 0) && (sepadata->ntotallpiters >= sepadata->nmaxtotallpiters) )
2714  *terminate = TRUE;
2715 
2716  return SCIP_OKAY;
2717 }
2718 
2719 /** Searches and tries to add Lagromory cuts */
2720 static
2722  SCIP* scip, /**< SCIP data structure */
2723  SCIP_SEPA* sepa, /**< separator */
2724  SCIP_SEPADATA* sepadata, /**< separator data structure */
2725  SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */
2726  int depth, /**< depth of the current node in the tree */
2727  SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */
2728  SCIP_RESULT* result /**< final result of the separation round */
2729  )
2730 {
2731  SCIP_ROW** generatedcurrroundcuts;
2732  SCIP_ROW** aggregatedcurrroundcuts;
2733  SCIP_Real* generatedcutefficacies;
2734  SCIP_Bool solfound;
2735  SCIP_SOL* softcutslpsol;
2736  SCIP_Real* softcutslpsolvals;
2737  SCIP_SOL* hardcutslpsol;
2738  SCIP_Real* hardcutslpsolvals;
2739  SCIP_Real* dualsol;
2740  SCIP_Real* dualvector;
2741  SCIP_Real* bestdualvector;
2742  SCIP_Real bestlagrangianval;
2743  SCIP_Real* origobjcoefs;
2744  SCIP_Real origobjoffset;
2745  SCIP_Real objval;
2746  SCIP_Real maxscale;
2747  SCIP_Longint maxdnom;
2748  SCIP_Bool cutoff;
2749  SCIP_Bool cutoff2;
2750  SCIP_Bool dualvecsdiffer;
2751  SCIP_Bool terminate;
2752  SCIP_COL** cols;
2753 
2754  int* ngeneratedcutsperiter;
2755  int bestdualvectorlen;
2756  int nbestdualupdates;
2757  int totaliternum;
2758  int* cutindsperm;
2759  int nprocessedcuts;
2760  int ncols;
2761  int nrows;
2762  int ngeneratedcurrroundcuts;
2763  int nselectedcurrroundcuts;
2764  int nselectedcuts;
2765  int naddedcurrroundcuts;
2766  int naggregatedcurrroundcuts;
2767  int nmaxgeneratedperroundcuts;
2768  int ncurrroundlpiters;
2769  int nsoftcuts;
2770  int nsoftcutsold;
2771  int maxdepth;
2772 
2773  assert(*result == SCIP_DIDNOTRUN);
2774  assert(sepadata != NULL);
2775  sepadata->ncalls = SCIPsepaGetNCalls(sepa);
2776  objval = SCIPinfinity(scip);
2777  bestlagrangianval = -SCIPinfinity(scip);
2778  bestdualvectorlen = 0;
2779  nbestdualupdates = 0;
2780  totaliternum = 0;
2781  ngeneratedcurrroundcuts = 0;
2782  naddedcurrroundcuts = 0;
2783  naggregatedcurrroundcuts = 0;
2784  ncurrroundlpiters = 0;
2785  nsoftcuts = 0;
2786  solfound = FALSE;
2787  cutoff = FALSE;
2788  cutoff2 = FALSE;
2789  dualvecsdiffer = FALSE;
2790  terminate = FALSE;
2791 
2792  SCIPdebugMsg(scip, "Separating cuts...\n");
2793 
2794  /* initialize the LP that will be used to solve the Lagrangian dual with fixed Lagrangian multipliers */
2795  SCIP_CALL( createLPWithSoftCuts(scip, sepadata) );
2796  /* create the LP that represents the node relaxation including all the generated cuts in this separator */
2797  SCIP_CALL( createLPWithHardCuts(scip, sepadata, NULL, 0) );
2798 
2799  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2800  nrows = SCIPgetNLPRows(scip);
2801 
2802  /* get the maximal number of cuts allowed in a separation round */
2803  nmaxgeneratedperroundcuts = ((depth == 0) ? sepadata->nmaxperroundcutsroot : sepadata->nmaxperroundcuts);
2804 
2805  /* allocate memory */
2806  SCIP_CALL( SCIPallocBufferArray(scip, &generatedcurrroundcuts, nmaxgeneratedperroundcuts) );
2807  SCIP_CALL( SCIPallocBufferArray(scip, &aggregatedcurrroundcuts, nmaxgeneratedperroundcuts) );
2808  SCIP_CALL( SCIPallocBufferArray(scip, &generatedcutefficacies, nmaxgeneratedperroundcuts) );
2809  SCIP_CALL( SCIPallocBufferArray(scip, &cutindsperm, nmaxgeneratedperroundcuts) );
2810  SCIP_CALL( SCIPallocCleanBufferArray(scip, &ngeneratedcutsperiter, sepadata->nmaxmainiters *
2811  (sepadata->nmaxsubgradientiters + 1)) );
2812  SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualsol, nrows + nmaxgeneratedperroundcuts) );
2813  SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualvector, nmaxgeneratedperroundcuts) );
2814  SCIP_CALL( SCIPallocCleanBufferArray(scip, &bestdualvector, nmaxgeneratedperroundcuts) );
2815  SCIP_CALL( SCIPallocBufferArray(scip, &softcutslpsolvals, ncols) );
2816  SCIP_CALL( SCIPcreateSol(scip, &softcutslpsol, NULL) );
2817  SCIP_CALL( SCIPallocBufferArray(scip, &hardcutslpsolvals, ncols) );
2818  SCIP_CALL( SCIPcreateSol(scip, &hardcutslpsol, NULL) );
2819  SCIP_CALL( SCIPallocBufferArray(scip, &origobjcoefs, ncols) );
2820 
2821  /* store current objective function */
2822  for( int i = 0; i < ncols; i++ )
2823  {
2824  origobjcoefs[i] = SCIPcolGetObj(cols[i]);
2825  }
2826  origobjoffset = SCIPgetTransObjoffset(scip);
2827 
2828  /* solve node LP relaxation to have an initial simplex tableau */
2829  SCIP_CALL( solveLagromoryLP(scip, sepadata, depth, origobjoffset, &solfound, softcutslpsol, softcutslpsolvals, &objval,
2830  &ncurrroundlpiters));
2831 
2832  /* generate initial cut pool */
2833  SCIP_CALL( generateInitCutPool(scip, sepa, sepadata, 0, softcutslpsol, softcutslpsolvals, nmaxgeneratedperroundcuts, allowlocal,
2834  generatedcurrroundcuts, generatedcutefficacies, ngeneratedcutsperiter, &ngeneratedcurrroundcuts, depth,
2835  &cutoff) );
2836 
2837  /* termination check */
2838  SCIP_CALL( checkMainLoopTermination(sepadata, cutoff, TRUE, ngeneratedcurrroundcuts, nsoftcuts,
2839  nmaxgeneratedperroundcuts, ncurrroundlpiters, depth, &terminate) );
2840 
2841  /* compute cuts for each integer col with fractional val */
2842  for( int i = 0; i < sepadata->nmaxmainiters && !SCIPisStopped(scip) && !terminate; ++i )
2843  {
2844  nsoftcutsold = nsoftcuts;
2845  nsoftcuts = ngeneratedcurrroundcuts;
2846  dualvecsdiffer = FALSE;
2847 
2848  /* solve Lagrangain dual */
2849  SCIP_CALL( solveLagrangianDual(scip, sepa, sepadata, softcutslpsol, softcutslpsolvals, i, ubparam, depth, allowlocal,
2850  nmaxgeneratedperroundcuts, origobjcoefs, origobjoffset, dualvector, &nsoftcuts, generatedcurrroundcuts,
2851  generatedcutefficacies, ngeneratedcutsperiter, &ngeneratedcurrroundcuts, &ncurrroundlpiters, &cutoff,
2852  &bestlagrangianval, bestdualvector, &bestdualvectorlen, &nbestdualupdates, &totaliternum) );
2853 
2854  /* @todo filter cuts before adding them to the new LP that was created based on the node relaxation? */
2855 
2856  /* update the LP representing the node relaxation by adding newly generated cuts */
2857  if( !cutoff && (ngeneratedcurrroundcuts - nsoftcutsold > 0) )
2858  {
2859  SCIP_CALL( createLPWithHardCuts(scip, sepadata, generatedcurrroundcuts, ngeneratedcurrroundcuts) );
2860 
2861  /* solve the LP and get relevant information */
2862  SCIP_CALL( solveLPWithHardCuts(scip, sepadata, &solfound, hardcutslpsol, hardcutslpsolvals) );
2863 
2864  /* if primal solution is found, then pass it to trysol heuristic */
2865  /* @note if trysol heuristic is was not present, then the solution found above, which can potentially be a MIP
2866  * primal feasible solution, will go to waste.
2867  */
2868  if( solfound && sepadata->heurtrysol != NULL )
2869  {
2870  SCIP_CALL( SCIPheurPassSolTrySol(scip, sepadata->heurtrysol, hardcutslpsol) );
2871  }
2872 
2873  /* if dual feasible, then fetch dual solution and reset Lagrangian multipliers based on it. otherwise, retain the
2874  * Lagrangian multipliers and simply initialize the new multipliers to zeroes. */
2875  if( SCIPlpiIsDualFeasible(sepadata->lpiwithhardcuts) )
2876  {
2877  SCIP_CALL( SCIPlpiGetSol(sepadata->lpiwithhardcuts, NULL, NULL, dualsol, NULL, NULL) );
2878  SCIP_CALL( updateDualVector(scip, sepadata, dualvector, &(dualsol[nrows]),
2879  ngeneratedcurrroundcuts, 0, -1, -1, 0.0, NULL, ngeneratedcurrroundcuts, TRUE, 0.0, 0.0, 0.0, 0.0, -1,
2880  &dualvecsdiffer, NULL) );
2881  }
2882  else
2883  {
2884  SCIP_CALL( updateDualVector(scip, sepadata, dualvector, dualvector, nsoftcuts, 0, -1, -1, 0.0, NULL,
2885  ngeneratedcurrroundcuts, TRUE, 0.0, 0.0, 0.0, 0.0, -1, &dualvecsdiffer, NULL) );
2886  }
2887  }
2888 
2889  /* termination check */
2890  SCIP_CALL( checkMainLoopTermination(sepadata, cutoff, dualvecsdiffer, ngeneratedcurrroundcuts, nsoftcuts,
2891  nmaxgeneratedperroundcuts, ncurrroundlpiters, depth, &terminate) );
2892  }
2893 
2894  /* set the maximal denominator in rational representation of gomory cut and the maximal scale factor to
2895  * scale resulting cut to integral values to avoid numerical instabilities
2896  */
2897  /**@todo find better but still stable gomory cut settings: look at dcmulti, gesa3, khb0525, misc06, p2756 */
2898  /* @note: above todo was copied from sepa_gomory.c. So, if gomory code is changed, same changes can be done here. */
2899  maxdepth = SCIPgetMaxDepth(scip);
2900  if( depth == 0 )
2901  {
2902  maxdnom = 1000;
2903  maxscale = 1000.0;
2904  }
2905  else if( depth <= maxdepth/4 )
2906  {
2907  maxdnom = 1000;
2908  maxscale = 1000.0;
2909  }
2910  else if( depth <= maxdepth/2 )
2911  {
2912  maxdnom = 100;
2913  maxscale = 100.0;
2914  }
2915  else
2916  {
2917  maxdnom = 10;
2918  maxscale = 10.0;
2919  }
2920 
2921  /* filter cuts by calling cut selection algorithms and add cuts accordingly */
2922  /* @todo an idea is to filter cuts after every main iter */
2923  /* @todo we can remove !cutoff criterion by adding forcedcuts */
2924  if( !cutoff && (ngeneratedcurrroundcuts > 0) )
2925  {
2926  if( SCIPisGE(scip, sepadata->cutsfilterfactor, 1.0) )
2927  {
2928  nselectedcurrroundcuts = ngeneratedcurrroundcuts;
2929  SCIP_CALL( addCuts(scip, sepadata, generatedcurrroundcuts, nselectedcurrroundcuts, maxdnom, maxscale,
2930  &naddedcurrroundcuts, &cutoff2) );
2931  cutoff = cutoff2;
2932  }
2933  else if( SCIPisPositive(scip, sepadata->cutsfilterfactor) )
2934  {
2935  nprocessedcuts = 0;
2936  for( int i = 0; i < sepadata->nmaxmainiters * (sepadata->nmaxsubgradientiters + 1); i++ )
2937  {
2938  if( ngeneratedcutsperiter[i] != 0 )
2939  {
2940  for( int j = 0; j < ngeneratedcutsperiter[i]; j++ )
2941  cutindsperm[j] = j + nprocessedcuts;
2942 
2943  /* sort cut efficacies by fractionality */
2944  SCIPsortDownRealInt(&generatedcutefficacies[nprocessedcuts], cutindsperm, ngeneratedcutsperiter[i]);
2945 
2946  nselectedcuts = (int)SCIPceil(scip, sepadata->cutsfilterfactor * ngeneratedcutsperiter[i]);
2947 
2948  SCIP_CALL( addCuts(scip, sepadata, &generatedcurrroundcuts[nprocessedcuts], nselectedcuts, maxdnom,
2949  maxscale, &naddedcurrroundcuts, &cutoff2) );
2950  cutoff = cutoff2;
2951 
2952  nprocessedcuts += ngeneratedcutsperiter[i];
2953  }
2954  }
2955  }
2956  }
2957  else if( ngeneratedcurrroundcuts > 0 )
2958  {
2959  nselectedcurrroundcuts = ngeneratedcurrroundcuts;
2960  SCIP_CALL( addCuts(scip, sepadata, generatedcurrroundcuts, nselectedcurrroundcuts, maxdnom, maxscale,
2961  &naddedcurrroundcuts, &cutoff2) );
2962  }
2963 
2964  /* add an aggregated cut based on best Lagrangian multipliers */
2965  if( sepadata->aggregatecuts && (ngeneratedcurrroundcuts > 0) && (bestdualvectorlen > 0) )
2966  {
2967  assert(bestdualvectorlen <= ngeneratedcurrroundcuts);
2968  SCIP_CALL( aggregateGeneratedCuts(scip, sepa, sepadata, generatedcurrroundcuts, bestdualvector, bestdualvectorlen,
2969  aggregatedcurrroundcuts, &naggregatedcurrroundcuts, &cutoff2) );
2970  cutoff = (!cutoff ? cutoff2 : cutoff);
2971  if( naggregatedcurrroundcuts > 0 )
2972  {
2973  SCIP_CALL( addCuts(scip, sepadata, aggregatedcurrroundcuts, naggregatedcurrroundcuts, maxdnom, maxscale,
2974  &naddedcurrroundcuts, &cutoff2) );
2975  cutoff = (!cutoff ? cutoff2 : cutoff);
2976  }
2977  }
2978 
2979  if( cutoff )
2980  {
2981  *result = SCIP_CUTOFF;
2982  }
2983  else if( naddedcurrroundcuts > 0 )
2984  {
2985  *result = SCIP_SEPARATED;
2986  }
2987  else
2988  {
2989  *result = SCIP_DIDNOTFIND;
2990  }
2991 
2992  /* delete the LP representing the Lagrangian dual problem with fixed Lagrangian multipliers */
2993  SCIP_CALL( deleteLPWithSoftCuts(scip, sepadata) );
2994 
2995  /* release the rows in aggregatedcurrroundcuts */
2996  for( int i = 0; i < naggregatedcurrroundcuts; ++i )
2997  {
2998  SCIP_CALL( SCIPreleaseRow(scip, &(aggregatedcurrroundcuts[i])) );
2999  }
3000  /* release the rows in generatedcurrroundcuts */
3001  for( int i = 0; i < ngeneratedcurrroundcuts; ++i )
3002  {
3003  SCIP_CALL( SCIPreleaseRow(scip, &(generatedcurrroundcuts[i])) );
3004  }
3005 
3006  for( int i = 0; i < sepadata->nmaxmainiters * (sepadata->nmaxsubgradientiters + 1); i++ )
3007  {
3008  ngeneratedcutsperiter[i] = 0;
3009  }
3010  for( int i = 0; i < nrows; i++ )
3011  {
3012  dualsol[i] = 0.0;
3013  }
3014  for( int i = 0; i < nmaxgeneratedperroundcuts; i++ )
3015  {
3016  dualsol[nrows + i] = 0.0;
3017  dualvector[i] = 0.0;
3018  bestdualvector[i] = 0.0;
3019  }
3020  /* free memory */
3021  SCIPfreeBufferArray(scip, &origobjcoefs);
3022  SCIPfreeBufferArray(scip, &hardcutslpsolvals);
3023  SCIP_CALL( SCIPfreeSol(scip, &hardcutslpsol) );
3024  SCIPfreeBufferArray(scip, &softcutslpsolvals);
3025  SCIP_CALL( SCIPfreeSol(scip, &softcutslpsol) );
3026  SCIPfreeCleanBufferArray(scip, &ngeneratedcutsperiter);
3027  SCIPfreeBufferArray(scip, &cutindsperm);
3028  SCIPfreeBufferArray(scip, &generatedcutefficacies);
3029  SCIPfreeBufferArray(scip, &aggregatedcurrroundcuts);
3030  SCIPfreeBufferArray(scip, &generatedcurrroundcuts);
3031  SCIPfreeCleanBufferArray(scip, &dualvector);
3032  SCIPfreeCleanBufferArray(scip, &bestdualvector);
3033  SCIPfreeCleanBufferArray(scip, &dualsol);
3034 
3035  return SCIP_OKAY;
3036 }
3037 
3038 /** creates separator data */
3039 static
3041  SCIP* scip, /**< SCIP data structure */
3042  SCIP_SEPADATA** sepadata /**< separator data structure */
3043  )
3044 {
3045  assert(scip != NULL);
3046  assert(sepadata != NULL);
3047 
3048  SCIP_CALL( SCIPallocBlockMemory(scip, sepadata) );
3049  BMSclearMemory(*sepadata);
3050 
3051  return SCIP_OKAY;
3052 }
3053 
3054 
3055 /*
3056  * Callback methods of separator
3057  */
3058 
3059 /** copy method for separator plugins (called when SCIP copies plugins) */
3060 static
3061 SCIP_DECL_SEPACOPY(sepaCopyLagromory)
3062 { /*lint --e{715}*/
3063  assert(scip != NULL);
3064  assert(sepa != NULL);
3065  assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3066 
3067  /* call inclusion method of constraint handler */
3069 
3070  return SCIP_OKAY;
3071 }
3072 
3073 
3074 /** destructor of separator to free user data (called when SCIP is exiting) */
3075 static
3076 SCIP_DECL_SEPAFREE(sepaFreeLagromory)
3077 {
3078  SCIP_SEPADATA* sepadata;
3079 
3080  sepadata = SCIPsepaGetData(sepa);
3081  assert(sepadata != NULL);
3082 
3083  SCIP_CALL( sepadataFree(scip, &sepadata) );
3084  SCIPsepaSetData(sepa, NULL);
3085 
3086  return SCIP_OKAY;
3087 }
3088 
3089 
3090 /** initialization method of separator (called after problem was transformed) */
3091 static
3092 SCIP_DECL_SEPAINIT(sepaInitLagromory)
3093 {
3094  SCIP_SEPADATA* sepadata;
3095 
3096  sepadata = SCIPsepaGetData(sepa);
3097  assert(scip != NULL);
3098  assert(sepadata != NULL);
3099 
3100  /* create and initialize random number generator */
3101  SCIP_CALL( SCIPcreateRandom(scip, &(sepadata->randnumgen), RANDSEED, TRUE) );
3102 
3103  /* find trysol heuristic */
3104  if ( sepadata->heurtrysol == NULL )
3105  {
3106  sepadata->heurtrysol = SCIPfindHeur(scip, "trysol");
3107  }
3108 
3109  return SCIP_OKAY;
3110 }
3111 
3112 /** deinitialization method of separator (called before transformed problem is freed) */
3113 static
3114 SCIP_DECL_SEPAEXIT(sepaExitLagromory)
3115 { /*lint --e{715}*/
3116  SCIP_SEPADATA* sepadata;
3117 
3118  sepadata = SCIPsepaGetData(sepa);
3119  assert(sepadata != NULL);
3120 
3121  SCIPfreeRandom(scip, &(sepadata->randnumgen));
3122 
3123  return SCIP_OKAY;
3124 }
3125 
3126 /** LP solution separation method of separator */
3127 static
3128 SCIP_DECL_SEPAEXECLP(sepaExeclpLagromory)
3129 {
3130  /*lint --e{715}*/
3131 
3132  SCIP_SEPADATA* sepadata;
3133  SCIP_SOL* bestsol;
3134  SCIP_COL** cols;
3135  SCIP_COL* col;
3136  SCIP_VAR* var;
3137  SCIP_Real ubparam;
3138  SCIP_Real dualdegeneracyrate;
3139  SCIP_Real varconsratio;
3140  SCIP_Real threshold1;
3141  SCIP_Real threshold2;
3142  int nrows;
3143  int ncols;
3144  int ncalls;
3145  int runnum;
3146 
3147  assert(sepa != NULL);
3148  assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3149  sepadata = SCIPsepaGetData(sepa);
3150  assert(sepadata != NULL);
3151 
3152  assert(result != NULL);
3153  *result = SCIP_DIDNOTRUN;
3154 
3155  assert(scip != NULL);
3156 
3157  ncalls = SCIPsepaGetNCallsAtNode(sepa);
3158  runnum = SCIPgetNRuns(scip);
3159  dualdegeneracyrate = 0.0;
3160  varconsratio = 0.0;
3161 
3162  /* only generate Lagromory cuts if we are not close to terminating */
3163  if( SCIPisStopped(scip) )
3164  return SCIP_OKAY;
3165 
3166  /* only call the separator starting sepadata->minrestart runs */
3167  if( (sepadata->minrestart >= 1) && (runnum < sepadata->minrestart + 1) )
3168  return SCIP_OKAY;
3169 
3170  /* only call the separator a given number of times at each node */
3171  if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
3172  || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
3173  return SCIP_OKAY;
3174 
3175  /* only generate cuts if an optimal LP solution is at hand */
3177  return SCIP_OKAY;
3178 
3179  /* only generate cuts if the LP solution is basic */
3180  if( ! SCIPisLPSolBasic(scip) )
3181  return SCIP_OKAY;
3182 
3183  /* get LP data */
3184  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
3185  assert(cols != NULL);
3186 
3187  nrows = SCIPgetNLPRows(scip);
3188 
3189  /* return if LP has no columns or no rows */
3190  if( ncols == 0 || nrows == 0 )
3191  return SCIP_OKAY;
3192 
3193  /* return if dual degeneracy metrics are below threshold values */
3194  threshold1 = sepadata->dualdegeneracyratethreshold;
3195  threshold2 = sepadata->varconsratiothreshold;
3196  SCIP_CALL( SCIPgetLPDualDegeneracy(scip, &dualdegeneracyrate, &varconsratio) );
3197  if( (dualdegeneracyrate < threshold1) && (varconsratio < threshold2) )
3198  return SCIP_OKAY;
3199 
3200  bestsol = SCIPgetBestSol(scip);
3201  ubparam = 0.0;
3202 
3203  /* if a MIP primal solution exists, use it to estimate the optimal value of the Lagrangian dual problem */
3204  if( bestsol != NULL )
3205  {
3206  for( int i = 0; i < ncols; ++i )
3207  {
3208  col = cols[i];
3209  assert(col != NULL);
3210 
3211  var = SCIPcolGetVar(col);
3212 
3213  ubparam += SCIPgetSolVal(scip, bestsol, var) * SCIPcolGetObj(col);
3214  }
3215 
3216  ubparam += SCIPgetTransObjoffset(scip);
3217  }
3218  /* if a MIP primal solution does not exist, then use the node relaxation's LP solutin to estimate the optimal value
3219  * of the Lagrangian dual problem
3220  */
3221  else
3222  {
3223  for( int i = 0; i < ncols; ++i )
3224  {
3225  col = cols[i];
3226  assert(col != NULL);
3227 
3228  ubparam += SCIPcolGetPrimsol(col) * SCIPcolGetObj(col);
3229  }
3230 
3231  ubparam += SCIPgetTransObjoffset(scip);
3232  ubparam *= SCIPisPositive(scip, ubparam) ? sepadata->ubparamposfactor : sepadata->ubparamnegfactor;
3233  }
3234 
3235  /* the main separation function of the separator */
3236  SCIP_CALL( separateCuts(scip, sepa, sepadata, ubparam, depth, (allowlocal && sepadata->allowlocal), result) );
3237 
3238  return SCIP_OKAY;
3239 }
3240 
3241 
3242 /*
3243  * separator specific interface methods
3244  */
3245 
3246 /** creates the Lagromory separator and includes it in SCIP */
3248  SCIP* scip /**< SCIP data structure */
3249  )
3250 {
3251  SCIP_SEPADATA* sepadata;
3252  SCIP_SEPA* sepa;
3253 
3254  /* create separator data */
3255  SCIP_CALL( sepadataCreate(scip, &sepadata) );
3256 
3257  sepadata->heurtrysol = NULL;
3258 
3259  /* include separator */
3261  SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpLagromory, NULL, sepadata) );
3262 
3263  assert(sepa != NULL);
3264 
3265  /* set non fundamental callbacks via setter functions */
3266  SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyLagromory) );
3267  SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeLagromory) );
3268  SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitLagromory) );
3269  SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitLagromory) );
3270 
3271  /* add separator parameters */
3272  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/away", "minimal integrality violation of a basis "
3273  "variable to try separation", &sepadata->away, FALSE, DEFAULT_AWAY, 0.0, 1.0, NULL, NULL) );
3274 
3275  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/rootlpiterlimitfactor", "factor w.r.t. root node LP "
3276  "iterations for maximal separating LP iterations in the root node (negative for no limit)",
3277  &sepadata->rootlpiterlimitfactor, TRUE, DEFAULT_ROOTLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) );
3278 
3279  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/totallpiterlimitfactor", "factor w.r.t. root node LP "
3280  "iterations for maximal separating LP iterations in the tree (negative for no limit)",
3281  &sepadata->totallpiterlimitfactor, TRUE, DEFAULT_TOTALLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) );
3282 
3283  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundlpiterlimitfactor", "factor w.r.t. root node LP "
3284  "iterations for maximal separating LP iterations per separation round (negative for no limit)",
3285  &sepadata->perroundlpiterlimitfactor, TRUE, DEFAULT_PERROUNDLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL,
3286  NULL) );
3287 
3288  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundcutsfactorroot", "factor w.r.t. number of integer "
3289  "columns for number of cuts separated per separation round in root node", &sepadata->perroundcutsfactorroot,
3291 
3292  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundcutsfactor", "factor w.r.t. number of integer "
3293  "columns for number of cuts separated per separation round at a non-root node",
3294  &sepadata->perroundcutsfactor, TRUE, DEFAULT_PERROUNDCUTSFACTOR, 0.0, SCIP_REAL_MAX, NULL, NULL) );
3295 
3296  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/totalcutsfactor", "factor w.r.t. number of integer "
3297  "columns for total number of cuts separated", &sepadata->totalcutsfactor, TRUE, DEFAULT_TOTALCUTSFACTOR, 0.0,
3298  SCIP_REAL_MAX, NULL, NULL) );
3299 
3300  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparaminit", "initial value of the mu parameter (factor "
3301  "for step length)", &sepadata->muparaminit, TRUE, DEFAULT_MUPARAMINIT, 0.0, 100.0, NULL, NULL) );
3302 
3303  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparamlb", "lower bound of the mu parameter (factor for "
3304  "step length)", &sepadata->muparamlb, TRUE, DEFAULT_MUPARAMLB, 0.0, 1.0, NULL, NULL) );
3305 
3306  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparamub", "upper bound of the mu parameter (factor for "
3307  "step length)", &sepadata->muparamub, TRUE, DEFAULT_MUPARAMUB, 1.0, 10.0, NULL, NULL) );
3308 
3309  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/mubacktrackfactor", "factor of mu while backtracking the "
3310  "mu parameter (factor for step length)", &sepadata->mubacktrackfactor, TRUE, DEFAULT_MUBACKTRACKFACTOR,
3311  0.0, 1.0, NULL, NULL) );
3312 
3313  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab1factor", "factor of mu parameter (factor for step "
3314  "length) for larger increment" , &sepadata->muslab1factor, TRUE, DEFAULT_MUSLAB1FACTOR, 0.0, SCIP_REAL_MAX, NULL,
3315  NULL));
3316 
3317  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab2factor", "factor of mu parameter (factor for step "
3318  "length) for smaller increment", &sepadata->muslab2factor, TRUE, DEFAULT_MUSLAB2FACTOR, 0.0, SCIP_REAL_MAX, NULL,
3319  NULL) );
3320 
3321  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab3factor", "factor of mu parameter (factor for step "
3322  "length) for reduction", &sepadata->muslab3factor, TRUE, DEFAULT_MUSLAB3FACTOR, 0.0, SCIP_REAL_MAX, NULL, NULL) );
3323 
3324  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/deltaslab1ub", "factor of delta deciding larger increment "
3325  "of mu parameter (factor for step length)", &sepadata->deltaslab1ub, TRUE, DEFAULT_DELTASLAB1UB, 0.0, 1.0,
3326  NULL, NULL) );
3327 
3328  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/deltaslab2ub", "factor of delta deciding smaller "
3329  "increment of mu parameter (factor for step length)", &sepadata->deltaslab2ub, TRUE, DEFAULT_DELTASLAB2UB,
3330  0.0, 1.0, NULL, NULL) );
3331 
3332  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/ubparamposfactor", "factor for positive upper bound used "
3333  "as an estimate for the optimal Lagrangian dual value", &sepadata->ubparamposfactor, TRUE, DEFAULT_UBPARAMPOSFACTOR,
3334  1.0, 100.0, NULL, NULL) );
3335 
3336  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/ubparamnegfactor", "factor for negative upper bound used "
3337  "as an estimate for the optimal Lagrangian dual value", &sepadata->ubparamnegfactor, TRUE, DEFAULT_UBPARAMNEGFACTOR,
3338  0.0, 1.0, NULL, NULL) );
3339 
3340  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perrootlpiterfactor", "factor w.r.t. root node LP "
3341  "iterations for iteration limit of each separating LP (negative for no limit)",
3342  &sepadata->perrootlpiterfactor, TRUE, DEFAULT_PERROOTLPITERFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) );
3343 
3344  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perlpiterfactor", "factor w.r.t. non-root node LP "
3345  "iterations for iteration limit of each separating LP (negative for no limit)", &sepadata->perlpiterfactor, TRUE,
3347 
3348  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/cutsfilterfactor", "fraction of generated cuts per "
3349  "explored basis to accept from separator", &sepadata->cutsfilterfactor, TRUE, DEFAULT_CUTSFILTERFACTOR, 0.0,
3350  1.0, NULL, NULL) );
3351 
3352  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusinit", "initial radius of the ball used in "
3353  "stabilization of Lagrangian multipliers", &sepadata->radiusinit, TRUE, DEFAULT_RADIUSINIT, 0.0, 1.0, NULL,
3354  NULL) );
3355 
3356  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusmax", "maximum radius of the ball used in "
3357  "stabilization of Lagrangian multipliers", &sepadata->radiusmax, TRUE, DEFAULT_RADIUSMAX, 0.0, SCIP_REAL_MAX,
3358  NULL, NULL) );
3359 
3360  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusmin", "minimum radius of the ball used in "
3361  "stabilization of Lagrangian multipliers", &sepadata->radiusmin, TRUE, DEFAULT_RADIUSMIN, 0.0, SCIP_REAL_MAX,
3362  NULL, NULL) );
3363 
3364  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/constant", "a constant for stablity center based "
3365  "stabilization of Lagrangian multipliers", &sepadata->constant, TRUE, DEFAULT_CONST, 2.0, SCIP_REAL_MAX,
3366  NULL, NULL) );
3367 
3368  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusupdateweight", "multiplier to evaluate cut "
3369  "violation score used for updating ball radius", &sepadata->radiusupdateweight, TRUE,
3370  DEFAULT_RADIUSUPDATEWEIGHT, 0.0, 1.0, NULL, NULL) );
3371 
3372  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/dualdegeneracyratethreshold", "minimum dual degeneracy "
3373  "rate for separator execution", &sepadata->dualdegeneracyratethreshold, FALSE,
3375 
3376  SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/varconsratiothreshold", "minimum variable-constraint "
3377  "ratio on optimal face for separator execution", &sepadata->varconsratiothreshold, FALSE,
3379 
3380  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/muparamconst", "is the mu parameter (factor for step "
3381  "length) constant?" , &sepadata->muparamconst, TRUE, DEFAULT_MUPARAMCONST, NULL, NULL) );
3382 
3383  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/separaterows", "separate rows with integral slack?",
3384  &sepadata->separaterows, TRUE, DEFAULT_SEPARATEROWS, NULL, NULL) );
3385 
3386  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/sortcutoffsol", "sort fractional integer columns"
3387  "based on fractionality?" , &sepadata->sortcutoffsol, TRUE, DEFAULT_SORTCUTOFFSOL, NULL, NULL) );
3388 
3389  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/sidetypebasis", "choose side types of row (lhs/rhs) "
3390  "based on basis information?", &sepadata->sidetypebasis, TRUE, DEFAULT_SIDETYPEBASIS, NULL, NULL) );
3391 
3392  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/dynamiccuts", "should generated cuts be removed from "
3393  "LP if they are no longer tight?", &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) );
3394 
3395  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/makeintegral", "try to scale all cuts to integral "
3396  "coefficients?", &sepadata->makeintegral, TRUE, DEFAULT_MAKEINTEGRAL, NULL, NULL) );
3397 
3398  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/forcecuts", "force cuts to be added to the LP?",
3399  &sepadata->forcecuts, TRUE, DEFAULT_FORCECUTS, NULL, NULL) );
3400 
3401  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/delayedcuts", "should cuts be added to the delayed cut "
3402  "pool", &sepadata->delayedcuts, TRUE, DEFAULT_DELAYEDCUTS, NULL, NULL) );
3403 
3404  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/allowlocal", "should locally valid cuts be generated?",
3405  &sepadata->allowlocal, TRUE, DEFAULT_ALLOWLOCAL, NULL, NULL) );
3406 
3407  SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/aggregatecuts", "aggregate all generated cuts using the "
3408  "Lagrangian multipliers?", &sepadata->aggregatecuts, TRUE, DEFAULT_AGGREGATECUTS, NULL, NULL) );
3409 
3410  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrounds", "maximal number of separation rounds per node "
3411  "(-1: unlimited)", &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
3412 
3413  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxroundsroot", "maximal number of separation rounds in "
3414  "the root node (-1: unlimited)", &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL,
3415  NULL) );
3416 
3417  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/perroundnmaxlpiters", "maximal number of separating LP "
3418  "iterations per separation round (-1: unlimited)", &sepadata->perroundnmaxlpiters, FALSE,
3419  DEFAULT_PERROUNDMAXLPITERS, -1, INT_MAX, NULL, NULL) );
3420 
3421  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxcutsperlp", "maximal number of cuts separated per "
3422  "Lagromory LP in the non-root node", &sepadata->nmaxcutsperlp, FALSE, DEFAULT_PERLPMAXCUTS, 0, INT_MAX, NULL,
3423  NULL));
3424 
3425  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxcutsperlproot", "maximal number of cuts separated per "
3426  "Lagromory LP in the root node", &sepadata->nmaxcutsperlproot, FALSE, DEFAULT_PERLPMAXCUTSROOT, 0, INT_MAX,
3427  NULL, NULL));
3428 
3429  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxmainiters", "maximal number of main loop iterations of "
3430  "the relax-and-cut algorithm", &sepadata->nmaxmainiters, TRUE, DEFAULT_MAXMAINITERS, 0, INT_MAX, NULL, NULL)
3431  );
3432 
3433  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxsubgradientiters", "maximal number of subgradient loop "
3434  "iterations of the relax-and-cut algorithm", &sepadata->nmaxsubgradientiters, TRUE,
3435  DEFAULT_MAXSUBGRADIENTITERS, 0, INT_MAX, NULL, NULL) );
3436 
3437  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/cutgenfreq", "frequency of subgradient iterations for "
3438  "generating cuts", &sepadata->cutgenfreq, TRUE, DEFAULT_CUTGENFREQ, 0, INT_MAX, NULL, NULL) );
3439 
3440  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/cutaddfreq", "frequency of subgradient iterations for "
3441  "adding cuts to objective function", &sepadata->cutaddfreq, TRUE, DEFAULT_CUTADDFREQ, 0, INT_MAX, NULL, NULL) );
3442 
3443  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxlagrangianvalsforavg", "maximal number of iterations "
3444  "for rolling average of Lagrangian value", &sepadata->nmaxlagrangianvalsforavg, TRUE,
3445  DEFAULT_MAXLAGRANGIANVALSFORAVG, 0, INT_MAX, NULL, NULL) );
3446 
3447  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxconsecitersformuupdate", "consecutive number of "
3448  "iterations used to determine if mu needs to be backtracked", &sepadata->nmaxconsecitersformuupdate, TRUE,
3449  DEFAULT_MAXCONSECITERSFORMUUPDATE, 0, INT_MAX, NULL, NULL) );
3450 
3451  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/projectiontype", "the ball into which the Lagrangian multipliers "
3452  "are projected for stabilization (0: no projection, 1: L1-norm ball projection, 2: L2-norm ball projection, 3: "
3453  "L_inf-norm ball projection)", &sepadata->projectiontype, TRUE, DEFAULT_PROJECTIONTYPE, 0, 3, NULL, NULL));
3454 
3455  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/stabilitycentertype", "type of stability center for "
3456  "taking weighted average of Lagrangian multipliers for stabilization (0: no weighted stabilization, 1: best "
3457  "Lagrangian multipliers)", &sepadata->stabilitycentertype, TRUE, DEFAULT_STABILITYCENTERTYPE, 0, 1, NULL,
3458  NULL) );
3459 
3460  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/optimalfacepriority", "priority of the optimal face for "
3461  "separator execution (0: low priority, 1: medium priority, 2: high priority)",
3462  &sepadata->optimalfacepriority, TRUE, DEFAULT_OPTIMALFACEPRIORITY, 0, 2, NULL, NULL) );
3463 
3464  SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/minrestart", "minimum restart round for separator "
3465  "execution (0: from beginning of the instance solving, >= n with n >= 1: from restart round n)",
3466  &sepadata->minrestart, TRUE, DEFAULT_MINRESTART, 0, INT_MAX, NULL, NULL) );
3467 
3468  return SCIP_OKAY;
3469 }
enum SCIP_Result SCIP_RESULT
Definition: type_result.h:61
void SCIPfreeRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen)
#define DEFAULT_FORCECUTS
static SCIP_RETCODE checkLagrangianDualTermination(SCIP_SEPADATA *sepadata, int nnewaddedsoftcuts, int nyettoaddsoftcuts, SCIP_Bool objvecsdiffer, int ngeneratedcurrroundcuts, int nmaxgeneratedperroundcuts, int ncurrroundlpiters, int depth, SCIP_Bool *terminate)
SCIP_Bool SCIPisFeasZero(SCIP *scip, SCIP_Real val)
void SCIPaggrRowFree(SCIP *scip, SCIP_AGGRROW **aggrrow)
Definition: cuts.c:1763
SCIP_RETCODE SCIPgetLPBInvRow(SCIP *scip, int r, SCIP_Real *coefs, int *inds, int *ninds)
Definition: scip_lp.c:714
SCIP_Real SCIPgetSolvingTime(SCIP *scip)
Definition: scip_timing.c:378
#define NULL
Definition: def.h:267
SCIP_RETCODE SCIPlpiFree(SCIP_LPI **lpi)
Definition: lpi_clp.cpp:643
primal heuristic that tries a given solution
SCIP_RETCODE SCIPcacheRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1635
SCIP_RETCODE SCIPlpiSetState(SCIP_LPI *lpi, BMS_BLKMEM *blkmem, const SCIP_LPISTATE *lpistate)
Definition: lpi_clp.cpp:3429
#define DEFAULT_PERROUNDCUTSFACTOR
static SCIP_RETCODE l1BallProjection(SCIP *scip, SCIP_Real *dualvector, int dualvectorlen, SCIP_Real radius)
SCIP_Longint SCIPgetNLPIterations(SCIP *scip)
static void updateSubgradient(SCIP *scip, SCIP_SOL *sol, SCIP_ROW **cuts, int ncuts, SCIP_Real *subgradient, SCIP_Real *dualvector, SCIP_Bool *subgradientzero, int *ncutviols, SCIP_Real *maxcutviol, int *nnzsubgradientdualprod, SCIP_Real *maxnzsubgradientdualprod)
#define DEFAULT_PERROUNDCUTSFACTORROOT
#define SEPA_PRIORITY
static SCIP_RETCODE addGMICutsAsSoftConss(SCIP_Real *dualvector, int ngeneratedcuts, int *naddedcuts, int *nnewaddedsoftcuts)
static SCIP_DECL_SEPAINIT(sepaInitLagromory)
SCIP_RETCODE SCIPflushRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1658
SCIP_RETCODE SCIPlpiGetSol(SCIP_LPI *lpi, SCIP_Real *objval, SCIP_Real *primsol, SCIP_Real *dualsol, SCIP_Real *activity, SCIP_Real *redcost)
Definition: lpi_clp.cpp:2788
#define DEFAULT_DELAYEDCUTS
SCIP_RETCODE SCIPgetRealParam(SCIP *scip, const char *name, SCIP_Real *value)
Definition: scip_param.c:307
#define SCIP_MAXSTRLEN
Definition: def.h:288
SCIP_RETCODE SCIPlpiSolvePrimal(SCIP_LPI *lpi)
Definition: lpi_clp.cpp:1805
static SCIP_RETCODE generateGMICuts(SCIP *scip, SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata, int mainiternum, int subgradientiternum, SCIP_SOL *sol, SCIP_Real *solvals, int nmaxgeneratedperroundcuts, SCIP_Bool allowlocal, SCIP_ROW **generatedcurrroundcuts, SCIP_Real *generatedcutefficacies, int ngeneratedcurrroundcuts, int *ngeneratednewcuts, int depth, SCIP_Bool *cutoff)
SCIP_RETCODE SCIPaddVarToRow(SCIP *scip, SCIP_ROW *row, SCIP_VAR *var, SCIP_Real val)
Definition: scip_lp.c:1701
int SCIProwGetNNonz(SCIP_ROW *row)
Definition: lp.c:17213
#define SQR(x)
Definition: def.h:214
SCIP_Bool SCIPisPositive(SCIP *scip, SCIP_Real val)
SCIP_RETCODE SCIPgetLPDualDegeneracy(SCIP *scip, SCIP_Real *degeneracy, SCIP_Real *varconsratio)
Definition: scip_lp.c:2792
#define QUAD_ARRAY_STORE(a, idx, x)
Definition: dbldblarith.h:55
static SCIP_RETCODE createLPWithSoftCuts(SCIP *scip, SCIP_SEPADATA *sepadata)
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static SCIP_RETCODE solveLagrangianDual(SCIP *scip, SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata, SCIP_SOL *sol, SCIP_Real *solvals, int mainiternum, SCIP_Real ubparam, int depth, SCIP_Bool allowlocal, int nmaxgeneratedperroundcuts, SCIP_Real *origobjcoefs, SCIP_Real origobjoffset, SCIP_Real *dualvector, int *nsoftcuts, SCIP_ROW **generatedcurrroundcuts, SCIP_Real *generatedcutefficacies, int *ngeneratedcutsperiter, int *ngeneratedcurrroundcuts, int *ncurrroundlpiters, SCIP_Bool *cutoff, SCIP_Real *bestlagrangianval, SCIP_Real *bestdualvector, int *bestdualvectorlen, int *nbestdualupdates, int *totaliternum)
const char * SCIProwGetName(SCIP_ROW *row)
Definition: lp.c:17351
#define DEFAULT_AGGREGATECUTS
#define DEFAULT_SORTCUTOFFSOL
#define DEFAULT_RADIUSMAX
SCIP_Bool SCIPisFeasNegative(SCIP *scip, SCIP_Real val)
SCIP_RETCODE SCIPgetLPI(SCIP *scip, SCIP_LPI **lpi)
Definition: scip_lp.c:985
SCIP_RETCODE SCIPaddDelayedPoolCut(SCIP *scip, SCIP_ROW *row)
Definition: scip_cut.c:641
int SCIPgetMaxDepth(SCIP *scip)
int SCIProwGetNLPNonz(SCIP_ROW *row)
Definition: lp.c:17227
static SCIP_RETCODE sepadataCreate(SCIP *scip, SCIP_SEPADATA **sepadata)
static SCIP_RETCODE generateInitCutPool(SCIP *scip, SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata, int mainiternum, SCIP_SOL *sol, SCIP_Real *solvals, int nmaxgeneratedperroundcuts, SCIP_Bool allowlocal, SCIP_ROW **generatedcurrroundcuts, SCIP_Real *generatedcutefficacies, int *ngeneratedcutsperiter, int *ngeneratedcurrroundcuts, int depth, SCIP_Bool *cutoff)
SCIP_Real SCIProwGetLhs(SCIP_ROW *row)
Definition: lp.c:17292
#define FALSE
Definition: def.h:94
SCIP_Bool SCIPcolIsIntegral(SCIP_COL *col)
Definition: lp.c:17072
SCIP_Real SCIPcolGetUb(SCIP_COL *col)
Definition: lp.c:16973
SCIP_Real SCIPcolGetObj(SCIP_COL *col)
Definition: lp.c:16953
#define SEPA_DESC
SCIP_Real SCIPinfinity(SCIP *scip)
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10877
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
#define TRUE
Definition: def.h:93
const char * SCIPsepaGetName(SCIP_SEPA *sepa)
Definition: sepa.c:743
static SCIP_DECL_SEPAFREE(sepaFreeLagromory)
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63
#define DEFAULT_CUTSFILTERFACTOR
SCIP_RETCODE SCIPlpiSetRealpar(SCIP_LPI *lpi, SCIP_LPPARAM type, SCIP_Real dval)
Definition: lpi_clp.cpp:3833
SCIP_Longint SCIPgetNRootFirstLPIterations(SCIP *scip)
static SCIP_RETCODE constructCutRow(SCIP *scip, SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata, int mainiternum, int subgradientiternum, int cutnnz, int *cutinds, SCIP_Real *cutcoefs, SCIP_Real cutefficacy, SCIP_Real cutrhs, SCIP_Bool cutislocal, int cutrank, SCIP_ROW **generatedcuts, SCIP_Real *generatedcutefficacies, int ngeneratedcurrroundcuts, int *ngeneratednewcuts, SCIP_Bool *cutoff)
static void updateBallRadius(SCIP *scip, SCIP_SEPADATA *sepadata, SCIP_Real maxviolscore, SCIP_Real maxviolscoreold, SCIP_Real nviolscore, SCIP_Real nviolscoreold, int nlpiters, SCIP_Real *ballradius)
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:108
static SCIP_RETCODE solveLagromoryLP(SCIP *scip, SCIP_SEPADATA *sepadata, int depth, SCIP_Real origobjoffset, SCIP_Bool *solfound, SCIP_SOL *sol, SCIP_Real *solvals, SCIP_Real *objval, int *ncurrroundlpiters)
#define DEFAULT_MAXROUNDS
SCIP_MESSAGEHDLR * SCIPgetMessagehdlr(SCIP *scip)
Definition: scip_message.c:88
static SCIP_RETCODE updateObjectiveVector(SCIP *scip, SCIP_Real *dualvector, SCIP_ROW **cuts, int ncuts, SCIP_Real *origobjcoefs, SCIP_Bool *objvecsdiffer)
void SCIPsortDownRealInt(SCIP_Real *realarray, int *intarray, int len)
SCIP_Real SCIPfeasFrac(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:136
#define QUAD_ASSIGN(a, constant)
Definition: dbldblarith.h:51
enum SCIP_LPSolStat SCIP_LPSOLSTAT
Definition: type_lp.h:51
#define DEFAULT_MAXSUBGRADIENTITERS
#define SCIPallocBlockMemory(scip, ptr)
Definition: scip_mem.h:89
SCIP_RETCODE SCIPgetLPColsData(SCIP *scip, SCIP_COL ***cols, int *ncols)
Definition: scip_lp.c:471
SCIP_RETCODE SCIPsetSepaCopy(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPACOPY((*sepacopy)))
Definition: scip_sepa.c:151
#define SCIPdebugMsg
Definition: scip_message.h:78
SCIP_RETCODE SCIPaddIntParam(SCIP *scip, const char *name, const char *desc, int *valueptr, SCIP_Bool isadvanced, int defaultvalue, int minvalue, int maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:83
#define DEFAULT_SEPARATEROWS
SCIP_RETCODE SCIPlpiCreate(SCIP_LPI **lpi, SCIP_MESSAGEHDLR *messagehdlr, const char *name, SCIP_OBJSEN objsen)
Definition: lpi_clp.cpp:531
#define DEFAULT_DELTASLAB1UB
SCIP_Real SCIPepsilon(SCIP *scip)
#define QUAD_TO_DBL(x)
Definition: dbldblarith.h:49
SCIP_RETCODE SCIPlpiAddCols(SCIP_LPI *lpi, int ncols, const SCIP_Real *obj, const SCIP_Real *lb, const SCIP_Real *ub, char **colnames, int nnonz, const int *beg, const int *ind, const SCIP_Real *val)
Definition: lpi_clp.cpp:758
#define DEFAULT_DELTASLAB2UB
#define MAKECONTINTEGRAL
SCIP_Real SCIPgetRowMaxActivity(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1956
#define DEFAULT_CUTGENFREQ
#define DEFAULT_AWAY
SCIP_SEPADATA * SCIPsepaGetData(SCIP_SEPA *sepa)
Definition: sepa.c:633
static SCIP_RETCODE stabilizeDualVector(SCIP *scip, SCIP_SEPADATA *sepadata, SCIP_Real *dualvector, int dualvectorlen, SCIP_Real *bestdualvector, int bestdualvectorlen, int nbestdualupdates, int subgradientiternum, int totaliternum, SCIP_Real maxviolscore, SCIP_Real maxviolscoreold, SCIP_Real nviolscore, SCIP_Real nviolscoreold, int nlpiters, SCIP_Real *ballradius)
SCIP_RETCODE SCIPheurPassSolTrySol(SCIP *scip, SCIP_HEUR *heur, SCIP_SOL *sol)
Definition: heur_trysol.c:252
#define SCIPallocCleanBufferArray(scip, ptr, num)
Definition: scip_mem.h:142
#define DEFAULT_ROOTLPITERLIMITFACTOR
SCIP_Bool SCIPisLPSolBasic(SCIP *scip)
Definition: scip_lp.c:667
#define SEPA_USESSUBSCIP
SCIP_RETCODE SCIPlpiAddRows(SCIP_LPI *lpi, int nrows, const SCIP_Real *lhs, const SCIP_Real *rhs, char **rownames, int nnonz, const int *beg, const int *ind, const SCIP_Real *val)
Definition: lpi_clp.cpp:914
SCIP_Real SCIPcolGetPrimsol(SCIP_COL *col)
Definition: lp.c:16996
#define DEFAULT_MAXCONSECITERSFORMUUPDATE
#define SEPA_FREQ
#define DEFAULT_MUSLAB2FACTOR
SCIP_HEUR * SCIPfindHeur(SCIP *scip, const char *name)
Definition: scip_heur.c:258
static void updateLagrangianValue(SCIP *scip, SCIP_Real objval, SCIP_Real *dualvector, SCIP_ROW **cuts, int ncuts, SCIP_Real *lagrangianval)
SCIP_RETCODE SCIPsolveDiveLP(SCIP *scip, int itlim, SCIP_Bool *lperror, SCIP_Bool *cutoff)
Definition: scip_lp.c:2678
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPsetSepaExit(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAEXIT((*sepaexit)))
Definition: scip_sepa.c:199
SCIP_RETCODE SCIPincludeSepaLagromory(SCIP *scip)
#define DEFAULT_MAXLAGRANGIANVALSFORAVG
SCIP_Bool SCIProwIsLocal(SCIP_ROW *row)
Definition: lp.c:17401
#define DEFAULT_MUPARAMUB
#define FIXINTEGRALRHS
int SCIPsepaGetNCallsAtNode(SCIP_SEPA *sepa)
Definition: sepa.c:880
SCIP_RETCODE SCIPlpiFreeState(SCIP_LPI *lpi, BMS_BLKMEM *blkmem, SCIP_LPISTATE **lpistate)
Definition: lpi_clp.cpp:3503
SCIP_Bool SCIPisEfficacious(SCIP *scip, SCIP_Real efficacy)
Definition: scip_cut.c:135
BMS_BLKMEM * SCIPblkmem(SCIP *scip)
Definition: scip_mem.c:57
static SCIP_RETCODE aggregateGeneratedCuts(SCIP *scip, SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata, SCIP_ROW **generatedcuts, SCIP_Real *bestdualvector, int bestdualvectorlen, SCIP_ROW **aggrcuts, int *naggrcuts, SCIP_Bool *cutoff)
#define QUAD(x)
Definition: dbldblarith.h:47
static SCIP_RETCODE checkMainLoopTermination(SCIP_SEPADATA *sepadata, SCIP_Bool cutoff, SCIP_Bool dualvecsdiffer, int ngeneratedcurrroundcuts, int nsoftcuts, int nmaxgeneratedperroundcuts, int ncurrroundlpiters, int depth, SCIP_Bool *terminate)
SCIP_Real SCIPcolGetLb(SCIP_COL *col)
Definition: lp.c:16963
SCIP_Bool SCIProwIsIntegral(SCIP_ROW *row)
Definition: lp.c:17391
void SCIPsepaSetData(SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata)
Definition: sepa.c:643
static SCIP_RETCODE updateDualVector(SCIP *scip, SCIP_SEPADATA *sepadata, SCIP_Real *dualvector1, SCIP_Real *dualvector2, int dualvector2len, int ndualvector2updates, int subgradientiternum, int totaliternum, SCIP_Real steplength, SCIP_Real *subgradient, int ncuts, SCIP_Bool backtrack, SCIP_Real maxviolscore, SCIP_Real maxviolscoreold, SCIP_Real nviolscore, SCIP_Real nviolscoreold, int nlpiters, SCIP_Bool *dualvecsdiffer, SCIP_Real *ballradius)
#define DEFAULT_MUSLAB1FACTOR
#define REALABS(x)
Definition: def.h:197
static SCIP_DECL_SEPACOPY(sepaCopyLagromory)
int SCIPgetNLPRows(SCIP *scip)
Definition: scip_lp.c:626
#define QUAD_ARRAY_LOAD(r, a, idx)
Definition: dbldblarith.h:54
#define SCIP_CALL(x)
Definition: def.h:380
static SCIP_RETCODE solveLPWithHardCuts(SCIP *scip, SCIP_SEPADATA *sepadata, SCIP_Bool *solfound, SCIP_SOL *sol, SCIP_Real *solvals)
SCIP_Real SCIProwGetRhs(SCIP_ROW *row)
Definition: lp.c:17302
SCIP_RETCODE SCIPaddRow(SCIP *scip, SCIP_ROW *row, SCIP_Bool forcecut, SCIP_Bool *infeasible)
Definition: scip_cut.c:250
#define DEFAULT_OPTIMALFACEPRIORITY
#define SCIPquadprecProdDD(r, a, b)
Definition: dbldblarith.h:58
SCIP_Bool SCIProwIsModifiable(SCIP_ROW *row)
Definition: lp.c:17411
SCIP_COL ** SCIProwGetCols(SCIP_ROW *row)
Definition: lp.c:17238
#define DEFAULT_TOTALLPITERLIMITFACTOR
SCIP_RETCODE SCIPincludeSepaBasic(SCIP *scip, SCIP_SEPA **sepa, const char *name, const char *desc, int priority, int freq, SCIP_Real maxbounddist, SCIP_Bool usessubscip, SCIP_Bool delay, SCIP_DECL_SEPAEXECLP((*sepaexeclp)), SCIP_DECL_SEPAEXECSOL((*sepaexecsol)), SCIP_SEPADATA *sepadata)
Definition: scip_sepa.c:109
SCIP_RETCODE SCIPcreateRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen, unsigned int initialseed, SCIP_Bool useglobalseed)
SCIP_Longint SCIPsepaGetNCalls(SCIP_SEPA *sepa)
Definition: sepa.c:860
#define DEFAULT_MINRESTART
SCIP_RETCODE SCIPcalcMIR(SCIP *scip, SCIP_SOL *sol, SCIP_Bool postprocess, SCIP_Real boundswitch, SCIP_Bool usevbds, SCIP_Bool allowlocal, SCIP_Bool fixintegralrhs, int *boundsfortrans, SCIP_BOUNDTYPE *boundtypesfortrans, SCIP_Real minfrac, SCIP_Real maxfrac, SCIP_Real scale, SCIP_AGGRROW *aggrrow, SCIP_Real *cutcoefs, SCIP_Real *cutrhs, int *cutinds, int *cutnnz, SCIP_Real *cutefficacy, int *cutrank, SCIP_Bool *cutislocal, SCIP_Bool *success)
Definition: cuts.c:3879
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:124
SCIP_RETCODE SCIPsetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var, SCIP_Real val)
Definition: scip_sol.c:1077
SCIP_Real * SCIProwGetVals(SCIP_ROW *row)
Definition: lp.c:17248
#define DEFAULT_MUBACKTRACKFACTOR
#define DEFAULT_TOTALCUTSFACTOR
#define DEFAULT_DUALDEGENERACYRATETHRESHOLD
#define DEFAULT_PERROUNDMAXLPITERS
#define SCIP_Bool
Definition: def.h:91
static SCIP_DECL_SEPAEXIT(sepaExitLagromory)
SCIP_LPSOLSTAT SCIPgetLPSolstat(SCIP *scip)
Definition: scip_lp.c:168
#define QUAD_ARRAY_SIZE(size)
Definition: dbldblarith.h:53
static SCIP_RETCODE separateCuts(SCIP *scip, SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata, SCIP_Real ubparam, int depth, SCIP_Bool allowlocal, SCIP_RESULT *result)
#define DEFAULT_MAXROUNDSROOT
#define DEFAULT_MUPARAMLB
#define DEFAULT_PROJECTIONTYPE
#define DEFAULT_PERLPMAXCUTS
lagromory separator
SCIP_Real SCIPlpiInfinity(SCIP_LPI *lpi)
Definition: lpi_clp.cpp:3919
#define DEFAULT_RADIUSUPDATEWEIGHT
SCIP_Bool SCIPlpiIsDualFeasible(SCIP_LPI *lpi)
Definition: lpi_clp.cpp:2609
SCIP_RETCODE SCIPlpiGetState(SCIP_LPI *lpi, BMS_BLKMEM *blkmem, SCIP_LPISTATE **lpistate)
Definition: lpi_clp.cpp:3389
int SCIPgetRowNumIntCols(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1886
SCIP_RETCODE SCIPaddPoolCut(SCIP *scip, SCIP_ROW *row)
Definition: scip_cut.c:361
#define MIN(x, y)
Definition: def.h:243
SCIP_RETCODE SCIPcreateEmptyRowSepa(SCIP *scip, SCIP_ROW **row, SCIP_SEPA *sepa, const char *name, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool removable)
Definition: scip_lp.c:1453
static void l2BallProjection(SCIP *scip, SCIP_Real *dualvector, int dualvectorlen, SCIP_Real radius)
static SCIP_RETCODE createLPWithHardCuts(SCIP *scip, SCIP_SEPADATA *sepadata, SCIP_ROW **cuts, int ncuts)
SCIP_RETCODE SCIPfreeSol(SCIP *scip, SCIP_SOL **sol)
Definition: scip_sol.c:841
#define SEPA_DELAY
#define DEFAULT_STABILITYCENTERTYPE
int SCIPgetNRuns(SCIP *scip)
#define DEFAULT_PERLPITERFACTOR
#define POSTPROCESS
SCIP_RETCODE SCIPgetLPBasisInd(SCIP *scip, int *basisind)
Definition: scip_lp.c:686
SCIP_COL ** SCIPgetLPCols(SCIP *scip)
Definition: scip_lp.c:506
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
#define DEFAULT_UBPARAMPOSFACTOR
SCIP_Real SCIPgetRowSolActivity(SCIP *scip, SCIP_ROW *row, SCIP_SOL *sol)
Definition: scip_lp.c:2144
#define BMSclearMemory(ptr)
Definition: memory.h:129
SCIP_Bool SCIPisCutNew(SCIP *scip, SCIP_ROW *row)
Definition: scip_cut.c:343
static void linfBallProjection(SCIP *scip, SCIP_Real *dualvector, int dualvectorlen, SCIP_Real radius)
#define SCIPquadprecSumQQ(r, a, b)
Definition: dbldblarith.h:67
#define DEFAULT_MAKEINTEGRAL
SCIP_Real SCIPrandomGetReal(SCIP_RANDNUMGEN *randnumgen, SCIP_Real minrandval, SCIP_Real maxrandval)
Definition: misc.c:10130
int SCIProwGetRank(SCIP_ROW *row)
Definition: lp.c:17381
#define SCIP_REAL_MAX
Definition: def.h:174
#define SEPA_MAXBOUNDDIST
#define MAXAGGRLEN(ncols)
#define DEFAULT_CUTADDFREQ
#define SCIP_LONGINT_FORMAT
Definition: def.h:165
SCIP_Real SCIProwGetConstant(SCIP_ROW *row)
Definition: lp.c:17258
SCIP_Real SCIPgetTransObjoffset(SCIP *scip)
Definition: scip_prob.c:1367
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition: scip_lp.c:1562
SCIP_Longint SCIPgetNNodeInitLPIterations(SCIP *scip)
#define SEPA_NAME
#define MAX(x, y)
Definition: def.h:239
#define DEFAULT_ALLOWLOCAL
SCIP_RETCODE SCIPsetSepaFree(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAFREE((*sepafree)))
Definition: scip_sepa.c:167
SCIP_RETCODE SCIPsetSepaInit(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAINIT((*sepainit)))
Definition: scip_sepa.c:183
SCIP_Real SCIPgetLPObjval(SCIP *scip)
Definition: scip_lp.c:247
SCIP_SOL * SCIPgetBestSol(SCIP *scip)
Definition: scip_sol.c:2169
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPaggrRowSumRows(SCIP *scip, SCIP_AGGRROW *aggrrow, SCIP_Real *weights, int *rowinds, int nrowinds, SCIP_Bool sidetypebasis, SCIP_Bool allowlocal, int negslack, int maxaggrlen, SCIP_Bool *valid)
Definition: cuts.c:2287
SCIP_RETCODE SCIPchgVarObjDive(SCIP *scip, SCIP_VAR *var, SCIP_Real newobj)
Definition: scip_lp.c:2378
#define DEFAULT_UBPARAMNEGFACTOR
SCIP_VAR * SCIPcolGetVar(SCIP_COL *col)
Definition: lp.c:17042
#define DEFAULT_DYNAMICCUTS
void SCIProwChgRank(SCIP_ROW *row, int rank)
Definition: lp.c:17534
static SCIP_RETCODE addCuts(SCIP *scip, SCIP_SEPADATA *sepadata, SCIP_ROW **cuts, int ncuts, SCIP_Longint maxdnom, SCIP_Real maxscale, int *naddedcuts, SCIP_Bool *cutoff)
#define DEFAULT_PERROOTLPITERFACTOR
#define RANDSEED
#define DEFAULT_VARCONSRATIOTHRESHOLD
static SCIP_RETCODE sepadataFree(SCIP *scip, SCIP_SEPADATA **sepadata)
#define DEFAULT_MAXMAINITERS
SCIP_Bool SCIPisFeasPositive(SCIP *scip, SCIP_Real val)
#define DEFAULT_MUPARAMINIT
#define DEFAULT_MUPARAMCONST
static void updateStepLength(SCIP *scip, SCIP_Real muparam, SCIP_Real ubparam, SCIP_Real lagrangianval, SCIP_Real *subgradient, int ncuts, SCIP_Real *steplength)
#define SCIP_Real
Definition: def.h:173
#define DEFAULT_MUSLAB3FACTOR
SCIP_Bool SCIPlpiIsPrimalFeasible(SCIP_LPI *lpi)
Definition: lpi_clp.cpp:2521
#define SCIPfreeCleanBufferArray(scip, ptr)
Definition: scip_mem.h:146
SCIP_Bool SCIPisStopped(SCIP *scip)
Definition: scip_general.c:718
SCIP_Real SCIPgetRowMinActivity(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1939
#define DEFAULT_CONST
static SCIP_RETCODE updateMuSteplengthParam(SCIP *scip, SCIP_SEPADATA *sepadata, int subgradientiternum, SCIP_Real ubparam, SCIP_Real *lagrangianvals, SCIP_Real bestlagrangianval, SCIP_Real avglagrangianval, SCIP_Real *muparam, SCIP_Bool *backtrack)
#define DEFAULT_RADIUSMIN
SCIP_RETCODE SCIPaggrRowCreate(SCIP *scip, SCIP_AGGRROW **aggrrow)
Definition: cuts.c:1731
#define SCIP_Longint
Definition: def.h:158
#define DEFAULT_PERROUNDLPITERLIMITFACTOR
SCIP_VARTYPE SCIPvarGetType(SCIP_VAR *var)
Definition: var.c:17585
SCIP_OBJSENSE SCIPgetObjsense(SCIP *scip)
Definition: scip_prob.c:1225
SCIP_RETCODE SCIPstartDive(SCIP *scip)
Definition: scip_lp.c:2242
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define BOUNDSWITCH
SCIP_Real SCIPsumepsilon(SCIP *scip)
SCIP_RETCODE SCIPgetLPRowsData(SCIP *scip, SCIP_ROW ***rows, int *nrows)
Definition: scip_lp.c:570
static SCIP_RETCODE deleteLPWithSoftCuts(SCIP *scip, SCIP_SEPADATA *sepadata)
struct BMS_BlkMem BMS_BLKMEM
Definition: memory.h:437
SCIP_RETCODE SCIPendDive(SCIP *scip)
Definition: scip_lp.c:2291
#define USEVBDS
SCIP_Real SCIPceil(SCIP *scip, SCIP_Real val)
static SCIP_RETCODE weightedDualVector(SCIP *scip, SCIP_SEPADATA *sepadata, SCIP_Real *dualvector, int dualvectorlen, SCIP_Real *stabilitycenter, int stabilitycenterlen, int nbestdualupdates, int totaliternum)
#define DEFAULT_SIDETYPEBASIS
int SCIPcolGetLPPos(SCIP_COL *col)
Definition: lp.c:17093
SCIP_RETCODE SCIPmakeRowIntegral(SCIP *scip, SCIP_ROW *row, SCIP_Real mindelta, SCIP_Real maxdelta, SCIP_Longint maxdnom, SCIP_Real maxscale, SCIP_Bool usecontvars, SCIP_Bool *success)
Definition: scip_lp.c:1844
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1217
#define DEFAULT_PERLPMAXCUTSROOT
SCIP_Real SCIPgetVarObjDive(SCIP *scip, SCIP_VAR *var)
Definition: scip_lp.c:2587
#define DEFAULT_RADIUSINIT
static SCIP_DECL_SEPAEXECLP(sepaExeclpLagromory)
SCIP_RETCODE SCIPaddRealParam(SCIP *scip, const char *name, const char *desc, SCIP_Real *valueptr, SCIP_Bool isadvanced, SCIP_Real defaultvalue, SCIP_Real minvalue, SCIP_Real maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:139
struct SCIP_SepaData SCIP_SEPADATA
Definition: type_sepa.h:52
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:57
SCIP_Real SCIPgetRowActivity(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:2104
SCIP_RETCODE SCIPcreateSol(SCIP *scip, SCIP_SOL **sol, SCIP_HEUR *heur)
Definition: scip_sol.c:184