From 4783a9778aab933c1e4e554b4168f061b553686d Mon Sep 17 00:00:00 2001 From: Benjamin Maier Date: Mon, 11 Mar 2019 18:12:12 +0100 Subject: [PATCH] more accurate rates in eq.-state-trajectory --- tacoma/flockwork.py | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/tacoma/flockwork.py b/tacoma/flockwork.py index eb37c8b..7206629 100644 --- a/tacoma/flockwork.py +++ b/tacoma/flockwork.py @@ -459,16 +459,21 @@ def flockwork_P_group_life_time_distributions_for_varying_alpha_beta(tau,max_gro # from equilibrium assumption k = P/(1-P) compute adjusted P new_P = k / (k+1) gamma = alpha + beta - new_gamma = tc.sample_a_function(t, gamma, new_t) + #new_gamma = tc.sample_a_function(t, gamma, new_t) + new_alpha = tc.sample_a_function(t, alpha, new_t) distros = [[] for m in range(min_group_size, max_group_size+1) ] ks = np.arange(N) # for every time point and adjusted P, compute the equilibrium group size distribution - for g_, P_ in zip(new_gamma, new_P): + for a_, P_, _k_ in zip(new_alpha, new_P, k): + if P_ > 0: + g_ = a_ / P_ + else: + g_ = 0 for m in range(min_group_size, max_group_size+1): - lambda_m = m*g_*(1-P_)+2*g_*P_*m*(N-m)/(N-1.0) + lambda_m = m*g_*(1-P_)+2*a_*m*(N-m)/(N-1.0) y = lambda_m * np.exp(-lambda_m*tau) distros[m-min_group_size].append(y) @@ -518,18 +523,24 @@ def flockwork_P_contact_time_distributions_for_varying_alpha_beta(tau,N,k_initia new_P = k / (k+1) gamma = alpha + beta new_gamma = tc.sample_a_function(t, gamma, new_t) + new_alpha = tc.sample_a_function(t, alpha, new_t) + distro_c = [] distro_ic = [] ks = np.arange(N) # for every time point and adjusted P, compute the equilibrium group size distribution - for g_, P_ in zip(new_gamma, new_P): + for a_, P_, _k_ in zip(new_alpha, new_P, k): + if P_ > 0: + g_ = a_ / P_ + else: + g_ = 0 p_k = degree_distribution(N,P_) _k_ = p_k.dot(ks) _k2_ = p_k.dot(ks**2) omega = 2*g_*(1-P_/(N-1)*_k2_/_k_) - lambda_1 = 2*g_*P_ + lambda_1 = 2*a_ this_distro_c = omega*np.exp(-tau*omega) this_distro_ic = lambda_1 * np.exp(-lambda_1*tau) distro_c.append(this_distro_c)