# Source code for tacoma.flockwork

# -*- coding: utf-8 -*-
"""
This module provides functions related to the flockwork
temporal network model.
"""
from copy import deepcopy

import numpy as np
from numpy import random

from scipy.integrate import ode
from scipy.integrate import simps
from scipy.special import gamma as Gamma

from tacoma.power_law_fitting import fit_power_law_clauset

from _tacoma import flockwork_P_varying_rates
from tacoma import _get_raw_temporal_network

import tacoma as tc

[docs]def flockwork_P_equilibrium_group_size_distribution(N, P):
"""Get the equilibrium group size distribution of a Flockwork-P model
given node number N and probability to reconnect P.

Parameters
----------
N : int
Number of nodes
P : float
Probability to reconnect

Returns
-------
numpy.ndarray
Group size distribution of this configuration. The m-th entry
of this array contains the expected number of nodes of groups
of size m.
"""

N = int(N)
P = float(P)

assert N > 2
assert P >= 0.
assert P <= 1.

if P == 0.:
dist = np.zeros((N+1,))
dist = N
elif P == 1.:
dist = np.zeros((N+1,))
dist[-1] = 1.
else:
dist = [N*(1.-P)]

N_fak = N - np.arange(1, N-1)
j_fak = ((N-1)-P * np.arange(1, N-1))
div = N_fak / j_fak
cum_product_div = np.cumprod(div)
for m in range(2, N):
#dist.append( (-1)**(m%2) * float(N)/float(m) * (P-1.) * np.prod(N_fak[1:m]/j_fak[1:m]) * P**(m-1) )
dist.append(float(N)/float(m) *
(1-P) * cum_product_div[m-2] * P**(m-1))

value = P
for j in range(1, N-1):
#value *= float(N-j-1) / ((P-N+1.) / P + (j-1))
value *= float(N-j-1) * P / ((N-1)-P*j)
#value *= P**(N-1)
dist.append(value)

dist = [0.] + dist

return np.array(dist)

[docs]def group_size_distribution_asymptotics(N, P, mmax=None, simple_pochhammer_approximation=True):
"""Get the asymptotic equilibrium group size distribution of a Flockwork-P model
given node number N and probability to reconnect P.

Parameters
----------
N : int
Number of nodes
P : float
Probability to reconnect
mmax : int, default = None
The maximum group size for which to calculate the asymptotics
simple_pochhammer_approximation : bool, default = True
Whether to use a simple first order Pochhammer approximation
or the Gamma-function approximation.

Returns
-------
ms : numpy.ndarray
group size vector
dist : numpy.ndarray
Asymptotic group size distribution
"""

N = int(N)
P = float(P)

assert N > 2
assert P >= 0.
assert P < 1.

if mmax is None:
mmax = int(N) // 2

ms = np.arange(1,mmax+1)

if P == 0.:
ms = np.arange(N+1)
dist = np.zeros((N+1,))
dist = N
else:
dist = [(1.-P)]

for m in ms[1:]:
if simple_pochhammer_approximation:
factor = ( 1 - P/(N-1)*(m-1)*(m-2)/2 )**(-1)
else:
factor = ((-1)**(m % 2) * (P*m/(N-1)/np.exp(1))**m * 2*np.pi / Gamma(-(N-1)/P) * m**(-(N-1)/P - 0.5))**(-1)
this = 1/m * (P/np.exp(1))**(m-1) * ((N-1)/(N-m))**(N-m+0.5) * factor

dist.append(this)

return ms, N*np.array(dist)

[docs]def flockwork_P_equilibrium_configuration(N,
P,
shuffle_nodes=True,
return_histogram=False,
seed=0,
shuffle_group_sizes=False,
dist=None,
use_binomial=True,
construct_edges=True,
):
"""Get an equilibrium configuration of a Flockwork-P model
given node number N and probability to reconnect P.

Parameters
----------
N : int
Number of nodes
P : float
Probability to reconnect
shuffle_nodes : bool, default : True
Shuffle the node order in which nodes are distributed
to groups. 'True' is recommended.
return_group_size_histogram : bool, default : False
Return a numpy array containing the counts of groups
of size :math:g.
seed : int, default : 0
The random seed. RNG is initialized randomly if seed = 0.
shuffle_group_sizes : bool, default : False
Shuffle the order of group sizes in which nodes are distributed
to groups. 'True' is recommended.
dist : np.array, default : None
Is computed from N and P if None is provided
use_binomial : bool, default : True
whether to draw from a Binomial instead of a Poisson
(conceptually more sound)
construct_edgelist : bool, default : True
if to actually construct the edgelist.

Returns
-------
:obj:list of :obj:tuple of int
edge list of equilibrium configuration
numpy.ndarray
group size counter of this configuration (only
if return_group_size_histogram is True)
"""

if seed > 0:
random.seed(seed)

if dist is None:
dist = flockwork_P_equilibrium_group_size_distribution(N, P)

dist = dist[1:]

total_group_number = sum(dist)
size_dist = dist / total_group_number
num_components = []
edges = []

if shuffle_nodes:
node_ints = random.permutation(N)
else:
node_ints = np.arange(N)

# in the beginning, theres N nodes left
nodes_left = N

# ... and no group counter contains any groups yet
C_m = np.zeros(N)

# while there are still nodes left to distribute
while nodes_left > 0:

# loop through group sizes in descending order.
# may contain all of the nodes left
if shuffle_group_sizes:
group_sizes = 1 + random.permutation(nodes_left)
else:
group_sizes = np.arange(1,nodes_left+1)[::-1]
for m in group_sizes:

if (nodes_left-m) < 0:
continue

# if the expected number of groups of this size is not zero
if dist[m-1] > 0. and nodes_left >= m:

#new_N_m = random.binomial(total_group_number,size_dist[m-1])
# dist carries the expected number of m-groups in equilibrium,
# so we draw a number of m groups from a Poisson distribution
# with mean dist[m-1]
if use_binomial:
max_draws = N // m
p = dist[m-1] / max_draws
new_C_m = random.binomial(max_draws,p)
else:
new_C_m = random.poisson(dist[m-1])

# if the new number of groups of size m is larger than the previously drawn number
if new_C_m > C_m[m-1]:
# accept the new number and add the difference in group count of this group size
delta_C_m = int(new_C_m - C_m[m-1])
elif nodes_left == 1 and m == 1:
# if there's only one node left, just add it to the loners
delta_C_m = 1
else:
delta_C_m = 0

for group_instance in range(delta_C_m):
# for u in xrange(nodes_left-1,nodes_left-m,-1):
#    for v in xrange(u-1,nodes_left-m-1,-1):
#        edges.append((node_ints[u],node_ints[v]))
if nodes_left - m < 0:
break

# add fully connected clusters to the edge set
if construct_edges:
if m > 1 and nodes_left-m >= 0:
edges.extend([tuple(sorted((int(node_ints[u]), int(node_ints[v]))))
for u in range(nodes_left-m, nodes_left-1)
for v in range(u+1, nodes_left)])
elif nodes_left - m < 0:
break

# remove the grouped nodes from the pool of remaining nodes
nodes_left -= m

# save this group as an accepted group
# and increment the group count
C_m[m-1] += 1

"""
# if there's no nodes left to distribute,
# N_m was chosen too large and we abandon
if nodes_left <= 0.:
break
"""

# go to next smaller group size

if return_histogram:
return edges, np.append(np.array([0.]), C_m)
else:
return edges

[docs]def flockwork_P_mean_degree_for_varying_rates(flockwork_P_params, N=None):
"""Compute the mean group size distribution for a Flockwork-P system with varying rates.

Parameters
----------
flockwork_P_params : :obj:dict
Contains all parameters necessary for a Flockwork-P simulation, especially the
time-dependent rewiring rate and time-dependent reconnection probability
N : int, default : None
If given, compute everything for N nodes, where N is different from
N in flockwork_P_params.

Returns
-------
t : numpy.array
An array of times at which the mean degree was evaluated
k : numpy.array
An array of mean degree values corresponding to the times in t.
"""

data = flockwork_P_params

if N is None:
N = data['N']

# compute the initial mean degree
k0 = 2*len(data['E']) / float(N)

# get arrays for rewiring
gamma = np.array(data['rewiring_rate'])
T, gamma = gamma[:, 0], gamma[:, 1]
T = np.append(T, data['tmax'])
gamma = np.append(gamma, gamma[-1])

P = np.array(data['P'])
P = np.append(P, P[-1])

# define the linear ODE describing the mean degree evolution
def dkdt(t, k, g_, P_):
return 2 * g_ * P_ - 2*k * (g_*(1-P_))

# intialize the integrator
r = ode(dkdt)
r.set_integrator('dopri5')

k = [k0]
new_t = [T]

i = 0
for t_, g_, P_ in zip(T[:-1], gamma[:-1], P[:-1]):

# for every interval set the initial mean degree
# and new parameters
r.set_initial_value(k[-1], t_)
r.set_f_params(g_, P_)

# for 10 time points within this interval,
# compute the expected mean degree
this_t = np.linspace(t_, T[i+1], 10)
for t__ in this_t[1:]:
new_t.append(t__)
this_k = r.integrate(t__)
k.append(this_k)

# increase interval
i += 1

return np.array(new_t), np.array(k)

[docs]def flockwork_P_mean_degree_for_varying_alpha_beta(N,k_initial,t,alpha,beta,tmax,sampling_points=10):
r"""Compute the mean group size distribution for a Flockwork-P system with varying rates.

Parameters
----------
N : int
Number of nodes
k_initial : float
initial mean degree
t : numpy.ndarray of float
time points at which :math:\alpha(t) and
:math:\beta(t) change
alpha : numpy.ndarray of float
active reconnection rate associated with
the time points in t
beta : numpy.ndarray of float
active disconnection rate associated with
the time points in t
tmax : float
final time
sampling_points : int, default : 10
how many points to sample in between two time points in t

Returns
-------
t : numpy.array
An array of times at which the mean degree was evaluated
k : numpy.array
An array of mean degree values corresponding to the times in t.
"""

# compute the initial mean degree
k0 = k_initial

# get arrays for rewiring
t = np.append(t, tmax)
alphas = np.append(alpha, alpha[-1])
betas = np.append(beta,beta[-1])

# define the linear ODE describing the mean degree evolution
def dkdt(t, k, a_, b_):
return 2 * a_ - 2*k * (b_)

# intialize the integrator
r = ode(dkdt)
r.set_integrator('dopri5')

k = [k0]
new_t = [t]

i = 0
for t_, a_, b_ in zip(t[:-1], alphas[:-1], betas[:-1]):

# for every interval set the initial mean degree
# and new parameters
r.set_initial_value(k[-1], t_)
r.set_f_params(a_, b_)

# for 10 time points within this interval,
# compute the expected mean degree
this_t = np.linspace(t_, t[i+1], sampling_points)
for t__ in this_t[1:]:
new_t.append(t__)
this_k = r.integrate(t__)
k.append(this_k)

# increase interval
i += 1

return np.array(new_t), np.array(k)

[docs]def flockwork_P_group_life_time_distributions_for_varying_alpha_beta(tau,max_group_size,N,k_initial,t,alpha,beta,tmax,min_group_size=2,sampling_points=10):
r"""Compute the mean group size distribution for a Flockwork-P system with varying rates.

Parameters
----------
tau : numpy.ndarray of float
durations for which to evaluate the probability density
max_group_size : int
until which group size the life time distribution should be computed
N : int
Number of nodes
k_initial : float
initial mean degree
t : numpy.ndarray of float
time points at which :math:\alpha(t) and
:math:\beta(t) change
alpha : numpy.ndarray of float
active reconnection rate associated with
the time points in t
beta : numpy.ndarray of float
active disconnection rate associated with
the time points in t
tmax : float
final time
min_group_size : int, default : 2
min group size the life time distribution should be computed for
sampling_points : int, default : 10
how many points to sample in between two time points in t

Returns
-------
P_taus : list of numpy.ndarray
list of Mean probability density of values at tau (contact duration)
"""

# estimate mean degree from integrating ODE
new_t, k = flockwork_P_mean_degree_for_varying_alpha_beta(N,k_initial,t,alpha,beta,tmax,sampling_points)

# 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_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 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*a_*m*(N-m)/(N-1.0)
y = lambda_m * np.exp(-lambda_m*tau)
distros[m-min_group_size].append(y)

mean_distros = []
for m in range(min_group_size, max_group_size+1):
dist = np.array(distros[m-min_group_size])
mean_distros.append( np.trapz(dist, x=new_t, axis=0) / (new_t[-1] - new_t) );

return mean_distros

[docs]def flockwork_P_contact_time_distributions_for_varying_alpha_beta(tau,N,k_initial,t,alpha,beta,tmax,sampling_points=10):
r"""Compute the mean group size distribution for a Flockwork-P system with varying rates.

Parameters
----------
tau : numpy.ndarray of float
durations for which to evaluate the probability density
N : int
Number of nodes
k_initial : float
initial mean degree
t : numpy.ndarray of float
time points at which :math:\alpha(t) and
:math:\beta(t) change
alpha : numpy.ndarray of float
active reconnection rate associated with
the time points in t
beta : numpy.ndarray of float
active disconnection rate associated with
the time points in t
tmax : float
final time
sampling_points : int, default : 10
how many points to sample in between two time points in t

Returns
-------
P_tau_c : numpy.array
Mean probability density of values at tau (contact duration)
P_tau_ic : numpy.array
Mean probability density of values at tau (inter-contact duration)
"""

# estimate mean degree from integrating ODE
new_t, k = flockwork_P_mean_degree_for_varying_alpha_beta(N,k_initial,t,alpha,beta,tmax,sampling_points)

# 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_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 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*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)
distro_ic.append(this_distro_ic)

# compute the mean group size distribution as a time integral over the
# group size distribution
distro_c = np.array(distro_c)
distro_ic = np.array(distro_ic)
mean_distro_c = np.trapz(distro_c, x=new_t, axis=0) / (new_t[-1] - new_t)
mean_distro_ic = np.trapz(distro_ic, x=new_t, axis=0) / (new_t[-1] - new_t)

return mean_distro_c, mean_distro_ic

[docs]def flockwork_P_mean_degree_and_group_size_distribution_for_varying_alpha_beta(N,k_initial,t,alpha,beta,tmax,sampling_points=10):
r"""Compute the mean group size distribution for a Flockwork-P system with varying rates.

Parameters
----------
N : int
Number of nodes
k_initial : float
initial mean degree
t : numpy.ndarray of float
time points at which :math:\alpha(t) and
:math:\beta(t) change
alpha : numpy.ndarray of float
active reconnection rate associated with
the time points in t
beta : numpy.ndarray of float
active disconnection rate associated with
the time points in t
tmax : float
final time
sampling_points : int, default : 10
how many points to sample in between two time points in t

Returns
-------
P_g : numpy.array
An array of times at which the mean degree was evaluated
P_k : numpy.array
An array of mean degree values corresponding to the times in t.
"""

# estimate mean degree from integrating ODE
new_t, k = flockwork_P_mean_degree_for_varying_alpha_beta(N,k_initial,t,alpha,beta,tmax,sampling_points=10)

# from equilibrium assumption k = P/(1-P) compute adjusted P
new_P = k / (k+1)

distro_groups = []
distro_degrees = []

# for every time point and adjusted P, compute the equilibrium group size distribution
for P_ in new_P:
this_distro_groups = flockwork_P_equilibrium_group_size_distribution(N, P_)
this_distro_degrees = degree_distribution(N, P_)
distro_groups.append(this_distro_groups[1:])
distro_degrees.append(this_distro_degrees)

# compute the mean group size distribution as a time integral over the
# group size distribution
distro_groups = np.array(distro_groups)
distro_degrees = np.array(distro_degrees)
mean_distro_groups = np.trapz(distro_groups, x=new_t, axis=0) / (new_t[-1] - new_t)
mean_distro_degrees = np.trapz(distro_degrees, x=new_t, axis=0) / (new_t[-1] - new_t)

return mean_distro_groups, mean_distro_degrees

[docs]def estimated_mean_group_size_distribution(temporal_network):
"""Compute the mean group size distribution for a temporal network under the assumption
that it can be described reasonably by a flockwork-P model.

Parameters
----------
temporal_network : :mod:edge_changes or :mod:edge_lists
A temporal network.

Returns
-------
mean_distribution : numpy.array
The average group size distribution of the temporal network which is closer to
to the _true_ group size distribution than measuring over the binned system.
The result is an array of length N with its i-th entry containing the mean number of
groups of size m = i + 1.
"""

new_t, k = tc.mean_degree(temporal_network)
N = temporal_network.N

# from equilibrium assumption k = P/(1-P) compute adjusted P
new_P = k / (k+1)

distro = []

# for every time point and adjusted P, compute the equilibrium group size distribution
for P_ in new_P:
this_distro = flockwork_P_equilibrium_group_size_distribution(N, P_)
distro.append(this_distro[1:])

# compute the mean group size distribution as a time integral over the
# group size distribution
distro = np.array(distro)
mean_distro = np.trapz(distro, x=new_t, axis=0) / (new_t[-1] - new_t)

return mean_distro

[docs]def flockwork_P_mean_group_size_distribution_from_mean_degree_distribution(flockwork_P_params, dk, N=None):
r"""Compute the mean group size distribution for a Flockwork-P system with varying rates from the
mean degree distribution which is fitted as :math:\left\langle k\right\rangle^{-\alpha},
hence this returns

.. math::

\left\langle N_m \right\rangle = \int dk P(k)  \times N_m( k/(k+1) )

Parameters
----------
flockwork_P_params : :obj:dict
Contains all parameters necessary for a Flockwork-P simulation, especially the
time-dependent rewiring rate and time-dependent reconnection probability
dk : float
resolution of <k>-space for the solution of the integral.
N : int, default : None
If given, compute everything for N nodes, where N is different from
N in flockwork_P_params.

Returns
-------
mean_distribution : numpy.array
An array of length N with its i-th entry containing the mean number of
groups of size m = i + 1.
"""

if N is None:
N = flockwork_P_params['N']

# estimate mean degree from integrating ODE
new_t, k = flockwork_P_mean_degree_for_varying_rates(flockwork_P_params, N)

kmin = 2.0 / N
ind = np.where(k >= kmin)
new_t = new_t[ind]
k = k[ind]

alpha, err, xmin = fit_power_law_clauset(k)
kmin = k.min()
kmax = k.max()

norm = 1/(1-alpha) * (kmax**(-alpha+1) - kmin**(-alpha+1))

def dist(k_): return k_**(-alpha) / norm

k = np.linspace(kmin, kmax, int((kmax-kmin) / dk) + 1)

# from equilibrium assumption k = P/(1-P) compute adjusted P
new_P = k / (k+1)

distro = []

# for every time point and adjusted P, compute the equilibrium group size distribution
for P_ in new_P:
this_distro = flockwork_P_equilibrium_group_size_distribution(N, P_)
distro.append(this_distro[1:])

distro = np.array(distro)

mean_distro = simps(dist(k)[:, None] * distro, x=k, axis=0)

return mean_distro

[docs]def naive_varying_rate_flockwork_simulation(N, t, reconnection_rates, disconnection_rates, tmax):
r"""
Do a naive simulation of a Flockwork systems where the reconnection and disconnection rate
vary over time as step functions. It is called naive because the rate change will not
be considered when evaluating the inter-event time at time points when the rates change.
I.e. at a rate change at time t, the new Flockwork model will be initiated as if the
last event happened at time t. This error is neglibile if the upper bounded mean inter-event time
at rate changes is
:math:\tau=(N(\alpha_{\mathrm{min}}+\beta_{\mathrm{min}})^{-1}\ll \Deltat_\mathrm{min}
where :math:\Deltat_\mathrm{min} is the minimal time between two rate changes.

Parameters
----------
N : int
number of nodes
t : numpy.ndarray of float
time points at which :math:\alpha and :math\beta change
reconnection_rates : numpy.ndarray of float
values of :math:\alpha for the corresponding time values in t.
disconnection_rates : numpy.ndarray of float
values of :math:\beta for the corresponding time values in t.
tmax : float
time at which the simulation is supposed to end.

Returns
-------
edge_changes : :class:_tacoma.edge_changes
The simulated Flockwork instance.
"""

if len(t) != len(reconnection_rates) or len(reconnection_rates) != len(disconnection_rates):
raise ValueError('The arrays t, reconnection_rates and disconnection_rates must have the same length' )

assert(t[-1]<tmax)
t = np.append(t,tmax)

all_networks = []
it = 0

result = None
for a_, b_ in zip(reconnection_rates, disconnection_rates):

dt = t[it+1] - t[it]
if a_ == 0.0 and b_ == 0.0:
P = 0
else:
P = a_ / (a_+b_)
g = (a_+b_)

if it == 0:
initial_edges = flockwork_P_equilibrium_configuration(N,P)

FW = tc.FlockworkPModel(initial_edges, N, g, P, save_temporal_network = True)
FW.simulate(dt)
this_result = FW.edge_changes

if it < len(reconnection_rates) - 1:
initial_edges = FW.get_current_edgelist()

all_networks.append(this_result)

it += 1

result = tc.concatenate(all_networks)

return result

[docs]def flockwork_P(N, P, t_run_total, rewiring_rate = 1.0, initial_edges = None, seed = 0, return_edge_changes_with_histograms=False):
r"""
Simulate a flockwork P-model where the disconnection rate per
node is :math:\gamma and reconnection probability is :math:P.
:func:tacoma.flockwork.flockwork_P_equilibrium_configuration
or just pass initial_edges = None to this function.
The total event rate is :math:N\gamma such that in
one unit of time there's an expected number of
:math:N events.

Parameters
----------
N : int
number of nodes
P : float
The reconnection probability. Has to be :math:0\leq P\leq 1
t_run_total : float
The total run time in units of :math:\gamma^{-1}.
rewiring_rate : float, default : 1.0
Event rate :math:\gamma per node per unit of time.
initial_edges : list of tuple of int, default : None
The initial state of the network as an edge list.
If None is provided, the initial state will be taken
from an equilibrium configuration generated with
:func:tacoma.flockwork.flockwork_P_equilibrium_configuration
seed : int, default : 0
The random seed.
return_edge_changes_with_histograms : bool, default : False
Instead of the converted :class:_tacoma.edge_changes,
return the original instance of
:class:_tacoma.edge_changes_with_histograms.

Returns
-------
:class:_tacoma.edge_changes
The simulated network. if return_edge_changes_with_histograms is True,
returns an instance of :class:_tacoma.edge_changes_with_histograms instead.

"""

if initial_edges is None:
initial_edges = flockwork_P_equilibrium_configuration(N, P)

FW = tc.FlockworkPModel(initial_edges, N, rewiring_rate, P, save_temporal_network = True)
FW.simulate(t_run_total)
this_result = FW.edge_changes

return this_result

[docs]def degree_distribution(N,P):
"""Get the equilibrium degree distribution of a Flockwork-P model
given node number N and probability to reconnect P.

Parameters
----------
N : int
Number of nodes
P : float
Probability to reconnect

Returns
-------
numpy.ndarray
Degree distribution of this configuration. The :math:k-th entry
of this array is the probability that a node has degree :math:k.
"""

C_m = flockwork_P_equilibrium_group_size_distribution(N, P)
P_k = np.arange(1,N+1) * C_m[1:] / N

return P_k

[docs]def degree_moment(N,P,m):
r"""Get the :math:m-th moment of the degree of an Flockwork-P model
equilibrium configuration
given node number N and probability to reconnect P.

Parameters
----------
N : int
Number of nodes
P : float
Probability to reconnect

Returns
-------
<k^m> : float
The :math:m-th moment of the degree distribution.
"""

P_k = degree_distribution(N, P)

return (np.arange(N)**m).dot(P_k)

[docs]def mean_degree(N,P,m):
r"""
Get the exact theoretical mean degree :math:\left< k\right>
for a Flockwork-P model.

Parameters
----------
N : int
Number of nodes
P : float
Probability to reconnect

Returns
-------
<k> : float
The mean degree.
"""

return degree_moment(N, P, 1)

[docs]def convert_to_edge_activity_parameters_plus_minus(N, P, gamma=1.0):
r"""
Convert the Flockwork-P parameters math:P and :math:\gamma
to the corresponding parameters in the edge activity model
:math:\omega^{+} and :math:\omega^{-}.

Parameters
----------
N : int
Number of nodes
P : float
Probability to reconnect
gamma : float, default = 1
Rewiring rate.

Returns
-------
omega_plus : float
The rate with which inactive edges are activated.
omega_minus : float
The rate with which active edges are deactivated.
"""

k = degree_moment(N, P, 1)
k2 = degree_moment(N, P, 2)

omega_minus = 2*gamma*(1-P/(N-1)*k2/k)
rho = k / (N-1)
omega = omega_minus * rho
omega_plus = omega / (1-rho)

return omega_plus, omega_minus

[docs]def convert_to_edge_activity_parameters(N, P, gamma=1.0):
r"""
Convert the Flockwork-P parameters math:P and :math:\gamma
to the corresponding parameters in the edge activity model
:math:\rho and :math:\omega.

Parameters
----------
N : int
Number of nodes
P : float
Probability to reconnect
gamma : float, default = 1
Rewiring rate.

Returns
-------
rho : float
The network density.
omega : float
The rate with which a link is either switched on or off.
"""

k = degree_moment(N, P, 1)
k2 = degree_moment(N, P, 2)

if k == 0.0:
rho = 0.0
omega = 2*gamma
else:

omega_minus = 2*gamma*(1-P/(N-1)*k2/k)
rho = k / (N-1)
omega = omega_minus * rho

return rho, omega

[docs]def convert_to_varying_edge_activity_parameters(N, reconnection_rates, disconnection_rates):
r"""
Convert the Flockwork-P parameters math:\alpha and :math:\beta
to the corresponding parameters in the edge activity model
:math:\rho and :math:\omega.

Parameters
----------
N : int
number of nodes
t : numpy.ndarray of float
time points at which :math:\alpha and :math\beta change
reconnection_rates : numpy.ndarray of float
values of :math:\alpha for the corresponding time values in t.
disconnection_rates : numpy.ndarray of float
values of :math:\beta for the corresponding time values in t.
tmax : float
time at which the simulation is supposed to end.

Returns
-------
rho : float
The network density.
omega : float
The rate with which a link is either switched on or off.
"""

rho = []
omega = []

if len(reconnection_rates) != len(disconnection_rates):
raise ValueError('The arrays reconnection_rates and disconnection_rates must have the same length' )

all_networks = []

for a_, b_ in zip(reconnection_rates, disconnection_rates):

if a_ == 0.0 and b_ == 0.0:
P = 0
else:
P = a_ / (a_+b_)
g = (a_+b_)

rho_, omega_ = convert_to_edge_activity_parameters(N, P, g)

rho.append(rho_)
omega.append(omega_)

return np.array(rho), np.array(omega)

if __name__ == "__main__":

N = 10
P = 0.5

dist = flockwork_P_equilibrium_group_size_distribution(N, P)

print(dist, sum([m * h for m, h in enumerate(dist)]))