CLAM-Development  1.4.0
BasicOps.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2004 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 #ifndef _BasicOps_
23 #define _BasicOps_
24 
25 #include <numeric>
26 #include <functional>
27 #include "StaticBool.hxx"
28 #include "DataTypes.hxx"
29 #include "Array.hxx"
30 #include "CLAM_Math.hxx"
31 #include "Order.hxx"
32 
33 
34 using std::accumulate;
35 using std::inner_product;
36 using std::mem_fun;
37 
38 
51 namespace CLAM {
52 
53 
54 
55 template <int o> struct Pow
56 {
57 public:
58  template <class T>
59  T operator () (const T& n) const {return n*next(n);}
60  Pow<o-1> next;
61 };
62 
63 template<> struct Pow<1>
64 {
65 public:
66  template<class T>
67  T operator() (const T& n) const {return n;}
68 };
69 
70 template<> struct Pow<0>
71 {
72 public:
73  template<class T>
74  T operator() (const T& n) const {return T(1.0);}
75 };
76 
77 
78 
80 template <int s,bool abs=false,class T=TData> class Power
81 {
82 public:
83  Power(){}
84  T operator() (const T& orig,const T& num)
85  {
86  return (*this)(orig,num,(StaticBool<abs>*)(0));
87  }
88  T operator() (const T& orig,const T& num,StaticFalse* doAbsolute)
89  {
90  return orig+mP(num);
91  }
92  T operator() (const T& orig,const T& num,StaticTrue* doAbsolute)
93  {
94  return orig+Abs(mP(num));
95  }
96 protected:
98 };
99 
100 
102 template<bool abs=false,class T=TData> class NoPowerTmpl
103 :public Power<1,abs,T>{};
105 template<bool abs=false,class T=TData> class SquareTmpl
106 :public Power<2,abs,T>{};
108 template<bool abs=false,class T=TData> class CubeTmpl
109 :public Power<3,abs,T>{};
110 
114 typedef CubeTmpl<> Cube;
115 
116 
118 template <int s,bool abs=false, class T=TData> class WeightedPower
119 {
120 public:
122  T operator() (const T& orig,const T& num)
123  {
124  return (*this)(orig,num,(StaticBool<abs>*)(0));
125  }
126  T operator() (const T& orig,const T& num,StaticFalse*)
127  {
128  return orig+mP(num)*i++;
129  }
130  T operator() (const T& orig,const T& num,StaticTrue*)
131  {
132  return orig+Abs(mP(num))*i++;
133  }
134 protected:
137 };
138 
139 
141 template <int s=1,class T=TData> class PoweredProduct
142 {
143 public:
145  T operator() (const T& i1,const T& i2)
146  {
147  return mP(i1)*i2;
148  }
149 protected:
151 };
152 
153 
155 template<bool abs=false,class T=TData> class WeightedNoPowerTmpl
156 :public WeightedPower<1,abs,T>{};
158 template<bool abs=false,class T=TData> class WeightedSquareTmpl
159 :public WeightedPower<2,abs,T>{};
161 template<bool abs=false,class T=TData> class WeightedCubeTmpl
162 :public WeightedPower<3,abs,T>{};
163 
168 
170 template <int s=1,bool abs=false,class T=TData,class U=TData> class BiasedPower
171 {
172 public:
173  BiasedPower(U imean):mean(imean){}
174  U operator() (const U& orig,const T& num)
175  {
176  return (*this)(orig,num,(StaticBool<abs>*)(0));
177  }
178  U operator() (const U& orig,const T& num,StaticFalse*)
179  {
180  return orig+mP(num-mean);
181  }
182  U operator() (const U& orig,const T& num,StaticTrue*)
183  {
184  return orig+mP(Abs(num)-mean);
185  }
186 public:
189 };
190 
191 
192 
194 template <class T=TData> class ProductTmpl
195 {
196 public:
197  T operator()(const T& orig,const T& num)
198  {
199  return orig*num;
200  }
201 };
202 
204 
209 {
210 public:
212  void Reset(){alreadyComputed=false;}
213  virtual ~BaseMemOp(){};
214 protected:
215  bool alreadyComputed;
216 };
217 
220 template <int s,bool abs=false,class T=TData> class PoweredSum:public BaseMemOp
221 {
222 public:
223  PoweredSum():memory((T)0.0){}
224  T operator()(const Array<T>& a,StaticTrue* useMemory=NULL)
225  {
226  if (alreadyComputed) return memory;
227  alreadyComputed=true;
228  return memory=(*this)(a,(StaticFalse*)(0));
229  }
231  {
232  return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(0.0),mP);
233  }
234 private:
235  T memory;
236  Power<s,abs,T> mP;
237 };
238 
239 
241 template<bool abs=false,class T=TData> class SumTmpl:public PoweredSum<1,abs,T>{};
243 template<bool abs=false,class T=TData> class SquaredSumTmpl:public PoweredSum<2,abs,T>{};
245 template<bool abs=false,class T=TData> class CubedSumTmpl:public PoweredSum<3,abs,T>{};
246 
247 typedef SumTmpl<> Sum;
250 
251 
253 template <class T=TData> class LogPlusTmpl
254 {
255 public:
256  T operator()(const T& orig,const T& num)
257  {
258  return orig+log(Abs(num));
259  }
260 };
261 
263 
264 // TODO: Remove it as it seems dupplicated
265 //typedef ProductTmpl<> Product;
266 
267 
270 template <class T=TData> class LogSumTmpl:public BaseMemOp
271 {
272 public:
273  LogSumTmpl():memory(0.0){}
275  {
276  if (alreadyComputed) return memory;
277  alreadyComputed=true;
278  return memory=(*this)(a,(StaticFalse*)(0));
279  }
281  {
282  return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(1.0),LogPlusTmpl<T>());
283 
284  }
285 private:
286  T memory;
287 };
288 
289 
292 template <class T=TData> class InnerProductTmpl:public BaseMemOp
293 {
294 public:
295  InnerProductTmpl():memory(0.0){}
297  {
298  if (alreadyComputed) return memory;
299  alreadyComputed=true;
300  return memory=(*this)(a,(StaticFalse*)(0));
301  }
303  {
304  return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(1.0),std::multiplies<T>());
305  }
306 private:
307  T memory;
308 };
309 
311 
314 template <int s, bool abs=false, class T=TData> class WeightedPoweredSum:public BaseMemOp
315 {
316 public:
317  WeightedPoweredSum():memory(0.0){}
320  {
321  if(alreadyComputed) return memory;
322  alreadyComputed=true;
323  return memory=(*this)(a,(StaticFalse*)(0));
324  }
327  {
328  return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(0.0),mWP);
329  }
330 private:
331  T memory;
333 
334 };
335 
336 
339 template <int s,bool abs=false, class T=TData> class CrossWeightedPoweredSum:public BaseMemOp
340 {
341 public:
342  CrossWeightedPoweredSum():memory(0.0){}
344  T operator()(const Array<T>& x,const Array<T>& y,StaticTrue* b=NULL)
345  {
346  if(alreadyComputed) return memory;
347  alreadyComputed=true;
348  return memory=(*this)(x,y,(StaticFalse*)(0));
349  }
351  T operator()(const Array<T>& x,const Array<T>& y,StaticFalse*)
352  {
353  return inner_product(x.GetPtr(),x.GetPtr()+x.Size(),y.GetPtr(),T(1.0),std::plus<T>(),PoweredProduct<s,T>());
354  }
355 private:
356  T memory;
357 
358 };
359 
360 
362 template<bool abs=false,class T=TData> class WeightedSumTmpl:public WeightedPoweredSum<1,abs,T>{};
364 template<bool abs=false,class T=TData> class WeightedSquaredSumTmpl:public WeightedPoweredSum<2,abs,T>{};
366 template<bool abs=false,class T=TData> class WeightedCubedSumTmpl:public WeightedPoweredSum<3,abs,T>{};
367 
371 
372 
375 template<int o, bool abs=false, class T=TData,class U=TData> class Moment:public BaseMemOp
376 {
377 public:
378  Moment():memory(0.0){}
381  {
382  return static_cast<U>(powSum(a))/a.Size();
383  }
386  {
387  if(alreadyComputed) return memory;
388  alreadyComputed=true;
389  return memory=(*this)(a,powSum,(StaticFalse*)(0));
390  }
393  {
394  return (*this)(a,mPs,(StaticFalse*)(0));
395  }
398  {
399  return (*this)(a,mPs,(StaticTrue*)(0));
400  }
401 
402 protected:
405 };
406 
409 template<int o, bool abs=false, class T=TData,class U=TData> class CenterOfGravity:public BaseMemOp
410 {
411 public:
415  {
416  U normFactor = PowSum( a );
417 
418  //MRJ: the zero case. I have set the tolerance to 1e-7, which is appropiate for single
419  //precision floating point numbers.
420  if ( normFactor < 1e-7 )
421  return (a.Size()%2==0)? a.Size()/2 : (a.Size()+1)/2;
422 
423  return static_cast<U>(wPowSum(a))/normFactor;
424  }
427  {
428  if(alreadyComputed) return memory;
429  alreadyComputed=true;
430  return memory=(*this)(a,(StaticFalse*)(0));
431  }
434  {
435  return (*this)(a,mWPS,mPS,(StaticFalse*)(0));
436  }
439  {
440  return (*this)(a,mWPS,mPS,(StaticTrue*)(0));
441  }
442 
443 protected:
447 };
448 
450 template<int o, bool abs=false, class T=TData,class U=TData> class CrossCenterOfGravity:public BaseMemOp
451 {
452 public:
456  {
457  return static_cast<U>(cwPowSum(a1,a2))/powSum(a1);
458  }
461  {
462  if(alreadyComputed) return memory;
463  alreadyComputed=true;
464  return memory=(*this)(a1,a2,cwPowSum,powSum,(StaticFalse*)(0));
465  }
467  U operator()(const Array<T>& a1,const Array<T>& a2,StaticFalse*)
468  {
469  return (*this)(a1,a2,mWPS,mPS,(StaticFalse*)(0));
470  }
472  U operator()(const Array<T>& a1,const Array<T>& a2,StaticTrue* b=NULL)
473  {
474  return (*this)(a1,a2,mWPS,mPS,(StaticTrue*)(0));
475  }
476 
477 protected:
481 };
482 
485 template<bool abs=false,class T=TData,class U=TData> class CentroidTmpl:public CenterOfGravity<1,abs,T,U>{};
488 template<bool abs=false,class T=TData,class U=TData> class MeanTmpl:public Moment<1,abs,T,U>{};
491 template<class T=TData> class EnergyTmpl:public SquaredSumTmpl<false,T>{};
492 
494 typedef MeanTmpl<> Mean;
496 
497 
500 template<class T=TData,class U=TData> class RMSTmpl:public BaseMemOp
501 {
502 public:
505  {
506  return CLAM_sqrt(sqrSum(a,(StaticFalse*)(0)));
507  }
510  {
511  if(!alreadyComputed)
512  {
513  memory=(*this)(a, sqrSum, (StaticFalse*)(0));
514  alreadyComputed=true;
515  }
516  return memory;
517  }
519  U operator()(const Array<T>& a,StaticTrue* useMemory=NULL)
520  {
521  return (*this)(a,mSS,(StaticTrue*)(0));
522  }
524  U operator()(const Array<T>& a,StaticFalse* useMemory)
525  {
526  return (*this)(a,mSS,(StaticFalse*)(0));
527  }
528 private:
529  U memory;
531 
532 };
533 
534 typedef RMSTmpl<> RMS;
535 
538 template<class T=TData,class U=TData> class GeometricMeanTmpl:public BaseMemOp
539 {
540 public:
541  U operator()(const Array<T>& a,LogSumTmpl<T>& inProd,StaticFalse * useMemory)
542  {
543  return exp(inProd(a,(StaticFalse*)(0))*1.0/(double)a.Size());
544  }
545  U operator()(const Array<T>& a,LogSumTmpl<T>& inProd,StaticTrue* useMemory=NULL)
546  {
547  if (alreadyComputed) return memory;
548  alreadyComputed=true;
549  return memory=(*this)(a, inProd, (StaticFalse*)(0));
550  }
551  U operator()(const Array<T>& a,StaticTrue* useMemory=NULL)
552  {
553  return (*this)(a,mIP,(StaticTrue*)(0));
554  }
555  U operator()(const Array<T>& a,StaticFalse* useMemory)
556  {
557  return (*this)(a,mIP,(StaticFalse*)(0));
558  }
559 
560 private:
561  U memory;
562  LogSumTmpl<T> mIP;
563 
564 };
565 
567 
570 template <int s,bool abs=false,class T=TData,class U=TData> class BiasedPoweredSum:public BaseMemOp
571 {
572 public:
573  BiasedPoweredSum():memory(0.0){}
574 
575  U operator()(const Array<T>& a, MeanTmpl<abs,T,U>& imean, StaticFalse* useMemory)
576  {
577  return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),U(),BiasedPower<s,abs,T,U>(imean(a)));
578  }
579 
580  U operator()(const Array<T>& a, MeanTmpl<abs,T,U>& imean, StaticTrue* useMemory=NULL)
581  {
582  if (alreadyComputed) return memory;
583  alreadyComputed=true;
584  return memory=(*this)(a,(StaticFalse*)(0));
585  }
586 
588  {
589  return (*this)(a,mMean,(StaticTrue*)(0));
590  }
592  {
593  return (*this)(a,mMean,(StaticFalse*)(0));
594  }
595 private:
596  U memory;
597  MeanTmpl<abs,T,U> mMean;
598 };
599 
603 template<int o,bool abs=false,class T=TData,class U=TData> class CentralMoment:public BaseMemOp
604 {
605 public:
608  {
609  if (alreadyComputed) return memory;
610  alreadyComputed=true;
611  return memory=(*this)(a,bps,(StaticFalse*)(0));
612  }
614  {
615  return static_cast<U>(bps(a))/a.Size();
616  }
618  {
619  return (*this)(a,mBPS, (StaticTrue*)(0));
620  }
622  {
623  return (*this)(a,mBPS, (StaticFalse*)(0));
624  }
625 
628  {
629  if (alreadyComputed) return memory;
630  alreadyComputed=true;
631  return memory=(*this)(a,moments,(StaticFalse*)(0));
632  }
634  {
635  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");
636  return (*this)(a,moments,(O<o>*)(0));
637  }
638 
640  {
641  return 0;
642  }
643 
645  {
646  // -m1² + m2
647  U m1 = (*(dynamic_cast<Moment<1,abs,T,U>*>(moments[0])))(a);
648  U m2 = (*(dynamic_cast<Moment<2,abs,T,U>*>(moments[1])))(a);
649  return (-1)*m1*m1 + m2;
650  }
651 
653  {
654  // 2*m1³ - 3*m1*m2 + m3 = ... 5 Mult
655  // m1*(2*m1² - 3*m2) + m3 ... 4 Mult
656  U m1 = (*(dynamic_cast<Moment<1,abs,T,U>*>(moments[0])))(a);
657  U m2 = (*(dynamic_cast<Moment<2,abs,T,U>*>(moments[1])))(a);
658  U m3 = (*(dynamic_cast<Moment<3,abs,T,U>*>(moments[2])))(a);
659  return m1*(2*m1*m1 - 3*m2) + m3;
660  }
661 
663  {
664  // -3*m1^4 + 6*m1²*m2 - 4*m1*m3 + m4 ... 9 Mult
665  // m1*(m1*((-3)*m1² + 6*m2) - 4*m3) + m4 ... 6 Mult
666  U m1 = (*(dynamic_cast<Moment<1,abs,T,U>*>(moments[0])))(a);
667  U m2 = (*(dynamic_cast<Moment<2,abs,T,U>*>(moments[1])))(a);
668  U m3 = (*(dynamic_cast<Moment<3,abs,T,U>*>(moments[2])))(a);
669  U m4 = (*(dynamic_cast<Moment<4,abs,T,U>*>(moments[3])))(a);
670  return m1*(m1*((-3)*m1*m1 + 6*m2) - 4*m3) + m4;
671  }
672 
674  {
675  // 4*u1^5 - 10*u1³*u2 + 10*u1²*u3 - 5*u1*u4+u5 = .... 14 Mult
676  // u1*(u1*(u1*(4*u1² - 10*u2) + 10*u3) - 5*u4) + u5 .... 8 Mult
677  U m1 = (*(dynamic_cast<Moment<1,abs,T,U>*>(moments[0])))(a);
678  U m2 = (*(dynamic_cast<Moment<2,abs,T,U>*>(moments[1])))(a);
679  U m3 = (*(dynamic_cast<Moment<3,abs,T,U>*>(moments[2])))(a);
680  U m4 = (*(dynamic_cast<Moment<4,abs,T,U>*>(moments[3])))(a);
681  U m5 = (*(dynamic_cast<Moment<5,abs,T,U>*>(moments[4])))(a);
682 
683  return m1*(m1*(m1*(4*m1*m1 - 10*m2) + 10*m3) - 5*m4) + m5;
684  }
685 
686 protected:
689 
690 };
691 
692 
695 template <bool abs=false, class T=TData, class U=TData> class StandardDeviationTmpl:public BaseMemOp
696 {
697 public:
699  U operator()(const Array<T>& a,CentralMoment<2,abs,T,U>& centralMoment2,bool useMemory=false)
700  {
701  if(!useMemory) return CLAM_sqrt(centralMoment2(a));
702 
703  if (alreadyComputed) return memory;
704  alreadyComputed=true;
705  return memory=(*this)(a,centralMoment2,false);
706  }
708  U operator()(const Array<T>& a,bool useMemory=false)
709  {
710  return (*this)(a,mCM2,useMemory);
711 
712  }
713 
714 
715 protected:
718 
719 };
720 
722 
723 
726 template <bool abs=false,class T=TData,class U=TData> class SkewTmpl:public BaseMemOp
727 {
728 
729 public:
732  U operator()(const Array<T>& data, StandardDeviationTmpl<abs,T,U>& std, CentralMoment<3,abs,T,U>& ctrMnt3, bool useMemory=false)
733  {
734  if(!useMemory) return MemorylessCompute(data, std, ctrMnt3);
735 
736  if (alreadyComputed) return memory;
737  alreadyComputed=true;
738  return memory=(*this)(data, std, ctrMnt3, false);
739  }
741  U operator()(const Array<T>& data, bool useMemory=false)
742  {
744  return (*this)(data, mSD, mCM3, useMemory);
745  }
746 
747 protected:
749  {
750  // When the values tend to be the same, skew tends to be 0 (simetric)
751  U tmpStd = CLAM_max(U(1e-10),std(data));
752  U tmpCentralMoment3 = ctrMnt3(data);
753  return tmpCentralMoment3/(tmpStd*tmpStd*tmpStd);
754  }
758 };
759 
760 typedef SkewTmpl<> Skew;
761 
764 template <bool abs=false,class T=TData,class U=TData> class KurtosisTmpl:public BaseMemOp
765 {
766 public:
769  U operator()(const Array<T>& a,CentralMoment<2,abs,T,U>& var,CentralMoment<4,abs,T,U>& ctrMnt4,bool useMemory=false)
770  {
771  if(!useMemory) return MemoryLessCompute(a, var, ctrMnt4);
772 
773  if (alreadyComputed) return memory;
774  alreadyComputed=true;
775  return memory=(*this)(a,var,ctrMnt4,false);
776  }
778  U operator()(const Array<T>& a,bool useMemory=false)
779  {
780  return (*this)(a,mCM2,mCM4,useMemory);
781  }
782 
783 protected:
785  {
786  U variance = var(a);
787  if (variance<U(1e-10)) return U(3.0);
788  U centerMoment4 = ctrMnt4(a);
789  return centerMoment4/(variance*variance);
790  }
794 };
795 
797 
798 
800 template <bool abs=false, class T=TData> class ComplexMin
801 {
802 public:
804  T operator() (const T& orig,const T& num)
805  {
806  return (*this)(orig,num,(StaticBool<abs>*)(0));
807  }
808  T operator() (const T& orig,const T& num,StaticFalse*)
809  {
810  return CLAM_min(orig,num);
811  }
812  T operator() (const T& orig,const T& num,StaticTrue*)
813  {
814  return CLAM_min(orig,Abs(num));
815  }
816 };
817 
820 template <bool abs=false,class T=TData> class ComplexMinElement:public BaseMemOp
821 {
822 public:
823  ComplexMinElement():memory((T)0.0){}
825  {
826  if (alreadyComputed) return memory;
827  alreadyComputed=true;
828  return memory=(*this)(a,(StaticFalse*)(0));
829  }
831  {
832  return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(a[0]),mM);
833  }
834 private:
835  T memory;
837 };
838 
840 template <bool abs=false, class T=TData> class ComplexMax
841 {
842 public:
844  T operator() (const T& orig,const T& num)
845  {
846  return (*this)(orig,num,(StaticBool<abs>*)(0));
847  }
848  T operator() (const T& orig,const T& num,StaticFalse*)
849  {
850  return CLAM_max(orig,num);
851  }
852  T operator() (const T& orig,const T& num,StaticTrue*)
853  {
854  return CLAM_max(orig,Abs(num));
855  }
856 };
857 
860 template <bool abs=false,class T=TData> class ComplexMaxElement:public BaseMemOp
861 {
862 public:
863  ComplexMaxElement():memory((T)0.0){}
865  {
866  if (alreadyComputed) return memory;
867  alreadyComputed=true;
868  return memory=(*this)(a,(StaticFalse*)(0));
869  }
871  {
872  return accumulate(a.GetPtr(),a.GetPtr()+a.Size(),T(a[0]),mM);
873  }
874 private:
875  T memory;
877 };
878 
879 
880 
881 };
882 
883 #endif// _BasicOps_
884