CLAM-Development
1.1
|
00001 /* 00002 * Copyright (c) 2004 MUSIC TECHNOLOGY GROUP (MTG) 00003 * UNIVERSITAT POMPEU FABRA 00004 * 00005 * 00006 * This program is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program; if not, write to the Free Software 00018 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00019 * 00020 */ 00021 00022 #ifndef _BasicOps_ 00023 #define _BasicOps_ 00024 00025 #include <numeric> 00026 #include <functional> 00027 #include "StaticBool.hxx" 00028 #include "DataTypes.hxx" 00029 #include "Array.hxx" 00030 #include "CLAM_Math.hxx" 00031 #include "Order.hxx" 00032 00033 00034 using std::accumulate; 00035 using std::inner_product; 00036 using std::mem_fun; 00037 00038 00051 namespace CLAM { 00052 00053 00054 00055 template <int o> struct Pow 00056 { 00057 public: 00058 template <class T> 00059 T operator () (const T& n) const {return n*next(n);} 00060 Pow<o-1> next; 00061 }; 00062 00063 template<> struct Pow<1> 00064 { 00065 public: 00066 template<class T> 00067 T operator() (const T& n) const {return n;} 00068 }; 00069 00070 template<> struct Pow<0> 00071 { 00072 public: 00073 template<class T> 00074 T operator() (const T& n) const {return T(1.0);} 00075 }; 00076 00077 00078 00080 template <int s,bool abs=false,class T=TData> class Power 00081 { 00082 public: 00083 Power(){} 00084 T operator() (const T& orig,const T& num) 00085 { 00086 return (*this)(orig,num,(StaticBool<abs>*)(0)); 00087 } 00088 T operator() (const T& orig,const T& num,StaticFalse* doAbsolute) 00089 { 00090 return orig+mP(num); 00091 } 00092 T operator() (const T& orig,const T& num,StaticTrue* doAbsolute) 00093 { 00094 return orig+Abs(mP(num)); 00095 } 00096 protected: 00097 Pow<s> mP; 00098 }; 00099 00100 00102 template<bool abs=false,class T=TData> class NoPowerTmpl 00103 :public Power<1,abs,T>{}; 00105 template<bool abs=false,class T=TData> class SquareTmpl 00106 :public Power<2,abs,T>{}; 00108 template<bool abs=false,class T=TData> class CubeTmpl 00109 :public Power<3,abs,T>{}; 00110 00112 typedef NoPowerTmpl<> NoPower; 00113 typedef SquareTmpl<> Square; 00114 typedef CubeTmpl<> Cube; 00115 00116 00118 template <int s,bool abs=false, class T=TData> class WeightedPower 00119 { 00120 public: 00121 WeightedPower():i(0){} 00122 T operator() (const T& orig,const T& num) 00123 { 00124 return (*this)(orig,num,(StaticBool<abs>*)(0)); 00125 } 00126 T operator() (const T& orig,const T& num,StaticFalse*) 00127 { 00128 return orig+mP(num)*i++; 00129 } 00130 T operator() (const T& orig,const T& num,StaticTrue*) 00131 { 00132 return orig+Abs(mP(num))*i++; 00133 } 00134 protected: 00135 TIndex i; 00136 Pow<s> mP; 00137 }; 00138 00139 00141 template <int s=1,class T=TData> class PoweredProduct 00142 { 00143 public: 00144 PoweredProduct(){} 00145 T operator() (const T& i1,const T& i2) 00146 { 00147 return mP(i1)*i2; 00148 } 00149 protected: 00150 Pow<s> mP; 00151 }; 00152 00153 00155 template<bool abs=false,class T=TData> class WeightedNoPowerTmpl 00156 :public WeightedPower<1,abs,T>{}; 00158 template<bool abs=false,class T=TData> class WeightedSquareTmpl 00159 :public WeightedPower<2,abs,T>{}; 00161 template<bool abs=false,class T=TData> class WeightedCubeTmpl 00162 :public WeightedPower<3,abs,T>{}; 00163 00165 typedef WeightedNoPowerTmpl<> WeightedNoPower; 00166 typedef WeightedSquareTmpl<> WeightedSquare; 00167 typedef WeightedCubeTmpl<> WeightedCube; 00168 00170 template <int s=1,bool abs=false,class T=TData,class U=TData> class BiasedPower 00171 { 00172 public: 00173 BiasedPower(U imean):mean(imean){} 00174 U operator() (const U& orig,const T& num) 00175 { 00176 return (*this)(orig,num,(StaticBool<abs>*)(0)); 00177 } 00178 U operator() (const U& orig,const T& num,StaticFalse*) 00179 { 00180 return orig+mP(num-mean); 00181 } 00182 U operator() (const U& orig,const T& num,StaticTrue*) 00183 { 00184 return orig+mP(Abs(num)-mean); 00185 } 00186 public: 00187 U mean; 00188 Pow<s> mP; 00189 }; 00190 00191 00192 00194 template <class T=TData> class ProductTmpl 00195 { 00196 public: 00197 T operator()(const T& orig,const T& num) 00198 { 00199 return orig*num; 00200 } 00201 }; 00202 00203 typedef ProductTmpl<> Product; 00204 00208 class BaseMemOp 00209 { 00210 public: 00211 BaseMemOp():alreadyComputed(false){} 00212 void Reset(){alreadyComputed=false;} 00213 virtual ~BaseMemOp(){}; 00214 protected: 00215 bool alreadyComputed; 00216 }; 00217 00220 template <int s,bool abs=false,class T=TData> class PoweredSum:public BaseMemOp 00221 { 00222 public: 00223 PoweredSum():memory((T)0.0){} 00224 T operator()(const Array<T>& a,StaticTrue* useMemory=NULL) 00225 { 00226 if (alreadyComputed) return memory; 00227 alreadyComputed=true; 00228 return memory=(*this)(a,(StaticFalse*)(0)); 00229 } 00230 T operator()(const Array<T>& a,StaticFalse*) 00231 { 00232 return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(0.0),mP); 00233 } 00234 private: 00235 T memory; 00236 Power<s,abs,T> mP; 00237 }; 00238 00239 00241 template<bool abs=false,class T=TData> class SumTmpl:public PoweredSum<1,abs,T>{}; 00243 template<bool abs=false,class T=TData> class SquaredSumTmpl:public PoweredSum<2,abs,T>{}; 00245 template<bool abs=false,class T=TData> class CubedSumTmpl:public PoweredSum<3,abs,T>{}; 00246 00247 typedef SumTmpl<> Sum; 00248 typedef SquaredSumTmpl<> SquaredSum; 00249 typedef CubedSumTmpl<> CubedSum; 00250 00251 00253 template <class T=TData> class LogPlusTmpl 00254 { 00255 public: 00256 T operator()(const T& orig,const T& num) 00257 { 00258 return orig+log(Abs(num)); 00259 } 00260 }; 00261 00262 typedef LogPlusTmpl<> LogSum; 00263 00264 // TODO: Remove it as it seems dupplicated 00265 //typedef ProductTmpl<> Product; 00266 00267 00270 template <class T=TData> class LogSumTmpl:public BaseMemOp 00271 { 00272 public: 00273 LogSumTmpl():memory(0.0){} 00274 T operator()(const Array<T>& a,StaticTrue* b=NULL) 00275 { 00276 if (alreadyComputed) return memory; 00277 alreadyComputed=true; 00278 return memory=(*this)(a,(StaticFalse*)(0)); 00279 } 00280 T operator()(const Array<T>& a,StaticFalse*) 00281 { 00282 return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(1.0),LogPlusTmpl<T>()); 00283 00284 } 00285 private: 00286 T memory; 00287 }; 00288 00289 00292 template <class T=TData> class InnerProductTmpl:public BaseMemOp 00293 { 00294 public: 00295 InnerProductTmpl():memory(0.0){} 00296 T operator()(const Array<T>& a,StaticTrue* b=NULL) 00297 { 00298 if (alreadyComputed) return memory; 00299 alreadyComputed=true; 00300 return memory=(*this)(a,(StaticFalse*)(0)); 00301 } 00302 T operator()(const Array<T>& a,StaticFalse*) 00303 { 00304 return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(1.0),std::multiplies<T>()); 00305 } 00306 private: 00307 T memory; 00308 }; 00309 00310 typedef InnerProductTmpl<> InnerProduct; 00311 00314 template <int s, bool abs=false, class T=TData> class WeightedPoweredSum:public BaseMemOp 00315 { 00316 public: 00317 WeightedPoweredSum():memory(0.0){} 00319 T operator()(const Array<T>& a,StaticTrue* b=NULL) 00320 { 00321 if(alreadyComputed) return memory; 00322 alreadyComputed=true; 00323 return memory=(*this)(a,(StaticFalse*)(0)); 00324 } 00326 T operator()(const Array<T>& a,StaticFalse*) 00327 { 00328 return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(0.0),mWP); 00329 } 00330 private: 00331 T memory; 00332 WeightedPower<s,abs,T> mWP; 00333 00334 }; 00335 00336 00339 template <int s,bool abs=false, class T=TData> class CrossWeightedPoweredSum:public BaseMemOp 00340 { 00341 public: 00342 CrossWeightedPoweredSum():memory(0.0){} 00344 T operator()(const Array<T>& x,const Array<T>& y,StaticTrue* b=NULL) 00345 { 00346 if(alreadyComputed) return memory; 00347 alreadyComputed=true; 00348 return memory=(*this)(x,y,(StaticFalse*)(0)); 00349 } 00351 T operator()(const Array<T>& x,const Array<T>& y,StaticFalse*) 00352 { 00353 return inner_product(x.GetPtr(),x.GetPtr()+x.Size(),y.GetPtr(),T(1.0),std::plus<T>(),PoweredProduct<s,T>()); 00354 } 00355 private: 00356 T memory; 00357 00358 }; 00359 00360 00362 template<bool abs=false,class T=TData> class WeightedSumTmpl:public WeightedPoweredSum<1,abs,T>{}; 00364 template<bool abs=false,class T=TData> class WeightedSquaredSumTmpl:public WeightedPoweredSum<2,abs,T>{}; 00366 template<bool abs=false,class T=TData> class WeightedCubedSumTmpl:public WeightedPoweredSum<3,abs,T>{}; 00367 00368 typedef WeightedSumTmpl<> WeightedSum; 00369 typedef WeightedSquaredSumTmpl<> WeightedSquaredSum; 00370 typedef WeightedCubedSumTmpl<> WeightedCubedSum; 00371 00372 00375 template<int o, bool abs=false, class T=TData,class U=TData> class Moment:public BaseMemOp 00376 { 00377 public: 00378 Moment():memory(0.0){} 00380 U operator()(const Array<T>& a,PoweredSum<o,abs,T>& powSum,StaticFalse*) 00381 { 00382 return static_cast<U>(powSum(a))/a.Size(); 00383 } 00385 U operator()(const Array<T>& a,PoweredSum<o,abs,T>& powSum,StaticTrue* b=NULL) 00386 { 00387 if(alreadyComputed) return memory; 00388 alreadyComputed=true; 00389 return memory=(*this)(a,powSum,(StaticFalse*)(0)); 00390 } 00392 U operator()(const Array<T>& a,StaticFalse*) 00393 { 00394 return (*this)(a,mPs,(StaticFalse*)(0)); 00395 } 00397 U operator()(const Array<T>& a,StaticTrue* b=NULL) 00398 { 00399 return (*this)(a,mPs,(StaticTrue*)(0)); 00400 } 00401 00402 protected: 00403 U memory; 00404 PoweredSum<o,abs,T> mPs; 00405 }; 00406 00409 template<int o, bool abs=false, class T=TData,class U=TData> class CenterOfGravity:public BaseMemOp 00410 { 00411 public: 00412 CenterOfGravity():memory(0.0){} 00414 U operator()(const Array<T>& a,WeightedPoweredSum<o,abs,T>& wPowSum,PoweredSum<o,abs,T>& PowSum,StaticFalse*) 00415 { 00416 U normFactor = PowSum( a ); 00417 00418 //MRJ: the zero case. I have set the tolerance to 1e-7, which is appropiate for single 00419 //precision floating point numbers. 00420 if ( normFactor < 1e-7 ) 00421 return (a.Size()%2==0)? a.Size()/2 : (a.Size()+1)/2; 00422 00423 return static_cast<U>(wPowSum(a))/normFactor; 00424 } 00426 U operator()(const Array<T>& a,WeightedPoweredSum<o,abs,T>& wPowSum,PoweredSum<o,abs,T>& PowSum,StaticTrue* b=NULL) 00427 { 00428 if(alreadyComputed) return memory; 00429 alreadyComputed=true; 00430 return memory=(*this)(a,(StaticFalse*)(0)); 00431 } 00433 U operator()(const Array<T>& a,StaticFalse*) 00434 { 00435 return (*this)(a,mWPS,mPS,(StaticFalse*)(0)); 00436 } 00438 U operator()(const Array<T>& a,StaticTrue* b=NULL) 00439 { 00440 return (*this)(a,mWPS,mPS,(StaticTrue*)(0)); 00441 } 00442 00443 protected: 00444 U memory; 00445 WeightedPoweredSum<o,abs,T> mWPS; 00446 PoweredSum<o,abs,T> mPS; 00447 }; 00448 00450 template<int o, bool abs=false, class T=TData,class U=TData> class CrossCenterOfGravity:public BaseMemOp 00451 { 00452 public: 00453 CrossCenterOfGravity():memory(0.0){} 00455 U operator()(const Array<T>& a1,const Array<T>& a2,CrossWeightedPoweredSum<o,abs,T>& cwPowSum,PoweredSum<o,abs,T>& powSum,StaticFalse*) 00456 { 00457 return static_cast<U>(cwPowSum(a1,a2))/powSum(a1); 00458 } 00460 U operator()(const Array<T>& a1,const Array<T>& a2,CrossWeightedPoweredSum<o,abs,T>& cwPowSum,PoweredSum<o,abs,T>& powSum,StaticTrue* b=NULL) 00461 { 00462 if(alreadyComputed) return memory; 00463 alreadyComputed=true; 00464 return memory=(*this)(a1,a2,cwPowSum,powSum,(StaticFalse*)(0)); 00465 } 00467 U operator()(const Array<T>& a1,const Array<T>& a2,StaticFalse*) 00468 { 00469 return (*this)(a1,a2,mWPS,mPS,(StaticFalse*)(0)); 00470 } 00472 U operator()(const Array<T>& a1,const Array<T>& a2,StaticTrue* b=NULL) 00473 { 00474 return (*this)(a1,a2,mWPS,mPS,(StaticTrue*)(0)); 00475 } 00476 00477 protected: 00478 U memory; 00479 CrossWeightedPoweredSum<o,abs,T> mWPS; 00480 PoweredSum<o,abs,T> mPS; 00481 }; 00482 00485 template<bool abs=false,class T=TData,class U=TData> class CentroidTmpl:public CenterOfGravity<1,abs,T,U>{}; 00488 template<bool abs=false,class T=TData,class U=TData> class MeanTmpl:public Moment<1,abs,T,U>{}; 00491 template<class T=TData> class EnergyTmpl:public SquaredSumTmpl<false,T>{}; 00492 00493 typedef CentroidTmpl<> Centroid; 00494 typedef MeanTmpl<> Mean; 00495 typedef EnergyTmpl<> Energy; 00496 00497 00500 template<class T=TData,class U=TData> class RMSTmpl:public BaseMemOp 00501 { 00502 public: 00504 U operator()(const Array<T>& a,SquaredSumTmpl<false,T>& sqrSum,StaticFalse* useMemory) 00505 { 00506 return CLAM_sqrt(sqrSum(a,(StaticFalse*)(0))); 00507 } 00509 U operator()(const Array<T>& a,SquaredSumTmpl<false,T>& sqrSum,StaticTrue* useMemory=NULL) 00510 { 00511 if(!alreadyComputed) 00512 { 00513 memory=(*this)(a, sqrSum, (StaticFalse*)(0)); 00514 alreadyComputed=true; 00515 } 00516 return memory; 00517 } 00519 U operator()(const Array<T>& a,StaticTrue* useMemory=NULL) 00520 { 00521 return (*this)(a,mSS,(StaticTrue*)(0)); 00522 } 00524 U operator()(const Array<T>& a,StaticFalse* useMemory) 00525 { 00526 return (*this)(a,mSS,(StaticFalse*)(0)); 00527 } 00528 private: 00529 U memory; 00530 SquaredSumTmpl<false,T> mSS; 00531 00532 }; 00533 00534 typedef RMSTmpl<> RMS; 00535 00538 template<class T=TData,class U=TData> class GeometricMeanTmpl:public BaseMemOp 00539 { 00540 public: 00541 U operator()(const Array<T>& a,LogSumTmpl<T>& inProd,StaticFalse * useMemory) 00542 { 00543 return exp(inProd(a,(StaticFalse*)(0))*1.0/(double)a.Size()); 00544 } 00545 U operator()(const Array<T>& a,LogSumTmpl<T>& inProd,StaticTrue* useMemory=NULL) 00546 { 00547 if (alreadyComputed) return memory; 00548 alreadyComputed=true; 00549 return memory=(*this)(a, inProd, (StaticFalse*)(0)); 00550 } 00551 U operator()(const Array<T>& a,StaticTrue* useMemory=NULL) 00552 { 00553 return (*this)(a,mIP,(StaticTrue*)(0)); 00554 } 00555 U operator()(const Array<T>& a,StaticFalse* useMemory) 00556 { 00557 return (*this)(a,mIP,(StaticFalse*)(0)); 00558 } 00559 00560 private: 00561 U memory; 00562 LogSumTmpl<T> mIP; 00563 00564 }; 00565 00566 typedef GeometricMeanTmpl<> GeometricMean; 00567 00570 template <int s,bool abs=false,class T=TData,class U=TData> class BiasedPoweredSum:public BaseMemOp 00571 { 00572 public: 00573 BiasedPoweredSum():memory(0.0){} 00574 00575 U operator()(const Array<T>& a, MeanTmpl<abs,T,U>& imean, StaticFalse* useMemory) 00576 { 00577 return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),U(),BiasedPower<s,abs,T,U>(imean(a))); 00578 } 00579 00580 U operator()(const Array<T>& a, MeanTmpl<abs,T,U>& imean, StaticTrue* useMemory=NULL) 00581 { 00582 if (alreadyComputed) return memory; 00583 alreadyComputed=true; 00584 return memory=(*this)(a,(StaticFalse*)(0)); 00585 } 00586 00587 U operator()(const Array<T>& a,StaticTrue* b=NULL) 00588 { 00589 return (*this)(a,mMean,(StaticTrue*)(0)); 00590 } 00591 U operator()(const Array<T>& a,StaticFalse*) 00592 { 00593 return (*this)(a,mMean,(StaticFalse*)(0)); 00594 } 00595 private: 00596 U memory; 00597 MeanTmpl<abs,T,U> mMean; 00598 }; 00599 00603 template<int o,bool abs=false,class T=TData,class U=TData> class CentralMoment:public BaseMemOp 00604 { 00605 public: 00606 CentralMoment():memory(){} 00607 U operator()(const Array<T>& a, BiasedPoweredSum<o,abs,T,U>& bps, StaticTrue* b=NULL) 00608 { 00609 if (alreadyComputed) return memory; 00610 alreadyComputed=true; 00611 return memory=(*this)(a,bps,(StaticFalse*)(0)); 00612 } 00613 U operator()(const Array<T>& a, BiasedPoweredSum<o,abs,T,U>& bps, StaticFalse*) 00614 { 00615 return static_cast<U>(bps(a))/a.Size(); 00616 } 00617 U operator()(const Array<T>& a, StaticTrue* b=NULL) 00618 { 00619 return (*this)(a,mBPS, (StaticTrue*)(0)); 00620 } 00621 U operator()(const Array<T>& a, StaticFalse*) 00622 { 00623 return (*this)(a,mBPS, (StaticFalse*)(0)); 00624 } 00625 00627 U operator()(const Array<T>& a, Array<BaseMemOp*>& moments, StaticTrue* b=NULL) 00628 { 00629 if (alreadyComputed) return memory; 00630 alreadyComputed=true; 00631 return memory=(*this)(a,moments,(StaticFalse*)(0)); 00632 } 00633 U operator()(const Array<T>& a, Array<BaseMemOp*>& moments, StaticFalse*) 00634 { 00635 CLAM_DEBUG_ASSERT(moments.Size()>=o,"Central Moment: you need as many raw moments as the order of the central moment you want to compute"); 00636 return (*this)(a,moments,(O<o>*)(0)); 00637 } 00638 00639 U operator()(const Array<T>& a, Array<BaseMemOp*>& moments, O<1>*) 00640 { 00641 return 0; 00642 } 00643 00644 U operator()(const Array<T>& a, Array<BaseMemOp*>& moments, O<2>*) 00645 { 00646 // -m1² + m2 00647 U m1 = (*(dynamic_cast<Moment<1,abs,T,U>*>(moments[0])))(a); 00648 U m2 = (*(dynamic_cast<Moment<2,abs,T,U>*>(moments[1])))(a); 00649 return (-1)*m1*m1 + m2; 00650 } 00651 00652 U operator()(const Array<T>& a, Array<BaseMemOp*>& moments, O<3>*) 00653 { 00654 // 2*m1³ - 3*m1*m2 + m3 = ... 5 Mult 00655 // m1*(2*m1² - 3*m2) + m3 ... 4 Mult 00656 U m1 = (*(dynamic_cast<Moment<1,abs,T,U>*>(moments[0])))(a); 00657 U m2 = (*(dynamic_cast<Moment<2,abs,T,U>*>(moments[1])))(a); 00658 U m3 = (*(dynamic_cast<Moment<3,abs,T,U>*>(moments[2])))(a); 00659 return m1*(2*m1*m1 - 3*m2) + m3; 00660 } 00661 00662 U operator()(const Array<T>& a, Array<BaseMemOp*>& moments, O<4>*) 00663 { 00664 // -3*m1^4 + 6*m1²*m2 - 4*m1*m3 + m4 ... 9 Mult 00665 // m1*(m1*((-3)*m1² + 6*m2) - 4*m3) + m4 ... 6 Mult 00666 U m1 = (*(dynamic_cast<Moment<1,abs,T,U>*>(moments[0])))(a); 00667 U m2 = (*(dynamic_cast<Moment<2,abs,T,U>*>(moments[1])))(a); 00668 U m3 = (*(dynamic_cast<Moment<3,abs,T,U>*>(moments[2])))(a); 00669 U m4 = (*(dynamic_cast<Moment<4,abs,T,U>*>(moments[3])))(a); 00670 return m1*(m1*((-3)*m1*m1 + 6*m2) - 4*m3) + m4; 00671 } 00672 00673 U operator()(const Array<T>& a, Array<BaseMemOp*>& moments, O<5>*) 00674 { 00675 // 4*u1^5 - 10*u1³*u2 + 10*u1²*u3 - 5*u1*u4+u5 = .... 14 Mult 00676 // u1*(u1*(u1*(4*u1² - 10*u2) + 10*u3) - 5*u4) + u5 .... 8 Mult 00677 U m1 = (*(dynamic_cast<Moment<1,abs,T,U>*>(moments[0])))(a); 00678 U m2 = (*(dynamic_cast<Moment<2,abs,T,U>*>(moments[1])))(a); 00679 U m3 = (*(dynamic_cast<Moment<3,abs,T,U>*>(moments[2])))(a); 00680 U m4 = (*(dynamic_cast<Moment<4,abs,T,U>*>(moments[3])))(a); 00681 U m5 = (*(dynamic_cast<Moment<5,abs,T,U>*>(moments[4])))(a); 00682 00683 return m1*(m1*(m1*(4*m1*m1 - 10*m2) + 10*m3) - 5*m4) + m5; 00684 } 00685 00686 protected: 00687 BiasedPoweredSum<o,abs,T,U> mBPS; 00688 U memory; 00689 00690 }; 00691 00692 00695 template <bool abs=false, class T=TData, class U=TData> class StandardDeviationTmpl:public BaseMemOp 00696 { 00697 public: 00698 StandardDeviationTmpl():memory(){} 00699 U operator()(const Array<T>& a,CentralMoment<2,abs,T,U>& centralMoment2,bool useMemory=false) 00700 { 00701 if(!useMemory) return CLAM_sqrt(centralMoment2(a)); 00702 00703 if (alreadyComputed) return memory; 00704 alreadyComputed=true; 00705 return memory=(*this)(a,centralMoment2,false); 00706 } 00708 U operator()(const Array<T>& a,bool useMemory=false) 00709 { 00710 return (*this)(a,mCM2,useMemory); 00711 00712 } 00713 00714 00715 protected: 00716 U memory; 00717 CentralMoment<2,abs,T,U> mCM2; 00718 00719 }; 00720 00721 typedef StandardDeviationTmpl<> StandardDeviation; 00722 00723 00726 template <bool abs=false,class T=TData,class U=TData> class SkewTmpl:public BaseMemOp 00727 { 00728 00729 public: 00730 SkewTmpl():memory(){} 00732 U operator()(const Array<T>& data, StandardDeviationTmpl<abs,T,U>& std, CentralMoment<3,abs,T,U>& ctrMnt3, bool useMemory=false) 00733 { 00734 if(!useMemory) return MemorylessCompute(data, std, ctrMnt3); 00735 00736 if (alreadyComputed) return memory; 00737 alreadyComputed=true; 00738 return memory=(*this)(data, std, ctrMnt3, false); 00739 } 00741 U operator()(const Array<T>& data, bool useMemory=false) 00742 { 00744 return (*this)(data, mSD, mCM3, useMemory); 00745 } 00746 00747 protected: 00748 U MemorylessCompute(const Array<T>& data, StandardDeviationTmpl<abs,T,U>& std, CentralMoment<3,abs,T,U>& ctrMnt3) 00749 { 00750 // When the values tend to be the same, skew tends to be 0 (simetric) 00751 U tmpStd = CLAM_max(U(1e-10),std(data)); 00752 U tmpCentralMoment3 = ctrMnt3(data); 00753 return tmpCentralMoment3/(tmpStd*tmpStd*tmpStd); 00754 } 00755 U memory; 00756 StandardDeviationTmpl<abs,T,U> mSD; 00757 CentralMoment<3,abs,T,U> mCM3; 00758 }; 00759 00760 typedef SkewTmpl<> Skew; 00761 00764 template <bool abs=false,class T=TData,class U=TData> class KurtosisTmpl:public BaseMemOp 00765 { 00766 public: 00767 KurtosisTmpl():memory(){} 00769 U operator()(const Array<T>& a,CentralMoment<2,abs,T,U>& var,CentralMoment<4,abs,T,U>& ctrMnt4,bool useMemory=false) 00770 { 00771 if(!useMemory) return MemoryLessCompute(a, var, ctrMnt4); 00772 00773 if (alreadyComputed) return memory; 00774 alreadyComputed=true; 00775 return memory=(*this)(a,var,ctrMnt4,false); 00776 } 00778 U operator()(const Array<T>& a,bool useMemory=false) 00779 { 00780 return (*this)(a,mCM2,mCM4,useMemory); 00781 } 00782 00783 protected: 00784 U MemoryLessCompute(const Array<T>& a,CentralMoment<2,abs,T,U>& var,CentralMoment<4,abs,T,U>& ctrMnt4) 00785 { 00786 U variance = var(a); 00787 if (variance<U(1e-10)) return U(3.0); 00788 U centerMoment4 = ctrMnt4(a); 00789 return centerMoment4/(variance*variance); 00790 } 00791 U memory; 00792 CentralMoment<2,abs,T,U> mCM2; 00793 CentralMoment<4,abs,T,U> mCM4; 00794 }; 00795 00796 typedef KurtosisTmpl<> Kurtosis; 00797 00798 00800 template <bool abs=false, class T=TData> class ComplexMin 00801 { 00802 public: 00803 ComplexMin(){} 00804 T operator() (const T& orig,const T& num) 00805 { 00806 return (*this)(orig,num,(StaticBool<abs>*)(0)); 00807 } 00808 T operator() (const T& orig,const T& num,StaticFalse*) 00809 { 00810 return CLAM_min(orig,num); 00811 } 00812 T operator() (const T& orig,const T& num,StaticTrue*) 00813 { 00814 return CLAM_min(orig,Abs(num)); 00815 } 00816 }; 00817 00820 template <bool abs=false,class T=TData> class ComplexMinElement:public BaseMemOp 00821 { 00822 public: 00823 ComplexMinElement():memory((T)0.0){} 00824 T operator()(const Array<T>& a,StaticTrue* b=NULL) 00825 { 00826 if (alreadyComputed) return memory; 00827 alreadyComputed=true; 00828 return memory=(*this)(a,(StaticFalse*)(0)); 00829 } 00830 T operator()(const Array<T>& a,StaticFalse*) 00831 { 00832 return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(a[0]),mM); 00833 } 00834 private: 00835 T memory; 00836 ComplexMin<abs,T> mM; 00837 }; 00838 00840 template <bool abs=false, class T=TData> class ComplexMax 00841 { 00842 public: 00843 ComplexMax(){} 00844 T operator() (const T& orig,const T& num) 00845 { 00846 return (*this)(orig,num,(StaticBool<abs>*)(0)); 00847 } 00848 T operator() (const T& orig,const T& num,StaticFalse*) 00849 { 00850 return CLAM_max(orig,num); 00851 } 00852 T operator() (const T& orig,const T& num,StaticTrue*) 00853 { 00854 return CLAM_max(orig,Abs(num)); 00855 } 00856 }; 00857 00860 template <bool abs=false,class T=TData> class ComplexMaxElement:public BaseMemOp 00861 { 00862 public: 00863 ComplexMaxElement():memory((T)0.0){} 00864 T operator()(const Array<T>& a,StaticTrue* b=NULL) 00865 { 00866 if (alreadyComputed) return memory; 00867 alreadyComputed=true; 00868 return memory=(*this)(a,(StaticFalse*)(0)); 00869 } 00870 T operator()(const Array<T>& a,StaticFalse*) 00871 { 00872 return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(a[0]),mM); 00873 } 00874 private: 00875 T memory; 00876 ComplexMax<abs,T> mM; 00877 }; 00878 00879 00880 00881 }; 00882 00883 #endif// _BasicOps_ 00884