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