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

Last change on this file since 956 was 956, checked in by pbelotti, 7 years ago

restored working version. Trying to fix compilation errors on const_iterators

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