Commit ce176b3a authored by Luuk Balkenende's avatar Luuk Balkenende

Added population size calculation and processing in simfunctions and processfunctions.py

parent 1bc4890e
......@@ -17,8 +17,8 @@ random_walk = False # Put on True if random walk is desired
T_dim = 10 # T * (Kb / eps) # Non-dimensional temperature
# Simulate configuration of polymers and their weights
end_end_array, polweight_array = calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk)
end_end_array, polweight_array, pol_size = calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk)
# Process configuration of polymers and their weights into figures (end2end distance against n_beads)
# and preform bootstrap to calculate errors.
exp_calc(polweight_array, end_end_array, n_beads)
exp_calc(polweight_array, end_end_array, n_beads, pol_size)
# Import libaries/functions
import numpy as np
from matplotlib import pyplot as plt
import random
from functions_perm import add_bead, limmits, prune, enrich, maxsize, stack_init, exp_calc
# Simulation parameters
n_beads = 100 # number of beads
n = 6 # number of angles
n_pol = 100 # number of polymers
alfa_up = 1.5 # alfa to compute upper limit
alfa_low = 0.5 # alfa to compute lwoer limit
cst_dep = 1/(0.75*n) # multiplication constant for L dependence perm
# Physical parameters
T_dim = 1 # T * (Kb / eps) # Non-dimensional temperature
# Initial global variables
end_end_array = np.zeros([n_beads, n_pol-1])
polweight_array = np.zeros([n_beads, n_pol-1])
numprune = 0
numenrich = 0
k = 0
fails = 0
# Loop until the desired amount of polymers is achieved.
while k < n_pol-1:
# Initialize stacks, which are used to store polymerts with lengt < n_beads
polweight_stack, end_end_stack, pos_stack = stack_init(n_beads)
j = 2
flag = 1
while flag == 1:
if j >= n_beads -1:
# if the polymer in the stack reaches n_beads, it is stored and discarded
polweight_array[:, k], end_end_array[:, k], polweight_stack, end_end_stack, pos_stack = maxsize(polweight_stack, end_end_stack, pos_stack)
k = k + 1
print('numenrich = ', numenrich )
print('numprune = ',numprune/2 )
print('k = ', k)
# if stack is empty, break and initialize stack again using stack_init
if polweight_stack.size == 0:
break
# if n_pol is reached break
if k == n_pol - 1:
break
# Add next bead to the first polymer in the stack and compute/store weight and end_end distance and position
pos_stack[:,:2], polweight, end_end , fails, j = add_bead(pos_stack[:,:2], n, T_dim, polweight_stack[:,0], fails, cst_dep)
polweight_stack[j,0] = polweight
end_end_stack[j,0] = end_end
# Deduce upper and lower limit for PERM, for first polymer in stack, using average of all polymers in stack + the final array
upper_lim, lower_lim = limmits(polweight_array, polweight_stack, alfa_up, j, alfa_low)
# Prune
if polweight_stack[j,0] < lower_lim:
numprune += 1
# Delete first polymer from stacks or increase weight by a factor 2, both with change 1/2.
polweight_stack, end_end_stack, pos_stack = prune(polweight_stack, end_end_stack, pos_stack, j)
# If stack is empty, leave loop, and to create new stack
if polweight_stack.size == 0:
flag = 0
#print('stack = 0')
# Enrich
elif polweight_stack[j,0] > upper_lim:
# Copy the first polymer into the stack while decreasing weight by a factor 2.
if j < n_beads -1:
polweight_stack, pos_stack, end_end_stack = enrich(polweight_stack, pos_stack, end_end_stack, j)
numenrich += 1
# Compute expetation value end distance polymers
exp_value = exp_calc(polweight_array, end_end_array)
#------------------------------------------------------- Plot/print values -------------------------------------------------------
# plt.figure()
# plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
#print(end_end_array)
#print(polweight_array)
print('prune = ', numprune/2)
print('enrich = ', numenrich)
print('size stack = ', numenrich - numprune/2 - n_pol)
a = np.arange(2,n_beads-1,1)
plt.figure()
plt.loglog(a, exp_value[2: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)
print('fails = ', fails)
#plt.figure()
# plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
print(polweight_stack.shape)
hist = np.count_nonzero(polweight_stack, axis=1)
bins = np.arange(0, n_beads, 1)
#histi = np.histogram(hist, n_beads)
plt.figure()
plt.hist(hist, bins)
plt.show()
print('joe')
\ No newline at end of file
import numpy as np
def add_bead(pos, n, T_dim, L, polweight):
# 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[L - 1, 0] + np.cos(phi)
pos_new[:, 1] = pos[L - 1, 1] + np.sin(phi)
# Calculate distance between new bead and beads
rx = np.zeros([n, L])
ry = np.zeros([n, L])
for i in range(n):
rx[i, :] = pos_new[i, 0] - pos[0:L, 0]
ry[i, :] = pos_new[i, 1] - pos[0:L, 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)
print(w_l)
#print(w_l)
polweight = polweight*w_l #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[L] = np.array([[pos[L - 1, 0] + np.cos(final_phi), pos[L - 1, 1] + np.sin(final_phi)]])
return pos, polweight
\ No newline at end of file
# Import libaries/functions
import numpy as np
from matplotlib import pyplot as plt
from new_functions import add_bead
# Simulation parameters
n_beads = 250 # number of beads
n = 6 # number of angles
n_pol = 10 # number of polymers
# Physical parameters
T_dim = 0.1 # T * (Kb / eps) # Non-dimensional temperature
polweight = 1
# Make n polymers
k = 0
while n_pol > k:
# Preallocate arrays
pos = np.zeros([n_beads, 2])
# Initialization
pos[1] = np.array([1, 0]) # bead 1: pos[0, 0]. bead 2: pos[1, 0]
#AddBead(Polymer,Polweight,L)
#Add third bead
L = 2
while L < n_beads:
#add next bead
pos, polweight = add_bead(pos, n, T_dim, L, polweight)
L = L +1
k = k+1
print(k)
plt.figure()
plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
plt.show()
\ No newline at end of file
......@@ -37,6 +37,8 @@ def exp_calc(polweight_array, end_end_array, n_beads):
#print('normal std = ', exp_value_std2.shape)
a = np.arange(0, n_beads - 1, 1)
b = np.arange(0, n_beads, 1)
exp_value[0:2] = [0, 1]
exp_value_mean[0:2] = [0, 1]
......@@ -50,6 +52,14 @@ def exp_calc(polweight_array, end_end_array, n_beads):
plt.grid(True, which="both", ls="-")
plt.ylim(1, 10000)
plt.xlim(1, 250)
plt.figure(3)
plt.title('Population size')
#plt.yscale("log")
plt.plot(b[2:], popsize[2:])
plt.grid(True, which="both", ls="-")
plt.ylabel('Population size')
plt.xlabel('N (Number of beads)')
plt.show()
......@@ -93,4 +103,11 @@ def bootstrap(polweight_array, end_end_array, n_beads):
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
return exp_value_mean, exp_value_std
def polsize(end_end_array):
popsize = np.count_nonzero(end_end_array, axis=1)
print(popsize)
print(popsize.shape)
return popsize
\ No newline at end of file
......@@ -40,6 +40,9 @@ def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk):
flag1 = 1
j = 2
m = 0
pol_size = np.zeros(n_beads)
pol_size[0] = n_beads
pol_size[1] = n_beads
# Calculate total polymers and their weights
while flag1 == 1:
......@@ -64,9 +67,11 @@ def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk):
flag1 = 0
print('stack equal to zero at number of beads = ', j)
j = j + 1
print('n_pol = ', polweight_stack.shape[1])
print(j)
pol_size[j] = polweight_stack.shape[1]
j = j + 1
# Append the two end2end matrices, and their weights
end_end_array = np.c_[end_end_array, end_end_stack]
......@@ -75,7 +80,7 @@ def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk):
print('n_enrich = ', n_enrich)
print('n_prune = ', n_prune)
return end_end_array, polweight_array
return end_end_array, polweight_array, pol_size
def allocate(n_beads, n_pol):
"""Allocates necessary arrays to store positions and weights of beads of polymers
......
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