Commit 1bc4890e authored by Luuk Balkenende's avatar Luuk Balkenende

Fixed divide by zero error in bootstrap

parent 6bcfa381
import numpy as np
import random
def add_bead(pos, n, T_dim, polweight, fails, cst_dep):
"""Adds one bead to polymer, according to the Rosenblaum method
Parameters:
----------
pos: array of size (n_beads, 2)
Position of n_beads beads in x,y direction
n: integer > 0
Number of possible angles
T_dim: float > 0
Non-dimensionalized temperature
j: integer > 1
Bead j is the j-th bead which is added to polymer
p_pol: GAAT WAARSCHIJNLIJK TOCH NOG VERANDEREN
Results:
--------
pos: array of size (n_beads, 2)
Position of n_beads beads in x,y direction
p_pol: GAAT WAARSCHIJNLIJK TOCH NOG VERANDEREN
"""
#compute first nonzero element of posstack to determine j
itemindex = np.where(pos[1::,0]==0)
j = int(itemindex[0][0]) + 1 # since first bead is not included
polweight = polweight[j-1]
# Standard phi, KAN NAAR MAIN.PY
phi_0 = np.linspace(0, 2 * np.pi - (2 / n) * np.pi, n)
phi = phi_0 + (1 - 2 * np.random.rand(n)) * np.pi / n
# Calculate new positions
pos_new = np.zeros([n, 2])
pos_new[:, 0] = pos[j - 1, 0] + np.cos(phi)
pos_new[:, 1] = pos[j - 1, 1] + np.sin(phi)
# Calculate distance between new bead and beads
rx = np.zeros([n, j])
ry = np.zeros([n, j])
for i in range(n):
rx[i, :] = pos_new[i, 0] - pos[0:j, 0]
ry[i, :] = pos_new[i, 1] - pos[0:j, 1]
r = np.sqrt(rx ** 2 + ry ** 2)
r[r == 0] = np.nan
E = np.sum(4 * 0.25*((0.8 / r) ** 12 - (0.8 / r) ** 6), axis=1)
#print(L)
w_jl = np.exp(-E / T_dim) #small weights (per angle)
w_l = np.sum(w_jl) #big weight per added bead (sum of angle weights)
if w_l == 0:
w_jl = np.ones(w_jl.shape)
w_l = np.sum(w_jl)
fails = fails+1
#print(w_l)
#print(w_l)
polweight = polweight*w_l *cst_dep #weight per polymer off by our algorithm
# Calculate new position of next bead
P = w_jl / w_l
final_phi = np.random.choice(phi, p=P)
pos[j] = np.array([[pos[j - 1, 0] + np.cos(final_phi), pos[j - 1, 1] + np.sin(final_phi)]])
end_end = np.sqrt(pos[j, 0]**2 + pos[j,1]**2)
return pos, polweight, end_end, fails, j
def enrich(polweight_stack, pos_stack, end_end_stack, j):
new_weight = 0.5 * polweight_stack[:, 0]
polweight_stack[:, 0] = new_weight
polweight_stack = np.c_[polweight_stack, polweight_stack[:,0]]
#polweight_stack = np.insert(polweight_stack, 1, new_weight, axis =1)
end_end_stack = np.c_[end_end_stack, end_end_stack[:,0]]
#end_end_stack = np.insert(end_end_stack, 1, end_end_stack[:,0], axis=1 )
pos_stack = np.c_[pos_stack, pos_stack[:,0]]
pos_stack = np.c_[pos_stack, pos_stack[:,1]]
#pos_stack = np.insert(pos_stack, 2, pos_stack[:,1], axis=1)
#pos_stack = np.insert(pos_stack, 2, pos_stack[:,0], axis=1)
return polweight_stack, pos_stack, end_end_stack
def limmits(polweight_matrix, polweight_stack, alfa_up, j, alfa_low):
# find all nonzero elements of stack & pw_matrix, of bead j, sum and average:
b = polweight_matrix[j,np.nonzero(polweight_matrix[j,:])]
c = polweight_stack[j,np.nonzero(polweight_stack[j,:])]
avg_weight = np.average(np.append(b,c))
upper = (alfa_up * avg_weight)/polweight_stack[2, 0]
lower = (alfa_low * avg_weight)/polweight_stack[2, 0]
return upper, lower
def prune(polweight_stack, end_end_stack, pos_stack, j ):
b = random.random()
if b < 0.5:
polweight_stack[:,0] = polweight_stack[:,0] * 2
else:
polweight_stack = np.delete(polweight_stack, (0), axis = 1)
end_end_stack = np.delete(end_end_stack, (0), axis = 1)
pos_stack = np.delete(pos_stack, (0), axis = 1)
pos_stack = np.delete(pos_stack, (0), axis = 1)
return polweight_stack, end_end_stack, pos_stack
def maxsize(polweight_stack, end_end_stack, pos_stack):
a = polweight_stack[:, 0]
a = a.flatten()
b = end_end_stack[:,0]
b = b.flatten()
polweight_stack = np.delete(polweight_stack, (0), axis = 1)
end_end_stack = np.delete(end_end_stack, (0), axis = 1)
pos_stack = np.delete(pos_stack, (0), axis = 1)
pos_stack = np.delete(pos_stack, (0), axis = 1)
return a, b, polweight_stack, end_end_stack, pos_stack
def stack_init(n_beads):
# initilize stacks to for storing data
polweight_stack = np.zeros([n_beads,1])
end_end_stack = np.zeros([n_beads, 1])
pos_stack = np.zeros([n_beads, 2])
polweight_stack[0], polweight_stack[1] = 1, 1
# Preallocate arrays
#pos = np.zeros([n_beads, 2])
# Initialization
pos_stack[1] = np.array([1, 0]) # bead 1: pos[0, 0]. bead 2: pos[1, 0]
#AddBead(Polymer,Polweight,L)
#Add third bead
return polweight_stack, end_end_stack, pos_stack
def exp_calc(weigths, exp_values):
end_end_times_weight = np.multiply(exp_values, weigths)
sum_end_end_times_weight = np.sum(end_end_times_weight, axis =1)
sum_weights = np.sum(weigths, axis=1)
exp_value = np.divide(sum_end_end_times_weight, sum_weights)
return exp_value
\ No newline at end of file
......@@ -24,16 +24,17 @@ def exp_calc(polweight_array, end_end_array, n_beads):
sum_weights = np.sum(polweight_array, axis=1)
exp_value = np.divide(sum_end_end_times_weight, sum_weights)
#foute manier van berekenen std, moet ook sws met kwadraten ofzoiets
std_weight_end = np.std(end_end_times_weight, axis=1)
std_weight = np.std(polweight_array, axis=1)
exp_value_std2 = np.divide(std_weight_end, std_weight)
exp_value_mean, exp_value_std = bootstrap(polweight_array, end_end_array, n_beads)
print('bootstrap = ', exp_value_std)
print('normal std = ', exp_value_std2)
print('bootstrap = ', exp_value_std.shape)
print('normal std = ', exp_value_std2.shape)
#print('bootstrap = ', exp_value_std)
#print('normal std = ', exp_value_std2)
#print('bootstrap = ', exp_value_std.shape)
#print('normal std = ', exp_value_std2.shape)
a = np.arange(0, n_beads - 1, 1)
exp_value[0:2] = [0, 1]
......@@ -43,32 +44,13 @@ def exp_calc(polweight_array, end_end_array, n_beads):
plt.title('bootstrap')
plt.xscale("log")
plt.yscale("log")
plt.errorbar(a, exp_value_mean[0:n_beads - 1] ** 2, xerr=0, yerr=np.sqrt(np.sqrt(a))*exp_value_std[:]**2, color="k")
plt.errorbar(a, exp_value_mean[0:n_beads - 1] ** 2, xerr=0, yerr=2*exp_value_std[:]**2, color="k")
plt.ylabel('R^2 (end to end distance)')
plt.xlabel('N (Number of beads)')
plt.grid(True, which="both", ls="-")
plt.ylim(1, 10000)
plt.xlim(1, 250)
plt.figure(2)
plt.title('normal std')
plt.xscale("log")
plt.yscale("log")
plt.errorbar(a, exp_value_mean[0:n_beads - 1] ** 2, xerr=0, yerr=exp_value_std2[1:] ** 2, color="k")
plt.ylabel('R^2 (end to end distance)')
plt.xlabel('N (Number of beads)')
plt.grid(True, which="both", ls="-")
plt.ylim(1, 10000)
plt.xlim(1, 250)
# plt.figure(2)
# plt.loglog(a, exp_value[0:n_beads - 1] ** 2, marker="x", markersize=7, color="k")
# plt.ylabel('R^2 (end to end distance)')
# plt.xlabel('N (Number of beads)')
# plt.grid(True, which="both", ls="-")
# plt.ylim(1, 10000)
# plt.xlim(1, 250)
plt.show()
def bootstrap(polweight_array, end_end_array, n_beads):
......@@ -94,16 +76,21 @@ def bootstrap(polweight_array, end_end_array, n_beads):
a = end_end_array.shape
nn = 100
exp_value = np.zeros((a[0]-1, nn))
for j in range(nn):
j = 0
while j < nn:
c = np.random.randint(a[1], size=(a[1]))
end_end_bootstrap = end_end_array[:, c]
polweight_bootstrap = polweight_array[:, c]
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 sum_weights_bootstrap == 0
# exp_value[:, j] =
exp_value[:, j] = np.divide(sum_end_end_times_bootstrap[0:n_beads - 1], sum_weights_bootstrap[0:n_beads - 1])
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])
j = j + 1
else:
j = j
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
\ No newline at end of file
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