Commit 6bcfa381 authored by Luuk Balkenende's avatar Luuk Balkenende

Added all function descriptions to functions

parent 338f704b
# Import libaries/functions
import numpy as np
from matplotlib import pyplot as plt
import random
from functions_perm import calc_pol, allocate, init, perm, add_bead, limits, prune, enrich, maxsize, stack_init, exp_calc
from simfunctions import calc_pol, allocate, init, perm, add_bead, limits, prune, enrich
from processfunctions import exp_calc, bootstrap
# Simulation parameters
n_beads = 100 # number of beads
n_beads = 250 # number of beads
n = 6 # number of angles
n_pol = 500 # number of polymers
n_pol = 300 # number of polymers
alfa_up = 2 # alfa to compute upper limit
alfa_low = 1 # alfa to compute lwoer limit
cst_dep = 1 / (0.75 * n) # multiplication constant for L dependence perm
alfa_low = 0.7 # alfa to compute lower limit
PERM = False # boolean
random_walk = False # boolean
# Physical parameters
T_dim = 100 # T * (Kb / eps) # Non-dimensional temperature
# Choose algorithm
PERM = True # Put on True if PERM algorithm is desired
random_walk = False # Put on True if random walk is desired
# Calculate configuration of polymers and their weights
end_end_array, polweight_array = calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, cst_dep, T_dim, PERM, random_walk)
# Physical parameters
T_dim = 10 # T * (Kb / eps) # Non-dimensional temperature
exp_value = exp_calc(polweight_array, end_end_array)
a = np.arange(2, n_beads - 1, 1)
# 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)
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)
plt.show()
\ No newline at end of file
# 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)
import numpy as np
from matplotlib import pyplot as plt
def exp_calc(polweight_array, end_end_array, n_beads):
"""Calculates the end to end distances with taking into account their weights, and
makes a figure of the mean values of those distances against number of beads.
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
n_beads: integer > 2
Number of maximum beats in a polymer
Results:
--------
Figure: end to end distance against number of beads
"""
end_end_times_weight = np.multiply(end_end_array, polweight_array)
sum_end_end_times_weight = np.sum(end_end_times_weight, axis=1)
sum_weights = np.sum(polweight_array, axis=1)
exp_value = np.divide(sum_end_end_times_weight, sum_weights)
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)
a = np.arange(0, n_beads - 1, 1)
exp_value[0:2] = [0, 1]
exp_value_mean[0:2] = [0, 1]
plt.figure(1)
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.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):
"""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
n_beads: integer > 2
Number of maximum beats in a polymer
Results:
--------
exp_value_mean: array of size (n_beads, 1)
Mean value of end-to-end distance per number of beads
exp_value_std: array of size (n_beads, 1)
Standard deviation of end-to-end distances per number of beads
"""
a = end_end_array.shape
nn = 100
exp_value = np.zeros((a[0]-1, nn))
for j in range(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])
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
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