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   presol_stuffing.c
26   	 * @ingroup DEFPLUGINS_PRESOL
27   	 * @brief  fix singleton continuous variables
28   	 * @author Dieter Weninger
29   	 *
30   	 * Investigate singleton continuous variables if one can be fixed at a bound.
31   	 *
32   	 * @todo enhancement from singleton continuous variables to continuous
33   	 *       variables with only one lock in a common row
34   	 */
35   	
36   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
37   	
38   	#include "blockmemshell/memory.h"
39   	#include "scip/presol_stuffing.h"
40   	#include "scip/pub_matrix.h"
41   	#include "scip/pub_message.h"
42   	#include "scip/pub_misc_sort.h"
43   	#include "scip/pub_presol.h"
44   	#include "scip/pub_var.h"
45   	#include "scip/scip_general.h"
46   	#include "scip/scip_mem.h"
47   	#include "scip/scip_message.h"
48   	#include "scip/scip_nlp.h"
49   	#include "scip/scip_numerics.h"
50   	#include "scip/scip_presol.h"
51   	#include "scip/scip_pricer.h"
52   	#include "scip/scip_prob.h"
53   	#include "scip/scip_probing.h"
54   	#include "scip/scip_var.h"
55   	#include <string.h>
56   	
57   	#define PRESOL_NAME            "stuffing"
58   	#define PRESOL_DESC            "fix redundant singleton continuous variables"
59   	#define PRESOL_PRIORITY             -100     /**< priority of the presolver (>= 0: before, < 0: after constraint handlers) */
60   	#define PRESOL_MAXROUNDS               0     /**< maximal number of presolving rounds the presolver participates in (-1: no limit) */
61   	#define PRESOL_TIMING           SCIP_PRESOLTIMING_EXHAUSTIVE /* timing of the presolver (fast, medium, or exhaustive) */
62   	
63   	/** type of fixing direction */
64   	enum Fixingdirection
65   	{
66   	   FIXATLB = -1,         /**< fix variable at lower bound */
67   	   NOFIX   =  0,         /**< do not fix variable */
68   	   FIXATUB =  1          /**< fix variable at upper bound */
69   	};
70   	typedef enum Fixingdirection FIXINGDIRECTION;
71   	
72   	/*
73   	 * Local methods
74   	 */
75   	
76   	/** try to fix singleton continuous variables */
77   	static
78   	SCIP_RETCODE singletonColumnStuffing(
79   	   SCIP*                 scip,               /**< SCIP main data structure */
80   	   SCIP_MATRIX*          matrix,             /**< matrix containing the constraints */
81   	   FIXINGDIRECTION*      varstofix,          /**< array holding fixing information */
82   	   int*                  nfixings            /**< number of possible fixings */
83   	   )
84   	{
85   	   SCIP_Real* valpnt;
86   	   SCIP_Real* colratios;
87   	   SCIP_Real* colcoeffs;
88   	   SCIP_Bool* rowprocessed;
89   	   int* rowpnt;
90   	   int* rowend;
91   	   int* colindices;
92   	   int* dummy;
93   	   SCIP_Bool* swapped;
94   	   SCIP_Real upperconst;
95   	   SCIP_Real lowerconst;
96   	   SCIP_Real rhs;
97   	   SCIP_Bool tryfixing;
98   	   int idx;
99   	   int col;
100  	   int row;
101  	   int fillcnt;
102  	   int k;
103  	   int nrows;
104  	   int ncols;
105  	
106  	   assert(scip != NULL);
107  	   assert(matrix != NULL);
108  	   assert(varstofix != NULL);
109  	   assert(nfixings != NULL);
110  	
111  	   nrows = SCIPmatrixGetNRows(matrix);
112  	   ncols = SCIPmatrixGetNColumns(matrix);
113  	
114  	   SCIP_CALL( SCIPallocBufferArray(scip, &colindices, ncols) );
115  	   SCIP_CALL( SCIPallocBufferArray(scip, &colratios, ncols) );
116  	   SCIP_CALL( SCIPallocBufferArray(scip, &colcoeffs, ncols) );
117  	   SCIP_CALL( SCIPallocBufferArray(scip, &dummy, ncols) );
118  	
119  	   SCIP_CALL( SCIPallocBufferArray(scip, &swapped, ncols) );
120  	   BMSclearMemoryArray(swapped, ncols);
121  	
122  	   SCIP_CALL( SCIPallocBufferArray(scip, &rowprocessed, nrows) );
123  	   BMSclearMemoryArray(rowprocessed, nrows);
124  	
125  	   for( col = 0; col < ncols; col++ )
126  	   {
127  	      /* consider only rows with minimal one continuous singleton column */
128  	      if( SCIPmatrixGetColNNonzs(matrix, col) == 1
129  	         && SCIPvarGetType(SCIPmatrixGetVar(matrix, col)) == SCIP_VARTYPE_CONTINUOUS
130  	         && SCIPvarGetNLocksUpType(SCIPmatrixGetVar(matrix, col), SCIP_LOCKTYPE_MODEL) == SCIPmatrixGetColNUplocks(matrix, col)
131  	         && SCIPvarGetNLocksDownType(SCIPmatrixGetVar(matrix, col), SCIP_LOCKTYPE_MODEL) == SCIPmatrixGetColNDownlocks(matrix, col) )
132  	      {
133  	         row = *(SCIPmatrixGetColIdxPtr(matrix, col));
134  	         if( rowprocessed[row] )
135  	            continue;
136  	
137  	         rowprocessed[row] = TRUE;
138  	
139  	         /* treat all >= rows from matrix, but internally we transform to <= relation */
140  	         if( SCIPmatrixIsRowRhsInfinity(matrix, row) )
141  	         {
142  	            fillcnt = 0;
143  	            tryfixing = TRUE;
144  	            upperconst = 0.0;
145  	            lowerconst = 0.0;
146  	            rhs = -SCIPmatrixGetRowLhs(matrix, row);
147  	
148  	            rowpnt = SCIPmatrixGetRowIdxPtr(matrix, row);
149  	            rowend = rowpnt + SCIPmatrixGetRowNNonzs(matrix, row);
150  	            valpnt = SCIPmatrixGetRowValPtr(matrix, row);
151  	
152  	            for( ; (rowpnt < rowend); rowpnt++, valpnt++ )
153  	            {
154  	               SCIP_Real coef;
155  	               SCIP_VAR* var;
156  	               SCIP_Real lb;
157  	               SCIP_Real ub;
158  	
159  	               coef = -(*valpnt);
160  	               idx = *rowpnt;
161  	               var = SCIPmatrixGetVar(matrix, idx);
162  	               lb = SCIPvarGetLbGlobal(var);
163  	               ub = SCIPvarGetUbGlobal(var);
164  	
165  	               /* we need to check if this is a singleton continuous variable and
166  	                * all constraints containing this variable are present inside
167  	                * the mixed integer linear matrix
168  	                */
169  	               if( SCIPmatrixGetColNNonzs(matrix, idx) == 1
170  	                  && (SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) + SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL)) == 1
171  	                  && SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
172  	               {
173  	                  if( SCIPisLT(scip, SCIPvarGetObj(var), 0.0) && SCIPisGT(scip, coef, 0.0) )
174  	                  {
175  	                     /* case 1: obj < 0 and coef > 0 */
176  	                     if( SCIPisInfinity(scip, -lb) )
177  	                     {
178  	                        tryfixing = FALSE;
179  	                        break;
180  	                     }
181  	
182  	                     upperconst += coef * lb;
183  	                     lowerconst += coef * lb;
184  	                     colratios[fillcnt] = SCIPvarGetObj(var) / coef;
185  	                     colindices[fillcnt] = idx;
186  	                     colcoeffs[fillcnt] = coef;
187  	                     fillcnt++;
188  	                  }
189  	                  else if( SCIPisGT(scip, SCIPvarGetObj(var), 0.0) && SCIPisLT(scip, coef, 0.0) )
190  	                  {
191  	                     /* case 2: obj > 0 and coef < 0 */
192  	                     if( SCIPisInfinity(scip, ub) )
193  	                     {
194  	                        tryfixing = FALSE;
195  	                        break;
196  	                     }
197  	
198  	                     /* multiply column by (-1) to become case 1.
199  	                      * now bounds are swapped: ub := -lb, lb := -ub
200  	                      */
201  	                     swapped[idx] = TRUE;
202  	                     upperconst += coef * ub;
203  	                     lowerconst += coef * ub;
204  	                     colratios[fillcnt] = SCIPvarGetObj(var) / coef;
205  	                     colindices[fillcnt] = idx;
206  	                     colcoeffs[fillcnt] = -coef;
207  	                     fillcnt++;
208  	                  }
209  	                  else if( SCIPisGE(scip, SCIPvarGetObj(var), 0.0) && SCIPisGE(scip, coef, 0.0) )
210  	                  {
211  	                     /* case 3: obj >= 0 and coef >= 0 is handled by duality fixing.
212  	                      *  we only consider the lower bound for the constants
213  	                      */
214  	                     if( SCIPisInfinity(scip, -lb) )
215  	                     {
216  	                        /* maybe unbounded */
217  	                        tryfixing = FALSE;
218  	                        break;
219  	                     }
220  	
221  	                     upperconst += coef * lb;
222  	                     lowerconst += coef * lb;
223  	                  }
224  	                  else
225  	                  {
226  	                     /* case 4: obj <= 0 and coef <= 0 is also handled by duality fixing.
227  	                      * we only consider the upper bound for the constants
228  	                      */
229  	                     assert(SCIPisLE(scip, SCIPvarGetObj(var), 0.0));
230  	                     assert(SCIPisLE(scip, coef, 0.0));
231  	
232  	                     if( SCIPisInfinity(scip, ub) )
233  	                     {
234  	                        /* maybe unbounded */
235  	                        tryfixing = FALSE;
236  	                        break;
237  	                     }
238  	
239  	                     upperconst += coef * ub;
240  	                     lowerconst += coef * ub;
241  	                  }
242  	               }
243  	               else
244  	               {
245  	                  /* consider contribution of discrete variables, non-singleton
246  	                   * continuous variables and variables with more than one lock
247  	                   */
248  	                  if( SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub) )
249  	                  {
250  	                     tryfixing = FALSE;
251  	                     break;
252  	                  }
253  	
254  	                  if( coef > 0 )
255  	                  {
256  	                     upperconst += coef * ub;
257  	                     lowerconst += coef * lb;
258  	                  }
259  	                  else
260  	                  {
261  	                     upperconst += coef * lb;
262  	                     lowerconst += coef * ub;
263  	                  }
264  	               }
265  	            }
266  	
267  	            if( tryfixing )
268  	            {
269  	               SCIPsortRealRealIntInt(colratios, colcoeffs, colindices, dummy, fillcnt);
270  	
271  	               /* verify which singleton continuous variable can be fixed */
272  	               for( k = 0; k < fillcnt; k++ )
273  	               {
274  	                  SCIP_VAR* var;
275  	                  SCIP_Real lb;
276  	                  SCIP_Real ub;
277  	                  SCIP_Real delta;
278  	
279  	                  idx = colindices[k];
280  	                  var = SCIPmatrixGetVar(matrix, idx);
281  	                  lb = SCIPvarGetLbGlobal(var);
282  	                  ub = SCIPvarGetUbGlobal(var);
283  	
284  	                  /* stop fixing if variable bounds are not finite */
285  	                  if( SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub) )
286  	                     break;
287  	
288  	                  assert(SCIPmatrixGetColNNonzs(matrix, idx) == 1);
289  	                  assert(SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS);
290  	                  assert((SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) + SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL)) == 1);
291  	                  assert(colcoeffs[k] >= 0);
292  	
293  	                  /* calculate the change in the row activities if this variable changes
294  	                   * its value from to its other bound
295  	                   */
296  	                  if( swapped[idx] )
297  	                     delta = -(lb - ub) * colcoeffs[k];
298  	                  else
299  	                     delta =  (ub - lb) * colcoeffs[k];
300  	
301  	                  assert(delta >= 0);
302  	
303  	                  if( SCIPisLE(scip, delta, rhs - upperconst) )
304  	                  {
305  	                     if( swapped[idx] )
306  	                        varstofix[idx] = FIXATLB;
307  	                     else
308  	                        varstofix[idx] = FIXATUB;
309  	
310  	                     (*nfixings)++;
311  	                  }
312  	                  else if( SCIPisLE(scip, rhs, lowerconst) )
313  	                  {
314  	                     if( swapped[idx] )
315  	                        varstofix[idx] = FIXATUB;
316  	                     else
317  	                        varstofix[idx] = FIXATLB;
318  	
319  	                     (*nfixings)++;
320  	                  }
321  	
322  	                  upperconst += delta;
323  	                  lowerconst += delta;
324  	               }
325  	            }
326  	         }
327  	      }
328  	   }
329  	
330  	   SCIPfreeBufferArray(scip, &rowprocessed);
331  	   SCIPfreeBufferArray(scip, &swapped);
332  	   SCIPfreeBufferArray(scip, &dummy);
333  	   SCIPfreeBufferArray(scip, &colcoeffs);
334  	   SCIPfreeBufferArray(scip, &colratios);
335  	   SCIPfreeBufferArray(scip, &colindices);
336  	
337  	   return SCIP_OKAY;
338  	}
339  	
340  	/*
341  	 * Callback methods of presolver
342  	 */
343  	
344  	/** copy method for constraint handler plugins (called when SCIP copies plugins) */
345  	static
346  	SCIP_DECL_PRESOLCOPY(presolCopyStuffing)
347  	{  /*lint --e{715}*/
348  	   assert(scip != NULL);
349  	   assert(presol != NULL);
350  	   assert(strcmp(SCIPpresolGetName(presol), PRESOL_NAME) == 0);
351  	
352  	   /* call inclusion method of presolver */
353  	   SCIP_CALL( SCIPincludePresolStuffing(scip) );
354  	
355  	   return SCIP_OKAY;
356  	}
357  	
358  	/** execution method of presolver */
359  	static
360  	SCIP_DECL_PRESOLEXEC(presolExecStuffing)
361  	{  /*lint --e{715}*/
362  	   SCIP_MATRIX* matrix;
363  	   SCIP_Bool initialized;
364  	   SCIP_Bool complete;
365  	   SCIP_Bool infeasible;
366  	
367  	   assert(result != NULL);
368  	   *result = SCIP_DIDNOTRUN;
369  	
370  	   if( (SCIPgetStage(scip) != SCIP_STAGE_PRESOLVING) || SCIPinProbing(scip) || SCIPisNLPEnabled(scip) )
371  	      return SCIP_OKAY;
372  	
373  	   if( SCIPgetNContVars(scip) == 0 || SCIPisStopped(scip) || SCIPgetNActivePricers(scip) > 0 )
374  	      return SCIP_OKAY;
375  	
376  	   if( !SCIPallowStrongDualReds(scip) )
377  	      return SCIP_OKAY;
378  	
379  	   *result = SCIP_DIDNOTFIND;
380  	
381  	   matrix = NULL;
382  	   SCIP_CALL( SCIPmatrixCreate(scip, &matrix, FALSE, &initialized, &complete, &infeasible,
383  	      naddconss, ndelconss, nchgcoefs, nchgbds, nfixedvars) );
384  	
385  	   /* if infeasibility was detected during matrix creation, return here */
386  	   if( infeasible )
387  	   {
388  	      if( initialized )
389  	         SCIPmatrixFree(scip, &matrix);
390  	
391  	      *result = SCIP_CUTOFF;
392  	      return SCIP_OKAY;
393  	   }
394  	
395  	   if( initialized )
396  	   {
397  	      FIXINGDIRECTION* varstofix;
398  	      int nfixings;
399  	      int ncols;
400  	
401  	      nfixings = 0;
402  	      ncols = SCIPmatrixGetNColumns(matrix);
403  	
404  	      SCIP_CALL( SCIPallocBufferArray(scip, &varstofix, ncols) );
405  	      BMSclearMemoryArray(varstofix, ncols);
406  	
407  	      SCIP_CALL( singletonColumnStuffing(scip, matrix, varstofix, &nfixings) );
408  	
409  	      if( nfixings > 0 )
410  	      {
411  	         int v;
412  	         int oldnfixedvars;
413  	
414  	         oldnfixedvars = *nfixedvars;
415  	
416  	         /* look for fixable variables */
417  	         for( v = ncols - 1; v >= 0; --v )
418  	         {
419  	            SCIP_Bool fixed;
420  	            SCIP_VAR* var;
421  	
422  	            var = SCIPmatrixGetVar(matrix, v);
423  	
424  	            if( varstofix[v] == FIXATLB )
425  	            {
426  	               SCIP_Real lb;
427  	
428  	               assert(SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == SCIPmatrixGetColNUplocks(matrix, v)
429  	                  && SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == SCIPmatrixGetColNDownlocks(matrix, v));
430  	
431  	               lb = SCIPvarGetLbGlobal(var);
432  	               assert(SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS);
433  	
434  	               /* avoid fixings to infinite values */
435  	               assert(!SCIPisInfinity(scip, -lb));
436  	
437  	               SCIPdebugMsg(scip, "Fix variable %s at lower bound %.15g\n", SCIPvarGetName(var), lb);
438  	
439  	               /* fix at lower bound */
440  	               SCIP_CALL( SCIPfixVar(scip, var, lb, &infeasible, &fixed) );
441  	               if( infeasible )
442  	               {
443  	                  SCIPdebugMsg(scip, " -> infeasible fixing\n");
444  	                  *result = SCIP_CUTOFF;
445  	
446  	                  break;
447  	               }
448  	               assert(fixed);
449  	               (*nfixedvars)++;
450  	            }
451  	            else if( varstofix[v] == FIXATUB )
452  	            {
453  	               SCIP_Real ub;
454  	
455  	               assert(SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == SCIPmatrixGetColNUplocks(matrix, v)
456  	                  && SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == SCIPmatrixGetColNDownlocks(matrix, v));
457  	
458  	               ub = SCIPvarGetUbGlobal(var);
459  	               assert(SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS);
460  	
461  	               /* avoid fixings to infinite values */
462  	               assert(!SCIPisInfinity(scip, ub));
463  	
464  	               SCIPdebugMsg(scip, "Fix variable %s at upper bound %.15g\n", SCIPvarGetName(var), ub);
465  	
466  	               /* fix at upper bound */
467  	               SCIP_CALL( SCIPfixVar(scip, var, ub, &infeasible, &fixed) );
468  	               if( infeasible )
469  	               {
470  	                  SCIPdebugMsg(scip, " -> infeasible fixing\n");
471  	                  *result = SCIP_CUTOFF;
472  	
473  	                  break;
474  	               }
475  	               assert(fixed);
476  	
477  	               (*nfixedvars)++;
478  	            }
479  	         }
480  	
481  	         if( *result != SCIP_CUTOFF && *nfixedvars > oldnfixedvars )
482  	            *result = SCIP_SUCCESS;
483  	      }
484  	
485  	      SCIPfreeBufferArray(scip, &varstofix);
486  	   }
487  	
488  	   SCIPmatrixFree(scip, &matrix);
489  	
490  	   return SCIP_OKAY;
491  	}
492  	
493  	/*
494  	 * presolver specific interface methods
495  	 */
496  	
497  	/** creates the stuffing presolver and includes it in SCIP */
498  	SCIP_RETCODE SCIPincludePresolStuffing(
499  	   SCIP*                 scip                /**< SCIP data structure */
500  	   )
501  	{
502  	   SCIP_PRESOL* presol;
503  	
504  	   /* include presolver */
505  	   SCIP_CALL( SCIPincludePresolBasic(scip, &presol, PRESOL_NAME, PRESOL_DESC, PRESOL_PRIORITY, PRESOL_MAXROUNDS,
506  	         PRESOL_TIMING, presolExecStuffing, NULL) );
507  	   SCIP_CALL( SCIPsetPresolCopy(scip, presol, presolCopyStuffing) );
508  	
509  	   return SCIP_OKAY;
510  	}
511  	
512  	/*lint --e{749}*/
513