statistics.h

/***************************************************************************
                          statistics.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2001 by Thies Jochimsen
    email                : jochimse@cns.mpg.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) {
  Data<T,N_rank> ensemble_copy(ensemble);

  Array<T,N_rank> diff(ensemble.shape());

  TinyVector<int,N_rank> index;
  for(int i=0; i<ensemble.numElements(); i++) {
    index=ensemble_copy.create_index(i);
    diff(index)=sum(fabs(ensemble-ensemble(index)));
  }

  return ensemble(minIndex(diff));
}


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(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;
  }
  result/=weightsum;

  return result;
}


#endif

---