00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <math.h>
00012
00013 #include "CouenneExprSin.hpp"
00014 #include "CouenneExprClone.hpp"
00015 #include "CouenneExprCos.hpp"
00016 #include "CouenneExprBSin.hpp"
00017 #include "CouenneExprMul.hpp"
00018
00019 namespace Couenne {
00020
00021 static const CouNumber
00022 pi = M_PI,
00023 pi2 = M_PI * 2.,
00024 pih = M_PI / 2.;
00025
00026
00027
00028 expression *exprSin::differentiate (int index) {
00029
00030 return new exprMul (new exprCos (new exprClone (argument_)),
00031 argument_ -> differentiate (index));
00032 }
00033
00034
00035
00036 void exprSin::getBounds (expression *&lb, expression *&ub) {
00037
00038 expression *xl, *xu;
00039
00040 argument_ -> getBounds (xl, xu);
00041
00042 lb = new exprLBSin (xl, xu);
00043 ub = new exprUBSin (new exprClone (xl), new exprClone (xu));
00044 }
00045
00046
00047 void exprSin::getBounds (CouNumber &lb, CouNumber &ub) {
00048
00049 CouNumber l, u;
00050 argument_ -> getBounds (l, u);
00051
00052 if ((u - l >= pi2) ||
00053 (floor (l/pi2 - 0.75) <
00054 floor (u/pi2 - 0.75)))
00055 lb = -1.;
00056 else lb = CoinMin (sin (l), sin (u));
00057
00058 if ((u - l >= pi2) ||
00059 (floor (l/pi2 - 0.25) <
00060 floor (u/pi2 - 0.25)))
00061 ub = 1.;
00062 else ub = CoinMax (sin (l), sin (u));
00063 }
00064
00066 bool trigImpliedBound (enum cou_trig type, int wind, int xind,
00067 CouNumber *l, CouNumber *u, t_chg_bounds *chg) {
00068
00069
00070
00071 CouNumber *xl = l + xind, wl = l [wind],
00072 *xu = u + xind, wu = u [wind];
00073
00074 bool tighter = false;
00075
00076 CouNumber fl, fu, iwl, iwu, displacement;
00077
00078 if (type == COU_SINE) {fl = sin (*xl); fu = sin (*xu); displacement = pih;}
00079 else {fl = cos (*xl); fu = cos (*xu); displacement = 0.;}
00080
00081 iwl = acos (wl);
00082 iwu = acos (wu);
00083
00084
00085
00086
00087
00088
00090
00091 if (wu < fl - COUENNE_EPS) {
00092
00093 CouNumber base = displacement + pi2 * floor ((*xl + pi - displacement) / pi2);
00094
00095
00096
00097 if (updateBound (-1, xl, base + iwu)) {
00098 tighter = true;
00099 chg [xind]. setLower (t_chg_bounds::CHANGED);
00100 }
00101 }
00102
00103 if (wu < fu - COUENNE_EPS) {
00104
00105 CouNumber base = displacement + pi2 * floor ((*xu + pi - displacement) / pi2);
00106
00107
00108
00109 if (updateBound (+1, xu, base - iwu)) {
00110 tighter = true;
00111 chg [xind]. setUpper (t_chg_bounds::CHANGED);
00112 }
00113 }
00114
00115 if (wl > fl + COUENNE_EPS) {
00116
00117 CouNumber base = displacement + pi2 * floor ((*xl - displacement) / pi2) + pi;
00118
00119
00120
00121 if (updateBound (-1, xl, base + pi - iwl)) {
00122 tighter = true;
00123 chg [xind]. setLower (t_chg_bounds::CHANGED);
00124 }
00125 }
00126
00127 if (wl > fu + COUENNE_EPS) {
00128
00129 CouNumber base = displacement + pi2 * floor ((*xu - displacement) / pi2) + pi;
00130
00131
00132
00133 if (updateBound (+1, xu, base - pi + iwl)) {
00134 tighter = true;
00135 chg [xind]. setUpper (t_chg_bounds::CHANGED);
00136 }
00137 }
00138
00139
00140
00141 return tighter;
00142 }
00143
00144
00146 void exprSin::closestFeasible (expression *varind, expression *vardep,
00147 CouNumber& left, CouNumber& right) const
00148 {
00149 CouNumber curr = (*varind)() - pih;
00150 int period = (int)(curr/pi2);
00151 CouNumber curr_noperiod = curr - pi2*period;
00152 CouNumber inv = acos((*vardep)());
00153
00154 if (curr_noperiod < inv) {
00155 left = pi2*period - inv;
00156 right = pi2*period + inv;
00157 }
00158 else if (curr_noperiod < pi2-inv) {
00159 left = pi2*period + inv;
00160 right = pi2*(period+1) - inv;
00161 }
00162 else {
00163 left = pi2*(period+1) - inv;
00164 right = pi2*(period+1) + inv;
00165 }
00166 left += pih;
00167 right += pih;
00168 }
00169
00170 }