1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24
25 /**@file sepa_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 */
186
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 */
231
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
355 SCIP_RETCODE createLPWithSoftCuts(
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
450 SCIP_RETCODE deleteLPWithSoftCuts(
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
467 SCIP_RETCODE createLPWithHardCuts(
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 */
517 if( SCIPgetObjsense(scip) == SCIP_OBJSENSE_MAXIMIZE )
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
667 SCIP_RETCODE sepadataFree(
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
704 SCIP_RETCODE updateMuSteplengthParam(
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;
724
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
895 void updateLagrangianValue(
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);
961
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
1010 SCIP_RETCODE l1BallProjection(
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 */
(1) Event cond_true: |
Condition "i < dualvectorlen", taking true branch. |
(3) Event loop_begin: |
Jumped back to beginning of loop. |
(4) Event cond_true: |
Condition "i < dualvectorlen", taking true branch. |
(6) Event loop_begin: |
Jumped back to beginning of loop. |
(7) Event cond_true: |
Condition "i < dualvectorlen", taking true branch. |
(9) Event loop_begin: |
Jumped back to beginning of loop. |
(10) Event cond_false: |
Condition "i < dualvectorlen", taking false branch. |
1032 for( int i = 1; i < dualvectorlen; i++ )
1033 {
1034 val += REALABS(dualvector[i]);
(2) Event loop: |
Jumping back to the beginning of the loop. |
(5) Event loop: |
Jumping back to the beginning of the loop. |
(8) Event loop: |
Jumping back to the beginning of the loop. |
(11) Event loop_end: |
Reached end of loop. |
1035 }
1036
1037 /* if the vector of Lagrangian multipliers lies outside the L1-norm ball, then do the projection */
(12) Event cond_true: |
Condition "val - radius > scip->set->num_epsilon", taking true branch. |
1038 if( SCIPisGT(scip, val, radius) )
1039 {
(13) Event cond_false: |
Condition "(temp1vals = BMSallocBufferMemoryArray_call(SCIPcleanbuffer(scip), (size_t)(ptrdiff_t)dualvectorlen, 8UL /* sizeof (*temp1vals) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/sepa_lagromory.c", 1040)) == NULL", taking false branch. |
(14) Event cond_false: |
Condition "(_restat_ = (((temp1vals = BMSallocBufferMemoryArray_call(SCIPcleanbuffer(scip), (size_t)(ptrdiff_t)dualvectorlen, 8UL /* sizeof (*temp1vals) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/sepa_lagromory.c", 1040)) == NULL) ? SCIP_NOMEMORY : SCIP_OKAY)) != SCIP_OKAY", taking false branch. |
(15) Event if_end: |
End of if statement. |
1040 SCIP_CALL( SCIPallocCleanBufferArray(scip, &temp1vals, dualvectorlen) );
(16) Event cond_false: |
Condition "(temp2vals = BMSallocBufferMemoryArray_call(SCIPcleanbuffer(scip), (size_t)(ptrdiff_t)dualvectorlen, 8UL /* sizeof (*temp2vals) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/sepa_lagromory.c", 1041)) == NULL", taking false branch. |
(17) Event cond_false: |
Condition "(_restat_ = (((temp2vals = BMSallocBufferMemoryArray_call(SCIPcleanbuffer(scip), (size_t)(ptrdiff_t)dualvectorlen, 8UL /* sizeof (*temp2vals) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/sepa_lagromory.c", 1041)) == NULL) ? SCIP_NOMEMORY : SCIP_OKAY)) != SCIP_OKAY", taking false branch. |
(18) Event if_end: |
End of if statement. |
1041 SCIP_CALL( SCIPallocCleanBufferArray(scip, &temp2vals, dualvectorlen) );
(19) Event cond_false: |
Condition "(temp1inds = BMSallocBufferMemoryArray_call(SCIPbuffer(scip), (size_t)(ptrdiff_t)dualvectorlen, 4UL /* sizeof (*temp1inds) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/sepa_lagromory.c", 1042)) == NULL", taking false branch. |
(20) Event cond_false: |
Condition "(_restat_ = (((temp1inds = BMSallocBufferMemoryArray_call(SCIPbuffer(scip), (size_t)(ptrdiff_t)dualvectorlen, 4UL /* sizeof (*temp1inds) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/sepa_lagromory.c", 1042)) == NULL) ? SCIP_NOMEMORY : SCIP_OKAY)) != SCIP_OKAY", taking false branch. |
(21) Event if_end: |
End of if statement. |
1042 SCIP_CALL( SCIPallocBufferArray(scip, &temp1inds, dualvectorlen) );
(22) Event cond_false: |
Condition "(temp2inds = BMSallocBufferMemoryArray_call(SCIPbuffer(scip), (size_t)(ptrdiff_t)dualvectorlen, 4UL /* sizeof (*temp2inds) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/sepa_lagromory.c", 1043)) == NULL", taking false branch. |
(23) Event cond_false: |
Condition "(_restat_ = (((temp2inds = BMSallocBufferMemoryArray_call(SCIPbuffer(scip), (size_t)(ptrdiff_t)dualvectorlen, 4UL /* sizeof (*temp2inds) */, "/sapmnt/home1/d029903/my_work/SCIPSoPlex_coverity/scip-soplex-9.0.0_rc1/scip/src/scip/sepa_lagromory.c", 1043)) == NULL) ? SCIP_NOMEMORY : SCIP_OKAY)) != SCIP_OKAY", taking false branch. |
(24) Event if_end: |
End of if statement. |
1043 SCIP_CALL( SCIPallocBufferArray(scip, &temp2inds, dualvectorlen) );
(25) Event cond_true: |
Condition "i < dualvectorlen", taking true branch. |
(27) Event loop_begin: |
Jumped back to beginning of loop. |
(28) Event cond_true: |
Condition "i < dualvectorlen", taking true branch. |
(30) Event loop_begin: |
Jumped back to beginning of loop. |
(31) Event cond_false: |
Condition "i < dualvectorlen", taking false branch. |
1044 for( int i = 0; i < dualvectorlen; i++ )
1045 {
1046 temp1inds[i] = -1;
1047 temp2inds[i] = -1;
(26) Event loop: |
Jumping back to the beginning of the loop. |
(29) Event loop: |
Jumping back to the beginning of the loop. |
(32) Event loop_end: |
Reached end of loop. |
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
(33) Event cond_true: |
Condition "i < dualvectorlen", taking true branch. |
(39) Event loop_begin: |
Jumped back to beginning of loop. |
(40) Event cond_true: |
Condition "i < dualvectorlen", taking true branch. |
(53) Event loop_begin: |
Jumped back to beginning of loop. |
(54) Event cond_true: |
Condition "i < dualvectorlen", taking true branch. |
(65) Event loop_begin: |
Jumped back to beginning of loop. |
(66) Event cond_false: |
Condition "i < dualvectorlen", taking false branch. |
1056 for( int i = 1; i < dualvectorlen; i++ )
1057 {
(34) Event cond_true: |
Condition "fabs(dualvector[i]) - pivotparam > scip->set->num_epsilon", taking true branch. |
(41) Event cond_true: |
Condition "fabs(dualvector[i]) - pivotparam > scip->set->num_epsilon", taking true branch. |
(55) Event cond_true: |
Condition "fabs(dualvector[i]) - pivotparam > scip->set->num_epsilon", taking true branch. |
1058 if( SCIPisGT(scip, REALABS(dualvector[i]), pivotparam) )
1059 {
1060 pivotparam += ((REALABS(dualvector[i]) - pivotparam) / (temp1len + 1));
(35) Event cond_true: |
Condition "pivotparam - (fabs(dualvector[i]) - radius) > scip->set->num_epsilon", taking true branch. |
(42) Event cond_false: |
Condition "pivotparam - (fabs(dualvector[i]) - radius) > scip->set->num_epsilon", taking false branch. |
(56) Event cond_false: |
Condition "pivotparam - (fabs(dualvector[i]) - radius) > scip->set->num_epsilon", taking false branch. |
1061 if( SCIPisGT(scip, pivotparam, REALABS(dualvector[i]) - radius) )
1062 {
1063 temp1vals[temp1len] = REALABS(dualvector[i]);
1064 temp1inds[temp1len] = i;
1065 temp1len++;
(36) Event if_fallthrough: |
Falling through to end of if statement. |
1066 }
1067 else
(43) Event else_branch: |
Reached else branch. |
(57) Event else_branch: |
Reached else branch. |
1068 {
(44) Event cond_true: |
Condition "j < temp1len", taking true branch. |
(46) Event loop_begin: |
Jumped back to beginning of loop. |
(47) Event cond_true: |
Condition "j < temp1len", taking true branch. |
(49) Event loop_begin: |
Jumped back to beginning of loop. |
(50) Event cond_false: |
Condition "j < temp1len", taking false branch. |
(58) Event cond_true: |
Condition "j < temp1len", taking true branch. |
(60) Event loop_begin: |
Jumped back to beginning of loop. |
(61) Event cond_false: |
Condition "j < temp1len", taking false branch. |
1069 for( int j = 0; j < temp1len; j++ )
1070 {
1071 temp2vals[temp2len + j] = temp1vals[j];
1072 temp2inds[temp2len + j] = temp1inds[j];
(45) Event loop: |
Jumping back to the beginning of the loop. |
(48) Event loop: |
Jumping back to the beginning of the loop. |
(51) Event loop_end: |
Reached end of loop. |
(59) Event loop: |
Jumping back to the beginning of the loop. |
(62) Event loop_end: |
Reached end of loop. |
1073 }
1074 temp2len += temp1len;
1075 temp1vals[0] = REALABS(dualvector[i]);
1076 temp1inds[0] = i;
1077 temp1len = 1;
1078 pivotparam = REALABS(dualvector[i]) - radius;
(37) Event if_end: |
End of if statement. |
1079 }
1080 }
(38) Event loop: |
Jumping back to the beginning of the loop. |
(52) Event loop: |
Jumping back to the beginning of the loop. |
(64) Event loop: |
Jumping back to the beginning of the loop. |
(67) Event loop_end: |
Reached end of loop. |
1081 }
1082
(68) Event cond_true: |
Condition "i < temp2len", taking true branch. |
(72) Event loop_begin: |
Jumped back to beginning of loop. |
(73) Event cond_false: |
Condition "i < temp2len", taking false branch. |
1083 for( int i = 0; i < temp2len; i++ )
1084 {
(69) Event cond_false: |
Condition "temp2vals[i] - pivotparam > scip->set->num_epsilon", taking false branch. |
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);
(70) Event if_end: |
End of if statement. |
1091 }
(71) Event loop: |
Jumping back to the beginning of the loop. |
(74) Event loop_end: |
Reached end of loop. |
1092 }
1093
1094 temp1changed = TRUE;
1095 ntemp1removed = 0;
(76) Event cond_true: |
Condition "temp1changed", taking true branch. |
1096 while( temp1changed )
1097 {
1098 temp1changed = FALSE;
1099
(77) Event cond_true: |
Condition "i < temp1len", taking true branch. |
1100 for( int i = 0; i < temp1len; i++ )
1101 {
(78) Event cond_true: |
Condition "temp1inds[i] >= 0", taking true branch. |
(79) Event cond_true: |
Condition "temp1vals[i] - pivotparam <= scip->set->num_epsilon", taking true branch. |
1102 if( (temp1inds[i] >= 0) && SCIPisLE(scip, temp1vals[i], pivotparam) )
1103 {
1104 temp1inds[i] = -1;
1105 temp1changed = TRUE;
1106 ntemp1removed++;
1107 assert(temp1len - ntemp1removed > 0);
(81) Event divide_by_zero: |
In expression "(pivotparam - temp1vals[i]) / (temp1len - ntemp1removed)", division by expression "temp1len - ntemp1removed" which may be zero has undefined behavior. |
Also see events: |
[assignment][assignment][incr] |
1108 pivotparam += ((pivotparam - temp1vals[i]) / (temp1len - ntemp1removed));
1109 }
1110 }
1111 }
1112
1113 for( int i = 0; i < dualvectorlen; i++ )
1114 {
1115 term = REALABS(dualvector[i]);
1116 val = MAX(term - pivotparam, 0.0);
1117
1118 if( SCIPisPositive(scip, dualvector[i]) )
1119 {
1120 dualvector[i] = val;
1121 }
1122 else if( SCIPisNegative(scip, dualvector[i]) )
1123 {
1124 dualvector[i] = -val;
1125 }
1126 }
1127
1128 /* free memory */
1129 for( int i = 0; i < dualvectorlen; i++ )
1130 {
1131 temp2vals[i] = 0.0;
1132 temp1vals[i] = 0.0;
1133 }
1134 SCIPfreeBufferArray(scip, &temp2inds);
1135 SCIPfreeBufferArray(scip, &temp1inds);
1136 SCIPfreeCleanBufferArray(scip, &temp2vals);
1137 SCIPfreeCleanBufferArray(scip, &temp1vals);
1138 }
1139
1140 return SCIP_OKAY;
1141 }
1142
1143 /** projection of Lagrangian multipliers onto L2-norm ball */
1144 static
1145 void l2BallProjection(
1146 SCIP* scip, /**< SCIP data structure */
1147 SCIP_Real* dualvector, /**< Lagrangian multipliers to be projected onto L2-norm vall */
1148 int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1149 SCIP_Real radius /**< radius of the L2-norm ball */
1150 )
1151 {
1152 SCIP_Real l2norm;
1153 SCIP_Real factor;
1154
1155 assert(!SCIPisNegative(scip, radius));
1156
1157 l2norm = 0.0;
1158 /* calculate the L2-norm of the Lagrangian multipliers */
1159 for( int i = 0; i < dualvectorlen; i++ )
1160 {
1161 l2norm += SQR(dualvector[i]);
1162 }
1163 l2norm = SQRT(l2norm);
1164 factor = radius/(1.0 + l2norm);
1165
1166 /* if the vector of Lagrangian multipliers is outside the L2-norm ball, then do the projection */
1167 if( SCIPisGT(scip, l2norm, radius) && SCIPisLT(scip, factor, 1.0) )
1168 {
1169 for( int i = 0; i < dualvectorlen; i++ )
1170 {
1171 dualvector[i] *= factor;
1172 }
1173 }
1174 }
1175
1176 /** projection of Lagrangian multipliers onto L_infinity-norm ball */
1177 static
1178 void linfBallProjection(
1179 SCIP* scip, /**< SCIP data structure */
1180 SCIP_Real* dualvector, /**< Lagrangian multipliers to be projected onto L_infinity-norm vall */
1181 int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1182 SCIP_Real radius /**< radius of the L_infinity-norm ball */
1183 )
1184 {
1185 assert(!SCIPisNegative(scip, radius));
1186
1187 /* if the vector of Lagrangian multipliers is outside the L_infinity-norm ball, then do the projection */
1188 for( int i = 0; i < dualvectorlen; i++ )
1189 {
1190 if( SCIPisLT(scip, dualvector[i], -radius) )
1191 {
1192 dualvector[i] = -radius;
1193 }
1194 else if( SCIPisGT(scip, dualvector[i], radius) )
1195 {
1196 dualvector[i] = radius;
1197 }
1198 }
1199 }
1200
1201 /** weighted Lagrangian multipliers based on a given vector as stability center */
1202 /* @todo calculate weight outside this function and pass it (so that this function becomes generic and independent of
1203 * the terminology related to best Lagrangian multipliers)
1204 */
1205 static
1206 SCIP_RETCODE weightedDualVector(
1207 SCIP* scip, /**< SCIP data structure */
1208 SCIP_SEPADATA* sepadata, /**< separator data structure */
1209 SCIP_Real* dualvector, /**< Lagrangian multipliers */
1210 int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1211 SCIP_Real* stabilitycenter, /**< stability center (i.e., core vector of Lagrangian multipliers) */
1212 int stabilitycenterlen, /**< length of the stability center */
1213 int nbestdualupdates, /**< number of best Lagrangian values found so far */
1214 int totaliternum /**< total number of iterations of the relax-and-cut algorithm performed so far */
1215 )
1216 {
1217 SCIP_Real constant;
1218 SCIP_Real weight;
1219 SCIP_Real alpha;
1220
1221 constant = MAX(2.0, sepadata->constant);
1222 /* weight factor from the literature on Dantzig-Wolfe decomposition stabilization schemes */
1223 weight = MIN(constant, (totaliternum + 1 + nbestdualupdates) / 2.0);
1224 alpha = 1.0 / weight;
1225
1226 assert(dualvectorlen >= stabilitycenterlen);
1227
1228 /* weighted Lagrangian multipliers */
1229 for( int i = 0; i < stabilitycenterlen; i++ )
1230 {
1231 dualvector[i] = alpha * dualvector[i] + (1 - alpha) * stabilitycenter[i];
1232 }
1233 for( int i = stabilitycenterlen; i < dualvectorlen; i++ )
1234 {
1235 dualvector[i] = alpha * dualvector[i];
1236 }
1237
1238 return SCIP_OKAY;
1239 }
1240
1241 /** stabilize Lagrangian multipliers */
1242 static
1243 SCIP_RETCODE stabilizeDualVector(
1244 SCIP* scip, /**< SCIP data structure */
1245 SCIP_SEPADATA* sepadata, /**< separator data structure */
1246 SCIP_Real* dualvector, /**< Lagrangian multipliers */
1247 int dualvectorlen, /**< length of the Lagrangian multipliers vector */
1248 SCIP_Real* bestdualvector, /**< best Lagrangian multipliers found so far */
1249 int bestdualvectorlen, /**< length of the best Lagrangian multipliers vector */
1250 int nbestdualupdates, /**< number of best Lagrangian values found so far */
1251 int subgradientiternum, /**< iteration number of the subgradient algorithm */
1252 int totaliternum, /**< total number of iterations of the relax-and-cut algorithm performed so far */
1253 SCIP_Real maxviolscore, /**< weighted average of maximum value of generated cut violations and maximum value of
1254 complementarity slackness violations, in the current iteration */
1255 SCIP_Real maxviolscoreold, /**< weighted average of maximum value of generated cut violations and maximum value of
1256 complementarity slackness violations, in the previous iteration */
1257 SCIP_Real nviolscore, /**< weighted average of number of generated cut violations and number of complementarity
1258 slackness violations, in the current iteration */
1259 SCIP_Real nviolscoreold, /**< weighted average of number of generated cut violations and number of complementarity
1260 slackness violations, in the previous iteration */
1261 int nlpiters, /**< number of LP iterations taken for solving the Lagrangian dual with fixed multipliers in
1262 current iteration */
1263 SCIP_Real* ballradius /**< norm ball radius */
1264 )
1265 {
1266 if( sepadata->projectiontype > 0 )
1267 {
1268 if( subgradientiternum >= 1 )
1269 {
1270 /* update the ball radius */
1271 updateBallRadius(scip, sepadata, maxviolscore, maxviolscoreold, nviolscore, nviolscoreold, nlpiters,
1272 ballradius);
1273 }
1274
1275 if( sepadata->projectiontype == 1 )
1276 {
1277 /* projection of Lagrangian multipliers onto L1-norm ball */
1278 SCIP_CALL( l1BallProjection(scip, dualvector, dualvectorlen, *ballradius) );
1279 }
1280 else if( sepadata->projectiontype == 2 )
1281 {
1282 /* projection of Lagrangian multipliers onto L2-norm ball */
1283 l2BallProjection(scip, dualvector, dualvectorlen, *ballradius);
1284 }
1285 else if( sepadata->projectiontype == 3 )
1286 {
1287 /* projection of Lagrangian multipliers onto L_inf-norm ball */
1288 linfBallProjection(scip, dualvector, dualvectorlen, *ballradius);
1289 }
1290 }
1291
1292 if( sepadata->stabilitycentertype == 1 )
1293 {
1294 /* weighted Lagrangian multipliers based on best Langrangian multipliers as stability center */
1295 SCIP_CALL( weightedDualVector(scip, sepadata, dualvector, dualvectorlen, bestdualvector,
1296 bestdualvectorlen, nbestdualupdates, totaliternum) );
1297 }
1298
1299 return SCIP_OKAY;
1300 }
1301
1302 /** update Lagrangian multipliers */
1303 static
1304 SCIP_RETCODE updateDualVector(
1305 SCIP* scip, /**< SCIP data structure */
1306 SCIP_SEPADATA* sepadata, /**< separator data structure */
1307 SCIP_Real* dualvector1, /**< Lagrangian multipliers vector to be updated */
1308 SCIP_Real* dualvector2, /**< Lagrangian multipliers vector used for backtracking */
1309 int dualvector2len, /**< length of the Lagrangian multipliers vector used for backtracking */
1310 int ndualvector2updates,/**< number of best Lagrangian values found so far */
1311 int subgradientiternum, /**< iteration number of the subgradient algorithm */
1312 int totaliternum, /**< total number of iterations of the relax-and-cut algorithm performed so far */
1313 SCIP_Real steplength, /**< step length used for updating Lagrangian multipliers */
1314 SCIP_Real* subgradient, /**< subgradient vector */
1315 int ncuts, /**< number of generated cuts so far */
1316 SCIP_Bool backtrack, /**< whether the Lagrangian multipliers need to be backtracked */
1317 SCIP_Real maxviolscore, /**< weighted average of maximum value of generated cut violations and maximum value of
1318 complementarity slackness violations, in the current iteration */
1319 SCIP_Real maxviolscoreold, /**< weighted average of maximum value of generated cut violations and maximum value of
1320 complementarity slackness violations, in the previous iteration */
1321 SCIP_Real nviolscore, /**< weighted average of number of generated cut violations and number of complementarity
1322 slackness violations, in the current iteration */
1323 SCIP_Real nviolscoreold, /**< weighted average of number of generated cut violations and number of complementarity
1324 slackness violations, in the previous iteration */
1325 int nlpiters, /**< number of LP iterations taken for solving the Lagrangian dual with fixed multipliers in
1326 current iteration */
1327 SCIP_Bool* dualvecsdiffer, /**< whether the updated Lagrangian multipliers differ from the old one */
1328 SCIP_Real* ballradius /**< norm ball radius */
1329 )
1330 {
1331 SCIP_Real* dualvector1copy;
1332
1333 assert(dualvector2len <= ncuts);
1334 assert(dualvecsdiffer != NULL);
1335
1336 *dualvecsdiffer = FALSE;
1337 /* @todo do allocation and free operations outside only once instead of every time this function is called? */
1338 /* copy of the Lagrangian multipliers to be used to check if the updated vector is different than this */
1339 SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualvector1copy, ncuts) );
1340 for( int i = 0; i < ncuts; i++ )
1341 {
1342 dualvector1copy[i] = dualvector1[i];
1343 }
1344
1345 /* if backtracking was not identified at the time of the mu parameter update, then update the Lagrangian multipliers
1346 * based on the given subgradient vector
1347 */
1348 if( !backtrack )
1349 {
1350 assert((subgradient != NULL) || (ncuts == 0));
1351 assert(subgradientiternum >= 0);
1352
1353 /* update Lagrangian multipliers */
1354 for( int i = 0; i < ncuts; i++ )
1355 {
1356 dualvector1[i] += steplength * subgradient[i];
1357 }
1358
1359 /* projection onto non-negative orthant */
1360 for( int i = 0; i < ncuts; i++ )
1361 {
1362 dualvector1[i] = MAX(dualvector1[i], 0.0);
1363 }
1364
1365 /* stabilization of Lagrangian multipliers */
1366 SCIP_CALL( stabilizeDualVector(scip, sepadata, dualvector1, ncuts, dualvector2, dualvector2len,
1367 ndualvector2updates, subgradientiternum, totaliternum, maxviolscore, maxviolscoreold, nviolscore,
1368 nviolscoreold, nlpiters, ballradius) );
1369
1370 /* projection onto non-negative orthant again in case stabilization changed some components negative*/
1371 for( int i = 0; i < ncuts; i++ )
1372 {
1373 dualvector1[i] = MAX(dualvector1[i], 0.0);
1374 }
1375 }
1376 /* if backtracking was identified at the time of the mu parameter update, then backtrack the Lagrangian multipliers
1377 * based on the given backtracking multipliers
1378 */
1379 else
1380 {
1381 for( int i = 0; i < dualvector2len; i++ )
1382 {
1383 dualvector1[i] = dualvector2[i];
1384 }
1385
1386 for( int i = dualvector2len; i < ncuts; i++ )
1387 {
1388 dualvector1[i] = 0.0;
1389 }
1390 }
1391
1392 /* identify if the vector of Lagrangian multipliers is indeed different after updating */
1393 for( int i = 0; i < ncuts; i++ )
1394 {
1395 if( !SCIPisEQ(scip, dualvector1[i], dualvector1copy[i]) )
1396 {
1397 *dualvecsdiffer = TRUE;
1398 break;
1399 }
1400 }
1401
1402 /* free memory */
1403 for( int i = 0; i < ncuts; i++ )
1404 {
1405 dualvector1copy[i] = 0.0;
1406 }
1407 SCIPfreeCleanBufferArray(scip, &dualvector1copy);
1408
1409 return SCIP_OKAY;
1410 }
1411
1412 /** check different termination criteria */
1413 /* @note: the criterion based on objvecsdiffer assumes deterministic solving process (i.e., we would get same LP solution
1414 * for "Lagrangian dual with fixed Lagrangian multipliers" when the objective vector remains the same across iterations).
1415 */
1416 /* @todo nlpssolved criterion? */
1417 static
1418 SCIP_RETCODE checkLagrangianDualTermination(
1419 SCIP_SEPADATA* sepadata, /**< separator data structure */
1420 int nnewaddedsoftcuts, /**< number of cuts that were recently penalized and added to the Lagrangian dual's
1421 objective function */
1422 int nyettoaddsoftcuts, /**< number of cuts that are yet to be penalized and added to the Lagrangian dual's
1423 objective function */
1424 SCIP_Bool objvecsdiffer, /**< whether the Lagrangian dual's objective function has changed */
1425 int ngeneratedcurrroundcuts, /**< number of cuts generated in the current separation round */
1426 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
1427 int ncurrroundlpiters, /**< number of separating LP iterations in the current separation round */
1428 int depth, /**< depth of the current node */
1429 SCIP_Bool* terminate /**< whether to terminate the subgradient algorithm loop */
1430 )
1431 {
1432 *terminate = FALSE;
1433
1434 /* check if no new cuts were added to the Lagrangian dual, no cuts are remaining to be added, and the objective
1435 * function of the Lagrangian dual with fixed multipliers had not changed from the previous iteration
1436 */
1437 if( (nnewaddedsoftcuts == 0) && (nyettoaddsoftcuts == 0) && !objvecsdiffer )
1438 *terminate = TRUE;
1439
1440 /* check if allowed number of cuts in this separation round have already been generated */
1441 if( ngeneratedcurrroundcuts >= nmaxgeneratedperroundcuts )
1442 *terminate = TRUE;
1443
1444 /* check if allowed number of cuts in the tree have already been generated */
1445 if( sepadata->ntotalcuts >= sepadata->nmaxtotalcuts )
1446 *terminate = TRUE;
1447
1448 /* check if allowed number of simplex iterations in this separation round have already been used up */
1449 if( (sepadata->nmaxperroundlpiters >= 0) && (ncurrroundlpiters >= sepadata->nmaxperroundlpiters) )
1450 *terminate = TRUE;
1451
1452 /* check if allowed number of simplex iterations in the root node have already been used up */
1453 if( (depth == 0) && (sepadata->nmaxrootlpiters >= 0) && (sepadata->nrootlpiters >= sepadata->nmaxrootlpiters) )
1454 *terminate = TRUE;
1455
1456 /* check if allowed number of simplex iterations in the tree have already been used up */
1457 if( (sepadata->nmaxtotallpiters >= 0) && (sepadata->ntotallpiters >= sepadata->nmaxtotallpiters) )
1458 *terminate = TRUE;
1459
1460 return SCIP_OKAY;
1461 }
1462
1463 /** solve the LP corresponding to the Lagrangian dual with fixed Lagrangian multipliers */
1464 static
1465 SCIP_RETCODE solveLagromoryLP(
1466 SCIP* scip, /**< SCIP data structure */
1467 SCIP_SEPADATA* sepadata, /**< separator data structure */
1468 int depth, /**< depth of the current node in the tree */
1469 SCIP_Real origobjoffset, /**< objective offset in the current node's relaxation */
1470 SCIP_Bool* solfound, /**< whether an LP optimal solution has been found */
1471 SCIP_SOL* sol, /**< data structure to store LP optimal solution, if found */
1472 SCIP_Real* solvals, /**< values of the LP optimal solution, if found */
1473 SCIP_Real* objval, /**< optimal objective value of the LP optimal solution, if found */
1474 int* ncurrroundlpiters /**< number of LP iterations taken for solving Lagrangian dual problems with fixed multipliers
1475 in the current separator round */
1476 )
1477 {
1478 SCIP_Real timelimit;
1479 SCIP_COL** cols;
1480 SCIP_COL* col;
1481 SCIP_VAR* var;
1482 SCIP_LPI* lpi;
1483 SCIP_Bool lperror;
1484 SCIP_Bool cutoff;
1485 SCIP_LPSOLSTAT stat;
1486 SCIP_Longint ntotallpiters;
1487 SCIP_Longint nlpiters;
1488 int ncols;
1489 int iterlimit;
1490
1491 assert(solfound != NULL);
1492 assert(sol != NULL);
1493 assert(solvals != NULL);
1494 assert(ncurrroundlpiters != NULL);
1495 assert(objval != NULL);
1496
1497 *solfound = FALSE;
1498 lperror = FALSE;
1499 cutoff = FALSE;
1500 iterlimit = -1;
1501 lpi = sepadata->lpiwithsoftcuts;
1502
1503 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
1504
1505 /* set time limit */
1506 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
1507 if( !SCIPisInfinity(scip, timelimit) )
1508 {
1509 timelimit -= SCIPgetSolvingTime(scip);
1510 if( timelimit <= 0.0 )
1511 {
1512 SCIPdebugMsg(scip, "skip Lagromory cut generation since no time left\n");
1513 goto TERMINATE;
1514 }
1515 /* @note: the following direct LPI call is being used because of the lack of an equivalent function call in
1516 * scip_lp.c (lpSetRealpar exists in lp.c though)
1517 */
1518 SCIP_CALL( SCIPlpiSetRealpar(lpi, SCIP_LPPAR_LPTILIM, timelimit) );
1519 }
1520
1521 /* find iteration limit */
1522 if( (depth == 0) &&
1523 (sepadata->perrootlpiterfactor >= 0.0 && !SCIPisInfinity(scip, sepadata->perrootlpiterfactor)) )
1524 {
1525 iterlimit = (int)(sepadata->perrootlpiterfactor * SCIPgetNRootFirstLPIterations(scip));
1526 }
1527 else if( (depth > 0) &&
1528 (sepadata->perlpiterfactor >= 0.0 && !SCIPisInfinity(scip, sepadata->perlpiterfactor)) )
1529 {
1530 iterlimit = (int)(sepadata->perlpiterfactor * SCIPgetNNodeInitLPIterations(scip));
1531 }
1532 if( sepadata->nmaxperroundlpiters >= 0 )
1533 {
1534 if( sepadata->nmaxperroundlpiters - *ncurrroundlpiters >= 0 )
1535 {
1536 if( iterlimit >= 0 )
1537 {
1538 iterlimit = MIN(iterlimit, sepadata->nmaxperroundlpiters - *ncurrroundlpiters);
1539 }
1540 else
1541 {
1542 iterlimit = sepadata->nmaxperroundlpiters - *ncurrroundlpiters;
1543 }
1544 }
1545 else
1546 {
1547 iterlimit = 0;
1548 }
1549 }
1550 /* @todo impose a finite iteration limit only when the dualvector changes from zero to non-zero for the first time because
1551 * many simplex pivots are performed in this case even with warm starting (compared to the case when the
1552 * dualvector changes from non-zero to non-zero).
1553 */
1554
1555 /* solve the LP with an iteration limit and get number of simplex iterations taken */
1556 ntotallpiters = SCIPgetNLPIterations(scip);
1557
1558 SCIP_CALL( SCIPsolveDiveLP(scip, iterlimit, &lperror, &cutoff) );
1559
1560 nlpiters = SCIPgetNLPIterations(scip) - ntotallpiters;
1561
1562 /* get the solution and objective value if optimal */
1563 stat = SCIPgetLPSolstat(scip);
1564 /* @todo is there any way to accept terminations due to iterlimit and timelimit as well? It is not possible
1565 * currently because primal sol is not saved in these cases.
1566 */
1567 /* @note: ideally, only primal feasibility is sufficient. But, there is no such option with SCIPgetLPSolstat. */
1568 if( stat == SCIP_LPSOLSTAT_OPTIMAL )
1569 {
1570 if( SCIPisLPSolBasic(scip) )
1571 {
1572 *solfound = TRUE;
1573
1574 /* update sol */
1575 for( int i = 0; i < ncols; ++i )
1576 {
1577 col = cols[i];
1578 assert(col != NULL);
1579
1580 var = SCIPcolGetVar(col);
1581
1582 solvals[i] = SCIPcolGetPrimsol(col);
1583 SCIP_CALL( SCIPsetSolVal(scip, sol, var, solvals[i]) );
1584 }
1585
1586 *objval = SCIPgetLPObjval(scip);
1587 *objval = *objval + origobjoffset;
1588 }
1589 }
1590
1591 /* update some statistics */
1592 if( depth == 0 )
1593 {
1594 sepadata->nrootlpiters += (int)nlpiters;
1595 }
1596 sepadata->ntotallpiters += (int)nlpiters;
1597 *ncurrroundlpiters += (int)nlpiters;
1598
1599 TERMINATE:
1600 return SCIP_OKAY;
1601 }
1602
1603 /** solve the LP corresponding to the node relaxation upon adding all the generated cuts */
1604 static
1605 SCIP_RETCODE solveLPWithHardCuts(
1606 SCIP* scip, /**< SCIP data structure */
1607 SCIP_SEPADATA* sepadata, /**< separator data structure */
1608 SCIP_Bool* solfound, /**< whether an LP optimal solution has been found */
1609 SCIP_SOL* sol, /**< data structure to store LP optimal solution, if found */
1610 SCIP_Real* solvals /**< values of the LP optimal solution, if found */
1611 )
1612 {
1613 SCIP_Real timelimit;
1614 SCIP_COL** cols;
1615 SCIP_COL* col;
1616 SCIP_VAR* var;
1617 int ncols;
1618
1619 assert(solfound != NULL);
1620 assert(sol != NULL);
1621 assert(solvals != NULL);
1622
1623 *solfound = FALSE;
1624
1625 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
1626
1627 /* set time limit */
1628 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
1629 if( !SCIPisInfinity(scip, timelimit) )
1630 {
1631 timelimit -= SCIPgetSolvingTime(scip);
1632 if( timelimit <= 0.0 )
1633 {
1634 SCIPdebugMsg(scip, "skip Lagromory cut generation since no time left\n");
1635 goto TERMINATE;
1636 }
1637 SCIP_CALL( SCIPlpiSetRealpar(sepadata->lpiwithhardcuts, SCIP_LPPAR_LPTILIM, timelimit) );
1638 }
1639
1640 /* solve the LP */
1641 SCIPlpiSolvePrimal(sepadata->lpiwithhardcuts); /*lint !e534*/
1642
1643 /* get the solution if primal feasible */
1644 if( SCIPlpiIsPrimalFeasible(sepadata->lpiwithhardcuts) )
1645 {
1646 *solfound = TRUE;
1647 SCIP_CALL( SCIPlpiGetSol(sepadata->lpiwithhardcuts, NULL, solvals, NULL, NULL, NULL) );
1648
1649 /* update sol */
1650 for( int i = 0; i < ncols; ++i )
1651 {
1652 col = cols[i];
1653 assert(col != NULL);
1654
1655 var = SCIPcolGetVar(col);
1656
1657 SCIP_CALL( SCIPsetSolVal(scip, sol, var, solvals[i]) );
1658 }
1659 }
1660
1661 TERMINATE:
1662 return SCIP_OKAY;
1663 }
1664
1665 /** construct a cut based on the input cut coefficients, sides, etc */
1666 static
1667 SCIP_RETCODE constructCutRow(
1668 SCIP* scip, /**< SCIP data structure */
1669 SCIP_SEPA* sepa, /**< pointer to the separator */
1670 SCIP_SEPADATA* sepadata, /**< separator data structure */
1671 int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */
1672 int subgradientiternum, /**< iteration number of the subgradient algorithm */
1673 int cutnnz, /**< number of nonzeros in cut */
1674 int* cutinds, /**< column indices in cut */
1675 SCIP_Real* cutcoefs, /**< cut cofficients */
1676 SCIP_Real cutefficacy, /**< cut efficacy */
1677 SCIP_Real cutrhs, /**< RHS of cut */
1678 SCIP_Bool cutislocal, /**< whether cut is local */
1679 int cutrank, /**< rank of cut */
1680 SCIP_ROW** generatedcuts, /**< array of generated cuts */
1681 SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut
1682 generations */
1683 int ngeneratedcurrroundcuts, /**< number of cuts generated until the previous basis in the current separation round */
1684 int* ngeneratednewcuts, /**< number of new cuts generated using the current basis */
1685 SCIP_Bool* cutoff /**< should the current node be cutoff? */
1686 )
1687 {
1688 SCIP_COL** cols;
1689 SCIP_Real minactivity;
1690 SCIP_Real maxactivity;
1691 SCIP_Real lhs;
1692 SCIP_Real rhs;
1693
1694 assert(scip != NULL);
1695 assert(mainiternum >= 0);
1696 assert(ngeneratednewcuts != NULL);
1697 assert(cutoff != NULL);
1698
1699 *cutoff = FALSE;
1700
1701 if( cutnnz == 0 && SCIPisFeasNegative(scip, cutrhs) ) /*lint !e644*/
1702 {
1703 SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", cutrhs);
1704 *cutoff = TRUE;
1705 return SCIP_OKAY;
1706 }
1707
1708 /* only take efficient cuts */
1709 if( SCIPisEfficacious(scip, cutefficacy) )
1710 {
1711 SCIP_ROW* cut;
1712 SCIP_VAR* var;
1713 char cutname[SCIP_MAXSTRLEN];
1714 int v;
1715
1716 /* construct cut name */
1717 if( subgradientiternum >= 0 )
1718 {
1719 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_%d" "_%d" "_%d", SCIPsepaGetName(sepa),
1720 sepadata->ncalls, mainiternum, subgradientiternum, ngeneratedcurrroundcuts + *ngeneratednewcuts);
1721 }
1722 else
1723 {
1724 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_%d" "_%d", SCIPsepaGetName(sepa),
1725 sepadata->ncalls, mainiternum, ngeneratedcurrroundcuts + *ngeneratednewcuts);
1726 }
1727
1728 /* create empty cut */
1729 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE,
1730 sepadata->dynamiccuts) );
1731
1732 /* set cut rank */
1733 SCIProwChgRank(cut, cutrank); /*lint !e644*/
1734
1735 /* cache the row extension and only flush them if the cut gets added */
1736 SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
1737
1738 cols = SCIPgetLPCols(scip);
1739
1740 /* collect all non-zero coefficients */
1741 for( v = 0; v < cutnnz; ++v )
1742 {
1743 var = SCIPcolGetVar(cols[cutinds[v]]);
1744 SCIP_CALL( SCIPaddVarToRow(scip, cut, var, cutcoefs[v]) );
1745 }
1746
1747 /* flush all changes before adding the cut */
1748 SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
1749
1750 if( SCIProwGetNNonz(cut) == 0 )
1751 {
1752 assert( SCIPisFeasNegative(scip, cutrhs) );
1753 SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", cutrhs);
1754 *cutoff = TRUE;
1755 return SCIP_OKAY;
1756 }
1757 else
1758 {
1759 /* gathering lhs and rhs again in case the separator is extended later to add other cuts/constraints that may
1760 * have non-inf lhs or inf rhs */
1761 lhs = SCIProwGetLhs(cut);
1762 rhs = SCIProwGetRhs(cut);
1763 assert(SCIPisInfinity(scip, -lhs));
1764 assert(!SCIPisInfinity(scip, rhs));
1765
1766 SCIPdebugMsg(scip, " -> %s cut <%s>: rhs=%f, eff=%f\n", "lagromory", cutname, cutrhs, cutefficacy);
1767
1768 /* check if the cut leads to infeasibility (i.e., *cutoff = TRUE) */
1769 {
1770 /* modifiable cuts cannot be declared infeasible, since we don't know all coefficients */
1771 if( SCIProwIsModifiable(cut) )
1772 *cutoff = FALSE;
1773
1774 /* check for activity infeasibility */
1775 minactivity = SCIPgetRowMinActivity(scip, cut);
1776 maxactivity = SCIPgetRowMaxActivity(scip, cut);
1777
1778 if( (!SCIPisInfinity(scip, rhs) && SCIPisFeasPositive(scip, minactivity - rhs)) ||
1779 (!SCIPisInfinity(scip, -lhs) && SCIPisFeasNegative(scip, maxactivity - lhs)) )
1780 {
1781 SCIPdebugMsg(scip, "cut <%s> is infeasible (sides=[%g,%g], act=[%g,%g])\n",
1782 SCIProwGetName(cut), lhs, rhs, minactivity, maxactivity);
1783 *cutoff = TRUE;
1784 }
1785 }
1786
1787 /* store the newly generated cut in an array and update some statistics */
1788 generatedcuts[ngeneratedcurrroundcuts + *ngeneratednewcuts] = cut;
1789 generatedcutefficacies[ngeneratedcurrroundcuts + *ngeneratednewcuts] = cutefficacy;
1790 ++(*ngeneratednewcuts);
1791 }
1792 }
1793
1794 return SCIP_OKAY;
1795 }
1796
1797 /** aggregated generated cuts based on the best Lagrangian multipliers */
1798 static
1799 SCIP_RETCODE aggregateGeneratedCuts(
1800 SCIP* scip, /**< SCIP data structure */
1801 SCIP_SEPA* sepa, /**< pointer to the separator */
1802 SCIP_SEPADATA* sepadata, /**< separator data structure */
1803 SCIP_ROW** generatedcuts, /**< cuts generated in the current separation round */
1804 SCIP_Real* bestdualvector, /**< best Lagrangian multipliers vector */
1805 int bestdualvectorlen, /**< length of the best Lagrangian multipliers vector */
1806 SCIP_ROW** aggrcuts, /**< aggregated cuts generated so far in the current separation round */
1807 int* naggrcuts, /**< number of aggregated cuts generated so far in the current separation round */
1808 SCIP_Bool* cutoff /**< should the current node be cutoff? */
1809 )
1810 {
1811 SCIP_Real* cutvals; /**< cut cofficients */
1812 SCIP_COL** cutcols;
1813 SCIP_COL** cols;
1814 SCIP_VAR* var;
1815 SCIP_ROW* cut;
1816 SCIP_ROW* aggrcut;
1817 SCIP_Real minactivity;
1818 SCIP_Real maxactivity;
1819 SCIP_Real cutlhs;
1820 SCIP_Real cutrhs;
1821 SCIP_Real cutconst;
1822 SCIP_Real aggrcutlhs;
1823 SCIP_Real aggrcutrhs;
1824 SCIP_Real aggrcutconst;
1825 SCIP_Real* aggrcutvals;
1826 SCIP_Real* aggrcutcoefs;
1827 SCIP_Real multiplier;
1828 SCIP_Bool aggrcutislocal;
1829 SCIP_Bool aggrindicator;
1830 SCIP_Real* tmpcutvals;
1831 SCIP_Real QUAD(quadterm);
1832 SCIP_Real QUAD(tmpcutconst);
1833 SCIP_Real QUAD(tmpcutrhs);
1834 SCIP_Real QUAD(quadprod);
1835 char aggrcutname[SCIP_MAXSTRLEN];
1836 int cutnnz; /**< number of nonzeros in cut */
1837 int aggrcutnnz;
1838 int* aggrcutinds; /**< column indices in cut */
1839 int aggrcutrank; /**< rank of cut */
1840 int cutrank;
1841 int ncols;
1842 int collppos;
1843 int nlocalcuts;
1844
1845 assert(scip != NULL);
1846 assert(naggrcuts != NULL);
1847 assert(cutoff != NULL);
1848
1849 *cutoff = FALSE;
1850 aggrcutlhs = -SCIPinfinity(scip);
1851 aggrcutnnz = 0;
1852 aggrcutrank = -1;
1853 nlocalcuts = 0;
1854 aggrindicator = FALSE;
1855
1856 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
1857
1858 /* allocate memory */
1859 SCIP_CALL( SCIPallocCleanBufferArray(scip, &tmpcutvals, QUAD_ARRAY_SIZE(ncols)) );
1860 QUAD_ASSIGN(tmpcutconst, 0.0);
1861 QUAD_ASSIGN(tmpcutrhs, 0.0);
1862 SCIP_CALL( SCIPallocCleanBufferArray(scip, &aggrcutvals, ncols) );
1863 SCIP_CALL( SCIPallocBufferArray(scip, &aggrcutcoefs, ncols) );
1864 SCIP_CALL( SCIPallocBufferArray(scip, &aggrcutinds, ncols) );
1865
1866 /* aggregate cuts based on the input Lagrangian multipliers */
1867 for( int i = 0; i < bestdualvectorlen; i++ )
1868 {
1869 multiplier = bestdualvector[i];
1870 if( SCIPisGE(scip, multiplier, 1e-4) )
1871 {
1872 cut = generatedcuts[i];
1873 cutnnz = SCIProwGetNNonz(cut);
1874 cutlhs = SCIProwGetLhs(cut);
1875 cutrhs = SCIProwGetRhs(cut);
1876 assert(SCIPisInfinity(scip, -cutlhs));
1877 assert(!SCIPisInfinity(scip, cutrhs));
1878 cutvals = SCIProwGetVals(cut);
1879 cutcols = SCIProwGetCols(cut);
1880 cutconst = SCIProwGetConstant(cut);
1881
1882 for( int j = 0; j < cutnnz; j++ )
1883 {
1884 collppos = SCIPcolGetLPPos(cutcols[j]);
1885 assert(collppos >= 0);
1886 assert(collppos <= ncols);
1887
1888 QUAD_ARRAY_LOAD(quadterm, tmpcutvals, collppos);
1889 SCIPquadprecProdDD(quadprod, multiplier, cutvals[j]);
1890 SCIPquadprecSumQQ(quadterm, quadterm, quadprod);
1891 QUAD_ARRAY_STORE(tmpcutvals, collppos, quadterm);
1892 }
1893
1894 SCIPquadprecProdDD(quadprod, multiplier, cutconst);
1895 SCIPquadprecSumQQ(tmpcutconst, tmpcutconst, quadprod);
1896 SCIPquadprecProdDD(quadprod, multiplier, cutrhs);
1897 SCIPquadprecSumQQ(tmpcutrhs, tmpcutrhs, quadprod);
1898
1899 cutrank = SCIProwGetRank(cut);
1900 aggrcutrank = MAX(aggrcutrank, cutrank);
1901 nlocalcuts += (SCIProwIsLocal(cut) ? 1 : 0);
1902 aggrindicator = TRUE;
1903 }
1904 }
1905 /* if at least one cut is local, then the aggregated cut is local */
1906 aggrcutislocal = (nlocalcuts > 0 ? TRUE : FALSE);
1907
1908 /* if the aggregation was successful, then create an empty row and build a cut */
1909 if( aggrindicator )
1910 {
1911 aggrcutconst = QUAD_TO_DBL(tmpcutconst);
1912 aggrcutrhs = QUAD_TO_DBL(tmpcutrhs);
1913
1914 /* build sparse representation of the aggregated cut */
1915 for( int i = 0; i < ncols; i++ )
1916 {
1917 QUAD_ARRAY_LOAD(quadterm, tmpcutvals, i);
1918 aggrcutvals[i] = QUAD_TO_DBL(quadterm);
1919 if( !SCIPisZero(scip, aggrcutvals[i]) )
1920 {
1921 aggrcutcoefs[aggrcutnnz] = aggrcutvals[i];
1922 aggrcutinds[aggrcutnnz] = i;
1923 aggrcutnnz++;
1924 }
1925 }
1926
1927 if( aggrcutnnz == 0 && SCIPisFeasNegative(scip, aggrcutrhs) ) /*lint !e644*/
1928 {
1929 SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", aggrcutrhs);
1930 *cutoff = TRUE;
1931 goto TERMINATE;
1932 }
1933
1934 /* construct aggregated cut name */
1935 (void) SCIPsnprintf(aggrcutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_aggregated", SCIPsepaGetName(sepa),
1936 sepadata->ncalls);
1937
1938 /* create empty cut */
1939 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &aggrcut, sepa, aggrcutname, aggrcutlhs - aggrcutconst, aggrcutrhs -
1940 aggrcutconst, aggrcutislocal, FALSE, sepadata->dynamiccuts) );
1941
1942 /* set cut rank */
1943 SCIProwChgRank(aggrcut, aggrcutrank); /*lint !e644*/
1944
1945 /* cache the row extension and only flush them if the cut gets added */
1946 SCIP_CALL( SCIPcacheRowExtensions(scip, aggrcut) );
1947
1948 /* collect all non-zero coefficients */
1949 for( int i = 0; i < aggrcutnnz; i++ )
1950 {
1951 var = SCIPcolGetVar(cols[aggrcutinds[i]]);
1952 SCIP_CALL( SCIPaddVarToRow(scip, aggrcut, var, aggrcutcoefs[i]) );
1953 }
1954
1955 /* flush all changes before adding the cut */
1956 SCIP_CALL( SCIPflushRowExtensions(scip, aggrcut) );
1957
1958 if( SCIProwGetNNonz(aggrcut) == 0 )
1959 {
1960 assert( SCIPisFeasNegative(scip, aggrcutrhs) );
1961 SCIPdebugMsg(scip, " -> Lagromory cut detected node infeasibility with cut 0 <= %g.\n", aggrcutrhs);
1962 *cutoff = TRUE;
1963 goto TERMINATE;
1964 }
1965 else
1966 {
1967 /* gathering lhs and rhs again in case the separator is extended later to add other cuts/constraints that may
1968 * have non-inf lhs or inf rhs */
1969 cutlhs = SCIProwGetLhs(aggrcut);
1970 cutrhs = SCIProwGetRhs(aggrcut);
1971 assert(SCIPisInfinity(scip, -cutlhs));
1972 assert(!SCIPisInfinity(scip, cutrhs));
1973
1974 SCIPdebugMsg(scip, " -> %s cut <%s>: rhs=%f\n", "lagromory", aggrcutname, aggrcutrhs);
1975
1976 /* check if the cut leads to infeasibility (i.e., *cutoff = TRUE) */
1977 {
1978 /* modifiable cuts cannot be declared infeasible, since we don't know all coefficients */
1979 if( SCIProwIsModifiable(aggrcut) )
1980 *cutoff = FALSE;
1981
1982 /* check for activity infeasibility */
1983 minactivity = SCIPgetRowMinActivity(scip, aggrcut);
1984 maxactivity = SCIPgetRowMaxActivity(scip, aggrcut);
1985
1986 if( (!SCIPisInfinity(scip, cutrhs) && SCIPisFeasPositive(scip, minactivity - cutrhs)) ||
1987 (!SCIPisInfinity(scip, -cutlhs) && SCIPisFeasNegative(scip, maxactivity - cutlhs)) )
1988 {
1989 SCIPdebugMsg(scip, "cut <%s> is infeasible (sides=[%g,%g], act=[%g,%g])\n",
1990 SCIProwGetName(aggrcut), cutlhs, cutrhs, minactivity, maxactivity);
1991 *cutoff = TRUE;
1992 }
1993 }
1994
1995 /* add the aggregated cut to a separate data structure */
1996 aggrcuts[*naggrcuts] = aggrcut;
1997 (*naggrcuts)++;
1998 }
1999
2000 QUAD_ASSIGN(quadterm, 0.0);
2001 for( int i = 0; i < ncols; i++ )
2002 {
2003 aggrcutvals[i] = 0.0;
2004 QUAD_ARRAY_STORE(tmpcutvals, i, quadterm);
2005 }
2006 }
2007
2008 TERMINATE:
2009 /* free memory */
2010 SCIPfreeBufferArray(scip, &aggrcutinds);
2011 SCIPfreeBufferArray(scip, &aggrcutcoefs);
2012 SCIPfreeCleanBufferArray(scip, &aggrcutvals);
2013 SCIPfreeCleanBufferArray(scip, &tmpcutvals);
2014
2015 return SCIP_OKAY;
2016 }
2017
2018 /** main method: LP solution separation method of separator */
2019 static
2020 SCIP_RETCODE generateGMICuts(
2021 SCIP* scip, /**< SCIP data structure */
2022 SCIP_SEPA* sepa, /**< pointer to the separator */
2023 SCIP_SEPADATA* sepadata, /**< separator data structure */
2024 int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */
2025 int subgradientiternum, /**< iteration number of the subgradient algorithm */
2026 SCIP_SOL* sol, /**< LP solution to be used for cut generation */
2027 SCIP_Real* solvals, /**< values of the LP solution to be used for cut generation */
2028 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
2029 SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */
2030 SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */
2031 SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut
2032 generations */
2033 int ngeneratedcurrroundcuts, /**< number of cuts generated until the previous basis in the current separation round */
2034 int* ngeneratednewcuts, /**< number of new cuts generated using the current basis */
2035 int depth, /**< depth of the current node in the tree */
2036 SCIP_Bool* cutoff /**< should the current node be cutoff? */
2037 )
2038 {
2039 SCIP_Real minfrac;
2040 SCIP_Real maxfrac;
2041 SCIP_Real frac;
2042 SCIP_Real* basisfrac;
2043 SCIP_Real* binvrow;
2044 SCIP_AGGRROW* aggrrow;
2045 SCIP_Bool success;
2046 SCIP_Real* cutcoefs;
2047 SCIP_Real cutrhs;
2048 SCIP_Real cutefficacy;
2049 SCIP_Bool cutislocal;
2050 SCIP_ROW** rows;
2051 SCIP_COL** cols;
2052 SCIP_VAR* var;
2053 SCIP_ROW* row;
2054 int cutrank;
2055 int cutnnz;
2056 int* cutinds;
2057 int* basisind;
2058 int* inds;
2059 int nrows;
2060 int ncols;
2061 int* basisperm;
2062 int k;
2063 int c;
2064 int ninds;
2065 int nmaxcutsperlp;
2066
2067 assert(ngeneratednewcuts != NULL);
2068
2069 minfrac = sepadata->away;
2070 maxfrac = 1.0 - sepadata->away;
2071 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
2072 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2073 *ngeneratednewcuts = 0;
2074 nmaxcutsperlp = ((depth == 0) ? sepadata->nmaxcutsperlproot : sepadata->nmaxcutsperlp);
2075
2076 /* allocate memory */
2077 SCIP_CALL( SCIPallocBufferArray(scip, &basisperm, nrows) );
2078 SCIP_CALL( SCIPallocBufferArray(scip, &basisfrac, nrows) );
2079 SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) );
2080 SCIP_CALL( SCIPallocBufferArray(scip, &inds, nrows) );
2081 SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
2082 SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, ncols) );
2083 SCIP_CALL( SCIPallocBufferArray(scip, &cutinds, ncols) );
2084 SCIP_CALL( SCIPaggrRowCreate(scip, &aggrrow) );
2085
2086 SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
2087
2088 /* check if the rows of the simplex tableau are suitable for cut generation and build an array of fractions */
2089 for( int i = 0; i < nrows; ++i )
2090 {
2091 frac = 0.0;
2092
2093 c = basisind[i];
2094
2095 basisperm[i] = i;
2096
2097 /* if the simplex tableau row corresponds to an LP column */
2098 if( c >= 0 )
2099 {
2100 assert(c < ncols);
2101
2102 var = SCIPcolGetVar(cols[c]);
2103 /* if the column is non-continuous one */
2104 if( SCIPvarGetType(var) != SCIP_VARTYPE_CONTINUOUS )
2105 {
2106 frac = SCIPfeasFrac(scip, solvals[c]);
2107 frac = MIN(frac, 1.0 - frac);
2108 }
2109 }
2110 /* if the simplex tableau row corresponds to an LP row and separation on rows is allowed */
2111 else if( sepadata->separaterows )
2112 {
2113 assert(0 <= -c-1);
2114 assert(-c-1 < nrows);
2115
2116 row = rows[-c-1];
2117 /* if the row is suitable for cut generation */
2118 if( SCIProwIsIntegral(row) && !SCIProwIsModifiable(row) )
2119 {
2120 frac = SCIPfeasFrac(scip, SCIPgetRowActivity(scip, row));
2121 frac = MIN(frac, 1.0 - frac);
2122 }
2123 }
2124
2125 if( frac >= minfrac )
2126 {
2127 /* slightly change fractionality to have random order for equal fractions */
2128 basisfrac[i] = frac + SCIPrandomGetReal(sepadata->randnumgen, -1e-6, 1e-6);
2129 }
2130 else
2131 {
2132 basisfrac[i] = 0.0;
2133 }
2134 }
2135
2136 /* if there is a need to sort the fractionalities */
2137 if( sepadata->sortcutoffsol )
2138 {
2139 /* sort basis indices by fractionality */
2140 SCIPsortDownRealInt(basisfrac, basisperm, nrows);
2141 }
2142
2143 /* for all basic columns belonging to integer variables, try to generate a GMI cut */
2144 for( int i = 0; i < nrows && !SCIPisStopped(scip) && !*cutoff; ++i )
2145 {
2146 if( (ngeneratedcurrroundcuts + *ngeneratednewcuts >= nmaxgeneratedperroundcuts) ||
2147 (sepadata->ntotalcuts + *ngeneratednewcuts >= sepadata->nmaxtotalcuts) ||
2148 (*ngeneratednewcuts >= nmaxcutsperlp) )
2149 break;
2150
2151 ninds = -1;
2152 cutefficacy = 0.0;
2153
2154 /* either break the loop or proceed to the next iteration if the fractionality is zero */
2155 if( basisfrac[i] == 0.0 )
2156 {
2157 if( sepadata->sortcutoffsol )
2158 break;
2159 else
2160 continue;
2161 }
2162
2163 k = basisperm[i];
2164
2165 /* get the row of B^-1 for this basic integer variable with fractional solution value and call aggregate function */
2166 SCIP_CALL( SCIPgetLPBInvRow(scip, k, binvrow, inds, &ninds) );
2167
2168 SCIP_CALL( SCIPaggrRowSumRows(scip, aggrrow, binvrow, inds, ninds, sepadata->sidetypebasis, allowlocal, 2,
2169 (int) MAXAGGRLEN(ncols), &success) );
2170
2171 if( !success )
2172 continue;
2173
2174 /* @todo Currently we are using the SCIPcalcMIR() function to compute the coefficients of the Gomory
2175 * cut. Alternatively, we could use the direct version (see thesis of Achterberg formula (8.4)) which
2176 * leads to cut a of the form \sum a_i x_i \geq 1. Rumor has it that these cuts are better.
2177 */
2178
2179 /* try to create GMI cut out of the aggregation row */
2180 SCIP_CALL( SCIPcalcMIR(scip, sol, POSTPROCESS, BOUNDSWITCH, USEVBDS, allowlocal, FIXINTEGRALRHS, NULL,
2181 NULL, minfrac, maxfrac, 1.0, aggrrow, cutcoefs, &cutrhs, cutinds, &cutnnz, &cutefficacy, &cutrank,
2182 &cutislocal, &success) );
2183
2184 if( success )
2185 {
2186 assert(allowlocal || !cutislocal); /*lint !e644*/
2187
2188 SCIP_CALL( constructCutRow(scip, sepa, sepadata, mainiternum, subgradientiternum, cutnnz, cutinds, cutcoefs,
2189 cutefficacy, cutrhs, cutislocal, cutrank, generatedcurrroundcuts, generatedcutefficacies,
2190 ngeneratedcurrroundcuts, ngeneratednewcuts, cutoff));
2191 }
2192 }
2193
2194 /* free temporary memory */
2195 SCIPfreeBufferArray(scip, &cutinds);
2196 SCIPfreeBufferArray(scip, &cutcoefs);
2197 SCIPfreeBufferArray(scip, &basisind);
2198 SCIPfreeBufferArray(scip, &inds);
2199 SCIPfreeBufferArray(scip, &binvrow);
2200 SCIPfreeBufferArray(scip, &basisfrac);
2201 SCIPfreeBufferArray(scip, &basisperm);
2202 SCIPaggrRowFree(scip, &aggrrow);
2203
2204 return SCIP_OKAY;
2205 }
2206
2207 /** update objective vector w.r.t. the fixed Lagrangian multipliers */
2208 static
2209 SCIP_RETCODE updateObjectiveVector(
2210 SCIP* scip, /**< SCIP data structure */
2211 SCIP_Real* dualvector, /**< Lagrangian multipliers vector */
2212 SCIP_ROW** cuts, /**< cuts generated so far in the current separation round */
2213 int ncuts, /**< number of cuts generated so far in the current separation round */
2214 SCIP_Real* origobjcoefs, /**< original objective function coefficients of the node linear relaxation */
2215 SCIP_Bool* objvecsdiffer /**< whether the updated objective function coefficients differ from the old ones */
2216 )
2217 {
2218 SCIP_Real* objvals;
2219 SCIP_Real* prod;
2220 SCIP_Real* oldobjvals;
2221 SCIP_Real* cutvals;
2222 SCIP_COL** cutcols;
2223 SCIP_COL** cols;
2224 SCIP_VAR* var;
2225 int cutnnz;
2226 int collppos;
2227 int ncols;
2228
2229 assert(objvecsdiffer != NULL);
2230 assert(ncuts > 0);
2231
2232 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2233
2234 /* allocate memory */
2235 SCIP_CALL( SCIPallocBufferArray(scip, &objvals, ncols) );
2236 SCIP_CALL( SCIPallocBufferArray(scip, &oldobjvals, ncols) );
2237 SCIP_CALL( SCIPallocCleanBufferArray(scip, &prod, ncols) );
2238 *objvecsdiffer = FALSE;
2239
2240 /* find the product of Lagrangian multipliers and cut coefficients */
2241 for( int i = 0; i < ncuts; i++ )
2242 {
2243 if( !SCIPisZero(scip, dualvector[i]) )
2244 {
2245 assert(!(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])) && SCIPisInfinity(scip, SCIProwGetRhs(cuts[i]))));
2246 assert(SCIPisInfinity(scip, -SCIProwGetLhs(cuts[i])));
2247 assert(!SCIPisInfinity(scip, SCIProwGetRhs(cuts[i])));
2248
2249 cutnnz = SCIProwGetNNonz(cuts[i]);
2250 assert(cutnnz <= ncols);
2251 cutvals = SCIProwGetVals(cuts[i]);
2252 cutcols = SCIProwGetCols(cuts[i]);
2253
2254 for( int j = 0; j < cutnnz; ++j )
2255 {
2256 collppos = SCIPcolGetLPPos(cutcols[j]);
2257 assert(collppos >= 0);
2258 assert(collppos <= ncols);
2259
2260 prod[collppos] += dualvector[i] * cutvals[j];
2261 }
2262 }
2263 }
2264
2265 /* change objective coefficients */
2266 for( int i = 0; i < ncols; i++ )
2267 {
2268 var = SCIPcolGetVar(cols[i]);
2269 oldobjvals[i] = SCIPgetVarObjDive(scip, var);
2270 objvals[i] = origobjcoefs[i] + prod[i];
2271
2272 SCIP_CALL( SCIPchgVarObjDive(scip, var, objvals[i]) );
2273
2274 /* identify if the updated objective vector is indeed different from the previous one */
2275 if( !(*objvecsdiffer) && !SCIPisEQ(scip, oldobjvals[i], objvals[i]) )
2276 *objvecsdiffer = TRUE;
2277 }
2278
2279 for( int i = 0; i < ncols; i++)
2280 {
2281 prod[i] = 0.0;
2282 }
2283
2284 /* free memory */
2285 SCIPfreeCleanBufferArray(scip, &prod);
2286 SCIPfreeBufferArray(scip, &oldobjvals);
2287 SCIPfreeBufferArray(scip, &objvals);
2288
2289 return SCIP_OKAY;
2290 }
2291
2292 /** add GMI cuts to the objective function of the Lagrangian dual problem by introducing new Lagrangian multipliers */
2293 static
2294 SCIP_RETCODE addGMICutsAsSoftConss(
2295 SCIP_Real* dualvector, /**< Lagrangian multipliers vector */
2296 int ngeneratedcuts, /**< number of cuts generated so far in the current separation round */
2297 int* naddedcuts, /**< number of cuts added so far in the current separation round to the Lagrangian dual problem
2298 upon penalization */
2299 int* nnewaddedsoftcuts /**< number of cuts added newly to the Lagrangian dual problem upon penalization */
2300 )
2301 {
2302 assert(*naddedcuts <= ngeneratedcuts);
2303
2304 /* set the initial penalty of the newly penalized cuts as zero */
2305 for( int i = *naddedcuts; i < ngeneratedcuts; i++ )
2306 dualvector[i] = 0.0;
2307
2308 *nnewaddedsoftcuts = ngeneratedcuts - *naddedcuts;
2309 *naddedcuts = ngeneratedcuts;
2310
2311 return SCIP_OKAY;
2312 }
2313
2314 /** solve the Lagrangian dual problem */
2315 static
2316 SCIP_RETCODE solveLagrangianDual(
2317 SCIP* scip, /**< SCIP data structure */
2318 SCIP_SEPA* sepa, /**< pointer to the separator */
2319 SCIP_SEPADATA* sepadata, /**< separator data structure */
2320 SCIP_SOL* sol, /**< data structure to store an LP solution upon solving a Lagrangian dual problem with
2321 fixed Lagrangian multipliers */
2322 SCIP_Real* solvals, /**< values of the LP solution obtained upon solving a Lagrangian dual problem with fixed
2323 Lagrangian multipliers */
2324 int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */
2325 SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */
2326 int depth, /**< depth of the current node in the tree */
2327 SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */
2328 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
2329 SCIP_Real* origobjcoefs, /**< original objective function coefficients of the node linear relaxation */
2330 SCIP_Real origobjoffset, /**< original objective function offset of the node linear relaxation */
2331 SCIP_Real* dualvector, /**< Lagrangian multipliers vector */
2332 int* nsoftcuts, /**< number of generated cuts that were penalized and added to the Lagrangian dual problem */
2333 SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */
2334 SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut
2335 generations */
2336 int* ngeneratedcutsperiter, /**< number of cuts generated per subgradient iteration in the current separation round */
2337 int* ngeneratedcurrroundcuts, /**< number of cuts generated so far in the current separation round */
2338 int* ncurrroundlpiters, /**< number of LP iterations taken for solving Lagrangian dual problems with fixed
2339 multipliers in the current separator round */
2340 SCIP_Bool* cutoff, /**< should the current node be cutoff? */
2341 SCIP_Real* bestlagrangianval, /**< best Lagrangian value found so far */
2342 SCIP_Real* bestdualvector, /**< Lagrangian multipliers corresponding to the best Lagrangian value found so far */
2343 int* bestdualvectorlen, /**< length of the Lagrangian multipliers corresponding to the best Lagrangian value
2344 found so far */
2345 int* nbestdualupdates, /**< number of best Lagrangian values found so far */
2346 int* totaliternum /**< total number of iterations of the relax-and-cut algorithm performed so far */
2347 )
2348 {
2349 SCIP_Real* subgradient;
2350 SCIP_Real muparam;
2351 SCIP_Real steplength;
2352 SCIP_Real objval;
2353 SCIP_Real lagrangianval;
2354 SCIP_Real* lagrangianvals;
2355 SCIP_Real avglagrangianval;
2356 SCIP_Real maxsoftcutviol;
2357 SCIP_Real maxnzsubgradientdualprod;
2358 SCIP_Real maxviolscore;
2359 SCIP_Real maxviolscoreold;
2360 SCIP_Real nviolscore;
2361 SCIP_Real nviolscoreold;
2362 SCIP_Real scoreweight;
2363 SCIP_Real ballradius;
2364 SCIP_Bool solfound;
2365 SCIP_Bool backtrack;
2366 SCIP_Bool terminate;
2367 SCIP_Bool subgradientzero;
2368 SCIP_Bool objvecsdiffer;
2369 SCIP_Bool dualvecsdiffer;
2370 SCIP_Bool solvelp;
2371 int ncurrroundlpiterslast;
2372 int nlpiters;
2373 int ngeneratednewcuts;
2374 int nnewaddedsoftcuts;
2375 int nsoftcutviols;
2376 int nnzsubgradientdualprod;
2377
2378 SCIP_CALL( SCIPallocBufferArray(scip, &lagrangianvals, sepadata->nmaxsubgradientiters) );
2379 SCIP_CALL( SCIPallocCleanBufferArray(scip, &subgradient, nmaxgeneratedperroundcuts) );
2380
2381 muparam = sepadata->muparaminit;
2382 steplength = 0.0;
2383 objval = 0.0;
2384 avglagrangianval = 0.0;
2385 maxsoftcutviol = 0.0;
2386 maxnzsubgradientdualprod = 0.0;
2387 maxviolscore = 0.0;
2388 nviolscore = 0.0;
2389 scoreweight = 1.0;
2390 ballradius = sepadata->radiusinit;
2391 ngeneratednewcuts = 0;
2392 nsoftcutviols = 0;
2393 nnzsubgradientdualprod = 0;
2394 terminate = FALSE;
2395 subgradientzero = FALSE;
2396 objvecsdiffer = FALSE;
2397 dualvecsdiffer = FALSE;
2398 solvelp = TRUE;
2399
2400 /* update objective vector based on input Lagrangian multipliers */
2401 if( *nsoftcuts > 0 )
2402 {
2403 SCIP_CALL( updateObjectiveVector(scip, dualvector, generatedcurrroundcuts, *nsoftcuts, origobjcoefs, &objvecsdiffer) );
2404 }
2405
2406 /* termination check */
2407 SCIP_CALL( checkLagrangianDualTermination(sepadata, -1, -1, FALSE, *ngeneratedcurrroundcuts,
2408 nmaxgeneratedperroundcuts, *ncurrroundlpiters, depth, &terminate) );
2409
2410 /* the subgradient algorithm loop */
2411 for( int i = 0; i < sepadata->nmaxsubgradientiters && !SCIPisStopped(scip) && !terminate; i++ )
2412 {
2413 solfound = FALSE;
2414 subgradientzero = FALSE;
2415 objvecsdiffer = FALSE;
2416 dualvecsdiffer = FALSE;
2417 nnewaddedsoftcuts = 0;
2418 scoreweight *= sepadata->radiusupdateweight;
2419
2420 ncurrroundlpiterslast = *ncurrroundlpiters;
2421 if( solvelp )
2422 {
2423 /* solve Lagrangian dual for fixed Lagrangian multipliers */
2424 SCIP_CALL( solveLagromoryLP(scip, sepadata, depth, origobjoffset, &solfound, sol, solvals, &objval,
2425 ncurrroundlpiters) );
2426 }
2427 nlpiters = *ncurrroundlpiters - ncurrroundlpiterslast;
2428
2429 /* if an optimal solution has been found, then generate cuts and do other operations */
2430 if( solfound )
2431 {
2432 /* generate GMI cuts if a new basis solution is found */
2433 if( (nlpiters >= 1) && (i % sepadata->cutgenfreq == 0) )
2434 {
2435 ngeneratednewcuts = 0;
2436 SCIP_CALL( generateGMICuts(scip, sepa, sepadata, mainiternum, i, sol, solvals,
2437 nmaxgeneratedperroundcuts, allowlocal, generatedcurrroundcuts, generatedcutefficacies,
2438 *ngeneratedcurrroundcuts, &ngeneratednewcuts, depth, cutoff));
2439 sepadata->ntotalcuts += ngeneratednewcuts;
2440 *ngeneratedcurrroundcuts += ngeneratednewcuts;
2441 ngeneratedcutsperiter[mainiternum * sepadata->nmaxsubgradientiters + i + 1] = ngeneratednewcuts;
2442 }
2443
2444 /* update subgradient, i.e., find the residuals of the penalized cuts, and determine various violations */
2445 updateSubgradient(scip, sol, generatedcurrroundcuts, *nsoftcuts, subgradient, dualvector, &subgradientzero,
2446 &nsoftcutviols, &maxsoftcutviol, &nnzsubgradientdualprod, &maxnzsubgradientdualprod);
2447
2448 /* calculate Lagrangian value for the fixed Lagrangian multipliers, and update best and avg values */
2449 updateLagrangianValue(scip, objval, dualvector, generatedcurrroundcuts, *nsoftcuts, &lagrangianval);
2450 if( SCIPisPositive(scip, lagrangianval - *bestlagrangianval) )
2451 {
2452 *bestlagrangianval = lagrangianval;
2453 for( int j = 0; j < *nsoftcuts; j++ )
2454 {
2455 bestdualvector[j] = dualvector[j];
2456 }
2457 *bestdualvectorlen = *nsoftcuts;
2458 (*nbestdualupdates)++;
2459 }
2460 lagrangianvals[i] = lagrangianval;
2461 if( i < sepadata->nmaxlagrangianvalsforavg )
2462 {
2463 avglagrangianval = (avglagrangianval * i + lagrangianval)/(i+1);
2464 }
2465 else
2466 {
2467 avglagrangianval = (avglagrangianval * sepadata->nmaxlagrangianvalsforavg -
2468 lagrangianvals[i - sepadata->nmaxlagrangianvalsforavg] +
2469 lagrangianval)/(sepadata->nmaxlagrangianvalsforavg);
2470 }
2471
2472 /* if the subgradient vector is non-zero, then update the mu parameter and the Lagrangian multipliers */
2473 if( !subgradientzero )
2474 {
2475 /* update mu param */
2476 SCIP_CALL( updateMuSteplengthParam(scip, sepadata, i, ubparam, lagrangianvals, *bestlagrangianval, avglagrangianval,
2477 &muparam, &backtrack) );
2478
2479 /* update step length */
2480 updateStepLength(scip, muparam, ubparam, lagrangianval, subgradient, *nsoftcuts, &steplength);
2481
2482 /* update scores to determine how to update the stabilization ball radius */
2483 maxviolscoreold = maxviolscore;
2484 nviolscoreold = nviolscore;
2485 maxviolscore = (1.0 - scoreweight) * maxsoftcutviol + scoreweight * maxnzsubgradientdualprod;
2486 nviolscore = (1.0 - scoreweight) * nsoftcutviols + scoreweight * nnzsubgradientdualprod;
2487
2488 /* update Lagrangian multipliers */
2489 SCIP_CALL( updateDualVector(scip, sepadata, dualvector, bestdualvector, *bestdualvectorlen,
2490 *nbestdualupdates, i, *totaliternum, steplength, subgradient, *nsoftcuts, backtrack, maxviolscore,
2491 maxviolscoreold, nviolscore, nviolscoreold, nlpiters, &dualvecsdiffer, &ballradius) );
2492
2493 /* update objective vector based on updated Lagrangian multipliers */
2494 if( dualvecsdiffer )
2495 {
2496 SCIP_CALL( updateObjectiveVector(scip, dualvector, generatedcurrroundcuts, *nsoftcuts, origobjcoefs, &objvecsdiffer) );
2497 }
2498 }
2499 /* if the subgradient vector if zero, then simply mark that the Lagrangian multipliers and the objective
2500 * function of the Lagrangian dual did not change */
2501 else
2502 {
2503 dualvecsdiffer = FALSE;
2504 objvecsdiffer = FALSE;
2505 }
2506
2507
2508 /* add generated GMI cuts to the objective function of the Lagrangian dual problem by introducing new
2509 * Lagrangian multipliers */
2510 if( (i % sepadata->cutaddfreq == 0) || (!dualvecsdiffer && !objvecsdiffer &&
2511 (*ngeneratedcurrroundcuts - *nsoftcuts > 0)) )
2512 {
2513 SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) );
2514 }
2515 }
2516 else
2517 {
2518 /* add any remaining generated GMI cuts to the objective function of the Lagrangian dual problem by introducing
2519 * new Lagrangian multipliers */
2520 if( (*ngeneratedcurrroundcuts - *nsoftcuts) > 0 )
2521 {
2522 SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) );
2523 }
2524
2525 solvelp = FALSE;
2526 }
2527
2528 /* termination check */
2529 SCIP_CALL( checkLagrangianDualTermination(sepadata, nnewaddedsoftcuts, *ngeneratedcurrroundcuts - *nsoftcuts,
2530 objvecsdiffer, *ngeneratedcurrroundcuts, nmaxgeneratedperroundcuts, *ncurrroundlpiters, depth,
2531 &terminate) );
2532
2533 (*totaliternum)++;
2534 }
2535
2536 /* add any remaining generated GMI cuts to the objective function of the Lagrangian dual problem by introducing new
2537 * Lagrangian multipliers */
2538 if( (*ngeneratedcurrroundcuts - *nsoftcuts) > 0 )
2539 {
2540 SCIP_CALL( addGMICutsAsSoftConss(dualvector, *ngeneratedcurrroundcuts, nsoftcuts, &nnewaddedsoftcuts) );
2541 }
2542
2543 /* free memory */
2544 for( int i = 0; i < nmaxgeneratedperroundcuts; i++)
2545 {
2546 subgradient[i] = 0.0;
2547 }
2548 SCIPfreeCleanBufferArray(scip, &subgradient);
2549 SCIPfreeBufferArray(scip, &lagrangianvals);
2550
2551 return SCIP_OKAY;
2552 }
2553
2554 /** generates initial cut pool before solving the Lagrangian dual */
2555 static
2556 SCIP_RETCODE generateInitCutPool(
2557 SCIP* scip, /**< SCIP data structure */
2558 SCIP_SEPA* sepa, /**< separator */
2559 SCIP_SEPADATA* sepadata, /**< separator data structure */
2560 int mainiternum, /**< iteration number of the outer loop of the relax-and-cut algorithm */
2561 SCIP_SOL* sol, /**< LP solution to be used for cut generation */
2562 SCIP_Real* solvals, /**< values of the LP solution to be used for cut generation */
2563 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
2564 SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */
2565 SCIP_ROW** generatedcurrroundcuts, /**< cuts generated in the current separation round */
2566 SCIP_Real* generatedcutefficacies, /**< array of generated cut efficacies w.r.t. the respective LP bases used for cut
2567 generations */
2568 int* ngeneratedcutsperiter, /**< number of cuts generated per subgradient iteration in the current separation round */
2569 int* ngeneratedcurrroundcuts, /**< number of cuts generated so far in the current separation round */
2570 int depth, /**< depth of the current node in the tree */
2571 SCIP_Bool* cutoff /**< should the current node be cutoff? */
2572 )
2573 {
2574 int ngeneratednewcuts;
2575
2576 ngeneratednewcuts = 0;
2577
2578 /* generate initial set of cuts */
2579 SCIP_CALL( generateGMICuts(scip, sepa, sepadata, mainiternum, -1, sol, solvals, nmaxgeneratedperroundcuts,
2580 allowlocal, generatedcurrroundcuts, generatedcutefficacies, *ngeneratedcurrroundcuts, &ngeneratednewcuts,
2581 depth, cutoff) );
2582
2583 /* update certain statistics */
2584 sepadata->ntotalcuts += ngeneratednewcuts;
2585 *ngeneratedcurrroundcuts += ngeneratednewcuts;
2586 ngeneratedcutsperiter[sepadata->nmaxsubgradientiters * mainiternum] = ngeneratednewcuts;
2587
2588 return SCIP_OKAY;
2589 }
2590
2591 /** add cuts to SCIP */
2592 static
2593 SCIP_RETCODE addCuts(
2594 SCIP* scip, /**< SCIP data structure */
2595 SCIP_SEPADATA* sepadata, /**< separator data structure */
2596 SCIP_ROW** cuts, /**< cuts generated so far in the current separation round */
2597 int ncuts, /**< number of cuts generated so far in the current separation round */
2598 SCIP_Longint maxdnom, /**< maximum denominator in the rational representation of cuts */
2599 SCIP_Real maxscale, /**< maximal scale factor to scale the cuts to integral values */
2600 int* naddedcuts, /**< number of cuts added to either global cutpool or sepastore */
2601 SCIP_Bool* cutoff /**< should the current node be cutoff? */
2602 )
2603 {
2604 SCIP_ROW* cut;
2605 SCIP_Bool madeintegral;
2606
2607 assert(scip != NULL);
2608 assert(sepadata != NULL);
2609 assert(naddedcuts != NULL);
2610 assert(cutoff != NULL);
2611
2612 *cutoff = FALSE;
2613 madeintegral = FALSE;
2614
2615 for( int i = 0; i < ncuts && !*cutoff; i++ )
2616 {
2617 cut = cuts[i];
2618
2619 if( SCIProwGetNNonz(cut) == 1 )
2620 {
2621 /* Add the bound change as cut to avoid that the LP gets modified. This would mean that the LP is not flushed
2622 * and the method SCIPgetLPBInvRow() fails; SCIP internally will apply this bound change automatically. */
2623 SCIP_CALL( SCIPaddRow(scip, cut, TRUE, cutoff) );
2624 ++(*naddedcuts);
2625 }
2626 else
2627 {
2628 if( sepadata->makeintegral && SCIPgetRowNumIntCols(scip, cut) == SCIProwGetNNonz(cut) )
2629 {
2630 /* try to scale the cut to integral values */
2631 SCIP_CALL( SCIPmakeRowIntegral(scip, cut, -SCIPepsilon(scip), SCIPsumepsilon(scip),
2632 maxdnom, maxscale, MAKECONTINTEGRAL, &madeintegral) );
2633
2634 /* if RHS = plus infinity (due to scaling), the cut is useless, so we are not adding it */
2635 if( madeintegral && SCIPisInfinity(scip, SCIProwGetRhs(cut)) )
2636 return SCIP_OKAY;
2637 }
2638
2639 if( SCIPisCutNew(scip, cut) )
2640 {
2641 /* add global cuts which are not implicit bound changes to the cut pool */
2642 if( !SCIProwIsLocal(cut) )
2643 {
2644 if( sepadata->delayedcuts )
2645 {
2646 SCIP_CALL( SCIPaddDelayedPoolCut(scip, cut) );
2647 }
2648 else
2649 {
2650 SCIP_CALL( SCIPaddPoolCut(scip, cut) );
2651 }
2652 }
2653 else
2654 {
2655 /* local cuts we add to the sepastore */
2656 SCIP_CALL( SCIPaddRow(scip, cut, sepadata->forcecuts, cutoff) );
2657 }
2658 ++(*naddedcuts);
2659 }
2660 }
2661 }
2662
2663 return SCIP_OKAY;
2664 }
2665
2666 /** check different termination criteria */
2667 /* @todo nlpssolved criterion? */
2668 static
2669 SCIP_RETCODE checkMainLoopTermination(
2670 SCIP_SEPADATA* sepadata, /**< separator data structure */
2671 SCIP_Bool cutoff, /**< should the current node be cutoff? */
2672 SCIP_Bool dualvecsdiffer, /**< whether the updated Lagrangian multipliers differ from the old one */
2673 int ngeneratedcurrroundcuts, /**< number of cuts generated in the current separation round */
2674 int nsoftcuts, /**< number of generated cuts that were penalized and added to the Lagrangian dual problem */
2675 int nmaxgeneratedperroundcuts,/**< maximal number of cuts allowed to generate per separation round */
2676 int ncurrroundlpiters, /**< number of LP iterations taken for solving Lagrangian dual problems with fixed
2677 multipliers in the current separator round */
2678 int depth, /**< depth of the current node in the tree */
2679 SCIP_Bool* terminate /**< whether to terminate the relax-and-cut algorithm */
2680 )
2681 {
2682 *terminate = FALSE;
2683
2684 /* check if the node has been identified to be cutoff */
2685 if( cutoff )
2686 *terminate = TRUE;
2687
2688 /* check if the Lagrangian multipliers do not differ from the previous iteration and no new cuts exist for penalizing */
2689 if( !dualvecsdiffer && (ngeneratedcurrroundcuts == nsoftcuts) )
2690 *terminate = TRUE;
2691
2692 /* check if allowed number of cuts in this separation round have already been generated */
2693 if( ngeneratedcurrroundcuts >= nmaxgeneratedperroundcuts )
2694 *terminate = TRUE;
2695
2696 /* check if allowed number of cuts in the tree have already been generated */
2697 if( sepadata->ntotalcuts >= sepadata->nmaxtotalcuts )
2698 *terminate = TRUE;
2699
2700 /* check if allowed number of simplex iterations in this separation round have already been used up */
2701 if( (sepadata->nmaxperroundlpiters >= 0) && (ncurrroundlpiters >= sepadata->nmaxperroundlpiters) )
2702 *terminate = TRUE;
2703
2704 /* check if allowed number of simplex iterations in the root node have already been used up */
2705 if( (depth == 0) && (sepadata->nmaxrootlpiters >= 0) && (sepadata->nrootlpiters >= sepadata->nmaxrootlpiters) )
2706 *terminate = TRUE;
2707
2708 /* check if allowed number of simplex iterations in the tree have already been used up */
2709 if( (sepadata->nmaxtotallpiters >= 0) && (sepadata->ntotallpiters >= sepadata->nmaxtotallpiters) )
2710 *terminate = TRUE;
2711
2712 return SCIP_OKAY;
2713 }
2714
2715 /** Searches and tries to add Lagromory cuts */
2716 static
2717 SCIP_RETCODE separateCuts(
2718 SCIP* scip, /**< SCIP data structure */
2719 SCIP_SEPA* sepa, /**< separator */
2720 SCIP_SEPADATA* sepadata, /**< separator data structure */
2721 SCIP_Real ubparam, /**< estimate of the optimal Lagrangian dual value */
2722 int depth, /**< depth of the current node in the tree */
2723 SCIP_Bool allowlocal, /**< should locally valid cuts be generated? */
2724 SCIP_RESULT* result /**< final result of the separation round */
2725 )
2726 {
2727 SCIP_ROW** generatedcurrroundcuts;
2728 SCIP_ROW** aggregatedcurrroundcuts;
2729 SCIP_Real* generatedcutefficacies;
2730 SCIP_Bool solfound;
2731 SCIP_SOL* softcutslpsol;
2732 SCIP_Real* softcutslpsolvals;
2733 SCIP_SOL* hardcutslpsol;
2734 SCIP_Real* hardcutslpsolvals;
2735 SCIP_Real* dualsol;
2736 SCIP_Real* dualvector;
2737 SCIP_Real* bestdualvector;
2738 SCIP_Real bestlagrangianval;
2739 SCIP_Real* origobjcoefs;
2740 SCIP_Real origobjoffset;
2741 SCIP_Real objval;
2742 SCIP_Real maxscale;
2743 SCIP_Longint maxdnom;
2744 SCIP_Bool cutoff;
2745 SCIP_Bool cutoff2;
2746 SCIP_Bool dualvecsdiffer;
2747 SCIP_Bool terminate;
2748 SCIP_COL** cols;
2749
2750 int* ngeneratedcutsperiter;
2751 int bestdualvectorlen;
2752 int nbestdualupdates;
2753 int totaliternum;
2754 int* cutindsperm;
2755 int nprocessedcuts;
2756 int ncols;
2757 int nrows;
2758 int ngeneratedcurrroundcuts;
2759 int nselectedcurrroundcuts;
2760 int nselectedcuts;
2761 int naddedcurrroundcuts;
2762 int naggregatedcurrroundcuts;
2763 int nmaxgeneratedperroundcuts;
2764 int ncurrroundlpiters;
2765 int nsoftcuts;
2766 int nsoftcutsold;
2767 int maxdepth;
2768
2769 assert(*result == SCIP_DIDNOTRUN);
2770 assert(sepadata != NULL);
2771 sepadata->ncalls = SCIPsepaGetNCalls(sepa);
2772 objval = SCIPinfinity(scip);
2773 bestlagrangianval = -SCIPinfinity(scip);
2774 bestdualvectorlen = 0;
2775 nbestdualupdates = 0;
2776 totaliternum = 0;
2777 ngeneratedcurrroundcuts = 0;
2778 naddedcurrroundcuts = 0;
2779 naggregatedcurrroundcuts = 0;
2780 ncurrroundlpiters = 0;
2781 nsoftcuts = 0;
2782 solfound = FALSE;
2783 cutoff = FALSE;
2784 cutoff2 = FALSE;
2785 dualvecsdiffer = FALSE;
2786 terminate = FALSE;
2787
2788 SCIPdebugMsg(scip, "Separating cuts...\n");
2789
2790 /* initialize the LP that will be used to solve the Lagrangian dual with fixed Lagrangian multipliers */
2791 SCIP_CALL( createLPWithSoftCuts(scip, sepadata) );
2792 /* create the LP that represents the node relaxation including all the generated cuts in this separator */
2793 SCIP_CALL( createLPWithHardCuts(scip, sepadata, NULL, 0) );
2794
2795 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2796 nrows = SCIPgetNLPRows(scip);
2797
2798 /* get the maximal number of cuts allowed in a separation round */
2799 nmaxgeneratedperroundcuts = ((depth == 0) ? sepadata->nmaxperroundcutsroot : sepadata->nmaxperroundcuts);
2800
2801 /* allocate memory */
2802 SCIP_CALL( SCIPallocBufferArray(scip, &generatedcurrroundcuts, nmaxgeneratedperroundcuts) );
2803 SCIP_CALL( SCIPallocBufferArray(scip, &aggregatedcurrroundcuts, nmaxgeneratedperroundcuts) );
2804 SCIP_CALL( SCIPallocBufferArray(scip, &generatedcutefficacies, nmaxgeneratedperroundcuts) );
2805 SCIP_CALL( SCIPallocBufferArray(scip, &cutindsperm, nmaxgeneratedperroundcuts) );
2806 SCIP_CALL( SCIPallocCleanBufferArray(scip, &ngeneratedcutsperiter, sepadata->nmaxmainiters *
2807 (sepadata->nmaxsubgradientiters + 1)) );
2808 SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualsol, nrows + nmaxgeneratedperroundcuts) );
2809 SCIP_CALL( SCIPallocCleanBufferArray(scip, &dualvector, nmaxgeneratedperroundcuts) );
2810 SCIP_CALL( SCIPallocCleanBufferArray(scip, &bestdualvector, nmaxgeneratedperroundcuts) );
2811 SCIP_CALL( SCIPallocBufferArray(scip, &softcutslpsolvals, ncols) );
2812 SCIP_CALL( SCIPcreateSol(scip, &softcutslpsol, NULL) );
2813 SCIP_CALL( SCIPallocBufferArray(scip, &hardcutslpsolvals, ncols) );
2814 SCIP_CALL( SCIPcreateSol(scip, &hardcutslpsol, NULL) );
2815 SCIP_CALL( SCIPallocBufferArray(scip, &origobjcoefs, ncols) );
2816
2817 /* store current objective function */
2818 for( int i = 0; i < ncols; i++ )
2819 {
2820 origobjcoefs[i] = SCIPcolGetObj(cols[i]);
2821 }
2822 origobjoffset = SCIPgetTransObjoffset(scip);
2823
2824 /* solve node LP relaxation to have an initial simplex tableau */
2825 SCIP_CALL( solveLagromoryLP(scip, sepadata, depth, origobjoffset, &solfound, softcutslpsol, softcutslpsolvals, &objval,
2826 &ncurrroundlpiters));
2827
2828 /* generate initial cut pool */
2829 SCIP_CALL( generateInitCutPool(scip, sepa, sepadata, 0, softcutslpsol, softcutslpsolvals, nmaxgeneratedperroundcuts, allowlocal,
2830 generatedcurrroundcuts, generatedcutefficacies, ngeneratedcutsperiter, &ngeneratedcurrroundcuts, depth,
2831 &cutoff) );
2832
2833 /* termination check */
2834 SCIP_CALL( checkMainLoopTermination(sepadata, cutoff, TRUE, ngeneratedcurrroundcuts, nsoftcuts,
2835 nmaxgeneratedperroundcuts, ncurrroundlpiters, depth, &terminate) );
2836
2837 /* compute cuts for each integer col with fractional val */
2838 for( int i = 0; i < sepadata->nmaxmainiters && !SCIPisStopped(scip) && !terminate; ++i )
2839 {
2840 nsoftcutsold = nsoftcuts;
2841 nsoftcuts = ngeneratedcurrroundcuts;
2842 dualvecsdiffer = FALSE;
2843
2844 /* solve Lagrangain dual */
2845 SCIP_CALL( solveLagrangianDual(scip, sepa, sepadata, softcutslpsol, softcutslpsolvals, i, ubparam, depth, allowlocal,
2846 nmaxgeneratedperroundcuts, origobjcoefs, origobjoffset, dualvector, &nsoftcuts, generatedcurrroundcuts,
2847 generatedcutefficacies, ngeneratedcutsperiter, &ngeneratedcurrroundcuts, &ncurrroundlpiters, &cutoff,
2848 &bestlagrangianval, bestdualvector, &bestdualvectorlen, &nbestdualupdates, &totaliternum) );
2849
2850 /* @todo filter cuts before adding them to the new LP that was created based on the node relaxation? */
2851
2852 /* update the LP representing the node relaxation by adding newly generated cuts */
2853 if( !cutoff && (ngeneratedcurrroundcuts - nsoftcutsold > 0) )
2854 {
2855 SCIP_CALL( createLPWithHardCuts(scip, sepadata, generatedcurrroundcuts, ngeneratedcurrroundcuts) );
2856
2857 /* solve the LP and get relevant information */
2858 SCIP_CALL( solveLPWithHardCuts(scip, sepadata, &solfound, hardcutslpsol, hardcutslpsolvals) );
2859
2860 /* if primal solution is found, then pass it to trysol heuristic */
2861 /* @note if trysol heuristic is was not present, then the solution found above, which can potentially be a MIP
2862 * primal feasible solution, will go to waste.
2863 */
2864 if( solfound && sepadata->heurtrysol != NULL )
2865 {
2866 SCIP_CALL( SCIPheurPassSolTrySol(scip, sepadata->heurtrysol, hardcutslpsol) );
2867 }
2868
2869 /* if dual feasible, then fetch dual solution and reset Lagrangian multipliers based on it. otherwise, retain the
2870 * Lagrangian multipliers and simply initialize the new multipliers to zeroes. */
2871 if( SCIPlpiIsDualFeasible(sepadata->lpiwithhardcuts) )
2872 {
2873 SCIP_CALL( SCIPlpiGetSol(sepadata->lpiwithhardcuts, NULL, NULL, dualsol, NULL, NULL) );
2874 SCIP_CALL( updateDualVector(scip, sepadata, dualvector, &(dualsol[nrows]),
2875 ngeneratedcurrroundcuts, 0, -1, -1, 0.0, NULL, ngeneratedcurrroundcuts, TRUE, 0.0, 0.0, 0.0, 0.0, -1,
2876 &dualvecsdiffer, NULL) );
2877 }
2878 else
2879 {
2880 SCIP_CALL( updateDualVector(scip, sepadata, dualvector, dualvector, nsoftcuts, 0, -1, -1, 0.0, NULL,
2881 ngeneratedcurrroundcuts, TRUE, 0.0, 0.0, 0.0, 0.0, -1, &dualvecsdiffer, NULL) );
2882 }
2883 }
2884
2885 /* termination check */
2886 SCIP_CALL( checkMainLoopTermination(sepadata, cutoff, dualvecsdiffer, ngeneratedcurrroundcuts, nsoftcuts,
2887 nmaxgeneratedperroundcuts, ncurrroundlpiters, depth, &terminate) );
2888 }
2889
2890 /* set the maximal denominator in rational representation of gomory cut and the maximal scale factor to
2891 * scale resulting cut to integral values to avoid numerical instabilities
2892 */
2893 /**@todo find better but still stable gomory cut settings: look at dcmulti, gesa3, khb0525, misc06, p2756 */
2894 /* @note: above todo was copied from sepa_gomory.c. So, if gomory code is changed, same changes can be done here. */
2895 maxdepth = SCIPgetMaxDepth(scip);
2896 if( depth == 0 )
2897 {
2898 maxdnom = 1000;
2899 maxscale = 1000.0;
2900 }
2901 else if( depth <= maxdepth/4 )
2902 {
2903 maxdnom = 1000;
2904 maxscale = 1000.0;
2905 }
2906 else if( depth <= maxdepth/2 )
2907 {
2908 maxdnom = 100;
2909 maxscale = 100.0;
2910 }
2911 else
2912 {
2913 maxdnom = 10;
2914 maxscale = 10.0;
2915 }
2916
2917 /* filter cuts by calling cut selection algorithms and add cuts accordingly */
2918 /* @todo an idea is to filter cuts after every main iter */
2919 /* @todo we can remove !cutoff criterion by adding forcedcuts */
2920 if( !cutoff && (ngeneratedcurrroundcuts > 0) )
2921 {
2922 if( SCIPisGE(scip, sepadata->cutsfilterfactor, 1.0) )
2923 {
2924 nselectedcurrroundcuts = ngeneratedcurrroundcuts;
2925 SCIP_CALL( addCuts(scip, sepadata, generatedcurrroundcuts, nselectedcurrroundcuts, maxdnom, maxscale,
2926 &naddedcurrroundcuts, &cutoff2) );
2927 cutoff = cutoff2;
2928 }
2929 else if( SCIPisPositive(scip, sepadata->cutsfilterfactor) )
2930 {
2931 nprocessedcuts = 0;
2932 for( int i = 0; i < sepadata->nmaxmainiters * (sepadata->nmaxsubgradientiters + 1); i++ )
2933 {
2934 if( ngeneratedcutsperiter[i] != 0 )
2935 {
2936 for( int j = 0; j < ngeneratedcutsperiter[i]; j++ )
2937 cutindsperm[j] = j + nprocessedcuts;
2938
2939 /* sort cut efficacies by fractionality */
2940 SCIPsortDownRealInt(&generatedcutefficacies[nprocessedcuts], cutindsperm, ngeneratedcutsperiter[i]);
2941
2942 nselectedcuts = (int)SCIPceil(scip, sepadata->cutsfilterfactor * ngeneratedcutsperiter[i]);
2943
2944 SCIP_CALL( addCuts(scip, sepadata, &generatedcurrroundcuts[nprocessedcuts], nselectedcuts, maxdnom,
2945 maxscale, &naddedcurrroundcuts, &cutoff2) );
2946 cutoff = cutoff2;
2947
2948 nprocessedcuts += ngeneratedcutsperiter[i];
2949 }
2950 }
2951 }
2952 }
2953 else if( ngeneratedcurrroundcuts > 0 )
2954 {
2955 nselectedcurrroundcuts = ngeneratedcurrroundcuts;
2956 SCIP_CALL( addCuts(scip, sepadata, generatedcurrroundcuts, nselectedcurrroundcuts, maxdnom, maxscale,
2957 &naddedcurrroundcuts, &cutoff2) );
2958 }
2959
2960 /* add an aggregated cut based on best Lagrangian multipliers */
2961 if( sepadata->aggregatecuts && (ngeneratedcurrroundcuts > 0) && (bestdualvectorlen > 0) )
2962 {
2963 assert(bestdualvectorlen <= ngeneratedcurrroundcuts);
2964 SCIP_CALL( aggregateGeneratedCuts(scip, sepa, sepadata, generatedcurrroundcuts, bestdualvector, bestdualvectorlen,
2965 aggregatedcurrroundcuts, &naggregatedcurrroundcuts, &cutoff2) );
2966 cutoff = (!cutoff ? cutoff2 : cutoff);
2967 if( naggregatedcurrroundcuts > 0 )
2968 {
2969 SCIP_CALL( addCuts(scip, sepadata, aggregatedcurrroundcuts, naggregatedcurrroundcuts, maxdnom, maxscale,
2970 &naddedcurrroundcuts, &cutoff2) );
2971 cutoff = (!cutoff ? cutoff2 : cutoff);
2972 }
2973 }
2974
2975 if( cutoff )
2976 {
2977 *result = SCIP_CUTOFF;
2978 }
2979 else if( naddedcurrroundcuts > 0 )
2980 {
2981 *result = SCIP_SEPARATED;
2982 }
2983 else
2984 {
2985 *result = SCIP_DIDNOTFIND;
2986 }
2987
2988 /* delete the LP representing the Lagrangian dual problem with fixed Lagrangian multipliers */
2989 SCIP_CALL( deleteLPWithSoftCuts(scip, sepadata) );
2990
2991 /* release the rows in aggregatedcurrroundcuts */
2992 for( int i = 0; i < naggregatedcurrroundcuts; ++i )
2993 {
2994 SCIP_CALL( SCIPreleaseRow(scip, &(aggregatedcurrroundcuts[i])) );
2995 }
2996 /* release the rows in generatedcurrroundcuts */
2997 for( int i = 0; i < ngeneratedcurrroundcuts; ++i )
2998 {
2999 SCIP_CALL( SCIPreleaseRow(scip, &(generatedcurrroundcuts[i])) );
3000 }
3001
3002 for( int i = 0; i < sepadata->nmaxmainiters * (sepadata->nmaxsubgradientiters + 1); i++ )
3003 {
3004 ngeneratedcutsperiter[i] = 0;
3005 }
3006 for( int i = 0; i < nrows; i++ )
3007 {
3008 dualsol[i] = 0.0;
3009 }
3010 for( int i = 0; i < nmaxgeneratedperroundcuts; i++ )
3011 {
3012 dualsol[nrows + i] = 0.0;
3013 dualvector[i] = 0.0;
3014 bestdualvector[i] = 0.0;
3015 }
3016 /* free memory */
3017 SCIPfreeBufferArray(scip, &origobjcoefs);
3018 SCIPfreeBufferArray(scip, &hardcutslpsolvals);
3019 SCIP_CALL( SCIPfreeSol(scip, &hardcutslpsol) );
3020 SCIPfreeBufferArray(scip, &softcutslpsolvals);
3021 SCIP_CALL( SCIPfreeSol(scip, &softcutslpsol) );
3022 SCIPfreeCleanBufferArray(scip, &ngeneratedcutsperiter);
3023 SCIPfreeBufferArray(scip, &cutindsperm);
3024 SCIPfreeBufferArray(scip, &generatedcutefficacies);
3025 SCIPfreeBufferArray(scip, &aggregatedcurrroundcuts);
3026 SCIPfreeBufferArray(scip, &generatedcurrroundcuts);
3027 SCIPfreeCleanBufferArray(scip, &dualvector);
3028 SCIPfreeCleanBufferArray(scip, &bestdualvector);
3029 SCIPfreeCleanBufferArray(scip, &dualsol);
3030
3031 return SCIP_OKAY;
3032 }
3033
3034 /** creates separator data */
3035 static
3036 SCIP_RETCODE sepadataCreate(
3037 SCIP* scip, /**< SCIP data structure */
3038 SCIP_SEPADATA** sepadata /**< separator data structure */
3039 )
3040 {
3041 assert(scip != NULL);
3042 assert(sepadata != NULL);
3043
3044 SCIP_CALL( SCIPallocBlockMemory(scip, sepadata) );
3045 BMSclearMemory(*sepadata);
3046
3047 return SCIP_OKAY;
3048 }
3049
3050
3051 /*
3052 * Callback methods of separator
3053 */
3054
3055 /** copy method for separator plugins (called when SCIP copies plugins) */
3056 static
3057 SCIP_DECL_SEPACOPY(sepaCopyLagromory)
3058 { /*lint --e{715}*/
3059 assert(scip != NULL);
3060 assert(sepa != NULL);
3061 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3062
3063 /* call inclusion method of constraint handler */
3064 SCIP_CALL( SCIPincludeSepaLagromory(scip) );
3065
3066 return SCIP_OKAY;
3067 }
3068
3069
3070 /** destructor of separator to free user data (called when SCIP is exiting) */
3071 static
3072 SCIP_DECL_SEPAFREE(sepaFreeLagromory)
3073 {
3074 SCIP_SEPADATA* sepadata;
3075
3076 sepadata = SCIPsepaGetData(sepa);
3077 assert(sepadata != NULL);
3078
3079 SCIP_CALL( sepadataFree(scip, &sepadata) );
3080 SCIPsepaSetData(sepa, NULL);
3081
3082 return SCIP_OKAY;
3083 }
3084
3085
3086 /** initialization method of separator (called after problem was transformed) */
3087 static
3088 SCIP_DECL_SEPAINIT(sepaInitLagromory)
3089 {
3090 SCIP_SEPADATA* sepadata;
3091
3092 sepadata = SCIPsepaGetData(sepa);
3093 assert(scip != NULL);
3094 assert(sepadata != NULL);
3095
3096 /* create and initialize random number generator */
3097 SCIP_CALL( SCIPcreateRandom(scip, &(sepadata->randnumgen), RANDSEED, TRUE) );
3098
3099 /* find trysol heuristic */
3100 if ( sepadata->heurtrysol == NULL )
3101 {
3102 sepadata->heurtrysol = SCIPfindHeur(scip, "trysol");
3103 }
3104
3105 return SCIP_OKAY;
3106 }
3107
3108 /** deinitialization method of separator (called before transformed problem is freed) */
3109 static
3110 SCIP_DECL_SEPAEXIT(sepaExitLagromory)
3111 { /*lint --e{715}*/
3112 SCIP_SEPADATA* sepadata;
3113
3114 sepadata = SCIPsepaGetData(sepa);
3115 assert(sepadata != NULL);
3116
3117 SCIPfreeRandom(scip, &(sepadata->randnumgen));
3118
3119 return SCIP_OKAY;
3120 }
3121
3122 /** LP solution separation method of separator */
3123 static
3124 SCIP_DECL_SEPAEXECLP(sepaExeclpLagromory)
3125 {
3126 /*lint --e{715}*/
3127
3128 SCIP_SEPADATA* sepadata;
3129 SCIP_SOL* bestsol;
3130 SCIP_COL** cols;
3131 SCIP_COL* col;
3132 SCIP_VAR* var;
3133 SCIP_Real ubparam;
3134 SCIP_Real dualdegeneracyrate;
3135 SCIP_Real varconsratio;
3136 SCIP_Real threshold1;
3137 SCIP_Real threshold2;
3138 int nrows;
3139 int ncols;
3140 int ncalls;
3141 int runnum;
3142
3143 assert(sepa != NULL);
3144 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3145 sepadata = SCIPsepaGetData(sepa);
3146 assert(sepadata != NULL);
3147
3148 assert(result != NULL);
3149 *result = SCIP_DIDNOTRUN;
3150
3151 assert(scip != NULL);
3152
3153 ncalls = SCIPsepaGetNCallsAtNode(sepa);
3154 runnum = SCIPgetNRuns(scip);
3155 dualdegeneracyrate = 0.0;
3156 varconsratio = 0.0;
3157
3158 /* only generate Lagromory cuts if we are not close to terminating */
3159 if( SCIPisStopped(scip) )
3160 return SCIP_OKAY;
3161
3162 /* only call the separator starting sepadata->minrestart runs */
3163 if( (sepadata->minrestart >= 1) && (runnum < sepadata->minrestart + 1) )
3164 return SCIP_OKAY;
3165
3166 /* only call the separator a given number of times at each node */
3167 if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
3168 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
3169 return SCIP_OKAY;
3170
3171 /* only generate cuts if an optimal LP solution is at hand */
3172 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
3173 return SCIP_OKAY;
3174
3175 /* only generate cuts if the LP solution is basic */
3176 if( ! SCIPisLPSolBasic(scip) )
3177 return SCIP_OKAY;
3178
3179 /* get LP data */
3180 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
3181 assert(cols != NULL);
3182
3183 nrows = SCIPgetNLPRows(scip);
3184
3185 /* return if LP has no columns or no rows */
3186 if( ncols == 0 || nrows == 0 )
3187 return SCIP_OKAY;
3188
3189 /* return if dual degeneracy metrics are below threshold values */
3190 threshold1 = sepadata->dualdegeneracyratethreshold;
3191 threshold2 = sepadata->varconsratiothreshold;
3192 SCIP_CALL( SCIPgetLPDualDegeneracy(scip, &dualdegeneracyrate, &varconsratio) );
3193 if( (dualdegeneracyrate < threshold1) && (varconsratio < threshold2) )
3194 return SCIP_OKAY;
3195
3196 bestsol = SCIPgetBestSol(scip);
3197 ubparam = 0.0;
3198
3199 /* if a MIP primal solution exists, use it to estimate the optimal value of the Lagrangian dual problem */
3200 if( bestsol != NULL )
3201 {
3202 for( int i = 0; i < ncols; ++i )
3203 {
3204 col = cols[i];
3205 assert(col != NULL);
3206
3207 var = SCIPcolGetVar(col);
3208
3209 ubparam += SCIPgetSolVal(scip, bestsol, var) * SCIPcolGetObj(col);
3210 }
3211
3212 ubparam += SCIPgetTransObjoffset(scip);
3213 }
3214 /* if a MIP primal solution does not exist, then use the node relaxation's LP solutin to estimate the optimal value
3215 * of the Lagrangian dual problem
3216 */
3217 else
3218 {
3219 for( int i = 0; i < ncols; ++i )
3220 {
3221 col = cols[i];
3222 assert(col != NULL);
3223
3224 ubparam += SCIPcolGetPrimsol(col) * SCIPcolGetObj(col);
3225 }
3226
3227 ubparam += SCIPgetTransObjoffset(scip);
3228 ubparam *= SCIPisPositive(scip, ubparam) ? sepadata->ubparamposfactor : sepadata->ubparamnegfactor;
3229 }
3230
3231 /* the main separation function of the separator */
3232 SCIP_CALL( separateCuts(scip, sepa, sepadata, ubparam, depth, (allowlocal && sepadata->allowlocal), result) );
3233
3234 return SCIP_OKAY;
3235 }
3236
3237
3238 /*
3239 * separator specific interface methods
3240 */
3241
3242 /** creates the Lagromory separator and includes it in SCIP */
3243 SCIP_RETCODE SCIPincludeSepaLagromory(
3244 SCIP* scip /**< SCIP data structure */
3245 )
3246 {
3247 SCIP_SEPADATA* sepadata;
3248 SCIP_SEPA* sepa;
3249
3250 /* create separator data */
3251 SCIP_CALL( sepadataCreate(scip, &sepadata) );
3252
3253 sepadata->heurtrysol = NULL;
3254
3255 /* include separator */
3256 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
3257 SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpLagromory, NULL, sepadata) );
3258
3259 assert(sepa != NULL);
3260
3261 /* set non fundamental callbacks via setter functions */
3262 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyLagromory) );
3263 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeLagromory) );
3264 SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitLagromory) );
3265 SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitLagromory) );
3266
3267 /* add separator parameters */
3268 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/away", "minimal integrality violation of a basis "
3269 "variable to try separation", &sepadata->away, FALSE, DEFAULT_AWAY, 0.0, 1.0, NULL, NULL) );
3270
3271 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/rootlpiterlimitfactor", "factor w.r.t. root node LP "
3272 "iterations for maximal separating LP iterations in the root node (negative for no limit)",
3273 &sepadata->rootlpiterlimitfactor, TRUE, DEFAULT_ROOTLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) );
3274
3275 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/totallpiterlimitfactor", "factor w.r.t. root node LP "
3276 "iterations for maximal separating LP iterations in the tree (negative for no limit)",
3277 &sepadata->totallpiterlimitfactor, TRUE, DEFAULT_TOTALLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) );
3278
3279 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundlpiterlimitfactor", "factor w.r.t. root node LP "
3280 "iterations for maximal separating LP iterations per separation round (negative for no limit)",
3281 &sepadata->perroundlpiterlimitfactor, TRUE, DEFAULT_PERROUNDLPITERLIMITFACTOR, -1.0, SCIP_REAL_MAX, NULL,
3282 NULL) );
3283
3284 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundcutsfactorroot", "factor w.r.t. number of integer "
3285 "columns for number of cuts separated per separation round in root node", &sepadata->perroundcutsfactorroot,
3286 TRUE, DEFAULT_PERROUNDCUTSFACTORROOT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
3287
3288 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perroundcutsfactor", "factor w.r.t. number of integer "
3289 "columns for number of cuts separated per separation round at a non-root node",
3290 &sepadata->perroundcutsfactor, TRUE, DEFAULT_PERROUNDCUTSFACTOR, 0.0, SCIP_REAL_MAX, NULL, NULL) );
3291
3292 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/totalcutsfactor", "factor w.r.t. number of integer "
3293 "columns for total number of cuts separated", &sepadata->totalcutsfactor, TRUE, DEFAULT_TOTALCUTSFACTOR, 0.0,
3294 SCIP_REAL_MAX, NULL, NULL) );
3295
3296 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparaminit", "initial value of the mu parameter (factor "
3297 "for step length)", &sepadata->muparaminit, TRUE, DEFAULT_MUPARAMINIT, 0.0, 100.0, NULL, NULL) );
3298
3299 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparamlb", "lower bound of the mu parameter (factor for "
3300 "step length)", &sepadata->muparamlb, TRUE, DEFAULT_MUPARAMLB, 0.0, 1.0, NULL, NULL) );
3301
3302 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muparamub", "upper bound of the mu parameter (factor for "
3303 "step length)", &sepadata->muparamub, TRUE, DEFAULT_MUPARAMUB, 1.0, 10.0, NULL, NULL) );
3304
3305 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/mubacktrackfactor", "factor of mu while backtracking the "
3306 "mu parameter (factor for step length)", &sepadata->mubacktrackfactor, TRUE, DEFAULT_MUBACKTRACKFACTOR,
3307 0.0, 1.0, NULL, NULL) );
3308
3309 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab1factor", "factor of mu parameter (factor for step "
3310 "length) for larger increment" , &sepadata->muslab1factor, TRUE, DEFAULT_MUSLAB1FACTOR, 0.0, SCIP_REAL_MAX, NULL,
3311 NULL));
3312
3313 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab2factor", "factor of mu parameter (factor for step "
3314 "length) for smaller increment", &sepadata->muslab2factor, TRUE, DEFAULT_MUSLAB2FACTOR, 0.0, SCIP_REAL_MAX, NULL,
3315 NULL) );
3316
3317 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/muslab3factor", "factor of mu parameter (factor for step "
3318 "length) for reduction", &sepadata->muslab3factor, TRUE, DEFAULT_MUSLAB3FACTOR, 0.0, SCIP_REAL_MAX, NULL, NULL) );
3319
3320 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/deltaslab1ub", "factor of delta deciding larger increment "
3321 "of mu parameter (factor for step length)", &sepadata->deltaslab1ub, TRUE, DEFAULT_DELTASLAB1UB, 0.0, 1.0,
3322 NULL, NULL) );
3323
3324 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/deltaslab2ub", "factor of delta deciding smaller "
3325 "increment of mu parameter (factor for step length)", &sepadata->deltaslab2ub, TRUE, DEFAULT_DELTASLAB2UB,
3326 0.0, 1.0, NULL, NULL) );
3327
3328 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/ubparamposfactor", "factor for positive upper bound used "
3329 "as an estimate for the optimal Lagrangian dual value", &sepadata->ubparamposfactor, TRUE, DEFAULT_UBPARAMPOSFACTOR,
3330 1.0, 100.0, NULL, NULL) );
3331
3332 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/ubparamnegfactor", "factor for negative upper bound used "
3333 "as an estimate for the optimal Lagrangian dual value", &sepadata->ubparamnegfactor, TRUE, DEFAULT_UBPARAMNEGFACTOR,
3334 0.0, 1.0, NULL, NULL) );
3335
3336 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perrootlpiterfactor", "factor w.r.t. root node LP "
3337 "iterations for iteration limit of each separating LP (negative for no limit)",
3338 &sepadata->perrootlpiterfactor, TRUE, DEFAULT_PERROOTLPITERFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) );
3339
3340 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/perlpiterfactor", "factor w.r.t. non-root node LP "
3341 "iterations for iteration limit of each separating LP (negative for no limit)", &sepadata->perlpiterfactor, TRUE,
3342 DEFAULT_PERLPITERFACTOR, -1.0, SCIP_REAL_MAX, NULL, NULL) );
3343
3344 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/cutsfilterfactor", "fraction of generated cuts per "
3345 "explored basis to accept from separator", &sepadata->cutsfilterfactor, TRUE, DEFAULT_CUTSFILTERFACTOR, 0.0,
3346 1.0, NULL, NULL) );
3347
3348 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusinit", "initial radius of the ball used in "
3349 "stabilization of Lagrangian multipliers", &sepadata->radiusinit, TRUE, DEFAULT_RADIUSINIT, 0.0, 1.0, NULL,
3350 NULL) );
3351
3352 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusmax", "maximum radius of the ball used in "
3353 "stabilization of Lagrangian multipliers", &sepadata->radiusmax, TRUE, DEFAULT_RADIUSMAX, 0.0, SCIP_REAL_MAX,
3354 NULL, NULL) );
3355
3356 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusmin", "minimum radius of the ball used in "
3357 "stabilization of Lagrangian multipliers", &sepadata->radiusmin, TRUE, DEFAULT_RADIUSMIN, 0.0, SCIP_REAL_MAX,
3358 NULL, NULL) );
3359
3360 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/constant", "a constant for stablity center based "
3361 "stabilization of Lagrangian multipliers", &sepadata->constant, TRUE, DEFAULT_CONST, 2.0, SCIP_REAL_MAX,
3362 NULL, NULL) );
3363
3364 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/radiusupdateweight", "multiplier to evaluate cut "
3365 "violation score used for updating ball radius", &sepadata->radiusupdateweight, TRUE,
3366 DEFAULT_RADIUSUPDATEWEIGHT, 0.0, 1.0, NULL, NULL) );
3367
3368 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/dualdegeneracyratethreshold", "minimum dual degeneracy "
3369 "rate for separator execution", &sepadata->dualdegeneracyratethreshold, FALSE,
3370 DEFAULT_DUALDEGENERACYRATETHRESHOLD, 0.0, 1.0, NULL, NULL) );
3371
3372 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/varconsratiothreshold", "minimum variable-constraint "
3373 "ratio on optimal face for separator execution", &sepadata->varconsratiothreshold, FALSE,
3374 DEFAULT_VARCONSRATIOTHRESHOLD, 1.0, SCIP_REAL_MAX, NULL, NULL) );
3375
3376 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/muparamconst", "is the mu parameter (factor for step "
3377 "length) constant?" , &sepadata->muparamconst, TRUE, DEFAULT_MUPARAMCONST, NULL, NULL) );
3378
3379 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/separaterows", "separate rows with integral slack?",
3380 &sepadata->separaterows, TRUE, DEFAULT_SEPARATEROWS, NULL, NULL) );
3381
3382 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/sortcutoffsol", "sort fractional integer columns"
3383 "based on fractionality?" , &sepadata->sortcutoffsol, TRUE, DEFAULT_SORTCUTOFFSOL, NULL, NULL) );
3384
3385 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/sidetypebasis", "choose side types of row (lhs/rhs) "
3386 "based on basis information?", &sepadata->sidetypebasis, TRUE, DEFAULT_SIDETYPEBASIS, NULL, NULL) );
3387
3388 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/dynamiccuts", "should generated cuts be removed from "
3389 "LP if they are no longer tight?", &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) );
3390
3391 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/makeintegral", "try to scale all cuts to integral "
3392 "coefficients?", &sepadata->makeintegral, TRUE, DEFAULT_MAKEINTEGRAL, NULL, NULL) );
3393
3394 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/forcecuts", "force cuts to be added to the LP?",
3395 &sepadata->forcecuts, TRUE, DEFAULT_FORCECUTS, NULL, NULL) );
3396
3397 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/delayedcuts", "should cuts be added to the delayed cut "
3398 "pool", &sepadata->delayedcuts, TRUE, DEFAULT_DELAYEDCUTS, NULL, NULL) );
3399
3400 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/allowlocal", "should locally valid cuts be generated?",
3401 &sepadata->allowlocal, TRUE, DEFAULT_ALLOWLOCAL, NULL, NULL) );
3402
3403 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/aggregatecuts", "aggregate all generated cuts using the "
3404 "Lagrangian multipliers?", &sepadata->aggregatecuts, TRUE, DEFAULT_AGGREGATECUTS, NULL, NULL) );
3405
3406 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrounds", "maximal number of separation rounds per node "
3407 "(-1: unlimited)", &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
3408
3409 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxroundsroot", "maximal number of separation rounds in "
3410 "the root node (-1: unlimited)", &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL,
3411 NULL) );
3412
3413 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/perroundnmaxlpiters", "maximal number of separating LP "
3414 "iterations per separation round (-1: unlimited)", &sepadata->perroundnmaxlpiters, FALSE,
3415 DEFAULT_PERROUNDMAXLPITERS, -1, INT_MAX, NULL, NULL) );
3416
3417 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxcutsperlp", "maximal number of cuts separated per "
3418 "Lagromory LP in the non-root node", &sepadata->nmaxcutsperlp, FALSE, DEFAULT_PERLPMAXCUTS, 0, INT_MAX, NULL,
3419 NULL));
3420
3421 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxcutsperlproot", "maximal number of cuts separated per "
3422 "Lagromory LP in the root node", &sepadata->nmaxcutsperlproot, FALSE, DEFAULT_PERLPMAXCUTSROOT, 0, INT_MAX,
3423 NULL, NULL));
3424
3425 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxmainiters", "maximal number of main loop iterations of "
3426 "the relax-and-cut algorithm", &sepadata->nmaxmainiters, TRUE, DEFAULT_MAXMAINITERS, 0, INT_MAX, NULL, NULL)
3427 );
3428
3429 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxsubgradientiters", "maximal number of subgradient loop "
3430 "iterations of the relax-and-cut algorithm", &sepadata->nmaxsubgradientiters, TRUE,
3431 DEFAULT_MAXSUBGRADIENTITERS, 0, INT_MAX, NULL, NULL) );
3432
3433 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/cutgenfreq", "frequency of subgradient iterations for "
3434 "generating cuts", &sepadata->cutgenfreq, TRUE, DEFAULT_CUTGENFREQ, 0, INT_MAX, NULL, NULL) );
3435
3436 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/cutaddfreq", "frequency of subgradient iterations for "
3437 "adding cuts to objective function", &sepadata->cutaddfreq, TRUE, DEFAULT_CUTADDFREQ, 0, INT_MAX, NULL, NULL) );
3438
3439 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxlagrangianvalsforavg", "maximal number of iterations "
3440 "for rolling average of Lagrangian value", &sepadata->nmaxlagrangianvalsforavg, TRUE,
3441 DEFAULT_MAXLAGRANGIANVALSFORAVG, 0, INT_MAX, NULL, NULL) );
3442
3443 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/nmaxconsecitersformuupdate", "consecutive number of "
3444 "iterations used to determine if mu needs to be backtracked", &sepadata->nmaxconsecitersformuupdate, TRUE,
3445 DEFAULT_MAXCONSECITERSFORMUUPDATE, 0, INT_MAX, NULL, NULL) );
3446
3447 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/projectiontype", "the ball into which the Lagrangian multipliers "
3448 "are projected for stabilization (0: no projection, 1: L1-norm ball projection, 2: L2-norm ball projection, 3: "
3449 "L_inf-norm ball projection)", &sepadata->projectiontype, TRUE, DEFAULT_PROJECTIONTYPE, 0, 3, NULL, NULL));
3450
3451 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/stabilitycentertype", "type of stability center for "
3452 "taking weighted average of Lagrangian multipliers for stabilization (0: no weighted stabilization, 1: best "
3453 "Lagrangian multipliers)", &sepadata->stabilitycentertype, TRUE, DEFAULT_STABILITYCENTERTYPE, 0, 1, NULL,
3454 NULL) );
3455
3456 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/optimalfacepriority", "priority of the optimal face for "
3457 "separator execution (0: low priority, 1: medium priority, 2: high priority)",
3458 &sepadata->optimalfacepriority, TRUE, DEFAULT_OPTIMALFACEPRIORITY, 0, 2, NULL, NULL) );
3459
3460 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/minrestart", "minimum restart round for separator "
3461 "execution (0: from beginning of the instance solving, >= n with n >= 1: from restart round n)",
3462 &sepadata->minrestart, TRUE, DEFAULT_MINRESTART, 0, INT_MAX, NULL, NULL) );
3463
3464 return SCIP_OKAY;
3465 }
3466