20 #ifndef _STATFACTORY_H_
21 #define _STATFACTORY_H_
32 #include <gsl/gsl_fit.h>
33 #include <gsl/gsl_errno.h>
34 #include <gsl/gsl_spline.h>
35 #include <gsl/gsl_rng.h>
36 #include <gsl/gsl_randist.h>
37 #include <gsl/gsl_cdf.h>
45 enum INTERPOLATION_TYPE {undefined=0,linear=1,polynomial=2,cspline=3,cspline_periodic=4,akima=5,akima_periodic=6};
47 enum DISTRIBUTION_TYPE {uniform=1,gaussian=2};
51 INTERPOLATION_TYPE getInterpolationType(
const std::string interpolationType){
52 std::map<std::string, INTERPOLATION_TYPE> m_interpMap;
54 return m_interpMap[interpolationType];
56 DISTRIBUTION_TYPE getDistributionType(
const std::string distributionType){
57 std::map<std::string, DISTRIBUTION_TYPE> m_distMap;
59 return m_distMap[distributionType];
61 static void allocAcc(gsl_interp_accel *&acc){
62 acc = gsl_interp_accel_alloc ();
65 static void getSpline(
const std::string type,
int size, gsl_spline *& spline){
66 std::map<std::string, INTERPOLATION_TYPE> m_interpMap;
68 switch(m_interpMap[type]){
70 spline=gsl_spline_alloc(gsl_interp_polynomial,size);
73 spline=gsl_spline_alloc(gsl_interp_cspline,size);
75 case(cspline_periodic):
76 spline=gsl_spline_alloc(gsl_interp_cspline_periodic,size);
79 spline=gsl_spline_alloc(gsl_interp_akima,size);
82 spline=gsl_spline_alloc(gsl_interp_akima_periodic,size);
86 spline=gsl_spline_alloc(gsl_interp_linear,size);
91 static void initSpline(gsl_spline *spline,
const double *x,
const double *y,
int size){
92 gsl_spline_init (spline, x, y, size);
94 static double evalSpline(gsl_spline *spline,
double x, gsl_interp_accel *acc){
95 return gsl_spline_eval (spline, x, acc);
98 static gsl_rng* getRandomGenerator(
unsigned long int theSeed){
99 const gsl_rng_type * T;
103 r = gsl_rng_alloc (gsl_rng_mt19937);
104 gsl_rng_set(r,theSeed);
107 void getNodataValues(std::vector<double>& nodatav)
const{nodatav=m_noDataValues;};
108 bool isNoData(
double value)
const{
109 if(m_noDataValues.empty())
112 return find(m_noDataValues.begin(),m_noDataValues.end(),value)!=m_noDataValues.end();
114 unsigned int pushNodDataValue(
double noDataValue){
115 if(find(m_noDataValues.begin(),m_noDataValues.end(),noDataValue)==m_noDataValues.end())
116 m_noDataValues.push_back(noDataValue);
117 return m_noDataValues.size();
119 unsigned int setNoDataValues(std::vector<double> vnodata){
120 m_noDataValues=vnodata;
121 return m_noDataValues.size();
123 double getRandomValue(
const gsl_rng* r,
const std::string type,
double a=0,
double b=1)
const{
124 std::map<std::string, DISTRIBUTION_TYPE> m_distMap;
127 switch(m_distMap[type]){
129 randValue = a+gsl_rng_uniform(r);
133 randValue = a+gsl_ran_gaussian(r, b);
140 template<
class T> T mymin(
const typename std::vector<T>& v)
const;
141 template<
class T> T mymax(
const typename std::vector<T>& v)
const;
142 template<
class T> T mymin(
const typename std::vector<T>& v, T minConstraint)
const;
143 template<
class T> T mymax(
const typename std::vector<T>& v, T maxConstraint)
const;
145 template<
class T>
typename std::vector<T>::const_iterator mymin(
const typename std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end)
const;
146 template<
class T>
typename std::vector<T>::iterator mymin(
const typename std::vector<T>& v,
typename std::vector<T>::iterator begin,
typename std::vector<T>::iterator end)
const;
147 template<
class T>
typename std::vector<T>::const_iterator mymin(
const typename std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, T minConstraint)
const;
148 template<
class T>
typename std::vector<T>::iterator mymin(
const typename std::vector<T>& v,
typename std::vector<T>::iterator begin,
typename std::vector<T>::iterator end, T minConstraint)
const;
149 template<
class T>
typename std::vector<T>::const_iterator mymax(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end)
const;
150 template<
class T>
typename std::vector<T>::iterator mymax(
const std::vector<T>& v,
typename std::vector<T>::iterator begin,
typename std::vector<T>::iterator end)
const;
151 template<
class T>
typename std::vector<T>::const_iterator mymax(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, T maxConstraint)
const;
152 template<
class T>
typename std::vector<T>::iterator mymax(
const std::vector<T>& v,
typename std::vector<T>::iterator begin,
typename std::vector<T>::iterator end, T maxConstraint)
const;
153 template<
class T>
typename std::vector<T>::const_iterator absmin(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end)
const;
154 template<
class T>
typename std::vector<T>::const_iterator absmax(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end)
const;
156 template<
class T>
void minmax(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, T& theMin, T& theMax)
const;
157 template<
class T> T sum(
const std::vector<T>& v)
const;
158 template<
class T>
double mean(
const std::vector<T>& v)
const;
159 template<
class T>
void eraseNoData(std::vector<T>& v)
const;
160 template<
class T> T median(
const std::vector<T>& v)
const;
161 template<
class T>
double var(
const std::vector<T>& v)
const;
162 template<
class T>
double moment(
const std::vector<T>& v,
int n)
const;
163 template<
class T>
double cmoment(
const std::vector<T>& v,
int n)
const;
164 template<
class T>
double skewness(
const std::vector<T>& v)
const;
165 template<
class T>
double kurtosis(
const std::vector<T>& v)
const;
166 template<
class T>
void meanVar(
const std::vector<T>& v,
double& m1,
double& v1)
const;
167 template<
class T1,
class T2>
void scale2byte(
const std::vector<T1>& input, std::vector<T2>& output,
unsigned char lbound=0,
unsigned char ubound=255)
const;
168 template<
class T>
void distribution(
const std::vector<T>& input,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, std::vector<double>& output,
int nbin, T &minimum, T &maximum,
double sigma=0,
const std::string &filename=
"")
const;
169 template<
class T>
void distribution(
const std::vector<T>& input, std::vector<double>& output,
int nbin,
double sigma=0,
const std::string &filename=
"")
const{distribution(input,input.begin(),input.end(),output,nbin,0,0,sigma,filename);};
170 template<
class T>
void distribution2d(
const std::vector<T>& inputX,
const std::vector<T>& inputY, std::vector< std::vector<double> >& output,
int nbin, T& minX, T& maxX, T& minY, T& maxY,
double sigma=0,
const std::string& filename=
"")
const;
171 template<
class T>
void cumulative (
const std::vector<T>& input,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, std::vector<int>& output,
int nbin, T &minimum, T &maximum)
const;
172 template<
class T>
void percentiles (
const std::vector<T>& input,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, std::vector<T>& output,
int nbin, T &minimum, T &maximum,
const std::string &filename=
"")
const;
173 template<
class T>
void signature(
const std::vector<T>& input,
double& k,
double& alpha,
double& beta,
double e)
const;
174 void signature(
double m1,
double m2,
double& k,
double& alpha,
double& beta,
double e)
const;
175 template<
class T>
void normalize(
const std::vector<T>& input, std::vector<double>& output)
const;
176 template<
class T>
void normalize_pct(std::vector<T>& input)
const;
177 template<
class T>
double rmse(
const std::vector<T>& x,
const std::vector<T>& y)
const;
178 template<
class T>
double correlation(
const std::vector<T>& x,
const std::vector<T>& y,
int delay=0)
const;
179 template<
class T>
double cross_correlation(
const std::vector<T>& x,
const std::vector<T>& y,
int maxdelay, std::vector<T>& z)
const;
180 template<
class T>
double linear_regression(
const std::vector<T>& x,
const std::vector<T>& y,
double &c0,
double &c1)
const;
181 template<
class T>
double linear_regression_err(
const std::vector<T>& x,
const std::vector<T>& y,
double &c0,
double &c1)
const;
182 template<
class T>
void interpolateUp(
const std::vector<double>& wavelengthIn,
const std::vector<T>& input,
const std::vector<double>& wavelengthOut,
const std::string& type, std::vector<T>& output,
bool verbose=
false)
const;
183 template<
class T>
void interpolateUp(
const std::vector<double>& wavelengthIn,
const std::vector< std::vector<T> >& input,
const std::vector<double>& wavelengthOut,
const std::string& type, std::vector< std::vector<T> >& output,
bool verbose=
false)
const;
186 template<
class T>
void interpolateUp(
const std::vector<T>& input, std::vector<T>& output,
int nbin)
const;
187 template<
class T>
void nearUp(
const std::vector<T>& input, std::vector<T>& output)
const;
188 template<
class T>
void interpolateUp(
double* input,
int dim, std::vector<T>& output,
int nbin);
189 template<
class T>
void interpolateDown(
const std::vector<T>& input, std::vector<T>& output,
int nbin)
const;
190 template<
class T>
void interpolateDown(
double* input,
int dim, std::vector<T>& output,
int nbin);
193 static void initMap(std::map<std::string, INTERPOLATION_TYPE>& m_interpMap){
195 m_interpMap[
"linear"]=linear;
196 m_interpMap[
"polynomial"]=polynomial;
197 m_interpMap[
"cspline"]=cspline;
198 m_interpMap[
"cspline_periodic"]=cspline_periodic;
199 m_interpMap[
"akima"]=akima;
200 m_interpMap[
"akima_periodic"]=akima_periodic;
202 static void initDist(std::map<std::string, DISTRIBUTION_TYPE>& m_distMap){
204 m_distMap[
"gaussian"]=gaussian;
205 m_distMap[
"uniform"]=uniform;
207 std::vector<double> m_noDataValues;
211 template<
class T>
inline typename std::vector<T>::const_iterator StatFactory::mymin(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end)
const
213 typename std::vector<T>::const_iterator tmpIt=begin;
214 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
222 template<
class T>
inline typename std::vector<T>::iterator StatFactory::mymin(
const std::vector<T>& v,
typename std::vector<T>::iterator begin,
typename std::vector<T>::iterator end)
const
224 typename std::vector<T>::iterator tmpIt=begin;
225 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
233 template<
class T>
inline typename std::vector<T>::const_iterator StatFactory::mymin(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, T minConstraint)
const
235 typename std::vector<T>::const_iterator tmpIt=v.end();
236 T minValue=minConstraint;
237 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
240 if((minConstraint<=*it)&&(*it<=minValue)){
248 template<
class T>
inline typename std::vector<T>::iterator StatFactory::mymin(
const std::vector<T>& v,
typename std::vector<T>::iterator begin,
typename std::vector<T>::iterator end, T minConstraint)
const
250 typename std::vector<T>::iterator tmpIt=v.end();
251 T minValue=minConstraint;
252 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
255 if((minConstraint<=*it)&&(*it<=minValue)){
263 template<
class T>
inline typename std::vector<T>::const_iterator StatFactory::mymax(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end)
const
265 typename std::vector<T>::const_iterator tmpIt=begin;
266 for (
typename std::vector<T>::iterator it = begin; it!=end; ++it){
275 template<
class T>
inline typename std::vector<T>::iterator StatFactory::mymax(
const std::vector<T>& v,
typename std::vector<T>::iterator begin,
typename std::vector<T>::iterator end)
const
277 typename std::vector<T>::iterator tmpIt=begin;
278 for (
typename std::vector<T>::iterator it = begin; it!=end; ++it){
287 template<
class T>
inline typename std::vector<T>::const_iterator StatFactory::mymax(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, T maxConstraint)
const
289 typename std::vector<T>::const_iterator tmpIt=v.end();
290 T maxValue=maxConstraint;
291 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
294 if((maxValue<=*it)&&(*it<=maxConstraint)){
302 template<
class T>
inline typename std::vector<T>::iterator StatFactory::mymax(
const std::vector<T>& v,
typename std::vector<T>::iterator begin,
typename std::vector<T>::iterator end, T maxConstraint)
const
304 typename std::vector<T>::iterator tmpIt=v.end();
305 T maxValue=maxConstraint;
306 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
309 if((maxValue<=*it)&&(*it<=maxConstraint)){
320 template<
class T>
inline T StatFactory::mymin(
const std::vector<T>& v)
const
322 T minValue=*(v.begin());
323 for (
typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
332 template<
class T>
inline T StatFactory::mymin(
const std::vector<T>& v, T minConstraint)
const
334 T minValue=minConstraint;
335 for (
typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
336 if((minConstraint<=*it)&&(*it<=minValue))
342 template<
class T>
inline T StatFactory::mymax(
const std::vector<T>& v)
const
344 T maxValue=*(v.begin());
345 for (
typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
354 template<
class T>
inline T StatFactory::mymax(
const std::vector<T>& v, T maxConstraint)
const
356 T maxValue=maxConstraint;
357 for (
typename std::vector<T>::const_iterator it = v.begin(); it!=v.end(); ++it){
360 if((maxValue<=*it)&&(*it<=maxConstraint))
366 template<
class T>
inline typename std::vector<T>::const_iterator StatFactory::absmax(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end)
const
368 typename std::vector<T>::const_iterator tmpIt=begin;
369 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
372 if(abs(*tmpIt)<abs(*it))
378 template<
class T>
inline typename std::vector<T>::const_iterator StatFactory::absmin(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end)
const
380 typename std::vector<T>::const_iterator tmpIt=begin;
381 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
384 if(abs(*tmpIt)>abs(*it))
389 template<
class T>
inline void StatFactory::minmax(
const std::vector<T>& v,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, T& theMin, T& theMax)
const
393 for (
typename std::vector<T>::const_iterator it = begin; it!=end; ++it){
403 template<
class T>
inline T StatFactory::sum(
const std::vector<T>& v)
const
405 typename std::vector<T>::const_iterator it;
407 for (it = v.begin(); it!= v.end(); ++it){
415 template<
class T>
inline double StatFactory::mean(
const std::vector<T>& v)
const
418 typename std::vector<T>::const_iterator it;
420 unsigned int validSize=0;
421 for (it = v.begin(); it!= v.end(); ++it){
428 return static_cast<double>(tmpSum)/validSize;
433 template<
class T>
inline void StatFactory::eraseNoData(std::vector<T>& v)
const
435 typename std::vector<T>::iterator it=v.begin();
444 template<
class T> T StatFactory::median(
const std::vector<T>& v)
const
446 std::vector<T> tmpV=v;
448 sort(tmpV.begin(),tmpV.end());
450 return tmpV[tmpV.size()/2];
452 return 0.5*(tmpV[tmpV.size()/2-1]+tmpV[tmpV.size()/2]);
455 template<
class T>
double StatFactory::var(
const std::vector<T>& v)
const
457 typename std::vector<T>::const_iterator it;
458 unsigned int validSize=0;
461 for (it = v.begin(); it!= v.end(); ++it){
486 template<
class T>
double StatFactory::moment(
const std::vector<T>& v,
int n)
const
489 unsigned int validSize=0;
490 typename std::vector<T>::const_iterator it;
493 for(it = v.begin(); it!= v.end(); ++it){
507 template<
class T>
double StatFactory::cmoment(
const std::vector<T>& v,
int n)
const
510 unsigned int validSize=0;
511 typename std::vector<T>::const_iterator it;
514 for(it = v.begin(); it!= v.end(); ++it){
524 template<
class T>
double StatFactory::skewness(
const std::vector<T>& v)
const
526 return cmoment(v,3)/pow(var(v),1.5);
529 template<
class T>
double StatFactory::kurtosis(
const std::vector<T>& v)
const
531 double m2=cmoment(v,2);
532 double m4=cmoment(v,4);
536 template<
class T>
void StatFactory::meanVar(
const std::vector<T>& v,
double& m1,
double& v1)
const
538 typename std::vector<T>::const_iterator it;
539 unsigned int validSize=0;
543 for (it = v.begin(); it!= v.end(); ++it){
566 template<
class T1,
class T2>
void StatFactory::scale2byte(
const std::vector<T1>& input, std::vector<T2>& output,
unsigned char lbound,
unsigned char ubound)
const
568 output.resize(input.size());
569 T1 minimum=mymin(input);
570 T1 maximum=mymax(input);
571 assert(maximum>minimum);
572 double scale=(ubound-lbound)/(maximum-minimum);
573 for (
int i=0;i<input.size();++i)
574 output[i]=scale*(input[i]-(minimum))+lbound;
577 template<
class T>
void StatFactory::distribution(
const std::vector<T>& input,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, std::vector<double>& output,
int nbin, T &minimum, T &maximum,
double sigma,
const std::string &filename)
const
581 minmax(input,begin,end,minValue,maxValue);
582 if(minimum<maximum&&minimum>minValue)
584 if(minimum<maximum&&maximum<maxValue)
591 if(maximum<=minimum){
592 std::ostringstream s;
593 s<<
"Error: could not calculate distribution (min>=max)";
597 assert(input.size());
598 if(output.size()!=nbin){
600 for(
int i=0;i<nbin;output[i++]=0);
602 typename std::vector<T>::const_iterator it;
603 for(it=begin;it!=end;++it){
616 for(
int ibin=0;ibin<nbin;++ibin){
617 double icenter=minimum+
static_cast<double>(maximum-minimum)*(ibin+0.5)/nbin;
618 double thePdf=gsl_ran_gaussian_pdf(*it-icenter, sigma);
619 output[ibin]+=thePdf;
626 else if(*it>minimum && *it<maximum)
627 theBin=
static_cast<int>(
static_cast<double>((nbin-1)*(*it)-minimum)/(maximum-minimum));
635 if(!filename.empty()){
636 std::ofstream outputfile;
637 outputfile.open(filename.c_str());
639 std::ostringstream s;
640 s<<
"Error opening distribution file , " << filename;
643 for(
int ibin=0;ibin<nbin;++ibin)
644 outputfile << minimum+static_cast<double>(maximum-minimum)*(ibin+0.5)/nbin <<
" " <<
static_cast<double>(output[ibin])/input.size() << std::endl;
649 template<
class T>
void StatFactory::distribution2d(
const std::vector<T>& inputX,
const std::vector<T>& inputY, std::vector< std::vector<double> >& output,
int nbin, T& minX, T& maxX, T& minY, T& maxY,
double sigma,
const std::string& filename)
const
651 assert(inputX.size());
652 assert(inputX.size()==inputY.size());
653 unsigned long int npoint=inputX.size();
655 minmax(inputX,inputX.begin(),inputX.end(),minX,maxX);
657 std::ostringstream s;
658 s<<
"Error: could not calculate distribution (minX>=maxX)";
662 minmax(inputY,inputY.begin(),inputY.end(),minY,maxY);
664 std::ostringstream s;
665 s<<
"Error: could not calculate distribution (minY>=maxY)";
670 for(
int i=0;i<nbin;++i){
671 output[i].resize(nbin);
672 for(
int j=0;j<nbin;++j)
677 for(
unsigned long int ipoint=0;ipoint<npoint;++ipoint){
678 if(inputX[ipoint]==maxX)
681 binX=
static_cast<int>(
static_cast<double>(inputX[ipoint]-minX)/(maxX-minX)*nbin);
682 if(inputY[ipoint]==maxY)
685 binY=
static_cast<int>(
static_cast<double>(inputY[ipoint]-minY)/(maxY-minY)*nbin);
687 assert(binX<output.size());
689 assert(binY<output[binX].size());
697 for(
int ibinX=0;ibinX<nbin;++ibinX){
698 double centerX=minX+
static_cast<double>(maxX-minX)*ibinX/nbin;
699 double pdfX=gsl_ran_gaussian_pdf(inputX[ipoint]-centerX, sigma);
700 for(
int ibinY=0;ibinY<nbin;++ibinY){
702 double centerY=minY+
static_cast<double>(maxY-minY)*ibinY/nbin;
703 double pdfY=gsl_ran_gaussian_pdf(inputY[ipoint]-centerY, sigma);
704 output[ibinX][binY]+=pdfX*pdfY;
709 ++output[binX][binY];
711 if(!filename.empty()){
712 std::ofstream outputfile;
713 outputfile.open(filename.c_str());
715 std::ostringstream s;
716 s<<
"Error opening distribution file , " << filename;
719 for(
int bin=0;bin<nbin;++bin){
720 for(
int bin=0;bin<nbin;++bin){
721 double value=
static_cast<double>(output[binX][binY])/npoint;
722 outputfile << (maxX-minX)*bin/(nbin-1)+minX <<
" " << (maxY-minY)*bin/(nbin-1)+minY <<
" " << value << std::endl;
729 template<
class T>
void StatFactory::percentiles (
const std::vector<T>& input,
typename std::vector<T>::const_iterator begin,
typename std::vector<T>::const_iterator end, std::vector<T>& output,
int nbin, T &minimum, T &maximum,
const std::string &filename)
const
732 minmax(input,begin,end,minimum,maximum);
737 assert(maximum>minimum);
739 assert(input.size());
741 std::vector<T> inputSort;
742 inputSort.assign(begin,end);
743 typename std::vector<T>::iterator vit=inputSort.begin();
744 while(vit!=inputSort.end()){
745 if(*vit<minimum||*vit>maximum)
746 inputSort.erase(vit);
750 std::sort(inputSort.begin(),inputSort.end());
751 vit=inputSort.begin();
752 std::vector<T> inputBin;
753 for(
int ibin=0;ibin<nbin;++ibin){
755 while(inputBin.size()<inputSort.size()/nbin&&vit!=inputSort.end()){
756 inputBin.push_back(*vit);
761 output[ibin]=mymax(inputBin);
764 if(!filename.empty()){
765 std::ofstream outputfile;
766 outputfile.open(filename.c_str());
768 std::ostringstream s;
769 s<<
"error opening distribution file , " << filename;
772 for(
int ibin=0;ibin<nbin;++ibin)
773 outputfile << ibin*100.0/nbin <<
" " << static_cast<double>(output[ibin])/input.size() << std::endl;
799 template<
class T>
void StatFactory::signature(
const std::vector<T>& input,
double&k,
double& alpha,
double& beta,
double e)
const
801 double m1=moment(input,1);
802 double m2=moment(input,2);
803 signature(m1,m2,k,alpha,beta,e);
806 template<
class T>
void StatFactory::normalize(
const std::vector<T>& input, std::vector<double>& output)
const{
807 double total=sum(input);
809 output.resize(input.size());
810 for(
int index=0;index<input.size();++index)
811 output[index]=input[index]/total;
817 template<
class T>
void StatFactory::normalize_pct(std::vector<T>& input)
const{
818 double total=sum(input);
820 typename std::vector<T>::iterator it;
821 for(it=input.begin();it!=input.end();++it)
822 *it=100.0*(*it)/total;
826 template<
class T>
double StatFactory::rmse(
const std::vector<T>& x,
const std::vector<T>& y)
const{
827 assert(x.size()==y.size());
830 for(
int isample=0;isample<x.size();++isample){
831 double e=x[isample]-y[isample];
837 template<
class T>
double StatFactory::correlation(
const std::vector<T>& x,
const std::vector<T>& y,
int delay)
const{
843 meanVar(x,meanX,varX);
844 meanVar(y,meanY,varY);
845 double denom = sqrt(varX*varY);
849 for (
int i=0;i<x.size();++i) {
851 if (j < 0 || j >= y.size())
854 assert(i>=0&&i<x.size());
855 assert(j>=0&&j<y.size());
856 sXY += (x[i] - meanX) * (y[j] - meanY);
859 double minSize=(x.size()<y.size())?x.size():y.size();
860 return(sXY / denom / (minSize-1));
866 template<
class T>
double StatFactory::cross_correlation(
const std::vector<T>& x,
const std::vector<T>& y,
int maxdelay, std::vector<T>& z)
const{
868 double sumCorrelation=0;
869 for (
int delay=-maxdelay;delay<maxdelay;delay++) {
870 z.push_back(correlation(x,y,delay));
871 sumCorrelation+=z.back();
873 return sumCorrelation;
876 template<
class T>
double StatFactory::linear_regression(
const std::vector<T>& x,
const std::vector<T>& y,
double &c0,
double &c1)
const{
877 assert(x.size()==y.size());
883 gsl_fit_linear(&(x[0]),1,&(y[0]),1,x.size(),&c0,&c1,&cov00,&cov01,&cov11,&sumsq);
884 return (1-sumsq/var(y)/(y.size()-1));
887 template<
class T>
double StatFactory::linear_regression_err(
const std::vector<T>& x,
const std::vector<T>& y,
double &c0,
double &c1)
const{
888 assert(x.size()==y.size());
894 gsl_fit_linear(&(x[0]),1,&(y[0]),1,x.size(),&c0,&c1,&cov00,&cov01,&cov11,&sumsq);
895 return sqrt((sumsq)/(y.size()));
901 template<
class T>
void StatFactory::interpolateUp(
const std::vector<double>& wavelengthIn,
const std::vector<T>& input,
const std::vector<double>& wavelengthOut,
const std::string& type, std::vector<T>& output,
bool verbose)
const{
902 assert(wavelengthIn.size());
903 assert(input.size()==wavelengthIn.size());
904 assert(wavelengthOut.size());
905 int nband=wavelengthIn.size();
907 gsl_interp_accel *acc;
910 getSpline(type,nband,spline);
912 assert(&(wavelengthIn[0]));
914 initSpline(spline,&(wavelengthIn[0]),&(input[0]),nband);
915 for(
int index=0;index<wavelengthOut.size();++index){
916 if(wavelengthOut[index]<*wavelengthIn.begin()){
917 output.push_back(*(input.begin()));
920 else if(wavelengthOut[index]>wavelengthIn.back()){
921 output.push_back(input.back());
934 double dout=evalSpline(spline,wavelengthOut[index],acc);
935 output.push_back(dout);
937 gsl_spline_free(spline);
938 gsl_interp_accel_free(acc);
972 template<
class T>
void StatFactory::interpolateUp(
const std::vector<T>& input, std::vector<T>& output,
int nbin)
const
974 assert(input.size());
977 int dim=input.size();
978 for(
int i=0;i<dim;++i){
980 double left=input[i];
982 double right=(i<dim-1)? input[i+1]:input[i];
983 deltaX=(right-left)/static_cast<double>(nbin);
984 for(
int x=0;x<nbin;++x){
985 output.push_back(left+x*deltaX);
989 output.push_back(input.back());
993 template<
class T>
void StatFactory::nearUp(
const std::vector<T>& input, std::vector<T>& output)
const
995 assert(input.size());
996 assert(output.size()>=input.size());
997 int dimInput=input.size();
998 int dimOutput=output.size();
1000 for(
int iin=0;iin<dimInput;++iin){
1001 for(
int iout=0;iout<dimOutput/dimInput;++iout){
1002 int indexOutput=iin*dimOutput/dimInput+iout;
1003 assert(indexOutput<output.size());
1004 output[indexOutput]=input[iin];
1009 template<
class T>
void StatFactory::interpolateUp(
double* input,
int dim, std::vector<T>& output,
int nbin)
1013 for(
int i=0;i<dim;++i){
1015 double left=input[i];
1017 double right=(i<dim-1)? input[i+1]:input[i];
1018 deltaX=(right-left)/static_cast<double>(nbin);
1019 for(
int x=0;x<nbin;++x){
1020 output.push_back(left+x*deltaX);
1024 output.push_back(input[dim-1]);
1028 template<
class T>
void StatFactory::interpolateDown(
const std::vector<T>& input, std::vector<T>& output,
int nbin)
const
1030 assert(input.size());
1033 int dim=input.size();
1035 output.push_back(input[0]);
1036 for(
int i=1;i<dim;++i){
1041 output.push_back(input[i]);
1046 template<
class T>
void StatFactory::interpolateDown(
double* input,
int dim, std::vector<T>& output,
int nbin)
1051 output.push_back(input[0]);
1052 for(
int i=1;i<dim;++i){
1057 output.push_back(input[i]);