source: trunk/Cbc/examples/modk.c @ 2496

Last change on this file since 2496 was 2474, checked in by unxusr, 7 months ago

fixed modk cuts

File size: 6.4 KB
Line 
1// $Id: crew.cpp 2469 2019-01-06 23:17:46Z unxusr $
2// Copyright (C) 2005, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7/*! \file modk.c
8  \brief Example of cut generator using callbacks in the CBC C API
9
10  Tries to generate mod-k cuts, i.e. for each constraint
11  rounds it down after multiplying it by {1 ... k-1} / k
12
13*/
14
15#include <stdio.h>
16#include <stdlib.h>
17#include <limits.h>
18#include <assert.h>
19#include <float.h>
20#include <ctype.h>
21#include <math.h>
22#include <string.h>
23#include <Cbc_C_Interface.h>
24
25#define MIN(a,b) ( (a<b) ? (a) : (b) )
26
27static int dblcmp( const void *pd1, const void *pd2 );
28static char dbl_is_integer( const double value );
29static int intcmp( const void *pd1, const void *pd2 );
30static int mcd( int n, const int values[] );
31static int dbl_as_int( double value );
32static double try_cut( void *osiSolver, int row, const double mult, int num, int den, 
33    int *nzCut, int *idx, double *coef, double *rhs, int *icoef );
34
35void cutcallback( void *osiSolver, void *osiCuts, void *appData )
36{
37  printf("entered callback\n");
38  int m = Osi_getNumRows( osiSolver );
39  int n = Osi_getNumCols( osiSolver );
40
41  // sorted coefs and rhs values
42  int *scoef = (int *) malloc( sizeof(int)*(n+1) );
43
44  // different integer positive coeffs (values for k)
45  int *dcoef = (int *) malloc( sizeof(int)*(n+1) );
46
47  // integer coefficients for computing the cut
48  int *icoef = (int*) malloc( sizeof(int)*n );
49
50  // double coefficients for storing the cut
51  double *coef = (double*) malloc(sizeof(double)*n);
52
53  // variables in this cut
54  int *cidx = (int*) malloc(sizeof(int)*n);
55
56  for ( int i=0 ; i<m ; ++i )
57  {
58    // checking different coefficients
59    int nz = Osi_getRowNz( osiSolver, i );
60
61    double rhs = Osi_getRowRHS( osiSolver, i );
62    if (!dbl_is_integer(rhs))
63      continue;
64
65    int irhs = dbl_as_int( rhs );
66
67    char allInt = 1;
68    const double *coeffs = Osi_getRowCoeffs( osiSolver, i );
69    for ( int j=0 ; (j<nz) ; ++j )
70    {
71      if (!dbl_is_integer(coeffs[j]))
72      {
73        allInt = 0;
74        break;
75      }
76      scoef[j] = abs(dbl_as_int(coeffs[j]));
77    }
78
79    if (!allInt)
80      continue;
81
82    scoef[nz++] = abs(irhs);
83
84    qsort( scoef, nz, sizeof(int), intcmp );
85
86
87    int prev = INT_MAX;
88    int ndiff = 0;
89
90    // checking for different values for k
91    for ( int j=0 ; (j<nz) ; ++j )
92    {
93      if (prev!=scoef[j])
94      {
95        prev = scoef[j];
96        dcoef[ndiff++] = scoef[j];
97      } // different coefs
98    } // all columns
99
100    if (ndiff<=1)
101      continue;
102
103    int nRem = 0;
104
105    // removing those which are divisors of larger values
106    for ( int j=0 ; (j<ndiff-1) ; ++j )
107    {
108      for ( int jl=j+1 ; (jl<ndiff) ; ++jl )
109      {
110        if (!dcoef[j])
111          continue;
112        if ( dcoef[jl] % dcoef[j] == 0 )
113        {
114          dcoef[j] = INT_MAX;
115          ++nRem;
116          break;
117        }
118      } // larger
119    } // all non-zeros
120
121    qsort( dcoef, ndiff, sizeof(int), intcmp );
122    ndiff -= nRem;
123
124    for ( int ik=0 ; ik<ndiff ; ++ ik )
125    {
126      int k = dcoef[ik];
127      for ( int num=1 ; (num<k) ; ++num )
128      {
129        int nzCut = 0;
130        double rhsCut = 0.0;
131
132        int nMult, startMult = 0;
133        double mults[] = { 1.0, -1.0 };
134
135        char sense = toupper(Osi_getRowSense( osiSolver, i ));
136
137        switch (sense)
138        {
139          case 'E':
140            {
141              nMult = 2;
142              break;
143            }
144          case 'G':
145            {
146              nMult = 1;
147              startMult = 1;
148              break;
149            }
150          case 'L':
151            {
152              nMult = 1;
153              startMult = 0;
154              break;
155            }
156          default:
157            {
158              fprintf( stderr, "unknow sense!");
159              abort();
160            }
161        }
162
163        for ( int iMult=startMult ; (iMult<nMult) ; ++iMult )
164        {
165          const double mult = mults[iMult];
166          const double viol = try_cut( osiSolver, i, mult, num, k, &nzCut,
167              cidx, coef, &rhsCut, icoef );
168          if ( viol<0.001 )
169            continue;
170
171          OsiCuts_addRowCut( osiCuts, nzCut, cidx, coef, 'L', rhsCut );
172        } // constraint in the form <=
173      } // every numerator < k
174    } // every denominator k
175  } // every row
176
177  free( scoef );
178  free( dcoef );
179  free( cidx );
180  free( coef );
181  free( icoef );
182}
183
184static int dblcmp( const void *pd1, const void *pd2 )
185{
186  const double *d1 = (const double *) pd1;
187  const double *d2 = (const double *) pd2;
188
189  if ((*d1)<(*d2))
190    return -1;
191  if ((*d1)>(*d2))
192    return 1;
193
194  return 0;
195}
196
197static char dbl_is_integer( const double value )
198{
199  double rv = floor( value+0.5 );
200
201  return fabs(rv-value) < 1e-6;
202}
203
204static int intcmp( const void *pi1, const void *pi2 )
205{
206  const int *i1 = (const int *) pi1;
207  const int *i2 = (const int *) pi2;
208
209  return ( (*i1) - (*i2) );
210}
211
212static int mcd( int n, const int values[] )
213{
214  //for ( int j=0 ; (j<n) ; ++j )
215  //
216  return 0;
217}
218
219static int dbl_as_int( double value )
220{
221  return (int) floor(value+0.5);
222}
223
224static double try_cut( void *osiSolver, int row, const double mult, int num, int den, 
225    int *nzCut, int *idx, double *coef, double *rhs, int *icoef )
226{
227  int nz = Osi_getRowNz( osiSolver, row );
228  const int *rIdx = Osi_getRowIndices( osiSolver, row );
229  const double *rCoef = Osi_getRowCoeffs( osiSolver, row );
230  int imult = dbl_as_int(mult);
231  int irhs = dbl_as_int( imult*Osi_getRowRHS( osiSolver, row ) );
232  for ( int j=0 ; j<nz ; ++j )
233    icoef[j] = dbl_as_int( imult*rCoef[j] );
234
235  const double *x = Osi_getColSolution( osiSolver );
236
237  double sumLHS = 0.0;
238
239  (*nzCut) = 0;
240  for ( int j=0 ; j<nz ; ++j )
241  {
242    int col = rIdx[j];
243    assert( col >= 0 && col<Osi_getNumCols(osiSolver) );
244
245    int cc = (icoef[j]*num)/den;
246    if (!cc)
247      continue;
248
249    idx[(*nzCut)] = col;
250    coef[(*nzCut)] = cc;
251
252    sumLHS += ((double)cc) * x[col];
253
254    (*nzCut)++;
255  }
256
257  *rhs = (irhs*num)/den;
258
259  return (*rhs) - sumLHS;
260}
261
262int main( int argc, char **argv )
263{
264  if (argc<2)
265  {
266    fprintf( stderr, "enter mps instance name.\n");
267    exit(1);
268  }
269
270  Cbc_Model *mip = Cbc_newModel();
271
272  Cbc_readMps( mip, argv[1] );
273
274  Cbc_setParameter( mip, "cuts", "off" );
275  Cbc_setParameter( mip, "gomory", "ifmove" );
276  Cbc_setParameter( mip, "preprocess", "off" );
277
278  Cbc_addCutCallback( mip,  cutcallback, "modk", NULL );
279
280
281  Cbc_solve( mip );
282
283  Cbc_deleteModel( mip );
284}
285
286/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
287*/
Note: See TracBrowser for help on using the repository browser.