statistics.h

/***************************************************************************
                          statistics.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2000-2021 by Thies Jochimsen
    email                : thies@jochimsen.de
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef STATISTICS_H
#define STATISTICS_H
 
#include<odindata/data.h>
#include<odindata/utils.h>
 
 
 
struct statisticResult {
 
  double min;
 
  double max;
 
  double mean;
 
  double stdev;
 
  double meandev;
 
  friend STD_ostream& operator << (STD_ostream& s, statisticResult stats);
 
};
 
 
 
template<typename T, int N_rank>
statisticResult statistics(const Array<T,N_rank>& ensemble, const Array<T,N_rank>* mask=0) {
  Log<OdinData> odinlog("","statistics");
  statisticResult result;
 
  result.min=0.0;
  result.max=0.0;
  result.mean=0.0;
  result.stdev=0.0;
  result.meandev=0.0;
 
 
  if(mask && !same_shape(ensemble,*mask)) {
    ODINLOG(odinlog,errorLog) << "size mismatch (ensemble.shape()=" << ensemble.shape() << ") != (mask.shape()=" << mask->shape() << ")" << STD_endl;
    return result;
  }
 
 
  int n=ensemble.numElements();
 
  Data<T,N_rank> ensemble_copy(ensemble);
 
  int nvals=0;
  TinyVector<int,N_rank> index;
  for(int i=0; i<n; i++) {
    index=ensemble_copy.create_index(i);
    bool include=true;
    if(mask && (*mask)(index)==0.0) include=false;
    if(include) {
      float val=ensemble(index);
      result.mean+=val;
      nvals++;
      if(i==0) {
        result.min=result.max=val;
      } else {
        if(val<result.min) result.min=val;
        if(val>result.max) result.max=val;
      }
    }
  }
  result.mean=secureDivision(result.mean,double(nvals));
 
  nvals=0;
  for(int i=0; i<n; i++) {
    index=ensemble_copy.create_index(i);
    bool include=true;
    if(mask && (*mask)(index)==0.0) include=false;
    if(include) {
      double diff=result.mean-ensemble(index);
      result.stdev+=diff*diff;
      nvals++;
    }
  }
  if(nvals>1) result.stdev=sqrt(result.stdev/double(nvals-1));
  else result.stdev=0.0;
 
  result.meandev=result.stdev/sqrt(double(nvals));
 
  return result;
}
 
 
template<typename T, int N_rank>
T median(const Array<T,N_rank>& ensemble, const Array<T,N_rank>* mask=0) {
  T result(0);
  Data<T,N_rank> ensemble_copy(ensemble);
  STD_list<T> vallist;
  int n=ensemble_copy.numElements();
  if(!n) return result;
  TinyVector<int,N_rank> index;
  for(int i=0; i<n; i++) {
    index=ensemble_copy.create_index(i);
    bool include=true;
    if(mask && (*mask)(index)==0.0) include=false;
    if(include) vallist.push_back(ensemble_copy(index));
  }
  vallist.sort();
 
  STD_vector<T> vec(list2vector(vallist));
 
  if(n%2) {
    result=vec[(n-1)/2]; // odd number of elements
  } else {
    result=0.5*(vec[n/2]+vec[n/2-1]); // even number of elements
  }
 
  return result;
}
 
 
template<typename T, int N_rank>
T weightmean(const Array<T,N_rank>& ensemble, const Array<float,N_rank>& weight) {
  Log<OdinData> odinlog("","weightmean");
  T result;
  result=0; // double-stage init for TinyVector
 
  Data<T,N_rank> ensemble_copy(ensemble);
 
  if(ensemble.shape()!=weight.shape()) {
    ODINLOG(odinlog,errorLog) << "size mismatch (ensemble.shape()=" << ensemble.shape() << ") != (weight.shape()=" << weight.shape() << ")" << STD_endl;
    return result;
  }
 
  float weightsum=0.0;
  TinyVector<int,N_rank> index;
  for(unsigned int i=0; i<ensemble_copy.numElements(); i++) {
    index=ensemble_copy.create_index(i);
    float w=fabs(weight(index));
    result+=w*ensemble_copy(index);
    weightsum+=w;
  }
  if(weightsum) result/=weightsum;
 
  return result;
}
 
 
#endif