source: branches/Couenne/Couenne/src/convex/operators/conv-exprInv.cpp @ 520

Last change on this file since 520 was 520, checked in by pbelotti, 13 years ago

limit power upper/lower bounds to get reasonable coefficients

File size: 3.1 KB
Line 
1/*
2 * Name:    conv-exprInv.cpp
3 * Author:  Pietro Belotti
4 * Purpose: convexification and bounding methods for the inverse operator
5 *
6 * This file is licensed under the Common Public License (CPL)
7 */
8
9#include <CouenneTypes.h>
10
11#include <exprInv.h>
12#include <exprClone.h>
13#include <exprConst.h>
14#include <exprMin.h>
15#include <exprOpp.h>
16#include <exprDiv.h>
17#include <exprSum.h>
18#include <exprMul.h>
19#include <exprPow.h>
20
21#include <CouenneProblem.h>
22#include <CouenneCutGenerator.h>
23
24// compute upper- and lower bounds of the expression w = 1/f(x) given
25// bounds of f(x)
26
27void exprInv::getBounds (expression *&lb, expression *&ub) {
28
29  expression *lba, *uba;
30  argument_ -> getBounds (lba, uba);
31
32  expression **all = new expression * [6];
33  all [0] = new exprConst (0);       all [1] = new exprConst (- COUENNE_INFINITY);  // l<0<u
34  all [2] = new exprOpp   (lba);     all [3] = new exprInv   (uba);                 // 0<l<u
35  all [4] = new exprClone (uba);     all [5] = new exprInv   (new exprClone (uba)); // l<u<0
36
37  lb = new exprMin (all, 6);
38
39  expression **alu = new expression * [6];
40  alu [0] = new exprConst (0);       alu [1] = new exprConst (COUENNE_INFINITY);   // l<0<u
41  alu [2] = new exprClone (all [2]); alu [3] = new exprInv (new exprClone (lba));  // 0<l<u
42  alu [4] = new exprClone (uba);     alu [5] = new exprInv (new exprClone (lba));  // l<u<0
43
44  ub = new exprMin (alu, 6);
45}
46
47
48#define MIN_DENOMINATOR 1e-10
49
50// generate convexification cut for constraint w = 1/x
51
52void exprInv::generateCuts (exprAux *aux, const OsiSolverInterface &si, 
53                            OsiCuts &cs, const CouenneCutGenerator *cg) {
54
55  // get bounds of numerator and denominator
56
57  expression *xle, *xue;
58
59  argument_ -> getBounds (xle, xue);
60
61  CouNumber l = (*xle) (), 
62            u = (*xue) ();
63
64  // if the denominator's bound interval has 0 as internal point,
65  // there is no convexification
66
67  if ((l < - COUENNE_EPS) && 
68      (u >   COUENNE_EPS)) 
69    return;
70
71  expression *wle, *wue;
72
73  aux -> getBounds (wle, wue);
74
75  CouNumber x;
76
77  int w_ind = aux       -> Index (), 
78      x_ind = argument_ -> Index ();
79
80  // choose sampling points. If unbounded, bound using a rule of thumb
81
82  int ns = cg -> nSamples ();
83
84  if      (l < - COUENNE_INFINITY + 1) l = ns * (u-1); // (-infinity, u] where u < 0
85  else if (u >   COUENNE_INFINITY - 1) u = ns * (l+1); // [l, +infinity) where l > 0
86
87  // make bounds nonzero
88
89  if (fabs (l) < COUENNE_EPS) {
90    /*l /= (COUENNE_EPS / MIN_DENOMINATOR);
91      if (fabs (l) < COUENNE_EPS) */
92    l = (l<0) ? - MIN_DENOMINATOR : MIN_DENOMINATOR;
93  }
94
95  if (fabs (u) < COUENNE_EPS) {
96    /*u /= (COUENNE_EPS / MIN_DENOMINATOR);
97      if (fabs (u) < COUENNE_EPS) */
98      u = (u<0) ? - MIN_DENOMINATOR : MIN_DENOMINATOR;
99  }
100
101  // bound
102  cg -> addEnvelope
103    (cs, (l > 0) ? +1 : -1, 
104     inv, oppInvSqr, w_ind, x_ind, 
105     (cg -> isFirst ()) ?
106       // place it somewhere in the interval (we won't care)
107       ((l > COUENNE_EPS) ? l : u) :
108       // not first call, gotta replace it where it gives deepest cut
109       powNewton ((*argument_) (), (*aux) (), inv, oppInvSqr, inv_dblprime),
110     l, u);
111
112  delete xle; 
113  delete xue;
114}
Note: See TracBrowser for help on using the repository browser.