diff --git a/docs/api/analysis.rst b/docs/api/analysis.rst new file mode 100644 index 0000000..0d6cacd --- /dev/null +++ b/docs/api/analysis.rst @@ -0,0 +1,18 @@ +================= +Analysis +================= + +.. currentmodule:: folie + + +Overdamped 1D models +========================= +.. autosummary:: + :toctree: generated/ + :template: function.rst + + analysis.free_energy_profile_1d + + analysis.mfpt_1d + + diff --git a/docs/conf.py b/docs/conf.py index be79685..c3e74f1 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -41,11 +41,12 @@ "nbsphinx", "numpydoc", "sphinx_gallery.gen_gallery", - "sphinx_gallery.load_style", + # "sphinx_gallery.load_style", "sphinx.ext.inheritance_diagram", "sphinx.ext.githubpages", "sphinx.ext.linkcode", "sphinx.ext.mathjax", + "sphinxcontrib.bibtex", ] numpydoc_show_class_members = False diff --git a/docs/index.rst b/docs/index.rst index 71082c3..9e291dd 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -83,6 +83,7 @@ Table of contents api/models api/functions api/estimation + api/analysis Indices and tables ================== diff --git a/docs/notebooks/estimation.ipynb b/docs/notebooks/estimation.ipynb index 5d92cd1..f56842a 100644 --- a/docs/notebooks/estimation.ipynb +++ b/docs/notebooks/estimation.ipynb @@ -143,9 +143,11 @@ "id": "550cf6a2", "metadata": {}, "source": [ - "But also obtain the free energy profile. The free energy profile is related to the force and the diffusion from\n", + "But also obtain the free energy profile. The free energy profile $V(x)$ is related to the force $F(x)$ and the diffusion $D(x)$ from\n", "\n", - "$$ F(x) = -D(x) \\nabla V(x) + \\mathrm{div} D(x)$$" + "$$ F(x) = -D(x) \\nabla V(x) + \\mathrm{div} D(x)$$\n", + "\n", + "The relation can then be inverted to obtain the free energy profile." ] }, { @@ -188,6 +190,43 @@ "ax.plot(xfa, pmf)" ] }, + { + "cell_type": "markdown", + "id": "ce0d9688-9d66-4231-92ee-1c4c053e13ec", + "metadata": {}, + "source": [ + "Since there is two well, we can also compute the mean first passage time to go from point x to 0." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "9944b6d8-5317-4ff2-a436-983d7bccffa7", + "metadata": {}, + "outputs": [ + { + "ename": "NameError", + "evalue": "name 'fl' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[1], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m x_mfpt, mfpt \u001b[38;5;241m=\u001b[39m \u001b[43mfl\u001b[49m\u001b[38;5;241m.\u001b[39manalysis\u001b[38;5;241m.\u001b[39mmfpt_1d(model_simu, \u001b[38;5;241m0\u001b[39m, [\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m20.0\u001b[39m, \u001b[38;5;241m50.0\u001b[39m], Npoints\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m500\u001b[39m)\n", + "\u001b[0;31mNameError\u001b[0m: name 'fl' is not defined" + ] + } + ], + "source": [ + "x_mfpt, mfpt = fl.analysis.mfpt_1d(model_simu, 0, [-20.0, 50.0], Npoints=500)\n", + "fig, ax = plt.subplots(1, 1)\n", + "# MFPT plot\n", + "ax.set_title(\"MFPT from x to 0\")\n", + "ax.grid()\n", + "ax.set_xlabel(\"$x$\")\n", + "ax.set_ylabel(\"$MFPT(x,0)$\")\n", + "ax.plot(x_mfpt, mfpt)" + ] + }, { "cell_type": "markdown", "id": "24ed7493", @@ -204,9 +243,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3.8", "language": "python", - "name": "python3" + "name": "python3.8" }, "language_info": { "codemirror_mode": { @@ -218,7 +257,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.8.18" } }, "nbformat": 4, diff --git a/docs/references.bib b/docs/references.bib index 2f47551..974ca1e 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -11,19 +11,38 @@ @article{Palacio-Rodriguez2022 doi = {10.1021/acs.jctc.2c00324}, note ={PMID: 35899416}, -URL = { - +URL = { + https://doi.org/10.1021/acs.jctc.2c00324 - - + + }, -eprint = { - +eprint = { + https://doi.org/10.1021/acs.jctc.2c00324 - - + + + +} } + +@article{Jungblut2016, + title = {Pathways to Self-Organization: {{Crystallization}} via Nucleation and Growth}, + shorttitle = {Pathways to Self-Organization}, + author = {Jungblut, S. and Dellago, C.}, + year = {2016}, + month = aug, + journal = {The European Physical Journal E}, + volume = {39}, + number = {8}, + pages = {77}, + issn = {1292-895X}, + doi = {10.1140/epje/i2016-16077-6}, + urldate = {2024-06-03}, + abstract = {Crystallization, a prototypical self-organization process during which a disordered state spontaneously transforms into a crystal characterized by a regular arrangement of its building blocks, usually proceeds by nucleation and growth. In the initial stages of the transformation, a localized nucleus of the new phase forms in the old one due to a random fluctuation. Most of these nuclei disappear after a short time, but rarely a crystalline embryo may reach a critical size after which further growth becomes thermodynamically favorable and the entire system is converted into the new phase. In this article, we will discuss several theoretical concepts and computational methods to study crystallization. More specifically, we will address the rare event problem arising in the simulation of nucleation processes and explain how to calculate nucleation rates accurately. Particular attention is directed towards discussing statistical tools to analyze crystallization trajectories and identify the transition mechanism.}, + langid = {english}, + keywords = {Flowing Matter: Liquids and Complex Fluids} } diff --git a/examples/toy_models/plot_1D_Double_Well.py b/examples/toy_models/plot_1D_Double_Well.py index aaa2fff..5be109f 100644 --- a/examples/toy_models/plot_1D_Double_Well.py +++ b/examples/toy_models/plot_1D_Double_Well.py @@ -18,17 +18,17 @@ # Plot of Free Energy and Force x_values = np.linspace(-7, 7, 100) -fig, axs = plt.subplots(1, 2) -axs[0].plot(x_values, free_energy(x_values)) -axs[1].plot(x_values, force_function(x_values.reshape(len(x_values), 1))) -axs[0].set_title("Potential") -axs[0].set_xlabel("$x$") -axs[0].set_ylabel("$V(x)$") -axs[0].grid() -axs[1].set_title("Force") -axs[1].set_xlabel("$x$") -axs[1].set_ylabel("$F(x)$") -axs[1].grid() +# fig, axs = plt.subplots(1, 2) +# axs[0].plot(x_values, free_energy(x_values)) +# axs[1].plot(x_values, force_function(x_values.reshape(len(x_values), 1))) +# axs[0].set_title("Potential") +# axs[0].set_xlabel("$x$") +# axs[0].set_ylabel("$V(x)$") +# axs[0].grid() +# axs[1].set_title("Force") +# axs[1].set_xlabel("$x$") +# axs[1].set_ylabel("$F(x)$") +# axs[1].grid() # Define model to simulate and type of simulator to use dt = 1e-3 @@ -45,12 +45,12 @@ time_steps = 10000 data = simulator.run(time_steps, q0, save_every=1) -# Plot resulting Trajectories -fig, axs = plt.subplots() -for n, trj in enumerate(data): - axs.plot(trj["x"]) - axs.set_title("Trajectory") - +# # Plot resulting Trajectories +# fig, axs = plt.subplots() +# for n, trj in enumerate(data): +# axs.plot(trj["x"]) +# axs.set_title("Trajectory") +# fig, axs = plt.subplots(1, 2) axs[0].set_title("Force") @@ -68,23 +68,27 @@ axs[1].plot(xfa, model_simu.diffusion(xfa.reshape(-1, 1)), label="Exact") trainforce = fl.functions.Polynomial(deg=3, coefficients=np.array([0, 0, 0, 0])) trainmodel = fl.models.Overdamped(force=trainforce, diffusion=fl.functions.Polynomial(deg=0, coefficients=np.asarray([0.9])), has_bias=False) -for name,marker, transitioncls in zip( +for name, marker, transitioncls in zip( ["Euler", "Ozaki", "ShojiOzaki", "Elerian", "Kessler", "Drozdov"], - ["x", "|",".","1","2","3"], + ["x", "|", ".", "1", "2", "3"], [ fl.EulerDensity, - fl.OzakiDensity, - fl.ShojiOzakiDensity, - fl.ElerianDensity, - fl.KesslerDensity, - fl.DrozdovDensity, ], ): estimator = fl.LikelihoodEstimator(transitioncls(trainmodel)) res = estimator.fit_fetch(data) - print(res.coefficients) - axs[0].plot(xfa, res.force(xfa.reshape(-1, 1)),marker=marker, label=name) - axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), marker=marker,label=name) + # print(res.coefficients) + axs[0].plot(xfa, res.force(xfa.reshape(-1, 1)), marker=marker, label=name) + axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), marker=marker, label=name) axs[0].legend() axs[1].legend() + +plt.figure() + +x_mfpt, mfpt = fl.analysis.mfpt_1d(model_simu, -5.0, [-10.0, 10.0], Npoints=500) +plt.plot(x_mfpt, mfpt) +x_mfpt, mfpt = fl.analysis.mfpt_1d(model_simu, 5.0, [-10.0, 10.0], Npoints=500) +plt.plot(x_mfpt, mfpt) +print(fl.analysis.mfpt_1d(model_simu, -5.0, [-10.0, 10.0], x_start=5.0, Npoints=500)) +print(fl.analysis.mfpt_1d(model_simu, 5.0, [-10.0, 10.0], x_start=-5.0, Npoints=500)) plt.show() diff --git a/folie/analysis/overdamped_1d.py b/folie/analysis/overdamped_1d.py index ebae799..167a5a2 100644 --- a/folie/analysis/overdamped_1d.py +++ b/folie/analysis/overdamped_1d.py @@ -1,4 +1,4 @@ -""" +""" Set of analysis methods focused on 1D overdamped models """ @@ -7,33 +7,89 @@ def free_energy_profile_1d(model, x): - """ - From the force and diffusion construct the free energy profile + r""" + From the force F(x) and diffusion D(x) construct the free energy profile V(x) using the formula + + .. math:: + F(x) = -D(x) \nabla V(x) + \mathrm{div} D(x) + """ x = x.ravel() diff_prime_val = model.diffusion.grad_x(x.reshape(-1, 1)).ravel() force_val = model.force(x.reshape(-1, 1)).ravel() diff_val = model.diffusion(x.reshape(-1, 1)).ravel() - diff_U = -1 * (force_val - 2 * diff_val * diff_prime_val) / diff_val ** 2 + diff_U = diff_U = -force_val / diff_val + diff_prime_val pmf = cumulative_trapezoid(diff_U, x, initial=0.0) return pmf - np.min(pmf) -def mfpt_1d(model, x_start: float, x_abs: float, Npoints=500, cumulative=False): - """ - Compute Mean first passage time from x0 to x1 +def mfpt_1d(model, x_end: float, x_range, Npoints=500, x_start=None): + r""" + Compute the mean first passage time from any point x within x_range to x_end, or from x_start to x_end if x_start is defined. + + It use numerical integration of the following formula for point from x_range[0] to x_end :cite:p:`Jungblut2016` + + .. math:: + MFPT(x,x_{end}) = \int_x^{x_{end}} \mathrm{d}y \frac{e^{\beta V(y)}}{D(y)} \int_{x\_range[0]}^y \mathrm{d} z e^{-\beta V(y)} + + and for point from x_end to x_range[1] + + .. math:: + + MFPT(x,x_{end}) = \int^x_{x_{end}} \mathrm{d}y \frac{e^{\beta V(y)}}{D(y)} \int^{x\_range[1]}_y \mathrm{d} z e^{-\beta V(y)} + + Parameters + ------------ + + model: + A fitted overdamped model + + x_end: float + The point to reach + + x_start: float, default to None + If this is not None it return the MFPT from x_start to x_end, otherwise it return the mean first passage from any point within x_range to x_end + + x_range: + A range of integration, It should be big enough to be able to compute the normalisation factor of the steady state probability. + + Npoints: int, default=500 + Number of point to use for + + References + -------------- + + .. footbibliography:: + """ - min_x = np.min(self.knots_diff) - x_range = np.linspace(min_x, x_abs, Npoints) + if x_start is not None: + if x_start < x_end: + int_range = np.linspace(x_range[0], x_end, Npoints) + prob_well = cumulative_trapezoid(np.exp(-free_energy_profile_1d(model, int_range)), int_range, initial=0.0) + else: + int_range = np.linspace(x_end, x_range[1], Npoints) + prob_well = cumulative_trapezoid(np.exp(-free_energy_profile_1d(model, int_range)), int_range, initial=0.0) + prob_well -= prob_well[-1] + + x_int = np.linspace(x_start, x_end, Npoints) + return np.trapz(np.exp(free_energy_profile_1d(model, x_int)) * np.interp(x_int, int_range, prob_well) / model.diffusion(x_int.reshape(-1, 1)).ravel(), x_int) - prob_well = cumulative_trapezoid(np.exp(-self.free_energy_profile(x_range), x_range), initial=0.0) - x_int = np.linspace(x_start, x_abs, Npoints) - if cumulative: - res = cumulative_trapezoid(np.exp(self.free_energy_profile(x_int)) * np.interp(x_int, x_range, prob_well) / self.diffusion(x_int, 0.0) ** 2, x_int, initial=0.0) - return res else: - return np.trapz(np.exp(self.free_energy_profile(x_int)) * np.interp(x_int, x_range, prob_well) / self.diffusion(x_int, 0.0) ** 2, x_int) + # Compute lower part + int_range = np.linspace(x_range[0], x_end, Npoints) + prob_well = cumulative_trapezoid(np.exp(-free_energy_profile_1d(model, int_range)), int_range, initial=0.0) + + x_int_lower = np.linspace(x_end, x_range[0], Npoints) + res_lower = -1 * cumulative_trapezoid(np.exp(free_energy_profile_1d(model, x_int_lower)) * np.interp(x_int_lower, int_range, prob_well) / model.diffusion(x_int_lower.reshape(-1, 1)).ravel(), x_int_lower, initial=0.0) + + int_range = np.linspace(x_end, x_range[1], Npoints) + prob_well = -1 * cumulative_trapezoid(np.exp(-free_energy_profile_1d(model, int_range)), int_range, initial=0.0) + prob_well -= prob_well[-1] + # return (int_range, prob_well) + x_int_upper = np.linspace(x_end, x_range[1], Npoints) + res_upper = cumulative_trapezoid(np.exp(free_energy_profile_1d(model, x_int_upper)) * np.interp(x_int_upper, int_range, prob_well) / model.diffusion(x_int_upper.reshape(-1, 1)).ravel(), x_int_upper, initial=0.0) + return np.hstack((x_int_lower[::-1], x_int_upper)), np.hstack((res_lower[::-1], res_upper)) diff --git a/folie/functions/base.py b/folie/functions/base.py index b9b6d42..6e0a2c0 100644 --- a/folie/functions/base.py +++ b/folie/functions/base.py @@ -225,6 +225,7 @@ class ModelOverlay(Function): output_shape: tuple or array The output shape of the term + """ def __init__(self, model, function_name, output_shape=None, **kwargs):