Help with Methodology (vectors' angles from Kramers-Moyal coefficients)

[apathanasiadis]

Hello, I came across your package in my effort to replicate the methodology from paper:
https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1000061
where in the Data Analysis section they use the Kramers-Moyal coefficients to get the drift coefficients of the vectors and subsequently the angle of neighbor vectors in the phase space in order to make conclusions about the dynamics of the flow in the state space.
I don't know if you would be interested in helping me a bit with the code. I actually want to do it for a 3d space, but even in 2d i don't quite understand how to get the angles from the output of the km function since i get this bins x bins x bins vector.
I hope I was clear. Thanks a lot in advance and apologies if my issue fall out of the scope of the issues section :)

[fmeirinhos]

I only glanced at the paper and the formula for the angle should actually be arccos((u . v)/(|u| |v|)).

Do you have any minimal example you can share with us? Note that the bins are just the discretization, per dimension, of the underlying real-space -- that is, the space where your vectors are defined.

[LRydin]

Hey @athantas95,
Nice to know that you are interested in this. I skimmed through the paper and, from what I can understand what you have in mind should be possible - and not too much of a nightmare. Let's start with a 2D example because 3D is simply harder to visualise.

You should have access to a 2D time series of your dynamical (stochastic) system. Let's call it `y'. You should just have to use

import kramersmoyal as kmc
# Y is a (2,N) time series
powers = [[0,0], [1,0], [0,1]]
kmc, edges = kmc.km(y, powers=2)

Now you are interested in the kmc[1,...] and kmc[2,...] which correspond to the drift in the first and the first in the second dimension (here in our 2D problem). You will need something like quiver from matplotlib to visualise the results. I always get terribly confused with the way quiver works, but I think something like this should do:

import matplotlib.pyplot as plt
import numpy as np

X, Y = np.meshgrid(edges[0], edges[1])
fig, ax = plt.subplots()
ax.quiver(X, Y, kmc[1,...], kmc[2,...], units='width')

Can you give that a try?

[fmeirinhos]

Building on @LRydin's answer, note that a power [m, n] enters the formula, as in the paper's Data Analysis, as

D_{x^m y^n} = \int (x - x')^m (y - y')^n P(...)

So to get D_x and D_y, you'd need the power [1,0] and [0,1], respectively. Note that these have index 1 and 2 respectively if powers=[[0,0], [1,0], [0,1]]. Hence,

import numpy as np
import kramersmoyal

y = np.random.randn(10000, 2) # Some fake data
D, edges = kramersmoyal.km(y, powers=[[0,0], [1,0], [0,1]], bins=[20,10])

D_x = D[1, ...]
D_y = D[2, ...]

To obtain a proper velocity vector we should stack these

D = np.stack((D_x, D_y))

Note that this "2d-vector" is is defined for all points of the underlying grid {x, y} and hence has shape (2, bins[0], bins[1]). To get the angle, you'd need a correct function

def angle_between(v1, v2):
    norms = np.linalg.norm(v1) * np.linalg.norm(v2)

    # Avoid division by zero and numerical instabilities with clamp
    if norms == 0:
        return 0
    return np.arccos(np.clip(np.dot(v1, v2) / norms, -1.0, 1.0))

So now the paper states that you are computing the angle of this velocity vector with its neighbours. This is hard-coded for 2d but trivial to generalise for 3 dimensions.

def calculate_angles(v):
    _, n_x, n_y = v.shape
    angles = np.zeros((n_x, n_y))

    # Iterate over each point in the lattice
    for x in range(n_x):
        for y in range(n_y):
            v_current = v[:, x, y]

            # Pre-determine the valid neighbor offsets (luckily the underlying grid is box-like)
            neighbor_offsets = [(dx, dy) for dx in [-1, 0, 1] for dy in [-1, 0, 1] 
                                if 0 <= x + dx < n_x and 0 <= y + dy < n_y and not (dx == 0 and dy == 0)]

            # Calculate angles with valid neighbors and SUM THEM
            for dx, dy in neighbor_offsets:
                v_neighbor = v[:, x + dx, y + dy]
                angles[x, y] += angle_between(v_current, v_neighbor)

    return angles

Hence,

theta_xy = calculate_angles(D)

To visualise

import matplotlib.pyplot as plt
X, Y = np.meshgrid(edges[0], edges[1])
plt.pcolormesh(X, Y, theta_xy.T)
plt.show()

To now solve your problem, you'd need just to adapt the powers to 3d, that is powers = [[0,0,0], [1,0,0], [0,1,0], [0,0,1] to obtain (D_x, D_y, D_z) and then generalise calculate_angles to all the different angles in a higher dimensional space

[apathanasiadis]

Hi again! Thanks a lot for the very quick responses and help. I tried the first snippet by @LRydin and I guess now I understand a bit better the output of the km function. I will now try the extended code by @fmeirinhos and let you know. Thanks a ton for your help :D

[fmeirinhos]

It just occurred to me that angle_between may not be 100% suitable since it's just the angle magnitude, defined in [0, pi], and since we are adding angles, it would be better to use a formula valid for [0, 2pi). I found out this formula, which works in 2d

def angle_between(v1, v2):
    angle = np.arctan2(v2[1], v2[0]) - np.arctan2(v1[1], v1[0])
    return np.mod(angle, 2 * np.pi)

Similarly, since we are adding so many angles, it's important the final result is mod 2\pi

def calculate_angles(v):
  # ...
  return np.mod(angles, 2*np.pi)

Please take the script with a rock of salt, there might be other "problems" lingering, we thought we would just give you a head start. Good luck!