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

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

moved include files to make them doxygenable. Introduced three-way branching, with fixed intervals for now. Added check for small bound interval within all generateCuts()

File size: 3.3 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  // special case: l and u are very close, replace function with
81  // linear term
82
83  if (fabs (u - l) < COUENNE_EPS) {
84
85    CouNumber x0 = 0.5 * (u+l);
86    cg -> createCut (cs, 2/x0, 0, w_ind, 1., x_ind, 1/(x0*x0));
87    return;
88  }
89
90  // choose sampling points. If unbounded, bound using a rule of thumb
91
92  int ns = cg -> nSamples ();
93
94  if      (l < - COUENNE_INFINITY) l = ns * (u-1); // (-infinity, u] where u < 0
95  else if (u >   COUENNE_INFINITY) u = ns * (l+1); // [l, +infinity) where l > 0
96
97  // make bounds nonzero
98
99  if (fabs (l) < COUENNE_EPS) {
100    /*l /= (COUENNE_EPS / MIN_DENOMINATOR);
101      if (fabs (l) < COUENNE_EPS) */
102    l = (l<0) ? - MIN_DENOMINATOR : MIN_DENOMINATOR;
103  }
104
105  if (fabs (u) < COUENNE_EPS) {
106    /*u /= (COUENNE_EPS / MIN_DENOMINATOR);
107      if (fabs (u) < COUENNE_EPS) */
108      u = (u<0) ? - MIN_DENOMINATOR : MIN_DENOMINATOR;
109  }
110
111  // bound
112  cg -> addEnvelope
113    (cs, (l > 0) ? +1 : -1, 
114     inv, oppInvSqr, w_ind, x_ind, 
115     (cg -> isFirst ()) ?
116       // place it somewhere in the interval (we won't care)
117       ((l > COUENNE_EPS) ? l : u) :
118       // not first call, gotta replace it where it gives deepest cut
119       powNewton ((*argument_) (), (*aux) (), inv, oppInvSqr, inv_dblprime),
120     l, u);
121
122  delete xle; 
123  delete xue;
124}
Note: See TracBrowser for help on using the repository browser.