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

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

hopefully fixed the last const_iterator problem

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