TruncMean.cxx

Go to the documentation of this file.

#include "TruncMean.h"

#include <algorithm>
#include <cmath>
#include <vector>

float TruncMean::CalcIterativeTruncMean(std::vector<float> v,
                                        const size_t& nmin,
                                        const size_t& nmax,
                                        const size_t& currentiteration,
                                        const size_t& lmin,
                                        const float& convergencelimit,
                                        const float& nsigma,
                                        const float& oldmed)
{

  auto const& mean = Mean(v);
  auto const& med = Median(v);
  auto const& rms = RMS(v);

  // if the vector length is below the lower limit -> return
  if (v.size() < lmin) return mean;

  // if we have passed the maximum number of iterations -> return
  if (currentiteration >= nmax) return mean;

  // if we passed the minimum number of iterations and the mean is close enough to the old value
  float fracdiff = fabs(med - oldmed) / oldmed;
  if ((currentiteration >= nmin) && (fracdiff < convergencelimit)) return mean;

  // if reached here it means we have to go on for another iteration

  // cutoff tails of distribution surrounding the mean
  // use erase-remove : https://en.wikipedia.org/wiki/Erase%E2%80%93remove_idiom
  // https://stackoverflow.com/questions/17270837/stdvector-removing-elements-which-fulfill-some-conditions
  v.erase(std::remove_if(v.begin(),
                         v.end(),
                         [med, nsigma, rms](const float& x) {
                           return ((x < (med - nsigma * rms)) || (x > (med + nsigma * rms)));
                         }), // lamdda condition for events to be removed
          v.end());

  return CalcIterativeTruncMean(
    v, nmin, nmax, lmin, currentiteration + 1, convergencelimit, nsigma, med);
}

void TruncMean::CalcTruncMeanProfile(const std::vector<float>& rr_v,
                                     const std::vector<float>& dq_v,
                                     std::vector<float>& dq_trunc_v,
                                     const float& nsigma)
{

  // how many points to sample
  int Nneighbor = (int)(_rad * 3 * 2);

  dq_trunc_v.clear();
  dq_trunc_v.reserve(rr_v.size());

  int Nmax = dq_v.size() - 1;

  for (size_t n = 0; n < dq_v.size(); n++) {

    // current residual range
    float rr = rr_v.at(n);

    int nmin = n - Nneighbor;
    int nmax = n + Nneighbor;

    if (nmin < 0) nmin = 0;
    if (nmax > Nmax) nmax = Nmax;

    // vector for local dq values
    std::vector<float> dq_local_v;

    for (int i = nmin; i < nmax; i++) {

      float dr = rr - rr_v[i];
      if (dr < 0) dr *= -1;

      if (dr > _rad) continue;

      dq_local_v.push_back(dq_v[i]);

    } // for all ticks we want to scan

    if (dq_local_v.size() == 0) {
      dq_trunc_v.push_back(dq_v.at(n));
      continue;
    }

    // calculate median and rms
    float median = Median(dq_local_v);
    float rms = RMS(dq_local_v);

    float truncated_dq = 0.;
    int npts = 0;
    for (auto const& dq : dq_local_v) {
      if ((dq < (median + rms * nsigma)) && (dq > (median - rms * nsigma))) {
        truncated_dq += dq;
        npts += 1;
      }
    }

    dq_trunc_v.push_back(truncated_dq / npts);
  } // for all values

  return;
}

float TruncMean::Mean(const std::vector<float>& v)
{

  float mean = 0.;
  for (auto const& n : v)
    mean += n;
  mean /= v.size();

  return mean;
}

float TruncMean::Median(const std::vector<float>& v)
{

  if (v.size() == 1) return v[0];

  std::vector<float> vcpy = v;

  std::sort(vcpy.begin(), vcpy.end());

  float median = vcpy[vcpy.size() / 2];

  return median;
}

float TruncMean::RMS(const std::vector<float>& v)
{

  float avg = 0.;
  for (auto const& val : v)
    avg += val;
  avg /= v.size();
  float rms = 0.;
  for (auto const& val : v)
    rms += (val - avg) * (val - avg);
  rms = sqrt(rms / (v.size() - 1));

  return rms;
}