/**
   \file Accumulator.h

   \author Darko Veberic
   \date 24 Aug 2008

   \version $Id$
*/

#ifndef _utl_Accumulator_h_
#define _utl_Accumulator_h_

#include <vector>
#include <deque>
#include <utl/Math.h>
#include <utl/SafeBoolCast.h>

#include <boost/shared_ptr.hpp>

#include <algorithm>


namespace utl {
  namespace Accumulator {

    template<typename T>
    class Min {
    public:
      typedef T Type;
      Min(const T first) : fMin(first) { }
      void operator()(const T x) { if (fMin > x) fMin = x; }
      T GetMin() const { return fMin; }
      void Clear(const T x) { fMin = x; }
    protected:
      Min() : fMin(T()) { }
      void SetFirst(const T x) { fMin = x; }
      void Clear() { fMin = T(); }
    private:
      T fMin;
    };


    template<typename T>
    class Max {
    public:
      typedef T Type;
      Max(const T first) : fMax(first) { }
      void operator()(const T x) { if (fMax < x) fMax = x; }
      T GetMax() const { return fMax; }
      void Clear(const T x) { fMax = x; }
    protected:
      Max() : fMax(T()) { }
      void Clear() { fMax = T(); }
    private:
      T fMax;
    };


    template<typename T>
    class MinMax : public Min<T>, public Max<T> {
    public:
      typedef T Type;
      MinMax(const T first) : Min<T>(first), Max<T>(first) { }
      MinMax(const T firstMin, const T firstMax) : Min<T>(firstMin), Max<T>(firstMax) { }
      void operator()(const T x) { Min<T>::operator()(x); Max<T>::operator()(x); }
      void Clear(const T x) { Clear(x, x); }
      void Clear(const T min, const T max) { Min<T>::Clear(min); Max<T>::Clear(max); }
    protected:
      MinMax() { }
      void Clear() { Min<T>::Clear(); Max<T>::Clear(); }
    };


    class AverageN {
    public:
      AverageN() : fX(0) { }
      void operator()(const double x) { fX += x; }
      double GetAverage(const int n) const { return fX / n; }
      double GetSum() const { return fX; }
      void Clear() { fX = 0; }
    protected:
      double fX;
    };


    class Average {
    public:
      Average() : fAvg(), fN(0) { }
      void operator()(const double x) { fAvg(x); ++fN; }
      double GetAverage() const { return fAvg.GetAverage(fN); }
      int GetN() const { return fN; }
      void Clear() { fAvg.Clear(); fN = 0; }
    private:
      AverageN fAvg;
      int fN;
    };


    template<typename T>
    class MinMaxAverage : public Min<T>, public Max<T>, public Average {
    public:
      typedef T Type;
      MinMaxAverage(const T first) : Min<T>(first), Max<T>(first) { Average::operator()(first); }
      void operator()(const T x) { Min<T>::operator()(x); Max<T>::operator()(x); Average::operator()(x); }
      void Clear(const T x) { Min<T>::Clear(x); Max<T>::Clear(x); Average::Clear(); Average::operator()(x); }
    protected:
      MinMaxAverage() { }
      void Clear() { Min<T>::Clear(); Max<T>::Clear(); Average::Clear(); }
    };


    //


    class Validator {
    public:
      Validator() : fIsValid(false) { }
      bool IsValid() const { return fIsValid; }
      void Validate(const bool state = true) { fIsValid = state; }
      void Clear() { fIsValid = false; }
    private:
      bool fIsValid;
    };

    /**
       This class enables the usage of the uninitialized versions
       (default empty constructors) of the accumulators that require
       the first value to be passed in the constructor.

       Example:
       \code
       Accumulator::Safe<Accumulator::Min<int> > min;

       if (min)
         cout << "min is initialized" << endl;
       else
         cout << "min is uninitialized" << endl;

       min(13);

       if (min)
         cout << "min is initialized" << endl;
       else
         cout << "min is uninitialized" << endl;

       cout << min.GetMin() << endl;
       \endcode

       This example prints
       \code
       min is uninitialized
       min is initialized
       13
       \endcode
    */
    template<typename AccumulatorT>
    class Safe : public AccumulatorT, public Validator,
                 public SafeBoolCast<Safe<AccumulatorT> > {
    public:
      void
      operator()(const typename AccumulatorT::Type x)
      {
        if (*this)
          AccumulatorT::operator()(x);
        else {
          AccumulatorT::Clear(x);
          Validate();
        }
      }

      bool BoolCast() const { return IsValid(); }

      void Clear() { AccumulatorT::Clear(); Validator::Clear(); }
    };


    //


    class SampleVarianceN : public AverageN {
    public:
      SampleVarianceN() : AverageN(), fX2(0) { }
      void operator()(const double x) { AverageN::operator()(x); fX2 += Sqr(x); }
      double GetVariance(const int n) const { return (fX2 - Sqr(fX)/n) / (n - 1); }
      double GetAverageError(const int n) const { return std::sqrt(GetVariance(n)/n); }
      double GetSumOfSquares() const { return fX2; }
      void Clear() { AverageN::Clear(); fX2 = 0; }
    protected:
      double fX2;
    };


    //


    class SampleStandardDeviationN : public SampleVarianceN {
    public:
      SampleStandardDeviationN() : SampleVarianceN() { }
      void operator()(const double x) { SampleVarianceN::operator()(x); }
      double GetStandardDeviation(const int n) const { return std::sqrt(GetVariance(n)); }
    };


    //


    class SampleVariance : public Average {
    public:
      SampleVariance() : Average(), fVariance() { }
      void operator()(const double x) { Average::operator()(x); fVariance(x); }
      double GetVariance() const { return fVariance.GetVariance(Average::GetN()); }
      double GetAverageError() const { return fVariance.GetAverageError(GetN()); }
      double GetSumOfSquares() const { return fVariance.GetSumOfSquares(); }
      void Clear() { Average::Clear(); fVariance.Clear(); }
    protected:
      SampleVarianceN fVariance;
    };


    /**
       \class SampleStandardDeviation Accumulator.h "utl/Accumulator.h"

       \brief Accumulates and calculates standard deviation

       Typical usage:
       \code
       vector<double> v;
       utl::Accumulator::SampleStandardDeviation sigma;
       for_each(v.begin(), v.end(), sigma);
       cout << sigma.GetStandardDeviation() << endl;
       \endcode
       or in within some loop
       \code
       vector<double> v;
       utl::Accumulator::SampleStandardDeviation sigma;
       for (vector<double>::const_iterator it = v.begin(); it != v.end(); ++it) {
         ...
         sigma(*it);
         ...
       }
       \endcode
    */
    class SampleStandardDeviation : public SampleVariance {
    public:
      SampleStandardDeviation() : SampleVariance() { }
      void operator()(const double x) { SampleVariance::operator()(x); }
      double GetStandardDeviation() const { return std::sqrt(GetVariance()); }
    };


    //


    class CorrelationCoefficient {
    public:
      unsigned int GetN() const { return fX.GetN(); }
      double GetAverageX() const { return fX.GetAverage(); }
      double GetAverageY() const { return fY.GetAverage(GetN()); }
      double GetAverageXY() const { return fXY.GetAverage(GetN()); }
      double GetStandardDeviationX() const { return fX.GetStandardDeviation(); }
      double GetStandardDeviationY() const { return fY.GetStandardDeviation(GetN()); }
      double GetCorrelationCoefficient() const
      { return GetN() * (GetAverageXY() - GetAverageX()*GetAverageY()) /
          ((GetN() - 1) * GetStandardDeviationX() * GetStandardDeviationY()); }
      void Clear() { fX.Clear(); fY.Clear(); fXY.Clear(); }
      void operator()(const double x, const double y) { fX(x); fY(y); fXY(x*y); }
    private:
      SampleStandardDeviation fX;
      SampleStandardDeviationN fY;
      AverageN fXY;
    };


    //


    class Threshold {
    public:
      Threshold(const double threshold, const double baseline = 0,
                const int n = 1, const int multiplicity = 0) :
        fN(n),
        fThreshold(n, threshold),
        fBaseline(n, baseline),
        fMultiplicity(multiplicity ? multiplicity : n)
      { }

      /*Threshold(const double threshold,
                  const int n = 1, const int multiplicity = 0) :
        fN(n),
        fThreshold(n, threshold),
        fBaseline(n, 0),
        fMultiplicity(multiplicity ? multiplicity : n)
      { }*/

      /*template<class Vector>
      Threshold(const Vector& threshold, const int n = 1, const int multiplicity = 0) :
        fN(n),
        fThreshold(&threshold[0], &threshold[n]),
        fBaseline(n, 0),
        fMultiplicity(multiplicity ? multiplicity : n)
      { }*/

      template<class Vector>
      Threshold(const Vector& threshold, const Vector& baseline,
                const int n = 1, const int multiplicity = 0) :
        fN(n),
        fThreshold(&threshold[0], &threshold[n]),
        fBaseline(&baseline[0], &baseline[n]),
        fMultiplicity(multiplicity ? multiplicity : n)
      { }

      template<class Vector>
      bool
      operator()(const Vector* const v)
        const
      {
        int coincidences = 0;
        for (int i = 0; i < fN; ++i)
          coincidences += ((v[i] - fBaseline[i]) >= fThreshold[i]);
        return (coincidences >= fMultiplicity);
      }

    private:
      const int fN;
      const std::vector<double> fThreshold;
      const std::vector<double> fBaseline;
      const int fMultiplicity;
    };


    class ThresholdLatch {
    public:
      ThresholdLatch(const double threshold, const double baseline,
                     const int n = 1, const int multiplicity = 0) :
        fThreshold(threshold, baseline, n, multiplicity ? multiplicity : n),
        fLatch(false)
      { }

      /*ThresholdLatch(const double threshold,
                       const int n = 1, const int multiplicity = 0) :
        fThreshold(threshold, n, multiplicity ? multiplicity : n),
        fLatch(false)
      { }*/

      /*template<class Vector>
      ThresholdLatch(const Vector& threshold,
                     const int n = 1, const int multiplicity = 0) :
        fThreshold(threshold, n, multiplicity ? multiplicity : n),
        fLatch(false)
      { }*/

      template<class Vector>
      ThresholdLatch(const Vector& threshold, const Vector& baseline,
                     const int n = 1, const int multiplicity = 0) :
        fThreshold(threshold, baseline, n, multiplicity ? multiplicity : n),
        fLatch(false)
      { }

      template<typename T>
      bool
      operator()(const T* const v)
      {
        const bool trigger = fThreshold(v);
        fLatch(trigger);
        return trigger;
      }

      bool GetLatch() const { return fLatch.GetMax(); }

      void Clear() { fLatch.Clear(false); }

    private:
      Threshold fThreshold;
      Max<bool> fLatch;
    };


    class TimeOverThreshold {
    public:
      /*TimeOverThreshold(const double threshold,
                          const int nBins, const int window,
                          const int n = 1, const int multiplicity = 0) :
        fThreshold(threshold, n, multiplicity ? multiplicity : n),
        fNBins(nBins),
        fWindow(window)
      { Clear(); }*/

      TimeOverThreshold(const double threshold, const double baseline,
                        const int nBins, const int window,
                        const int n = 1, const int multiplicity = 0) :
        fThreshold(threshold, baseline, n, multiplicity ? multiplicity : n),
        fNBins(nBins),
        fWindow(window)
      { Clear(); }

      /*template<class Vector>
      TimeOverThreshold(const Vector& threshold,
                        const int nBins, const int window,
                        const int n = 1, const int multiplicity = 0) :
        fThreshold(threshold, n, multiplicity ? multiplicity : n),
        fNBins(nBins),
        fWindow(window)
      { Clear(); }*/

      template<class Vector>
      TimeOverThreshold(const Vector& threshold, const Vector& baseline,
                        const int nBins, const int window,
                        const int n = 1, const int multiplicity = 0) :
        fThreshold(threshold, baseline, n, multiplicity ? multiplicity : n),
        fNBins(nBins),
        fWindow(window)
      { Clear(); }

      template<typename T>
      bool
      operator()(const T* const v)
      {
        // fifo
        const bool trigger = fThreshold(v);
        fTriggersInWindow.push_front(trigger);
        fTriggerCount += trigger - fTriggersInWindow.back();
        fTriggersInWindow.pop_back();
        return fTriggerCount >= fNBins;
      }

      void
      Clear()
      {
        fTriggerCount = 0;
        fTriggersInWindow.clear();
        fTriggersInWindow.resize(fWindow, false);
      }

    private:
      Threshold fThreshold;
      int fNBins;
      int fWindow;
      int fTriggerCount;
      std::deque<bool> fTriggersInWindow;
    };


    class TimeOverThresholdDeconvoluted {
    public:
      /*TimeOverThresholdDeconvoluted(const double threshold, const double decayBins,
                                      const int nBins, const int window,
                                      const int n = 1, const int multiplicity = 0) :
        fN(n),
        fDeconvolutionFactor(std::exp(-1/decayBins)),
        fMultiplicity(multiplicity),
        fBaseline(n, 0),
        fPrevious(n, 0),
        fCurrent(n, 0),
        fTimeOverThreshold(n)
      {
        for (int i = 0; i < n; ++i)
          fTimeOverThreshold[i].reset(new TimeOverThreshold(threshold, nBins, window, 1, 1));
        Clear();
      }*/

      TimeOverThresholdDeconvoluted(const double threshold, const double baseline,
                                    const double decayBins,
                                    const int nBins, const int window,
                                    const int n = 1, const int multiplicity = 0) :
        fN(n),
        fDeconvolutionFactor(std::exp(-1/decayBins)),
        fMultiplicity(multiplicity),
        fBaseline(n, baseline),
        fPrevious(n, 0),
        fCurrent(n, 0),
        fTimeOverThreshold(n)
      {
        // ToT triggers have zero baseline because they are subtracted in this class
        for (int i = 0; i < n; ++i)
          fTimeOverThreshold[i].reset(new TimeOverThreshold(threshold, 0., nBins, window, 1, 1));
        Clear();
      }

      /*template<class Vector>
      TimeOverThresholdDeconvoluted(const Vector& threshold, const double decayBins,
                                    const int nBins, const int window,
                                    const int n = 1, const int multiplicity = 0) :
        fN(n),
        fDeconvolutionFactor(std::exp(-1/decayBins)),
        fMultiplicity(multiplicity),
        fBaseline(n, 0),
        fPrevious(n, 0),
        fCurrent(n, 0),
        fTimeOverThreshold(n)
      {
        for (int i = 0; i < n; ++i)
          fTimeOverThreshold[i].reset(new TimeOverThreshold(threshold, nBins, window, 1, 1));
        Clear();
      }*/

      template<class Vector>
      TimeOverThresholdDeconvoluted(const Vector& threshold, const Vector& baseline, const double decayBins,
                                    const int nBins, const int window,
                                    const int n = 1, const int multiplicity = 0) :
        fN(n),
        fDeconvolutionFactor(std::exp(-1/decayBins)),
        fMultiplicity(multiplicity),
        fBaseline(&baseline[0], &baseline[n]),
        fPrevious(n, 0),
        fCurrent(n, 0),
        fTimeOverThreshold(n)
      {
        // ToT triggers have zero baseline because they are subtracted in this class
        Vector zeroBaseline;
        for (unsigned int i = 0; i != fBaseline.size(); ++i)
          zeroBaseline[i] = 0;

        for (int i = 0; i < n; ++i)
          fTimeOverThreshold[i].reset(new TimeOverThreshold(threshold, zeroBaseline, nBins, window, 1, 1));
        Clear();
      }

      template<typename T>
      bool
      operator()(const T* const v)
      {
        for (int i = 0; i < fN; ++i) {
          fCurrent[i] = (v[i] - fBaseline[i] - fDeconvolutionFactor*fPrevious[i]) / (1 - fDeconvolutionFactor);
          fPrevious[i] = v[i] - fBaseline[i];
        }

        int coincidences = 0;
        for (int i = 0; i < fN; ++i)
          coincidences += (*fTimeOverThreshold[i])(&fCurrent[i]);

        return (coincidences >= fMultiplicity);
      }

      void
      Clear()
      {
        for (int i = 0; i < fN; ++i) {
          fPrevious[i] = 0;
          fCurrent[i] = 0;
        }
        for (unsigned int i = 0; i < fTimeOverThreshold.size(); ++i)
          fTimeOverThreshold[i]->Clear();
      }

    private:
      const int fN;
      const double fDeconvolutionFactor;
      const int fMultiplicity;

      const std::vector<double> fBaseline;

      std::vector<double> fPrevious;
      std::vector<double> fCurrent;

      std::vector<boost::shared_ptr<TimeOverThreshold> > fTimeOverThreshold;
    };


    class MultiplicityOfPositiveSteps {
    public:
      MultiplicityOfPositiveSteps(const int low, const int high, const int nVeto,
                                  const int nBin, const int window, const int n = 1,
                                  const int multiplicity = 0) :
        fN(n),
        fMultiplicity(multiplicity ? multiplicity : n),
        //fPosSteps(0),
        //fJump(false),
        fNVeto(nVeto),
        fPrevious(n, 0),
        fJumpCount(n, 0),
        fWindow(window),
        fOneTrace(fWindow, false),
        fJumpsInWindow(n, fOneTrace),
        fCumulSteps(n, 0),
        fNSatisfy(0),
        fLVeto(n,0),
        fVeto(n,0),
        fHigh(high),
        fLow(low),
        fMin(nBin)
      { Clear(); }

      /*Threshold(const double threshold,
                  const int n = 1, const int multiplicity = 0) :
        fN(n),
        fThreshold(n, threshold),
        fBaseline(n, 0),
        fMultiplicity(multiplicity ? multiplicity : n)
      { }*/

      /*template<class Vector>
      Threshold(const Vector& threshold, const int n = 1, const int multiplicity = 0) :
        fN(n),
        fThreshold(&threshold[0], &threshold[n]),
        fBaseline(n, 0),
        fMultiplicity(multiplicity ? multiplicity : n)
      { }*/

      template<class Vector>
      MultiplicityOfPositiveSteps(const int low, const int high, const int nVeto,
                                  const int nBin, const int window, const int n = 1,
                                  const int multiplicity = 0) :
        fN(n),
        fMultiplicity(multiplicity ? multiplicity : n),
        //fPosSteps(0),
        //fJump(false),
        fNVeto(nVeto),
        fPrevious(n, 0),
        fJumpCount(n, 0),
        fJumpCountSum(0),
        fWindow(window),
        fOneTrace(fWindow, false),
        fJumpsInWindow(n, fOneTrace),
        fCumulSteps(n, 0),
        fNSatisfy(0),
        fLVeto(n,0),
        fVeto(n,0),
        fHigh(high),
        fLow(low),
        fMin(nBin)
      { Clear(); }

      template<typename T>
      bool
      operator()(const T* const v)
        //const
      {
        for (int i = 0; i < fN; ++i) {
          if (fPrevious[i] == -1) {
            fPrevious[i] = v[i];
            return false;
          }
          int isJump = 0;
          const int step = v[i] - fPrevious[i];
          if (step > 0)
            fCumulSteps[i] += step;
          else {
            isJump += (fCumulSteps[i] >= fLow && fCumulSteps[i] <= fHigh && fVeto[i] <= 0);
            if (fNVeto >= 0 && fCumulSteps[i] >= 16) {
              fLVeto[i] = 2;
              if (fCumulSteps[i] >= 32)
                ++fLVeto[i];
              if (fCumulSteps[i] >= 64)
                ++fLVeto[i];
              if (fCumulSteps[i] >= 128)
                ++fLVeto[i];
              if (fCumulSteps[i] >= 256)
                ++fLVeto[i];
              if (fCumulSteps[i] >= 512)
                ++fLVeto[i];
              fVeto[i] = std::max(fVeto[i], fNVeto + fLVeto[i]);
            }
            fCumulSteps[i] = 0;
          }

          if (fVeto[i] > 0)
            --fVeto[i];

          fJumpsInWindow[i].push_front(isJump);
          fJumpCount[i] += isJump - int(fJumpsInWindow[i].back());
          fNSatisfy += (fJumpCount[i] > fMin);
          fJumpsInWindow[i].pop_back();
          fPrevious[i] = v[i];
        }
        return (fNSatisfy >= fMultiplicity);
      }

      void
      Clear()
      {
        for (int i = 0; i < fN; ++i) {
          fPrevious[i] = -1;
          fJumpCount[i] = 0;
          fLVeto[i] = 0;
          fVeto[i] = 0;
          fJumpsInWindow[i].clear();
          fJumpsInWindow[i].resize(fWindow, false);
        }
        fJumpsInWindow.resize(fN, fOneTrace);
      }

    private:
      const int fN;
      const int fMultiplicity;
      int fNVeto;
      std::vector<int> fPrevious;
      std::vector<int> fJumpCount;
      int fJumpCountSum;
      int fPosSteps;
      int fWindow;
      std::deque<bool> fOneTrace;
      std::vector<std::deque<bool> > fJumpsInWindow;
      std::vector<int> fCumulSteps;
      int fNSatisfy;
      std::vector<int> fLVeto;
      std::vector<int> fVeto;
      const int fHigh;
      const int fLow;
      const int fMin;
    };


  } // namespace Accumulator
} // namespace utl


#endif