From 3a1364acc7b5dcbebfd9fda2d0dde1d93af5aaf2 Mon Sep 17 00:00:00 2001 From: pvanvelde Date: Wed, 8 May 2019 12:14:04 +0200 Subject: [PATCH] Change implementation of the PERM algorithm by storing the data of low weight polymers --- functions.py | 128 +++++++++++++++++++++++++++++++++++++++------------ main.py | 94 +++++++++++++++++++------------------ 2 files changed, 147 insertions(+), 75 deletions(-) diff --git a/functions.py b/functions.py index d21428d..848dd81 100644 --- a/functions.py +++ b/functions.py @@ -1,7 +1,6 @@ import numpy as np -import sys - -def add_bead(pos, n, T_dim, j, p_pol, fails): +import random +def add_bead(pos, n, T_dim, polweight, fails, cst_dep,j): """Adds one bead to polymer, according to the Rosenblaum method Parameters: @@ -22,10 +21,13 @@ def add_bead(pos, n, T_dim, j, p_pol, fails): 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 - flag = 0 + # Calculate new positions pos_new = np.zeros([n, 2]) pos_new[:, 0] = pos[j - 1, 0] + np.cos(phi) @@ -39,35 +41,103 @@ def add_bead(pos, n, T_dim, j, p_pol, fails): ry[i, :] = pos_new[i, 1] - pos[0:j, 1] r = np.sqrt(rx ** 2 + ry ** 2) r[r == 0] = np.nan - - # Calculate energy of different angles - E = np.sum(4 * ((1 / r) ** 12 - (1 / r) ** 6), axis=1) - - # Calculate (normalized) weights - w = np.exp(-E / T_dim) - w[np.isnan(w)] = 0 - w_tot = np.sum(w) - - # Stop execution if beads can only be placed (almost) on top of other beads - if w_tot < 1e-3: # BEETJE ME GEPIELD - #print('Help, I have nowhere to go') - #print('w_tot = ', w_tot) - flag = 1 - fails = fails + 1 + 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 / w_tot + 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 - # Calculate weights of one polymer to account for deviation from the theoretical value - w_idx = np.where(phi == final_phi) - w_jl = w[w_idx] - # print(w_jl) + +def enrich(polweight_stack, pos_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) + k = k + 1 + - w_l = w_tot - #print(w_tot) - p_pol = p_pol * w_jl / w_l - end_end = np.sqrt(pos[j, 0]**2 + pos[j,1]**2) + return polweight_stack, pos_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, k): + b = random.random() + if b < 0.5: + polweight_stack[:,k] = polweight_stack[:,k] * 2 + else: + polweight_stack = np.delete(polweight_stack, (k), axis = 1) + pos_stack = np.delete(pos_stack, [2*k, 2*k+1], axis = 1) + k = k-1 + + return pos_stack, polweight_stack, k + + + +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): + # 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, polweight_stack = 1, 1 + # Preallocate arrays + #pos = np.zeros([n_beads, 2]) + # Initialization + pos_stack = 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 - return pos, p_pol, flag, end_end, fails \ No newline at end of file +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 diff --git a/main.py b/main.py index 5215912..cbea062 100644 --- a/main.py +++ b/main.py @@ -1,64 +1,66 @@ # Import libaries/functions import numpy as np from matplotlib import pyplot as plt -import time +import random +from functionsperm2 import add_bead, limits, prune, enrich, maxsize, stack_init, exp_calc -from functions import add_bead -fails =0 -start = time.time() # Simulation parameters -n_beads = 250 #number of beads -n = 6 #number of angles -n_pol = 10000 #number of polymers +n_beads = 250 # number of beads +n = 6 # number of angles +n_pol = 100 # number of polymers +alfa_up = 100 # alfa to compute upper limit +alfa_low = 0 # 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 = 10 # T * (Kb / eps) # Non-dimensional temperature -p = np.zeros([n_beads, n_pol]) -end_end_matrix = np.zeros([n_beads, n_pol]) -# Make n polymers -for k in range(n_pol): - # Preallocate arrays - pos = np.zeros([n_beads, 2]) - - # Initialization - pos = np.array([1, 0]) #bead 1: pos[0, 0]. bead 2: pos[1, 0] - - # Make polymer of lenght n_beads - #w_jl = 1 # start value of weight w_j - #w_l = 1 # start value of weight W - p_pol = 1 - for j in range(2,n_beads): +m = np.arange(0, 2*n_pol, 2) +pos_stack = np.zeros([n_beads, n_pol*2]) +pos_stack[1, m] = 1 +#print(pos_stack) +polweight_stack = np.zeros([n_beads,n_pol]) +polweight_stack[0:2,:] = 1 +end_end_stack = np.zeros([n_beads,n_pol]) - pos, p_pol, flag, end_end, fails = add_bead(pos, n, T_dim, j, p_pol, fails) - if flag == 1: - #print('unsucesfull polymer') - #p[:, k] = 0 - k = k - 1 # throw polymer away - break - - if flag == 0: # keep polymer - p[j,k] = p_pol #Weight per polymer - end_end_matrix[j,k] = end_end +enrich = 0 +prune = 0 +flag = 1 -#print('P = ', p) -# Plot polymer -#plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" ) +fails = 0 +j = 2 +while flag == 1: + k=0 + while k < int(pos_stack.shape/2) : + pos_stack[:, 2*k:2*k+2], polweight_stack[j, k], end_end_stack[j,k], fails = add_bead(pos_stack[:, 2*k:2*k+2], n, T_dim, polweight_stack[j-1,k], fails, cst_dep,j) + + uplim, lowlim = limits(polweight_stack, alfa_up, k, alfa_low, j) + + if polweight_stack[j, k] < lowlim: + prune = (prune + 1)/2 + pos_stack, polweight_stack, k = prune(polweight_stack, pos_stack, k) -end = time.time() + elif polweight_stack[j, k] > uplim: + enrich = enrich + 1 + polweight_stack, pos_stack, k = enrich(polweight_stack, pos_stack, k) + + k = k +1 + if pos_stack.size == 0 or j== n_beads-1 : + flag = 0 + print('stack equal to zero at number of beads = ', j ) + j=j+1 -print('fails = ', fails) -average_end_end = np.average(end_end_matrix[2::], axis=1, weights=p[2::]) -a = np.arange(2,n_beads,1) - -plt.loglog(a, average_end_end**2, marker= "x", markersize=7, color="k" ) +exp_value = exp_calc(polweight_stack, end_end_stack) +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('time = ', end-start) -plt.show() - - +print(polweight_stack.shape) +print('enrich = ', enrich) +print('prune = ', prune) +plt.show() \ No newline at end of file -- 2.21.0