00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <math.h>
00012 #ifndef M_PI
00013 # define M_PI 3.14159265358979323846
00014 #endif
00015 #ifndef M_PI_2
00016 # define M_PI_2 1.57079632679489661923
00017 #endif
00018
00019 #include "OsiSolverInterface.hpp"
00020 #include "CouenneTypes.hpp"
00021 #include "CouenneCutGenerator.hpp"
00022 #include "exprSin.hpp"
00023 #include "exprCos.hpp"
00024 #include "exprAux.hpp"
00025
00026 #define NEW_TRIG
00027
00028 #ifndef NEW_TRIG
00030 int addHexagon (const CouenneCutGenerator *, // cut generator that has called us
00031 OsiCuts &,
00032 enum cou_trig,
00033 expression *,
00034 expression *);
00035 #endif
00036
00038 int trigEnvelope (const CouenneCutGenerator *, OsiCuts &,
00039 expression *, expression *, enum cou_trig);
00040
00041
00043
00044 void exprSin::generateCuts (expression *w,
00045 OsiCuts &cs, const CouenneCutGenerator *cg,
00046 t_chg_bounds *chg, int wind,
00047 CouNumber lbw, CouNumber ubw) {
00048
00049
00050
00051
00052
00053
00054
00055
00056 #ifdef NEW_TRIG
00057 if (trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_SINE) == 0)
00058 #else
00059 if (addHexagon (cg, cs, COU_SINE, w, w -> Image () -> Argument()) == 0)
00060 #endif
00061 {
00062
00063 }
00064 }
00065
00066
00068
00069 void exprCos::generateCuts (expression *w,
00070 OsiCuts &cs, const CouenneCutGenerator *cg,
00071 t_chg_bounds *chg, int wind,
00072 CouNumber lbw, CouNumber ubw) {
00073
00074
00075
00076
00077
00078
00079
00080
00081 #ifdef NEW_TRIG
00082 if (trigEnvelope (cg, cs, w, w -> Image () -> Argument (), COU_COSINE) == 0)
00083 #else
00084 if (addHexagon (cg, cs, COU_COSINE, w, w -> Image () -> Argument()) == 0)
00085 #endif
00086 {
00087
00088 }
00089 }
00090
00091
00093 int bayEnvelope (const CouenneCutGenerator *, OsiCuts &, int, int,
00094 CouNumber, CouNumber, CouNumber, bool &, bool &);
00095
00096
00098
00099 int trigEnvelope (const CouenneCutGenerator *cg,
00100 OsiCuts &cs,
00101 expression *w,
00102 expression *arg,
00103 enum cou_trig which_trig) {
00104
00105 CouNumber lb, ub;
00106 arg -> getBounds (lb, ub);
00107
00108
00109 CouNumber displ = (which_trig == COU_COSINE) ? M_PI_2 : 0.;
00110
00111 int ncuts = 0,
00112 xi = arg -> Index (),
00113 wi = w -> Index ();
00114
00115 if (fabs (ub - lb) < COUENNE_EPS) {
00116
00117 CouNumber x0 = 0.5 * (ub+lb), f, fp;
00118
00119 if (which_trig == COU_SINE) {f = sin (x0); fp = cos (x0);}
00120 else {f = cos (x0); fp = -sin (x0);}
00121
00122 return cg -> createCut (cs, f - fp*x0, 0, wi, 1., xi, -fp);
00123 }
00124
00125
00126
00127 bool skip_up = false,
00128 skip_dn = false;
00129
00130 if (lb > -COUENNE_INFINITY) ncuts += bayEnvelope (cg, cs, wi, xi, lb, ub, displ, skip_up, skip_dn);
00131 if (ub < COUENNE_INFINITY) ncuts += bayEnvelope (cg, cs, wi, xi, ub, lb, displ, skip_up, skip_dn);
00132
00133 return ncuts;
00134 }
00135
00136
00137
00138
00139
00140
00141 int bayEnvelope (const CouenneCutGenerator *cg,
00142 OsiCuts &cs,
00143 int wi,
00144 int xi,
00145 CouNumber x0,
00146 CouNumber x1,
00147 CouNumber displacement,
00148 bool &skip_up,
00149 bool &skip_dn) {
00150
00151 CouNumber tpt,
00152 rx0 = modulo (x0 + displacement, 2*M_PI),
00153 rx1 = rx0 + x1 - x0,
00154 base = x0 - rx0,
00155 sinrx0 = sin (rx0), zero;
00156
00157 int ncuts = 0,
00158 up = (rx0 < M_PI) ? +1 : -1,
00159 left = (x0 < x1) ? +1 : -1;
00160
00161
00162 zero = (up>0) ? 0. : M_PI;
00163
00164 bool *s0, *s1;
00165
00166 if (up>0) {s0 = &skip_up; s1 = &skip_dn;}
00167 else {s0 = &skip_dn; s1 = &skip_up;}
00168
00169 if (left * (modulo (rx0, M_PI) - M_PI_2) < 0) {
00170
00171
00172
00173
00174
00175
00176 ncuts += cg -> addTangent (cs, wi, xi, x0, sin (rx0), cos (rx0), -up);
00177
00178
00179 CouNumber extr0 = .75 * M_PI * (left+1) - M_PI_2 * up;
00180
00181
00182 if ((left * (rx1 - M_PI * ((left - up) / 2 + 1)) <= 0) ||
00183 (left * (rx1 - (tpt = trigNewton
00184 (rx0, extr0, extr0 + M_PI_2))) <= 0)) {
00185 if (!*s1)
00186 *s1 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), up)) > 0);
00187 } else ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base+tpt, sin (tpt), up);
00188 }
00189 else {
00190
00191
00192
00193
00194
00195 if (left * (rx1 - (4*left - up + 2) * M_PI_2) < 0) {
00196 CouNumber cosrx0 = cos (rx0);
00197 if (up * (sin (rx1) - sinrx0 - cosrx0 * (rx1-rx0)) < 0)
00198
00199 ncuts += cg -> addTangent (cs, wi, xi, x0, sinrx0, cosrx0, -up);
00200 else {
00201 CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
00202 tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
00203 if (left * (rx1 - tpt) < 0) {
00204 if (!*s0)
00205 *s0 = ((ncuts += cg->addSegment (cs, wi, xi, x0, sin(rx0), x1, sin(rx1), -up)) > 0);
00206 } else ncuts += cg->addSegment (cs, wi, xi, x0, sin(rx0), base+tpt, sin(tpt), -up);
00207 }
00208 } else {
00209 CouNumber searchpt = M_PI_2 * (2 + 3*left - up);
00210 tpt = trigNewton (rx0, searchpt, searchpt + left * M_PI_2);
00211 ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), -up);
00212 }
00213
00214
00215 if ((left * (rx1 - (zero + M_PI)) < 0) ||
00216 (left * (rx1 - (tpt = trigNewton (rx0, (2 + left - up) * M_PI_2,
00217 (2 + 2*left - up) * M_PI_2))) < 0)) {
00218 if (!*s1)
00219 *s1 = ((ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), x1, sin (rx1), up)) > 0);
00220 } else ncuts += cg -> addSegment (cs, wi, xi, x0, sin (rx0), base + tpt, sin (tpt), up);
00221 }
00222
00223 return ncuts;
00224 }
00225
00226
00227 #ifndef NEW_TRIG
00228
00229
00231
00232 int addHexagon (const CouenneCutGenerator *cg,
00233 OsiCuts &cs,
00234 enum cou_trig tt,
00235 expression *aux,
00236 expression *arg) {
00237
00238
00239 CouNumber lb, ub;
00240 arg -> getBounds (lb, ub);
00241
00242 int ncuts = 0,
00243 x_ind = arg -> Index (),
00244 w_ind = aux -> Index ();
00245
00246 if (fabs (ub - lb) < COUENNE_EPS) {
00247
00248 CouNumber x0 = 0.5 * (ub+lb), f, fp;
00249
00250 if (tt == COU_SINE) {f = sin (x0); fp = cos (x0);}
00251 else {f = cos (x0); fp = -sin (x0);}
00252
00253 return cg -> createCut (cs, f - fp*x0, 0, w_ind, 1., x_ind, -fp);
00254 }
00255
00256
00257
00258
00259
00260 if (lb > -COUENNE_INFINITY) {
00261 if (tt == COU_SINE) {
00262 ncuts += cg -> createCut (cs, sin (lb) - lb, -1, w_ind, 1., x_ind, -1.);
00263 ncuts += cg -> createCut (cs, sin (lb) + lb, +1, w_ind, 1., x_ind, 1.);
00264 }
00265 else {
00266 ncuts += cg -> createCut (cs, cos (lb) - lb, -1, w_ind, 1., x_ind, -1.);
00267 ncuts += cg -> createCut (cs, cos (lb) + lb, +1, w_ind, 1., x_ind, 1.);
00268 }
00269 }
00270
00271
00272 if (ub < COUENNE_INFINITY) {
00273 if (tt == COU_SINE) {
00274 ncuts += cg -> createCut (cs, sin (ub) - ub, +1, w_ind, 1., x_ind, -1.);
00275 ncuts += cg -> createCut (cs, sin (ub) + ub, -1, w_ind, 1., x_ind, 1.);
00276 }
00277 else {
00278 ncuts += cg -> createCut (cs, cos (ub) - ub, +1, w_ind, 1., x_ind, -1.);
00279 ncuts += cg -> createCut (cs, cos (ub) + ub, -1, w_ind, 1., x_ind, 1.);
00280 }
00281 }
00282
00283 return ncuts;
00284 }
00285 #endif