/*************************************************************************** * Copyright (C) 2006 by Nicolas PASCAL * * pascal@icare-pc12 * * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the * * Free Software Foundation, Inc., * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ #ifndef STATISTIC_H #define STATISTIC_H #include #include /** * @file statistic.h contain some statistics computing methods */ /** * @brief compute the mean value of a vector * @param first the first element of the vector to consider * @param last the last one * @return the computed mean */ template inline T mean(Iter_T first, Iter_T last) { T sum=accumulate(first,last,T()); size_t nb_val=distance(first,last); return sum/nb_val; } /** * Make a generic bilinear interpolation * @param x * @param y * @param x1 * @param x2 * @param y1 * @param y2 * @param f11 f(x1,y1) * @param f12 f(x1,y2) * @param f21 f(x2,y1) * @param f22 f(x2,y2) * @return f(x,y) */ inline double bilinear_interpolation( const double &x, const double &y, const double &x1, const double &y1, const double &x2, const double &y2, const double &f11, const double &f12, const double &f21, const double &f22) { double u = (x-x1)/(x2-x1); double v = (y-y1)/(y2-y1); return (1.-u)*(1.-v)*f11 + (1.-u)*(v)*f12 + (u)*(1.-v)*f21 + (u)*(v)*f22; }; /** * Make a generic linear interpolation * @param x * @param x1 * @param x2 * @param f1 f(x1) * @param f2 f(x2) * @return f(x) */ inline double linear_interpolation( const double x, const double x1, const double x2, const double f1, const double f2 ) { return f1+(x-x1)*(f2-f1)/(x2-x1); }; /** * @brief linear interpolation the @a val1 value * This method supports either classic C arrays and STL's vector * @param v1 the input abscissa values. Must be sorted either in ascending or descending order. Doesn't support STL vectors * @param v2 the output abscissa values * @param val1 the input value (in v1 coordinates) * @param val2 the value interpolated in v2 corresponding to val1 */ template inline void interpolate( Iter_T1 v1_start, Iter_T1 v1_end, Iter_T2 v2_start, Iter_T2 v2_end, T1 &val1, T2 &val2) { int size_v1=distance(v1_start, v1_end); int size_v2=distance(v2_start, v2_end); //--- Check the input parameters validity if (size_v1<2 || size_v2<2 || size_v1!=size_v2) { string msg="the abscissa vectors do not have a valid size : size(v1)="+MyTools::to_string(size_v1)+" size(v2)="+MyTools::to_string(size_v2); InterpolationError e(__FILE__,__LINE__,msg); throw e; } //--- find the 2 values in v1 that "bounds" val1. If val1 is less that the min value of v1 or gretaer than the max value of v1, an extrapolatino either of an interpolation int v1_up_bound_idx=-1,v1_low_bound_idx=-1; // indexes of the values that bound val1 in v1 Iter_T1 up_bound=(Iter_T1)(0); if (*(v1_start)<*(v1_start+1)) // v1 sorted in increasing order up_bound=upper_bound(v1_start,v1_end,val1,less()); else // v1 sorted in decreasing order up_bound=upper_bound(v1_start,v1_end,val1,greater()); if (up_bound==v1_start) { // val1 is lower than the min value of v1 -> extrapolate v1_up_bound_idx=1; v1_low_bound_idx=0; } else if (up_bound==v1_end) { // val1 is upper than the max value of v1 -> extrapolate v1_up_bound_idx=size_v1-1; v1_low_bound_idx=size_v1-2; } else { // general case v1_up_bound_idx=distance(v1_start,up_bound); v1_low_bound_idx=v1_up_bound_idx-1; } //--- find the 2 values that correspond to v1_upper_val and v1_lower_val in v2 T1 v1_upper_val=*(v1_start+v1_up_bound_idx); T1 v1_lower_val=*(v1_start+v1_low_bound_idx); T2 v2_upper_val=*(v2_start+v1_up_bound_idx); T2 v2_lower_val=*(v2_start+v1_low_bound_idx); //--- do the interpolation val2=(T2)(linear_interpolation((double)(val1),(double)(v1_lower_val),(double)(v1_upper_val),(double)(v2_lower_val),(double)(v2_upper_val))); } /** * @brief linear interpolation of the values of a vector. * This method supports either classic C arrays and STL's vector * @param v1_start start of input abscissa values. Must be sorted either in ascending or descending order * @param v1_end end of input abscissa values * @param v2_start start of output abscissa values * @param v2_end end of output abscissa values * @param vin_start start of the input data to interpolate * @param vin_end end of the input data to interpolate * @param vout_start start of the interpolated output data */ template inline void interpolate( Iter_T1 v1_start, Iter_T1 v1_end, Iter_T2 v2_start, Iter_T2 v2_end, Iter_T1 vin_start, Iter_T1 vin_end, Iter_T2 vout_start) { for (Iter_T1 p_inval=vin_start;p_inval!=vin_end;++p_inval ) // interpolate all values of vin and store them in vout interpolate(v1_start, v1_end, v2_start, v2_end,*p_inval,*(vout_start+distance(vin_start,p_inval))); } template /** * @brief search the value that is the nearest to @a val in a sorted vector * If @a val is situed exactly between 2 values, the upper index is returned * @warning The method does not support vectors that contain equal value. In other mathematical words, the vector has to be in strict monotonic order * This method supports either classic C arrays (v_start and v_end will be pointers) and STL's vector (use iterators in the case) * @param v_start vector start position (pointer or iterator) * @param v_end vector end position (pointer that points just after the last value to check or iterator) * @param val the value we are searching the nearest position * @return the index of the nearest value in the vector. In range [0,vector_size-1]. Return -1 in case of problem : empty vector, v_start or v_end NULL */ long get_nearest_val_idx( const NumType* v_start, const NumType* v_end, const NumType &val) { if (v_start==NULL || v_end==NULL || v_end-v_start<=1) return -1; int nb_val=v_end-v_start; // number of values of the vector /* if the vector has only one value, the nearest one is at index 0 */ if (nb_val==0) return 0; assert(nb_val>1); /* search for the insertion position of val */ NumType* low=NULL; if (*(v_start+1)>*(v_start)) // vector is sorted in ascending order low=lower_bound(const_cast(v_start),const_cast(v_end),val,less()); else low=lower_bound(const_cast(v_start),const_cast(v_end),val,greater()); /* search for the nearest value */ int low_idx=low-v_start; // index of the lower bound if (low_idx==0) return 0; else if (low_idx==(nb_val)) return nb_val-1; else { // compute distance to the bounding values NumType lower_val=*(v_start+low_idx-1); NumType lower_val_dist=fabs(val-lower_val); NumType upper_val=*(v_start+low_idx); NumType upper_val_dist=fabs(val-upper_val); if (lower_val_dist