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