AutoCorr.h

#ifndef UTIL_AUTO_CORR_H
#define UTIL_AUTO_CORR_H

/*
* Util Package - C++ Utilities for Scientific Computation
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <util/param/ParamComposite.h>    // Base class
#include <util/containers/RingBuffer.h>   // member
#include <util/containers/DArray.h>       // member

// Needed for implementation
#include <util/accumulators/setToZero.h>
#include <util/accumulators/product.h>
#include <util/space/Vector.h>
#include <util/format/Int.h>
#include <util/format/Dbl.h>
#include <util/format/write.h>

#include <complex>

using std::complex;

namespace Util
{

   template <typename Data, typename Product>
   class AutoCorr : public ParamComposite
   {

   public:

      AutoCorr();

      ~AutoCorr();

      void clear();

      void readParameters(std::istream& in);

      void setParam(int bufferCapacity);

      virtual void loadParameters(Serializable::IArchive& ar);

      virtual void save(Serializable::OArchive& ar);

      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);

      void sample(Data value);

      void output(std::ostream& out);

      int bufferCapacity() const;

      int nSample() const;

      Data average() const;

      double corrTime() const;

      Product autoCorrelation(int t) const;

   private:

      // Ring buffer containing a sequence of stored Data values.
      RingBuffer<Data>  buffer_;

      // Array in which corr[j] = sum of values of <x(i-j), x(i)>
      DArray<Product>  corr_;

      // Array in which nCorr[i] = number of values added to corr[i]
      DArray<int>      nCorr_;

      // Sum of all previous values of x(t)
      Data        sum_;

      // Physical capacity (# of elements) of buffer, corr, and nCorr
      int         bufferCapacity_;

      // Total number of previous values of x(t)
      int         nSample_;

      void allocate();

      bool isValid();

   };

   /*
   * Default constructor.
   */
   template <typename Data, typename Product>
   AutoCorr<Data, Product>::AutoCorr()
    : buffer_(),
      corr_(),
      nCorr_(),
      bufferCapacity_(0),
      nSample_(0)
   {
      setClassName("AutoCorr");
      setToZero(sum_);
   }

   /*
   * Destructor
   */
   template <typename Data, typename Product>
   AutoCorr<Data, Product>::~AutoCorr()
   {}

   /*
   * Read buffer capacity and allocate all required memory.
   */
   template <typename Data, typename Product>
   void AutoCorr<Data, Product>::readParameters(std::istream& in)
   {
      read<int>(in, "capacity", bufferCapacity_);
      allocate();
   }

   /*
   * Set buffer capacity and initialize.
   */
   template <typename Data, typename Product>
   void AutoCorr<Data, Product>::setParam(int bufferCapacity)
   {
      bufferCapacity_ = bufferCapacity;
      allocate();
   }

   /*
   * Load state from an archive.
   */
   template <typename Data, typename Product>
   void AutoCorr<Data, Product>::loadParameters(Serializable::IArchive& ar)
   {
      loadParameter<int>(ar, "capacity", bufferCapacity_);
      ar & buffer_;
      ar & corr_;
      ar & nCorr_;
      ar & sum_;
      ar & nSample_;
      isValid();
   }

   /*
   * Save state to an archive.
   */
   template <typename Data, typename Product>
   void AutoCorr<Data, Product>::save(Serializable::OArchive& ar)
   {  ar & *this; }

   /*
   * Set previously allocated to initial empty state.
   */
   template <typename Data, typename Product>
   void AutoCorr<Data, Product>::clear()
   {
      setToZero(sum_);
      nSample_ = 0;

      if (bufferCapacity_ > 0) {
         for (int i=0; i < bufferCapacity_; ++i) {
            setToZero(corr_[i]);
            nCorr_[i] = 0;
         }
         buffer_.clear();
      }
   }

   /*
   * Allocate arrays and CyclicBuffer, and initialize.
   */
   template <typename Data, typename Product>
   void AutoCorr<Data, Product>::allocate()
   {
      if (bufferCapacity_ > 0) {
         corr_.allocate(bufferCapacity_);
         nCorr_.allocate(bufferCapacity_);
         buffer_.allocate(bufferCapacity_);
      }
      clear();
   }

   /*
   * Are capacities consistent?
   */
   template <typename Data, typename Product>
   bool AutoCorr<Data, Product>::isValid()
   {
      bool valid = true;
      if (bufferCapacity_ != corr_.capacity()) valid = false;
      if (bufferCapacity_ != nCorr_.capacity()) valid = false;
      if (bufferCapacity_ != buffer_.capacity()) valid = false;
      if (!valid) {
         UTIL_THROW("Invalid AutoCorr");
      }
      return valid;
   }

   /*
   * Sample a single value from a time sequence.
   */
   template <typename Data, typename Product>
   void AutoCorr<Data, Product>::sample(Data value)
   {
      ++nSample_;
      sum_ += value;
      buffer_.append(value);
      for (int i=0; i < buffer_.size(); ++i) {
         corr_[i] += product(buffer_[i], value);
         ++nCorr_[i];
      };
   }

   /*
   * Return capacity of history buffer.
   */
   template <typename Data, typename Product>
   int AutoCorr<Data, Product>::bufferCapacity() const
   { return bufferCapacity_; }

   /*
   * Return number of values sampled thus far.
   */
   template <typename Data, typename Product>
   int AutoCorr<Data, Product>::nSample() const
   { return nSample_; }

   /*
   * Return average of sampled values.
   */
   template <typename Data, typename Product>
   Data AutoCorr<Data, Product>::average() const
   {
      Data ave  = sum_;
      ave /= double(nSample_);
      return ave;
   }

   /*
   * Final output
   */
   template <typename Data, typename Product>
   void AutoCorr<Data, Product>::output(std::ostream& outFile)
   {
      Data    ave;
      Product autocorr;
      Product aveSq;

      // Calculate and output average of sampled values
      ave  = sum_;
      ave /= double(nSample_);
      aveSq = product(ave, ave);

      // Calculate and output autocorrelation
      for (int i = 0; i < buffer_.size(); ++i) {
         autocorr = corr_[i]/double(nCorr_[i]);
         autocorr = autocorr - aveSq;
         outFile << Int(i);
         write<Product>(outFile, autocorr);
         outFile << std::endl;
      }

   }

   /*
   *  Return correlation time in unit of sampling interval
   */
   template <typename Data, typename Product>
   double AutoCorr<Data, Product>::corrTime() const
   {
      Data    ave;
      Product aveSq, variance, autocorr, sum;
      int     size;

      size = buffer_.size();

      // Calculate average sampled values
      ave   = sum_/double(nSample_);
      ave  /= double(nSample_);
      aveSq = product(ave, ave);

      // Calculate variance sampled values
      variance = corr_[0]/double(nCorr_[0]);
      variance = variance - aveSq;

      setToZero(sum);
      for (int i = 1; i < size/2; ++i) {
         autocorr = corr_[i]/double(nCorr_[i]);
         autocorr = autocorr - aveSq;
         sum     += autocorr;
      }
      sum = sum/variance;

      return sum;
   }

   /*
   * Return autocorrelation at a given lag time index
   *
   * \param t the lag time index
   */
   template <typename Data, typename Product>
   Product AutoCorr<Data, Product>::autoCorrelation(int t) const
   {
      Data    ave;
      Product autocorr;
      Product aveSq;

      // Calculate average of sampled values
      ave  = sum_;
      ave /= double(nSample_);
      aveSq = product(ave, ave);

      // Calculate and return autocorrelation
      assert(t < buffer_.size());
      autocorr = corr_[t]/double(nCorr_[t]);
      autocorr = autocorr - aveSq;

      return autocorr;
   }

   /*
   * Serialize this AutoCorr.
   */
   template <typename Data, typename Product>
   template <class Archive>
   void AutoCorr<Data, Product>::serialize(Archive& ar,
                                           const unsigned int version)
   {
      ar & bufferCapacity_;
      ar & buffer_;
      ar & corr_;
      ar & nCorr_;
      ar & sum_;
      ar & nSample_;
      isValid();
   }

}
#endif