MeanSqDispArray.h

#ifndef UTIL_MEAN_SQ_DISP_ARRAY_H
#define UTIL_MEAN_SQ_DISP_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.
*/

// Needed in header
#include <util/param/ParamComposite.h>
#include <util/containers/RingBuffer.h>
#include <util/containers/Array.h>

// Needed in implementation
#include <util/accumulators/MeanSqDispArray.h>
#include <util/accumulators/setToZero.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 Data>
   class MeanSqDispArray : public ParamComposite 
   {
   
   public:
      
      /* 
      * Constructor.
      */
      MeanSqDispArray();

      /* 
      * Destructor.
      */
      ~MeanSqDispArray();

      void readParameters(std::istream& in);
   
      void setParam(int ensembleCapacity, 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 setNEnsemble(int nEnsemble);
   
      void clear();
   
      void sample(const Array<Data>& values);
   
      void output(std::ostream& out);
  
      int bufferCapacity() 
      {  return bufferCapacity_; }
  
      int nEnsemble() 
      {  return nEnsemble_; }
  
      int nSample() 
      {  return nSample_; }
   
   private:

      DArray< RingBuffer<Data> >  buffers_; 
   
      DArray<double>  sqDiffSums_;
   
      DArray<int>  nValues_;
   
      int  ensembleCapacity_;
      
      int  bufferCapacity_;
      
      int  nEnsemble_;
   
      int  nSample_;

      void allocate();
   
      double sqDiff(const Data& data1, const Data& data2);
   
   };


   /*
   * Default constructor.
   */
   template <typename Data>
   MeanSqDispArray<Data>::MeanSqDispArray() 
    : buffers_(),
      sqDiffSums_(),
      nValues_(),
      ensembleCapacity_(0),
      bufferCapacity_(0),
      nEnsemble_(0),
      nSample_(0)
   {  setClassName("MeanSqDispArray"); }

   /*
   * Default destructor.
   */
   template <typename Data>
   MeanSqDispArray<Data>::~MeanSqDispArray() 
   {}
 
   /*
   * Read parameters from file.
   */
   template <typename Data>
   void MeanSqDispArray<Data>::readParameters(std::istream& in)
   {
      read<int>(in, "ensembleCapacity", ensembleCapacity_);
      read<int>(in, "bufferCapacity", bufferCapacity_);
      allocate();
      nEnsemble_ = ensembleCapacity_;
   }
   
   /*
   * Set parameters and initialize.
   */
   template <typename Data> 
   void 
   MeanSqDispArray<Data>::setParam(int ensembleCapacity, int bufferCapacity)
   {
      ensembleCapacity_ = ensembleCapacity;
      bufferCapacity_ = bufferCapacity;
      allocate();
      // Set number of sequences to maximum capacity as a default.
      nEnsemble_  = ensembleCapacity;
   }

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

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

   /*
   * Serialize this MeanSqDispArray.
   */
   template <typename Data>
   template <class Archive>
   void MeanSqDispArray<Data>::serialize(Archive& ar, const unsigned int version)
   {
      ar & ensembleCapacity_;
      ar & bufferCapacity_;
      ar & nEnsemble_;
      ar & nValues_;
      ar & nSample_;
      ar & buffers_; 
      ar & sqDiffSums_;
   }

   /*
   * Save internal state to archive.
   */
   template <typename Data>
   void MeanSqDispArray<Data>::save(Serializable::OArchive &ar)
   { ar & *this; }
   
   /* 
   * Set previously allocated to initial empty state.
   */
   template <typename Data>
   void MeanSqDispArray<Data>::clear()
   {   
      nSample_ = 0;

      if (bufferCapacity_ > 0) {

         int i;
         for (i = 0; i < bufferCapacity_; ++i) {
            setToZero(sqDiffSums_[i]);
            nValues_[i] = 0;
         }
   
         for (i = 0; i < ensembleCapacity_; ++i) {
            buffers_[i].clear();
         }

      }

   }
   
   /*
   * Allocate arrays and an array of CyclicBuffer objects (private method).
   */
   template <typename Data>
   void MeanSqDispArray<Data>::allocate()
   { 
      if (bufferCapacity_ > 0) { 
 
         // Allocate accumulator arrays
         sqDiffSums_.allocate(bufferCapacity_);
         nValues_.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>
   inline void MeanSqDispArray<Data>::sample(const Array<Data>& values)
   {
      int i, j;

      ++nSample_;

      for (i = 0; i < nEnsemble_; ++i) {
         buffers_[i].append(values[i]);
      }
      
      int bufferSize = buffers_[0].size();
      for (j = 0; j < bufferSize; ++j) {
         for (i = 0; i < nEnsemble_; ++i) {
            sqDiffSums_[j] += sqDiff(buffers_[i][j], values[i]);
         }
         ++nValues_[j];
      };
   }
  
   template <>
   inline double 
   MeanSqDispArray<int>::sqDiff(const int& data1, const int& data2) 
   {
      int diff; 
      diff = data1 - data2;
      return double(diff*diff);
   }

   template <>
   inline double 
   MeanSqDispArray<double>::sqDiff(const double& data1, const double& data2) 
   {
      double diff; 
      diff = data1 - data2;
      return diff*diff;
   }

   template <>
   inline double 
   MeanSqDispArray<Vector>::sqDiff(const Vector& data1, const Vector& data2) 
   {
      Vector diff;
      diff.subtract(data1, data2);
      return diff.square();
   }

   /*
   * Final output
   */
   template <typename Data>
   void MeanSqDispArray<Data>::output(std::ostream& out) 
   {
      double msd;
   
      // Calculate and output mean-squared difference
      int bufferSize = buffers_[0].size();
      for (int i = 0; i < bufferSize; ++i) {
         msd = sqDiffSums_[i]/double(nValues_[i]*nEnsemble_);
         out << Int(i) << Dbl(msd) << std::endl;
      }
      
   }
   
}
#endif