Commit 20f7a781 authored by Pieter van Velde's avatar Pieter van Velde

Small changes in commenting in main.py, simfunctions.py and processfunction.py

parent 12f945e4
import numpy as np
import random
def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, cst_dep, T_dim, PERM, random_walk):
# Allocate necessary arrays
pos_stack, pos_array, polweight_stack, polweight_array, end_end_stack,\
end_end_array = allocate(n_beads, n_pol)
# Initialise position of first two beads and their weighting factor
pos_stack, polweight_stack = init(n_pol, pos_stack, polweight_stack)
# Initialise some variables
n_enrich = 0 # amount of times enrichment was done
n_prune = 0 # amount of times pruning was done
flag1 = 1
j = 2
m = 0
# Calculate total polymers and their weights
while flag1 == 1:
k = 0
flag2 = 1
# Calculate next bead of polymer
while flag2 == 1:
# Calculate next bead of every polymer, and apply PERM
pos_stack, pos_array, polweight_stack, polweight_array, \
end_end_stack, end_end_array, j, k, m, n_prune, n_enrich = \
perm(pos_stack, pos_array, polweight_stack, polweight_array,
end_end_stack, end_end_array, j, k, m, n, n_prune, n_enrich,
T_dim, cst_dep, alfa_up, alfa_low, PERM, random_walk)
k = k + 1
if k >= int(pos_stack.shape[1] / 2):
flag2 = 0
if pos_stack.size == 0 or j == n_beads - 1:
flag1 = 0
print('stack equal to zero at number of beads = ', j)
j = j + 1
print('n_pol = ', polweight_stack.shape[1])
print(j)
# Append the two end2end matrices, and their weights
end_end_array = np.c_[end_end_array, end_end_stack]
polweight_array = np.c_[polweight_array, polweight_stack]
print('n_enrich = ', n_enrich)
print('n_prune = ', n_prune)
return end_end_array, polweight_array
def allocate(n_beads, n_pol):
pos_stack = np.zeros([n_beads, n_pol * 2])
pos_array = np.zeros([n_beads, 2])
polweight_stack = np.zeros([n_beads, n_pol])
polweight_array = np.zeros([n_beads])
end_end_stack = np.zeros([n_beads, n_pol])
end_end_array = np.zeros([n_beads, 1])
return pos_stack, pos_array, polweight_stack, polweight_array, end_end_stack, end_end_array
def init(n_pol, pos_stack, polweight_stack):
m = np.arange(0, 2 * n_pol, 2)
pos_stack[1, m] = 1
polweight_stack[0:2, :] = 1
return pos_stack, polweight_stack
def perm(pos_stack, pos_array, polweight_stack, polweight_array, end_end_stack, end_end_array,
j, k, m, n, n_prune, n_enrich, T_dim, cst_dep, alfa_up, alfa_low, PERM, random_walk):
# Avoid PERM and random_walk are true together
if random_walk == True and PERM == True:
print('ERROR: random_walk and PERM cannot be True both, change this in simulations parameters.')
exit()
# Add next bead to polymer
pos_stack[:, 2 * k:2 * k + 2], polweight_stack[j, k], end_end_stack[j, k] \
= add_bead(pos_stack[:, 2 * k:2 * k + 2], n, T_dim, polweight_stack[j - 1, k], cst_dep, j, random_walk)
# Calculate upper and lower limits for PERM algorithm
if PERM == True:
uplim, lowlim = limits(polweight_stack, alfa_up, k, alfa_low, j)
# Conditions for pruning
if polweight_stack[j, k] < lowlim:
n_prune = n_prune + 1 # (prune1 + 1) / 2
pos_stack, polweight_stack, pos_array, polweight_array, end_end_stack, end_end_array, k, m \
= prune(polweight_stack, pos_stack, pos_array, polweight_array,
end_end_stack, end_end_array, k, m)
# Conditions for enriching
elif polweight_stack[j, k] > uplim:
n_enrich = n_enrich + 1
polweight_stack, pos_stack, end_end_stack, k \
= enrich(polweight_stack, pos_stack, end_end_stack, k)
return pos_stack, pos_array, polweight_stack, polweight_array, \
end_end_stack, end_end_array, j, k, m, n_prune, n_enrich
def add_bead(pos, n, T_dim, polweight, cst_dep, j, random_walk):
"""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
"""
# 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)
if random_walk == True:
w_jl = np.ones(n)
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) *10**-50
w_l = np.sum(w_jl)
print('fail')
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
def enrich(polweight_stack, pos_stack, end_end_stack, k):
polweight_stack[:, k] = 0.5 * polweight_stack[:, k]
polweight_stack = np.insert(polweight_stack, k, polweight_stack[:, k], axis=1)
pos_stack = np.insert(pos_stack, 2 * k, pos_stack[:, 2 * k + 1], axis=1)
pos_stack = np.insert(pos_stack, 2 * k, pos_stack[:, 2 * k + 1], axis=1)
end_end_stack = np.insert(end_end_stack, k, end_end_stack[:,k], axis=1)
k = k + 1
return polweight_stack, pos_stack, end_end_stack, k
def limits(polweight_stack, alfa_up, k, alfa_low, j):
# find all nonzero elements of stack & pw_matrix, of bead j, sum and average:
avg_weight = np.average(polweight_stack[j, :k + 1])
weight3 = polweight_stack[2, k]
upper = (alfa_up * avg_weight) / weight3
lower = (alfa_low * avg_weight) / weight3
return upper, lower
def prune(polweight_stack, pos_stack, pos_array, polweight_array, end_end_stack, end_end_array, k, m):
b = random.random()
# print('b =', b)
if b < 0.5:
polweight_stack[:, k] = polweight_stack[:, k] * 2
else:
polweight_array = np.c_[polweight_array, polweight_stack[:, k]]
#polweight_array[:, m] = polweight_stack[:, k]
polweight_stack = np.delete(polweight_stack, k, axis=1)
pos_array = np.c_[pos_array, pos_stack[:, 2*k]]
pos_array = np.c_[pos_array, pos_stack[:, 2*k+1]]
# pos_array[:, 2*m] = pos_stack[:, 2*k]
# pos_array[:, 2*m+1] = pos_stack[:, 2*k+1]
pos_stack = np.delete(pos_stack, [2 * k, 2 * k + 1], axis=1)
end_end_array = np.c_[end_end_array, end_end_stack[:, k]]
end_end_stack = np.delete(end_end_stack, k, axis=1)
k = k - 1
m = m + 1
return pos_stack, polweight_stack, pos_array, polweight_array, end_end_stack, end_end_array, k, m
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):
# initialize 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
# Import libaries/functions
from simfunctions import calc_pol, allocate, init, perm, add_bead, limits, prune, enrich
from processfunctions import exp_calc, bootstrap
# Import functions
from simfunctions import calc_pol
from processfunctions import exp_calc
# Simulation parameters
n_beads = 250 # number of beads
n = 6 # number of angles
n_pol = 1000 # number of polymers
alfa_up = 2 # alfa to compute upper limit
alfa_low = 1.2 # alfa to compute lower limit
n_pol = 200 # number of polymers
alfa_up = 2 # alfa to compute upper limit for PERM
alfa_low = 1.2 # alfa to compute lower limit for PERM
# Physical parameters
T_dim = 10 # T * (Kb / eps) # Non-dimensional temperature
# Choose algorithm
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 = 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)
# Process configuration of polymers and their weights into figures (end2end distance against n_beads)
# Process configuration of polymers and their weights into figures
# and preform bootstrap to calculate errors.
exp_calc(polweight_array, end_end_array, n_beads, pos_array)
......@@ -20,30 +20,32 @@ def exp_calc(polweight_array, end_end_array, n_beads, pos_array):
--------
Figure: end to end distance against number of beads
"""
# Compute mean values and standard deviation using bootstrap.
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
# Compute population size
popsize = polsize(end_end_array)
# Plot data
font = {'family': 'DejaVu Sans',
'weight': 'normal',
'size': 17}
matplotlib.rc('font', **font)
# Create fit for end to end distance
a = np.arange(0, n_beads - 1, 1)
q = np.polyfit(np.log(a[1::]), np.log(exp_value_mean[1::] ** 2), 1)
print(q)
print('Fitting parameters for end2end distance -> ', q)
xx = np.linspace(0, 250, 1000)
yy = np.exp(q[1]) * xx ** q[0]
plt.figure(1)
# Plot end to end distance
plt.figure(1)
plt.xscale("log")
plt.yscale("log")
plt.plot(xx, yy, 'r--', label="Fit", linewidth=3)
......@@ -62,25 +64,23 @@ def exp_calc(polweight_array, end_end_array, n_beads, pos_array):
plt.ylim(top=10000)
plt.savefig('endend', bbox_inches='tight', dpi=300)
###########################################################################
# Create fit for gyraton radius.
q1 = np.polyfit(np.log(a[1::]), np.log(rad_gyr_mean[1::]), 1)
print(q1)
print('Fitting parameters for gyration radius -> ', q1)
xx1 = np.linspace(0, 250, 1000)
yy1 = np.exp(q1[1]) * xx ** q1[0]
# Plot gyraton radius.
plt.figure(2)
plt.xscale("log")
plt.yscale("log")
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.grid(True, which="both", ls=":")
plt.legend(loc='best')
plt.xlim(1, 250)
......@@ -88,26 +88,25 @@ def exp_calc(polweight_array, end_end_array, n_beads, pos_array):
plt.ylim(top=1000)
plt.savefig('gyros', bbox_inches='tight', dpi=300)
# Plot population size
plt.figure(3)
plt.title('Population size')
# plt.yscale("log")
plt.plot(b[2:], popsize[2:])
plt.plot(np.arange(2, n_beads, 1), popsize[2:])
plt.grid(True, which="both", ls="-")
plt.ylabel('Population size')
plt.xlabel('N (Number of beads)')
plt.show()
def bootstrap(polweight_array, end_end_array, n_beads):
def bootstrap(polweight_array, exp_value_array, n_beads):
"""Preforms bootstrap onto end-to-end distances of polymers to calculate the error
Parameters:
----------
polweight_array: array of size (n_beads, n_poll + amount of enriching - 0.5 times amount of pruning)
Array with all weights of all beads of all polymers
end_end_array: array of size (n_beads, n_poll + amount of enriching - 0.5 times amount of pruning)
Distance of bead to center bead, for every bead of every simulated polymer
exp_value_array: array of size (n_beads, n_poll + amount of enriching - 0.5 times amount of pruning)
Array at which the bootstrap is applied
n_beads: integer > 2
Number of maximum beats in a polymer
......@@ -118,12 +117,10 @@ 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
a = exp_value_array.shape
exp_value = np.zeros((a[0] - 1, nn)) # Allocate observable
j = 0
m = 0
......@@ -132,29 +129,27 @@ def bootstrap(polweight_array, end_end_array, n_beads):
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]
exp_value_bootstrap = exp_value_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)
exp_value_times_bootstrap = np.multiply(exp_value_bootstrap, polweight_bootstrap)
sum_exp_value_times_bootstrap = np.sum(exp_value_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],
exp_value[:, j] = np.divide(sum_exp_value_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:
elif np.count_nonzero(exp_value_array[n_beads - 1, :]) == 0:
print('No polymer reached ', n_beads, ' beads, try again, or lower variable n_beads.')
exit()
elif m > 100:
elif m > 500:
print('Not enough polymers reached ', n_beads, ' beads, try again, or adjust variable n_beads.')
exit()
else:
m = m + 1
print('fail')
# print(sum_weights_bootstrap)
m = m + 1
# Calculate mean and standard deviation of observable
exp_value_mean = np.mean(exp_value, axis=1)
......@@ -163,15 +158,36 @@ def bootstrap(polweight_array, end_end_array, n_beads):
def polsize(end_end_array):
# HIER NOG EEN DOCSTRING DINGENS!!
"""Calulates the population size of polymers for different lenghts
Parameters:
----------
end_end_array: array of size (n_beads, n_poll + amount of enriching - 0.5 times amount of pruning)
Distance of bead to center bead, for every bead of every simulated polymer.
Results:
--------
popsize: array of size (n_beads)
population size of polymers for different lenghts
"""
# Count nonzero elements
popsize = np.count_nonzero(end_end_array, axis=1)
return popsize
def radiusgyration(pos_array):
# HIER NOG EEN DOCSTRING DINGENS!
"""Calulates the radius of gyration.
Parameters:
----------
pos_array: array of size (n_beads, ?)
Array with positions of the beads of polymers which stopped growing
Results:
--------
rad_gyr: array of size (n_beads, n_poll + amount of enriching - 0.5 times amount of pruning)
Array containing the radius of gyration for all polymers at different lenghts
"""
# Preparation for calculation can start (allocation, initialization)
pos_array = np.delete(pos_array, [0, 1], axis=1)
......
......@@ -65,10 +65,10 @@ def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk):
# Loop over all n_beads, and abort while loop if all polymers have their final length
if pos_stack.size == 0 or j == n_beads - 1:
flag1 = 0
print('stack equal to zero at number of beads = ', j)
print('n_pol = ', polweight_stack.shape[1])
print(j)
print('stack equal to zero at number of beads = ', j + 1)
if PERM == True:
print('number of polymers = ', polweight_stack.shape[1])
print('polymers have lenght: ', j+1, 'beads')
j = j + 1
# Append the two end2end matrices, the position matrices, and their weights, an
......@@ -375,10 +375,11 @@ def limits(polweight_stack, alfa_up, k, alfa_low, j):
Number to indicate when to prune
"""
# find all nonzero elements of stack & pw_matrix, of bead j, sum and average:
# Find all nonzero elements of stack & pw_matrix, of bead j, sum and average:
avg_weight = np.average(polweight_stack[j, :k+1])
weight3 = polweight_stack[2, k] # Average weights of first 3 beads
# Limit computation
upper = (alfa_up * avg_weight) / weight3
lower = (alfa_low * avg_weight) / weight3
......
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