Commit 0e1d0ad2 authored by Olaf's avatar Olaf

fix error estimation energy plot

parent 167e55db
# import
import numpy as np
from scipy import optimize
def chi(A,x):
"""
function to determine chi_A(t)
last modified: 20/03/19
Parameters
----------
A : (1,n_max)-array
input data
x : (1,n_max-1)-array
x data
Returns
-------
chi : (1,x)-array
autocorrelation function
"""
# Process input
n_max = np.size(A)
A_tile = np.tile(A,[n_max,1]) # (n_max,n_max)-array
# Calculate scalefactor N
N = 1/(n_max-x) # (1,x)-array
# Calculate chi
# # chi_1(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n)A(n+t)
AA = np.convolve(A,np.flip(A,0))[0:n_max][::-1]
chi_1 = np.multiply(N,AA) # (1,n_max)-array
# # chi_2(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n)
A_triu_2 = np.triu(np.flip(A_tile,1))
A_2 = np.sum(A_triu_2,1)
chi_2 = np.multiply(N,A_2) # (1,n_max)-array
# # chi_3(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n+1)
A_triu_3 = np.triu(A_tile)
A_3 = np.sum(A_triu_3,1)
chi_3 = np.multiply(N,A_3) # (1,n_max)-array
# # chi = chi_1 + chi_2*chi_3
chi = chi_1 - np.multiply(chi_2,chi_3) # (1,n_max)-array
# return
return chi
def test_func(x, a, b):
"""
exponential function to fit chi
last modified: 26/03/2019
Parameters
----------
x : (1,N) array
x_data
a : (float)
amplitude
b : (float)
tau
Returns
-------
a*np.exp(-x/b) : (1,N) array
function
"""
return a*np.exp(-x/b)
def tau(y_data):
"""
function to find correlation time of the simulation (tau) by fitting chi to exponential test function
last modified: 26/03/2019
Parameters
----------
y_data : (1,N) array
data to which function is fitted
Returns
-------
tau : (float)
correlation time of the simulation
"""
# data
x_data = np.array(range(len(y_data)))
# fit and determine tau
[amplitude, tau], params_covariance = optimize.curve_fit(test_func, x_data,y_data)
# return
return tau
def calc_std_mean(tau,data):
"""
function to calculate standard deviation (std) and mean (mean) of data using data blocking
last modified: 20/03/19
Parameters
----------
tau : float
correlation time of the simulation
data : (1,N)-array
data from which std and mean are calculated
Returns
-------
data_std : float
standard deviation of data
data_mean : float
mean of data
"""
# define length of data block
m = 2*np.abs(np.ceil(tau))
# select multiple of m datasets of data
data_select = data[:int(len(data)-len(data)%m)]
# reshape data to (n,m)-array
n = len(data_select)/m
data_reshape= data_select.reshape(int(n),int(m))
# average in data_blocks
data_j = np.mean(data_reshape,1)
# determine std and mean from average data_blocks
data_std = np.std(data_j)
data_mean= np.mean(data_j)
# return
return data_std, data_mean
def error_combine(data):
"""
overhead function to calculate standard deviation and mean using datablocking
last modified: 20/03/2019
Parameters
----------
data : (1,N)-array
data from which std and mean are calculated
Returns
-------
data_std : (float)
standard deviation of data
data_mean : (float)
mean of data
"""
# process input
n_max = np.size(data)
x_data= np.linspace(0,n_max-1,n_max)
# calculate chi
chi_output = chi(data,x_data)
# calculate tau
tau_output = tau(chi_output)
# calculate std and mean
[data_std, data_mean] = calc_std_mean(tau_output,data)
# return
return data_std, data_mean
\ No newline at end of file
......@@ -137,9 +137,9 @@ def plot_energy(bodies, settings):
E_kin += np.asarray(bodies[i].E_kin)
E_pot += np.asarray(bodies[i].E_pot)
E_kin_std, E_kin_mean = error_combine(E_kin)
E_pot_std, E_pot_mean = error_combine(E_pot)
E_std,E_mean = error_combine(E_kin + E_pot)
E_kin_std = np.std(E_kin)
E_pot_std = np.std(E_pot)
E_std,E_mean = np.std(E_kin + E_pot)
fig = plt.figure()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment