provide updated files for those parts needed by cthyb.

README.md

+This library is intended to provide tools for the error analysis of Monte Carlo data

analysis_common.h

+/***************************************************************************
+ *   Copyright (C) 2009-2016 by Florian Goth   *
+ *   fgoth@physik.uni-wuerzburg.de   *
+ *                                                                         *
+ *   All rights reserved.                                                  *
+ *                                                                         *
+ *   Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: *
+ *     * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. *
+ *     * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. *
+ *     * Neither the name of the <ORGANIZATION> nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. *
+ *                                                                         *
+ *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS   *
+ *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT     *
+ *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR *
+ *   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR *
+ *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
+ *   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,   *
+ *   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR    *
+ *   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
+ *   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING  *
+ *   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS    *
+ *   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.          *
+ ***************************************************************************/
+#ifndef ANALYSIS_COMMON_H
+#define ANALYSIS_COMMON_H
+
+namespace mc_analysis_common
+{
+
+#define GCC_VERSION (__GNUC__ * 10000 \
+                     + __GNUC_MINOR__ * 100 \
+                     + __GNUC_PATCHLEVEL__)
+
+#if(GCC_Version < 40300)
+#include <vector>
+#endif
+
+/**
+The following helper classes are for properly cleaning up memory in the case that the used container has a virtual destructor.
+Starting with the helper functions available gcc 4.3 we can preoperly deal with it.
+For older versions we just catch the common case of the std::vector.
+*/
+template <bool isvirtual, class T>
+struct Destroy
+{
+    static inline void destroy(const T&) {}
+};
+
+template<class T>
+struct Destroy<false, T>
+{
+    static inline void destroy(T& t)
+    {
+        t.clear();
+        t.reserve(0);
+    }
+};
+
+#if(GCC_VERSION < 40300)
+template <class T>
+struct isnotVector
+{
+    enum
+    {
+        Ret = true
+    };
+};
+
+template <class T>
+struct isnotVector<std::vector<T> >
+{
+    enum
+    {
+        Ret = false
+    };
+};
+#endif
+
+template <class T>
+struct CallClear
+{
+    static inline void destroy(T& t)
+    {
+#if (GCC_VERSION >= 40300)
+        Destroy<__has_virtual_destructor(T), T>::destroy(t);
+#else
+        Destroy<isnotVector<T>::Ret, T>::destroy(t);
+#endif
+    }
+};
+
+};
+
+#endif

analysis_functions.h

+# ifndef MC_ANALYSIS_ANALYSIS_FUNCTIONS_H
+# define MC_ANALYSIS_ANALYSIS_FUNCTIONS_H // include guard
+# include<vector>
+# include "helper_classes.h"
+# include <cassert>
+
+namespace mc_analysis
+{
+        /**
+        This function calculates the expectation value of the samples stored in the container cont
+        @param cont the container that stores the samples
+        @return the average value of the stored samples
+        */
+        template<typename T>
+		typename T::value_type mean(const T& cont)
+        {
+                 typename T::value_type retval(cont[0]);
+                 for (unsigned int k = 1; k < cont.size(); ++k)
+                        retval += cont[k];
+                 return retval/static_cast<typename ToScalar<T>::RetType >(cont.size());
+        }
+
+	template <typename Scalar> 
+		Scalar mean_signed(const std::vector< Scalar > & signdata__ , const std::vector< Scalar > & sign__)
+	{
+		Scalar average_sign(0.), average_data(0.);
+		average_data=mean(signdata__);
+		average_sign=mean(sign__);
+		return average_data/average_sign;		
+	}
+
+	template <typename T>
+		typename T::value_type standard_deviation(const T& cont)
+		{
+                  typedef typename T::value_type VType;
+			assert(cont.size() > 1);
+			VType mean_ = mean(cont);
+			VType sum (Cwise<VType>::product(cont[0] - mean_, cont[0] - mean_) );
+			for(unsigned int i = 1; i < cont.size(); i++) sum += Cwise<VType>::product(cont[i] - mean_, cont[i] - mean_);
+			return Cwise<VType>::sqrt( sum/static_cast<typename ToScalar<VType>::RetType> (cont.size() - 1.0) );
+		}
+
+
+	template <typename Scalar> errordata_signed<Scalar> jackknife_signed(const std::vector<Scalar> & signdata__, const std::vector<Scalar> & sign__)
+		/**
+		 *	\brief Calculates the jackknife error, mean and average sign of the data.
+		 *	\param signdata__ a vector containing the measured data already multiplied by the sign in 
+		 *			  the Monte Carlo simulation
+		 *	\param sign__	  a vector containing the signs. Note that the sizes of the two vectors have
+		 *			  to match!
+		 */
+	{
+		errordata_signed<Scalar> errd;
+		Scalar average_sign, average_data;
+		average_data=mean(signdata__);
+		average_sign=mean(sign__);
+		Scalar length= static_cast<Scalar>(signdata__.size());
+		std::vector<Scalar> jackdata(signdata__.size());
+		// Fill Jackknife bins
+		for(unsigned int i=0; i< jackdata.size(); i++)
+		{
+			jackdata[i] = (length*average_data - signdata__[i])/(length*average_sign - sign__[i]);	
+		}
+		Scalar jackmean(0.);
+		for(unsigned int i=0; i< signdata__.size();i++)
+		{
+			jackmean += jackdata[i];
+		}
+		jackmean *= 1.0/length;
+
+		Scalar mean, error, sign;
+		mean= average_data/average_sign;
+		for(unsigned int i=0; i< signdata__.size(); i++)
+		{
+			error+= (jackmean - jackdata[i])*(jackmean - jackdata[i]);
+		}
+		error *= (length -1.)/length;
+		error = std::sqrt(error);
+
+		errd.set_mean(mean);
+		errd.set_error(error);
+		errd.set_sign(average_sign);
+		return errd;
+	}
+
+
+
+	template <class DataCont, class SignCont> errordata<typename DataCont::value_type> jackknife_signed(const DataCont & signdata__, const SignCont & sign__)
+		/**
+		 *	\brief Calculates the jackknife error, mean and average sign of the data.
+		 *	\param signdata__ a vector containing the measured data already multiplied by the sign in 
+		 *			  the Monte Carlo simulation
+		 *	\param sign__	  a vector containing the signs. Note that the sizes of the two vectors have
+		 *			  to match!
+		 */
+	{
+		// DEBUG
+		//std::cout << "\n\n\n Performing Jackknife analysis" << std::endl;
+		// END DEBUG
+		typedef typename DataCont::value_type ValT;
+		typedef typename ToScalar<ValT>::RetType ScalT;
+		errordata<ValT> errd;
+		ValT average_data(mean(signdata__));
+		ScalT average_sign(mean(sign__));
+		ScalT length= static_cast<ScalT>(signdata__.size());
+		std::vector<ValT> jackdata(signdata__.size(), signdata__[0]);
+		// Fill Jackknife bins
+		for(unsigned int i=0; i< jackdata.size(); i++)
+		{
+			jackdata[i] = (length*average_data - signdata__[i])/(length*average_sign - sign__[i]);	
+		}
+		ValT jackmean(jackdata[0]);
+		for(unsigned int i=1; i< signdata__.size();i++)
+		{
+			jackmean += jackdata[i];
+		}
+		jackmean *= 1.0/length;
+
+		ValT mean(average_data/average_sign);
+		ValT error((jackmean - jackdata[0])*(jackmean - jackdata[0]));
+		for(unsigned int i=1; i< signdata__.size(); i++)
+		{
+			error+= (jackmean - jackdata[i])*(jackmean - jackdata[i]);
+		}
+		error *= (length -1.)/length;
+		error = std::sqrt(error);
+
+		errd.set_mean(mean);
+		errd.set_error(error);
+		//DEBUG
+		//std::cout << "\tResults:\n";
+		//std::cout << "Mean.size " << mean.size() <<std::endl;
+		//std::cout << "\n\n\n" <<std::endl;
+		//END DEBUG
+		return errd;
+	}
+
+
+
+
+	template <typename Scalar, typename Func> errordata_signed<typename Func::res_t> jack_eval_signed(const std::vector<Scalar> & signdata__, const std::vector<Scalar> & sign__, const Func & f) 
+	{
+		std::vector<typename Func::res_t> jackdata(signdata__.size());
+		errordata_signed<typename Func::res_t> errd;
+		Scalar average_sign, average_data;
+		average_data=mean(signdata__);
+		average_sign=mean(sign__);
+		Scalar length= static_cast<Scalar>(signdata__.size());
+
+		for(unsigned int i=0; i<signdata__.size(); i++)
+		{
+			Scalar x_J; 
+			x_J=(length*average_data - signdata__[i])/(length*average_sign - signdata__[i]);
+			jackdata[i] = f(x_J);	
+		}
+
+		typename Func::res_t jackmean(0.);
+		for(unsigned int i=0; i< signdata__.size();i++)
+		{
+			jackmean += jackdata[i];
+		}
+		jackmean *= 1.0/static_cast<typename Func::res_t>(signdata__.size());
+
+		typename Func::res_t fun_mean(0.), fun_error(0.);
+		fun_mean= f(average_data/average_sign);
+		for(unsigned int i=0; i< signdata__.size(); i++)
+		{
+			fun_error+= (jackmean - jackdata[i])*(jackmean - jackdata[i]);
+		}
+		fun_error *= static_cast<typename Func::res_t>(signdata__.size() -1)/static_cast<typename Func::res_t>(signdata__.size());
+		fun_error = std::sqrt(fun_error);
+
+		errd.set_mean(fun_mean);
+		errd.set_error(fun_error);
+		errd.set_sign(average_sign);
+		return errd;
+	
+	}
+
+
+	template <typename Scalar> Scalar autocorrelation_function(const std::vector<Scalar> data, const unsigned int displacement)
+		/**
+		 * 	\brief Calculate the normalized Autocorrelation function A(displacement). It is defined by:
+		 *
+		 * 	\f[A(disp)= \frac{( <data[i] data[i+disp]> - <data[i]>^2 )}{( <data[i]^2> - <data[i]>^2 )} \f]
+		 * 	\param data	vector containing the data. The size of the vector should be much larger 
+		 * 			than the displacement.
+		 * 	\param displacement
+		 */
+	{
+		Scalar sum_disp(0.), sum_mean(0.), sum_sq(0.);	
+		for(unsigned int i=0; i< data.size()-displacement; i++)
+		{
+			sum_disp += data[i]*data[i+displacement];
+			sum_mean += data[i];
+			sum_sq += data[i]*data[i];	
+		}	
+		Scalar length=static_cast<Scalar>(data.size() - displacement);
+		return (sum_disp - 1.0/(length)*sum_mean*sum_mean )/( sum_sq - 1.0/(length)* sum_mean*sum_mean);
+	}
+
+        /**
+        A stream operator to output the elements of a valarray
+        @param out the stream in which to enter the elements
+        @param rhs the valarray to output
+        @return the modified stream
+        */
+        template <class T>
+        inline std::ostream& operator<<(std::ostream& out, const std::valarray<T>& rhs )
+        {
+               for (unsigned int i = 0; i < rhs.size(); ++i)
+                         out<<rhs[i]<<std::endl;
+               return out;
+        }
+}
+
+
+# endif // include guard

array_observable.h

+# ifndef MONTE_CARLO_TOOLS_ARRAY_OBSERVABLE_H
+# define MONTE_CARLO_TOOLS_ARRAY_OBSERVABLE_H
+# include "errordata.h"
+# include "analysis_functions.h"
+# include "analysis_common.h"
+namespace mc_analysis
+{
+
+	template <typename Scalar, class Cont = std::vector<std::valarray<Scalar> > >  class array_observable : public Cont
+	{
+		private:
+			std::valarray< std::valarray<Scalar> > data;///< the binned data
+			unsigned int num_bins;///< the number of bins
+			unsigned int binsize;///< how many measurements make up a in
+                        using Cont::size;//make the size function of the base class inaccesible to the outside world
+                        using Cont::operator[];
+
+
+
+		public:
+			inline array_observable(void);
+                        inline virtual ~array_observable();
+			inline std::size_t get_num_data_points(void){return Cont::size();};
+			inline std::valarray<Scalar> mean(void);
+			inline std::valarray<Scalar> sd(void); // This should be implemented by a general function, as mean
+	//		std::valarray<Scalar> standard_error(void);
+			inline errordata<std::valarray<Scalar> > jackknife(void);
+                        /**
+                        You call this function to rebin the measurements to the supplied binsize
+                        @param binsize__ how many measurements make up a bin
+                        */
+			inline void rebin_data(unsigned int binsize__);
+	//		void fractional_rebinning(unsigned int num_bins__);
+			inline void print_data(void);
+
+			template <typename Func> inline errordata<typename Func::res_t> jack_eval(Func & f);
+                        using Cont::value_type;
+                        /**
+                        these access operators are handy for the jackknife analysis and provide access to the binned data
+                        @param k the index of the Element
+                        @return a reference to the element
+                        */
+                        inline std::valarray<Scalar>& operator[] (std::size_t k) {return data[k];}
+                        inline const std::valarray<Scalar>& operator[] (std::size_t k) const {return data[k];}
+                        inline std::size_t size() const {return data.size();}
+			inline void clear_tmp(void);
+			/**
+			 * \brief delete all data
+			 */
+			inline void reset(void){clear_tmp(); this->clear();}
+	};
+
+	template <typename Scalar, class Cont> 
+        array_observable<Scalar, Cont>::~array_observable(void)
+	{
+            mc_analysis_common::CallClear<Cont>::destroy(static_cast<Cont&>(*this));
+	}
+
+	template <typename Scalar, class Cont> array_observable<Scalar, Cont>::array_observable(void)
+	{
+		num_bins=0;
+		binsize=0;
+		data.resize(0);	
+	}
+
+
+	template <typename Scalar, class Cont> std::valarray<Scalar> array_observable<Scalar, Cont>::mean(void)
+	{
+		return mc_analysis::mean(data);	
+	}
+
+	template <typename Scalar, class Cont> std::valarray<Scalar> array_observable<Scalar,Cont>::sd(void)
+	{
+/*		std::valarray<Scalar> sum(data[0].size());
+		sum=0.;
+		std::valarray<Scalar> mean(data[0].size());
+		mean=this->mean();
+		for(unsigned int i=0; i<num_bins; i++) sum += (data[i]-mean)*(data[i]-mean);
+		return std::sqrt( sum/static_cast<Scalar> (num_bins -1.) );
+		*/
+		return standard_deviation(data);
+	}
+
+
+
+	template <typename Scalar, class Cont> errordata<std::valarray<Scalar> > array_observable<Scalar, Cont>::jackknife(void)
+		/**
+		 *	\brief this function performs a simple jackknife analysis for the data
+		 *	it calculates the jackknife mean and the standard error.
+		 */
+	{
+		// For this special case, the Jackknife mean is the same as the arithmetic average of the data
+		errordata<std::valarray<Scalar> > errd;
+		errd.set_mean(this->mean());
+		// The Jackknife error for this case also reduces to the standard error:
+		errd.set_error(this->sd()/std::sqrt( static_cast<Scalar>(num_bins) ));
+		return errd;		
+	}
+
+
+
+
+	template <typename Scalar, class Cont> void array_observable<Scalar, Cont>::rebin_data(unsigned int binsize__)
+		/**
+		 *	\brief This function has to be called before using the data analysis functions. Rebinning is mandatory
+		 *	as it fills the data vector. The original data remains intact in the class
+		 *	\param binsize__ Specifies, how many data points should be averaged to give a new bin
+		 */
+	{
+		this->binsize=binsize__;
+		this->num_bins= Cont::size()/binsize__;//how many functions do we have
+		data.resize(num_bins, Cont::operator[](0));//make place for num_bins functions
+		for(unsigned int i=0;i<this->num_bins;i++)
+		{
+			std::valarray<Scalar> sum(Cont::operator[]( i * this->binsize  ));
+			for(unsigned int j = 1; j < (this->binsize); j++)
+			{
+				sum += Cont::operator[]( i*this->binsize + j );
+			}
+			data[i] = sum;
+			data[i] *= 1./static_cast<double>(this->binsize); 
+		}
+	}
+
+
+	template <typename Scalar, class Cont> void  array_observable<Scalar, Cont>::print_data(void)
+	{
+		for(unsigned int i=0; i< data.size() ; i++)
+		{
+			for(unsigned int j=0; j< data[i].size(); j++)
+			{
+				std::cout << data[i][j] << "\t";			
+			}
+			std::cout << std::endl;
+		}
+		std::cout << std::endl;
+	}
+
+
+
+	template <typename Scalar, class Cont> template <typename Func> errordata<typename Func::res_t> array_observable<Scalar, Cont>::jack_eval(Func & f)
+		/**
+		 * \brief This function calculates the mean and error of the data after applying a function
+		 * specified in a Function object.
+		 * 
+		 * \param F is a functor defined by the class Functor. You need to overload the operator().
+		 * and to specify the typedef res_t for the result/return type.
+		 *
+		 * \warning Rebinning before calling this function is absolutely required!
+		 */
+	{
+		// Generate Jackknife bins
+		Scalar length = static_cast<Scalar>(this->num_bins); 
+		std::valarray<Scalar> mean( this->mean() );
+		typename Func::res_t fun_mean=f(mean);
+
+		std::valarray<typename Func::res_t> jackdata(fun_mean,this->num_bins); // Init the jackknife bins with f(mean), to make them large enough
+
+
+		for(unsigned int i=0; i < this->num_bins; i++)
+		{
+			std::valarray<Scalar> x_J(mean.size());
+			x_J = (length * mean - data[i]);
+			x_J *= 1./(length - 1.0);
+			jackdata[i] = f(x_J);
+		}
+
+		typename Func::res_t jackmean(jackdata[0]);
+		for(unsigned int i=1; i< this->num_bins;i++)
+		{
+			jackmean += jackdata[i];
+		}
+		jackmean *= 1.0/(length);
+		
+
+		typename Func::res_t fun_error= (jackmean - jackdata[0])*(jackmean - jackdata[0]);
+
+		for(unsigned int i=1; i< this->num_bins; i++)
+		{
+			fun_error+= (jackmean - jackdata[i])*(jackmean - jackdata[i]);
+		}
+		fun_error *= (length-1.)/length;
+		fun_error = std::sqrt(fun_error);
+		
+		errordata<typename Func::res_t> errd;
+		errd.set_mean(fun_mean);
+		errd.set_error(fun_error);
+		return errd;
+	}
+
+	/**
+	 *	\brief Strip all temporary data from the array observable class
+	 *	
+	 */
+	template <typename Scalar, class Cont> void array_observable<Scalar, Cont>::clear_tmp(void)
+	{
+	    data.resize(0);
+            num_bins = 0;
+            binsize = 0;
+            return;
+	}
+}
+# endif

errordata.h

+# ifndef MONTE_CARLO_TOOLS_ERRORDATA_H // include guard
+# define MONTE_CARLO_TOOLS_ERRORDATA_H
+#include <valarray>
+#include <iostream>
+namespace mc_analysis
+{
+
+template <typename R> class errordata
+{
+private:
+    R mean;
+    R error;
+    R bias;///< the bias of the jackknife procedure. This is e.g useful for estimating the quality of the binsize. WARNING: not all functions of the mc_analysis class set this
+public:
+    /**
+    A constructor for the error data of a scalar observable
+    @param m this sets the estimation for the mean value of observable
+    @param e this set the estimation for the standard error of the observable
+    @param b this set the estimate for the bias
+    */
+    errordata(const R& m, const R& e, const R& b) : mean(m), error(e), bias(b) {}
+    /**
+    Two argument version of the constructor. The bias is set to zero.
+    @param m this sets the estimation for the mean value of observable
+    @param e this set the estimation for the standard error of the observable
+    */
+    errordata(const R& m, const R& e) : mean(m), error(e), bias(0) {}
+    void set_mean(R mean__)
+    {
+        mean=mean__;
+    };
+    void set_error(R error__)
+    {
+        error=error__;
+    };
+    R get_mean(void) const
+    {
+        return mean;
+    }
+    R get_error(void) const
+    {
+        return error;
+    }
+    /**
+    An accessor to the bias of the observable
+    @return the bias of the observable
+    */
+    R get_bias() const
+    {
+        return bias;
+    }
+    typedef R value_type;
+//    size_t size() const {return }
+};
+
+template <typename R>
+inline std::ostream& operator<<(std::ostream& out, const errordata<R>& arg )
+{
+    out<<arg.get_mean()<<" +- "<<arg.get_error()<<std::endl;
+    out<<"bias: "<<arg.get_bias();
+    return out;
+}
+
+// Specialization for valarray:
+template <typename Scalar> class errordata<std::valarray<Scalar> >
+{
+public:
+    typedef std::valarray<Scalar> CovType;
+    void setCov(std::valarray<Scalar>& rhs) {
+      cov.resize(rhs.size());
+      cov = rhs;      
+    }
+    const CovType& getCov() const {return cov;}
+    /**
+    A constructor for the error data of a scalar observable
+    @param m this sets the estimation for the mean value of observable
+    @param e this set the estimation for the standard error of the observable
+    @param b this set the estimate for the bias
+    */
+    errordata(const std::valarray<Scalar>& m, const std::valarray<Scalar>& e, const std::valarray<Scalar>& b) : mean(m), error(e), bias(b) {}
+    void set_mean(std::valarray<Scalar> mean_ )
+    {
+        mean.resize( mean_.size() );
+        mean=mean_;
+    }
+    void set_error(std::valarray<Scalar> error_ )
+    {
+        error.resize( error_.size() );
+        error=error_;
+    }
+    std::valarray<Scalar> get_mean(void) const
+    {
+        return mean;
+    }