Merge branch 'develop' into feature/rwf

This commit is contained in:
Simon Kuberski 2022-02-09 14:16:54 +01:00
commit c939e49f73
6 changed files with 47 additions and 55 deletions

View file

@ -39,7 +39,7 @@ class Corr:
region indentified for this correlator. region indentified for this correlator.
""" """
if isinstance(data_input, np.ndarray): # Input is an array of Corrs if isinstance(data_input, np.ndarray):
# This only works, if the array fulfills the conditions below # This only works, if the array fulfills the conditions below
if not len(data_input.shape) == 2 and data_input.shape[0] == data_input.shape[1]: if not len(data_input.shape) == 2 and data_input.shape[0] == data_input.shape[1]:
@ -95,7 +95,6 @@ class Corr:
# An undefined timeslice is represented by the None object # An undefined timeslice is represented by the None object
self.content = [None] * padding[0] + self.content + [None] * padding[1] self.content = [None] * padding[0] + self.content + [None] * padding[1]
self.T = len(self.content) self.T = len(self.content)
self.prange = prange self.prange = prange
self.gamma_method() self.gamma_method()
@ -160,9 +159,6 @@ class Corr:
raise Exception("Vectors are of wrong shape!") raise Exception("Vectors are of wrong shape!")
if normalize: if normalize:
vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r) vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r)
# if (not (0.95 < vector_r @ vector_r < 1.05)) or (not (0.95 < vector_l @ vector_l < 1.05)):
# print("Vectors are normalized before projection!")
newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content] newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content]
else: else:

View file

@ -35,6 +35,38 @@ def load_object(path):
return pickle.load(file) return pickle.load(file)
def pseudo_Obs(value, dvalue, name, samples=1000):
"""Generate an Obs object with given value, dvalue and name for test purposes
Parameters
----------
value : float
central value of the Obs to be generated.
dvalue : float
error of the Obs to be generated.
name : str
name of the ensemble for which the Obs is to be generated.
samples: int
number of samples for the Obs (default 1000).
"""
if dvalue <= 0.0:
return Obs([np.zeros(samples) + value], [name])
else:
for _ in range(100):
deltas = [np.random.normal(0.0, dvalue * np.sqrt(samples), samples)]
deltas -= np.mean(deltas)
deltas *= dvalue / np.sqrt((np.var(deltas) / samples)) / np.sqrt(1 + 3 / samples)
deltas += value
res = Obs(deltas, [name])
res.gamma_method(S=2, tau_exp=0)
if abs(res.dvalue - dvalue) < 1e-10 * dvalue:
break
res._value = float(value)
return res
def gen_correlated_data(means, cov, name, tau=0.5, samples=1000): def gen_correlated_data(means, cov, name, tau=0.5, samples=1000):
""" Generate observables with given covariance and autocorrelation times. """ Generate observables with given covariance and autocorrelation times.

View file

@ -1301,7 +1301,7 @@ def correlate(obs_a, obs_b):
Keep in mind to only correlate primary observables which have not been reweighted Keep in mind to only correlate primary observables which have not been reweighted
yet. The reweighting has to be applied after correlating the observables. yet. The reweighting has to be applied after correlating the observables.
Currently only works if ensembles are identical. This is not really necessary. Currently only works if ensembles are identical (this is not strictly necessary).
""" """
if sorted(obs_a.names) != sorted(obs_b.names): if sorted(obs_a.names) != sorted(obs_b.names):
@ -1461,38 +1461,6 @@ def covariance(obs1, obs2, correlation=False, **kwargs):
return dvalue return dvalue
def pseudo_Obs(value, dvalue, name, samples=1000):
"""Generate a pseudo Obs with given value, dvalue and name
Parameters
----------
value : float
central value of the Obs to be generated.
dvalue : float
error of the Obs to be generated.
name : str
name of the ensemble for which the Obs is to be generated.
samples: int
number of samples for the Obs (default 1000).
"""
if dvalue <= 0.0:
return Obs([np.zeros(samples) + value], [name])
else:
for _ in range(100):
deltas = [np.random.normal(0.0, dvalue * np.sqrt(samples), samples)]
deltas -= np.mean(deltas)
deltas *= dvalue / np.sqrt((np.var(deltas) / samples)) / np.sqrt(1 + 3 / samples)
deltas += value
res = Obs(deltas, [name])
res.gamma_method(S=2, tau_exp=0)
if abs(res.dvalue - dvalue) < 1e-10 * dvalue:
break
res._value = float(value)
return res
def import_jackknife(jacks, name, idl=None): def import_jackknife(jacks, name, idl=None):
"""Imports jackknife samples and returns an Obs """Imports jackknife samples and returns an Obs

View file

@ -1,7 +1,6 @@
import numpy as np import numpy as np
import pyerrors as pe import pyerrors as pe
import pytest import pytest
import time
np.random.seed(0) np.random.seed(0)

View file

@ -1,7 +1,5 @@
import autograd.numpy as np import autograd.numpy as np
import os import os
import random
import string
import copy import copy
import pyerrors as pe import pyerrors as pe
import pytest import pytest
@ -142,7 +140,7 @@ def test_overloading_vectorization():
assert [o.value for o in b / a] == [o.value for o in [b / p for p in a]] assert [o.value for o in b / a] == [o.value for o in [b / p for p in a]]
def test_gamma_method(): def test_gamma_method_standard_data():
for data in [np.tile([1, -1], 1000), for data in [np.tile([1, -1], 1000),
np.random.rand(100001), np.random.rand(100001),
np.zeros(1195), np.zeros(1195),
@ -285,7 +283,7 @@ def test_covariance_symmetry():
assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps) assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps)
def test_gamma_method(): def test_gamma_method_uncorrelated():
# Construct pseudo Obs with random shape # Construct pseudo Obs with random shape
value = np.random.normal(5, 10) value = np.random.normal(5, 10)
dvalue = np.abs(np.random.normal(0, 1)) dvalue = np.abs(np.random.normal(0, 1))

View file

@ -1,14 +1,14 @@
import os,sys,inspect import os
import sys
import inspect
import pyerrors as pe
import pyerrors.input.sfcf as sfin
import shutil
current_dir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) current_dir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
parent_dir = os.path.dirname(current_dir) parent_dir = os.path.dirname(current_dir)
sys.path.insert(0, parent_dir) sys.path.insert(0, parent_dir)
import pyerrors as pe
import pyerrors.input.openQCD as qcdin
import pyerrors.input.sfcf as sfin
import shutil
from time import sleep
def build_test_environment(env_type, cfgs, reps): def build_test_environment(env_type, cfgs, reps):
if env_type == "o": if env_type == "o":
@ -26,7 +26,6 @@ def build_test_environment(env_type, cfgs, reps):
def clean_test_environment(env_type, cfgs, reps): def clean_test_environment(env_type, cfgs, reps):
if env_type == "o": if env_type == "o":
for i in range(1,reps): for i in range(1,reps):