Commit 4772ea07 authored by Florian Goth's avatar Florian Goth
Browse files

add code for calculation of covariance matrix.

parent dad70634
......@@ -69,6 +69,44 @@ private:
uint len;
};
template <typename T>
class ReVecFunction
{
public:
inline std::valarray<T> operator()(std::valarray<std::complex<T> >& dat)
{
auto denom = dat[len-1].real()*dat[len-1].real() + dat[len-1].imag()*dat[len-1].imag();
std::valarray<T> retval(len);
for(uint k = 0; k < (len - 1); ++k)
{
retval[k] = (dat[k].real()*dat[len-1].real() + dat[k].imag()*dat[len-1].imag())/denom;
}
return retval;
}
typedef std::valarray<T> res_t;
ReVecFunction(uint l) : len(l) {}
private:
uint len;
};
template <typename T>
class ImVecFunction
{
public:
inline std::valarray<T> operator() (std::valarray<std::complex<T> >& dat)
{
auto denom = dat[len-1].real()*dat[len-1].real() + dat[len-1].imag()*dat[len-1].imag();
std::valarray<T> retval(len);
for(uint k = 0; k < (len - 1); ++k)
retval[k] = (dat[k].imag() * dat[len-1].real() - dat[len-1].imag()*dat[k].real())/denom;
return retval;
}
typedef std::valarray<T> res_t;
ImVecFunction(uint l) : len(l) {}
private:
uint len;
};
template <typename T>
class PlainSign
{
......@@ -83,6 +121,23 @@ private:
uint len;
};
template <typename T>
class PlainVecSign
{
public:
inline T operator() (T* dat)
{
std::valarray<T> retval(len);
for(uint k = 0; k < len-1; ++k)
retval[k] = dat[k]/dat[len-1];
return retval;
}
typedef std::valarray<T> res_t;
PlainVecSign(uint l) : len(l) {}
private:
uint len;
};
template <typename T>
class PlainSign<std::complex<T> >
{//this is supposed to be used if we mix a complex observable and a real sign
......@@ -97,6 +152,23 @@ private:
uint len;
};
template <typename T>
class PlainVecSign<std::complex<T> >
{//this is supposed to be used if we mix a complex observable and a real sign
public:
inline std::valarray<std::complex<T> > operator() (std::complex<T>* dat)
{
std::valarray<std::complex<T> > retval(len);
for(uint k = 0; k < len-1; ++k)
retval[k] = dat[k]/ dat[len-1];
return retval;
}
typedef std::valarray<std::complex<T> > res_t;
PlainVecSign< std::complex< T > >(uint l) : len(l) {}
private:
uint len;
};
template <typename FPType>
inline mc_analysis::errordata<std::complex<FPType> > handlesign(FileVector<std::complex<FPType> >& a, FileVector<std::complex<FPType> >& b)
{
......@@ -126,8 +198,8 @@ inline mc_analysis::errordata<std::complex<FPType> > handlesign(FileVector<std::
cont[0][k] = a[k];
cont[1][k] = b[k];
}
mc_analysis::errordata<typename ReFunction<FPType>::res_t> ed = mc_analysis::jackknife(ReFunction<FPType>(2), cont, 2);
mc_analysis::errordata<typename ReFunction<FPType>::res_t> edim = mc_analysis::jackknife(ImFunction<FPType>(2), cont, 2);
mc_analysis::errordata<typename ReFunction<FPType>::res_t> ed = mc_analysis::jackknife(ReFunction<FPType>(2), (&cont[0]), 2);
mc_analysis::errordata<typename ImFunction<FPType>::res_t> edim = mc_analysis::jackknife(ImFunction<FPType>(2), (&cont[0]), 2);
return mc_analysis::errordata<std::complex<FPType> >(std::complex<FPType>(ed.get_mean(), edim.get_mean()), std::complex<FPType>(ed.get_error(), edim.get_error()), std::complex<FPType>(ed.get_bias(), edim.get_bias()));
}
......@@ -187,13 +259,76 @@ inline mc_analysis::errordata<FPType> handlesign(FileVector<FPType>& a, std::val
template <typename FPType>
mc_analysis::errordata<std::complex<FPType> > handlesign(const std::valarray<std::complex<FPType> >& a, const std::valarray<std::complex<FPType> >& b)
{
typedef std::complex<FPType> SignType;
const std::valarray<SignType> *const cont[2] = {&a, &b};
const std::valarray<std::complex<FPType> > *const cont[2] = {&a, &b};
mc_analysis::errordata<typename ReFunction<FPType>::res_t> ed = mc_analysis::jackknife (ReFunction<FPType>(2), cont, 2);
mc_analysis::errordata<typename ImFunction<FPType>::res_t> edim = mc_analysis::jackknife(ImFunction<FPType>(2), cont, 2);
return mc_analysis::errordata<std::complex<FPType> >(std::complex<FPType>(ed.get_mean(), edim.get_mean()), std::complex<FPType>(ed.get_error(), edim.get_error()), std::complex<FPType>(ed.get_bias(), edim.get_bias()));
}
template <typename FPType>
mc_analysis::errordata<std::valarray<std::complex<FPType> > > handlesign(const std::valarray<std::valarray<std::complex<FPType> > >& a, const std::valarray<FPType>& b)
{
std::valarray<std::complex<FPType> > sign(b.size());
for(uint k = 0; k < b.size(); ++k)
sign[k] = b[k];
const std::valarray<std::complex<FPType> > ** cont = new const std::valarray<std::complex<FPType> > *[a.size() + 1];
for(uint k = 0; k < a.size(); ++k)
cont[k] = &(a[k]);
cont[a.size()] = &sign;
mc_analysis::errordata<typename PlainVecSign<std::complex<FPType> >::res_t> ed = mc_analysis::jackknife(PlainVecSign<std::complex<FPType> >(a.size() + 1), cont, a.size() + 1);
return ed;
}
template <typename FPType>
mc_analysis::errordata<std::valarray<std::complex<FPType> > > handlesign(const std::valarray<std::valarray<std::complex<FPType> > >& a, bool covariance)
{
uint funnr = a[0].size();
mc_analysis::errordata<typename ReVecFunction<FPType>::res_t> ed = mc_analysis::vecjackknife(ReVecFunction<FPType>(funnr), a, funnr, covariance);
mc_analysis::errordata<typename ImVecFunction<FPType>::res_t> edim = mc_analysis::vecjackknife(ImVecFunction<FPType >(funnr), a, funnr, covariance);
std::valarray<std::complex<FPType> > retval(funnr-1);
std::valarray<std::complex<FPType> > retvalerr(funnr-1);
std::valarray<std::complex<FPType> > retvalbias(funnr-1);
std::valarray<std::complex<FPType> > retvalcov((funnr - 1) * (funnr - 1));
for(uint y = 0; y < (funnr-1); ++y)
{
retval[y] = std::complex<FPType>(ed.get_mean()[y], edim.get_mean()[y]);
retvalerr[y] = std::complex<FPType>(ed.get_error()[y], edim.get_error()[y]);
retvalbias[y] = std::complex<FPType>(ed.get_bias()[y], edim.get_bias()[y]);
if(covariance)
{
for(uint j = 0; j < (funnr - 1); ++j)
{
retvalcov[y*(funnr - 1) + j] = std::complex<FPType>(ed.getCov()[y*funnr + j], edim.getCov()[y*funnr + j]);
}
}
}
mc_analysis::errordata<std::valarray<std::complex<FPType> > > finalretval(retval, retvalerr, retvalbias);
if(covariance)
finalretval.setCov(retvalcov);
return finalretval;
}
template <typename FPType>
mc_analysis::errordata<std::valarray<std::complex<FPType> > > handlesign(const FileVector<std::valarray<std::complex<FPType> > >& a, const std::valarray<std::complex<FPType> >& b)
{
const std::valarray<std::complex<FPType> > *const* cont = new const std::valarray<std::complex<FPType> > *const[a.size() + 1];
for(uint k = 0; k < a.size(); ++k)
cont[k] = &(a[k]);
cont[a.size()] = &b;
mc_analysis::errordata<typename ReVecFunction<FPType>::res_t> ed = mc_analysis::jackknife (ReFunction<FPType>(a.size() + 1), cont, a.size() + 1);
mc_analysis::errordata<typename ImVecFunction<FPType>::res_t> edim = mc_analysis::jackknife(ImFunction<FPType>(a.size() + 1), cont, a.size() + 1);
std::valarray<std::complex<FPType> > retval(a.size());
std::valarray<std::complex<FPType> > retvalerr(a.size());
std::valarray<std::complex<FPType> > retvalbias(a.size());
for(uint k = 0; k < a.size(); ++k)
{
retval[k] = std::complex<FPType>(ed.get_mean()[k], edim.get_mean()[k]);
retvalerr[k] = std::complex<FPType>(ed.get_error()[k], edim.get_error()[k]);
retvalbias[k] = std::complex<FPType>(ed.get_bias()[k], edim.get_bias()[k]);
}
return mc_analysis::errordata<std::complex<FPType> >(retval, retvalerr, retvalbias);
}
template<class ObsType, typename SignType>
class Error_Analysis
{
......@@ -322,10 +457,11 @@ struct Point
@param delta_s the resolution of the function
*/
template <class Cont, class SignType, typename FPType>
static inline void writeVectorFunctiontoFile(Cont& container, unsigned int len, const valarray<SignType>& avsign, std::string& nameRe, std::string& nameIm, unsigned int slices, FPType delta_s)
static inline void writeVectorFunctiontoFile(Cont& container, unsigned int len, const valarray<SignType>& avsign, std::string& nameRe, std::string& nameIm, unsigned int slices, FPType delta_s, bool covariance)
{
std::ofstream fileRe(nameRe.c_str());
std::ofstream fileIm(nameIm.c_str());
std::ofstream covfile((nameRe+"Cov").c_str());
//the template parameter Cont is usually a filevector
typedef typename Cont::value_type::value_type::value_type ScalarType;
std::ofstream xdep[5];
......@@ -337,61 +473,55 @@ static inline void writeVectorFunctiontoFile(Cont& container, unsigned int len,
container.sync();
std::vector<Point<FPType> > points;
points.reserve(4*slices);
unsigned int nrbins = container.size();//get nr. of bins
uint availmem = sysconf (_SC_AVPHYS_PAGES) * sysconf(_SC_PAGESIZE);//some approximation to the amount of available memory
uint memoryperpoint = sizeof(ScalarType)*nrbins;
int possibleparallelism = std::min(availmem / memoryperpoint, slices);
std::valarray<ScalarType>* temp = NULL;
bool memcritical = (possibleparallelism == 0);
std::cout<<"Possible parallelism: "<<possibleparallelism<<std::endl;
if(!memcritical)
unsigned int nroffunctions = container.size();//get nr. of bins
std::valarray<std::valarray<ScalarType> > temp(nroffunctions + 1);//allocate a valarray for one function
for(uint j = 0; j < (nroffunctions+1); ++j)
temp[j].resize(slices+1);
for(uint k = 0; k < len; ++k) // for every point
{
try
for (uint j = 0; j < nroffunctions; ++j) //get every bin
{
temp = new std::valarray<ScalarType>[possibleparallelism];
for(int j = 0; j < possibleparallelism; ++j)
temp[j].resize(nrbins);
auto ret = container(j, k);
for(uint i = 0; i < ret.size(); ++i)
temp[j][i] = ret[i];
temp[j][ret.size()] = avsign[j];
}
catch(std::bad_alloc& e)
mc_analysis::errordata<std::valarray<ScalarType> > edp(handlesign(temp, covariance));//do sign - analysis
for(uint l = 0; l < slices; ++l)
{
delete [] temp;
memcritical = true;
}
}
if(memcritical)
{//memory is critical. let's try it anyway and see if the kernel frees memory.
possibleparallelism = 1;
temp = new std::valarray<ScalarType>(nrbins);
}
// std::valarray<ScalarType> temp(nrbins);//allocate a valarray for one function
for (unsigned int k = 0; k < len; ++k)//for every index of the function that we stored
{
for (unsigned int j = 0; j < slices; j += possibleparallelism)//for all functionpoints
{
int upperlimit = possibleparallelism;
if(j + possibleparallelism >= slices)
upperlimit = slices - j;//determine the remaining indices
for (unsigned int l = 0; l < nrbins; ++l)//copy function points
{
auto tempdatafromfile = container(l, k);
for(int j1 = 0; j1 < upperlimit; ++j1)
temp[j1][l] = tempdatafromfile[j+j1];
// temp[l] = container(l, k)[j];
}
for(int j1 = 0; j1 < upperlimit; ++j1)
{
mc_analysis::errordata<ScalarType> edp(handlesign(temp[j1], avsign));//do sign - analysis
//output data
points.push_back(Point<FPType>((j+j1) * delta_s, k, real(edp.get_mean()), real(edp.get_error())));
fileRe<<(j+j1) * delta_s<<" "<<real(edp.get_mean())<<" "<<real(edp.get_error())<<std::endl;
fileIm<<(j+j1) * delta_s<<" "<<imag(edp.get_mean())<<" "<<imag(edp.get_error())<<std::endl;
}
}
points.push_back(Point<FPType>(l * delta_s, k, real(edp.get_mean()[l]), real(edp.get_error()[l])));
fileRe<<l * delta_s<<" "<<real(edp.get_mean()[l])<<" "<<real(edp.get_error()[l])<<std::endl;
fileIm<<l * delta_s<<" "<<imag(edp.get_mean()[l])<<" "<<imag(edp.get_error()[l])<<std::endl;
}
//output xmgrace separators
fileRe<<"&"<<std::endl;
fileIm<<"&"<<std::endl;
if(covariance)
{
for(uint k = 0; k < slices; ++k)
{
for(uint j = 0; j < slices; ++j)
covfile<<real(edp.getCov()[k*slices + j])<<" ";
covfile<<std::endl;
}
}
delete [] temp;
}
// for (unsigned int k = 0; k < len; ++k)//for every function that we stored
// {
// for (unsigned int j = 0; j < slices; ++j)//for all functionpoints
// {
// for (unsigned int l = 0; l < nroffunctions; ++l)//copy function points
// temp[l] = container(l,k)[j];
// mc_analysis::errordata<FPType> edp(handlesign(temp, avsign));//do sign - analysis
// //output data
// points.push_back(Point<FPType>(j * delta_s, k, edp.get_mean(), edp.get_error()));
// file<<j * delta_s<<" "<<edp.get_mean()<<" "<<edp.get_error()<<std::endl;
// }
// //output xmgrace separators
// file<<"&"<<std::endl;
// }
//now let's generate the transformed data
sort(points.begin(), points.end());
#pragma omp parallel for schedule(dynamic)
......@@ -410,6 +540,102 @@ static inline void writeVectorFunctiontoFile(Cont& container, unsigned int len,
return;
}
// /**
// @param fileRe the file to which we write the realpart
// @param fileIm the file to which we write the imaginarypart
// @param len the length of the Vector
// @param slices the number of points the function has
// @param delta_s the resolution of the function
// */
// template <class Cont, class SignType, typename FPType>
// static inline void writeVectorFunctiontoFile(Cont& container, unsigned int len, const valarray<SignType>& avsign, std::string& nameRe, std::string& nameIm, unsigned int slices, FPType delta_s)
// {
// std::ofstream fileRe(nameRe.c_str());
// std::ofstream fileIm(nameIm.c_str());
// //the template parameter Cont is usually a filevector
// typedef typename Cont::value_type::value_type::value_type ScalarType;
// std::ofstream xdep[5];
// for (unsigned int i = 0; i < 5; ++i)
// {
// unsigned int dt = 2 << i;
// xdep[i].open((nameRe + "_" + toString(dt)).c_str());
// }
// container.sync();
// std::vector<Point<FPType> > points;
// points.reserve(4*slices);
// unsigned int nrbins = container.size();//get nr. of bins
// uint availmem = sysconf (_SC_AVPHYS_PAGES) * sysconf(_SC_PAGESIZE);//some approximation to the amount of available memory
// uint memoryperpoint = sizeof(ScalarType)*nrbins;
// int possibleparallelism = std::min(availmem / memoryperpoint, slices);
// std::valarray<ScalarType>* temp = NULL;
// bool memcritical = (possibleparallelism == 0);
// std::cout<<"Possible parallelism: "<<possibleparallelism<<std::endl;
// if(!memcritical)
// {
// try
// {
// temp = new std::valarray<ScalarType>[possibleparallelism];
// for(int j = 0; j < possibleparallelism; ++j)
// temp[j].resize(nrbins);
// }
// catch(std::bad_alloc& e)
// {
// delete [] temp;
// memcritical = true;
// }
// }
// if(memcritical)
// {//memory is critical. let's try it anyway and see if the kernel frees memory.
// possibleparallelism = 1;
// temp = new std::valarray<ScalarType>(nrbins);
// }
// // std::valarray<ScalarType> temp(nrbins);//allocate a valarray for one function
// for (unsigned int k = 0; k < len; ++k)//for every index of the function that we stored
// {
// for (unsigned int j = 0; j < slices; j += possibleparallelism)//for all functionpoints
// {
// int upperlimit = possibleparallelism;
// if(j + possibleparallelism >= slices)
// upperlimit = slices - j;//determine the remaining indices
// for (unsigned int l = 0; l < nrbins; ++l)//copy function points
// {
// auto tempdatafromfile = container(l, k);
// for(int j1 = 0; j1 < upperlimit; ++j1)
// temp[j1][l] = tempdatafromfile[j+j1];
// // temp[l] = container(l, k)[j];
// }
// for(int j1 = 0; j1 < upperlimit; ++j1)
// {
// mc_analysis::errordata<ScalarType> edp(handlesign(temp[j1], avsign));//do sign - analysis
// //output data
// points.push_back(Point<FPType>((j+j1) * delta_s, k, real(edp.get_mean()), real(edp.get_error())));
// fileRe<<(j+j1) * delta_s<<" "<<real(edp.get_mean())<<" "<<real(edp.get_error())<<std::endl;
// fileIm<<(j+j1) * delta_s<<" "<<imag(edp.get_mean())<<" "<<imag(edp.get_error())<<std::endl;
// }
// }
// //output xmgrace separators
// fileRe<<"&"<<std::endl;
// fileIm<<"&"<<std::endl;
// }
// delete [] temp;
// //now let's generate the transformed data
// sort(points.begin(), points.end());
// #pragma omp parallel for schedule(dynamic)
// for(size_t i = 0; i < 5; ++i)
// {
// unsigned int inc = 2 << i;
// for(unsigned int k = 0; k < slices; k += inc)
// {
// for(unsigned int j = 0; j < len; ++j)
// {
// xdep[i]<<j<<" "<<points[k*len+j].y<<" "<<points[k*len+j].dy<<std::endl;
// }
// xdep[i]<<"&"<<std::endl;
// }
// }
// return;
// }
/**
@param name the file to which we write the data
@param len the length of the Vector
......@@ -502,7 +728,7 @@ void Error_Analysis<std::valarray<std::valarray<T> >, SignType>::analyze(const s
{
std::string outRe(p + obs.myname + string("_Re"));//create file for outputting the real part
std::string outIm(p + obs.myname + string("_Im"));//create file for outputting the imaginary part
writeVectorFunctiontoFile(obs.fv, obs.tensorindices, avsign.signcache, outRe, outIm, obs.functionpoints, obs.delta_s);//well, write to files...
writeVectorFunctiontoFile(obs.fv, obs.tensorindices, avsign.signcache, outRe, outIm, obs.functionpoints, obs.delta_s, obs.covariance);//well, write to files...
}
template<>
......
# ifndef MC_ANALYSIS_ANALYSIS_FUNCTIONS_H
# define MC_ANALYSIS_ANALYSIS_FUNCTIONS_H // include guard
template <class T>
struct ToScalar
{
typedef T RetType;
};
template <class T>
struct ToScalar<std::valarray<T> >
{
typedef typename ToScalar<T>::RetType RetType;
};
namespace mc_analysis
{
/**
......@@ -13,7 +26,7 @@ namespace mc_analysis
typename T::value_type retval(cont[0]);
for (unsigned int k = 1; k < cont.size(); ++k)
retval += cont[k];
return retval/static_cast<typename T::value_type>(cont.size());
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__)
......
......@@ -64,11 +64,13 @@ inline std::ostream& operator<<(std::ostream& out, const errordata<R>& arg )
// Specialization for valarray:
template <typename Scalar> class errordata<std::valarray<Scalar> >
{
private:
std::valarray<Scalar> mean;
std::valarray<Scalar> error;
std::valarray<Scalar> bias;
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
......@@ -89,12 +91,11 @@ public:
std::valarray<Scalar> get_mean(void)
{
return mean;
};
}
std::valarray<Scalar> get_error(void)
{
return error;
};
}
/**
An accessor to the bias of the observable
@return the bias of the observable
......@@ -109,9 +110,15 @@ public:
this->mean=errd__.mean;
this->error.resize(errd__.error.size() );
this->error=errd__.error;
this->cov.resize(errd__.cov.size());
this->cov=errd__.cov;
return *this;
}
private:
std::valarray<Scalar> mean;
std::valarray<Scalar> error;
std::valarray<Scalar> bias;
CovType cov;
};
......
......@@ -23,6 +23,7 @@
***************************************************************************/
#ifndef JACKKNIFE_H
#include <cmath>
#include <omp.h>
#include "errordata.h"
#include "analysis_functions.h"
namespace mc_analysis
......@@ -39,7 +40,7 @@ T expectationvalue(const T *const cont, unsigned int len)
T retval(cont[0]);
for (unsigned int k = 1; k < len; ++k)
retval += cont[k];
return retval/static_cast<T>(len);
return retval/static_cast<typename ToScalar<T>::RetType>(len);
}
/**
......@@ -78,7 +79,6 @@ Cont the Container that is used for the Storage of the data
template <typename Func, typename Cont>
errordata<typename Func::res_t > jackknife(Func func, const Cont *const data, unsigned int numcont)
{
bool covariance = false;
typedef typename Access_Trait<Cont>::ContainerElementType ElementType;
typedef typename Func::res_t RetType;
unsigned int len = Access_Trait<Cont>::deref(data[0]).size();//data[0].size();//we assume all datasets have the same length
......@@ -112,7 +112,75 @@ errordata<typename Func::res_t > jackknife(Func func, const Cont *const data, un
delete [] jackbin;
fun_error *= static_cast<double>(len - 1)/static_cast<double>(len);
fun_error = std::sqrt(fun_error);
return errordata<RetType>(fun_mean, fun_error, (fun_mean - jackmean) * RetType(len - 1.0));
return errordata<RetType>(fun_mean, fun_error, (fun_mean - jackmean) * RetType(len - 1));
}
/**
This function calculates the function func from the given data and performs a jackknife analysis afterwards
The template parameters are as follows:
Func the Function-Object that represents the Function we want to calculate from the data
Cont the Container that is used for the Storage of the data
@param func the function we want to evaluate
@param data an array of Containers with the stored samples
@param numcont the number of arrays pointed to by data
*/
template <typename Func, typename Cont>
errordata<typename Func::res_t > vecjackknife(Func func, const Cont& data, unsigned int numcont, bool covariance)
{
typedef typename Cont::value_type::value_type ElementType;
typedef typename Func::res_t RetType;
unsigned int len = data[0].size();//we assume all datasets have the same length. this is the number of functions not the number of bins!!!
auto mean = mc_analysis::mean(data);
//mean now contains the averages of each dataset
auto x_J = mean;
RetType* jackbin = new RetType[data.size()];
auto norm = ElementType(1.0/data.size());
//calculate jackknife samples of the function
for (unsigned int k = 0; k < data.size(); ++k)
{
x_J = norm *(ElementType(data.size()) * mean - data[k]);
jackbin[k] = func(x_J);
}
//calculate the jackknife average of the jackknife samples
RetType jackmean = expectationvalue(jackbin, data.size());
RetType fun_mean = func(mean);
RetType fun_error = (jackmean - jackbin[0]) * (jackmean - jackbin[0]);
//calculate the error
for (unsigned int k = 1; k < data.size(); ++k)
{
RetType temp = jackmean - jackbin[k];
fun_error += temp * temp;
}
auto fac = typename ToScalar<RetType>::RetType(data.size() - 1);
fun_error *= fac/static_cast<double>(data.size());
fun_error = std::sqrt(fun_error);
errordata<RetType> retval(fun_mean, fun_error, (fun_mean - jackmean) * fac);
//let's do the covariance matrix
if(covariance)
{
#ifdef _OPENMP
double start = omp_get_wtime();
#endif
std::valarray<typename ToScalar<RetType>::RetType> cov(numcont*numcont);
for(uint y = 0; y < numcont; ++y)
{
for(uint x = 0; x < numcont; ++x)
{
for(uint k = 0; k < data.size(); ++k)
{
cov[y * numcont + x] += (jackbin[k][x] - fun_mean[x]) * (jackbin[k][y] - fun_mean[y]);
}
cov[y * numcont + x] *= (data.size() - 1.0)/data.size();
}
}
#ifdef _OPENMP
std::cout<<"calculation took "<<omp_get_wtime() - start<<" s."<<std::endl;
#endif
retval.setCov(cov);
}
delete [] jackbin;
return retval;
}
}
#endif
Supports Markdown
0% or