Commit 8ef5e54a authored by Pieter van Velde's avatar Pieter van Velde

Comment and changes weight calculation in functions_perm and main_perm

parent f4ed33e2
......@@ -66,8 +66,8 @@ def add_bead(pos, n, T_dim, polweight, fails, cst_dep):
def enrich(polweight_stack, pos_stack, end_end_stack, j):
new_weight = 0.5 * polweight_stack[j, 0]
polweight_stack[j, 0] = new_weight
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]]
......@@ -105,7 +105,7 @@ def limmits(polweight_matrix, polweight_stack, alfa_up, j, alfa_low):
def prune(polweight_stack, end_end_stack, pos_stack, j ):
b = random.random()
if b < 0.5:
polweight_stack[j,0] = polweight_stack[j,0] * 2
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)
......@@ -128,3 +128,24 @@ def maxsize(polweight_stack, end_end_stack, pos_stack):
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
......@@ -2,112 +2,96 @@
import numpy as np
from matplotlib import pyplot as plt
import random
from functions_perm import add_bead, limmits, prune, enrich, maxsize
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 = 1000 # number of polymers
alfa_up = 2 # alfa to compute upper limit
alfa_low = 1.0 # alfa to compute lwoer limit
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 = 100 # T * (Kb / eps) # Non-dimensional temperature
T_dim = 1 # T * (Kb / eps) # Non-dimensional temperature
polweight_stack = np.zeros([n_beads,1])
end_end_stack = np.zeros([n_beads, 1])
pos_stack = np.zeros([n_beads, 2])
# 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
# Make n polymers
k = 0
fails = 0
# Loop until the desired amount of polymers is achieved.
while k < n_pol-1:
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
# 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 -2:
# print('pos = ', pos_stack )
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 )
print('numprune = ',numprune/2 )
print('k = ', k)
# print('pol = ', polweight_stack,)
# print( 'end = ', end_end_stack)
# print('pos = ', pos_stack )
# 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
pos_stack[:,:2], polweight, end_end, fails, j = add_bead(pos_stack[:,:2], n, T_dim, polweight_stack[:,0], fails, cst_dep)
end_end_stack[j,0] = end_end
# 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)
#print('prune')
# 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)
#print('enrich')
numenrich += 1
if j == n_beads:
print('error maybe')
#print('size stack = ', numenrich - numprune)
# 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)
print('prune = ', numprune/2)
print('enrich = ', numenrich)
print('size stack = ', numenrich - numprune - n_pol)
end_end_times_weight = np.multiply(end_end_array, polweight_array)
print('size stack = ', numenrich - numprune/2 - n_pol)
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)
a = np.arange(2,n_beads-1,1)
plt.figure()
......@@ -128,3 +112,5 @@ bins = np.arange(0, n_beads, 1)
plt.figure()
plt.hist(hist, bins)
plt.show()
print('joe')
\ 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