Commit 12f945e4 authored by Luuk Balkenende's avatar Luuk Balkenende

Adjusted radius of gyration to be fit for bootstrap function in processfunctions.py

parent 247457d7
......@@ -5,16 +5,16 @@ from processfunctions import exp_calc, bootstrap
# Simulation parameters
n_beads = 250 # number of beads
n = 6 # number of angles
n_pol = 10000 # number of polymers
n_pol = 1000 # number of polymers
alfa_up = 2 # alfa to compute upper limit
alfa_low = 0.7 # alfa to compute lower limit
alfa_low = 1.2 # alfa to compute lower limit
# Choose algorithm
PERM = False # Put on True if PERM algorithm is desired
PERM = False # Put on True if PERM algorithm is desired
random_walk = False # Put on True if random walk is desired
# Physical parameters
T_dim = 10 # T * (Kb / eps) # Non-dimensional temperature
T_dim = 0.2 # T * (Kb / eps) # Non-dimensional temperature
# Simulate configuration of polymers and their weights
end_end_array, polweight_array, pos_array = calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk)
......
......@@ -22,39 +22,39 @@ def exp_calc(polweight_array, end_end_array, n_beads, pos_array):
"""
exp_value_mean, exp_value_std = bootstrap(polweight_array, end_end_array, n_beads)
rad_gyr = radiusgyration(pos_array)
rad_gyr_mean, rad_gyr_std = bootstrap(polweight_array, rad_gyr, n_beads)
popsize = polsize(end_end_array)
a = np.arange(0, n_beads - 1, 1)
b = np.arange(0, n_beads, 1)
exp_value_mean[0:2] = [0, 1]
rad_gyr_mean[0] = 0
font = {'family' : 'DejaVu Sans',
'weight' : 'normal',
'size' : 17}
font = {'family': 'DejaVu Sans',
'weight': 'normal',
'size': 17}
matplotlib.rc('font', **font)
q = np.polyfit(np.log(a[1::]), np.log(exp_value_mean[1::]**2), 1)
q = np.polyfit(np.log(a[1::]), np.log(exp_value_mean[1::] ** 2), 1)
print(q)
xx = np.linspace(0, 250 , 1000)
yy = np.exp(q[1])*xx**q[0]
xx = np.linspace(0, 250, 1000)
yy = np.exp(q[1]) * xx ** q[0]
plt.figure(1)
plt.xscale("log")
plt.yscale("log")
plt.plot( xx, yy, 'r--', label="Fit", linewidth =3)
plt.plot(xx, yy, 'r--', label="Fit", linewidth=3)
plt.xlabel(r'$N_{beads}$')
plt.ylabel(r'$R^2$ $[\sigma^2]$')
plt.xscale("log")
plt.yscale("log")
plt.errorbar(a, exp_value_mean[0:n_beads - 1] ** 2, fmt='x', xerr=None, yerr=2*exp_value_std*exp_value_mean, label='Data',color="k", capsize=3, capthick=1, markersize=3)
plt.errorbar(a, exp_value_mean[0:n_beads - 1] ** 2, fmt='x', xerr=None, yerr=2 * exp_value_std * exp_value_mean,
label='Data', color="k", capsize=3, capthick=1, markersize=3)
plt.grid(True, which="both", ls=":")
plt.legend(loc='best')
plt.xlim(1, 250)
......@@ -65,34 +65,32 @@ def exp_calc(polweight_array, end_end_array, n_beads, pos_array):
###########################################################################
q1 = np.polyfit(np.log(a[1::]), np.log(rad_gyr_mean[1::]), 1)
print(q1)
xx1 = np.linspace(0, 250 , 1000)
yy1 = np.exp(q1[1])*xx**q1[0]
xx1 = np.linspace(0, 250, 1000)
yy1 = np.exp(q1[1]) * xx ** q1[0]
plt.figure(2)
plt.xscale("log")
plt.yscale("log")
plt.plot( xx1, yy1, 'r--', label="Fit", linewidth =3)
plt.plot(xx1, yy1, 'r--', label="Fit", linewidth=3)
plt.xlabel(r'$N_{beads}$')
plt.ylabel(r'$R_g^2$ $[\sigma^2]$')
plt.xscale("log")
plt.yscale("log")
plt.errorbar(a, rad_gyr_mean, fmt='x', xerr=None, yerr=rad_gyr_std, label='Data',color="k", capsize=3, capthick=1, markersize=3)
plt.errorbar(a, rad_gyr_mean, fmt='x', xerr=None, yerr=rad_gyr_std, label='Data', color="k", capsize=3, capthick=1,
markersize=3)
plt.grid(True, which="both", ls=":")
plt.legend(loc='best')
plt.xlim(1, 250)
plt.ylim(bottom=0.1)
plt.ylim(top=1000)
plt.savefig('gyros', bbox_inches='tight', dpi=300)
plt.figure(3)
plt.title('Population size')
#plt.yscale("log")
# plt.yscale("log")
plt.plot(b[2:], popsize[2:])
plt.grid(True, which="both", ls="-")
plt.ylabel('Population size')
......@@ -100,6 +98,7 @@ def exp_calc(polweight_array, end_end_array, n_beads, pos_array):
plt.show()
def bootstrap(polweight_array, end_end_array, n_beads):
"""Preforms bootstrap onto end-to-end distances of polymers to calculate the error
......@@ -119,76 +118,92 @@ def bootstrap(polweight_array, end_end_array, n_beads):
exp_value_std: array of size (n_beads, 1)
Standard deviation of end-to-end distances per number of beads
"""
# HIER EIGENLIJK DE NAAM NOG VERANDEREN VAN END TO END NAAR OBSERVABLE (WANT WE GEBRUIKEN OOK GYROS HIERO)
# Initialize/allocate values
nn = 100 # Amount of bootstraps (n)
a = end_end_array.shape
nn = 100
exp_value = np.zeros((a[0]-1, nn))
exp_value = np.zeros((a[0] - 1, nn)) # Allocate observable
j = 0
m = 0
while j < nn:
# Make random array of size with values ranging from 0 to max(index variable)
c = np.random.randint(a[1], size=(a[1]))
# Get random values from variables using random array c
end_end_bootstrap = end_end_array[:, c]
polweight_bootstrap = polweight_array[:, c]
# Calculate true observable values
end_end_times_bootstrap = np.multiply(end_end_bootstrap, polweight_bootstrap)
sum_end_end_times_bootstrap = np.sum(end_end_times_bootstrap, axis=1)
sum_weights_bootstrap = np.sum(polweight_bootstrap, axis=1)
# If statements which warns the user when data is not sufficient
if np.count_nonzero(sum_weights_bootstrap) == n_beads:
exp_value[:, j] = np.divide(sum_end_end_times_bootstrap[0:n_beads - 1], sum_weights_bootstrap[0:n_beads - 1])
exp_value[:, j] = np.divide(sum_end_end_times_bootstrap[0:n_beads - 1],
sum_weights_bootstrap[0:n_beads - 1])
j = j + 1
elif np.count_nonzero(end_end_array[n_beads - 1, :]) == 0:
print('No polymer reached ', n_beads, ' beads, try again, or lower variable n_beads.')
exit()
elif m > 100:
print('Not enough polymers reached ', n_beads, ' beads, try again, or adjust variable n_beads.')
exit()
else:
j = j
m = m + 1
print('fail')
# print(sum_weights_bootstrap)
# Calculate mean and standard deviation of observable
exp_value_mean = np.mean(exp_value, axis=1)
exp_value_std = np.std(exp_value, axis=1)
return exp_value_mean, exp_value_std
def polsize(end_end_array):
# HIER NOG EEN DOCSTRING DINGENS!!
popsize = np.count_nonzero(end_end_array, axis=1)
print(popsize)
print(popsize.shape)
return popsize
def radiusgyration(pos_array):
# HIER NOG EEN DOCSTRING DINGENS!
# Preparation for calculation can start (allocation, initialization)
pos_array = np.delete(pos_array, [0, 1], axis=1)
#pos_array = np.delete(pos_array, 0, axis=0)
[nx, ny] = pos_array.shape
print('nx',nx)
#print('ny',ny)
#print('pos_array', pos_array)
imax = nx
jmax = ny
rad_gyr = np.zeros([imax, int(jmax/2)])
[imax, jmax] = pos_array.shape
rad_gyr = np.zeros([imax, int(jmax / 2)])
j = 0
# Loop over all polymers
while j < jmax:
i = 1
imax = np.count_nonzero(pos_array[:, j], axis=0)
#print('imax = ', imax)
#print('j = ', j)
# Loop over all beads
while i < imax:
# Calculate radius of gyration
N = i + 1
xk = pos_array[:i+1, j]
#print(xk)
#exit()
yk = pos_array[:i+1, j+1]
xk = pos_array[:i + 1, j]
yk = pos_array[:i + 1, j + 1]
xmean = np.sum(xk)/N
ymean = np.sum(yk)/N
xmean = np.sum(xk) / N
ymean = np.sum(yk) / N
xr = xk - xmean
yr = yk - ymean
rad_gyr[i, int(j/2)] = 1/N * np.sum(np.power(xr, 2) + np.power(yr, 2))
#print('rad_gyr ij', rad_gyr[i, int(j/2)])
#print('i = ', i, ' j = ', j)
i = i + 1
rad_gyr[i, int(j / 2)] = 1 / N * np.sum(np.power(xr, 2) + np.power(yr, 2))
i = i + 1
j = j + 2
print('shape gyros = ', rad_gyr.shape)
#rad_gyr = np.insert(rad_gyr, 0, rad_gyr[0, :], axis=0)
# Adjust array to get it fit for bootstrap function
rad_gyr_zeros = np.zeros(rad_gyr.shape[0])
rad_gyr = np.insert(rad_gyr, 0, rad_gyr_zeros, axis=1)
print('shape gyros = ', rad_gyr.shape)
return rad_gyr
\ No newline at end of file
This diff is collapsed.
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