conv-exprTrilinear-gencuts.cpp
Go to the documentation of this file.
1 /* $Id: conv-exprTrilinear-gencuts.cpp 1055 2014-01-31 18:27:50Z fmargot $
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 
21 //#define DEBUG
22 using namespace Couenne;
23 
24 #define EPSILONT 1.e-6
25 
26 //typedef CouNumber double;
27 
29 void 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 
41 void 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 
62  int *ibnd;
63  ibnd = new int[3];
64  ibnd[0] = varIndices[0]; ibnd[1] = varIndices[1]; ibnd[2] = varIndices[2];
65 #ifdef DEBUG
66 std::cout << "ibnd[0] =" << ibnd[0] << " ibnd[1] =" << ibnd[1] << " ibnd[2] =" << ibnd[2] << std::endl;
67 std::cout << "vlb[ibnd[0]] =" << vlb[ibnd[0]] << " vub[ibnd[0]] =" << vub[ibnd[0]] << std::endl;
68 std::cout << "vlb[ibnd[1]] =" << vlb[ibnd[1]] << " vub[ibnd[1]] =" << vub[ibnd[1]] << std::endl;
69 std::cout << "vlb[ibnd[2]] =" << vlb[ibnd[2]] << " vub[ibnd[2]] =" << vub[ibnd[2]] << std::endl;
70 #endif
71 
72  // compute the 6 permutations of the 3 variables
73  permutation3(ind,ibnd);
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) {
203 #ifdef DEBUG
204  std::cout << " -- case 1 --" << std::endl;
205 #endif
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) {
245 #ifdef DEBUG
246  std::cout << " -- case 2 --" << std::endl;
247 #endif
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 
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 
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) {
332 #ifdef DEBUG
333  std::cout << " -- case 3 --" << std::endl;
334 #endif
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]) {
375 #ifdef DEBUG
376 std::cout << "else if " << std::endl;
377 #endif
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 {
409 #ifdef DEBUG
410 std::cout << "else " << std::endl;
411 #endif
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]]
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]]) ))
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] ) {
476 
477 #ifdef DEBUG
478 std::cout << "2 - if " << std::endl;
479 #endif
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++) {
488 
489  cutIndices [ii][0] = v1;
490  cutIndices [ii][1] = v2;
491  cutIndices [ii][2] = v3;
492  cutIndices [ii][3] = v4;
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]) {
506 #ifdef DEBUG
507 std::cout << "2 - else if" << std::endl;
508 #endif
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 
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
542 std::cout << "2 - another else if" << std::endl;
543 std::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];
563 
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 
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;
588  bnd[last+4] = (-(theta1*xU1)/(xL1-xU1)) - xL1*xL2*xU3 - xL1*xU2*xL3 + xU1*xL2*xL3;
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) {
603 #ifdef DEBUG
604  std::cout << " -- case 4 --" << std::endl;
605 #endif
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) {
644 #ifdef DEBUG
645  std::cout << " -- case 5 --" << std::endl;
646  std::cout << "v1 = " << v1 << " v2 =" << v2 << " v3 =" << v3 << std::endl;
647 #endif
648 
649  defcons_size = 12;
650  prepareVectors (defcons_size);
651 
652  // compute the permutations of the 3 variables
653  ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
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];
656  int i, flagg=0, idx=0;
657  i = 0;
658  while(i < 2 && flagg == 0) {
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];
671 #ifdef DEBUG
672  std::cout << "v1 = " << v1 << " v2 =" << v2 << " v3 =" << v3 << std::endl;
673 #endif
674 
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];
678 
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 
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]) {
690 #ifdef DEBUG
691  std::cout << " -- 5 if --" << std::endl;
692 #endif
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 {
707 #ifdef DEBUG
708  std::cout << " -- 5 else --" << std::endl;
709 #endif
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) {
741 #ifdef DEBUG
742  std::cout << " -- case 6 --" << std::endl;
743 #endif
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) {
782 #ifdef DEBUG
783  std::cout << " -- case 7 --" << std::endl;
784 #endif
785 
786  defcons_size = 12;
787  prepareVectors (defcons_size);
788 
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])) {
791 #ifdef DEBUG
792  std::cout << " -- epsilonT --" << std::endl;
793 #endif
794  } else {
795  // compute the 6 permutations of the 3 variables
796  ibnd[0] = v1; ibnd[1] = v2; ibnd[2] = v3;
797  permutation3(ind,ibnd);
798  int i, flagg=0, idx=0;
799  i = 0;
800  while(i < 6 && flagg == 0) {
801  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]]
802  <= vlb[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vub[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]] &&
803  vub[ind[i][0]]*vlb[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vub[ind[i][1]]*vub[ind[i][2]]
804  <= vub[ind[i][0]]*vub[ind[i][1]]*vlb[ind[i][2]] + vlb[ind[i][0]]*vlb[ind[i][1]]*vub[ind[i][2]])
805  {
806  idx = i; // store the index of the permutation satisfying the condition
807  flagg = 1; // condition is satisfied
808  }
809  i++;
810  }
811  if (flagg==0) {
812  std::cout << "ERROR!!!" << std::endl; exit(0);
813  }
814 
815  v1 = ind[idx][0]; v2 = ind[idx][1]; v3 = ind[idx][2];
816  }
817 
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];
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 
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 
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) {
860 #ifdef DEBUG
861  std::cout << " -- case 8 --" << std::endl;
862 #endif
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 
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
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) {
952 #ifdef DEBUG
953  std::cout << " -- case 9 --" << std::endl;
954 #endif
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 
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 
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) {
1047 #ifdef DEBUG
1048  std::cout << " -- case 10 --" << std::endl;
1049 #endif
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]]
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]])
1063  {
1064  idx = i; // store the index of the permutation satisfying the condition
1065  flagg = 1; // condition is satisfied
1066  }
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];
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 
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
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  }
1131 #ifdef DEBUG
1132 std::cout << ii << ") cutLb =" << cutLb[ii] << " " << "cutUb = " << cutUb[ii] << std::endl;
1133 #endif
1134  }
1135 
1136  for (int i=0; i<6; i++)
1137  delete [] ind [i];
1138 
1139  delete [] ibnd;
1140  delete [] bnd;
1141  delete [] ind;
1142 }
1143 
1144 
1145 // generate cuts for trilinear expressions
1146 void 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 [4];
1154 
1155  enum auxSign sign = cg -> Problem () -> Var (w -> Index ()) -> sign ();
1156 
1157  for (int i=0; i<3; i++)
1158  varInd [i] = args [i] -> Index ();
1159 
1160  varInd [3] = w -> Index ();
1161 
1162  std::vector <std::vector <int> > cutIndices;
1163  std::vector <std::vector <double> > cutCoeff;
1164  std::vector <double> cutLb, cutUb;
1165 
1166  // These cuts rely on three non-fixed variables: if l_i=u_i for i in
1167  // 1,2,3, there is a NaN due to division by (u_i - l_i)
1168 
1169  int n_var_fixed = 0, isFixed [3] = {0,0,0};
1170  double fixed_prod = 1.;
1171 
1172  for (int i=0; i<3; ++i) {
1173 
1174  double lb, ub;
1175 
1176  ArgList () [i] -> getBounds (lb, ub);
1177 
1178  if (fabs (ub - lb) < COUENNE_EPS) {
1179  isFixed [i] = 1;
1180  ++ n_var_fixed;
1181  fixed_prod *= (.5*(lb+ub));
1182  }
1183  }
1184 
1185  if (n_var_fixed) {
1186 
1187  if ((fixed_prod < COUENNE_EPS) || (n_var_fixed == 3)) {
1188 
1189  if (!(cg->createCut (cs,
1190  (sign == expression::AUX_LEQ) ? - COIN_DBL_MAX : fixed_prod,
1191  (sign == expression::AUX_GEQ) ? COIN_DBL_MAX : fixed_prod,
1192  w -> Index (), 1))) {
1193 
1194  cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "exprTriLin: variable should be fixed but cut can't be added: ");
1195  if (cg -> Problem () -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_CONVEXIFYING)) w -> print ();
1196  cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "\n");
1197  }
1198 
1199  } else {
1200 
1201  // We have either one or two fixed variables
1202 
1203  switch (n_var_fixed) {
1204 
1205  case 1: // this is w = f1 * v2 * v3, where one variable is
1206  // fixed. Revert to exprMul
1207 
1208  {
1209  int
1210  xi = (!isFixed [0]) ? 0 : (!isFixed [1]) ? 1 : 2,
1211  yi = (!isFixed [2]) ? 2 : (!isFixed [1]) ? 1 : 0;
1212 
1213  assert (xi != 2);
1214  assert (yi != 0);
1215 
1216  xi = varInd [xi];
1217  yi = varInd [yi];
1218 
1219  double
1220  *lb = cg -> Problem () -> Lb (),
1221  *ub = cg -> Problem () -> Ub (),
1222  xl = lb [xi],
1223  xu = ub [xi],
1224  yl = lb [yi],
1225  yu = ub [yi],
1226  wl = lb [varInd [3]],
1227  wu = ub [varInd [3]];
1228 
1229  unifiedProdCuts (cg, cs,
1230  xi, (*(arglist_ [0])) (), xl, xu,
1231  yi, (*(arglist_ [1])) (), yl, yu,
1232  varInd [3], (*w) (), wl, wu,
1233  chg, sign);
1234  }
1235 
1236  break;
1237 
1238  case 2:
1239 
1240  { // easy... cut is w = (f1 * f2) * v3, where f* are fixed and v* is remaining
1241 
1242  int varInd = (!isFixed [0]) ? 0 : (!isFixed [1]) ? 1 : 2;
1243 
1244  double
1245  lb = ((sign == expression::AUX_LEQ) ? - COIN_DBL_MAX : 0),
1246  ub = ((sign == expression::AUX_GEQ) ? COIN_DBL_MAX : 0);
1247 
1248  if (!(cg->createCut (cs, lb, ub, w -> Index (), 1, ArgList () [varInd] -> Index (), -fixed_prod))) {
1249  cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "exprTriLin: variable should be fixed but cut can't be added: ");
1250  if (cg -> Problem () -> Jnlst () -> ProduceOutput (Ipopt::J_ERROR, J_CONVEXIFYING)) w -> print ();
1251  cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "\n");
1252  }
1253  }
1254 
1255  break;
1256 
1257  default:
1258  cg -> Problem () -> Jnlst () -> Printf (Ipopt::J_ERROR, J_CONVEXIFYING, "exprTriLin: Error, there should be one or two fixed variables");
1259  break;
1260  }
1261  }
1262 
1263  return;
1264  }
1265 
1266  TriLinCuts (cg -> Problem () -> Lb (),
1267  cg -> Problem () -> Ub (),
1268  varInd,
1269  cutIndices, cutCoeff,
1270  cutLb, cutUb);
1271 
1272  // sanity check on returned vectors
1273  assert (cutIndices.size () == cutCoeff.size () &&
1274  cutIndices.size () == cutLb.size () &&
1275  cutIndices.size () == cutUb.size ());
1276 
1277  //printf ("trilinear cuts:\n");
1278 
1279  for (int i = (int) cutIndices.size (); i--;) {
1280 
1281  // Fix right hand sides: all cuts have coefficients of w equal to
1282  // one, but they might be inequality-type auxiliaries.
1283 
1284  exprAux *waux = dynamic_cast <exprAux *> (w);
1285 
1286  if (waux) {
1287  if (waux -> sign () == expression::AUX_LEQ) cutLb [i] = - COUENNE_INFINITY;
1288  else if (waux -> sign () == expression::AUX_GEQ) cutUb [i] = COUENNE_INFINITY;
1289  }
1290 
1291  if ((cutLb [i] > - COUENNE_INFINITY/10) ||
1292  (cutUb [i] < COUENNE_INFINITY/10)) {
1293 
1294  int
1295  size = (int) cutIndices [i].size (),
1296  *ind = new int [size];
1297 
1298  double *coe = new double [size];
1299 
1300  int cardCut = 0;
1301  for(int fmi=0; fmi<4; fmi++) {
1302  if(fabs(cutCoeff[i][fmi]) > 1e-8) {
1303  ind[cardCut] = cutIndices[i][fmi];
1304  coe[cardCut] = cutCoeff[i][fmi];
1305  cardCut++;
1306  }
1307  }
1308 
1309  OsiRowCut cut (cutLb [i], cutUb [i], 4, cardCut, ind, coe);
1310  //cut.print ();
1311 
1312  if (cg -> Problem () -> bestSol ()) {
1313 
1314  // check validity of cuts by verifying they don't cut the
1315  // optimum in the node
1316 
1317  double *sol = cg -> Problem () -> bestSol ();
1318  const double
1319  *lb = cg -> Problem () -> Lb (),
1320  *ub = cg -> Problem () -> Ub ();
1321 
1322  int nVars = cg -> Problem () -> nVars ();
1323 
1324  bool optIn = true;
1325 
1326  for (int i=0; i< nVars; i++)
1327  if ((sol [i] < lb [i] - COUENNE_EPS) ||
1328  (sol [i] > ub [i] + COUENNE_EPS)) {
1329  optIn = false;
1330  break;
1331  }
1332 
1333  if (optIn) {
1334 
1335  for (unsigned int i=0; i<cutIndices.size (); i++) {
1336 
1337  double chs = 0.;
1338 
1339  for (unsigned int j=0; j<cutIndices[i].size(); j++)
1340  chs += cutCoeff [i] [j] * sol [cutIndices [i] [j]];
1341 
1342  if ((chs < cutLb [i] - COUENNE_EPS) ||
1343  (chs > cutUb [i] + COUENNE_EPS)) {
1344 
1345  printf ("cut %d violates optimum:\n", i);
1346 
1347  if (cutLb [i] > -COUENNE_INFINITY) printf ("%g <= ", cutLb [i]);
1348  for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g x%d ", cutCoeff [i] [j], cutIndices [i] [j]); printf ("\n = ");
1349  for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g *%g ", cutCoeff [i] [j], sol [cutIndices [i] [j]]); printf ("\n = ");
1350  for (unsigned int j=0; j<cutIndices[i].size(); j++) printf ("%+g ", cutCoeff [i] [j] * sol [cutIndices [i] [j]]); printf (" = %g", chs);
1351  if (cutUb [i] < COUENNE_INFINITY) printf (" <= %g", cutUb [i]);
1352  printf ("\n");
1353 
1354  }
1355  }
1356  }
1357  }
1358 
1359  cs.insert (cut);
1360  }
1361  }
1362  //delete [] varInd;
1363 }
Couenne::CouenneCutGenerator
Cut Generator for linear convexifications.
Definition: CouenneCutGenerator.hpp:49
Couenne::exprOp::print
virtual void print(std::ostream &out=std::cout, bool=false) const
I/O.
Definition: exprOp.cpp:42
Couenne::unifiedProdCuts
void unifiedProdCuts(const CouenneCutGenerator *, OsiCuts &, int, CouNumber, CouNumber, CouNumber, int, CouNumber, CouNumber, CouNumber, int, CouNumber, CouNumber, CouNumber, t_chg_bounds *, enum expression::auxSign)
unified convexification of products and divisions
Definition: unifiedProdCuts.cpp:96
CouenneExprTrilinear.hpp
permutation3
void permutation3(int **ind, int *ibnd)
Definition: conv-exprTrilinear-gencuts.cpp:29
CouenneExprMul.hpp
Couenne::t_chg_bounds
status of lower/upper bound of a variable, to be checked/modified in bound tightening ...
Definition: CouenneTypes.hpp:66
Couenne::exprTrilinear::getBounds
virtual void getBounds(expression *&, expression *&)
Get lower and upper bound of an expression (if any)
Definition: conv-exprTrilinear.cpp:20
Couenne::J_CONVEXIFYING
const Ipopt::EJournalCategory J_CONVEXIFYING(Ipopt::J_USER3)
TriLinCuts
void TriLinCuts(double *vlb, double *vub, int *varIndices, std::vector< std::vector< int > > &cutIndices, std::vector< std::vector< double > > &cutCoeff, std::vector< double > &cutLb, std::vector< double > &cutUb)
generate convexification cuts for constraint w = x*y*z
Definition: conv-exprTrilinear-gencuts.cpp:41
Couenne::exprTrilinear::generateCuts
void generateCuts(expression *w, OsiCuts &cs, const CouenneCutGenerator *cg, t_chg_bounds *=NULL, int=-1, CouNumber=-COUENNE_INFINITY, CouNumber=COUENNE_INFINITY)
generate equality between *this and *w
Definition: conv-exprTrilinear-gencuts.cpp:1146
j
static char * j
Definition: OSdtoa.cpp:3622
EPSILONT
#define EPSILONT
Definition: conv-exprTrilinear-gencuts.cpp:24
Couenne::expression::Image
virtual expression * Image() const
return pointer to corresponding expression (for auxiliary variables only)
Definition: CouenneExpression.hpp:100
e
void fint fint fint real fint real real real real real real real real real * e
Definition: BonBqpdSolver.cpp:27
Couenne::expression::Index
virtual int Index() const
Return index of variable (only valid for exprVar and exprAux)
Definition: CouenneExpression.hpp:73
Couenne::expression::auxSign
auxSign
&quot;sign&quot; of the constraint defining an auxiliary.
Definition: CouenneExpression.hpp:55
CouenneCutGenerator.hpp
Couenne::CouenneCutGenerator::createCut
int createCut(OsiCuts &, CouNumber, CouNumber, int, CouNumber, int=-1, CouNumber=0., int=-1, CouNumber=0., bool=false) const
create cut and check violation. Insert and return status
Definition: createCuts.cpp:32
COUENNE_EPS
#define COUENNE_EPS
Definition: CouennePrecisions.hpp:19
int
static int
Definition: OSdtoa.cpp:2173
CouenneExprAux.hpp
Couenne::exprOp::ArgList
expression ** ArgList() const
return argument list
Definition: CouenneExprOp.hpp:66
Couenne::exprOp::arglist_
expression ** arglist_
argument list is an array of pointers to other expressions
Definition: CouenneExprOp.hpp:35
Couenne::CouNumber
double CouNumber
main number type in Couenne
Definition: CouenneTypes.hpp:100
COUENNE_INFINITY
#define COUENNE_INFINITY
Definition: CouennePrecisions.hpp:32
Couenne::exprAux
Auxiliary variable.
Definition: CouenneExprAux.hpp:32
Couenne::expression
Expression base class.
Definition: CouenneExpression.hpp:48
prepareVectors
#define prepareVectors(a)
CouenneTypes.hpp
CouenneProblem.hpp
w
void fint fint fint real fint real real real real real real real real * w
Definition: BonBqpdSolver.cpp:27
Generated by   doxygen 1.8.5