00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <algorithm>
00037 #include <climits>
00038
00039 #include <gecode/float/linear.hh>
00040 #include <gecode/float/rel.hh>
00041
00042 namespace Gecode { namespace Float { namespace Linear {
00043
00044 void
00045 estimate(Term* t, int n, FloatVal c, FloatNum& l, FloatNum &u) {
00046 FloatVal est = c;
00047 for (int i=n; i--; )
00048 est += t[i].a * t[i].x.domain();
00049 FloatNum min = est.min();
00050 FloatNum max = est.max();
00051 if (min < Limits::min)
00052 min = Limits::min;
00053 if (min > Limits::max)
00054 min = Limits::max;
00055 l = min;
00056 if (max < Limits::min)
00057 max = Limits::min;
00058 if (max > Limits::max)
00059 max = Limits::max;
00060 u = max;
00061 }
00062
00063 forceinline bool
00064 overflow(Term* t, int n, FloatVal c) {
00065 FloatVal est = c;
00066 for (int i=n; i--; )
00067 est += t[i].a * t[i].x.domain();
00068 FloatNum min = est.min();
00069 FloatNum max = est.max();
00070 return ((min < Limits::min) || (min > Limits::max) ||
00071 (max < Limits::min) || (max > Limits::max));
00072 }
00073
00075 class TermLess {
00076 public:
00077 forceinline bool
00078 operator ()(const Term& a, const Term& b) {
00079 return before(a.x,b.x);
00080 }
00081 };
00082
00084 FloatView
00085 extend(Home home, Region& r, Term*& t, int& n) {
00086 FloatNum min, max;
00087 estimate(t,n,0.0,min,max);
00088 FloatVar x(home,min,max);
00089 Term* et = r.alloc<Term>(n+1);
00090 for (int i=n; i--; )
00091 et[i] = t[i];
00092 et[n].a=-1.0; et[n].x=x;
00093 t=et; n++;
00094 return x;
00095 }
00096
00097
00102 template<class View>
00103 forceinline void
00104 post_nary(Home home,
00105 ViewArray<View>& x, ViewArray<View>& y, FloatRelType frt,
00106 FloatVal c) {
00107 switch (frt) {
00108 case FRT_EQ:
00109 GECODE_ES_FAIL((Eq<View,View >::post(home,x,y,c)));
00110 break;
00111 case FRT_LQ:
00112 GECODE_ES_FAIL((Lq<View,View >::post(home,x,y,c)));
00113 break;
00114 default: GECODE_NEVER;
00115 }
00116 }
00117
00118 void
00119 dopost(Home home, Term* t, int n, FloatRelType frt, FloatVal c) {
00120 Limits::check(c,"Float::linear");
00121
00122 for (int i=n; i--; ) {
00123 if ((t[i].a.min() < 0.0) && (t[i].a.max() > 0.0))
00124 throw ValueMixedSign("Float::linear[coefficient]");
00125 if (t[i].x.assigned()) {
00126 c -= t[i].a * t[i].x.val();
00127 t[i]=t[--n];
00128 }
00129 }
00130
00131 if ((c < Limits::min) || (c > Limits::max) || overflow(t, n, c))
00132 throw OutOfLimits("Float::linear");
00133
00134
00135
00136
00137
00138 {
00139
00140 TermLess tl;
00141 Support::quicksort<Term,TermLess>(t,n,tl);
00142
00143
00144 int i = 0;
00145 int j = 0;
00146 while (i < n) {
00147 Limits::check(t[i].a,"Float::linear");
00148 FloatVal a = t[i].a;
00149 FloatView x = t[i].x;
00150 while ((++i < n) && same(t[i].x,x)) {
00151 a += t[i].a;
00152 Limits::check(a,"Float::linear");
00153 }
00154 if (a != 0.0) {
00155 t[j].a = a; t[j].x = x; j++;
00156 }
00157 }
00158 n = j;
00159 }
00160
00161 Term *t_p, *t_n;
00162 int n_p, n_n;
00163
00164
00165
00166
00167
00168 if (n > 0) {
00169 int i = 0;
00170 int j = n-1;
00171 while (true) {
00172 while ((t[j].a < 0) && (--j >= 0)) ;
00173 while ((t[i].a > 0) && (++i < n)) ;
00174 if (j <= i) break;
00175 std::swap(t[i],t[j]);
00176 }
00177 t_p = t; n_p = i;
00178 t_n = t+n_p; n_n = n-n_p;
00179 } else {
00180 t_p = t; n_p = 0;
00181 t_n = t; n_n = 0;
00182 }
00183
00184
00185
00186
00187
00188 for (int i=n_n; i--; )
00189 t_n[i].a = -t_n[i].a;
00190
00191 if (frt == FRT_GQ) {
00192 frt = FRT_LQ;
00193 std::swap(n_p,n_n); std::swap(t_p,t_n); c = -c;
00194 }
00195
00196 if (n == 0) {
00197 switch (frt) {
00198 case FRT_EQ: if (!c.in(0.0)) home.fail(); break;
00199 case FRT_LQ: if (c.max() < 0.0) home.fail(); break;
00200 default: GECODE_NEVER;
00201 }
00202 return;
00203 }
00204
00205
00206
00207
00208
00209 bool is_unit = true;
00210 for (int i=n; i--; )
00211 if (!(t[i].a == 1.0)) {
00212 is_unit = false;
00213 break;
00214 }
00215
00216 if (is_unit) {
00217
00218 ViewArray<FloatView> x(home,n_p);
00219 for (int i = n_p; i--; )
00220 x[i] = t_p[i].x;
00221 ViewArray<FloatView> y(home,n_n);
00222 for (int i = n_n; i--; )
00223 y[i] = t_n[i].x;
00224 post_nary<FloatView>(home,x,y,frt,c);
00225 } else {
00226
00227 ViewArray<ScaleView> x(home,n_p);
00228 for (int i = n_p; i--; )
00229 x[i] = ScaleView(t_p[i].a,t_p[i].x);
00230 ViewArray<ScaleView> y(home,n_n);
00231 for (int i = n_n; i--; )
00232 y[i] = ScaleView(t_n[i].a,t_n[i].x);
00233 post_nary<ScaleView>(home,x,y,frt,c);
00234 }
00235 }
00236
00237 void
00238 post(Home home, Term* t, int n, FloatRelType frt, FloatVal c) {
00239 Region re;
00240 switch (frt) {
00241 case FRT_EQ: case FRT_LQ: case FRT_GQ:
00242 break;
00243 case FRT_NQ: case FRT_LE: case FRT_GR:
00244 rel(home, extend(home,re,t,n), frt, c);
00245 frt=FRT_EQ; c=0.0;
00246 break;
00247 default:
00248 throw UnknownRelation("Float::linear");
00249 }
00250 dopost(home, t, n, frt, c);
00251 }
00252
00253 void
00254 post(Home home, Term* t, int n, FloatRelType frt, FloatVal c, Reify r) {
00255 Region re;
00256 rel(home, extend(home,re,t,n), frt, c, r);
00257 dopost(home, t, n, FRT_EQ, 0.0);
00258 }
00259
00260 }}}
00261
00262
00263