diff --git a/Ben_Manuscripts/stochastic_transport/Supporting_Information.pdf b/Ben_Manuscripts/stochastic_transport/Supporting_Information.pdf index 10dead2f..fe2a097c 100644 Binary files a/Ben_Manuscripts/stochastic_transport/Supporting_Information.pdf and b/Ben_Manuscripts/stochastic_transport/Supporting_Information.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/Supporting_Information.tex b/Ben_Manuscripts/stochastic_transport/Supporting_Information.tex index 3c22478c..7e893990 100644 --- a/Ben_Manuscripts/stochastic_transport/Supporting_Information.tex +++ b/Ben_Manuscripts/stochastic_transport/Supporting_Information.tex @@ -364,6 +364,51 @@ more negatively correlated than they should be.}\label{fig:hurst_correction} \end{figure} + \section{Verifying Markovianity}\label{section:markov_validation} + + We verified the Markovianity of our transition matrix, $T$, in two ways. First we + ensured that the process satisfied detailed balance: + \begin{equation} + T_{i,j}P_i(t=\infty) = T_{j,i}P_j(t=\infty) + \end{equation} + where $P$ is the equilibrium distribution of states. This implies that the number + of transitions from state $i$ to $j$ and from state $j$ to $i$ should be equal. Graphical + representations of the count matrices show that this is true in Figure~\ref{fig:counts}. + + Second, we ensured that the transition matrix did not change on coarser time scales. + In Figures~\ref{fig:counts} and~\ref{fig:transitions}, we show that increasing the + length of time between samples does not change the properties of the count or + probability transition matrices. + + \begin{figure} + \centering + \begin{subfigure}{\textwidth} + \includegraphics[width=\textwidth]{URE_counts.pdf} + \caption{Urea}\label{fig:URE_counts} + \end{subfigure} + \begin{subfigure}{\textwidth} + \includegraphics[width=\textwidth]{GCL_counts.pdf} + \caption{Ethylene Glycol}\label{fig:GCL_counts} + \end{subfigure} + \caption{The number of transitions from state $i$ to $j$ and $j$ to $i$ are very close + indicating that our process obeys detailed balance. Detailed balance is conserved for + different sized time steps.}\label{fig:counts} + \end{figure} + + \begin{figure} + \centering + \begin{subfigure}{\textwidth} + \includegraphics[width=\textwidth]{URE_transitions.pdf} + \caption{Urea}\label{fig:URE_transitions} + \end{subfigure} + \begin{subfigure}{\textwidth} + \includegraphics[width=\textwidth]{GCL_transitions.pdf} + \caption{Ethylene Glycol}\label{fig:GCL_transitions} + \end{subfigure} + \caption{As the timestep between observations increases, the probability + transition matrix does not change signficantly.}\label{fig:transitions} + \end{figure} + \section{2-state MSDDM parameters}\label{section:simple_msddm_params} We used a simple 2-state model to illustrate how the self-transition probability diff --git a/Ben_Manuscripts/stochastic_transport/figures/ACH_hop_acf.pdf b/Ben_Manuscripts/stochastic_transport/figures/ACH_hop_acf.pdf new file mode 100644 index 00000000..0f908358 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/ACH_hop_acf.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/GCL_hop_acf.pdf b/Ben_Manuscripts/stochastic_transport/figures/GCL_hop_acf.pdf new file mode 100644 index 00000000..06fe5c7a Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/GCL_hop_acf.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/GCL_msddm.pdf b/Ben_Manuscripts/stochastic_transport/figures/GCL_msddm.pdf new file mode 100644 index 00000000..7f8aabc7 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/GCL_msddm.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/GCL_msddm_zeroH.pdf b/Ben_Manuscripts/stochastic_transport/figures/GCL_msddm_zeroH.pdf new file mode 100644 index 00000000..e06419d1 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/GCL_msddm_zeroH.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/MET_hop_acf.pdf b/Ben_Manuscripts/stochastic_transport/figures/MET_hop_acf.pdf new file mode 100644 index 00000000..82c21c00 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/MET_hop_acf.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/MET_msddm_zeroH.pdf b/Ben_Manuscripts/stochastic_transport/figures/MET_msddm_zeroH.pdf new file mode 100644 index 00000000..ece506c5 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/MET_msddm_zeroH.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/T_sensitivity.pdf b/Ben_Manuscripts/stochastic_transport/figures/T_sensitivity.pdf new file mode 100644 index 00000000..1b8d060e Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/T_sensitivity.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/URE_hop_acf.pdf b/Ben_Manuscripts/stochastic_transport/figures/URE_hop_acf.pdf new file mode 100644 index 00000000..17228788 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/URE_hop_acf.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/URE_msddm.pdf b/Ben_Manuscripts/stochastic_transport/figures/URE_msddm.pdf new file mode 100644 index 00000000..6ea2c8d0 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/URE_msddm.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/URE_msddm_zeroH.pdf b/Ben_Manuscripts/stochastic_transport/figures/URE_msddm_zeroH.pdf new file mode 100644 index 00000000..b85fbb72 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/URE_msddm_zeroH.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/URE_powerlaw.pdf b/Ben_Manuscripts/stochastic_transport/figures/URE_powerlaw.pdf new file mode 100644 index 00000000..655ae832 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/URE_powerlaw.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/URE_trajectories.pdf b/Ben_Manuscripts/stochastic_transport/figures/URE_trajectories.pdf new file mode 100644 index 00000000..8b695ea9 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/URE_trajectories.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison.pdf b/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison.pdf new file mode 100644 index 00000000..9d52104b Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous.pdf b/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous.pdf new file mode 100644 index 00000000..5c279863 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous_GCL.pdf b/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous_GCL.pdf new file mode 100644 index 00000000..1f48e9e0 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous_GCL.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous_URE.pdf b/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous_URE.pdf new file mode 100644 index 00000000..2de9e6aa Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/gaussian_levy_comparison_anomalous_URE.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/hop_acf.pdf b/Ben_Manuscripts/stochastic_transport/figures/hop_acf.pdf new file mode 100644 index 00000000..20f6ef5b Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/hop_acf.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/msddm_acf.pdf b/Ben_Manuscripts/stochastic_transport/figures/msddm_acf.pdf new file mode 100644 index 00000000..cd60bed2 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/msddm_acf.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/figures/state_probabilities.pdf b/Ben_Manuscripts/stochastic_transport/figures/state_probabilities.pdf new file mode 100644 index 00000000..0c8a0f05 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/figures/state_probabilities.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/ACH_equilibration.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/ACH_equilibration.pdf new file mode 100644 index 00000000..3387e02d Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/ACH_equilibration.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_counts.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_counts.pdf new file mode 100644 index 00000000..29866235 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_counts.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_equilibration.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_equilibration.pdf new file mode 100644 index 00000000..7ae0fdcc Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_equilibration.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_transitions.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_transitions.pdf new file mode 100644 index 00000000..728d13cc Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/GCL_transitions.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/MET_equilibration.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/MET_equilibration.pdf new file mode 100644 index 00000000..070fbd21 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/MET_equilibration.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_counts.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_counts.pdf new file mode 100644 index 00000000..f55335a4 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_counts.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_equilibration.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_equilibration.pdf new file mode 100644 index 00000000..6520ebae Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_equilibration.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_transitions.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_transitions.pdf new file mode 100644 index 00000000..611b794d Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/URE_transitions.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/flm_autocovariance.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/flm_autocovariance.pdf new file mode 100644 index 00000000..47e0cc64 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/flm_autocovariance.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/hurst_autocorrelation.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/hurst_autocorrelation.pdf new file mode 100644 index 00000000..e5a78ae4 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/hurst_autocorrelation.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/hurst_correction.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/hurst_correction.pdf new file mode 100644 index 00000000..16498900 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/hurst_correction.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/truncation_correction.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/truncation_correction.pdf new file mode 100644 index 00000000..79c7ca20 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/truncation_correction.pdf differ diff --git a/Ben_Manuscripts/stochastic_transport/supporting_figures/truncation_cutoff.pdf b/Ben_Manuscripts/stochastic_transport/supporting_figures/truncation_cutoff.pdf new file mode 100644 index 00000000..0d0737c2 Binary files /dev/null and b/Ben_Manuscripts/stochastic_transport/supporting_figures/truncation_cutoff.pdf differ diff --git a/LLC_Membranes/analysis/markov_state_dependent_dynamics.py b/LLC_Membranes/analysis/markov_state_dependent_dynamics.py index eeb3d90f..69d1d88f 100755 --- a/LLC_Membranes/analysis/markov_state_dependent_dynamics.py +++ b/LLC_Membranes/analysis/markov_state_dependent_dynamics.py @@ -5,6 +5,9 @@ import numpy as np import mdtraj as md import matplotlib.pyplot as plt +import matplotlib.colors as colors +import matplotlib.cbook as cbook +from mpl_toolkits.axes_grid1 import AxesGrid from LLC_Membranes.analysis import hbonds, coordination_number from LLC_Membranes.llclib import file_rw, physical, topology, timeseries, fitting_functions from LLC_Membranes.timeseries.fractional_levy_motion import FLM @@ -295,14 +298,14 @@ def _determine_state_sequence(self): print("Done!") - def _make_transition_matrix(self, start=1): + def _make_transition_matrix(self, start=1, step=1): """ Estimate a transition matrix based on the state sequence """ self.transition_matrix = np.zeros([self.nstates, self.nstates]) self.count_matrix = np.zeros_like(self.transition_matrix, dtype=int) - for t in range(start, self.nT): # start at frame 1. May need to truncate more as equilibration + for t in range(start, self.nT, step): # start at frame 1. May need to truncate more as equilibration transitioned_from = self.state_sequence[t - 1, :] transitioned_to = self.state_sequence[t, :] for pair in zip(transitioned_from, transitioned_to): @@ -311,6 +314,73 @@ def _make_transition_matrix(self, start=1): # normalize so rows sum to unity self.transition_matrix = (self.count_matrix.T / self.count_matrix.sum(axis=1)).T + def verify_markovinity(self, timesteps=(1, 2, 4), cmap='bwr', dt=0.5, show=False, save=False): + """ Verify that the transition matrix is Markovian + + The condition of detailed balance is satisfied if the following holds: + .. math:: + + T_{i,j}P_j(t=\infty) = T_{j, i}P_j(t=\infty) + + In other words, the off-diagonal entries of the transition matrix are equal. + + A further test is to show that the transition matrix is independent of the time step. + + :param timesteps: timesteps at which to create transition matrix. A timestep of x would add to the transition + matrix every x frames. + + :type timesteps: tuple + """ + + start = 0.5 + end = self.nstates + start + extent = [start, end, start, end] + x_positions = np.linspace(start, end - 1, self.nstates) + y_positions = np.linspace(start, end - 1, self.nstates) + + fig = plt.figure(figsize=(15, 6)) + fig2 = plt.figure(figsize=(15, 6)) + grid = AxesGrid(fig, 111, nrows_ncols=(1, 3), cbar_mode='none', cbar_location='top', cbar_pad=0.2) + grid2 = AxesGrid(fig2, 111, nrows_ncols=(1, 3), cbar_mode='none', cbar_location='top', cbar_pad=0.2) + + for s, ax in enumerate(grid): + + self._make_transition_matrix(step=timesteps[s]) + + T = self.transition_matrix[::-1, :].T + C = self.count_matrix[::-1, :].T + + ax.imshow(T.T, extent=extent, origin='lower', cmap=cmap) + grid2.axes_all[s].imshow(C.T, extent=extent, origin='lower', cmap=cmap, + norm=colors.LogNorm(self.count_matrix.min(), self.count_matrix.max())) + + for yndx, y in enumerate(y_positions): + for xndx, x in enumerate(x_positions): + + label = T[xndx, yndx] + counts = C[xndx, yndx] + ax.text(x + start, y + start, '%.2f' % label, color='black', ha='center', va='center', weight='bold') + grid2.axes_all[s].text(x + start, y + start, '%d' % counts, color='black', ha='center', va='center', weight='bold') + + ax.tick_params(labelsize=14) + ax.set_xticks(np.arange(1, self.nstates + 1)) + ax.set_title('Timestep = %.1f ns' % (dt * timesteps[s]), fontsize=14) + ax.set_yticklabels(np.arange(9, 0, -1)) + + grid2.axes_all[s].tick_params(labelsize=14) + grid2.axes_all[s].set_xticks(np.arange(1, self.nstates + 1)) + grid2.axes_all[s].set_title('Timestep = %.1f ns' % (dt * timesteps[s]), fontsize=14) + grid2.axes_all[s].set_yticklabels(np.arange(9, 0, -1)) + + fig.tight_layout() + fig2.tight_layout() + + if save: + fig.savefig('%s_transitions.pdf' % self.residue) + fig2.savefig('%s_counts.pdf' % self.residue) + if show: + plt.show() + def measure_state_emissions(self, dim=2, fit_function=levy): """ Measure observations as a function of state. @@ -883,8 +953,10 @@ def plot_msd(self, cutoff=0.4, dt=0.5, save=False, savename='msd.pdf', label='MS # print(states.state_sequence) # exit() - # states._make_transition_matrix(1, end=100) - # + states.verify_markovinity() + + #states._make_transition_matrix(1, end=100) + total = states.count_matrix.sum(axis=0) p = total / total.sum() # probability of being in state i at any time t w, v = np.linalg.eig(states.transition_matrix.T) @@ -892,6 +964,7 @@ def plot_msd(self, cutoff=0.4, dt=0.5, save=False, savename='msd.pdf', label='MS print(w[0]) print(v[:, 0] / v[:, 0].sum()) print(p) + plt.show() exit() ############ Cut-off Justification Plot ################## # alpha, mu, sigma = states.fit_params[-1]