<jupyter_start><jupyter_text>Symbolic System Create a symbolic three-state system:<jupyter_code>import markoviandynamics as md
sym_system = md.SymbolicDiscreteSystem(3)<jupyter_output><empty_output><jupyter_text>Get the symbolic equilibrium distribution:<jupyter_code>sym_system.equilibrium()<jupyter_output><empty_output><jupyter_text>Create a symbolic three-state system with potential energy barriers:<jupyter_code>sym_system = md.SymbolicDiscreteSystemArrhenius(3)<jupyter_output><empty_output><jupyter_text>It's the same object as the previous one, only with additional symbolic barriers:<jupyter_code>sym_system.B_ij<jupyter_output><empty_output><jupyter_text>We can assing values to the free parameters in the equilibrium distribution:<jupyter_code>sym_system.equilibrium(energies=[0, 0.1, 1])
sym_system.equilibrium(energies=[0, 0.1, 1], temperature=1.5)<jupyter_output><empty_output><jupyter_text>and create multiple equilibrium points by assigning temperature sequence:<jupyter_code>import numpy as np
temperature_range = np.linspace(0.01, 10, 300)
equilibrium_line = sym_system.equilibrium([0, 0.1, 1], temperature_range)

equilibrium_line.shape<jupyter_output><empty_output><jupyter_text>Symbolic rate matrix Create a symbolic rate matrix with Arrhenius process transitions:<jupyter_code>sym_rate_matrix = md.SymbolicRateMatrixArrhenius(3)

sym_rate_matrix<jupyter_output><empty_output><jupyter_text>Energies and barriers can be substituted at once:<jupyter_code>energies = [0, 0.1, 1]
barriers = [[0, 0.11, 1.1], 
            [0.11, 0, 10], 
            [1.1, 10, 0]]

sym_rate_matrix.subs_symbols(energies, barriers)
sym_rate_matrix.subs_symbols(energies, barriers, temperature=2.5)<jupyter_output><empty_output><jupyter_text>A symbolic rate matrix can be also lambdified (transform to lambda function):<jupyter_code>rate_matrix_lambdified = sym_rate_matrix.lambdify()<jupyter_output><empty_output><jupyter_text>The parameters of this function are the free symbols in the rate matrix:<jupyter_code>rate_matrix_lambdified.__code__.co_varnames<jupyter_output><empty_output><jupyter_text>They are positioned in ascending order. First the temperature, then the energies and the barriers. Sequence of rate matrices can be created by calling this function with a sequence for each parameter.  Dynamics We start by computing an initial probability distribution by assigning the energies and temperature:<jupyter_code>p_initial = sym_system.equilibrium(energies, 0.5)

p_initial<jupyter_output><empty_output><jupyter_text>Trajectory - evolve by a fixed rate matrix Compute the rate matrix by substituting free symbols:<jupyter_code>rate_matrix = md.rate_matrix_arrhenius(energies, barriers, 1.2)

rate_matrix<jupyter_output><empty_output><jupyter_text>Create trajectory of probability distributions in time:<jupyter_code>import numpy as np

# Create time sequence
t_range = np.linspace(0, 5, 100)

trajectory = md.evolve(p_initial, rate_matrix, t_range)

trajectory.shape
import matplotlib.pyplot as plt
%matplotlib inline

for i in [0, 1, 2]:
    plt.plot(t_range, trajectory[i,0,:], label='$p_{}(t)$'.format(i + 1))
plt.xlabel('$t$')
plt.legend()<jupyter_output><empty_output><jupyter_text>Trajectory - evolve by a time-dependent rate matrix Create a temperature sequence in time:<jupyter_code>temperature_time = 1.4 + np.sin(4. * t_range)<jupyter_output><empty_output><jupyter_text>Create a rate matrix as a function of the temperature sequence:<jupyter_code># Array of stacked rate matrices that corresponds to ``temperature_time``
rate_matrix_time = md.rate_matrix_arrhenius(energies, barriers, temperature_time)

rate_matrix_time.shape
crazy_trajectory = md.evolve(p_initial, rate_matrix_time, t_range)

crazy_trajectory.shape
for i in [0, 1, 2]:
    plt.plot(t_range, crazy_trajectory[i,0,:], label='$p_{}(t)$'.format(i + 1))
plt.xlabel('$t$')
plt.legend()<jupyter_output><empty_output><jupyter_text>Diagonalize the rate matrix Calculate the eigenvalues, left and right eigenvectors:<jupyter_code>U, eigenvalues, V = md.eigensystem(rate_matrix)

U.shape, eigenvalues.shape, V.shape<jupyter_output><empty_output><jupyter_text>The eigenvalues are in descending order (the eigenvectors are ordered accordingly):<jupyter_code>eigenvalues<jupyter_output><empty_output><jupyter_text>We can also compute the eigensystem for multiple rate matrices at once (or evolution of a rate matrix, i.e., `rate_matrix_time`):<jupyter_code>U, eigenvalues, V = md.eigensystem(rate_matrix_time)

U.shape, eigenvalues.shape, V.shape<jupyter_output><empty_output><jupyter_text>Decompose to rate matrix eigenvectors A probability distribution, in general, can be decomposed to the right eigenvectors of the rate matrix:$$\left|p\right\rangle = a_1\left|v_1\right\rangle + a_2\left|v_2\right\rangle + a_3\left|v_3\right\rangle$$where $a_i$ is the coefficient of the i'th right eigenvector $\left|v_i\right\rangle$. A rate matrix that satisfies detailed balance has its first eigenvector as the equilibrium distribution $\left|\pi\right\rangle$. Therefore, *markovian-dynamics* normalizes $a_1$ to $1$ and decompose a probability distribution to $$\left|p\right\rangle = \left|\pi\right\rangle + a_2\left|v_2\right\rangle + a_3\left|v_3\right\rangle$$ Decompose ``p_initial``:<jupyter_code>md.decompose(p_initial, rate_matrix)<jupyter_output><empty_output><jupyter_text>We can decompose also multiple points and/or by multiple rate matrices. For example, decompose multiple points:<jupyter_code>first_decomposition = md.decompose(equilibrium_line, rate_matrix)

first_decomposition.shape
for i in [0, 1, 2]:
    plt.plot(temperature_range, first_decomposition[i,:], label='$a_{}(T)$'.format(i + 1))
plt.xlabel('$T$')
plt.legend()<jupyter_output><empty_output><jupyter_text>or decompose a trajectory:<jupyter_code>second_decomposition = md.decompose(trajectory, rate_matrix)

second_decomposition.shape
for i in [0, 1, 2]:
    plt.plot(t_range, second_decomposition[i,0,:], label='$a_{}(t)$'.format(i + 1))
plt.xlabel('$t$')
plt.legend()<jupyter_output><empty_output><jupyter_text>Decompose single point using multiple rate matrices:<jupyter_code>third_decomposition = md.decompose(p_initial, rate_matrix_time)

third_decomposition.shape
for i in [0, 1, 2]:
    plt.plot(t_range, third_decomposition[i,0,:], label='$a_{}(t)$'.format(i + 1))
plt.legend()<jupyter_output><empty_output><jupyter_text>Decompose, for every time $t$, the corresponding point $\left|p(t)\right\rangle$ using the temporal rate matrix $R(t)$<jupyter_code>fourth_decomposition = md.decompose(trajectory, rate_matrix_time)

fourth_decomposition.shape
for i in [0, 1, 2]:
    plt.plot(t_range, fourth_decomposition[i,0,:], label='$a_{}(t)$'.format(i + 1))
plt.legend()<jupyter_output><empty_output><jupyter_text>Plotting the 2D probability simplex for three-state system The probability space of a three-state system is a three dimensional space. However, the normalization constraint $\sum_{i}p_i=1$ together with $0 < p_i \le 1$, form a 2D triangular plane in which all of the possible probability points reside. We'll start by importing the plotting module:<jupyter_code>import markoviandynamics.plotting.plotting2d as plt2d

# Use latex rendering
plt2d.latex()<jupyter_output><empty_output><jupyter_text>Plot the probability plane:<jupyter_code>plt2d.figure(figsize=(7, 5.5))
plt2d.equilibrium_line(equilibrium_line)
plt2d.legend()<jupyter_output><empty_output><jupyter_text>We can plot many objects on the probability plane, such as trajectories, points, and eigenvectors of the rate matrix:<jupyter_code># Final equilibrium point
p_final = sym_system.equilibrium(energies, 1.2)
plt2d.figure(focus=True, figsize=(7, 5.5))
plt2d.equilibrium_line(equilibrium_line)

# Plot trajectory
plt2d.plot(trajectory, c='r', label=r'$\left|p(t)\right>$')

# Initial & final points
plt2d.point(p_initial, c='k', label=r'$\left|p_0\right>$')
plt2d.point(p_final, c='r', label=r'$\left|\pi\right>$')

# Eigenvectors
plt2d.eigenvectors(md.eigensystem(rate_matrix), kwargs_arrow={'zorder': 1})

plt2d.legend()<jupyter_output><empty_output><jupyter_text>Plot multiple trajectories at once:<jupyter_code># Create temperature sequence
temperature_range = np.logspace(np.log10(0.01), np.log10(10), 50)

# Create the equilibrium line points
equilibrium_line = sym_system.equilibrium(energies, temperature_range)

# Create a trajectory for every point on ``equilibrium_line``
equilibrium_line_trajectory = md.evolve(equilibrium_line, rate_matrix, t_range)
plt2d.figure(focus=True, figsize=(7, 5))
plt2d.equilibrium_line(equilibrium_line)
plt2d.plot(equilibrium_line_trajectory, c='g', alpha=0.2)
plt2d.point(p_final, c='r', label=r'$\left|\pi\right>$')
plt2d.legend()
# Create a trajectory for every point on ``equilibrium_line``
equilibrium_line_crazy_trajectory = md.evolve(equilibrium_line, rate_matrix_time, t_range)
plt2d.figure(focus=True, figsize=(7, 5))
plt2d.equilibrium_line(equilibrium_line)
plt2d.plot(equilibrium_line_crazy_trajectory, c='r', alpha=0.1)
plt2d.text(p_final, r'Text $\alpha$', delta_x=0.05)
plt2d.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1143956043956044, 'fraction_numerical': 0.01769230769230769, 'mean_word_length': 8.073778664007977, 'pattern_counts': {'":': 0, '<': 110, '<?xml version=': 0, '>': 113, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '7224', 'license': 'MIT', 'n_tokens_mistral': 3323, 'n_tokens_neox': 3125, 'n_words': 691, 'path': 'notebooks/usage.ipynb', 'repo_name': 'amitg7/markovian-dynamics', 'size': '9100'}