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
00038 #ifndef MAT_INTERVAL
00039 #define MAT_INTERVAL
00040 #include <math.h>
00041 #include <iomanip>
00042 namespace mat {
00043
00044
00045 template<typename Treal>
00046 class Interval {
00047 public:
00048 explicit Interval(Treal low = 1, Treal upp = -1)
00049 : lowerBound(low), upperBound(upp) {
00050 }
00051 inline bool empty() const {return lowerBound > upperBound;}
00052
00053 static Interval intersect(Interval const & A, Interval const & B) {
00054 if (A.empty())
00055 return A;
00056 else if (B.empty())
00057 return B;
00058 else
00059 return Interval(A.lowerBound > B.lowerBound ?
00060 A.lowerBound : B.lowerBound,
00061 A.upperBound < B.upperBound ?
00062 A.upperBound : B.upperBound);
00063 }
00064 inline void intersect(Interval const & other) {
00065 if (this->empty())
00066 return;
00067 else if (other.empty()) {
00068 *this = other;
00069 return;
00070 }
00071 else {
00072 this->lowerBound = this->lowerBound > other.lowerBound ?
00073 this->lowerBound : other.lowerBound;
00074 this->upperBound = this->upperBound < other.upperBound ?
00075 this->upperBound : other.upperBound;
00076 return;
00077 }
00078 }
00079
00080 inline void intersect_always_non_empty(Interval const & other) {
00081 if (this->empty()) {
00082 *this = other;
00083 return;
00084 }
00085 if (other.empty())
00086 return;
00087
00088 Interval intersection = intersect(*this, other);
00089 if ( !intersection.empty() ) {
00090 *this = intersection;
00091 return;
00092 }
00093
00094 Treal tmp_val;
00095 if ( this->lowerBound > other.upperBound )
00096 tmp_val = (this->lowerBound + other.upperBound) / 2;
00097 else if ( this->upperBound < other.lowerBound )
00098 tmp_val = (this->upperBound + other.lowerBound) / 2;
00099 else
00100 assert(0);
00101 this->lowerBound = tmp_val;
00102 this->upperBound = tmp_val;
00103 return;
00104 }
00105
00109 inline Treal length() const {
00110 if (empty())
00111 return 0.0;
00112 else
00113 return upperBound - lowerBound;
00114 }
00115 inline Treal midPoint() const {
00116 assert(!empty());
00117 return (upperBound + lowerBound) / 2;
00118 }
00119 inline bool cover(Treal const value) const {
00120 if (empty())
00121 return false;
00122 else
00123 return (value <= upperBound && value >= lowerBound);
00124 }
00125 inline bool overlap(Interval const & other) const {
00126 Interval tmp = intersect(*this, other);
00127 return !tmp.empty();
00128 }
00129
00133 inline void increase(Treal const value) {
00134 if (empty())
00135 throw Failure("Interval<Treal>::increase(Treal const) : "
00136 "Attempt to increase empty interval.");
00137 lowerBound -= value;
00138 upperBound += value;
00139 }
00140 inline void decrease(Treal const value) {
00141 lowerBound += value;
00142 upperBound -= value;
00143 }
00144 inline Treal low() const {return lowerBound;}
00145 inline Treal upp() const {return upperBound;}
00146
00147 #if 0
00148 inline Interval<Treal> operator*(Interval<Treal> const & other) const {
00149 assert(lowerBound >= 0);
00150 assert(other.lowerBound >= 0);
00151 return Interval<Treal>(lowerBound * other.lowerBound,
00152 upperBound * other.upperBound);
00153 }
00154 #endif
00155
00156 inline Interval<Treal> operator*(Treal const & value) const {
00157 if (value < 0)
00158 return Interval<Treal>(upperBound * value, lowerBound * value);
00159 else
00160 return Interval<Treal>(lowerBound * value, upperBound * value);
00161 }
00162
00163
00164 inline Interval<Treal> operator-(Interval<Treal> const & other) const {
00165 return *this + (-1.0 * other);
00166
00167
00168 }
00169
00170 inline Interval<Treal> operator+(Interval<Treal> const & other) const {
00171 return Interval<Treal>(lowerBound + other.lowerBound,
00172 upperBound + other.upperBound);
00173 }
00174 inline Interval<Treal> operator/(Treal const & value) const {
00175 if (value < 0)
00176 return Interval<Treal>(upperBound / value, lowerBound / value);
00177 else
00178 return Interval<Treal>(lowerBound / value, upperBound / value);
00179 }
00180 inline Interval<Treal> operator-(Treal const & value) const {
00181 return Interval<Treal>(lowerBound - value, upperBound - value);
00182 }
00183 inline Interval<Treal> operator+(Treal const & value) const {
00184 return Interval<Treal>(lowerBound + value, upperBound + value);
00185 }
00186
00187
00188
00189 void puriStep(int poly);
00190 void invPuriStep(int poly);
00191
00192
00193 void puriStep(int poly, Treal alpha);
00194 void invPuriStep(int poly, Treal alpha);
00195 protected:
00196 Treal lowerBound;
00197 Treal upperBound;
00198 private:
00199 };
00200
00201 #if 0
00202
00203 template<typename Treal>
00204 inline Interval<Treal> operator-(Treal const value,
00205 Interval<Treal> const & other) {
00206 return Interval<Treal>(value - other.lowerBound,
00207 value - other.upperBound);
00208 }
00209 template<typename Treal>
00210 inline Interval<Treal> operator+(Treal const value,
00211 Interval<Treal> const & other) {
00212 return Interval<Treal>(value + other.lowerBound,
00213 value + other.upperBound);
00214 }
00215 #endif
00216
00217 #if 1
00218 template<typename Treal>
00219 Interval<Treal> sqrtInt(Interval<Treal> const & other) {
00220 if (other.empty())
00221 return other;
00222 else {
00223 assert(other.low() >= 0);
00224 return Interval<Treal>(template_blas_sqrt(other.low()), template_blas_sqrt(other.upp()));
00225 }
00226 }
00227 #endif
00228
00229 template<typename Treal>
00230 void Interval<Treal>::puriStep(int poly) {
00231 if (upperBound > 2.0 || lowerBound < -1.0)
00232 throw Failure("Interval<Treal>::puriStep(int) : It is assumed here "
00233 "that the interval I is within [-1.0, 2.0]");
00234 Interval<Treal> oneInt(-1.0,1.0);
00235 Interval<Treal> zeroInt(0.0,2.0);
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 bool nonEmptyIntervalInZeroToOne = false;
00248 if(upperBound > lowerBound && lowerBound >= 0 && upperBound <= 1)
00249 nonEmptyIntervalInZeroToOne = true;
00250 if (poly) {
00251
00252 *this = Interval<Treal>::intersect(*this,oneInt);
00253
00254 upperBound = 2 * upperBound - upperBound * upperBound;
00255 lowerBound = 2 * lowerBound - lowerBound * lowerBound;
00256 if(nonEmptyIntervalInZeroToOne && upperBound < lowerBound) {
00257
00258 Treal midPoint = (lowerBound + upperBound) / 2;
00259 upperBound = lowerBound = midPoint;
00260 }
00261 }
00262 else {
00263 *this = Interval<Treal>::intersect(*this,zeroInt);
00264
00265 upperBound = upperBound * upperBound;
00266 lowerBound = lowerBound * lowerBound;
00267 }
00268 }
00269
00270 template<typename Treal>
00271 void Interval<Treal>::invPuriStep(int poly) {
00272 if (upperBound > 2.0 || lowerBound < -1.0)
00273 throw Failure("Interval<Treal>::puriStep(int) : It is assumed here "
00274 "that the interval I is within [-1.0, 1.0]");
00275 Interval<Treal> oneInt(-1.0,1.0);
00276 Interval<Treal> zeroInt(0.0,2.0);
00277 if (poly) {
00278
00279 *this = Interval<Treal>::intersect(*this,oneInt);
00280
00281 upperBound = 1.0 - template_blas_sqrt(1.0 - upperBound);
00282 lowerBound = 1.0 - template_blas_sqrt(1.0 - lowerBound);
00283 }
00284 else {
00285 *this = Interval<Treal>::intersect(*this,zeroInt);
00286
00287 upperBound = template_blas_sqrt(upperBound);
00288 lowerBound = template_blas_sqrt(lowerBound);
00289 }
00290 }
00291
00292 #if 1
00293 template<typename Treal>
00294 void Interval<Treal>::puriStep(int poly, Treal alpha) {
00295 if (upperBound > 2.0 || lowerBound < -1.0)
00296 throw Failure("Interval<Treal>::puriStep(int, real) : It is assumed here "
00297 "that the interval I is within [-1.0, 2.0]");
00298 Interval<Treal> oneInt(-1.0,1.0);
00299 Interval<Treal> zeroInt(0.0,2.0);
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311 bool nonEmptyIntervalInZeroToOne = false;
00312 if(upperBound > lowerBound && lowerBound >= 0 && upperBound <= 1)
00313 nonEmptyIntervalInZeroToOne = true;
00314 if (poly) {
00315
00316 *this = Interval<Treal>::intersect(*this,oneInt);
00317
00318 upperBound = 2 * alpha * upperBound - alpha * alpha * upperBound * upperBound;
00319 lowerBound = 2 * alpha * lowerBound - alpha * alpha * lowerBound * lowerBound;
00320 if(nonEmptyIntervalInZeroToOne && upperBound < lowerBound) {
00321
00322 Treal midPoint = (lowerBound + upperBound) / 2;
00323 upperBound = lowerBound = midPoint;
00324 }
00325 }
00326 else {
00327 *this = Interval<Treal>::intersect(*this,zeroInt);
00328
00329 upperBound = ( alpha * upperBound + (1-alpha) ) * ( alpha * upperBound + (1-alpha) );
00330 lowerBound = ( alpha * lowerBound + (1-alpha) ) * ( alpha * lowerBound + (1-alpha) );
00331 }
00332 }
00333
00334 template<typename Treal>
00335 void Interval<Treal>::invPuriStep(int poly, Treal alpha) {
00336 if (upperBound > 2.0 || lowerBound < -1.0)
00337 throw Failure("Interval<Treal>::invPuriStep(int, real) : It is assumed here "
00338 "that the interval I is within [-1.0, 2.0]");
00339 Interval<Treal> oneInt(-1.0,1.0);
00340 Interval<Treal> zeroInt(0.0,2.0);
00341 if (poly) {
00342
00343 *this = Interval<Treal>::intersect(*this,oneInt);
00344
00345 upperBound = ( 1.0 - template_blas_sqrt(1.0 - upperBound) ) / alpha;
00346 lowerBound = ( 1.0 - template_blas_sqrt(1.0 - lowerBound) ) / alpha;
00347 }
00348 else {
00349 *this = Interval<Treal>::intersect(*this,zeroInt);
00350
00351 upperBound = ( template_blas_sqrt(upperBound) - (1-alpha) ) / alpha;
00352 lowerBound = ( template_blas_sqrt(lowerBound) - (1-alpha) ) / alpha;
00353 }
00354 }
00355
00356 #endif
00357
00358 template<typename Treal>
00359 std::ostream& operator<<(std::ostream& s,
00360 Interval<Treal> const & in) {
00361 if (in.empty())
00362 s<<"[empty]";
00363 else
00364 s<<"["<<in.low()<<", "<<in.upp()<<"]";
00365 return s;
00366 }
00367
00368 }
00369 #endif