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

Last change on this file since 939 was 939, checked in by pbelotti, 8 years ago

cosmetic changes

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