/home/pascal/depot/filedata/src/statistic.h

Go to the documentation of this file.

/***************************************************************************
 *   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 <numeric>
#include <algorithm>

template<class T, class Iter_T>
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;
}

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;

};
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);
};

template<class T1, class Iter_T1,class T2, class Iter_T2>
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<T1>());
  else
    // v1 sorted in decreasing order
    up_bound=upper_bound(v1_start,v1_end,val1,greater<T1>());

  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)));
}

template<class Iter_T1,class Iter_T2>
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<typename NumType>
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<NumType*>(v_start),const_cast<NumType*>(v_end),val,less<NumType>());
  else
    low=lower_bound(const_cast<NumType*>(v_start),const_cast<NumType*>(v_end),val,greater<NumType>());
  /* 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<upper_val_dist) return low_idx-1;
    else return low_idx;
  }
}

#endif