libodin-dev 1.8.8-2ubuntu1 / usr / include / odindata / statistics.h

/***************************************************************************
                          statistics.h  -  description
                             -------------------
    begin                : Fri Apr 6 2001
    copyright            : (C) 2000-2014 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>


/**
  * @addtogroup odindata
  * @{
  */

///////////////////////////////////////////////////////////

/**
  * Results of a statistical analysis
  *
  */
struct statisticResult {

/**
  * The minimum value
  */
  double min;

/**
  * The minimum value
  */
  double max;

/**
  * The mean value
  */
  double mean;

/**
  * The standard deviation
  */
  double stdev;

/**
  * The standard deviation of the mean value
  */
  double meandev;

/**
  *  Prints statistics to ostream
  */
  friend STD_ostream& operator << (STD_ostream& s, statisticResult stats);

};


///////////////////////////////////////////////////////////

/**
  * Returns a statistical description of 'ensemble'. If mask is non-zero,
  * only values with mask!=0 are considered.
  *
  */
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;
}

///////////////////////////////////////////////////////////

/**
  * Returns the median of 'ensemble'.
  *
  */
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;
}

///////////////////////////////////////////////////////////

/**
  * Returns the weighted mean of 'ensemble' using 'weight'
  * for the weighting.
  *
  */
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;
  }
  if(weightsum) result/=weightsum;

  return result;
}


/** @}
  */

#endif