source: trunk/Couenne/src/cut/sdpcuts/CutGen.cpp @ 945

Last change on this file since 945 was 945, checked in by stefan, 8 years ago

generateCuts in Cgl >= 0.58 will not be const anymore

  • Property svn:keywords set to Author Date Id Revision
File size: 14.5 KB
Line 
1/* $Id: CutGen.cpp 945 2013-04-06 20:25:21Z stefan $
2 *
3 * Name:    CutGen.cpp
4 * Authors: Andrea Qualizza
5 *          Pietro Belotti
6 * Purpose: Generation of all cuts
7 *
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11#include "stdlib.h"
12#include "stdio.h"
13#include "math.h"
14
15#include "CoinTime.hpp"
16#include "CoinHelperFunctions.hpp"
17
18#include "CouenneExprVar.hpp"
19#include "CouenneExprClone.hpp"
20#include "operators/CouenneExprMul.hpp"
21#include "CouenneProblem.hpp"
22#include "CouenneMatrix.hpp"
23#include "CouenneSdpCuts.hpp"
24
25#include "dsyevx_wrapper.hpp"
26
27//#define DEBUG
28
29const bool WISE_SPARSIFY = true;
30
31#define SPARSIFY_MAX_CARD  10000
32#define SPARSIFY_NEW_NZ_THRESHOLD 0.70
33
34#define EV_TOL 1e-13
35
36using namespace Couenne;
37
38
39// sdpcut separator -- all minors
40void CouenneSdpCuts::generateCuts (const OsiSolverInterface &si, OsiCuts &cs, 
41                                   const CglTreeInfo info)
42#if CGL_VERSION_MAJOR == 0 && CGL_VERSION_MINOR <= 57
43  const
44#endif
45  {
46
47  // only at root node once
48
49  if ((info . level + info . pass > 4)) return;
50
51  problem_ -> domain () -> push (&si, &cs);
52
53  for (std::vector <CouenneExprMatrix *>::const_iterator
54         minor  = minors_. begin ();
55       minor   != minors_. end   (); ++minor)
56
57    genCutSingle (*minor, si, cs, info);
58
59  problem_ -> domain () -> pop ();
60}
61
62
63// sdpcut separator -- one minor at a time
64void CouenneSdpCuts::genCutSingle (CouenneExprMatrix * const & minor, 
65                                   const OsiSolverInterface &si, 
66                                   OsiCuts &cs, 
67                                   const CglTreeInfo info) const {
68#ifdef DEBUG
69  printf ("Generating cut on minor -> ");
70  minor -> print ();
71#endif
72
73  std::vector <expression *> &varInd = minor -> varIndices ();
74
75  int
76    n = (int) (minor -> size ()),
77    m,
78    nVecs = (numEigVec_ < 0) ? n : numEigVec_,
79
80    **indA  = new int * [n], // contains indices of matrix A below
81    *indMap = new int [problem_ -> nVars ()],
82
83    compressed_index = 0;
84
85  double 
86    *A = new double [n * n];
87
88  // Fill in inverse map
89  for (std::vector <expression *>::const_iterator
90         i  = varInd . begin ();
91       i   != varInd . end   (); ++i)
92    indMap [(*i) -> Index ()] = compressed_index++;
93
94  // ---------------------------------------------------------------------
95
96  for (int i=0; i<n; ++i)
97    CoinFillN (indA [i] = new int [n], n, -2);
98
99  // build index and value matrices
100
101  // Get minors_
102
103  CoinFillN (A, n * n, 0.);
104
105#ifdef DEBUG
106  printf ("Components\n\n");
107#endif
108
109  int indI = 0;
110
111  for (std::set <std::pair <int, CouenneSparseVector *> >::iterator
112         i  = minor -> getRows () . begin ();
113       i   != minor -> getRows () . end   (); ++i, ++indI) {
114
115#ifdef DEBUG
116    printf ("row %d [var %d]: ", indI, i -> first); fflush (stdout);
117#endif
118
119    int
120      majInd = (i -> first == problem_ -> nVars ()) ? n-1 : indMap [i -> first],
121      indJ   = 0;
122
123    for (std::set <CouenneScalar *>::iterator
124           j  = (*i) . second -> getElements () . begin ();
125         j   != (*i) . second -> getElements () . end   (); ++j, ++indJ) {
126
127#ifdef DEBUG
128      printf ("[r%d,v%d,%d,%g,", indJ, 
129              (*j) -> getIndex (), 
130              (*j) -> getElem () -> Index (), 
131              (*((*j) -> getElem ()))  ());
132      (*j) -> getElem () -> print (); printf ("] ");
133      fflush (stdout);
134#endif
135
136      int minInd = ((*j) -> getIndex () == problem_ -> nVars ()) ? (n-1) : (indMap [(*j) -> getIndex ()]);
137
138      expression *Elem = (*j) -> getElem ();
139
140      A [majInd * n + minInd] = (*Elem) (); // evaluate variable at this position of matrix
141      //A [indJ * n + indI] = (*Elem) ();
142
143      indA [majInd] [minInd] = Elem -> Index ();
144      //indA [indJ] [indI] =
145    }
146
147#ifdef DEBUG
148    printf ("\n");
149#endif
150  }
151
152  delete [] indMap;
153
154  // fill in non-existing auxiliaries in X (if any) with their hypothetical product
155
156  for   (int i=0; i<n-1; ++i) // get to n-1 since varIndices not defined afterward
157    for (int j=i; j<n-1; ++j)
158      if (indA [i] [j] == -2) // never happens on border row/column
159        A [i * n + j] = 
160        A [j * n + i] = 
161          (*(varInd [i])) () * 
162          (*(varInd [j])) ();
163
164#ifdef DEBUG
165  for (int i=0; i<n; ++i) {
166    for (int j=0; j<n; ++j)
167      printf ("[%4d,%7.2g] ", indA [i][j], A [i * n + j]);
168    printf ("\n");
169  }
170#endif
171
172  double
173    *Acopy = useSparsity_ ? CoinCopyOfArray (A, n * n) : NULL,
174    *w = NULL,
175    *z = NULL;
176
177  //  printf ("calling dsyevx NONSPARSE\n");
178
179  //------------------------------------------------------------------------------------------------------
180  dsyevx_interface (n, A, m, w, z, EV_TOL, 
181                    -COIN_DBL_MAX, onlyNegEV_ ? 0. : COIN_DBL_MAX, 
182                    1, numEigVec_ < 0 ? n : numEigVec_);
183  //------------------------------------------------------------------------------------------------------
184
185  if (m < nVecs) 
186    nVecs = m;
187
188  double
189    **work_ev = new double * [m];
190
191  for (int i=0; i < nVecs; i++) {
192
193    if (w [i] >= 0) {
194      nVecs = i; // updated nVecs
195      break;
196    }
197
198    work_ev [i] = z + (i*n);
199
200#ifdef SCALE_EIGENV
201    double scaling_factor = sqrt (n);
202
203    for (int j=0; j<n; j++)
204      work_ev [i] [j] *= scaling_factor;
205#endif
206  }
207
208  // Generate sdp (possibly dense) cuts //////////////////////////////////////////////////////////
209
210  for (int i=0; i<nVecs; i++)
211    genSDPcut (si, cs, minor, work_ev [i], work_ev [i], indA);
212
213  int
214    wise_evdec_num    = 0,
215    card_sparse_v_mat = 0,
216    min_nz;
217
218  double **sparse_v_mat = NULL;
219
220  if (useSparsity_) {
221
222    sparse_v_mat = new double*[SPARSIFY_MAX_CARD];
223    for (int i=0; i<SPARSIFY_MAX_CARD; i++)
224      sparse_v_mat[i] = new double [n];
225
226    min_nz = ceil (n * SPARSIFY_NEW_NZ_THRESHOLD);
227    card_sparse_v_mat = 0;
228
229    sparsify2 (n, A, sparse_v_mat, &card_sparse_v_mat, min_nz, &wise_evdec_num);
230
231    for (int k=0; k < card_sparse_v_mat; k++)
232      genSDPcut (si, cs, minor, sparse_v_mat [k], sparse_v_mat [k], indA);
233
234    ////////////////////////////////////////////
235
236    for (int i=0; i<nVecs; ++i) {
237
238      card_sparse_v_mat = 0;
239      double *v = work_ev[i];
240
241      sparsify (WISE_SPARSIFY, i, w [i], v, n, Acopy, sparse_v_mat, &card_sparse_v_mat, &wise_evdec_num);
242
243      for (int k=0; k < card_sparse_v_mat; k++) {
244
245        genSDPcut (si, cs, minor, sparse_v_mat [k], sparse_v_mat [k], indA);
246
247        if (useSparsity_)
248          additionalSDPcuts (si, cs, minor, n, Acopy, sparse_v_mat[k], indA);
249      }
250    }
251  }
252
253  for (int i=0; i<n; ++i)
254    delete [] indA [i];
255  delete [] indA;
256
257  if (useSparsity_) {
258
259    for (int i=0; i < SPARSIFY_MAX_CARD; i++)
260      delete [] sparse_v_mat[i];
261
262    delete [] sparse_v_mat;
263    delete [] Acopy;
264  }
265
266  delete [] z;
267  delete [] w;
268  delete [] A;
269  delete [] work_ev;
270}
271
272
273/************************************************************************/
274void CouenneSdpCuts::genSDPcut (const OsiSolverInterface &si,
275                                OsiCuts &cs, 
276                                CouenneExprMatrix *XX,
277                                double *v1, double *v2, 
278                                int **indA) const {
279  int
280    nterms   = 0,
281    n        = (int) (XX -> size ()),
282    nvars    = problem_ -> nVars (),
283    N        = n * n,
284    *ind     = new int [N],
285    *inverse = new int [nvars];
286
287  double
288    *coeff = new double [N],
289    *xtraC = new double [nvars],
290    rhs    = 0.;
291
292  std::vector <expression *> &varIndices = XX -> varIndices ();
293
294  CoinFillN (xtraC,   nvars, 0.);
295  CoinFillN (inverse, nvars, -1);
296
297  // for (int i=0; i<n; ++i) {
298  //   if (fabs (v1 [i]) < COUENNE_EPS) v1 [i] = 0;
299  //   if (fabs (v2 [i]) < COUENNE_EPS) v2 [i] = 0;
300  // }
301
302#ifdef DEBUG
303  printf ("Solution: (");
304  for (int i=0; i<n; i++) {
305    if (i) printf (",");
306    printf ("%g", v1 [i]);
307  }
308  printf (")\n");
309#endif
310
311  // ASSUMPTION: matrix is symmetric
312
313  bool numerics_flag = false;
314
315  // coefficients for X_ij
316  for (int i=0; (i<n) && !numerics_flag; i++)
317
318    for (int j=i; j<n; j++) {
319
320      double coeff0 = v1 [i] * v2 [j] + v1 [j] * v2 [i];
321
322      if (fabs (coeff0) < 1e-21) continue; // why 1e-21? Cbc/Clp related
323
324      int index = indA [i] [j];
325
326#ifdef DEBUG
327      printf ("{%d,%g} ", index, coeff0);
328#endif
329
330      // Three cases:
331      //
332      // 1) index >= 0: corresponds to variable of the problem,
333      //    proceed as normal
334      //
335      // 2) index == -1: constant, subtract it from right-hand side
336      //
337      // 3) index < -1: this variable (x_ij = x_i * x_j) is not in the
338      //    problem, and must be replaced somehow.
339
340#ifdef DEBUG
341      if (index == -1) printf ("found constant: %g\n", (i==j) ? coeff0 / 2 : coeff0);
342#endif
343      if (index == -1)
344
345        rhs -= ((i==j) ? coeff0 / 2 : coeff0);
346
347      else if (index < -1) { // only happens in non-bordered matrix
348
349        expression
350          *Xi = XX -> varIndices () [i],
351          *Xj = XX -> varIndices () [j];
352
353        double 
354          valXi = (*Xi) (),   
355          valXj = (*Xj) (),
356          li, lj,
357          ui, uj,
358          L, U;
359
360        Xi -> getBounds (li, ui);
361        Xj -> getBounds (lj, uj);
362
363#ifdef DEBUG
364        printf ("expression: x%d [%g,%g] * x%d [%g,%g]\n", 
365                Xi -> Index (), li, ui, 
366                Xj -> Index (), lj, uj);
367#endif
368
369        if (i==j) {
370
371          // This variable (x_ii = x_i ^ 2) is not in the
372          // problem. Replace it with its upper envelope
373          //
374          // X_ii <= li^2 + (ui^2 - li^2) / (ui-li) * (xi-li) = (ui+li) * xi - li*ui
375
376          if ((fabs (li) > COUENNE_INFINITY) || 
377              (fabs (ui) > COUENNE_INFINITY)) { 
378
379            // term would be X_ii <= inf, which adds inf to the
380            // inequality
381
382            numerics_flag = true;
383            break;
384          }
385
386          xtraC [varIndices [i] -> Index ()] += coeff0 / 2 * (li + ui); 
387#ifdef DEBUG
388          printf ("adding %g=%g*(%g+%g) (sq) to xtrac[%d=varInd[%d]]\n", 
389                  coeff0 / 2 * (li + ui), coeff0, li, ui, varIndices [i] -> Index (), i);
390#endif
391          rhs += coeff0 / 2 * (li * ui);
392
393        } else {
394
395          // This variable (x_ij = x_i * x_j) is not in the problem,
396          // and must be replaced somehow. Apply Fourier-Motzkin
397          // elimination using bounds and McCormick inequalities on
398          // the fictitious variable y_ij = x_i x_j:
399          //
400          //    y_ij >= L = min {l_i*l_j, u_i*u_j, l_i*u_j, u_i*l_j}
401          //    y_ij >= l_j x_i + l_i x_j - l_i l_j
402          //    y_ij >= u_j x_i + u_i x_j - u_i u_j
403          //
404          //    y_ij <= U = max {l_i*l_j, u_i*u_j, l_i*u_j, u_i*l_j}
405          //    y_ij <= u_j x_i + l_i x_j - u_i l_j
406          //    y_ij <= l_j x_i + u_i x_j - l_i u_j
407
408          exprMul Xij (new exprClone (Xi), 
409                       new exprClone (Xj));
410
411          Xij . getBounds (L, U);
412
413          double 
414            rhsMll = lj * valXi + li * valXj - lj * li,
415            rhsMuu = uj * valXi + ui * valXj - uj * ui;
416
417          if (coeff0 < 0) {
418
419            if (L >= CoinMax (rhsMll, rhsMuu)) 
420
421              rhs -= coeff0 * L; // L is the tightest (greatest)
422
423            else if (rhsMuu > rhsMll) {  // rhsMuu is the tightest
424
425              if ((fabs (ui) > COUENNE_INFINITY) || 
426                  (fabs (uj) > COUENNE_INFINITY)) { 
427                numerics_flag = true;
428                break;
429              }
430
431              xtraC [varIndices [i] -> Index ()] += coeff0 * uj;
432              xtraC [varIndices [j] -> Index ()] += coeff0 * ui;
433              rhs += coeff0 * uj * ui;
434
435#ifdef DEBUG
436              printf ("adding (%g,%g) = %g*(%g,%g) (rhsMuu) to xtrac[%d=varInd[%d]] and xtrac[%d=varInd[%d]]\n", 
437                      coeff0 * uj, coeff0 * ui, coeff0, uj, ui, varIndices [i] -> Index (), i, varIndices [j] -> Index (), j);
438#endif
439            } else { // rhsMll is the tightest
440
441              if ((fabs (li) > COUENNE_INFINITY) || 
442                  (fabs (lj) > COUENNE_INFINITY)) { 
443
444                numerics_flag = true;
445                break;
446              }
447
448              xtraC [varIndices [i] -> Index ()] += coeff0 * lj;
449              xtraC [varIndices [j] -> Index ()] += coeff0 * li;
450              rhs += coeff0 * lj * li;
451
452#ifdef DEBUG
453              printf ("adding (%g,%g) = %g*(%g,%g) (rhsMll) to xtrac[%d=varInd[%d]] and xtrac[%d=varInd[%d]]\n", 
454                      coeff0 * lj, coeff0 * li, coeff0, lj, li, varIndices [i] -> Index (), i, varIndices [j] -> Index (), j);
455#endif
456            }
457
458          } else { // coeff0 > 0
459
460            double
461              rhsMlu = lj * valXi + ui * valXj - lj * ui,
462              rhsMul = uj * valXi + li * valXj - uj * li;
463
464            if (U <= CoinMin (rhsMlu, rhsMul)) 
465
466              rhs -= coeff0 * U; // U is the tightest (smallest)
467
468            else if (rhsMul < rhsMlu) { // rhsMul is the tightest
469
470              if ((fabs (li) > COUENNE_INFINITY) || 
471                  (fabs (uj) > COUENNE_INFINITY)) { 
472
473                numerics_flag = true;
474                break;
475              }
476
477              xtraC [varIndices [i] -> Index ()] += coeff0 * uj;
478              xtraC [varIndices [j] -> Index ()] += coeff0 * li;
479              rhs += coeff0 * uj * li;
480
481#ifdef DEBUG
482              printf ("adding (%g,%g) = %g*(%g,%g) (rhsMul) to xtrac[%d=varInd[%d]] and xtrac[%d=varInd[%d]]\n", 
483                      coeff0 * uj, coeff0 * li, coeff0, uj, li, varIndices [i] -> Index (), i, varIndices [j] -> Index (), j);
484#endif
485
486            } else { // rhsMlu is the tightest
487
488              if ((fabs (ui) > COUENNE_INFINITY) || 
489                  (fabs (lj) > COUENNE_INFINITY)) { 
490
491                numerics_flag = true;
492                break;
493              }
494
495              xtraC [varIndices [i] -> Index ()] += coeff0 * lj;
496              xtraC [varIndices [j] -> Index ()] += coeff0 * ui;
497              rhs += coeff0 * lj * ui;
498
499#ifdef DEBUG
500              printf ("adding (%g,%g) = %g*(%g,%g) (rhsMlu) to xtrac[%d=varInd[%d]] and xtrac[%d=varInd[%d]]\n", 
501                      coeff0 * lj, coeff0 * ui, coeff0, lj, ui, varIndices [i] -> Index (), i, varIndices [j] -> Index (), j);
502#endif
503            }
504          }
505        }
506
507        ///////////////////////////////////////////////////////////
508
509      } else {
510
511#ifdef DEBUG
512        printf ("normal term: %g x_%d [terms:%d]\n", (i==j) ? (0.5 * coeff0) : (coeff0), index, nterms);
513#endif
514
515        if      (inverse [index] >= 0)
516          coeff [inverse [index]] += (i==j) ? (0.5 * coeff0) : (coeff0);
517        else {
518
519          coeff   [nterms]   = (i==j) ? (0.5 * coeff0) : (coeff0);
520          // if (inverse [index] >= 0)
521          //   printf ("duplicate index at nterms=%d: inverse[%d] = %d, value %g\n",
522          //      nterms, index, inverse [index], coeff [nterms]);
523          inverse [index]    = nterms;
524          ind     [nterms++] = index;
525        }
526      }
527#ifdef DEBUG
528      printf ("%d terms so far\n", nterms);
529#endif
530    }
531
532  if (!numerics_flag)
533    for (std::vector <expression *>::iterator
534         i  = varIndices . begin (); 
535       i   != varIndices . end   (); ++i) {
536
537    int varInd = (*i) -> Index ();
538
539    if ((inverse [varInd] >= 0) && 
540        (fabs (xtraC [varInd]) > 1e-15)) {
541#ifdef DEBUG
542      printf ("now adding %g to coeff [inv [%d] = %d]\n", xtraC [varInd], varInd, inverse [varInd]);
543#endif
544      coeff [inverse [varInd]] += xtraC [varInd];
545    } else if (fabs (xtraC [varInd]) > COUENNE_EPS) { // independent variable not appearing so far
546
547      coeff   [nterms]   = xtraC [varInd];
548      inverse [varInd]   = nterms; // probably useless
549      ind     [nterms++] = varInd;
550    }
551  }
552
553  delete [] inverse;
554  delete [] xtraC;
555
556  if (!numerics_flag) {
557
558    OsiRowCut *cut = new OsiRowCut;
559    cut -> setRow (nterms, ind, coeff);
560    cut -> setLb (rhs);
561
562    if (nterms > 0) {
563
564#ifdef DEBUG
565      printf ("SDP: separating ");
566      cut -> print ();
567#endif
568
569      CoinAbsFltEq treatAsSame (COUENNE_EPS);
570      cs.insertIfNotDuplicate (*cut, treatAsSame);
571
572      double violation = 0.;
573
574      if (problem_ -> bestSol () && ((violation = cut -> violated (problem_ -> bestSol ())) > 0.)) {
575
576        printf ("Cut violates optimal solution by %g\n", violation);
577        cut -> print ();
578      }
579    }
580
581    delete cut;
582  }
583
584  delete [] ind;
585  delete [] coeff;
586}
Note: See TracBrowser for help on using the repository browser.