1 /* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */ 2 /* 3 * Main authors: 4 * Christian Schulte <schulte@gecode.org> 5 * Vincent Barichard <Vincent.Barichard@univ-angers.fr> 6 * 7 * Copyright: 8 * Christian Schulte, 2002 9 * Vincent Barichard, 2012 10 * 11 * This file is part of Gecode, the generic constraint 12 * development environment: 13 * http://www.gecode.org 14 * 15 * Permission is hereby granted, free of charge, to any person obtaining 16 * a copy of this software and associated documentation files (the 17 * "Software"), to deal in the Software without restriction, including 18 * without limitation the rights to use, copy, modify, merge, publish, 19 * distribute, sublicense, and/or sell copies of the Software, and to 20 * permit persons to whom the Software is furnished to do so, subject to 21 * the following conditions: 22 * 23 * The above copyright notice and this permission notice shall be 24 * included in all copies or substantial portions of the Software. 25 * 26 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 27 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 28 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 29 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 30 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 31 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 32 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 33 * 34 */ 35 36 #include <algorithm> 37 #include <climits> 38 39 #include <gecode/float/linear.hh> 40 #include <gecode/float/rel.hh> 41 42 namespace Gecode { namespace Float { namespace Linear { 43 44 void estimate(Term * t,int n,FloatVal c,FloatNum & l,FloatNum & u)45 estimate(Term* t, int n, FloatVal c, FloatNum& l, FloatNum &u) { 46 FloatVal est = c; 47 for (int i=n; i--; ) 48 est += t[i].a * t[i].x.domain(); 49 FloatNum min = est.min(); 50 FloatNum max = est.max(); 51 if (min < Limits::min) 52 min = Limits::min; 53 if (min > Limits::max) 54 min = Limits::max; 55 l = min; 56 if (max < Limits::min) 57 max = Limits::min; 58 if (max > Limits::max) 59 max = Limits::max; 60 u = max; 61 } 62 63 forceinline bool overflow(Term * t,int n,FloatVal c)64 overflow(Term* t, int n, FloatVal c) { 65 FloatVal est = c; 66 for (int i=n; i--; ) 67 est += t[i].a * t[i].x.domain(); 68 FloatNum min = est.min(); 69 FloatNum max = est.max(); 70 return ((min < Limits::min) || (min > Limits::max) || 71 (max < Limits::min) || (max > Limits::max)); 72 } 73 74 /// Sort linear terms by view 75 class TermLess { 76 public: 77 forceinline bool operator ()(const Term & a,const Term & b)78 operator ()(const Term& a, const Term& b) { 79 return a.x < b.x; 80 } 81 }; 82 83 /// Extend terms by adding view for result 84 FloatView extend(Home home,Region & r,Term * & t,int & n)85 extend(Home home, Region& r, Term*& t, int& n) { 86 FloatNum min, max; 87 estimate(t,n,0.0,min,max); 88 FloatVar x(home,min,max); 89 Term* et = r.alloc<Term>(n+1); 90 for (int i=n; i--; ) 91 et[i] = t[i]; 92 et[n].a=-1.0; et[n].x=x; 93 t=et; n++; 94 return x; 95 } 96 97 98 /** 99 * \brief Posting n-ary propagators 100 * 101 */ 102 template<class View> 103 forceinline void post_nary(Home home,ViewArray<View> & x,ViewArray<View> & y,FloatRelType frt,FloatVal c)104 post_nary(Home home, 105 ViewArray<View>& x, ViewArray<View>& y, FloatRelType frt, 106 FloatVal c) { 107 switch (frt) { 108 case FRT_EQ: 109 GECODE_ES_FAIL((Eq<View,View >::post(home,x,y,c))); 110 break; 111 case FRT_LQ: 112 GECODE_ES_FAIL((Lq<View,View >::post(home,x,y,c))); 113 break; 114 default: GECODE_NEVER; 115 } 116 } 117 118 void dopost(Home home,Term * t,int n,FloatRelType frt,FloatVal c)119 dopost(Home home, Term* t, int n, FloatRelType frt, FloatVal c) { 120 Limits::check(c,"Float::linear"); 121 122 for (int i=n; i--; ) { 123 if ((t[i].a.min() < 0.0) && (t[i].a.max() > 0.0)) 124 throw ValueMixedSign("Float::linear[coefficient]"); 125 if (t[i].x.assigned()) { 126 c -= t[i].a * t[i].x.val(); 127 t[i]=t[--n]; 128 } 129 } 130 131 if ((c < Limits::min) || (c > Limits::max) || overflow(t, n, c)) 132 throw OutOfLimits("Float::linear"); 133 134 /* 135 * Join coefficients for aliased variables: 136 * 137 */ 138 { 139 // Group same variables 140 TermLess tl; 141 Support::quicksort<Term,TermLess>(t,n,tl); 142 143 // Join adjacent variables 144 int i = 0; 145 int j = 0; 146 while (i < n) { 147 Limits::check(t[i].a,"Float::linear"); 148 FloatVal a = t[i].a; 149 FloatView x = t[i].x; 150 while ((++i < n) && (t[i].x == x)) { 151 a += t[i].a; 152 Limits::check(a,"Float::linear"); 153 } 154 if (a != 0.0) { 155 t[j].a = a; t[j].x = x; j++; 156 } 157 } 158 n = j; 159 } 160 161 Term *t_p, *t_n; 162 int n_p, n_n; 163 164 /* 165 * Partition into positive/negative coefficents 166 * 167 */ 168 if (n > 0) { 169 int i = 0; 170 int j = n-1; 171 while (true) { 172 while ((t[j].a < 0) && (--j >= 0)) ; 173 while ((t[i].a > 0) && (++i < n)) ; 174 if (j <= i) break; 175 std::swap(t[i],t[j]); 176 } 177 t_p = t; n_p = i; 178 t_n = t+n_p; n_n = n-n_p; 179 } else { 180 t_p = t; n_p = 0; 181 t_n = t; n_n = 0; 182 } 183 184 /* 185 * Make all coefficients positive 186 * 187 */ 188 for (int i=n_n; i--; ) 189 t_n[i].a = -t_n[i].a; 190 191 if (frt == FRT_GQ) { 192 frt = FRT_LQ; 193 std::swap(n_p,n_n); std::swap(t_p,t_n); c = -c; 194 } 195 196 if (n == 0) { 197 switch (frt) { 198 case FRT_EQ: if (!c.in(0.0)) home.fail(); break; 199 case FRT_LQ: if (c.max() < 0.0) home.fail(); break; 200 default: GECODE_NEVER; 201 } 202 return; 203 } 204 205 /* 206 * Test for unit coefficients only 207 * 208 */ 209 bool is_unit = true; 210 for (int i=n; i--; ) 211 if (!(t[i].a == 1.0)) { 212 is_unit = false; 213 break; 214 } 215 216 if (is_unit) { 217 // Unit coefficients 218 ViewArray<FloatView> x(home,n_p); 219 for (int i = n_p; i--; ) 220 x[i] = t_p[i].x; 221 ViewArray<FloatView> y(home,n_n); 222 for (int i = n_n; i--; ) 223 y[i] = t_n[i].x; 224 post_nary<FloatView>(home,x,y,frt,c); 225 } else { 226 // Arbitrary coefficients 227 ViewArray<ScaleView> x(home,n_p); 228 for (int i = n_p; i--; ) 229 x[i] = ScaleView(t_p[i].a,t_p[i].x); 230 ViewArray<ScaleView> y(home,n_n); 231 for (int i = n_n; i--; ) 232 y[i] = ScaleView(t_n[i].a,t_n[i].x); 233 post_nary<ScaleView>(home,x,y,frt,c); 234 } 235 } 236 237 void post(Home home,Term * t,int n,FloatRelType frt,FloatVal c)238 post(Home home, Term* t, int n, FloatRelType frt, FloatVal c) { 239 Region re; 240 switch (frt) { 241 case FRT_EQ: case FRT_LQ: case FRT_GQ: 242 break; 243 case FRT_NQ: case FRT_LE: case FRT_GR: 244 rel(home, extend(home,re,t,n), frt, c); 245 frt=FRT_EQ; c=0.0; 246 break; 247 default: 248 throw UnknownRelation("Float::linear"); 249 } 250 dopost(home, t, n, frt, c); 251 } 252 253 void post(Home home,Term * t,int n,FloatRelType frt,FloatVal c,Reify r)254 post(Home home, Term* t, int n, FloatRelType frt, FloatVal c, Reify r) { 255 Region re; 256 rel(home, extend(home,re,t,n), frt, c, r); 257 dopost(home, t, n, FRT_EQ, 0.0); 258 } 259 260 }}} 261 262 // STATISTICS: float-post 263 264