AutoCorrArray.h

#ifndef UTIL_AUTO_CORR_ARRAY_H
#define UTIL_AUTO_CORR_ARRAY_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/DArray.h>          // member template
#include <util/containers/RingBuffer.h>      // member template parameter

#include <util/accumulators/setToZero.h>
#include <util/accumulators/product.h>
#include <util/containers/Array.h>      
#include <util/space/Vector.h>
#include <util/format/Int.h>
#include <util/format/Dbl.h>

#include <complex>
using std::complex;

namespace Util
{

   template <typename> class Array;
   
   template <typename Data, typename Product>
   class AutoCorrArray : public ParamComposite
   {
   
   public:
   
      AutoCorrArray();

      ~AutoCorrArray();
  
      virtual void readParameters(std::istream& in);
   
      void setParam(int ensembleCapacity, int bufferCapacity);
   
      virtual void loadParameters(Serializable::IArchive &ar);
   
      virtual void save(Serializable::OArchive &ar);
   
      void setNEnsemble(int nEnsemble);
   
      void clear();
   
      void sample(const Array<Data>& values);
   
      void output(std::ostream& out);
  
      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);

      int bufferCapacity() const;
  
      int nEnsemble() const;
  
      int nSample() const;
  
      Data average() const;
  
      double corrTime() const; 
      
   private:
   
      DArray< RingBuffer<Data> > buffers_; 
   
      DArray<Product> corr_;
   
      DArray<int> nCorr_;
   
      Data sum_;
   
      int  ensembleCapacity_;
 
      int  bufferCapacity_;
      
      int  nEnsemble_;
   
      int  nSample_;

      void allocate();
   
   };

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

   /*
   * Destructor.
   */
   template <typename Data, typename Product>
   AutoCorrArray<Data, Product>::~AutoCorrArray() 
   {}
 
   /*
   * Read parameters from file.
   */
   template <typename Data, typename Product>
   void AutoCorrArray<Data, Product>::readParameters(std::istream& in)
   {
      read<int>(in, "ensembleCapacity", ensembleCapacity_);
      read<int>(in, "bufferCapacity",   bufferCapacity_);
      nEnsemble_ = ensembleCapacity_;
      allocate();
   }
   
   /*
   * Set parameters and initialize.
   */
   template <typename Data, typename Product>
   void AutoCorrArray<Data, Product>::setParam(int ensembleCapacity, int bufferCapacity)
   {
      ensembleCapacity_ = ensembleCapacity;
      bufferCapacity_ = bufferCapacity;
      allocate();
      nEnsemble_ = ensembleCapacity;
   }

   /*
   * Load internal state from archive.
   */
   template <typename Data, typename Product>
   void AutoCorrArray<Data, Product>::loadParameters(Serializable::IArchive &ar)
   {
      loadParameter<int>(ar, "ensembleCapacity", ensembleCapacity_);
      loadParameter<int>(ar, "bufferCapacity",   bufferCapacity_);
      ar & nEnsemble_;
      ar & buffers_; 
      ar & corr_;
      ar & nCorr_;
      ar & sum_;
      ar & nSample_;
   }

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

   /*
   * Set or reset nEnsemble.
   */
   template <typename Data, typename Product>
   void AutoCorrArray<Data, Product>::setNEnsemble(int nEnsemble)
   {
      if (ensembleCapacity_ == 0) {
         UTIL_THROW("No memory has been allocated: ensembleCapacity_ == 0"); 
      }
      if (nEnsemble > ensembleCapacity_) {
         UTIL_THROW("nEnsemble > ensembleCapacity_");      
      }
      nEnsemble_ = nEnsemble;
   }

   /* 
   * Set accumulator to initial empty state.
   */
   template <typename Data, typename Product>
   void AutoCorrArray<Data, Product>::clear()
   {   
      setToZero(sum_);
      nSample_ = 0;
      if (bufferCapacity_ > 0) {
         int i;
         for (i = 0; i < bufferCapacity_; ++i) {
            setToZero(corr_[i]);
            nCorr_[i] = 0;
         }
         for (i = 0; i < ensembleCapacity_; ++i) {
            buffers_[i].clear();
         }
      }
   }
   
   /*
   * Allocate arrays and CyclicBuffer (private method).
   */
   template <typename Data, typename Product>
   void AutoCorrArray<Data, Product>::allocate()
   { 
      if (bufferCapacity_ > 0) { 
         // Allocate autocorrelation accumulators
         corr_.allocate(bufferCapacity_);
         nCorr_.allocate(bufferCapacity_);
   
         // Allocate buffers
         buffers_.allocate(ensembleCapacity_);
         for (int i=0; i < ensembleCapacity_; ++i) {
            buffers_[i].allocate(bufferCapacity_);
         }
      }
      clear();
   }
   
   /*
   * Sample a single value from a time sequence.
   */
   template <typename Data, typename Product>
   void AutoCorrArray<Data, Product>::sample(const Array<Data>& values)
   {
      int i, j;
      ++nSample_;

      for (i = 0; i < nEnsemble_; ++i) {
         sum_ += values[i];
         buffers_[i].append(values[i]);
      }
      
      int bufferSize = buffers_[0].size();
      for (j=0; j < bufferSize; ++j) {
         for (i = 0; i < nEnsemble_; ++i) {
            corr_[j] += product(buffers_[i][j], values[i]);
         }
         ++nCorr_[j];
      };
   }
   
   /*
   * Return capacity of history buffer for each sequence.
   */
   template <typename Data, typename Product>
   int AutoCorrArray<Data, Product>::bufferCapacity() const
   {  return bufferCapacity_; }

   /*
   * Return number of independent sequences.
   */
   template <typename Data, typename Product>
   int AutoCorrArray<Data, Product>::nEnsemble() const
   {  return nEnsemble_; }

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

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

   /*
   * Final output.
   */
   template <typename Data, typename Product>
   void AutoCorrArray<Data, Product>::output(std::ostream& outFile) 
   {
      Product autocorr;
      // Product aveSq;
      // Data    ave = average();
   
      // Calculate and output autocorrelation
      // aveSq = product(ave, ave);
      int bufferSize = buffers_[0].size();
      for (int i = 0; i < bufferSize; ++i) {
         autocorr = corr_[i]/double(nCorr_[i]*nEnsemble_);
         //autocorr = autocorr - aveSq;
         //outFile << Int(i) << Dbl(autocorr) << Int(nCorr_[i]) << std::endl;
         outFile << Int(i) << Dbl(autocorr) << std::endl;
      }
   }
   
   /*
   *  Return correlation time in unit of sampling interval  
   */
   template <typename Data, typename Product>
   double AutoCorrArray<Data, Product>::corrTime() const
   {
      Data    ave;
      Product aveSq, variance, autocorr, sum;
      int     bufferSize = buffers_[0].size();
   
      // Calculate average of sampled values
      ave  = sum_;
      ave /= double(nSample_*nEnsemble_);
      aveSq = product(ave, ave);
      variance = corr_[0]/double(nCorr_[0]*nEnsemble_);
      variance = variance - aveSq;
   
      // Sum over autocorrelation function
      setToZero(sum);
      for (int i = 1; i < bufferSize/2; ++i) {
         autocorr = corr_[i]/double(nCorr_[i]*nEnsemble_);
         autocorr = autocorr - aveSq;
         sum += autocorr;
      }
      sum /= variance;
      return sum; 
   }
  
   /*
   * Serialize this AutoCorrArray.
   */
   template <typename Data, typename Product>
   template <class Archive>
   void AutoCorrArray<Data, Product>::serialize(Archive& ar, 
                                                const unsigned int version)
   {
      ar & ensembleCapacity_;
      ar & bufferCapacity_;
      ar & nEnsemble_;
      ar & buffers_; 
      ar & corr_;
      ar & nCorr_;
      ar & sum_;
      ar & nSample_;
   }

}
#endif