Improve doc and modify formulas for analysis

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
+
+

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

docs/index.rst

@@ -83,6 +83,7 @@ Table of contents
    api/models
    api/functions
    api/estimation
+   api/analysis
 
 Indices and tables
 ==================

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,

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}
 }

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()

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))

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):