source: branches/Couenne/Couenne/src/convex/operators/conv-exprSinCos.cpp @ 582

Last change on this file since 582 was 582, checked in by andreasw, 12 years ago

minor changes to make code compile with MSVC++; reran autotools

File size: 9.0 KB
Line 
1/*
2 * Name:    conv-exprSinCos.cpp
3 * Author:  Pietro Belotti
4 * Purpose: convexification methods for sines and cosines
5 *
6 * This file is licensed under the Common Public License (CPL)
7 */
8
9#include <math.h>
10#ifndef M_PI
11# define M_PI 3.14159265358979323846
12#endif
13#ifndef M_PI_2
14# define M_PI_2 1.57079632679489661923
15#endif
16
17#include <OsiSolverInterface.hpp>
18#include <CouenneTypes.h>
19#include <CouenneCutGenerator.h>
20#include <exprSin.h>
21#include <exprCos.h>
22#include <exprAux.h>
23
24#define NEW_TRIG
25
26/// convex cuts for sine or cosine
27int trigEnvelope (const CouenneCutGenerator *, OsiCuts &,
28                   exprAux *, expression *, enum cou_trig);
29
30
31/// add four cuts with slope 1 and -1
32int addHexagon (const CouenneCutGenerator *, // cut generator that has called us
33                 OsiCuts &,      // cut set to be enriched
34                 enum cou_trig,  // sine or cosine
35                 exprAux *,      // auxiliary variable
36                 expression *);  // argument of cos/sin (should be a variable)
37
38
39/// generate convexification cut for constraint w = sin (this)
40
41void exprSin::generateCuts (exprAux *w, const OsiSolverInterface &si, 
42                            OsiCuts &cs, const CouenneCutGenerator *cg,
43                            t_chg_bounds *chg) {
44
45  int wi = w -> Index ();
46
47  if (chg && !(cg -> isFirst ()) && 
48      (chg [wi].lower == UNCHANGED) && 
49      (chg [wi].upper == UNCHANGED))
50    return;
51
52#ifdef NEW_TRIG
53  if (trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_SINE) == 0)
54#else
55  if (addHexagon (cg, cs, COU_SINE, w, w -> Image () -> Argument()) == 0)
56#endif
57    {
58
59    }
60}
61
62
63/// generate convexification cut for constraint w = cos (this)
64
65void exprCos::generateCuts (exprAux *w, const OsiSolverInterface &si, 
66                            OsiCuts &cs, const CouenneCutGenerator *cg,
67                            t_chg_bounds *chg) {
68
69  int wi = w -> Index ();
70
71  if (chg && !(cg -> isFirst ()) && 
72      (chg [wi].lower == UNCHANGED) && 
73      (chg [wi].upper == UNCHANGED))
74    return;
75
76#ifdef NEW_TRIG
77  if (trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_COSINE) == 0) 
78#else
79  if (addHexagon (cg, cs, COU_COSINE, w, w -> Image () -> Argument()) == 0)
80#endif
81    {
82
83    }
84}
85
86
87/// add lateral edges of the hexagon providing
88
89int addHexagon (const CouenneCutGenerator *cg, // cut generator that has called us
90                 OsiCuts &cs,       // cut set to be enriched
91                 enum cou_trig tt,  // sine or cosine
92                 exprAux *aux,      // auxiliary variable
93                 expression *arg) { // argument of cos/sin (should be a variable)
94
95  // AW 2007-06-11: The following doesn't compile with MSVC++ because
96  // sin and cos are ambiguous
97  //unary_function fn = (tt == COU_SINE) ? sin : cos;
98
99  // retrieve argument bounds
100  expression *lbe, *ube;
101  arg -> getBounds (lbe, ube);
102
103  CouNumber lb = (*lbe) (), 
104            ub = (*ube) ();
105
106  delete lbe;
107  delete ube;
108
109  int ncuts = 0,
110    x_ind = arg -> Index (),
111    w_ind = aux -> Index ();
112
113  if (fabs (ub - lb) < COUENNE_EPS) {
114
115    CouNumber x0 = 0.5 * (ub+lb), f, fp;
116
117    if (tt == COU_SINE) {f = sin (x0); fp =  cos (x0);}
118    else                {f = cos (x0); fp = -sin (x0);}
119
120    return cg -> createCut (cs, f - fp*x0, 0, w_ind, 1., x_ind, -fp);
121  }
122
123  // add  /    \ envelope
124  //      \    /
125
126  // left
127  if (lb > -COUENNE_INFINITY) { // if not unbounded
128    if (tt == COU_SINE) {
129      ncuts += cg -> createCut (cs, sin (lb) - lb, -1, w_ind, 1., x_ind, -1.); // up:  w-x <= f lb - lb
130      ncuts += cg -> createCut (cs, sin (lb) + lb, +1, w_ind, 1., x_ind,  1.); // dn:  w+x >= f lb + lb
131    }
132    else {
133      ncuts += cg -> createCut (cs, cos (lb) - lb, -1, w_ind, 1., x_ind, -1.); // up:  w-x <= f lb - lb
134      ncuts += cg -> createCut (cs, cos (lb) + lb, +1, w_ind, 1., x_ind,  1.); // dn:  w+x >= f lb + lb
135    }
136  }
137
138  // right
139  if (ub <  COUENNE_INFINITY) { // if not unbounded
140    if (tt == COU_SINE) {
141      ncuts += cg -> createCut (cs, sin (ub) - ub, +1, w_ind, 1., x_ind, -1.); // dn: w - x >= f ub - ub
142      ncuts += cg -> createCut (cs, sin (ub) + ub, -1, w_ind, 1., x_ind,  1.); // up: w + x <= f ub + ub
143    }
144    else {
145      ncuts += cg -> createCut (cs, cos (ub) - ub, +1, w_ind, 1., x_ind, -1.); // dn: w - x >= f ub - ub
146      ncuts += cg -> createCut (cs, cos (ub) + ub, -1, w_ind, 1., x_ind,  1.); // up: w + x <= f ub + ub
147    }
148  }
149
150  return ncuts;
151}
152
153
154/// normalize angle within [0,b] (typically, pi or 2pi)
155inline CouNumber modulo (register CouNumber a, register CouNumber b)
156  {return a - b * floor (a/b);}
157
158
159/// restrict to quarter of the interval [0,2pi]
160int bayEnvelope (const CouenneCutGenerator *, OsiCuts &, int, int, 
161                 CouNumber, CouNumber, CouNumber, bool &, bool &);
162
163
164/// real linearization of sine/cosine
165
166int trigEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
167                   OsiCuts &cs,                   // cut set to be enriched
168                   exprAux *w,
169                   expression *arg,
170                   enum cou_trig which_trig) {
171
172  expression *lbe, *ube;
173
174  arg -> getBounds (lbe, ube);
175
176  CouNumber lb = (*lbe) (), 
177            ub = (*ube) (),
178            // if cosine, scale variables to pretend this is a sine problem
179            displ = (which_trig == COU_COSINE) ? M_PI_2 : 0.;
180
181  delete lbe;
182  delete ube;
183
184  int ncuts = 0,
185    xi = arg -> Index (),
186    wi = w   -> Index ();
187
188  if (fabs (ub - lb) < COUENNE_EPS) {
189
190    CouNumber x0 = 0.5 * (ub+lb), f, fp;
191
192    if (which_trig == COU_SINE) {f = sin (x0); fp =  cos (x0);}
193    else                        {f = cos (x0); fp = -sin (x0);}
194
195    return cg -> createCut (cs, f - fp*x0, 0, wi, 1., xi, -fp);
196  }
197
198  // true if, in the first call (lb), a lower/upper chord was added
199  // --> no such chord must be generated in the second call (ub)
200  bool skip_up = false, 
201       skip_dn = false;
202
203  if (lb > -COUENNE_INFINITY) ncuts += bayEnvelope (cg, cs, wi, xi, lb, ub, displ, skip_up, skip_dn);
204  if (ub <  COUENNE_INFINITY) ncuts += bayEnvelope (cg, cs, wi, xi, ub, lb, displ, skip_up, skip_dn);
205
206  return ncuts;
207}
208
209
210//                             __
211// study single bay ( \__/ or /  \ ) of the trigonometric function
212//
213
214int bayEnvelope (const CouenneCutGenerator *cg, // cut generator that has called us
215                 OsiCuts &cs,                   // cut set to be enriched
216                 int wi,
217                 int xi,
218                 CouNumber x0,           // starting point
219                 CouNumber x1,           // other bound
220                 CouNumber displacement, 
221                 bool &skip_up, 
222                 bool &skip_dn) {
223
224  CouNumber tpt,
225    rx0  = modulo (x0 + displacement, 2*M_PI),
226    rx1  = rx0 + x1 - x0,
227    base = x0 - rx0,
228    sinrx0 = sin (rx0), zero;
229
230  int ncuts = 0,
231    up   = (rx0 < M_PI) ? +1 : -1,
232    left = (x0  < x1)   ? +1 : -1;
233
234  // starting point of the current bay
235  zero = (up>0) ? 0. : M_PI;
236
237  bool *s0, *s1;
238
239  if (up>0) {s0 = &skip_up; s1 = &skip_dn;}
240  else      {s0 = &skip_dn; s1 = &skip_up;}
241
242  if (left * (modulo (rx0, M_PI) - M_PI_2) < 0) { 
243
244    // after  flex (i.e., at \_ or /~ ) for left  bound,
245    // before flex (i.e., at _/ or ~\ ) for right bound
246
247    // out of the "belly": tangent. If on upper bay we consider the
248    // lower half-plane, and viceversa --> use -up
249    ncuts += cg -> addTangent (cs, wi, xi, x0, sin (rx0), cos (rx0), -up);
250
251    // leftmost extreme to search for tangent point
252    CouNumber extr0 = .75 * M_PI * (left+1) - M_PI_2 * up; 
253
254    // in:
255    if ((left * (rx1 - M_PI * ((left - up) / 2 + 1)) <= 0) ||   // if rx1 in same "belly", or
256        (left * (rx1 - (tpt = trigNewton
257                        (rx0, extr0, extr0 + M_PI_2))) <= 0)) { // before closest leaning point
258      if (!*s1) // -> chord, if not already added in previous call
259        *s1 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1,  sin (rx1), up)) > 0);
260    } else      ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
261  }
262  else {
263
264    // after  stationary point (i.e., _/ or ~\ ) for left  bound,
265    // before stationary point (i.e., /~ or \_ ) for right bound
266 
267    //    if (left * (rx1 - left * (zero + 5*M_PI_2)) < 0) {
268    if (left * (rx1 - (4*left - up + 2) * M_PI_2) < 0) {
269      CouNumber cosrx0 = cos (rx0);
270      if (up * (sin (rx1) - sinrx0 - cosrx0 * (rx1-rx0)) < 0) 
271        // (b,sinb) below tangent --> tangent
272        ncuts += cg -> addTangent (cs, wi, xi, x0, sinrx0, cosrx0, -up);
273      else {    // up: either chord or leaning plane
274        CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
275        tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
276        if (left * (rx1 - tpt) < 0) {
277          if (!*s0)
278            *s0 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), -up)) > 0);
279        } else      ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
280      }
281    } else {
282      CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
283      tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
284      ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
285    }
286
287    // down: other chord or leaning plane
288    if ((left * (rx1 - (zero + M_PI)) < 0) || 
289        (left * (rx1 - (tpt = trigNewton (rx0, (2 +   left - up) * M_PI_2, 
290                                               (2 + 2*left - up) * M_PI_2))) < 0)) {
291      if (!*s1) 
292        *s1 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), up)) > 0);
293    } else      ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
294  }
295
296  return ncuts;
297}
Note: See TracBrowser for help on using the repository browser.