BasicOps.hxx

Go to the documentation of this file.
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 
Generated by  doxygen 1.6.3