Skip to content

GitLab

Commit 989ac7b1 authored by Lukas Weber's avatar Lukas Weber
Browse files

evalable api overhaul&taking over all rng bookkeeping from mc

parent 1415004c
Hide whitespace changes
Inline Side-by-side
Showing with 121 additions and 765 deletions
evalable.cpp
View file @ 989ac7b1
......@@ -76,9 +76,7 @@ void evalable :: get_statistics(ostream& os)
}
void evalable :: jackknife(vector<observable*>& v,
void (*f) (double&, vector <valarray<double>* >&))
{
void evalable :: jackknife(vector<observable*>& v, evalablefunc f) {
mean_v.clear();
error_v.clear();
bins_=v[0]->bins();
......@@ -108,7 +106,7 @@ void evalable :: jackknife(vector<observable*>& v,
*jack[i]-=v[i]->samples[slice(k*(v[i]->vector_length()),v[i]->vector_length(),1)];
*jack[i]/=static_cast<double>(nmo);
}
f(ali[0],jack);
ali[0] = f(jack);
if (k) me+=ali; else {me=ali;}
}
jm=me;
......@@ -117,7 +115,7 @@ void evalable :: jackknife(vector<observable*>& v,
*jack[i]=*sum[i];
*jack[i]/=static_cast<double>(n);
}
f(ali[0],jack);
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);
......@@ -128,7 +126,7 @@ void evalable :: jackknife(vector<observable*>& v,
*jack[i]-=v[i]->samples[slice(k*(v[i]->vector_length()),v[i]->vector_length(),1)];
*jack[i]/=static_cast<double>(nmo);
}
f(ali[0],jack);
ali[0] = f(jack);
ali-=jm;
ali*=ali;
if (k) me+=ali; else me=ali;
......@@ -142,74 +140,8 @@ void evalable :: jackknife(vector<observable*>& v,
}
}
void evalable :: jackknife(vector<observable*>& v,
void (*f) (double&, vector <valarray<double>* >&, double*),double* p)
{
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);
}
f(ali[0],jack,p);
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);
}
f(ali[0],jack,p);
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[0],jack,p);
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 :: vectorjackknife(vector<observable*>& v,
void (*f) (valarray<double>&, vector <valarray<double>* >&))
void evalable :: jackknife(vector<observable*>& v, vectorevalablefunc f)
{
mean_v.clear();
error_v.clear();
......@@ -346,206 +278,7 @@ void evalable :: vectorjackknife(vector<observable*>& v,
// delete [] binneddata1;
}
void evalable :: vectorjackknife(vector<observable*>& v,
void (*f) (valarray<double>&, vector <valarray<double>* >&, double*),double* p)
{
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,p);
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,p);
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,p);
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, void* obj,
void (*f) (void* ,double&, vector <valarray<double>* >&))
{
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);
}
f(obj,ali[0],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);
}
f(obj,ali[0],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(obj,ali[0],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 :: vectorjackknife(vector<observable*>& v, void* obj,
void (*f) (void*, valarray<double>&, vector <valarray<double>* >&))
{
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(obj,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(obj,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(obj,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];
}
}
string evalable::toString() {
std::stringstream ss;
......
evalable.h
View file @ 989ac7b1
......@@ -7,6 +7,9 @@
#include "types.h"
#include "observable.h"
typedef std::function<double(vector<valarray<double>*>&)> evalablefunc;
typedef std::function<void(valarray<double>&, vector<valarray<double>*>&)> vectorevalablefunc;
class evalable
{
public:
......@@ -31,13 +34,9 @@ class evalable
double variance(uint);
double autocorrelationtime(uint); //need to implement jackknife binning still
void jackknife(vector<observable*>& , void (*f) (double&, vector <valarray<double>* >&));
void jackknife(vector<observable*>& , void (*f) (double&, vector <valarray<double>* >&,double*), double*);
void vectorjackknife(vector<observable*>& , void (*f) (valarray<double>&, vector <valarray<double>* >&));
void vectorjackknife(vector<observable*>& , void (*f) (valarray<double>&, vector <valarray<double>* >&,double*), double*);
void jackknife(vector<observable*>&, evalablefunc f);
void jackknife(vector<observable*>&, vectorevalablefunc f);
void jackknife(vector<observable*>& , void* , void (*f) (void*, double&, vector <valarray<double>* >&));
void vectorjackknife(vector<observable*>& , void*, void (*f) (void*, valarray<double>&, vector <valarray<double>* >&));
void reset() {mean_v.clear();error_v.clear();}
......
mc.cpp
View file @ 989ac7b1
#include "mc.h"
mc :: mc (string dir) {
abstract_mc :: abstract_mc (string dir) {
param_init(dir);
therm=param.value_or_default<int>("THERMALIZATION",10000);
}
mc :: ~mc() {
abstract_mc :: ~abstract_mc() {
random_clear();
}
void mc::random_init() {
void abstract_mc::random_init() {
if (param.defined("SEED"))
rng = new randomnumbergenerator(param.value_of<luint>("SEED"));
else
rng = new randomnumbergenerator();
have_random=true;
}
void mc::param_init(string dir) {
void abstract_mc::param_init(string dir) {
param.read_file(dir);
}
void mc::random_write(odump& d) {
void abstract_mc::random_write(odump& d) {
rng->write(d);
}
void mc::seed_write(string fn) {
void abstract_mc::seed_write(string fn) {
ofstream s;
s.open(fn.c_str());
s << rng->seed()<<endl;s.close();
}
void mc::random_read(idump& d) {
void abstract_mc::random_read(idump& d) {
rng = new randomnumbergenerator();
have_random=true;
rng->read(d);
}
void mc::random_clear() {
if(have_random) {
delete rng;
have_random=false;
}
void abstract_mc::random_clear() {
delete rng;
rng = 0;
}
double mc::random01() {
double abstract_mc::random01() {
return rng->d();
}
void abstract_mc::_init() {
random_init();
init();
}
void abstract_mc::_do_update() {
sweep++;
do_update();
}
void abstract_mc::_write(std::string dir) {
odump d(dir+"dump");
random_write(d);
d.write(sweep);
write(d);
d.close();
seed_write(dir+"seed");
dir+="sweeps";
ofstream f; f.open(dir.c_str());
f << ( (is_thermalized()) ? sweep-therm : 0 ) << " " << sweep <<endl;
f.close();
}
bool abstract_mc::_read(std::string dir) {
idump d(dir+"dump");
if (!d)
return false;
random_read(d);
d.read(sweep);
if(!read(d))
return false;
d.close();
return true;
}
void abstract_mc::_write_output(std::string filename) {
ofstream f(filename.c_str());
write_output(f);
f << "PARAMETERS" << endl;
param.get_all(f);
measure.get_statistics(f);
f.close();
}
bool abstract_mc::is_thermalized() {
return sweep>therm;
}
mc.h
View file @ 989ac7b1
......@@ -10,37 +10,45 @@
using namespace std;
class mc
class abstract_mc
{
private:
bool have_random = false;
randomnumbergenerator * rng = 0;
void param_init(string dir);
void random_clear();
protected:
void random_write(odump& d);
void seed_write(string fn);
void random_read(idump& d);
parser param;
public:
double random01();
//int sweep;
void random_init();
int sweep = 0;
int therm = 0;
protected:
parser param;
virtual void init() = 0;
virtual void write(odump &out) = 0;
virtual bool read(idump &in) = 0;
virtual void write_output(std::ofstream &f) = 0;
virtual void do_update() = 0;
public:
double random01();
virtual void do_measurement() = 0;
virtual void write(string) = 0;
virtual bool read(string) = 0;
virtual void write_output(string) = 0;
// these functions do a little more, like taking care of the
// random number generator state, then call the child class versions.
void _init();
void _write(std::string dir);
bool _read(std::string dir);
void _write_output(std::string filename);
void _do_update();
virtual bool is_thermalized() = 0;
bool is_thermalized();
measurements measure;
mc(string dir);
virtual ~mc();
abstract_mc(string dir);
virtual ~abstract_mc();
};
measurements.cpp
View file @ 989ac7b1
#include "measurements.h"
#include <sstream>
void measurements :: add_observable(std::string name)
{
observable *o = new observable(name,1);
obs_v.push_back(o);
eo[name]=obs_v.size()-1;
tag[name]=1;
}
void measurements :: add_observable(std::string name, luint bin_length)
{
observable *o = new observable(name, 1,bin_length);
obs_v.push_back(o);
eo[name]=obs_v.size()-1;
tag[name]=1;
}
void measurements :: add_observable(std::string name, luint bin_length, luint initial_length)
{
observable *o = new observable(name, 1, bin_length, initial_length);
obs_v.push_back(o);
eo[name]=obs_v.size()-1;
tag[name]=1;
}
void measurements :: add_observable(std::string name, luint bin_length, luint initial_length, luint resizefactor, luint resizeoffset)
{
observable *o = new observable(name, 1, bin_length, initial_length, resizefactor, resizeoffset);
......@@ -34,39 +10,6 @@ void measurements :: add_observable(std::string name, luint bin_length, luint in
}
void measurements :: add_vectorobservable(std::string name)
{
observable *o = new observable(name);
obs_v.push_back(o);
eo[name]=obs_v.size()-1;
tag[name]=1;
}
void measurements :: add_vectorobservable(std::string name, luint vector_length)
{
observable *o = new observable(name, vector_length);
obs_v.push_back(o);
eo[name]=obs_v.size()-1;
tag[name]=1;
}
void measurements :: add_vectorobservable(std::string name, luint vector_length, luint bin_length)
{
observable *o = new observable(name, vector_length,bin_length);
obs_v.push_back(o);
eo[name]=obs_v.size()-1;
tag[name]=1;
}
void measurements :: add_vectorobservable(std::string name, luint vector_length, luint bin_length, luint initial_length)
{
observable *o = new observable(name, vector_length, bin_length, initial_length);
obs_v.push_back(o);
eo[name]=obs_v.size()-1;
tag[name]=1;
}
void measurements :: add_vectorobservable(std::string name, luint vector_length, luint bin_length, luint initial_length, luint resizefactor, luint resizeoffset)
{
observable *o = new observable(name, vector_length, bin_length, initial_length, resizefactor, resizeoffset);
......@@ -76,9 +19,7 @@ void measurements :: add_vectorobservable(std::string name, luint vector_length,
}