Commit f4ed33e2 authored by Pieter van Velde's avatar Pieter van Velde

Implemented PERM algorithm in main_perm and functions_perm

parent 51bb0d44
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
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[j, 0]
polweight_stack[j, 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[j,0] = polweight_stack[j,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
......@@ -7,9 +7,9 @@ from functions_perm import add_bead, limmits, prune, enrich, maxsize
# Simulation parameters
n_beads = 100 # number of beads
n = 6 # number of angles
n_pol = 100 # number of polymers
n_pol = 1000 # number of polymers
alfa_up = 2 # alfa to compute upper limit
alfa_low = 1.2 # alfa to compute lwoer limit
alfa_low = 1.0 # alfa to compute lwoer limit
cst_dep = 1/(0.75*n) # multiplication constant for L dependence perm
# Physical parameters
......@@ -43,14 +43,15 @@ while k < n_pol-1:
j = 2
flag = 1
while flag == 1:
if j >= n_beads -2:
# print('pos = ', pos_stack )
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('size stack = ', numenrich - numprune)
print('numenrich = ', numenrich )
print('numprune = ',numprune )
print('k = ', k)
# print('pol = ', polweight_stack,)
# print( 'end = ', end_end_stack)
......@@ -72,17 +73,17 @@ while k < n_pol-1:
if polweight_stack[j,0] < lower_lim:
numprune += 1
polweight_stack, end_end_stack, pos_stack = prune(polweight_stack, end_end_stack, pos_stack, j)
#print('prune')
if polweight_stack.size == 0:
flag = 0
#print('stack = 0')
numprune += 1
elif polweight_stack[j,0] > upper_lim:
if j < n_beads -1:
polweight_stack, pos_stack, end_end_stack = enrich(polweight_stack,j, pos_stack, end_end_stack)
polweight_stack, pos_stack, end_end_stack = enrich(polweight_stack, pos_stack, end_end_stack, j)
#print('enrich')
numenrich += 1
if j == n_beads:
......@@ -94,8 +95,8 @@ while k < n_pol-1:
# plt.figure()
# plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
print(end_end_array)
print(polweight_array)
#print(end_end_array)
#print(polweight_array)
print('prune = ', numprune)
print('enrich = ', numenrich)
print('size stack = ', numenrich - numprune - n_pol)
......@@ -120,4 +121,10 @@ 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()
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