analysis

Angle calculation

---

source

dataset_angles

 dataset_angles (trajs:list)

Given a set of trajectories, calculate all angles between displacements

	Type	Details
trajs	list	set of trajectories from which to calculate angles
Returns	list	list of angles between displacements

---

source

get_angle

 get_angle (a:tuple, b:tuple, c:tuple)

Calculates the angle between the segments generate by three points

	Type	Details
a	tuple	2d position point A
b	tuple	2d position point B
c	tuple	2d position point C
Returns	tuple	angle between segments AB and BC points

MSD based analysis

---

source

msd_analysis

 msd_analysis ()

Contains mean squared displacement (MSD) based methods to analyze trajectories.

---

source

msd_analysis.tamsd

 msd_analysis.tamsd (trajs:numpy.ndarray, t_lags:numpy.ndarray, dim=1)

Calculates the time average mean squared displacement (TA-MSD) of a trajectory at various time lags,

	Type	Default	Details
trajs	ndarray		Set of trajectories of dimenions NxTxD (N: number of trajectories, T: lenght, D: dimension)
t_lags	ndarray		Time lags used for the TA-MSD
dim	int	1	Dimension of the trajectories (currently only 1 and 2 supported)
Returns	np.array		TA-MSD of each trayectory / t_lag

---

source

msd_analysis.get_diff_coeff

 msd_analysis.get_diff_coeff (trajs:numpy.ndarray, t_lags:list=None)

Calculates the diffusion coefficient of a trajectory by means of the linear fitting of the TA-MSD.

	Type	Default	Details
trajs	ndarray
t_lags	list	None	Time lags used for the TA-MSD.
Returns	np.array		Diffusion coefficient of the given trajectory.

Here we show an example for the calculation of the diffusion coefficient of a 2D Brownian motion trajectory. We create trajectories from displacements of variance \(\sigma =1\), which results in a diffusion coefficient \(D=0.5\).

pos = np.cumsum(np.random.randn(1070, 100, 2), axis = 1)
D = msd_analysis().get_diff_coeff(pos)

plt.hist(D, bins = 30);
plt.axvline(np.mean(D), c = 'k', label = f'Mean of predictions = {np.round(np.mean(D), 2)}')                                                                  
plt.axvline(0.5, c = 'r', ls = '--', label = 'Expected')
plt.legend()

<matplotlib.legend.Legend>

---

source

msd_analysis.get_exponent

 msd_analysis.get_exponent (trajs:numpy.ndarray, t_lags:list=None)

Calculates the diffusion coefficient of a trajectory by means of the linear fitting of the TA-MSD.

	Type	Default	Details
trajs	ndarray
t_lags	list	None	Time lags used for the TA-MSD.
Returns	np.array		Diffusion coefficient of the given trajectory.

To showcase this function, we generate fractional brownian motion trajectories with \(\alpha = 0.5\) and calculate their exponent:

trajs, _ = models_phenom().single_state(N = 1001, T = 100, alphas = 0.5, dim = 2)
alpha = msd_analysis().get_exponent(trajs.transpose(1,0,2))

plt.hist(alpha, bins = 30);
plt.axvline(np.mean(alpha), c = 'k', label = f'Mean of predictions = {np.round(np.mean(alpha), 2)}')
plt.axvline(0.5, c = 'r', ls = '--', label = 'Expected')
plt.legend()

(1001, 100, 2)

<matplotlib.legend.Legend>

Velocity Autocorrelation Function (VACF)

---

source

vacf

 vacf (trajs, delta_t:int|list|numpy.ndarray=1,
       taus:bool|list|numpy.ndarray=None)

Calculates the velocity autocorrelation function for the given set of trajectories.

	Type	Default	Details
trajs	np.array		NxT matrix containing N trajectories of length T.
delta_t	int | list | numpy.ndarray	1	If not None, the vacf is calculated in the demanded time lags.
taus	bool | list | numpy.ndarray	None	Time windows at wich the vacf is calculated.
Returns	np.array		VACF of the given trajectories and the given time windows.

We show here an example of the VACF for FBM trajectories at various time lages, showing that they all coincide (as expected for this diffusion model).

deltats = np.arange(1, 5).tolist()
taus = np.arange(0, 100)
trajs, _ = models_phenom().single_state(N = 200, T = 200, alphas = 0.5)
trajs = trajs.transpose(1, 0, 2)[:,:,0]

for deltat in deltats:  
    v = vacf(trajs, deltat, taus)        
    plt.plot(taus/deltat, v.flatten(), 'o-', alpha = 0.4)    
plt.xlim(-1, 10)
plt.ylabel('VACF'); plt.xlabel(r'$\tau / \delta$')

Text(0.5, 0, '$\\tau / \\delta$')

Convex hull analysis

---

source

CH_changepoints

 CH_changepoints (trajs, tau:int=10, metric:{'area','volume'}='volume')

Computes the changes points a multistate trajectory based on the Convex Hull approach proposed in PRE 96 (022144), 2017.

	Type	Default	Details
trajs	np.array		NxT matrix containing N trajectories of length T.
tau	int	10	Time window over which the CH is calculated.
metric	{‘area’, ‘volume’}	volume	Calculate change points w.r.t. area or volume of CH.
Returns	list		Change points of the given trajectory.

We showcase the use of the convex hull in a Brownian motion trajectory with two distinct diffusion coefficients, one 10 times the other:

# Generate trajectories and plot change points
T = 100; on = 40; off = 60; tau = 5
traj = np.random.randn(T, 2)
traj[on:off, :] = traj[on:off, :]*10
traj = traj.cumsum(0)
plt.axvline(on-tau, c = 'k')
plt.axvline(off-tau, c = 'k', label = 'True change points')

# Calculate variable Sd 
Sd = np.zeros(traj.shape[0]-2*tau)
for k in range(traj.shape[0]-2*tau):       
    Sd[k] = ConvexHull(traj[k:(k+2*tau)], ).volume  

# Compute change points both with volume and area
CPs = CH_changepoints([traj], tau = tau)[0].flatten()-tau
CPs_a = CH_changepoints([traj], tau = tau, metric = 'area')[0].flatten()-tau


# Plot everything
label_cp = 'CH Volume'
for cp in CPs:
    plt.axvline(cp, c = 'g', alpha = 0.8, ls = '--', label = label_cp)
    label_cp = ''
label_cp = 'CH Area'    
for cp in CPs_a:
    plt.axvline(cp, alpha = 0.8, ls = '--', c = 'orange', label = label_cp)
    label_cp = ''
    
plt.plot(Sd, '-o')
plt.axhline(Sd.mean(), label = 'CH Volume mean', c = 'g',)
plt.legend()
plt.xlabel('n'); plt.ylabel(r'$S_d(n)$')

Text(0, 0.5, '$S_d(n)$')

Cramér-Rao lower bound

Expressions of the CRLB for the estimation of the diffusion coefficient for trajectories without noise, for dimension 1 and 2, at different trajectory lengths

---

source

CRLB_D

 CRLB_D (T:int, dim:int=1)

Calculates the bound for S(D)/D, i.e. ratio between the standard deviation and the expected value of D This holds for x->0 only (i.e. no noise)! See PRE 85, 061916 (2012) for full equation.

	Type	Default	Details
T	int		Length of the trajectory
dim	int	1	Dimension of the trajectoy
Returns	float		Cramér-Rao bound

Here we show how for Brownian motion trajectories, the MSD with is optimal and converges to the CRLB.

#eval: false
Ts = np.arange(10, 100) # Trajectory lengths
Ds = np.arange(0.5, 2.1, 0.5) # We will test on different diffusion coefficients

# Saving the standard deviations
std_ds = np.zeros((len(Ds), len(Ts)))

for idxD, D in enumerate(tqdm(Ds)):
    
    # Generate the BM trajectories with given parameters
    trajs = (np.sqrt(2*D)*np.random.randn(int(1e4), Ts[-1], 1)).cumsum(1) 

    for idxL, T in enumerate(tqdm(Ts)):
        
        std_ds[idxD, idxL] = np.array(msd_analysis().get_diff_coeff(trajs = trajs[:, :T], t_lags = [1,2])).std()

fig, axs = plt.subplots(1,4, figsize = (4*3,3), tight_layout = True)

for idx, (ax, std_d, D) in enumerate(zip(axs.flatten(), std_ds, Ds)):
    ax.plot(Ts, std_d/D, label = r'MSD$_{1-2}$')
    ax.scatter(Ts[::3], CRLB_D(Ts)[::3], label = 'Cramer-Rao bound', alpha = 0.9, facecolor = 'w', edgecolor = 'C1')
    ax.set_title(fr'$D = {D}$')
axs[0].legend()
plt.setp(axs[0], ylabel = 'std(D) / D')
plt.setp(axs, xlabel = 'Traj. length');

Here is now the same calculation for 2D trajectories, showing that we get the same convergence!

fig, axs = plt.subplots(1,4, figsize = (4*3,3), tight_layout = True)

for idx, (ax, std_d, D) in enumerate(zip(axs.flatten(), std_ds, Ds)):
    ax.plot(Ts, std_d/D, label = r'MSD$_{1-2}$', c = 'C2')
    ax.scatter(Ts[::3], CRLB_D(Ts, dim = 2)[::3], label = 'Cramer-Rao bound', alpha = 0.9, facecolor = 'w', edgecolor = 'C3')
    ax.set_title(fr'$D = {D}$')
axs[0].legend()
plt.setp(axs[0], ylabel = 'std(D) / D')
plt.setp(axs, xlabel = 'Traj. length');