# source:trunk/Couenne/src/convex/operators/conv-exprTrilinear-gencuts.cpp@487

Last change on this file since 487 was 487, checked in by pbelotti, 3 years ago

adopting new build system in trunk

• Property svn:keywords set to `Author Date Id Revision`
File size: 55.3 KB
RevLine
[458]1/* \$Id\$
2 *
3 * Name:       conv-exprTrilinear-gencuts.cpp.cpp
4 * Source:     GNU C++
5 * Author:     Sonia Cafieri
6 * Purpose:    generate inequalities defining the convex envelope of a
7 *             trilinear monomial
8 * History:    Nov 2010 work started
9 */
10
11#include "CouenneCutGenerator.hpp"
12
13#include "CouenneTypes.hpp"
14#include "CouenneExprMul.hpp"
15#include "CouenneExprTrilinear.hpp"
16#include "CouenneProblem.hpp"
17#include "CouenneExprAux.hpp"
18
19#include <vector>
20
[472]21//#define DEBUG
[458]22using namespace Couenne;
23
[470]24#define EPSILONT 1.e-6
25
[458]26//typedef CouNumber double;
27
28//// permutations of 3 elements
29void permutation3(int **ind,int *ibnd)
30{
31  ind[0][0] = ibnd[0]; ind[0][1] = ibnd[1]; ind[0][2] = ibnd[2];
32  ind[1][0] = ibnd[0]; ind[1][1] = ibnd[2]; ind[1][2] = ibnd[1];
33  ind[2][0] = ibnd[1]; ind[2][1] = ibnd[0]; ind[2][2] = ibnd[2];
34  ind[3][0] = ibnd[1]; ind[3][1] = ibnd[2]; ind[3][2] = ibnd[0];
35  ind[4][0] = ibnd[2]; ind[4][1] = ibnd[0]; ind[4][2] = ibnd[1];
36  ind[5][0] = ibnd[2]; ind[5][1] = ibnd[1]; ind[5][2] = ibnd[0];
37}
38
39
40/// generate convexification cuts for constraint w = x*y*z
41void TriLinCuts (double *vlb, double *vub, int *varIndices,
42                 std::vector <std::vector <int> >    &cutIndices,
43                 std::vector <std::vector <double> > &cutCoeff,
44                 std::vector <double>                &cutLb,
45                 std::vector <double>                &cutUb) {
46
47  // var indices
48  int v1, v2, v3, v4;
49  // number of cuts
50  int defcons_size = 20;
51  // bounds on cuts
52  double *bnd = new double [12];
53
54  CouNumber cf = 1.;
55
56  int **ind;
57  ind = new int*[6];
58  for(int i=0; i<6; i++) {
59    ind[i] = new int[6];
60  }
61
[469]62  int *ibnd;
63  ibnd = new int[3];
64  ibnd[0] = varIndices[0]; ibnd[1] = varIndices[1]; ibnd[2] = varIndices[2];
65#ifdef DEBUG
66std::cout << "ibnd[0] =" << ibnd[0] << "  ibnd[1] =" << ibnd[1] << "  ibnd[2] =" << ibnd[2] << std::endl;
67std::cout << "vlb[ibnd[0]] =" << vlb[ibnd[0]] << "  vub[ibnd[0]] =" << vub[ibnd[0]] << std::endl;
68std::cout << "vlb[ibnd[1]] =" << vlb[ibnd[1]] << "  vub[ibnd[1]] =" << vub[ibnd[1]] << std::endl;
69std::cout << "vlb[ibnd[2]] =" << vlb[ibnd[2]] << "  vub[ibnd[2]] =" << vub[ibnd[2]] << std::endl;
70#endif
71
[458]72  // compute the 6 permutations of the 3 variables
[469]73  permutation3(ind,ibnd);
[458]74
75  int i, flag=0, idx=0;
76  i = 0;
77  while(i < 6 && flag == 0) {
78    if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] <=0 && vub[ind[i][2]] >=0) {
79      idx = i;   // store the index of the permutation satisfying the condition
80      flag = 1;  // this is case 1
81    }
82    i++;
83  }
84  i = 0;
85  while(i < 6 && flag == 0) {
86    if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0
87       && vub[ind[i][1]] >=0 && vub[ind[i][2]] >=0) {
88      idx = i;   // store the index of the permutation satisfying the condition
89      flag = 2;  // this is case 2
90    }
91    i++;
92  }
93  i = 0;
94  while(i < 6 && flag == 0) {
95    if(vlb[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0
96       && vub[ind[i][0]] >=0 && vub[ind[i][1]] >=0 && vub[ind[i][2]] >=0) {
97      idx = i;   // store the index of the permutation satisfying the condition
98      flag = 3;  // this is case 3
99    }
100    i++;
101  }
102  i = 0;
103  while(i < 6 && flag == 0) {
104    if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0
105       && vub[ind[i][1]] >=0 && vub[ind[i][2]] <=0) {
106      idx = i;   // store the index of the permutation satisfying the condition
107      flag = 4;  // this is case 4
108    }
109    i++;
110  }
111  i = 0;
112  while(i < 6 && flag == 0) {
113    if(vlb[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vlb[ind[i][2]] <=0
114       && vub[ind[i][0]] >=0 && vub[ind[i][1]] >=0 && vub[ind[i][2]] <=0) {
115      idx = i;   // store the index of the permutation satisfying the condition
116      flag = 5;  // this is case 5
117    }
118    i++;
119  }
120  i = 0;
121  while(i < 6 && flag == 0) {
122    if(vlb[ind[i][0]] <=0 && vub[ind[i][0]] >=0 && vub[ind[i][1]] <=0 && vub[ind[i][2]] <=0) {
123      idx = i;   // store the index of the permutation satisfying the condition
124      flag = 6;  // this is case 6
125    }
126    i++;
127  }
128  i = 0;
129  while(i < 6 && flag == 0) {
130    if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] >=0) {
131      idx = i;   // store the index of the permutation satisfying the condition
132      flag = 7;  // this is case 7
133    }
134    i++;
135  }
136  i = 0;
137  while(i < 6 && flag == 0) {
138    if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] >=0 && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
139      idx = i;   // store the index of the permutation satisfying the condition
140      flag = 8;  // this is case 8
141    }
142    i++;
143  }
144  i = 0;
145  while(i < 6 && flag == 0) {
146    if(vlb[ind[i][0]] >=0 && vlb[ind[i][1]] <=0 && vub[ind[i][1]] <=0
147       && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
148      idx = i;   // store the index of the permutation satisfying the condition
149      flag = 9;  // this is case 9
150    }
151    i++;
152  }
153  i = 0;
154  while(i < 6 && flag == 0) {
155    if(vlb[ind[i][0]] <=0 && vub[ind[i][0]] <=0 && vlb[ind[i][1]] <=0 && vub[ind[i][1]] <=0
156       && vlb[ind[i][2]] <=0 && vub[ind[i][2]] <=0) {
157      idx = i;   // store the index of the permutation satisfying the condition
158      flag = 10;  // this is case 10
159    }
160    i++;
161  }
162
163  if (flag==0) {
164    std::cout << "ERROR: case not implemented" << std::endl;
165    exit(0);
166  }
167
168  // var indices
169  v1 = ind[idx][0];
170  v2 = ind[idx][1];
171  v3 = ind[idx][2];
172  v4 = varIndices [3];
173
174  // lower and upper bound on variables
175  double xL1 = cf*vlb[v1];  double  xU1 = cf*vub[v1];
176  double xL2 = vlb[v2];  double xU2 = vub[v2];
177  double xL3 = vlb[v3];  double xU3 = vub[v3];
178
179#define prepareVectors(a) {                             \
180                                                        \
181    int size = (int) (cutIndices.size ());              \
182                                                        \
183    for (int i=0; i<a; i++) {                           \
184                                                        \
185      cutIndices. push_back (std::vector <int>    ());  \
186      cutCoeff.   push_back (std::vector <double> ());  \
187                                                        \
188      cutLb.      push_back (-COUENNE_INFINITY);        \
189      cutUb.      push_back ( COUENNE_INFINITY);        \
190                                                        \
191      for (int j=0; j<4; j++) {                         \
192                                                        \
193        cutIndices [size+i].push_back (-1);             \
194        cutCoeff   [size+i].push_back (0.);             \
195      }                                                 \
196    }                                                   \
197}
198
199  /*----------------------------------------------------------------------------------------*/
200
201  // case 1
202  if(flag == 1) {
[469]203#ifdef DEBUG
204    std::cout << " -- case 1 --" << std::endl;
205#endif
[458]206
207    double theta  = xL1*xU2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
208    double theta1 = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
209
210    defcons_size = 12;
211
212    prepareVectors(defcons_size);
213
214    for(int ii = 0; ii < defcons_size; ii++) {
215
216      cutIndices [ii][0] = v1;
217      cutIndices [ii][1] = v2;
218      cutIndices [ii][2] = v3;
219      cutIndices [ii][3] = v4;
220
221      cutCoeff [ii][3] = 1.;
222    }
223
224    cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xU1*xU3; cutCoeff [0][2] = -xU1*xU2;  bnd[0] = - 2.*xU1*xU2*xU3;
225    cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3;
226    cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xL1*xU2*xL3 - xL1*xL2*xL3;
227    cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xU1*xL2;  bnd[3] = - xU1*xL2*xU3 - xU1*xL2*xL3;
228    cutCoeff [4][0] = -xL2*xL3; cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xL1*xL2;  bnd[4] = - xU1*xL2*xL3 - xL1*xL2*xL3;
229    cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -(theta/(xU3-xL3));
230    bnd[5] = (-(theta*xL3)/(xU3-xL3)) - xL1*xU2*xU3 - xU1*xL2*xU3 + xU1*xU2*xL3;
231
232    cutCoeff [6][0] = -xU2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2;  bnd[6] = - 2.*xU1*xU2*xL3;
233    cutCoeff [7][0] = -xL2*xL3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xL2;  bnd[7] = - xU1*xL2*xU3 - xU1*xL2*xL3;
234    cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xL1*xL2;  bnd[8] = - xL1*xU2*xU3 - xL1*xL2*xU3;
235    cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
236    cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xU1*xU3; cutCoeff [10][2] = -xL1*xL2;  bnd[10] = - xU1*xL2*xU3 - xL1*xL2*xU3;
237    cutCoeff [11][0] = -xL2*xL3; cutCoeff [11][1] = -xL1*xL3; cutCoeff [11][2] = -(theta1/(xL3-xU3)); cutCoeff [11][3] = 1.;
238    bnd[11] = (-(theta1*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
239  } // end if case 1
240
241  /*----------------------------------------------------------------------------------------*/
242
243  // case 2
244  if(flag == 2) {
[469]245#ifdef DEBUG
246    std::cout << " -- case 2 --" << std::endl;
247#endif
[458]248
249    defcons_size = 12;
250
251    prepareVectors (defcons_size);
252
253    // compute the 6 permutations of the 3 variables
254    ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
255    permutation3(ind,ibnd);
256    int i, flagg=0, idx=0;
257    i = 0;
258    while(i < 6 && flagg == 0) {
259      if(vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
260         <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
261        {
262          idx = i;   // store the index of the permutation satisfying the condition
263          flagg = 1;  // condition is satisfied
264        }
265      i++;
266    }
267    if (flagg==0) {
268      std::cout << "ERROR!!!" << std::endl; exit(0);
269    }
270    v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
271
272    double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
273    double xL2(vlb[v2]); double xU2(vub[v2]);
274    double xL3(vlb[v3]); double xU3(vub[v3]);
275
[469]276    for(int ii = 0; ii < defcons_size; ii++) {
277
278      cutIndices [ii][0] = v1;
279      cutIndices [ii][1] = v2;
280      cutIndices [ii][2] = v3;
281      cutIndices [ii][3] = v4;
282      cutCoeff [ii][3] = 1.;
283    }
284
[458]285    double theta1 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
286    double theta2 = xU1*xL2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xU2*xU3;
287
288    cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xU1*xU3; cutCoeff [0][2] = -xU1*xU2;  bnd[0] = - 2.*xU1*xU2*xU3;
289    cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
290    cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xU2;  bnd[2] = - xL1*xU2*xL3 - xL1*xU2*xU3;
291    cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xL1*xU2*xL3 - xL1*xL2*xL3;
292    cutCoeff [4][0] = -xL2*xU3; cutCoeff [4][1] = -(theta1/(xL2-xU2)); cutCoeff [4][2] = -xL1*xL2;
293    bnd[4] = (-(theta1*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
294    cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -(theta2/(xU3-xL3));
295    bnd[5] = (-(theta2*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xU1*xU2*xL3;
296
297    if (  vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
298          >= vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]) {
299
300      double theta1c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
301      double theta2c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xL2*xU3;
302
303      cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xU1*xU3; cutCoeff [6][2] = -xU1*xL2;  bnd[6] = - 2.*xU1*xL2*xU3;
304      cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
305      cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xL1*xL2;  bnd[8] = - xL1*xU2*xU3 - xL1*xL2*xU3;
306      cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
307      cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -xL1*xL3; cutCoeff [10][2] = -(theta1c/(xL3-xU3));
308      bnd[10] = (-(theta1c*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
309      cutCoeff [11][0] = -xL2*xL3; cutCoeff [11][1] = -(theta2c/(xL2-xU2)); cutCoeff [11][2] = -xL1*xL2;
310      bnd[11] = (-(theta2c*xU2)/(xL2-xU2)) - xU1*xL2*xL3 - xL1*xL2*xU3 + xU1*xU2*xU3;
311
312    } else {
313      double theta1c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
314      double theta2c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
315
316      cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xU1*xU3; cutCoeff [6][2] = -xU1*xL2;  bnd[6] = - 2.*xU1*xL2*xU3;
317      cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
318      cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2;  bnd[8] = - xL1*xL2*xL3 - xL1*xU2*xL3;
319      cutCoeff [9][0] = -xL2*xL3; cutCoeff [9][1] = -xL1*xU3; cutCoeff [9][2] = -xL1*xL2;  bnd[9] = - xL1*xL2*xL3 - xL1*xL2*xU3;
320      cutCoeff [10][0] = -xU2*xU3; cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -(theta1c/(xU3-xL3));
321      bnd[10] = (-(theta1c*xL3)/(xU3-xL3)) - xU1*xU2*xU3 - xL1*xL2*xU3 + xU1*xL2*xL3;
322      cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
323      bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
324    }
325
326  } // end if case 2
327
328  /*----------------------------------------------------------------------------------------*/
329
330  // case 3
331  if(flag == 3) {
[469]332#ifdef DEBUG
333    std::cout << " -- case 3 --" << std::endl;
334#endif
[458]335
336    int last;
337
338    if(vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3]
339       <= vlb[v1]*vlb[v2]*vlb[v3] + 2.*vub[v1]*vub[v2]*vub[v3] &&
340       vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]
341       <= vlb[v1]*vub[v2]*vub[v3] + 2.*vub[v1]*vlb[v2]*vlb[v3] &&
342       vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
343       <= vub[v1]*vlb[v2]*vub[v3] + 2.*vlb[v1]*vub[v2]*vlb[v3] &&
344       vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3]
345       <= vub[v1]*vub[v2]*vlb[v3] + 2.*vlb[v1]*vlb[v2]*vub[v3] ) {
346
347      double theta3x = 0.5*(xL1*xU2*xU3 + xL1*xL2*xL3 - xU1*xU2*xL3 - xU1*xL2*xU3)/(xL1-xU1);
348      double theta3y = 0.5*(xU1*xL2*xU3 + xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xU2*xU3)/(xL2-xU2);
349      double theta3z = 0.5*(xU1*xU2*xL3 + xL1*xL2*xL3 - xU1*xL2*xU3 - xL1*xU2*xU3)/(xL3-xU3);
350      double theta3c = xL1*xL2*xL3 - theta3x*xL1 - theta3y*xL2 - theta3z*xL3;
351
352      defcons_size = 5;
353      prepareVectors (defcons_size);
354
355      for(int ii = 0; ii < 5; ii++) {
356
357        cutIndices [ii][0] = v1;
358        cutIndices [ii][1] = v2;
359        cutIndices [ii][2] = v3;
360        cutIndices [ii][3] = v4;
361
362        cutCoeff [ii][3] = 1.;
363      }
364
365      cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
366      cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2;  bnd[1] = - 2.*xL1*xL2*xU3;
367      cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - 2.*xU1*xU2*xU3;
368      cutCoeff [3][0] = -xL2*xL3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xU1*xL2;  bnd[3] = - 2.*xU1*xL2*xL3;
369      cutCoeff [4][0] = -theta3x; cutCoeff [4][1] = -theta3y; cutCoeff [4][2] = -theta3z;  bnd[4] = theta3c;
370
371      last=4;
372
373    } else if (vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] + vlb[v1]*vub[v2]*vub[v3]
374               >= vlb[v1]*vlb[v2]*vlb[v3] + 2.*vub[v1]*vub[v2]*vub[v3]) {
[469]375#ifdef DEBUG
376std::cout << "else if " << std::endl;
377#endif
[458]378
379      double theta1 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
380      double theta2 = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
381      double theta3 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
382
383      defcons_size = 6;
384      prepareVectors (defcons_size);
385
386      for(int ii = 0; ii < 6; ii++) {
387
388        cutIndices [ii][0] = v1;
389        cutIndices [ii][1] = v2;
390        cutIndices [ii][2] = v3;
391        cutIndices [ii][3] = v4;
392
393        cutCoeff [ii][3] = 1.;
394      }
395
396      cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
397      cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2;  bnd[1] = - 2.*xL1*xL2*xU3;
398      cutCoeff [2][0] = -xL2*xL3; cutCoeff [2][1] = -xU1*xL3; cutCoeff [2][2] = -xU1*xL2;  bnd[2] = - 2.*xU1*xL2*xL3;
399      cutCoeff [3][0] = -(theta1/(xU1-xL1)); cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2;
400      bnd[3] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
401      cutCoeff [4][0] = -xU2*xU3; cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -(theta2/(xU3-xL3));
402      bnd[4] = (-(theta2*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xL1*xL2*xL3;
403      cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -(theta3/(xU2-xL2)); cutCoeff [5][2] = -xU1*xU2;
404      bnd[5] = (-(theta3*xL2)/(xU2-xL2)) - xU1*xU2*xL3 - xL1*xU2*xU3 + xL1*xL2*xL3;
405
406      last=5;
407
408    } else {
[469]409#ifdef DEBUG
410std::cout << "else " << std::endl;
411#endif
[458]412      // compute the 6 permutations of the 3 variables
413      ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
414      permutation3(ind,ibnd);
415      int i, flagg=0, idx=0;
416      i = 0;
417      while(i < 6 && flagg == 0) {
418        if (((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
[469]419             >= vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] + 2.*vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]]) ))
[458]420          {
421            idx = i;   // store the index of the permutation satisfying the condition
422            flagg = 1;  // condition is satisfied
423          }
424        i++;
425      }
426      if (flagg==0) {
427        std::cout << "ERROR!!!" << std::endl; exit(0);
428      }
429      v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
430
431      double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
432      double xL2(vlb[v2]); double xU2(vub[v2]);
433      double xL3(vlb[v3]); double xU3(vub[v3]);
434
435      //} else if (vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]
436      //         >= vlb[v1]*vub[v2]*vub[v3] + 2.*vub[v1]*vlb[v2]*vlb[v3]) {
437
438      defcons_size = 6;
439      prepareVectors (defcons_size);
440
441      for(int ii = 0; ii < 6; ii++) {
442
443        cutIndices [ii][0] = v1;
444        cutIndices [ii][1] = v2;
445        cutIndices [ii][2] = v3;
446        cutIndices [ii][3] = v4;
447
448        cutCoeff [ii][3] = 1.;
449      }
450
451      double theta1 = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
452      double theta2 = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
453      double theta3 = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xU2*xL3;
454
455      cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
456      cutCoeff [1][0] = -xL2*xU3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xL2;  bnd[1] = - 2.*xL1*xL2*xU3;
457      cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - 2.*xU1*xU2*xU3;
458      cutCoeff [3][0] = -xL2*xL3; cutCoeff [3][1] = -(theta1/(xL2-xU2)); cutCoeff [3][2] = -xU1*xL2;
459      bnd[3] = (-(theta1*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
460      cutCoeff [4][0] = -(theta2/(xU1-xL1)); cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xU1*xL2;
461      bnd[4] = (-(theta2*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
462      cutCoeff [5][0] = -xL2*xL3; cutCoeff [5][1] = -xU1*xL3; cutCoeff [5][2] = -(theta3/(xL3-xU3));
463      bnd[5] = (-(theta3*xU3)/(xL3-xU3)) - xL1*xL2*xL3 - xU1*xU2*xL3 + xL1*xU2*xU3;
464
465      last=5;
466    }
467
468    if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
469       >= vub[v1]*vub[v2]*vub[v3] + 2.*vlb[v1]*vlb[v2]*vlb[v3] &&
470       vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] + vub[v1]*vub[v2]*vub[v3]
471       >= vub[v1]*vlb[v2]*vlb[v3] + 2.*vlb[v1]*vub[v2]*vub[v3] &&
472       vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] + vub[v1]*vub[v2]*vub[v3]
473       >= vlb[v1]*vub[v2]*vlb[v3] + 2.*vub[v1]*vlb[v2]*vub[v3] &&
474       vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
475       >= vlb[v1]*vlb[v2]*vub[v3] + 2.*vub[v1]*vub[v2]*vlb[v3] ) {
[469]476
477#ifdef DEBUG
478std::cout << "2 - if " << std::endl;
479#endif
[458]480      double theta3x = 0.5*(xU1*xL2*xL3 + xU1*xU2*xU3 - xL1*xL2*xU3 - xL1*xU2*xL3)/(xU1-xL1);
481      double theta3y = 0.5*(xL1*xU2*xL3 + xU1*xU2*xU3 - xL1*xL2*xU3 - xU1*xL2*xL3)/(xU2-xL2);
482      double theta3z = 0.5*(xL1*xL2*xU3 + xU1*xU2*xU3 - xL1*xU2*xL3 - xU1*xL2*xL3)/(xU3-xL3);
483      double theta3c = xU1*xU2*xU3 - theta3x*xU1 - theta3y*xU2 - theta3z*xU3;
484
485      prepareVectors (5);
486
487      for(int ii = last+1; ii <= last+5; ii++) {
[469]488
489        cutIndices [ii][0] = v1;
490        cutIndices [ii][1] = v2;
491        cutIndices [ii][2] = v3;
492        cutIndices [ii][3] = v4;
[458]493        cutCoeff [ii][3] = 1.;
494      }
495
496      cutCoeff [last+1][0] = -xL2*xL3; cutCoeff [last+1][1] = -xL1*xL3; cutCoeff [last+1][2] = -xL1*xL2;  bnd[last+1] = - 2.*xL1*xL2*xL3;
497      cutCoeff [last+2][0] = -xU2*xL3; cutCoeff [last+2][1] = -xU1*xL3; cutCoeff [last+2][2] = -xU1*xU2;  bnd[last+2] = - 2.*xU1*xU2*xL3;
498      cutCoeff [last+3][0] = -xL2*xU3; cutCoeff [last+3][1] = -xU1*xU3; cutCoeff [last+3][2] = -xU1*xL2;  bnd[last+3] = - 2.*xU1*xL2*xU3;
499      cutCoeff [last+4][0] = -xU2*xU3; cutCoeff [last+4][1] = -xL1*xU3; cutCoeff [last+4][2] = -xL1*xU2;  bnd[last+4] = - 2.*xL1*xU2*xU3;
500      cutCoeff [last+5][0] = -theta3x; cutCoeff [last+5][1] = -theta3y; cutCoeff [last+5][2] = -theta3z;  bnd[last+5] = theta3c;
501
502      defcons_size = last+6;
503
504    } else if (vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
505               <= vub[v1]*vub[v2]*vub[v3] + 2.*vlb[v1]*vlb[v2]*vlb[v3]) {
[469]506#ifdef DEBUG
507std::cout << "2 - else if" << std::endl;
508#endif
[458]509
510      double theta1 = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
511      double theta2 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
512      double theta3 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xU1*xL2*xL3;
513
514      prepareVectors (6);
515
516      for(int ii = last+1; ii <= last+6; ii++) {
517
518        cutIndices [ii][0] = v1;
519        cutIndices [ii][1] = v2;
520        cutIndices [ii][2] = v3;
521        cutIndices [ii][3] = v4;
522
523        cutCoeff [ii][3] = 1.;
524      }
525
526      cutCoeff [last+1][0] = -xU2*xL3; cutCoeff [last+1][1] = -xU1*xL3; cutCoeff [last+1][2] = -xU1*xU2;  bnd[last+1] = - 2.*xU1*xU2*xL3;
527      cutCoeff [last+2][0] = -xL2*xU3; cutCoeff [last+2][1] = -xU1*xU3; cutCoeff [last+2][2] = -xU1*xL2;  bnd[last+2] = - 2.*xU1*xL2*xU3;
528      cutCoeff [last+3][0] = -xU2*xU3; cutCoeff [last+3][1] = -xL1*xU3; cutCoeff [last+3][2] = -xL1*xU2;  bnd[last+3] = - 2.*xL1*xU2*xU3;
529      cutCoeff [last+4][0] = -xL2*xL3; cutCoeff [last+4][1] = -xL1*xL3; cutCoeff [last+4][2] = -(theta1/(xL3-xU3));
530      bnd[last+4] = (-(theta1*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
531      cutCoeff [last+5][0] = -(theta2/(xL1-xU1)); cutCoeff [last+5][1] = -xL1*xL3; cutCoeff [last+5][2] = -xL1*xL2;
532      bnd[last+5] = (-(theta2*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xU2*xU3;
533      cutCoeff [last+6][0] = - xL2*xL3; cutCoeff [last+6][1] = -(theta3/(xL2-xU2)); cutCoeff [last+6][2] = -xL1*xL2;
534      bnd[last+6] = (-(theta3*xU2)/(xL2-xU2)) - xL1*xL2*xU3 - xU1*xL2*xL3 + xU1*xU2*xU3;
535
536      defcons_size = last+7;
537
[469]538    } else //if (vlb[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] + vub[v1]*vub[v2]*vub[v3]
539           //    <= vub[v1]*vlb[v2]*vlb[v3] + 2.*vlb[v1]*vub[v2]*vub[v3])
540      {
541#ifdef DEBUG
542std::cout << "2 - another else if" << std::endl;
543std::cout << "v1 = " << v1 << " v2 =" << v2 << "  v3 =" << v3 << std::endl;
544#endif
545      // compute the 6 permutations of the 3 variables
546      //ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
547      permutation3(ind,ibnd);
548      int i, flagg=0, idx=0;
549      i = 0;
550      while(i < 6 && flagg == 0) {
551        if (vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
552               <= vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + 2.*vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]])
553          {
554            idx = i;   // store the index of the permutation satisfying the condition
555            flagg = 1;  // condition is satisfied
556          }
557        i++;
558      }
559      if (flagg==0) {
560        std::cout << "ERROR!!!" << std::endl; exit(0);
561      }
562      v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
[458]563
[469]564      double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
565      double xL2(vlb[v2]); double xU2(vub[v2]);
566      double xL3(vlb[v3]); double xU3(vub[v3]);
567
[458]568      double theta1 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
569      double theta2 = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
570      double theta3 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
571
572      prepareVectors (6);
573
574      for(int ii = last+1; ii <= last+6; ii++) {
575
576        cutIndices [ii][0] = v1;
577        cutIndices [ii][1] = v2;
578        cutIndices [ii][2] = v3;
579        cutIndices [ii][3] = v4;
580
581        cutCoeff [ii][3] = 1.;
582      }
583
584      cutCoeff [last+1][0] = -xL2*xL3; cutCoeff [last+1][1] = -xL1*xL3; cutCoeff [last+1][2] = -xL1*xL2;  bnd[last+1] = - 2.*xL1*xL2*xL3;
585      cutCoeff [last+2][0] = -xU2*xL3; cutCoeff [last+2][1] = -xU1*xL3; cutCoeff [last+2][2] = -xU1*xU2;  bnd[last+2] = - 2.*xU1*xU2*xL3;
586      cutCoeff [last+3][0] = -xL2*xU3; cutCoeff [last+3][1] = -xU1*xU3; cutCoeff [last+3][2] = -xU1*xL2;  bnd[last+3] = - 2.*xU1*xL2*xU3;
587      cutCoeff [last+4][0] = -(theta1/(xL1-xU1)); cutCoeff [last+4][1] = -xL1*xU3; cutCoeff [last+4][2] = -xL1*xU2;
[469]588      bnd[last+4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
[458]589      cutCoeff [last+5][0] = -xU2*xU3; cutCoeff [last+5][1] = -(theta2/(xU2-xL2)); cutCoeff [last+5][2] = -xL1*xU2;
590      bnd[last+5] = (-(theta2*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
591      cutCoeff [last+6][0] = -xU2*xU3; cutCoeff [last+6][1] = -xL1*xU3; cutCoeff [last+6][2] = -(theta3/(xU3-xL3));
592      bnd[last+6] = (-(theta3*xL3)/(xU3-xL3)) - xL1*xL2*xU3 - xU1*xU2*xU3 + xU1*xL2*xL3;
593
594      defcons_size = last+7;
595    }
596
597  } // end if case 3
598
599  /*----------------------------------------------------------------------------------------*/
600
601  // case 4
602  if(flag == 4) {
[469]603#ifdef DEBUG
604    std::cout << " -- case 4 --" << std::endl;
605#endif
[458]606
607    double theta  = xU1*xL2*xU3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xL2*xL3;
608    double theta1 = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
609
610    defcons_size = 12;
611
612    prepareVectors (defcons_size);
613
614    for(int ii = 0; ii < defcons_size; ii++) {
615      cutIndices [ii][0] = v1;
616      cutIndices [ii][1] = v2;
617      cutIndices [ii][2] = v3;
618      cutIndices [ii][3] = v4;
619      cutCoeff [ii][3] = 1.;
620    }
621
622    cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xL2;  bnd[0] = - 2.*xU1*xL2*xL3;
623    cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3;
624    cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xL1*xU2*xL3 - xL1*xL2*xL3;
625    cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2;  bnd[3] = - xU1*xL2*xU3 - xU1*xU2*xU3;
626    cutCoeff [4][0] = -xU2*xU3; cutCoeff [4][1] = -xL1*xU3; cutCoeff [4][2] = -xU1*xU2;  bnd[4] = - xL1*xU2*xU3 - xU1*xU2*xU3;
627    cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
628    bnd[5] = (-(theta*xU2)/(xL2-xU2)) - xU1*xL2*xU3 - xL1*xL2*xL3 + xU1*xU2*xL3;
629
630    cutCoeff [6][0] = -xU2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2;  bnd[6] = - 2.*xU1*xU2*xL3;
631    cutCoeff [7][0] = -xU2*xU3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xL2;  bnd[7] = - xU1*xU2*xU3 - xU1*xL2*xU3;
632    cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2;  bnd[8] = - xL1*xU2*xL3 - xL1*xL2*xL3;
633    cutCoeff [9][0] = -xL2*xL3; cutCoeff [9][1] = -xL1*xU3; cutCoeff [9][2] = -xL1*xL2;  bnd[9] = - xL1*xL2*xU3 - xL1*xL2*xL3;
634    cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xU1*xL2;  bnd[10] = - xL1*xL2*xU3 - xU1*xL2*xU3;
635    cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta1/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
636    bnd[11] = (-(theta1*xL2)/(xU2-xL2)) - xL1*xU2*xL3 - xU1*xU2*xU3 + xU1*xL2*xL3;
637
638  } // end if case 4
639
640  /*----------------------------------------------------------------------------------------*/
641
642  // case 5
643  if(flag == 5) {
[469]644#ifdef DEBUG
645    std::cout << " -- case 5 --" << std::endl;
646 std::cout << "v1 = " << v1 << " v2 =" << v2 << "  v3 =" << v3 << std::endl;
647#endif
[458]648
649    defcons_size = 12;
650    prepareVectors (defcons_size);
651
[469]652    // compute the permutations of the 3 variables
[458]653    ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
[469]654   ind[0][0] = ibnd[0]; ind[0][1] = ibnd[1]; ind[0][2] = ibnd[2];
655   ind[1][0] = ibnd[1]; ind[1][1] = ibnd[0]; ind[1][2] = ibnd[2];
[458]656    int i, flagg=0, idx=0;
657    i = 0;
[469]658    while(i < 2 && flagg == 0) {
[458]659      if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
660         >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
661        {
662          idx = i;   // store the index of the permutation satisfying the condition
663          flagg = 1;  // condition is satisfied
664        }
665      i++;
666    }
667    if (flagg==0) {
668      std::cout << "ERROR!!!" << std::endl; exit(0);
669    }
670    v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
[469]671#ifdef DEBUG
672 std::cout << "v1 = " << v1 << " v2 =" << v2 << "  v3 =" << v3 << std::endl;
673#endif
[458]674
[469]675    double xL1 = cf*vlb[v1]; double xU1 = cf*vub[v1];
676    double xL2 = vlb[v2]; double xU2 = vub[v2];
677    double xL3 = vlb[v3]; double xU3 = vub[v3];
[458]678
[469]679    for(int ii = 0; ii < defcons_size; ii++) {
680
681      cutIndices [ii][0] = v1;
682      cutIndices [ii][1] = v2;
683      cutIndices [ii][2] = v3;
684      cutIndices [ii][3] = v4;
685      cutCoeff [ii][3] = 1.;
686    }
687
[458]688    if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
689       <= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]) {
[469]690#ifdef DEBUG
691    std::cout << " -- 5 if --" << std::endl;
692#endif
[458]693
694      double theta1 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
695      double theta2 = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
696
697      cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
698      cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
699      cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xL1*xU2*xU3 - xL1*xL2*xU3;
700      cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xU1*xL2*xU3 - xL1*xL2*xU3;
701      cutCoeff [4][0] = -(theta1/(xU1-xL1)); cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -xU1*xU2;
702      bnd[4] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
703      cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -(theta2/(xU2-xL2)); cutCoeff [5][2] = -xU1*xU2;
704      bnd[5] = (-(theta2*xL2)/(xU2-xL2)) - xU1*xU2*xL3 - xL1*xU2*xU3 + xL1*xL2*xL3;
705
706    } else {
[469]707#ifdef DEBUG
708    std::cout << " -- 5 else --" << std::endl;
709#endif
[458]710      double theta1 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xL1*xU2*xU3;
711      double theta2 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
712
713      cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xL3;
714      cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
715      cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - xU1*xL2*xU3 - xU1*xU2*xU3;
716      cutCoeff [3][0] = -xU2*xU3; cutCoeff [3][1] = -xL1*xU3; cutCoeff [3][2] = -xU1*xU2;  bnd[3] = - xL1*xU2*xU3 - xU1*xU2*xU3;
717      cutCoeff [4][0] = -(theta1/(xL1-xU1)); cutCoeff [4][1] = -xL1*xU3; cutCoeff [4][2] = -xL1*xL2;
718      bnd[4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xL3;
719      cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta2/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
720      bnd[5] = (-(theta2*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
721    }
722
723    double theta1c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
724    double theta2c = xU1*xU2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
725
726    cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xL3;
727    cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
728    cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
729    cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2;  bnd[9] = - xU1*xU2*xU3 - xU1*xL2*xU3;
730    cutCoeff [10][0] = -(theta1c/(xL1-xU1)); cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xL1*xU2;
731    bnd[10] = (-(theta1c*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
732    cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
733    bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
734
735  } // end if case 5
736
737  /*----------------------------------------------------------------------------------------*/
738
739  // case 6
740  if(flag == 6) {
[469]741#ifdef DEBUG
742    std::cout << " -- case 6 --" << std::endl;
743#endif
[458]744
745    double theta = xU1*xU2*xL3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xL2*xU3;
746    double theta1 = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
747
748    defcons_size = 12;
749    prepareVectors (defcons_size);
750
751    for (int ii = 0; ii < defcons_size; ii++) {
752
753      cutIndices [ii][0] = v1;
754      cutIndices [ii][1] = v2;
755      cutIndices [ii][2] = v3;
756      cutIndices [ii][3] = v4;
757
758      cutCoeff [ii][3] = 1.;
759    }
760
761    cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xL2;  bnd[0] = - 2.*xU1*xL2*xL3;
762    cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xL1*xL3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xU3 - xL1*xU2*xL3;
763    cutCoeff [2][0] = -xU2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xL1*xU2*xU3 - xL1*xL2*xU3;
764    cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xU1*xL2*xU3 - xL1*xL2*xU3;
765    cutCoeff [4][0] = -xU2*xL3; cutCoeff [4][1] = -xL1*xL3; cutCoeff [4][2] = -xU1*xU2;  bnd[4] = - xL1*xU2*xL3 - xU1*xU2*xL3;
766    cutCoeff [5][0] = -(theta/(xU1-xL1)); cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xU2;
767    bnd[5] = (-(theta*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xL2*xL3;
768
769    cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xL3;
770    cutCoeff [7][0] = -xL2*xU3; cutCoeff [7][1] = -xL1*xU3; cutCoeff [7][2] = -xU1*xL2;  bnd[7] = - xL1*xL2*xU3 - xU1*xL2*xU3;
771    cutCoeff [8][0] = -xU2*xU3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xU1*xU2*xU3 - xU1*xL2*xU3;
772    cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xL3; cutCoeff [9][2] = -xU1*xU2;  bnd[9] = - xU1*xU2*xU3 - xU1*xU2*xL3;
773    cutCoeff [10][0] = -xU2*xL3; cutCoeff [10][1] = -xU1*xL3; cutCoeff [10][2] = -xL1*xU2;  bnd[10] = - xL1*xU2*xL3 - xU1*xU2*xL3;
774    cutCoeff [11][0] = -(theta1/(xL1-xU1)); cutCoeff [11][1] = -xL1*xU3; cutCoeff [11][2] = -xL1*xU2;
775    bnd[11] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
776  } // end if case 6
777
778  /*----------------------------------------------------------------------------------------*/
779
780  // case 7
781  if(flag == 7) {
[469]782#ifdef DEBUG
783    std::cout << " -- case 7 --" << std::endl;
784#endif
[458]785
[469]786    defcons_size = 12;
[458]787    prepareVectors (defcons_size);
788
[471]789    if ((vlb[v1]<=EPSILONT && vlb[v2]<=EPSILONT && vlb[v3]<=EPSILONT) ||
790        (vlb[v1]==vlb[v2] && vlb[v1]==vlb[v3] && vub[v1]==vub[v2] && vub[v1]==vub[v3])) {
[470]791      int idx = 0;
792#ifdef DEBUG
[471]793      std::cout << " -- epsilonT --" << std::endl;
[470]794#endif
795    } else {
[458]796    // compute the 6 permutations of the 3 variables
797    ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
798    permutation3(ind,ibnd);
799    int i, flagg=0, idx=0;
800    i = 0;
801    while(i < 6 && flagg == 0) {
802      if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
803         <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
804         vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
805         <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
806        {
807          idx = i;   // store the index of the permutation satisfying the condition
808          flagg = 1;  // condition is satisfied
809        }
810      i++;
811    }
812    if (flagg==0) {
813      std::cout << "ERROR!!!" << std::endl; exit(0);
814    }
[470]815    }
[458]816    v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
817
[469]818    double xL1 = cf*vlb[v1]; double xU1 = cf*vub[v1];
819    double xL2 = vlb[v2]; double xU2 = vub[v2];
820    double xL3 = vlb[v3]; double xU3 = vub[v3];
[458]821
822    //if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
823    // <= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] &&
824    // vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
825    // <= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]) {
826
827    double theta1 = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
828    double theta2 = xL1*xL2*xU3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xU2*xL3;
829
[469]830    for(int ii = 0; ii < defcons_size; ii++) {
831      cutIndices [ii][0] = v1;
832      cutIndices [ii][1] = v2;
833      cutIndices [ii][2] = v3;
834      cutIndices [ii][3] = v4;
835
836      cutCoeff [ii][3] = 1.;
837    }
838
[458]839    cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xL2;  bnd[0] = - 2.*xL1*xL2*xL3;
840    cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xU1*xU3; cutCoeff [1][2] = -xU1*xU2;  bnd[1] = - 2.*xU1*xU2*xU3;
841    cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xL1*xU3; cutCoeff [2][2] = -xU1*xL2;  bnd[2] = - xL1*xL2*xU3 - xU1*xL2*xU3;
842    cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xU1*xL3; cutCoeff [3][2] = -xL1*xU2;  bnd[3] = - xU1*xU2*xL3 - xL1*xU2*xL3;
843    cutCoeff [4][0] = -(theta1/(xU1-xL1)); cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xU1*xL2;
844    bnd[4] = (-(theta1*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
845    cutCoeff [5][0] = -(theta2/(xL1-xU1)); cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -xL1*xU2;
846    bnd[5] = (-(theta2*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
847    //}
848    cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xU1*xL3; cutCoeff [6][2] = -xU1*xU2;  bnd[6] = - xU1*xU2*xL3 - xU1*xL2*xL3;
849    cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xL1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - xU1*xU2*xL3 - xL1*xU2*xL3;
850    cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xU1*xL2*xU3 - xU1*xL2*xL3;
851    cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
852    cutCoeff [10][0] = -xL2*xU3; cutCoeff [10][1] = -xU1*xU3; cutCoeff [10][2] = -xL1*xL2;  bnd[10] = - xU1*xL2*xU3 - xL1*xL2*xU3;
853    cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -xL1*xU3; cutCoeff [11][2] = -xL1*xL2;  bnd[11] = - xL1*xU2*xU3 - xL1*xL2*xU3;
854  } // end if case 7
855
856  /*----------------------------------------------------------------------------------------*/
857
858  // case 8
859  if(flag == 8) {
[469]860#ifdef DEBUG
861    std::cout << " -- case 8 --" << std::endl;
862#endif
[458]863
864    defcons_size = 12;
865    prepareVectors (defcons_size);
866
867    // compute the 6 permutations of the 3 variables
868    ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
869    permutation3(ind,ibnd);
870    int i, flagg=0, idx=0;
871    i = 0;
872    while(i < 6 && flagg == 0) {
873      if(vub[ind[i][2]] <=0  &&
874         ((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
875           >= vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] &&
876           vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
877           >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) ||
878          (vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
879           >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])) )
880        {
881          idx = i;   // store the index of the permutation satisfying the condition
882          flagg = 1;  // condition is satisfied
883        }
884      i++;
885    }
886    if (flagg==0) {
887      std::cout << "ERROR!!!" << std::endl; exit(0);
888    }
889    v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
890
891    double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
892    double xL2(vlb[v2]); double xU2(vub[v2]);
893    double xL3(vlb[v3]); double xU3(vub[v3]);
894
[469]895    for(int ii = 0; ii < defcons_size; ii++) {
896
897      cutIndices [ii][0] = v1;
898      cutIndices [ii][1] = v2;
899      cutIndices [ii][2] = v3;
900      cutIndices [ii][3] = v4;
901      cutCoeff [ii][3] = 1.;
902    }
903
904    // compute the 6 permutations of the 3 variables
[458]905    //if(vub[v3]<=0) {
906    cutCoeff [0][0] = -xU2*xL3; cutCoeff [0][1] = -xL1*xL3; cutCoeff [0][2] = -xL1*xL2;  bnd[0] = - xL1*xU2*xL3 - xL1*xL2*xL3;
907    cutCoeff [1][0] = -xU2*xL3; cutCoeff [1][1] = -xL1*xU3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xL3 - xL1*xU2*xU3;
908    cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xL3; cutCoeff [2][2] = -xU1*xL2;  bnd[2] = - xU1*xL2*xU3 - xU1*xL2*xL3;
909    cutCoeff [3][0] = -xL2*xU3; cutCoeff [3][1] = -xU1*xU3; cutCoeff [3][2] = -xU1*xU2;  bnd[3] = - xU1*xL2*xU3 - xU1*xU2*xU3;
910    cutCoeff [4][0] = -xL2*xL3; cutCoeff [4][1] = -xU1*xL3; cutCoeff [4][2] = -xL1*xL2;  bnd[4] = - xU1*xL2*xL3 - xL1*xL2*xL3;
911    cutCoeff [5][0] = -xU2*xU3; cutCoeff [5][1] = -xL1*xU3; cutCoeff [5][2] = -xU1*xU2;  bnd[5] = - xU1*xU2*xU3 - xL1*xU2*xU3;
912
913    if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
914       >= vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3] &&
915       vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
916       >= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
917
918      double theta1c = xU1*xL2*xL3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
919      double theta2c = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xL1*xU2*xU3;
920
921      cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xU3;
922      cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
923      cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xU1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xU1*xL2*xU3 - xU1*xL2*xL3;
924      cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xL1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xU3 - xL1*xU2*xL3;
925      cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -xL1*xL3; cutCoeff [10][2] = -(theta1c/(xL3-xU3));
926      bnd[10] = (-(theta1c*xU3)/(xL3-xU3)) - xU1*xL2*xL3 - xL1*xU2*xL3 + xU1*xU2*xU3;
927      cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -xU1*xU3; cutCoeff [11][2] = -(theta2c/(xU3-xL3));
928      bnd[11] = (-(theta2c*xL3)/(xU3-xL3)) - xU1*xL2*xU3 - xL1*xU2*xU3 + xL1*xL2*xL3;
929
930    } else if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
931              >= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
932
933      double theta1c = xL1*xL2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
934      double theta2c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xU1*xU2*xU3;
935
936      cutCoeff [6][0] = -xL2*xU3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xU3;
937      cutCoeff [7][0] = -xU2*xL3; cutCoeff [7][1] = -xU1*xL3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xL3;
938      cutCoeff [8][0] = -xL2*xL3; cutCoeff [8][1] = -xL1*xL3; cutCoeff [8][2] = -xL1*xU2;  bnd[8] = - xL1*xL2*xL3 - xL1*xU2*xL3;
939      cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2;  bnd[9] = - xU1*xL2*xU3 - xU1*xU2*xU3;
940      cutCoeff [10][0] = -xL2*xL3; cutCoeff [10][1] = -(theta1c/(xL2-xU2)); cutCoeff [10][2] = -xU1*xL2;
941      bnd[10] = (-(theta1c*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
942      cutCoeff [11][0] = -xU2*xU3; cutCoeff [11][1] = -(theta2c/(xU2-xL2)); cutCoeff [11][2] = -xL1*xU2;
943      bnd[11] = (-(theta2c*xL2)/(xU2-xL2)) - xL1*xU2*xL3 - xU1*xU2*xU3 + xU1*xL2*xL3;
944    }
945    //}
946  } // end if case 8
947
948  /*----------------------------------------------------------------------------------------*/
949
950  // case 9
951  if(flag == 9) {
[469]952#ifdef DEBUG
953    std::cout << " -- case 9 --" << std::endl;
954#endif
[458]955
956    defcons_size = 12;
957    prepareVectors (defcons_size);
958
959    // compute the 6 permutations of the 3 variables
960    ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
961    permutation3(ind,ibnd);
962    int i, flagg=0, idx=0;
963    i = 0;
964    while(i < 6 && flagg == 0) {
965      if(vlb[ind[i][0]] >=0  &&
966         ((vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
967           <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
968           vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
969           <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) ||
970          (vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
971           <= vlb[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]] &&
972           vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]
973           <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]]) ))
974        {
975          idx = i;   // store the index of the permutation satisfying the condition
976          flagg = 1;  // condition is satisfied
977        }
978      i++;
979    }
980    if (flagg==0) {
981      std::cout << "ERROR!!!" << std::endl; exit(0);
982    }
983    v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
984
985    double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
986    double xL2(vlb[v2]); double xU2(vub[v2]);
987    double xL3(vlb[v3]); double xU3(vub[v3]);
988
[469]989    for (int ii = 0; ii < defcons_size; ii++) {
990
991      cutIndices [ii][0] = v1;
992      cutIndices [ii][1] = v2;
993      cutIndices [ii][2] = v3;
994      cutIndices [ii][3] = v4;
995      cutCoeff [ii][3] = 1.;
996    }
997
[458]998    //if(vlb[v1]>=0) {
999    if(vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
1000       <= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3] &&
1001       vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3]
1002       <= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
1003
1004      double theta1 = xL1*xL2*xU3 - xU1*xU2*xU3 - xL1*xL2*xL3 + xL1*xU2*xL3;
1005      double theta2 = xU1*xL2*xU3 - xL1*xL2*xL3 - xU1*xU2*xU3 + xU1*xU2*xL3;
1006
1007      cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xL1*xU3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xU3;
1008      cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
1009      cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xL1*xL2;  bnd[2] = - xU1*xL2*xU3 - xL1*xL2*xU3;
1010      cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xU1*xU2;  bnd[3] = - xU1*xU2*xL3 - xL1*xU2*xL3;
1011      cutCoeff [4][0] = -(theta1/(xL1-xU1)); cutCoeff [4][1] = -xL1*xL3; cutCoeff [4][2] = -xL1*xL2;
1012      bnd[4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xU2*xU3;
1013      cutCoeff [5][0] = -(theta2/(xU1-xL1)); cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xU2;
1014      bnd[5] = (-(theta2*xL1)/(xU1-xL1)) - xU1*xL2*xU3 - xU1*xU2*xL3 + xL1*xL2*xL3;
1015
1016    } else if(vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
1017              <= vlb[v1]*vlb[v2]*vlb[v3] + vub[v1]*vub[v2]*vub[v3] &&
1018              vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]
1019              <= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3]) {
1020
1021      double theta1 = xU1*xU2*xU3 - xL1*xL2*xU3 - xU1*xU2*xL3 + xL1*xU2*xL3;
1022      double theta2 = xL1*xL2*xL3 - xU1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xU3;
1023
1024      cutCoeff [0][0] = -xU2*xU3; cutCoeff [0][1] = -xL1*xU3; cutCoeff [0][2] = -xL1*xU2;  bnd[0] = - 2.*xL1*xU2*xU3;
1025      cutCoeff [1][0] = -xL2*xL3; cutCoeff [1][1] = -xU1*xL3; cutCoeff [1][2] = -xU1*xL2;  bnd[1] = - 2.*xU1*xL2*xL3;
1026      cutCoeff [2][0] = -xL2*xU3; cutCoeff [2][1] = -xU1*xU3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - xU1*xL2*xU3 - xU1*xU2*xU3;
1027      cutCoeff [3][0] = -xU2*xL3; cutCoeff [3][1] = -xL1*xL3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xL1*xL2*xL3 - xL1*xU2*xL3;
1028      cutCoeff [4][0] = -xU2*xL3; cutCoeff [4][1] = -(theta1/(xU2-xL2)); cutCoeff [4][2] = -xU1*xU2;
1029      bnd[4] = (-(theta1*xL2)/(xU2-xL2)) - xU1*xU2*xU3 - xL1*xU2*xL3 + xL1*xL2*xU3;
1030      cutCoeff [5][0] = -xL2*xU3; cutCoeff [5][1] = -(theta2/(xL2-xU2)); cutCoeff [5][2] = -xL1*xL2;
1031      bnd[5] = (-(theta2*xU2)/(xL2-xU2)) - xL1*xL2*xL3 - xU1*xL2*xU3 + xU1*xU2*xL3;
1032    }
1033
1034    cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xU3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - xL1*xL2*xU3 - xL1*xL2*xL3;
1035    cutCoeff [7][0] = -xL2*xL3; cutCoeff [7][1] = -xL1*xL3; cutCoeff [7][2] = -xL1*xU2;  bnd[7] = - xL1*xU2*xL3 - xL1*xL2*xL3;
1036    cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
1037    cutCoeff [9][0] = -xU2*xU3; cutCoeff [9][1] = -xU1*xU3; cutCoeff [9][2] = -xU1*xL2;  bnd[9] = - xU1*xU2*xU3 - xU1*xL2*xU3;
1038    cutCoeff [10][0] = -xU2*xU3; cutCoeff [10][1] = -xU1*xL3; cutCoeff [10][2] = -xU1*xU2;  bnd[10] = - xU1*xU2*xU3 - xU1*xU2*xL3;
1039    cutCoeff [11][0] = -xU2*xL3; cutCoeff [11][1] = -xU1*xL3; cutCoeff [11][2] = -xL1*xU2;  bnd[11] = - xL1*xU2*xL3 - xU1*xU2*xL3;
1040    //}
1041  } // end if case 9
1042
1043  /*----------------------------------------------------------------------------------------*/
1044
1045  // case 10
1046  if(flag == 10) {
[469]1047#ifdef DEBUG
1048    std::cout << " -- case 10 --" << std::endl;
1049#endif
[458]1050
1051    defcons_size = 12;
1052    prepareVectors (defcons_size);
1053
1054    ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
1055    permutation3(ind,ibnd);
1056    int i, flagg=0, idx=0;
1057    i = 0;
1058    while(i < 6 && flagg == 0) {
1059      if(vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
1060         >= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
1061         vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
[469]1062         >= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
[458]1063        {
1064          idx = i;   // store the index of the permutation satisfying the condition
1065          flagg = 1;  // condition is satisfied
1066        }
[469]1067      i++;
1068    }
1069    if (flagg==0) {
1070       std::cout << "ERROR!!!" << std::endl; exit(0);
1071    }
1072    v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
[458]1073
1074      double xL1(cf*vlb[v1]); double xU1(cf*vub[v1]);
1075      double xL2(vlb[v2]); double xU2(vub[v2]);
1076      double xL3(vlb[v3]); double xU3(vub[v3]);
1077
[469]1078      for(int ii = 0; ii < defcons_size; ii++) {
1079
1080        cutIndices [ii][0] = v1;
1081        cutIndices [ii][1] = v2;
1082        cutIndices [ii][2] = v3;
1083        cutIndices [ii][3] = v4;
1084        cutCoeff [ii][3] = 1.;
1085      }
1086
1087    // compute the 6 permutations of the 3 variables
[458]1088      cutCoeff [0][0] = -xL2*xL3; cutCoeff [0][1] = -xU1*xL3; cutCoeff [0][2] = -xU1*xU2;  bnd[0] = - xU1*xU2*xL3 - xU1*xL2*xL3;
1089      cutCoeff [1][0] = -xU2*xU3; cutCoeff [1][1] = -xL1*xL3; cutCoeff [1][2] = -xL1*xU2;  bnd[1] = - xL1*xU2*xU3 - xL1*xU2*xL3;
1090      cutCoeff [2][0] = -xU2*xL3; cutCoeff [2][1] = -xL1*xL3; cutCoeff [2][2] = -xU1*xU2;  bnd[2] = - xU1*xU2*xL3 - xL1*xU2*xL3;
1091      cutCoeff [3][0] = -xU2*xU3; cutCoeff [3][1] = -xL1*xU3; cutCoeff [3][2] = -xL1*xL2;  bnd[3] = - xL1*xU2*xU3 - xL1*xL2*xU3;
1092      cutCoeff [4][0] = -xL2*xU3; cutCoeff [4][1] = -xU1*xU3; cutCoeff [4][2] = -xL1*xL2;  bnd[4] = - xU1*xL2*xU3 - xL1*xL2*xU3;
1093      cutCoeff [5][0] = -xL2*xL3; cutCoeff [5][1] = -xU1*xU3; cutCoeff [5][2] = -xU1*xL2;  bnd[5] = - xU1*xL2*xU3 - xU1*xL2*xL3;
1094
1095      //if(vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
1096      // >= vlb[v1]*vub[v2]*vlb[v3] + vub[v1]*vlb[v2]*vub[v3] &&
1097      // vub[v1]*vlb[v2]*vlb[v3] + vlb[v1]*vub[v2]*vub[v3]
1098      // >= vub[v1]*vub[v2]*vlb[v3] + vlb[v1]*vlb[v2]*vub[v3]) {
1099
1100      double theta1c = xL1*xU2*xL3 - xU1*xL2*xL3 - xL1*xU2*xU3 + xL1*xL2*xU3;
1101      double theta2c = xU1*xU2*xL3 - xL1*xU2*xU3 - xU1*xL2*xL3 + xU1*xL2*xU3;
1102
1103      cutCoeff [6][0] = -xL2*xL3; cutCoeff [6][1] = -xL1*xL3; cutCoeff [6][2] = -xL1*xL2;  bnd[6] = - 2.*xL1*xL2*xL3;
1104      cutCoeff [7][0] = -xU2*xU3; cutCoeff [7][1] = -xU1*xU3; cutCoeff [7][2] = -xU1*xU2;  bnd[7] = - 2.*xU1*xU2*xU3;
1105      cutCoeff [8][0] = -xL2*xU3; cutCoeff [8][1] = -xL1*xU3; cutCoeff [8][2] = -xU1*xL2;  bnd[8] = - xL1*xL2*xU3 - xU1*xL2*xU3;
1106      cutCoeff [9][0] = -xU2*xL3; cutCoeff [9][1] = -xU1*xL3; cutCoeff [9][2] = -xL1*xU2;  bnd[9] = - xL1*xU2*xL3 - xU1*xU2*xL3;
1107      cutCoeff [10][0] = -(theta1c/(xL1-xU1)); cutCoeff [10][1] = -xL1*xU3; cutCoeff [10][2] = -xL1*xU2;
1108      bnd[10] = (-(theta1c*xU1)/(xL1-xU1)) - xL1*xU2*xL3 - xL1*xL2*xU3 + xU1*xL2*xL3;
1109      cutCoeff [11][0] = -(theta2c/(xU1-xL1)); cutCoeff [11][1] = -xU1*xL3; cutCoeff [11][2] = -xU1*xL2;
1110      bnd[11] = (-(theta2c*xL1)/(xU1-xL1)) - xU1*xU2*xL3 - xU1*xL2*xU3 + xL1*xU2*xU3;
1111
1112  } // end if case 10
1113
1114  /*----------------------------------------------------------------------------------------*/
1115
1116  // get lower and upper bound on the constraints to be added to the problem
1117  for (int ii = 0; ii < defcons_size; ii++) {
1118
1119    //printf ("flag = %d\n", flag);
1120
1121    if (ii < defcons_size/2) {
1122
1123      if (cf > 0) {cutLb[ii] =  bnd[ii];           cutUb[ii] = COUENNE_INFINITY;}
1124      if (cf < 0) {cutLb[ii] = -COUENNE_INFINITY;  cutUb[ii] = bnd[ii];}
1125
1126    } else {
1127
1128      if (cf > 0) {cutLb[ii] = -COUENNE_INFINITY;  cutUb[ii] = bnd[ii];}
1129      if (cf < 0) {cutLb[ii] =  bnd[ii];           cutUb[ii] = COUENNE_INFINITY;}
1130    }
[469]1131#ifdef DEBUG
1132std::cout << ii << ") cutLb =" << cutLb[ii] << " " << "cutUb = " << cutUb[ii] << std::endl;
1133#endif
[458]1134  }
[462]1135
1136  for (int i=0; i<6; i++)
1137    delete [] ind [i];
1138
[469]1139  delete [] ibnd;
[462]1140  delete [] bnd;
1141  delete [] ind;
[458]1142}
1143
1144
1145// generate cuts for trilinear expressions
1146void exprTrilinear::generateCuts (expression *w,
1147                                  OsiCuts &cs, const CouenneCutGenerator *cg,
1148                                  t_chg_bounds *chg, int wind,
1149                                  CouNumber lbw, CouNumber ubw) {
1150
1151  expression **args = w -> Image () -> ArgList ();
1152
1153  int *varInd = new int [4];
1154
1155  for (int i=0; i<3; i++)
1156    varInd [i] = args [i] -> Index ();
1157
1158  varInd [3] = w -> Index ();
1159
1160  std::vector <std::vector <int> >    cutIndices;
1161  std::vector <std::vector <double> > cutCoeff;
1162  std::vector <double>                cutLb, cutUb;
1163
1164  TriLinCuts (cg -> Problem () -> Lb (),
1165              cg -> Problem () -> Ub (),
1166              varInd,
1167              cutIndices, cutCoeff,
1168              cutLb, cutUb);
1169
1170  // sanity check on returned vectors
1171  assert (cutIndices.size () == cutCoeff.size () &&
1172          cutIndices.size () == cutLb.size    () &&
1173          cutIndices.size () == cutUb.size    ());
1174
[469]1175  //printf ("trilinear cuts:\n");
[458]1176
[459]1177  for (int i = (int) cutIndices.size (); i--;) {
[458]1178
1179    int
1180       size = (int) cutIndices [i].size (),
1181      *ind  = new int [size];
1182
1183    double *coe = new double [size];
1184
1185    std::copy (cutIndices [i].begin (), cutIndices [i].end (), ind);
1186    std::copy (cutCoeff   [i].begin (), cutCoeff   [i].end (), coe);
1187
1188    OsiRowCut cut (cutLb [i], cutUb [i], 4, 4, ind, coe);
[469]1189    //cut.print ();
[458]1190
[462]1191    delete [] ind;
1192    delete [] coe;
1193
[459]1194    if (cg -> Problem () -> bestSol ()) {
1195
1196      // check validity of cuts by verifying they don't cut the
1197      // optimum in the node
1198
1199      double *sol = cg -> Problem () -> bestSol ();
1200      const double
1201        *lb = cg -> Problem () -> Lb (),
1202        *ub = cg -> Problem () -> Ub ();
1203
1204      int nVars = cg -> Problem () -> nVars ();
1205
1206      bool optIn = true;
1207
1208      for (int i=0; i< nVars; i++)
1209        if ((sol [i] < lb [i] - COUENNE_EPS) ||
1210            (sol [i] > ub [i] + COUENNE_EPS)) {
1211          optIn = false;
1212          break;
1213        }
1214
1215      if (optIn) {
1216
1217        for (unsigned int i=0; i<cutIndices.size (); i++) {
1218
1219          double chs = 0.;
1220
1221          for (unsigned int j=0; j<cutIndices[i].size(); j++)
1222            chs += cutCoeff [i] [j] * sol [cutIndices [i] [j]];
1223
1224          if ((chs < cutLb [i] - COUENNE_EPS) ||
1225              (chs > cutUb [i] + COUENNE_EPS)) {
1226
[461]1227            printf ("cut %d violates optimum:\n", i);
[459]1228
1229            if (cutLb [i] > -COUENNE_INFINITY) printf ("%g <= ", cutLb [i]);
[461]1230            for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g x%d ", cutCoeff [i] [j],       cutIndices [i] [j]);  printf ("\n = ");
1231            for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g *%g ", cutCoeff [i] [j],  sol [cutIndices [i] [j]]); printf ("\n = ");
1232            for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g ",     cutCoeff [i] [j] * sol [cutIndices [i] [j]]); printf (" = %g", chs);
[459]1233            if (cutUb [i] <  COUENNE_INFINITY) printf (" <= %g", cutUb [i]);
1234            printf ("\n");
1235
1236          }
1237        }
1238      }
1239    }
1240
[458]1241    cs.insert (cut);
1242  }
1243
1244  delete [] varInd;
1245}
Note: See TracBrowser for help on using the repository browser.