Skip to content
Snippets Groups Projects
Select Git revision
  • bbce02aa2fe47f7b2910fad7cf07db158cb1f083
  • main default protected
  • Coverage
  • gitlab-ci-docker
  • test-install-requirementsshell-executor
  • 1.0
6 results

ValidationSimplifiedSystem.py

Blame
  • Forked from Jan Habscheid / fxdgm
    Source project has a limited visibility.
    Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    evalable.cpp 7.79 KiB
    #include "evalable.h"
    #include <sstream>
    
    void evalable :: mean(valarray<double>& v) {
    	v=mean_v[0];
    }
    
    void evalable :: error(valarray<double>& v) {
    	v=error_v[error_v.size()-1];
    }
    
    void evalable :: variance(valarray<double>& v) {
    	v=error_v[error_v.size()-1];
    	v*=v;
    	v*=2.;
    }
    
    void evalable :: autocorrelationtime(valarray<double>& v) {
    	v=error_v[error_v.size()-1];
    	v/=error_v[0];
    	v*=v;
    	v*=0.5;
    }
    
    double evalable :: mean(uint j) {
    	return mean_v[0][j];
    }
    
    double evalable :: error(uint j) {
    	return error_v[error_v.size()-1][j];
    }
    
    double evalable :: naiveerror(uint j) {
             return error_v[0][j];
    }
    
    
    double evalable :: variance(uint j) {
    	return 2.*pow(error_v[error_v.size()-1][j],2.);
    }
    
    double evalable :: autocorrelationtime(uint j) {
    	return 0.5*pow(error_v[error_v.size()-1][j]/error_v[0][j],2.);
    }
    
    void evalable :: get_statistics(ostream& os)
    {
    	streamsize origprec=os.precision();
    	for (uint j=0;j<vector_length_;++j) {
     		os << "Name = " << name_;
    		if (vector_length_>1) os<<"["<<j<<"]";
    		os<<"\n";
     		os << "Bins = " << bins() << "\n";
     		os << "BinLength = " << bin_length() << "\n";
     		os.precision((int)(fabs((log(error_v[error_v.size()-1][j])/log(10.))))+3);
     		os << "Mean  = " << mean_v[0][j]<<"\n";
     		os.precision(3);
     		os << "Error = " << error_v[error_v.size()-1][j] <<"\n";
     		os << "Variance = " << pow(error_v[error_v.size()-1][j],2.)*bins() <<"\n";
    // not done 	os << "Autocorrelationtime = " << 0.5*pow(error_v[error_v.size()-1][j]/error_v[0][j],2.) <<"\n";
    // 		luint length=1;
    // 		luint n=bins();
    // 		for (uint i=0;i<error_v.size();++i) {
    // 			os << "  ReBinLength = "<< length
    // 			   << "  Bins = " << n;
    // 			os.precision((int)(fabs((log(error_v[i][j])/log(10.))))+3);
    // 			os << "  Mean = " << mean_v[i][j];
    // 			os.precision(3);
    // 			os << "  Error = "<< error_v[i][j] 
    // 					<<"\n";
    // 			length*=binning_base;
    // 			n/=binning_base;
    // 		}
    	}	
     	os.precision(origprec);
    
    }
    
    void evalable :: jackknife(vector<observable*>& v, evalablefunc f) {
     	mean_v.clear();
     	error_v.clear();
    	bins_=v[0]->bins();
    	for (uint i=0;i<v.size();++i) if (v[i]->bins()<bins_) bins_=v[i]->bins();
    	bin_length_=v[0]->bin_length();
    	binning_base=v[0]->binning_base();
    	luint n=bins_;
    	luint nmo=n-1;
    	vector<valarray<double>* > sum;
    	vector<valarray<double>* > jack;
    	for (uint i=0;i<v.size();++i) {
    		valarray<double>* nsump=new valarray<double>(0.,v[i]->vector_length());
    		valarray<double>* njackp=new valarray<double>(0.,v[i]->vector_length());	
    		sum.push_back(nsump);
    		jack.push_back(njackp);
    	}
    	for (uint i=0;i<v.size();++i)
    		for (uint k=0;k<n;++k) {
    			*sum[i]+=v[i]->samples[slice(k*(v[i]->vector_length()),v[i]->vector_length(),1)];
    		}
    	valarray<double> me(0.,1);
    	valarray<double> jm(0.,1);
    	valarray<double> ali(0.,1);
    	for (uint k=0;k<n;++k) {
    		for (uint i=0;i<v.size();++i) {
    			*jack[i]=*sum[i];
    			*jack[i]-=v[i]->samples[slice(k*(v[i]->vector_length()),v[i]->vector_length(),1)];
    			*jack[i]/=static_cast<double>(nmo);
    		}
    		ali[0] = f(jack);
    		if (k) me+=ali; else {me=ali;}
    	}
    	jm=me;
    	jm/=static_cast<double>(n);
    	for (uint i=0;i<v.size();++i) {
    		*jack[i]=*sum[i];
    		*jack[i]/=static_cast<double>(n);
    	}
    	ali[0] = f(jack);
    	for (uint j=0;j<me.size();++j)
    		me[j]=n*ali[j]-nmo*me[j]/n;
    	mean_v.push_back(me);
    	vector_length_=me.size();
    	for (uint k=0;k<n;++k) {
    		for (uint i=0;i<v.size();++i) {
    			*jack[i]=*sum[i];
    			*jack[i]-=v[i]->samples[slice(k*(v[i]->vector_length()),v[i]->vector_length(),1)];
    			*jack[i]/=static_cast<double>(nmo);
    		}
    		ali[0] = f(jack);
    		ali-=jm;
    		ali*=ali;
    		if (k) me+=ali; else me=ali;
    	}	
    	me*=static_cast<double>(nmo)/static_cast<double>(bins_);
    	me=sqrt(me);
    	error_v.push_back(me);
    	for (uint i=0;i<v.size();++i) {
    		delete sum[i];
    		delete jack[i];
    	}
    }
    
    
    void evalable :: jackknife(vector<observable*>& v, vectorevalablefunc f) 
    {
     	mean_v.clear();
     	error_v.clear();
    	bins_=v[0]->bins();
    	for (uint i=0;i<v.size();++i) if (v[i]->bins()<bins_) bins_=v[i]->bins();
    	bin_length_=v[0]->bin_length();
    	binning_base=v[0]->binning_base();
    	luint n=bins_;
    	luint nmo=n-1;
    	vector<valarray<double>* > sum;
    	vector<valarray<double>* > jack;
    	for (uint i=0;i<v.size();++i) {
    		valarray<double>* nsump=new valarray<double>(0.,v[i]->vector_length());
    		valarray<double>* njackp=new valarray<double>(0.,v[i]->vector_length());	
    		sum.push_back(nsump);
    		jack.push_back(njackp);
    	}
    	for (uint i=0;i<v.size();++i)
    		for (uint k=0;k<n;++k) {
    			*sum[i]+=v[i]->samples[slice(k*(v[i]->vector_length()),v[i]->vector_length(),1)];
    		}
    	valarray<double> me;
    	valarray<double> jm;
    	valarray<double> ali;
    	for (uint k=0;k<n;++k) {
    		for (uint i=0;i<v.size();++i) {
    			*jack[i]=*sum[i];
    			*jack[i]-=v[i]->samples[slice(k*(v[i]->vector_length()),v[i]->vector_length(),1)];
    			*jack[i]/=static_cast<double>(nmo);
    		}
    		f(ali,jack);
    		if (k) me+=ali; else {me.resize(ali.size());me=ali;}
    	}
    	jm.resize(me.size());
    	jm=me;
    	jm/=static_cast<double>(n);
    	for (uint i=0;i<v.size();++i) {
    		*jack[i]=*sum[i];
    		*jack[i]/=static_cast<double>(n);
    	}
    	f(ali,jack);
    	for (uint j=0;j<me.size();++j)
    		me[j]=n*ali[j]-nmo*me[j]/n;
    	mean_v.push_back(me);
    	vector_length_=me.size();
    	for (uint k=0;k<n;++k) {
    		for (uint i=0;i<v.size();++i) {
    			*jack[i]=*sum[i];
    			*jack[i]-=v[i]->samples[slice(k*(v[i]->vector_length()),v[i]->vector_length(),1)];
    			*jack[i]/=static_cast<double>(nmo);
    		}
    		f(ali,jack);
    		ali-=jm;
    		ali*=ali;
    		if (k) me+=ali; else me=ali;
    	}	
    	me*=static_cast<double>(nmo)/static_cast<double>(bins_);
    	me=sqrt(me);
    	error_v.push_back(me);
    	for (uint i=0;i<v.size();++i) {
    		delete sum[i];
    		delete jack[i];
    	}
    
    // 
    // 	luint max_exp=luint(log(double(n))/log(double(binning_base))-7);
    // 	luint max_l=luint(pow(double(binning_base),double(max_exp)));
    // 	double* binneddata1 = new double[n/binning_base];
    // 	luint length=binning_base;
    // 	n/=binning_base;
    // 	nmo=n-1;	
    // 	luint index=0;
    // 	for (luint j=0;j<n;++j)	{
    // 		double ali1=0;
    // 		for (luint k=0;k<binning_base;++k) {
    // 			ali1+=o1->samples[index];
    // 			++index;
    // 		}
    // 		binneddata1[j]=ali1/binning_base;		
    // 	}
    // 	sum1=0;
    // 	for (uint j=0;j<=nmo;++j){
    // 		sum1+=binneddata1[j];
    // 	}
    // 	double meanb=0;
    // 	for (uint j=0;j<=nmo;++j){
    // 		double J1=(sum1-binneddata1[j])/nmo;
    // 		meanb+=f(J1);
    // 	}
    // 	meanb=n*f(sum1/n)-nmo*meanb/n;
    // 	double varb=0;
    // 	for (uint j=0;j<=nmo;++j){
    // 		double J1=(sum1-binneddata1[j])/nmo;
    // 		varb+=pow(f(J1)-meanb,2.);
    // 	}
    // 	varb*=((double)nmo)/((double)n);
    // 	mean_v.push_back(meanb);
    // 	error_v.push_back(sqrt(varb));
    // 	while (length<max_l) {		
    // 		length*=binning_base;
    // 		n/=binning_base;
    // 		nmo=n-1;	
    // 		index=0;
    // 		for (luint j=0;j<n;++j)	{
    // 			double ali1=0;
    // 			for (luint k=0;k<binning_base;++k) {
    // 				ali1+=binneddata1[index];
    // 				++index;
    // 			}
    // 			binneddata1[j]=ali1/binning_base;		
    // 		}
    // 		sum1=0;
    // 		for (uint j=0;j<=nmo;++j){
    // 			sum1+=binneddata1[j];
    // 		}
    // 		meanb=0;
    // 		for (uint j=0;j<=nmo;++j){
    // 			double J1=(sum1-binneddata1[j])/nmo;
    // 			meanb+=f(J1);
    // 		}
    // 		meanb=n*f(sum1/n)-nmo*meanb/n;
    // 		varb=0;
    // 		for (uint j=0;j<=nmo;++j){
    // 			double J1=(sum1-binneddata1[j])/nmo;
    // 			varb+=pow(f(J1)-meanb,2.);
    // 		}
    // 		varb*=((double)nmo)/((double)n);
    // 		mean_v.push_back(meanb);
    // 		error_v.push_back(sqrt(varb));
    // 		
    // 	}
    // 	error_=error_v[error_v.size()-1];
    // 	autocorrelationtime_=0.5*pow(error_/error_v[0],2.);
    // 	delete [] binneddata1;
    }
    
    
    
    string evalable::toString() {
    	std::stringstream ss;
    
    	get_statistics(ss);
    /*	ss << "name: " << name_ << "\n";
    	ss << "bins: " << bins_ << "\n";
    	ss << "vector length: " << vector_length_ << "\n";
    	ss << "bin length: " << bin_length_ << "\n";
    	ss << "binning base: " << binning_base << "\n";
    	ss << "--------------------\n";
    	for (size_t i = 0; i < mean_v.size(); ++i) {
    		valarray<double> *v = &mean_v[i],
    				*e = &error_v[i];
    		for (size_t j = 0; j < v->size(); ++j)
    			ss << (*v)[j] << "\t" << (*e)[j] << "\n";
    	}*/
    	return ss.str();
    }