2021-10-15 12:11:06 +01:00
import warnings
2021-01-08 14:46:37 +01:00
import pickle
import numpy as np
import autograd . numpy as anp # Thinly-wrapped numpy
from autograd import jacobian
import matplotlib . pyplot as plt
import numdifftools as nd
2021-10-29 12:13:04 +01:00
from itertools import groupby
2021-11-29 12:15:27 +01:00
from . covobs import Covobs
2021-01-08 14:46:37 +01:00
class Obs :
""" Class for a general observable.
Instances of Obs are the basic objects of a pyerrors error analysis .
They are initialized with a list which contains arrays of samples for
different ensembles / replica and another list of same length which contains
the names of the ensembles / replica . Mathematical operations can be
performed on instances . The result is another instance of Obs . The error of
an instance can be computed with the gamma_method . Also contains additional
methods for output and visualization of the error calculation .
Attributes
- - - - - - - - - -
2021-11-01 14:41:57 +00:00
S_global : float
Standard value for S ( default 2.0 )
S_dict : dict
Dictionary for S values . If an entry for a given ensemble
exists this overwrites the standard value for that ensemble .
tau_exp_global : float
Standard value for tau_exp ( default 0.0 )
2021-11-12 14:06:06 +00:00
tau_exp_dict : dict
2021-11-01 14:41:57 +00:00
Dictionary for tau_exp values . If an entry for a given ensemble exists
this overwrites the standard value for that ensemble .
N_sigma_global : float
Standard value for N_sigma ( default 1.0 )
2021-11-12 14:06:06 +00:00
N_sigma_dict : dict
Dictionary for N_sigma values . If an entry for a given ensemble exists
this overwrites the standard value for that ensemble .
2021-01-08 14:46:37 +01:00
"""
2021-11-30 11:08:30 +01:00
__slots__ = [ ' names ' , ' shape ' , ' r_values ' , ' deltas ' , ' N ' , ' _value ' , ' _dvalue ' ,
' ddvalue ' , ' reweighted ' , ' S ' , ' tau_exp ' , ' N_sigma ' ,
' e_dvalue ' , ' e_ddvalue ' , ' e_tauint ' , ' e_dtauint ' ,
' e_windowsize ' , ' e_rho ' , ' e_drho ' , ' e_n_tauint ' , ' e_n_dtauint ' ,
' idl ' , ' is_merged ' , ' tag ' , ' covobs ' , ' __dict__ ' ]
2021-01-08 14:46:37 +01:00
S_global = 2.0
S_dict = { }
tau_exp_global = 0.0
tau_exp_dict = { }
N_sigma_global = 1.0
2021-11-12 14:06:06 +00:00
N_sigma_dict = { }
2021-10-28 17:46:09 +02:00
filter_eps = 1e-10
2021-01-08 14:46:37 +01:00
2021-11-29 12:15:27 +01:00
def __init__ ( self , samples , names , idl = None , means = None , covobs = None , * * kwargs ) :
2021-11-01 14:41:57 +00:00
""" Initialize Obs object.
2021-11-07 21:44:22 +00:00
Parameters
2021-11-01 14:41:57 +00:00
- - - - - - - - - -
samples : list
list of numpy arrays containing the Monte Carlo samples
names : list
2021-11-15 11:13:12 +00:00
list of strings labeling the individual samples
2021-11-01 14:41:57 +00:00
idl : list , optional
list of ranges or lists on which the samples are defined
means : list , optional
list of mean values for the case that the mean values were
already subtracted from the samples
"""
2021-01-08 14:46:37 +01:00
2021-11-30 11:08:30 +01:00
if means is None and samples is not None :
2021-10-21 16:06:29 +01:00
if len ( samples ) != len ( names ) :
raise Exception ( ' Length of samples and names incompatible. ' )
2021-11-12 15:26:53 +00:00
if idl is not None :
if len ( idl ) != len ( names ) :
raise Exception ( ' Length of idl incompatible with samples and names. ' )
2021-10-21 16:06:29 +01:00
if len ( names ) != len ( set ( names ) ) :
2021-11-12 15:26:53 +00:00
raise Exception ( ' names are not unique. ' )
2021-10-21 16:06:29 +01:00
if not all ( isinstance ( x , str ) for x in names ) :
raise TypeError ( ' All names have to be strings. ' )
if min ( len ( x ) for x in samples ) < = 4 :
2021-11-12 15:26:53 +00:00
raise Exception ( ' Samples have to have at least 5 entries. ' )
2021-11-30 13:32:50 +01:00
2021-11-30 11:08:30 +01:00
if names :
2021-11-29 12:15:27 +01:00
self . names = sorted ( names )
2021-11-30 11:08:30 +01:00
else :
self . names = [ ]
2021-01-08 14:46:37 +01:00
self . shape = { }
self . r_values = { }
self . deltas = { }
2021-11-29 12:15:27 +01:00
if covobs is None :
self . covobs = { }
else :
self . covobs = covobs
2021-10-21 15:51:17 +01:00
2021-10-28 17:46:09 +02:00
self . idl = { }
2021-11-30 11:08:30 +01:00
if samples is not None :
2021-11-29 12:15:27 +01:00
if idl is not None :
for name , idx in sorted ( zip ( names , idl ) ) :
if isinstance ( idx , range ) :
self . idl [ name ] = idx
elif isinstance ( idx , ( list , np . ndarray ) ) :
dc = np . unique ( np . diff ( idx ) )
if np . any ( dc < 0 ) :
raise Exception ( " Unsorted idx for idl[ %s ] " % ( name ) )
if len ( dc ) == 1 :
self . idl [ name ] = range ( idx [ 0 ] , idx [ - 1 ] + dc [ 0 ] , dc [ 0 ] )
else :
self . idl [ name ] = list ( idx )
2021-10-28 17:46:09 +02:00
else :
2021-11-29 12:15:27 +01:00
raise Exception ( ' incompatible type for idl[ %s ]. ' % ( name ) )
else :
for name , sample in sorted ( zip ( names , samples ) ) :
self . idl [ name ] = range ( 1 , len ( sample ) + 1 )
if means is not None :
for name , sample , mean in sorted ( zip ( names , samples , means ) ) :
self . shape [ name ] = len ( self . idl [ name ] )
if len ( sample ) != self . shape [ name ] :
raise Exception ( ' Incompatible samples and idx for %s : %d vs. %d ' % ( name , len ( sample ) , self . shape [ name ] ) )
self . r_values [ name ] = mean
self . deltas [ name ] = sample
else :
for name , sample in sorted ( zip ( names , samples ) ) :
self . shape [ name ] = len ( self . idl [ name ] )
if len ( sample ) != self . shape [ name ] :
raise Exception ( ' Incompatible samples and idx for %s : %d vs. %d ' % ( name , len ( sample ) , self . shape [ name ] ) )
self . r_values [ name ] = np . mean ( sample )
self . deltas [ name ] = sample - self . r_values [ name ]
self . is_merged = { }
self . N = sum ( list ( self . shape . values ( ) ) )
self . _value = 0
if means is None :
for name in self . names :
self . _value + = self . shape [ name ] * self . r_values [ name ]
self . _value / = self . N
2021-10-21 15:51:17 +01:00
else :
2021-11-29 12:15:27 +01:00
self . _value = 0
self . is_merged = { }
self . N = 0
2021-01-08 14:46:37 +01:00
2021-10-21 16:15:02 +01:00
self . _dvalue = 0.0
2021-01-08 14:46:37 +01:00
self . ddvalue = 0.0
2021-10-29 11:56:58 +01:00
self . reweighted = False
2021-01-08 14:46:37 +01:00
2021-10-18 13:28:47 +01:00
self . tag = None
2021-10-21 16:15:02 +01:00
@property
def value ( self ) :
return self . _value
@property
def dvalue ( self ) :
return self . _dvalue
2021-11-01 17:14:30 +00:00
@property
def e_names ( self ) :
return sorted ( set ( [ o . split ( ' | ' ) [ 0 ] for o in self . names ] ) )
2021-11-30 11:08:30 +01:00
@property
def cov_names ( self ) :
return sorted ( set ( [ o for o in self . covobs . keys ( ) ] ) )
@property
def mc_names ( self ) :
return sorted ( set ( [ o . split ( ' | ' ) [ 0 ] for o in self . names if o not in self . cov_names ] ) )
2021-11-01 17:14:30 +00:00
@property
def e_content ( self ) :
res = { }
for e , e_name in enumerate ( self . e_names ) :
res [ e_name ] = sorted ( filter ( lambda x : x . startswith ( e_name + ' | ' ) , self . names ) )
if e_name in self . names :
res [ e_name ] . append ( e_name )
return res
2021-01-08 14:46:37 +01:00
def gamma_method ( self , * * kwargs ) :
2021-11-17 13:42:04 +00:00
""" Estimate the error and related properties of the Obs.
2021-01-08 14:46:37 +01:00
2021-11-07 21:44:22 +00:00
Parameters
- - - - - - - - - -
2021-11-01 14:41:57 +00:00
S : float
2021-11-17 13:42:04 +00:00
specifies a custom value for the parameter S ( default 2.0 ) .
If set to 0 it is assumed that the data exhibits no
autocorrelation . In this case the error estimates coincides
with the sample standard error .
2021-11-01 14:41:57 +00:00
tau_exp : float
positive value triggers the critical slowing down analysis
2021-11-17 13:42:04 +00:00
( default 0.0 ) .
2021-11-01 14:41:57 +00:00
N_sigma : float
number of standard deviations from zero until the tail is
2021-11-17 13:42:04 +00:00
attached to the autocorrelation function ( default 1 ) .
2021-11-01 14:41:57 +00:00
fft : bool
determines whether the fft algorithm is used for the computation
of the autocorrelation function ( default True )
2021-01-08 14:46:37 +01:00
"""
2021-11-01 17:14:30 +00:00
e_content = self . e_content
2021-01-08 14:46:37 +01:00
self . e_dvalue = { }
self . e_ddvalue = { }
self . e_tauint = { }
self . e_dtauint = { }
self . e_windowsize = { }
self . e_n_tauint = { }
self . e_n_dtauint = { }
e_gamma = { }
self . e_rho = { }
self . e_drho = { }
2021-10-21 16:15:02 +01:00
self . _dvalue = 0
2021-10-28 17:46:09 +02:00
self . ddvalue = 0
2021-01-08 14:46:37 +01:00
self . S = { }
self . tau_exp = { }
2021-11-12 14:06:06 +00:00
self . N_sigma = { }
2021-01-08 14:46:37 +01:00
if kwargs . get ( ' fft ' ) is False :
fft = False
else :
fft = True
2021-11-12 14:00:57 +00:00
def _parse_kwarg ( kwarg_name ) :
2021-11-12 13:55:03 +00:00
if kwarg_name in kwargs :
tmp = kwargs . get ( kwarg_name )
2021-01-08 14:46:37 +01:00
if isinstance ( tmp , ( int , float ) ) :
if tmp < 0 :
2021-11-12 14:06:06 +00:00
raise Exception ( kwarg_name + ' has to be larger or equal to 0. ' )
2021-01-08 14:46:37 +01:00
for e , e_name in enumerate ( self . e_names ) :
2021-11-12 13:55:03 +00:00
getattr ( self , kwarg_name ) [ e_name ] = tmp
2021-01-08 14:46:37 +01:00
else :
2021-11-12 13:55:03 +00:00
raise TypeError ( kwarg_name + ' is not in proper format. ' )
else :
for e , e_name in enumerate ( self . e_names ) :
if e_name in getattr ( Obs , kwarg_name + ' _dict ' ) :
getattr ( self , kwarg_name ) [ e_name ] = getattr ( Obs , kwarg_name + ' _dict ' ) [ e_name ]
else :
getattr ( self , kwarg_name ) [ e_name ] = getattr ( Obs , kwarg_name + ' _global ' )
2021-01-08 14:46:37 +01:00
2021-11-12 14:00:57 +00:00
_parse_kwarg ( ' S ' )
_parse_kwarg ( ' tau_exp ' )
2021-11-12 14:06:06 +00:00
_parse_kwarg ( ' N_sigma ' )
2021-01-08 14:46:37 +01:00
2021-11-30 11:08:30 +01:00
for e , e_name in enumerate ( self . mc_names ) :
r_length = [ ]
for r_name in e_content [ e_name ] :
if isinstance ( self . idl [ r_name ] , range ) :
r_length . append ( len ( self . idl [ r_name ] ) )
else :
r_length . append ( ( self . idl [ r_name ] [ - 1 ] - self . idl [ r_name ] [ 0 ] + 1 ) )
e_N = np . sum ( [ self . shape [ r_name ] for r_name in e_content [ e_name ] ] )
w_max = max ( r_length ) / / 2
e_gamma [ e_name ] = np . zeros ( w_max )
self . e_rho [ e_name ] = np . zeros ( w_max )
self . e_drho [ e_name ] = np . zeros ( w_max )
for r_name in e_content [ e_name ] :
e_gamma [ e_name ] + = self . _calc_gamma ( self . deltas [ r_name ] , self . idl [ r_name ] , self . shape [ r_name ] , w_max , fft )
gamma_div = np . zeros ( w_max )
for r_name in e_content [ e_name ] :
gamma_div + = self . _calc_gamma ( np . ones ( ( self . shape [ r_name ] ) ) , self . idl [ r_name ] , self . shape [ r_name ] , w_max , fft )
e_gamma [ e_name ] / = gamma_div [ : w_max ]
if np . abs ( e_gamma [ e_name ] [ 0 ] ) < 10 * np . finfo ( float ) . tiny : # Prevent division by zero
self . e_tauint [ e_name ] = 0.5
self . e_dtauint [ e_name ] = 0.0
self . e_dvalue [ e_name ] = 0.0
self . e_ddvalue [ e_name ] = 0.0
self . e_windowsize [ e_name ] = 0
continue
2021-01-08 14:46:37 +01:00
2021-11-30 11:08:30 +01:00
self . e_rho [ e_name ] = e_gamma [ e_name ] [ : w_max ] / e_gamma [ e_name ] [ 0 ]
self . e_n_tauint [ e_name ] = np . cumsum ( np . concatenate ( ( [ 0.5 ] , self . e_rho [ e_name ] [ 1 : ] ) ) )
# Make sure no entry of tauint is smaller than 0.5
self . e_n_tauint [ e_name ] [ self . e_n_tauint [ e_name ] < = 0.5 ] = 0.5 + np . finfo ( np . float64 ) . eps
# hep-lat/0306017 eq. (42)
self . e_n_dtauint [ e_name ] = self . e_n_tauint [ e_name ] * 2 * np . sqrt ( np . abs ( np . arange ( w_max ) + 0.5 - self . e_n_tauint [ e_name ] ) / e_N )
self . e_n_dtauint [ e_name ] [ 0 ] = 0.0
2021-01-08 14:46:37 +01:00
2021-11-30 11:08:30 +01:00
def _compute_drho ( i ) :
tmp = self . e_rho [ e_name ] [ i + 1 : w_max ] + np . concatenate ( [ self . e_rho [ e_name ] [ i - 1 : : - 1 ] , self . e_rho [ e_name ] [ 1 : w_max - 2 * i ] ] ) - 2 * self . e_rho [ e_name ] [ i ] * self . e_rho [ e_name ] [ 1 : w_max - i ]
self . e_drho [ e_name ] [ i ] = np . sqrt ( np . sum ( tmp * * 2 ) / e_N )
2021-01-08 14:46:37 +01:00
2021-11-30 11:08:30 +01:00
_compute_drho ( 1 )
if self . tau_exp [ e_name ] > 0 :
texp = self . tau_exp [ e_name ]
# if type(self.idl[e_name]) is range: # scale tau_exp according to step size
# texp /= self.idl[e_name].step
# Critical slowing down analysis
if w_max / / 2 < = 1 :
raise Exception ( " Need at least 8 samples for tau_exp error analysis " )
for n in range ( 1 , w_max / / 2 ) :
_compute_drho ( n + 1 )
if ( self . e_rho [ e_name ] [ n ] - self . N_sigma [ e_name ] * self . e_drho [ e_name ] [ n ] ) < 0 or n > = w_max / / 2 - 2 :
# Bias correction hep-lat/0306017 eq. (49) included
self . e_tauint [ e_name ] = self . e_n_tauint [ e_name ] [ n ] * ( 1 + ( 2 * n + 1 ) / e_N ) / ( 1 + 1 / e_N ) + texp * np . abs ( self . e_rho [ e_name ] [ n + 1 ] ) # The absolute makes sure, that the tail contribution is always positive
self . e_dtauint [ e_name ] = np . sqrt ( self . e_n_dtauint [ e_name ] [ n ] * * 2 + texp * * 2 * self . e_drho [ e_name ] [ n + 1 ] * * 2 )
# Error of tau_exp neglected so far, missing term: self.e_rho[e_name][n + 1] ** 2 * d_tau_exp ** 2
self . e_dvalue [ e_name ] = np . sqrt ( 2 * self . e_tauint [ e_name ] * e_gamma [ e_name ] [ 0 ] * ( 1 + 1 / e_N ) / e_N )
self . e_ddvalue [ e_name ] = self . e_dvalue [ e_name ] * np . sqrt ( ( n + 0.5 ) / e_N )
self . e_windowsize [ e_name ] = n
break
else :
if self . S [ e_name ] == 0.0 :
2021-11-17 13:42:04 +00:00
self . e_tauint [ e_name ] = 0.5
self . e_dtauint [ e_name ] = 0.0
2021-11-30 11:08:30 +01:00
self . e_dvalue [ e_name ] = np . sqrt ( e_gamma [ e_name ] [ 0 ] / ( e_N - 1 ) )
self . e_ddvalue [ e_name ] = self . e_dvalue [ e_name ] * np . sqrt ( 0.5 / e_N )
2021-11-17 13:42:04 +00:00
self . e_windowsize [ e_name ] = 0
2021-11-30 11:08:30 +01:00
else :
# Standard automatic windowing procedure
tau = self . S [ e_name ] / np . log ( ( 2 * self . e_n_tauint [ e_name ] [ 1 : ] + 1 ) / ( 2 * self . e_n_tauint [ e_name ] [ 1 : ] - 1 ) )
g_w = np . exp ( - np . arange ( 1 , w_max ) / tau ) - tau / np . sqrt ( np . arange ( 1 , w_max ) * e_N )
for n in range ( 1 , w_max ) :
if n < w_max / / 2 - 2 :
_compute_drho ( n + 1 )
if g_w [ n - 1 ] < 0 or n > = w_max - 1 :
self . e_tauint [ e_name ] = self . e_n_tauint [ e_name ] [ n ] * ( 1 + ( 2 * n + 1 ) / e_N ) / ( 1 + 1 / e_N ) # Bias correction hep-lat/0306017 eq. (49)
self . e_dtauint [ e_name ] = self . e_n_dtauint [ e_name ] [ n ]
2021-11-17 13:42:04 +00:00
self . e_dvalue [ e_name ] = np . sqrt ( 2 * self . e_tauint [ e_name ] * e_gamma [ e_name ] [ 0 ] * ( 1 + 1 / e_N ) / e_N )
self . e_ddvalue [ e_name ] = self . e_dvalue [ e_name ] * np . sqrt ( ( n + 0.5 ) / e_N )
self . e_windowsize [ e_name ] = n
break
2021-11-29 12:15:27 +01:00
2021-11-30 13:32:50 +01:00
self . _dvalue + = self . e_dvalue [ e_name ] * * 2
self . ddvalue + = ( self . e_dvalue [ e_name ] * self . e_ddvalue [ e_name ] ) * * 2
2021-01-08 14:46:37 +01:00
2021-11-30 11:08:30 +01:00
for e_name in self . cov_names :
self . e_dvalue [ e_name ] = np . sqrt ( self . covobs [ e_name ] . errsq ( ) )
self . e_ddvalue [ e_name ] = 0
self . _dvalue + = self . e_dvalue [ e_name ] * * 2
2021-01-08 14:46:37 +01:00
2021-11-30 13:32:50 +01:00
self . _dvalue = np . sqrt ( self . _dvalue )
2021-11-01 17:14:30 +00:00
if self . _dvalue == 0.0 :
2021-01-08 14:46:37 +01:00
self . ddvalue = 0.0
else :
2021-11-30 13:32:50 +01:00
self . ddvalue = np . sqrt ( self . ddvalue ) / self . _dvalue
2021-11-01 17:14:30 +00:00
return
2021-01-08 14:46:37 +01:00
2021-11-12 11:50:14 +00:00
def _calc_gamma ( self , deltas , idx , shape , w_max , fft ) :
2021-11-07 21:44:22 +00:00
""" Calculate Gamma_ {AA} from the deltas, which are defined on idx.
idx is assumed to be a contiguous range ( possibly with a stepsize != 1 )
Parameters
- - - - - - - - - -
2021-11-12 11:50:14 +00:00
deltas : list
List of fluctuations
idx : list
2021-11-15 11:13:12 +00:00
List or range of configurations on which the deltas are defined .
2021-11-12 11:50:14 +00:00
shape : int
2021-11-15 11:13:12 +00:00
Number of configurations in idx .
2021-11-12 11:50:14 +00:00
w_max : int
Upper bound for the summation window .
fft : bool
determines whether the fft algorithm is used for the computation
of the autocorrelation function .
2021-11-07 21:44:22 +00:00
"""
gamma = np . zeros ( w_max )
2021-11-12 11:41:23 +00:00
deltas = _expand_deltas ( deltas , idx , shape )
2021-11-07 21:44:22 +00:00
new_shape = len ( deltas )
if fft :
max_gamma = min ( new_shape , w_max )
# The padding for the fft has to be even
padding = new_shape + max_gamma + ( new_shape + max_gamma ) % 2
gamma [ : max_gamma ] + = np . fft . irfft ( np . abs ( np . fft . rfft ( deltas , padding ) ) * * 2 ) [ : max_gamma ]
else :
for n in range ( w_max ) :
if new_shape - n > = 0 :
gamma [ n ] + = deltas [ 0 : new_shape - n ] . dot ( deltas [ n : new_shape ] )
return gamma
2021-11-05 14:49:43 +00:00
def details ( self , ens_content = True ) :
2021-11-12 11:59:40 +00:00
""" Output detailed properties of the Obs.
Parameters
- - - - - - - - - -
ens_content : bool
print details about the ensembles and replica if true .
"""
2021-11-12 16:24:13 +00:00
if self . tag is not None :
print ( " Description: " , self . tag )
2021-11-16 11:58:27 +00:00
if not hasattr ( self , ' e_dvalue ' ) :
print ( ' Result \t %3.8e ' % ( self . value ) )
2021-01-08 14:46:37 +01:00
else :
2021-11-16 11:58:27 +00:00
if self . value == 0.0 :
percentage = np . nan
else :
2021-11-30 13:32:50 +01:00
percentage = np . abs ( self . _dvalue / self . value ) * 100
print ( ' Result \t %3.8e +/- %3.8e +/- %3.8e ( %3.3f %% ) ' % ( self . value , self . _dvalue , self . ddvalue , percentage ) )
2021-11-03 10:20:48 +00:00
if len ( self . e_names ) > 1 :
print ( ' Ensemble errors: ' )
2021-11-30 11:08:30 +01:00
for e_name in self . mc_names :
if len ( self . e_names ) > 1 :
print ( ' ' , e_name , ' \t %3.8e +/- %3.8e ' % ( self . e_dvalue [ e_name ] , self . e_ddvalue [ e_name ] ) )
if self . tau_exp [ e_name ] > 0 :
print ( ' t_int \t %3.8e +/- %3.8e tau_exp = %3.2f , N_sigma = %1.0i ' % ( self . e_tauint [ e_name ] , self . e_dtauint [ e_name ] , self . tau_exp [ e_name ] , self . N_sigma [ e_name ] ) )
2021-11-03 10:20:48 +00:00
else :
2021-11-30 11:08:30 +01:00
print ( ' t_int \t %3.8e +/- %3.8e S = %3.2f ' % ( self . e_tauint [ e_name ] , self . e_dtauint [ e_name ] , self . S [ e_name ] ) )
for e_name in self . cov_names :
print ( ' ' , e_name , ' \t %3.8e ' % ( self . e_dvalue [ e_name ] ) )
2021-11-03 10:23:31 +00:00
if ens_content is True :
2021-11-03 10:25:42 +00:00
if len ( self . e_names ) == 1 :
print ( self . N , ' samples in ' , len ( self . e_names ) , ' ensemble: ' )
else :
print ( self . N , ' samples in ' , len ( self . e_names ) , ' ensembles: ' )
2021-11-12 16:21:17 +00:00
my_string_list = [ ]
for key , value in sorted ( self . e_content . items ( ) ) :
2021-11-29 12:15:27 +01:00
if key not in self . covobs :
my_string = ' ' + " \u00B7 Ensemble ' " + key + " ' "
if len ( value ) == 1 :
my_string + = f ' : { self . shape [ value [ 0 ] ] } configurations '
if isinstance ( self . idl [ value [ 0 ] ] , range ) :
my_string + = f ' (from { self . idl [ value [ 0 ] ] . start } to { self . idl [ value [ 0 ] ] [ - 1 ] } ' + int ( self . idl [ value [ 0 ] ] . step != 1 ) * f ' in steps of { self . idl [ value [ 0 ] ] . step } ' + ' ) '
else :
my_string + = ' (irregular range) '
2021-11-12 16:21:17 +00:00
else :
2021-11-29 12:15:27 +01:00
sublist = [ ]
for v in value :
my_substring = ' ' + " \u00B7 Replicum ' " + v [ len ( key ) + 1 : ] + " ' "
my_substring + = f ' : { self . shape [ v ] } configurations '
if isinstance ( self . idl [ v ] , range ) :
my_substring + = f ' (from { self . idl [ v ] . start } to { self . idl [ v ] [ - 1 ] } ' + int ( self . idl [ v ] . step != 1 ) * f ' in steps of { self . idl [ v ] . step } ' + ' ) '
else :
my_substring + = ' (irregular range) '
sublist . append ( my_substring )
my_string + = ' \n ' + ' \n ' . join ( sublist )
2021-11-12 16:21:17 +00:00
else :
2021-11-29 12:15:27 +01:00
my_string = ' ' + " \u00B7 Covobs ' " + key + " ' "
2021-11-12 16:21:17 +00:00
my_string_list . append ( my_string )
print ( ' \n ' . join ( my_string_list ) )
2021-01-08 14:46:37 +01:00
2021-11-12 14:25:51 +00:00
def print ( self , level = 1 ) :
warnings . warn ( " Method ' print ' renamed to ' details ' " , DeprecationWarning )
self . details ( level > 1 )
2021-01-08 14:46:37 +01:00
2021-10-20 13:12:58 +01:00
def is_zero_within_error ( self , sigma = 1 ) :
2021-10-26 10:00:18 +01:00
""" Checks whether the observable is zero within ' sigma ' standard errors.
2021-10-20 13:12:58 +01:00
2021-11-12 11:59:40 +00:00
Parameters
- - - - - - - - - -
sigma : int
Number of standard errors used for the check .
2021-10-20 13:12:58 +01:00
Works only properly when the gamma method was run .
"""
2021-11-30 13:32:50 +01:00
return self . is_zero ( ) or np . abs ( self . value ) < = sigma * self . _dvalue
2021-10-15 11:47:41 +01:00
2021-11-18 10:51:46 +00:00
def is_zero ( self , rtol = 1.e-5 , atol = 1.e-8 ) :
""" Checks whether the observable is zero within a given tolerance.
Parameters
- - - - - - - - - -
rtol : float
Relative tolerance ( for details see numpy documentation ) .
atol : float
Absolute tolerance ( for details see numpy documentation ) .
"""
2021-12-01 15:14:56 +01:00
return ( np . isclose ( 0.0 , self . value , rtol , atol ) and all ( np . allclose ( 0.0 , delta , rtol , atol ) for delta in self . deltas . values ( ) ) and all ( np . allclose ( 0.0 , delta . grad , rtol , atol ) for delta in self . covobs . values ( ) ) )
2021-10-17 12:28:59 +01:00
2021-09-30 12:43:28 +01:00
def plot_tauint ( self , save = None ) :
2021-11-12 11:59:40 +00:00
""" Plot integrated autocorrelation time for each ensemble.
Parameters
- - - - - - - - - -
save : str
saves the figure to a file named ' save ' if .
"""
2021-11-16 11:58:27 +00:00
if not hasattr ( self , ' e_dvalue ' ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' Run the gamma method first. ' )
2021-09-30 12:43:28 +01:00
fig = plt . figure ( )
2021-11-30 11:08:30 +01:00
for e , e_name in enumerate ( self . mc_names ) :
2021-09-30 12:43:28 +01:00
plt . xlabel ( r ' $W$ ' )
plt . ylabel ( r ' $ \ tau_ \ mathrm {int} $ ' )
2021-01-08 14:46:37 +01:00
length = int ( len ( self . e_n_tauint [ e_name ] ) )
if self . tau_exp [ e_name ] > 0 :
base = self . e_n_tauint [ e_name ] [ self . e_windowsize [ e_name ] ]
x_help = np . arange ( 2 * self . tau_exp [ e_name ] )
2021-10-11 12:22:58 +01:00
y_help = ( x_help + 1 ) * np . abs ( self . e_rho [ e_name ] [ self . e_windowsize [ e_name ] + 1 ] ) * ( 1 - x_help / ( 2 * ( 2 * self . tau_exp [ e_name ] - 1 ) ) ) + base
2021-01-08 14:46:37 +01:00
x_arr = np . arange ( self . e_windowsize [ e_name ] + 1 , self . e_windowsize [ e_name ] + 1 + 2 * self . tau_exp [ e_name ] )
2021-09-30 12:43:28 +01:00
plt . plot ( x_arr , y_help , ' C ' + str ( e ) , linewidth = 1 , ls = ' -- ' , marker = ' , ' )
2021-01-08 14:46:37 +01:00
plt . errorbar ( [ self . e_windowsize [ e_name ] + 2 * self . tau_exp [ e_name ] ] , [ self . e_tauint [ e_name ] ] ,
2021-09-30 12:43:28 +01:00
yerr = [ self . e_dtauint [ e_name ] ] , fmt = ' C ' + str ( e ) , linewidth = 1 , capsize = 2 , marker = ' o ' , mfc = plt . rcParams [ ' axes.facecolor ' ] )
2021-01-08 14:46:37 +01:00
xmax = self . e_windowsize [ e_name ] + 2 * self . tau_exp [ e_name ] + 1.5
2021-10-11 12:22:58 +01:00
label = e_name + r ' , $ \ tau_ \ mathrm {exp} $= ' + str ( np . around ( self . tau_exp [ e_name ] , decimals = 2 ) )
2021-01-08 14:46:37 +01:00
else :
2021-10-11 12:22:58 +01:00
label = e_name + ' , S= ' + str ( np . around ( self . S [ e_name ] , decimals = 2 ) )
2021-01-08 14:46:37 +01:00
xmax = max ( 10.5 , 2 * self . e_windowsize [ e_name ] - 0.5 )
2021-09-30 12:43:28 +01:00
plt . errorbar ( np . arange ( length ) , self . e_n_tauint [ e_name ] [ : ] , yerr = self . e_n_dtauint [ e_name ] [ : ] , linewidth = 1 , capsize = 2 , label = label )
plt . axvline ( x = self . e_windowsize [ e_name ] , color = ' C ' + str ( e ) , alpha = 0.5 , marker = ' , ' , ls = ' -- ' )
plt . legend ( )
2021-01-08 14:46:37 +01:00
plt . xlim ( - 0.5 , xmax )
2021-09-30 12:43:28 +01:00
plt . ylim ( bottom = 0.0 )
plt . draw ( )
if save :
fig . savefig ( save )
2021-01-08 14:46:37 +01:00
def plot_rho ( self ) :
""" Plot normalized autocorrelation function time for each ensemble. """
2021-11-16 11:58:27 +00:00
if not hasattr ( self , ' e_dvalue ' ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' Run the gamma method first. ' )
2021-11-30 11:08:30 +01:00
for e , e_name in enumerate ( self . mc_names ) :
2021-01-08 14:46:37 +01:00
plt . xlabel ( ' W ' )
plt . ylabel ( ' rho ' )
length = int ( len ( self . e_drho [ e_name ] ) )
plt . errorbar ( np . arange ( length ) , self . e_rho [ e_name ] [ : length ] , yerr = self . e_drho [ e_name ] [ : ] , linewidth = 1 , capsize = 2 )
2021-09-28 12:44:58 +01:00
plt . axvline ( x = self . e_windowsize [ e_name ] , color = ' r ' , alpha = 0.25 , ls = ' -- ' , marker = ' , ' )
2021-01-08 14:46:37 +01:00
if self . tau_exp [ e_name ] > 0 :
plt . plot ( [ self . e_windowsize [ e_name ] + 1 , self . e_windowsize [ e_name ] + 1 + 2 * self . tau_exp [ e_name ] ] ,
[ self . e_rho [ e_name ] [ self . e_windowsize [ e_name ] + 1 ] , 0 ] , ' k- ' , lw = 1 )
xmax = self . e_windowsize [ e_name ] + 2 * self . tau_exp [ e_name ] + 1.5
2021-10-11 12:22:58 +01:00
plt . title ( ' Rho ' + e_name + r ' , tau \ _exp= ' + str ( np . around ( self . tau_exp [ e_name ] , decimals = 2 ) ) )
2021-01-08 14:46:37 +01:00
else :
xmax = max ( 10.5 , 2 * self . e_windowsize [ e_name ] - 0.5 )
plt . title ( ' Rho ' + e_name + ' , S= ' + str ( np . around ( self . S [ e_name ] , decimals = 2 ) ) )
plt . plot ( [ - 0.5 , xmax ] , [ 0 , 0 ] , ' k-- ' , lw = 1 )
plt . xlim ( - 0.5 , xmax )
2021-09-30 12:43:28 +01:00
plt . draw ( )
2021-01-08 14:46:37 +01:00
def plot_rep_dist ( self ) :
""" Plot replica distribution for each ensemble with more than one replicum. """
2021-11-16 11:58:27 +00:00
if not hasattr ( self , ' e_dvalue ' ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' Run the gamma method first. ' )
2021-11-30 11:08:30 +01:00
for e , e_name in enumerate ( self . mc_names ) :
2021-01-08 14:46:37 +01:00
if len ( self . e_content [ e_name ] ) == 1 :
print ( ' No replica distribution for a single replicum ( ' , e_name , ' ) ' )
continue
r_length = [ ]
sub_r_mean = 0
for r , r_name in enumerate ( self . e_content [ e_name ] ) :
r_length . append ( len ( self . deltas [ r_name ] ) )
sub_r_mean + = self . shape [ r_name ] * self . r_values [ r_name ]
e_N = np . sum ( r_length )
sub_r_mean / = e_N
arr = np . zeros ( len ( self . e_content [ e_name ] ) )
for r , r_name in enumerate ( self . e_content [ e_name ] ) :
arr [ r ] = ( self . r_values [ r_name ] - sub_r_mean ) / ( self . e_dvalue [ e_name ] * np . sqrt ( e_N / self . shape [ r_name ] - 1 ) )
plt . hist ( arr , rwidth = 0.8 , bins = len ( self . e_content [ e_name ] ) )
2021-10-18 10:08:48 +01:00
plt . title ( ' Replica distribution ' + e_name + ' (mean=0, var=1) ' )
2021-10-15 15:46:59 +01:00
plt . draw ( )
2021-01-08 14:46:37 +01:00
2021-10-28 17:46:09 +02:00
def plot_history ( self , expand = True ) :
2021-11-12 11:59:40 +00:00
""" Plot derived Monte Carlo history for each ensemble
Parameters
- - - - - - - - - -
expand : bool
show expanded history for irregular Monte Carlo chains ( default : True ) .
"""
2021-11-30 11:08:30 +01:00
for e , e_name in enumerate ( self . mc_names ) :
2021-10-11 12:22:58 +01:00
plt . figure ( )
2021-01-08 14:46:37 +01:00
r_length = [ ]
2021-10-28 17:46:09 +02:00
tmp = [ ]
2021-01-08 14:46:37 +01:00
for r , r_name in enumerate ( self . e_content [ e_name ] ) :
2021-10-28 17:46:09 +02:00
if expand :
2021-11-16 12:27:34 +00:00
tmp . append ( _expand_deltas ( self . deltas [ r_name ] , list ( self . idl [ r_name ] ) , self . shape [ r_name ] ) + self . r_values [ r_name ] )
2021-10-28 17:46:09 +02:00
else :
tmp . append ( self . deltas [ r_name ] + self . r_values [ r_name ] )
r_length . append ( len ( tmp [ - 1 ] ) )
2021-01-08 14:46:37 +01:00
e_N = np . sum ( r_length )
x = np . arange ( e_N )
y = np . concatenate ( tmp , axis = 0 )
plt . errorbar ( x , y , fmt = ' . ' , markersize = 3 )
plt . xlim ( - 0.5 , e_N - 0.5 )
plt . title ( e_name )
2021-10-15 15:46:59 +01:00
plt . draw ( )
2021-01-08 14:46:37 +01:00
def plot_piechart ( self ) :
""" Plot piechart which shows the fractional contribution of each
ensemble to the error and returns a dictionary containing the fractions . """
2021-11-16 11:58:27 +00:00
if not hasattr ( self , ' e_dvalue ' ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' Run the gamma method first. ' )
2021-11-30 13:32:50 +01:00
if self . _dvalue == 0.0 :
2021-01-08 14:46:37 +01:00
raise Exception ( ' Error is 0.0 ' )
labels = self . e_names
2021-12-01 15:14:56 +01:00
sizes = [ self . e_dvalue [ name ] * * 2 for name in labels ] / self . _dvalue * * 2
2021-01-08 14:46:37 +01:00
fig1 , ax1 = plt . subplots ( )
ax1 . pie ( sizes , labels = labels , startangle = 90 , normalize = True )
ax1 . axis ( ' equal ' )
2021-10-15 15:46:59 +01:00
plt . draw ( )
2021-01-08 14:46:37 +01:00
return dict ( zip ( self . e_names , sizes ) )
def dump ( self , name , * * kwargs ) :
""" Dump the Obs to a pickle file ' name ' .
2021-11-07 21:44:22 +00:00
Parameters
- - - - - - - - - -
2021-11-12 11:59:40 +00:00
name : str
name of the file to be saved .
2021-11-07 21:44:22 +00:00
path : str
specifies a custom path for the file ( default ' . ' )
2021-01-08 14:46:37 +01:00
"""
if ' path ' in kwargs :
file_name = kwargs . get ( ' path ' ) + ' / ' + name + ' .p '
else :
file_name = name + ' .p '
with open ( file_name , ' wb ' ) as fb :
pickle . dump ( self , fb )
2021-11-15 14:11:45 +00:00
def export_jackknife ( self ) :
""" Export jackknife samples from the Obs
Returns
- - - - - - -
2021-11-15 14:18:26 +00:00
numpy . ndarray
2021-11-15 14:11:45 +00:00
Returns a numpy array of length N + 1 where N is the number of samples
for the given ensemble and replicum . The zeroth entry of the array contains
the mean value of the Obs , entries 1 to N contain the N jackknife samples
derived from the Obs . The current implementation only works for observables
2021-11-15 14:38:21 +00:00
defined on exactly one ensemble and replicum . The derived jackknife samples
should agree with samples from a full jackknife analysis up to O ( 1 / N ) .
2021-11-15 14:11:45 +00:00
"""
if len ( self . names ) != 1 :
raise Exception ( " ' export_jackknife ' is only implemented for Obs defined on one ensemble and replicum. " )
name = self . names [ 0 ]
full_data = self . deltas [ name ] + self . r_values [ name ]
n = full_data . size
mean = np . mean ( full_data )
tmp_jacks = np . zeros ( n + 1 )
tmp_jacks [ 0 ] = self . value
2021-11-17 16:25:50 +00:00
tmp_jacks [ 1 : ] = ( n * mean - full_data ) / ( n - 1 )
2021-11-15 14:11:45 +00:00
return tmp_jacks
2021-10-16 11:41:44 +01:00
def __float__ ( self ) :
2021-10-18 09:34:10 +01:00
return float ( self . value )
2021-10-16 11:41:44 +01:00
2021-01-08 14:46:37 +01:00
def __repr__ ( self ) :
2021-10-16 11:21:19 +01:00
return ' Obs[ ' + str ( self ) + ' ] '
def __str__ ( self ) :
2021-11-30 13:32:50 +01:00
if self . _dvalue == 0.0 :
2021-10-16 11:02:54 +01:00
return str ( self . value )
2021-11-30 13:32:50 +01:00
fexp = np . floor ( np . log10 ( self . _dvalue ) )
2021-01-08 14:46:37 +01:00
if fexp < 0.0 :
2021-11-30 13:32:50 +01:00
return ' { : {form} }( {:2.0f} ) ' . format ( self . value , self . _dvalue * 10 * * ( - fexp + 1 ) , form = ' . ' + str ( - int ( fexp ) + 1 ) + ' f ' )
2021-01-08 14:46:37 +01:00
elif fexp == 0.0 :
2021-11-30 13:32:50 +01:00
return ' {:.1f} ( {:1.1f} ) ' . format ( self . value , self . _dvalue )
2021-01-08 14:46:37 +01:00
else :
2021-11-30 13:32:50 +01:00
return ' {:.0f} ( {:2.0f} ) ' . format ( self . value , self . _dvalue )
2021-10-16 11:02:54 +01:00
2021-01-08 14:46:37 +01:00
# Overload comparisons
def __lt__ ( self , other ) :
return self . value < other
2021-10-18 13:07:33 +01:00
def __le__ ( self , other ) :
return self . value < = other
2021-01-08 14:46:37 +01:00
def __gt__ ( self , other ) :
return self . value > other
2021-10-18 13:07:33 +01:00
def __ge__ ( self , other ) :
return self . value > = other
def __eq__ ( self , other ) :
return ( self - other ) . is_zero ( )
def __ne__ ( self , other ) :
return not ( self - other ) . is_zero ( )
2021-01-08 14:46:37 +01:00
# Overload math operations
def __add__ ( self , y ) :
if isinstance ( y , Obs ) :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] + x [ 1 ] , [ self , y ] , man_grad = [ 1 , 1 ] )
else :
if isinstance ( y , np . ndarray ) :
return np . array ( [ self + o for o in y ] )
2021-09-29 15:07:35 +01:00
elif y . __class__ . __name__ == ' Corr ' :
return NotImplemented
2021-01-08 14:46:37 +01:00
else :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] + y , [ self ] , man_grad = [ 1 ] )
2021-10-11 12:22:58 +01:00
2021-01-08 14:46:37 +01:00
def __radd__ ( self , y ) :
return self + y
def __mul__ ( self , y ) :
if isinstance ( y , Obs ) :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] * x [ 1 ] , [ self , y ] , man_grad = [ y . value , self . value ] )
else :
if isinstance ( y , np . ndarray ) :
return np . array ( [ self * o for o in y ] )
2021-10-17 12:23:27 +01:00
elif isinstance ( y , complex ) :
return CObs ( self * y . real , self * y . imag )
2021-09-29 15:07:35 +01:00
elif y . __class__ . __name__ == ' Corr ' :
return NotImplemented
2021-01-08 14:46:37 +01:00
else :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] * y , [ self ] , man_grad = [ y ] )
def __rmul__ ( self , y ) :
return self * y
def __sub__ ( self , y ) :
if isinstance ( y , Obs ) :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] - x [ 1 ] , [ self , y ] , man_grad = [ 1 , - 1 ] )
else :
if isinstance ( y , np . ndarray ) :
return np . array ( [ self - o for o in y ] )
2021-09-29 15:07:35 +01:00
elif y . __class__ . __name__ == ' Corr ' :
return NotImplemented
2021-01-08 14:46:37 +01:00
else :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] - y , [ self ] , man_grad = [ 1 ] )
def __rsub__ ( self , y ) :
return - 1 * ( self - y )
def __neg__ ( self ) :
return - 1 * self
def __truediv__ ( self , y ) :
if isinstance ( y , Obs ) :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] / x [ 1 ] , [ self , y ] , man_grad = [ 1 / y . value , - self . value / y . value * * 2 ] )
else :
if isinstance ( y , np . ndarray ) :
return np . array ( [ self / o for o in y ] )
2021-09-29 15:07:35 +01:00
elif y . __class__ . __name__ == ' Corr ' :
2021-09-27 16:39:33 +01:00
return NotImplemented
2021-01-08 14:46:37 +01:00
else :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] / y , [ self ] , man_grad = [ 1 / y ] )
def __rtruediv__ ( self , y ) :
if isinstance ( y , Obs ) :
return derived_observable ( lambda x , * * kwargs : x [ 0 ] / x [ 1 ] , [ y , self ] , man_grad = [ 1 / self . value , - y . value / self . value * * 2 ] )
else :
if isinstance ( y , np . ndarray ) :
return np . array ( [ o / self for o in y ] )
2021-11-03 11:53:26 +00:00
elif y . __class__ . __name__ == ' Corr ' :
return NotImplemented
2021-01-08 14:46:37 +01:00
else :
return derived_observable ( lambda x , * * kwargs : y / x [ 0 ] , [ self ] , man_grad = [ - y / self . value * * 2 ] )
def __pow__ ( self , y ) :
if isinstance ( y , Obs ) :
return derived_observable ( lambda x : x [ 0 ] * * x [ 1 ] , [ self , y ] )
else :
return derived_observable ( lambda x : x [ 0 ] * * y , [ self ] )
def __rpow__ ( self , y ) :
if isinstance ( y , Obs ) :
return derived_observable ( lambda x : x [ 0 ] * * x [ 1 ] , [ y , self ] )
else :
return derived_observable ( lambda x : y * * x [ 0 ] , [ self ] )
def __abs__ ( self ) :
return derived_observable ( lambda x : anp . abs ( x [ 0 ] ) , [ self ] )
# Overload numpy functions
def sqrt ( self ) :
return derived_observable ( lambda x , * * kwargs : np . sqrt ( x [ 0 ] ) , [ self ] , man_grad = [ 1 / 2 / np . sqrt ( self . value ) ] )
def log ( self ) :
return derived_observable ( lambda x , * * kwargs : np . log ( x [ 0 ] ) , [ self ] , man_grad = [ 1 / self . value ] )
def exp ( self ) :
return derived_observable ( lambda x , * * kwargs : np . exp ( x [ 0 ] ) , [ self ] , man_grad = [ np . exp ( self . value ) ] )
def sin ( self ) :
return derived_observable ( lambda x , * * kwargs : np . sin ( x [ 0 ] ) , [ self ] , man_grad = [ np . cos ( self . value ) ] )
def cos ( self ) :
return derived_observable ( lambda x , * * kwargs : np . cos ( x [ 0 ] ) , [ self ] , man_grad = [ - np . sin ( self . value ) ] )
def tan ( self ) :
return derived_observable ( lambda x , * * kwargs : np . tan ( x [ 0 ] ) , [ self ] , man_grad = [ 1 / np . cos ( self . value ) * * 2 ] )
def arcsin ( self ) :
return derived_observable ( lambda x : anp . arcsin ( x [ 0 ] ) , [ self ] )
def arccos ( self ) :
return derived_observable ( lambda x : anp . arccos ( x [ 0 ] ) , [ self ] )
def arctan ( self ) :
return derived_observable ( lambda x : anp . arctan ( x [ 0 ] ) , [ self ] )
def sinh ( self ) :
return derived_observable ( lambda x , * * kwargs : np . sinh ( x [ 0 ] ) , [ self ] , man_grad = [ np . cosh ( self . value ) ] )
def cosh ( self ) :
return derived_observable ( lambda x , * * kwargs : np . cosh ( x [ 0 ] ) , [ self ] , man_grad = [ np . sinh ( self . value ) ] )
def tanh ( self ) :
return derived_observable ( lambda x , * * kwargs : np . tanh ( x [ 0 ] ) , [ self ] , man_grad = [ 1 / np . cosh ( self . value ) * * 2 ] )
def arcsinh ( self ) :
return derived_observable ( lambda x : anp . arcsinh ( x [ 0 ] ) , [ self ] )
def arccosh ( self ) :
return derived_observable ( lambda x : anp . arccosh ( x [ 0 ] ) , [ self ] )
def arctanh ( self ) :
return derived_observable ( lambda x : anp . arctanh ( x [ 0 ] ) , [ self ] )
def sinc ( self ) :
return derived_observable ( lambda x : anp . sinc ( x [ 0 ] ) , [ self ] )
2021-10-16 12:07:40 +01:00
class CObs :
2021-10-18 13:34:02 +01:00
""" Class for a complex valued observable. """
2021-10-21 10:59:53 +01:00
__slots__ = [ ' _real ' , ' _imag ' , ' tag ' ]
2021-10-18 13:34:02 +01:00
2021-10-16 12:07:40 +01:00
def __init__ ( self , real , imag = 0.0 ) :
2021-10-21 10:59:53 +01:00
self . _real = real
self . _imag = imag
2021-10-18 13:31:57 +01:00
self . tag = None
2021-10-16 12:07:40 +01:00
2021-10-21 10:59:53 +01:00
@property
def real ( self ) :
return self . _real
@property
def imag ( self ) :
return self . _imag
2021-10-16 12:07:40 +01:00
def gamma_method ( self , * * kwargs ) :
2021-10-26 17:12:25 +01:00
""" Executes the gamma_method for the real and the imaginary part. """
2021-10-17 11:34:53 +01:00
if isinstance ( self . real , Obs ) :
self . real . gamma_method ( * * kwargs )
if isinstance ( self . imag , Obs ) :
self . imag . gamma_method ( * * kwargs )
2021-10-16 12:07:40 +01:00
2021-10-18 17:36:49 +01:00
def is_zero ( self ) :
2021-10-26 17:12:25 +01:00
""" Checks whether both real and imaginary part are zero within machine precision. """
2021-10-18 17:36:49 +01:00
return self . real == 0.0 and self . imag == 0.0
2021-10-18 14:48:16 +01:00
def conjugate ( self ) :
return CObs ( self . real , - self . imag )
2021-10-16 12:07:40 +01:00
def __add__ ( self , other ) :
2021-10-25 15:10:18 +01:00
if isinstance ( other , np . ndarray ) :
return other + self
elif hasattr ( other , ' real ' ) and hasattr ( other , ' imag ' ) :
2021-10-16 12:07:40 +01:00
return CObs ( self . real + other . real ,
self . imag + other . imag )
else :
return CObs ( self . real + other , self . imag )
def __radd__ ( self , y ) :
return self + y
def __sub__ ( self , other ) :
2021-10-25 15:10:18 +01:00
if isinstance ( other , np . ndarray ) :
return - 1 * ( other - self )
elif hasattr ( other , ' real ' ) and hasattr ( other , ' imag ' ) :
2021-10-16 12:07:40 +01:00
return CObs ( self . real - other . real , self . imag - other . imag )
else :
return CObs ( self . real - other , self . imag )
def __rsub__ ( self , other ) :
return - 1 * ( self - other )
def __mul__ ( self , other ) :
2021-10-25 15:10:18 +01:00
if isinstance ( other , np . ndarray ) :
return other * self
elif hasattr ( other , ' real ' ) and hasattr ( other , ' imag ' ) :
2021-10-25 09:41:34 +01:00
if all ( isinstance ( i , Obs ) for i in [ self . real , self . imag , other . real , other . imag ] ) :
return CObs ( derived_observable ( lambda x , * * kwargs : x [ 0 ] * x [ 1 ] - x [ 2 ] * x [ 3 ] ,
[ self . real , other . real , self . imag , other . imag ] ,
man_grad = [ other . real . value , self . real . value , - other . imag . value , - self . imag . value ] ) ,
derived_observable ( lambda x , * * kwargs : x [ 2 ] * x [ 1 ] + x [ 0 ] * x [ 3 ] ,
[ self . real , other . real , self . imag , other . imag ] ,
man_grad = [ other . imag . value , self . imag . value , other . real . value , self . real . value ] ) )
2021-10-25 09:51:30 +01:00
elif getattr ( other , ' imag ' , 0 ) != 0 :
2021-10-25 09:41:34 +01:00
return CObs ( self . real * other . real - self . imag * other . imag ,
self . imag * other . real + self . real * other . imag )
2021-10-25 09:51:30 +01:00
else :
return CObs ( self . real * other . real , self . imag * other . real )
2021-10-16 12:07:40 +01:00
else :
2021-10-25 09:51:30 +01:00
return CObs ( self . real * other , self . imag * other )
2021-10-16 12:07:40 +01:00
def __rmul__ ( self , other ) :
return self * other
def __truediv__ ( self , other ) :
2021-10-25 15:10:18 +01:00
if isinstance ( other , np . ndarray ) :
return 1 / ( other / self )
elif hasattr ( other , ' real ' ) and hasattr ( other , ' imag ' ) :
2021-10-16 12:07:40 +01:00
r = other . real * * 2 + other . imag * * 2
return CObs ( ( self . real * other . real + self . imag * other . imag ) / r , ( self . imag * other . real - self . real * other . imag ) / r )
else :
return CObs ( self . real / other , self . imag / other )
2021-10-25 10:44:51 +01:00
def __rtruediv__ ( self , other ) :
r = self . real * * 2 + self . imag * * 2
if hasattr ( other , ' real ' ) and hasattr ( other , ' imag ' ) :
return CObs ( ( self . real * other . real + self . imag * other . imag ) / r , ( self . real * other . imag - self . imag * other . real ) / r )
else :
return CObs ( self . real * other / r , - self . imag * other / r )
2021-10-16 12:07:40 +01:00
def __abs__ ( self ) :
return np . sqrt ( self . real * * 2 + self . imag * * 2 )
def __neg__ ( other ) :
return - 1 * other
def __eq__ ( self , other ) :
return self . real == other . real and self . imag == other . imag
def __str__ ( self ) :
2021-10-18 14:48:16 +01:00
return ' ( ' + str ( self . real ) + int ( self . imag > = 0.0 ) * ' + ' + str ( self . imag ) + ' j) '
2021-10-16 12:07:40 +01:00
def __repr__ ( self ) :
return ' CObs[ ' + str ( self ) + ' ] '
2021-11-12 11:41:23 +00:00
def _expand_deltas ( deltas , idx , shape ) :
""" Expand deltas defined on idx to a regular, contiguous range, where holes are filled by 0.
If idx is of type range , the deltas are not changed
Parameters
- - - - - - - - - -
deltas : list
List of fluctuations
idx : list
List or range of configs on which the deltas are defined .
shape : int
Number of configs in idx .
"""
if isinstance ( idx , range ) :
return deltas
else :
ret = np . zeros ( idx [ - 1 ] - idx [ 0 ] + 1 )
for i in range ( shape ) :
ret [ idx [ i ] - idx [ 0 ] ] = deltas [ i ]
return ret
2021-11-11 10:50:22 +00:00
def _merge_idx ( idl ) :
2021-10-28 17:46:09 +02:00
""" Returns the union of all lists in idl
Parameters
- - - - - - - - - -
2021-11-07 21:44:22 +00:00
idl : list
List of lists or ranges .
2021-10-28 17:46:09 +02:00
"""
2021-10-29 12:13:04 +01:00
# Use groupby to efficiently check whether all elements of idl are identical
2021-10-29 15:52:18 +01:00
try :
g = groupby ( idl )
if next ( g , True ) and not next ( g , False ) :
return idl [ 0 ]
except :
pass
2021-10-29 12:13:04 +01:00
2021-10-28 17:46:09 +02:00
if np . all ( [ type ( idx ) is range for idx in idl ] ) :
if len ( set ( [ idx [ 0 ] for idx in idl ] ) ) == 1 :
idstart = min ( [ idx . start for idx in idl ] )
idstop = max ( [ idx . stop for idx in idl ] )
idstep = min ( [ idx . step for idx in idl ] )
return range ( idstart , idstop , idstep )
return list ( set ( ) . union ( * idl ) )
2021-11-11 10:50:22 +00:00
def _expand_deltas_for_merge ( deltas , idx , shape , new_idx ) :
2021-10-28 17:46:09 +02:00
""" Expand deltas defined on idx to the list of configs that is defined by new_idx.
New , empy entries are filled by 0. If idx and new_idx are of type range , the smallest
common divisor of the step sizes is used as new step size .
Parameters
- - - - - - - - - -
2021-11-01 14:41:57 +00:00
deltas : list
List of fluctuations
idx : list
List or range of configs on which the deltas are defined .
Has to be a subset of new_idx .
shape : list
Number of configs in idx .
new_idx : list
List of configs that defines the new range .
2021-10-28 17:46:09 +02:00
"""
if type ( idx ) is range and type ( new_idx ) is range :
if idx == new_idx :
return deltas
ret = np . zeros ( new_idx [ - 1 ] - new_idx [ 0 ] + 1 )
for i in range ( shape ) :
ret [ idx [ i ] - new_idx [ 0 ] ] = deltas [ i ]
return np . array ( [ ret [ new_idx [ i ] - new_idx [ 0 ] ] for i in range ( len ( new_idx ) ) ] )
2021-11-17 16:44:54 +01:00
def _filter_zeroes ( deltas , idx , eps = Obs . filter_eps ) :
2021-10-28 17:46:09 +02:00
""" Filter out all configurations with vanishing fluctuation such that they do not
2021-11-17 16:44:54 +01:00
contribute to the error estimate anymore . Returns the new deltas and
idx according to the filtering .
2021-10-28 17:46:09 +02:00
A fluctuation is considered to be vanishing , if it is smaller than eps times
the mean of the absolute values of all deltas in one list .
Parameters
- - - - - - - - - -
2021-11-17 16:44:54 +01:00
deltas : list
List of fluctuations
idx : list
List or ranges of configs on which the deltas are defined .
2021-11-07 21:44:22 +00:00
eps : float
Prefactor that enters the filter criterion .
2021-10-28 17:46:09 +02:00
"""
2021-11-17 16:44:54 +01:00
new_deltas = [ ]
new_idx = [ ]
maxd = np . mean ( np . fabs ( deltas ) )
for i in range ( len ( deltas ) ) :
if abs ( deltas [ i ] ) > eps * maxd :
new_deltas . append ( deltas [ i ] )
new_idx . append ( idx [ i ] )
if new_idx :
return np . array ( new_deltas ) , new_idx
2021-11-16 14:26:18 +01:00
else :
2021-11-17 16:44:54 +01:00
return deltas , idx
2021-10-28 17:46:09 +02:00
2021-01-08 14:46:37 +01:00
def derived_observable ( func , data , * * kwargs ) :
""" Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
Parameters
- - - - - - - - - -
2021-11-01 14:41:57 +00:00
func : object
arbitrary function of the form func ( data , * * kwargs ) . For the
automatic differentiation to work , all numpy functions have to have
the autograd wrapper ( use ' import autograd.numpy as anp ' ) .
data : list
list of Obs , e . g . [ obs1 , obs2 , obs3 ] .
num_grad : bool
if True , numerical derivatives are used instead of autograd
( default False ) . To control the numerical differentiation the
kwargs of numdifftools . step_generators . MaxStepGenerator
can be used .
man_grad : list
manually supply a list or an array which contains the jacobian
of func . Use cautiously , supplying the wrong derivative will
not be intercepted .
2021-01-08 14:46:37 +01:00
Notes
- - - - -
For simple mathematical operations it can be practical to use anonymous
functions . For the ratio of two observables one can e . g . use
new_obs = derived_observable ( lambda x : x [ 0 ] / x [ 1 ] , [ obs1 , obs2 ] )
"""
data = np . asarray ( data )
raveled_data = data . ravel ( )
2021-10-17 12:23:27 +01:00
# Workaround for matrix operations containing non Obs data
for i_data in raveled_data :
if isinstance ( i_data , Obs ) :
first_name = i_data . names [ 0 ]
first_shape = i_data . shape [ first_name ]
2021-11-05 14:41:22 +00:00
first_idl = i_data . idl [ first_name ]
2021-10-17 12:23:27 +01:00
break
for i in range ( len ( raveled_data ) ) :
if isinstance ( raveled_data [ i ] , ( int , float ) ) :
2021-11-05 14:41:22 +00:00
raveled_data [ i ] = Obs ( [ raveled_data [ i ] + np . zeros ( first_shape ) ] , [ first_name ] , idl = [ first_idl ] )
2021-10-17 12:23:27 +01:00
2021-11-29 12:15:27 +01:00
allcov = { }
for o in raveled_data :
2021-11-30 11:08:30 +01:00
for name in o . cov_names :
2021-11-29 12:15:27 +01:00
if name in allcov :
if not np . array_equal ( allcov [ name ] , o . covobs [ name ] . cov ) :
raise Exception ( ' Inconsistent covariance matrices for %s ! ' % ( name ) )
else :
allcov [ name ] = o . covobs [ name ] . cov
2021-01-08 14:46:37 +01:00
n_obs = len ( raveled_data )
new_names = sorted ( set ( [ y for x in [ o . names for o in raveled_data ] for y in x ] ) )
2021-11-17 17:19:24 +01:00
is_merged = { name : ( len ( list ( filter ( lambda o : o . is_merged . get ( name , False ) is True , raveled_data ) ) ) > 0 ) for name in new_names }
2021-10-29 16:42:01 +01:00
reweighted = len ( list ( filter ( lambda o : o . reweighted is True , raveled_data ) ) ) > 0
2021-10-28 17:46:09 +02:00
new_idl_d = { }
for name in new_names :
idl = [ ]
for i_data in raveled_data :
tmp = i_data . idl . get ( name )
2021-01-08 14:46:37 +01:00
if tmp is not None :
2021-10-28 17:46:09 +02:00
idl . append ( tmp )
2021-11-11 10:50:22 +00:00
new_idl_d [ name ] = _merge_idx ( idl )
2021-11-17 16:44:54 +01:00
if not is_merged [ name ] :
is_merged [ name ] = ( 1 != len ( set ( [ len ( idx ) for idx in [ * idl , new_idl_d [ name ] ] ] ) ) )
2021-10-28 17:46:09 +02:00
2021-10-18 10:00:15 +01:00
if data . ndim == 1 :
values = np . array ( [ o . value for o in data ] )
else :
values = np . vectorize ( lambda x : x . value ) ( data )
2021-01-08 14:46:37 +01:00
new_values = func ( values , * * kwargs )
multi = 0
if isinstance ( new_values , np . ndarray ) :
multi = 1
new_r_values = { }
for name in new_names :
tmp_values = np . zeros ( n_obs )
for i , item in enumerate ( raveled_data ) :
tmp = item . r_values . get ( name )
if tmp is None :
tmp = item . value
tmp_values [ i ] = tmp
if multi > 0 :
tmp_values = np . array ( tmp_values ) . reshape ( data . shape )
new_r_values [ name ] = func ( tmp_values , * * kwargs )
if ' man_grad ' in kwargs :
deriv = np . asarray ( kwargs . get ( ' man_grad ' ) )
if new_values . shape + data . shape != deriv . shape :
raise Exception ( ' Manual derivative does not have correct shape. ' )
elif kwargs . get ( ' num_grad ' ) is True :
if multi > 0 :
raise Exception ( ' Multi mode currently not supported for numerical derivative ' )
options = {
' base_step ' : 0.1 ,
' step_ratio ' : 2.5 ,
' num_steps ' : None ,
' step_nom ' : None ,
' offset ' : None ,
' num_extrap ' : None ,
' use_exact_steps ' : None ,
' check_num_steps ' : None ,
' scale ' : None }
for key in options . keys ( ) :
kwarg = kwargs . get ( key )
if kwarg is not None :
options [ key ] = kwarg
2021-10-11 12:22:58 +01:00
tmp_df = nd . Gradient ( func , order = 4 , * * { k : v for k , v in options . items ( ) if v is not None } ) ( values , * * kwargs )
2021-01-08 14:46:37 +01:00
if tmp_df . size == 1 :
deriv = np . array ( [ tmp_df . real ] )
else :
deriv = tmp_df . real
else :
deriv = jacobian ( func ) ( values , * * kwargs )
final_result = np . zeros ( new_values . shape , dtype = object )
for i_val , new_val in np . ndenumerate ( new_values ) :
new_deltas = { }
2021-11-29 12:15:27 +01:00
new_grad = { }
2021-01-08 14:46:37 +01:00
for j_obs , obs in np . ndenumerate ( data ) :
for name in obs . names :
2021-11-30 11:08:30 +01:00
if name in obs . cov_names :
2021-11-29 12:15:27 +01:00
if name in new_grad :
new_grad [ name ] + = deriv [ i_val + j_obs ] * obs . covobs [ name ] . grad
else :
new_grad [ name ] = deriv [ i_val + j_obs ] * obs . covobs [ name ] . grad
else :
new_deltas [ name ] = new_deltas . get ( name , 0 ) + deriv [ i_val + j_obs ] * _expand_deltas_for_merge ( obs . deltas [ name ] , obs . idl [ name ] , obs . shape [ name ] , new_idl_d [ name ] )
2021-11-30 13:32:50 +01:00
new_covobs = { name : Covobs ( 0 , allcov [ name ] , name , grad = new_grad [ name ] ) for name in new_grad }
2021-01-08 14:46:37 +01:00
new_samples = [ ]
2021-10-21 15:51:17 +01:00
new_means = [ ]
2021-10-28 17:46:09 +02:00
new_idl = [ ]
2021-11-29 12:15:27 +01:00
new_names_obs = [ ]
2021-11-17 16:44:54 +01:00
for name in new_names :
2021-11-29 12:15:27 +01:00
if name not in new_covobs :
if is_merged [ name ] :
filtered_deltas , filtered_idl_d = _filter_zeroes ( new_deltas [ name ] , new_idl_d [ name ] )
else :
filtered_deltas = new_deltas [ name ]
filtered_idl_d = new_idl_d [ name ]
new_samples . append ( filtered_deltas )
new_idl . append ( filtered_idl_d )
new_means . append ( new_r_values [ name ] [ i_val ] )
new_names_obs . append ( name )
final_result [ i_val ] = Obs ( new_samples , new_names_obs , means = new_means , idl = new_idl )
for name in new_covobs :
final_result [ i_val ] . names . append ( name )
final_result [ i_val ] . shape [ name ] = 1
final_result [ i_val ] . idl [ name ] = [ ]
final_result [ i_val ] . covobs = new_covobs
2021-10-21 16:15:02 +01:00
final_result [ i_val ] . _value = new_val
2021-10-28 17:46:09 +02:00
final_result [ i_val ] . is_merged = is_merged
final_result [ i_val ] . reweighted = reweighted
2021-01-08 14:46:37 +01:00
if multi == 0 :
final_result = final_result . item ( )
return final_result
2021-11-11 10:50:22 +00:00
def _reduce_deltas ( deltas , idx_old , idx_new ) :
2021-10-28 17:46:09 +02:00
""" Extract deltas defined on idx_old on all configs of idx_new.
Parameters
- - - - - - - - - -
2021-11-07 21:44:22 +00:00
deltas : list
List of fluctuations
idx_old : list
List or range of configs on which the deltas are defined
idx_new : list
List of configs for which we want to extract the deltas .
Has to be a subset of idx_old .
2021-10-28 17:46:09 +02:00
"""
if not len ( deltas ) == len ( idx_old ) :
2021-11-15 11:13:12 +00:00
raise Exception ( ' Length of deltas and idx_old have to be the same: %d != %d ' % ( len ( deltas ) , len ( idx_old ) ) )
2021-10-28 17:46:09 +02:00
if type ( idx_old ) is range and type ( idx_new ) is range :
if idx_old == idx_new :
return deltas
shape = len ( idx_new )
ret = np . zeros ( shape )
oldpos = 0
for i in range ( shape ) :
if oldpos == idx_old [ i ] :
raise Exception ( ' idx_old and idx_new do not match! ' )
pos = - 1
for j in range ( oldpos , len ( idx_old ) ) :
if idx_old [ j ] == idx_new [ i ] :
pos = j
break
if pos < 0 :
2021-11-11 10:50:22 +00:00
raise Exception ( ' Error in _reduce_deltas: Config %d not in idx_old ' % ( idx_new [ i ] ) )
2021-10-28 17:46:09 +02:00
ret [ i ] = deltas [ j ]
return np . array ( ret )
2021-01-08 14:46:37 +01:00
def reweight ( weight , obs , * * kwargs ) :
2021-10-28 17:46:09 +02:00
""" Reweight a list of observables.
Parameters
- - - - - - - - - -
2021-11-01 14:41:57 +00:00
weight : Obs
Reweighting factor . An Observable that has to be defined on a superset of the
configurations in obs [ i ] . idl for all i .
obs : list
list of Obs , e . g . [ obs1 , obs2 , obs3 ] .
2021-11-02 14:19:00 +00:00
all_configs : bool
if True , the reweighted observables are normalized by the average of
the reweighting factor on all configurations in weight . idl and not
on the configurations in obs [ i ] . idl .
2021-10-28 17:46:09 +02:00
"""
2021-01-08 14:46:37 +01:00
result = [ ]
for i in range ( len ( obs ) ) :
2021-11-30 11:08:30 +01:00
if len ( obs [ i ] . cov_names ) :
raise Exception ( ' Error: Not possible to reweight an Obs that contains covobs! ' )
2021-01-08 14:46:37 +01:00
if sorted ( weight . names ) != sorted ( obs [ i ] . names ) :
raise Exception ( ' Error: Ensembles do not fit ' )
for name in weight . names :
2021-10-28 17:46:09 +02:00
if not set ( obs [ i ] . idl [ name ] ) . issubset ( weight . idl [ name ] ) :
raise Exception ( ' obs[ %d ] has to be defined on a subset of the configs in weight.idl[ %s ]! ' % ( i , name ) )
2021-01-08 14:46:37 +01:00
new_samples = [ ]
2021-10-28 17:46:09 +02:00
w_deltas = { }
2021-01-08 14:46:37 +01:00
for name in sorted ( weight . names ) :
2021-11-11 10:50:22 +00:00
w_deltas [ name ] = _reduce_deltas ( weight . deltas [ name ] , weight . idl [ name ] , obs [ i ] . idl [ name ] )
2021-10-28 17:46:09 +02:00
new_samples . append ( ( w_deltas [ name ] + weight . r_values [ name ] ) * ( obs [ i ] . deltas [ name ] + obs [ i ] . r_values [ name ] ) )
tmp_obs = Obs ( new_samples , sorted ( weight . names ) , idl = [ obs [ i ] . idl [ name ] for name in sorted ( weight . names ) ] )
if kwargs . get ( ' all_configs ' ) :
new_weight = weight
else :
new_weight = Obs ( [ w_deltas [ name ] + weight . r_values [ name ] for name in sorted ( weight . names ) ] , sorted ( weight . names ) , idl = [ obs [ i ] . idl [ name ] for name in sorted ( weight . names ) ] )
2021-01-08 14:46:37 +01:00
2021-10-28 17:46:09 +02:00
result . append ( derived_observable ( lambda x , * * kwargs : x [ 0 ] / x [ 1 ] , [ tmp_obs , new_weight ] , * * kwargs ) )
2021-10-29 11:56:58 +01:00
result [ - 1 ] . reweighted = True
2021-10-28 17:46:09 +02:00
result [ - 1 ] . is_merged = obs [ i ] . is_merged
2021-01-08 14:46:37 +01:00
return result
def correlate ( obs_a , obs_b ) :
""" Correlate two observables.
2021-11-07 21:44:22 +00:00
Parameters
- - - - - - - - - -
2021-11-03 14:01:52 +00:00
obs_a : Obs
First observable
obs_b : Obs
Second observable
2021-01-08 14:46:37 +01:00
Keep in mind to only correlate primary observables which have not been reweighted
yet . The reweighting has to be applied after correlating the observables .
Currently only works if ensembles are identical . This is not really necessary .
"""
if sorted ( obs_a . names ) != sorted ( obs_b . names ) :
raise Exception ( ' Ensembles do not fit ' )
2021-11-30 11:08:30 +01:00
if len ( obs_a . cov_names ) or len ( obs_b . cov_names ) :
raise Exception ( ' Error: Not possible to correlate Obs that contain covobs! ' )
2021-01-08 14:46:37 +01:00
for name in obs_a . names :
if obs_a . shape [ name ] != obs_b . shape [ name ] :
raise Exception ( ' Shapes of ensemble ' , name , ' do not fit ' )
2021-11-05 12:19:23 +00:00
if obs_a . idl [ name ] != obs_b . idl [ name ] :
raise Exception ( ' idl of ensemble ' , name , ' do not fit ' )
2021-01-08 14:46:37 +01:00
2021-10-29 11:56:58 +01:00
if obs_a . reweighted is True :
2021-10-15 12:11:06 +01:00
warnings . warn ( " The first observable is already reweighted. " , RuntimeWarning )
2021-10-29 11:56:58 +01:00
if obs_b . reweighted is True :
2021-10-15 12:11:06 +01:00
warnings . warn ( " The second observable is already reweighted. " , RuntimeWarning )
2021-01-08 14:46:37 +01:00
new_samples = [ ]
2021-11-05 12:41:13 +00:00
new_idl = [ ]
2021-01-08 14:46:37 +01:00
for name in sorted ( obs_a . names ) :
new_samples . append ( ( obs_a . deltas [ name ] + obs_a . r_values [ name ] ) * ( obs_b . deltas [ name ] + obs_b . r_values [ name ] ) )
2021-11-05 12:41:13 +00:00
new_idl . append ( obs_a . idl [ name ] )
2021-01-08 14:46:37 +01:00
2021-11-05 12:41:13 +00:00
o = Obs ( new_samples , sorted ( obs_a . names ) , idl = new_idl )
2021-11-17 16:44:54 +01:00
o . is_merged = { name : ( obs_a . is_merged . get ( name , False ) or obs_b . is_merged . get ( name , False ) ) for name in o . names }
2021-10-29 11:56:58 +01:00
o . reweighted = obs_a . reweighted or obs_b . reweighted
2021-10-28 17:46:09 +02:00
return o
2021-01-08 14:46:37 +01:00
def covariance ( obs1 , obs2 , correlation = False , * * kwargs ) :
""" Calculates the covariance of two observables.
covariance ( obs , obs ) is equal to obs . dvalue * * 2
The gamma method has to be applied first to both observables .
If abs ( covariance ( obs1 , obs2 ) ) > obs1 . dvalue * obs2 . dvalue , the covariance
is constrained to the maximum value in order to make sure that covariance
matrices are positive semidefinite .
2021-11-07 21:44:22 +00:00
Parameters
- - - - - - - - - -
2021-11-12 11:59:40 +00:00
obs1 : Obs
First Obs
obs2 : Obs
Second Obs
2021-11-07 21:44:22 +00:00
correlation : bool
if true the correlation instead of the covariance is
returned ( default False )
2021-01-08 14:46:37 +01:00
"""
2021-11-15 16:39:50 +01:00
if set ( obs1 . names ) . isdisjoint ( set ( obs2 . names ) ) :
return 0.
2021-01-08 14:46:37 +01:00
for name in sorted ( set ( obs1 . names + obs2 . names ) ) :
2021-10-28 20:07:14 +02:00
if ( obs1 . shape . get ( name ) != obs2 . shape . get ( name ) ) and ( obs1 . shape . get ( name ) is not None ) and ( obs2 . shape . get ( name ) is not None ) :
raise Exception ( ' Shapes of ensemble ' , name , ' do not fit ' )
2021-11-11 10:50:22 +00:00
if ( 1 != len ( set ( [ len ( idx ) for idx in [ obs1 . idl [ name ] , obs2 . idl [ name ] , _merge_idx ( [ obs1 . idl [ name ] , obs2 . idl [ name ] ] ) ] ] ) ) ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' Shapes of ensemble ' , name , ' do not fit ' )
2021-10-26 15:43:30 +01:00
if not hasattr ( obs1 , ' e_names ' ) or not hasattr ( obs2 , ' e_names ' ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' The gamma method has to be applied to both Obs first. ' )
dvalue = 0
2021-11-30 11:08:30 +01:00
for e_name in obs1 . mc_names :
2021-01-08 14:46:37 +01:00
if e_name not in obs2 . e_names :
continue
gamma = 0
r_length = [ ]
for r_name in obs1 . e_content [ e_name ] :
if r_name not in obs2 . e_content [ e_name ] :
continue
r_length . append ( len ( obs1 . deltas [ r_name ] ) )
gamma + = np . sum ( obs1 . deltas [ r_name ] * obs2 . deltas [ r_name ] )
e_N = np . sum ( r_length )
tau_combined = ( obs1 . e_tauint [ e_name ] + obs2 . e_tauint [ e_name ] ) / 2
dvalue + = gamma / e_N * ( 1 + 1 / e_N ) / e_N * 2 * tau_combined
2021-11-30 11:08:30 +01:00
for e_name in obs1 . cov_names :
if e_name not in obs2 . cov_names :
continue
dvalue + = float ( np . dot ( np . transpose ( obs1 . covobs [ e_name ] . grad ) , np . dot ( obs1 . covobs [ e_name ] . cov , obs2 . covobs [ e_name ] . grad ) ) )
2021-01-08 14:46:37 +01:00
if np . abs ( dvalue / obs1 . dvalue / obs2 . dvalue ) > 1.0 :
dvalue = np . sign ( dvalue ) * obs1 . dvalue * obs2 . dvalue
if correlation :
dvalue = dvalue / obs1 . dvalue / obs2 . dvalue
return dvalue
def covariance2 ( obs1 , obs2 , correlation = False , * * kwargs ) :
""" Alternative implementation of the covariance of two observables.
covariance ( obs , obs ) is equal to obs . dvalue * * 2
The gamma method has to be applied first to both observables .
If abs ( covariance ( obs1 , obs2 ) ) > obs1 . dvalue * obs2 . dvalue , the covariance
is constrained to the maximum value in order to make sure that covariance
matrices are positive semidefinite .
Keyword arguments
- - - - - - - - - - - - - - - - -
correlation - - if true the correlation instead of the covariance is
returned ( default False )
"""
2021-10-28 17:46:09 +02:00
def expand_deltas ( deltas , idx , shape , new_idx ) :
""" Expand deltas defined on idx to a contiguous range [new_idx[0], new_idx[-1]].
New , empy entries are filled by 0. If idx and new_idx are of type range , the smallest
common divisor of the step sizes is used as new step size .
Parameters
- - - - - - - - - -
deltas - - List of fluctuations
idx - - List or range of configs on which the deltas are defined .
Has to be a subset of new_idx .
shape - - Number of configs in idx .
new_idx - - List of configs that defines the new range .
"""
if type ( idx ) is range and type ( new_idx ) is range :
if idx == new_idx :
return deltas
ret = np . zeros ( new_idx [ - 1 ] - new_idx [ 0 ] + 1 )
for i in range ( shape ) :
ret [ idx [ i ] - new_idx [ 0 ] ] = deltas [ i ]
return ret
def calc_gamma ( deltas1 , deltas2 , idx1 , idx2 , new_idx , w_max ) :
gamma = np . zeros ( w_max )
deltas1 = expand_deltas ( deltas1 , idx1 , len ( idx1 ) , new_idx )
deltas2 = expand_deltas ( deltas2 , idx2 , len ( idx2 ) , new_idx )
new_shape = len ( deltas1 )
max_gamma = min ( new_shape , w_max )
# The padding for the fft has to be even
padding = new_shape + max_gamma + ( new_shape + max_gamma ) % 2
gamma [ : max_gamma ] + = ( np . fft . irfft ( np . fft . rfft ( deltas1 , padding ) * np . conjugate ( np . fft . rfft ( deltas2 , padding ) ) ) [ : max_gamma ] + np . fft . irfft ( np . fft . rfft ( deltas2 , padding ) * np . conjugate ( np . fft . rfft ( deltas1 , padding ) ) ) [ : max_gamma ] ) / 2.0
return gamma
2021-01-08 14:46:37 +01:00
2021-11-15 16:39:50 +01:00
if set ( obs1 . names ) . isdisjoint ( set ( obs2 . names ) ) :
return 0.
2021-10-26 15:43:30 +01:00
if not hasattr ( obs1 , ' e_names ' ) or not hasattr ( obs2 , ' e_names ' ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' The gamma method has to be applied to both Obs first. ' )
dvalue = 0
e_gamma = { }
e_dvalue = { }
e_n_tauint = { }
e_rho = { }
2021-11-30 11:08:30 +01:00
for e_name in obs1 . mc_names :
2021-01-08 14:46:37 +01:00
2021-11-30 11:08:30 +01:00
if e_name not in obs2 . mc_names :
2021-01-08 14:46:37 +01:00
continue
2021-10-28 17:46:09 +02:00
idl_d = { }
2021-01-08 14:46:37 +01:00
r_length = [ ]
for r_name in obs1 . e_content [ e_name ] :
2021-10-28 18:07:09 +02:00
if r_name not in obs2 . e_content [ e_name ] :
continue
2021-11-11 10:50:22 +00:00
idl_d [ r_name ] = _merge_idx ( [ obs1 . idl [ r_name ] , obs2 . idl [ r_name ] ] )
2021-11-12 09:41:17 +00:00
if isinstance ( idl_d [ r_name ] , range ) :
2021-10-28 17:46:09 +02:00
r_length . append ( len ( idl_d [ r_name ] ) )
else :
r_length . append ( ( idl_d [ r_name ] [ - 1 ] - idl_d [ r_name ] [ 0 ] + 1 ) )
2021-01-08 14:46:37 +01:00
2021-10-28 18:07:09 +02:00
if not r_length :
return 0.
2021-10-28 19:53:15 +02:00
2021-01-08 14:46:37 +01:00
w_max = max ( r_length ) / / 2
e_gamma [ e_name ] = np . zeros ( w_max )
for r_name in obs1 . e_content [ e_name ] :
if r_name not in obs2 . e_content [ e_name ] :
continue
2021-10-28 17:46:09 +02:00
e_gamma [ e_name ] + = calc_gamma ( obs1 . deltas [ r_name ] , obs2 . deltas [ r_name ] , obs1 . idl [ r_name ] , obs2 . idl [ r_name ] , idl_d [ r_name ] , w_max )
2021-01-08 14:46:37 +01:00
2021-10-28 17:46:09 +02:00
if np . all ( e_gamma [ e_name ] == 0.0 ) :
2021-01-08 14:46:37 +01:00
continue
e_shapes = [ ]
for r_name in obs1 . e_content [ e_name ] :
e_shapes . append ( obs1 . shape [ r_name ] )
2021-10-28 17:46:09 +02:00
gamma_div = np . zeros ( w_max )
e_N = 0
for r_name in obs1 . e_content [ e_name ] :
if r_name not in obs2 . e_content [ e_name ] :
continue
gamma_div + = calc_gamma ( np . ones ( obs1 . shape [ r_name ] ) , np . ones ( obs2 . shape [ r_name ] ) , obs1 . idl [ r_name ] , obs2 . idl [ r_name ] , idl_d [ r_name ] , w_max )
e_N + = np . sum ( np . ones_like ( idl_d [ r_name ] ) )
e_gamma [ e_name ] / = gamma_div [ : w_max ]
2021-01-08 14:46:37 +01:00
e_rho [ e_name ] = e_gamma [ e_name ] [ : w_max ] / e_gamma [ e_name ] [ 0 ]
e_n_tauint [ e_name ] = np . cumsum ( np . concatenate ( ( [ 0.5 ] , e_rho [ e_name ] [ 1 : ] ) ) )
# Make sure no entry of tauint is smaller than 0.5
e_n_tauint [ e_name ] [ e_n_tauint [ e_name ] < 0.5 ] = 0.500000000001
window = max ( obs1 . e_windowsize [ e_name ] , obs2 . e_windowsize [ e_name ] )
# Bias correction hep-lat/0306017 eq. (49)
e_dvalue [ e_name ] = 2 * ( e_n_tauint [ e_name ] [ window ] + obs1 . tau_exp [ e_name ] * np . abs ( e_rho [ e_name ] [ window + 1 ] ) ) * ( 1 + ( 2 * window + 1 ) / e_N ) * e_gamma [ e_name ] [ 0 ] / e_N
dvalue + = e_dvalue [ e_name ]
2021-11-30 11:08:30 +01:00
for e_name in obs1 . cov_names :
if e_name not in obs2 . cov_names :
continue
dvalue + = float ( np . dot ( np . transpose ( obs1 . covobs [ e_name ] . grad ) , np . dot ( obs1 . covobs [ e_name ] . cov , obs2 . covobs [ e_name ] . grad ) ) )
2021-01-08 14:46:37 +01:00
if np . abs ( dvalue / obs1 . dvalue / obs2 . dvalue ) > 1.0 :
dvalue = np . sign ( dvalue ) * obs1 . dvalue * obs2 . dvalue
if correlation :
dvalue = dvalue / obs1 . dvalue / obs2 . dvalue
return dvalue
def covariance3 ( obs1 , obs2 , correlation = False , * * kwargs ) :
""" Another alternative implementation of the covariance of two observables.
covariance2 ( obs , obs ) is equal to obs . dvalue * * 2
Currently only works if ensembles are identical .
The gamma method has to be applied first to both observables .
If abs ( covariance2 ( obs1 , obs2 ) ) > obs1 . dvalue * obs2 . dvalue , the covariance
is constrained to the maximum value in order to make sure that covariance
matrices are positive semidefinite .
Keyword arguments
- - - - - - - - - - - - - - - - -
correlation - - if true the correlation instead of the covariance is
returned ( default False )
plot - - if true , the integrated autocorrelation time for each ensemble is
plotted .
"""
for name in sorted ( set ( obs1 . names + obs2 . names ) ) :
2021-10-28 20:07:14 +02:00
if ( obs1 . shape . get ( name ) != obs2 . shape . get ( name ) ) and ( obs1 . shape . get ( name ) is not None ) and ( obs2 . shape . get ( name ) is not None ) :
raise Exception ( ' Shapes of ensemble ' , name , ' do not fit ' )
2021-11-11 10:50:22 +00:00
if ( 1 != len ( set ( [ len ( idx ) for idx in [ obs1 . idl [ name ] , obs2 . idl [ name ] , _merge_idx ( [ obs1 . idl [ name ] , obs2 . idl [ name ] ] ) ] ] ) ) ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' Shapes of ensemble ' , name , ' do not fit ' )
2021-10-26 15:43:30 +01:00
if not hasattr ( obs1 , ' e_names ' ) or not hasattr ( obs2 , ' e_names ' ) :
2021-01-08 14:46:37 +01:00
raise Exception ( ' The gamma method has to be applied to both Obs first. ' )
tau_exp = [ ]
S = [ ]
for e_name in sorted ( set ( obs1 . e_names + obs2 . e_names ) ) :
t_1 = obs1 . tau_exp . get ( e_name )
t_2 = obs2 . tau_exp . get ( e_name )
if t_1 is None :
t_1 = 0
if t_2 is None :
t_2 = 0
tau_exp . append ( max ( t_1 , t_2 ) )
S_1 = obs1 . S . get ( e_name )
S_2 = obs2 . S . get ( e_name )
if S_1 is None :
S_1 = Obs . S_global
if S_2 is None :
S_2 = Obs . S_global
S . append ( max ( S_1 , S_2 ) )
check_obs = obs1 + obs2
check_obs . gamma_method ( tau_exp = tau_exp , S = S )
if kwargs . get ( ' plot ' ) :
check_obs . plot_tauint ( )
check_obs . plot_rho ( )
cov = ( check_obs . dvalue * * 2 - obs1 . dvalue * * 2 - obs2 . dvalue * * 2 ) / 2
if np . abs ( cov / obs1 . dvalue / obs2 . dvalue ) > 1.0 :
cov = np . sign ( cov ) * obs1 . dvalue * obs2 . dvalue
if correlation :
cov = cov / obs1 . dvalue / obs2 . dvalue
return cov
def pseudo_Obs ( value , dvalue , name , samples = 1000 ) :
""" Generate a pseudo Obs with given value, dvalue and name
2021-11-12 11:59:40 +00:00
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 ) .
2021-01-08 14:46:37 +01:00
"""
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
2021-10-21 16:15:02 +01:00
res . _value = float ( value )
2021-01-08 14:46:37 +01:00
return res
2021-11-18 11:17:20 +00:00
def import_jackknife ( jacks , name , idl = None ) :
2021-11-17 16:06:26 +00:00
""" Imports jackknife samples and returns an Obs
Parameters
- - - - - - - - - -
jacks : numpy . ndarray
numpy array containing the mean value as zeroth entry and
the N jackknife samples as first to Nth entry .
name : str
name of the ensemble the samples are defined on .
"""
length = len ( jacks ) - 1
prj = ( np . ones ( ( length , length ) ) - ( length - 1 ) * np . identity ( length ) )
samples = jacks [ 1 : ] @ prj
2021-11-18 11:17:20 +00:00
new_obs = Obs ( [ samples ] , [ name ] , idl = idl )
2021-11-17 16:06:26 +00:00
new_obs . _value = jacks [ 0 ]
return new_obs
2021-01-08 14:46:37 +01:00
def merge_obs ( list_of_obs ) :
""" Combine all observables in list_of_obs into one new observable
2021-11-07 21:44:22 +00:00
Parameters
- - - - - - - - - -
list_of_obs : list
list of the Obs object to be combined
2021-01-08 14:46:37 +01:00
It is not possible to combine obs which are based on the same replicum
"""
replist = [ item for obs in list_of_obs for item in obs . names ]
if ( len ( replist ) == len ( set ( replist ) ) ) is False :
2021-10-11 12:22:58 +01:00
raise Exception ( ' list_of_obs contains duplicate replica: %s ' % ( str ( replist ) ) )
2021-11-30 11:08:30 +01:00
if any ( [ len ( o . cov_names ) for o in list_of_obs ] ) :
raise Exception ( ' Not possible to merge data that contains covobs! ' )
2021-01-08 14:46:37 +01:00
new_dict = { }
2021-10-28 17:46:09 +02:00
idl_dict = { }
2021-01-08 14:46:37 +01:00
for o in list_of_obs :
new_dict . update ( { key : o . deltas . get ( key , 0 ) + o . r_values . get ( key , 0 )
2021-10-11 12:22:58 +01:00
for key in set ( o . deltas ) | set ( o . r_values ) } )
2021-10-28 17:46:09 +02:00
idl_dict . update ( { key : o . idl . get ( key , 0 ) for key in set ( o . deltas ) } )
2021-01-08 14:46:37 +01:00
2021-10-28 17:46:09 +02:00
names = sorted ( new_dict . keys ( ) )
o = Obs ( [ new_dict [ name ] for name in names ] , names , idl = [ idl_dict [ name ] for name in names ] )
2021-11-17 16:44:54 +01:00
o . is_merged = { name : np . any ( [ oi . is_merged . get ( name , False ) for oi in list_of_obs ] ) for name in o . names }
2021-10-28 17:46:09 +02:00
o . reweighted = np . max ( [ oi . reweighted for oi in list_of_obs ] )
return o
2021-11-29 12:15:27 +01:00
2021-11-30 11:08:30 +01:00
def cov_Obs ( means , cov , name , grad = None ) :
""" Create an Obs based on mean(s) and a covariance matrix
2021-11-29 12:15:27 +01:00
Parameters
- - - - - - - - - -
2021-11-30 11:08:30 +01:00
mean : list of floats or float
N mean value ( s ) of the new Obs
2021-11-29 12:15:27 +01:00
cov : list or array
2021-11-30 11:08:30 +01:00
2 d ( NxN ) Covariance matrix , 1 d diagonal entries or 0 d covariance
2021-11-29 12:15:27 +01:00
name : str
identifier for the covariance matrix
grad : list or array
Gradient of the Covobs wrt . the means belonging to cov .
"""
2021-11-30 11:08:30 +01:00
def covobs_to_obs ( co ) :
""" Make an Obs out of a Covobs
2021-11-29 12:15:27 +01:00
2021-11-30 11:08:30 +01:00
Parameters
- - - - - - - - - -
co : Covobs
Covobs to be embedded into the Obs
"""
o = Obs ( None , None , empty = True )
o . _value = co . value
o . names . append ( co . name )
o . covobs [ co . name ] = co
o . _dvalue = np . sqrt ( co . errsq ( ) )
o . shape [ co . name ] = 1
o . idl [ co . name ] = [ ]
return o
2021-11-29 12:15:27 +01:00
ol = [ ]
2021-11-30 11:08:30 +01:00
if isinstance ( means , ( float , int ) ) :
means = [ means ]
2021-11-29 12:15:27 +01:00
for i in range ( len ( means ) ) :
ol . append ( covobs_to_obs ( Covobs ( means [ i ] , cov , name , pos = i , grad = grad ) ) )
if ol [ 0 ] . covobs [ name ] . N != len ( means ) :
raise Exception ( ' You have to provide %d mean values! ' % ( ol [ 0 ] . N ) )
2021-11-30 11:08:30 +01:00
if len ( ol ) == 1 :
return ol [ 0 ]
2021-11-29 12:15:27 +01:00
return ol