Remove OzakiDensity and ShojiOzakiDensity

benchmark/benchmark_likelihood.py

@@ -19,7 +19,7 @@ def data(request):
 
 
 @pytest.mark.parametrize("data", ["numpy", "dask"], indirect=True)
-@pytest.mark.parametrize("transitioncls", [fl.EulerDensity, fl.OzakiDensity, fl.ShojiOzakiDensity, fl.ElerianDensity, fl.KesslerDensity, fl.DrozdovDensity])
+@pytest.mark.parametrize("transitioncls", [fl.EulerDensity, fl.ElerianDensity, fl.KesslerDensity, fl.DrozdovDensity])
 def test_likelihood_functions(data, request, benchmark, transitioncls):
     fun = fl.functions.Linear()
     model = fl.models.Overdamped(fun, fun.copy())

docs/api/estimation.rst

@@ -36,10 +36,6 @@ Overdamped Transition density
 
    EulerDensity
 
-   OzakiDensity
-
-   ShojiOzakiDensity
-
    ElerianDensity
 
    KesslerDensity

docs/for_developper.rst

@@ -25,7 +25,7 @@ Likelihood estimation: using Transition density
 
 The list of implemented transition density and their relations
 
-.. inheritance-diagram:: folie.EulerDensity folie.OzakiDensity folie.ShojiOzakiDensity folie.ElerianDensity folie.KesslerDensity folie.DrozdovDensity
+.. inheritance-diagram:: folie.EulerDensity folie.ElerianDensity folie.KesslerDensity folie.DrozdovDensity
     :top-classes: folie.estimation.transitionDensity.TransitionDensity
     :parts: 1
 

examples/plot_biasedOU.py

@@ -51,12 +51,10 @@
 axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), "--", label=name)
 
 for name, marker, transitioncls in zip(
-    ["Euler", "Ozaki", "ShojiOzaki", "Elerian", "Kessler", "Drozdov"],
-    ["+", "x", "P", "1", "2", "3"],
+    ["Euler", "Elerian", "Kessler", "Drozdov"],
+    ["+", "1", "2", "3"],
     [
         fl.EulerDensity,
-        fl.OzakiDensity,
-        fl.ShojiOzakiDensity,
         fl.ElerianDensity,
         fl.KesslerDensity,
         fl.DrozdovDensity,

examples/plot_likelihood.py

@@ -39,11 +39,9 @@
 
 
 for name, transitioncls in zip(
-    ["Euler", "Ozaki", "ShojiOzaki", "Elerian", "Kessler", "Drozdov"],
+    ["Euler", "Elerian", "Kessler", "Drozdov"],
     [
         fl.EulerDensity,
-        fl.OzakiDensity,
-        fl.ShojiOzakiDensity,
         fl.ElerianDensity,
         fl.KesslerDensity,
         fl.DrozdovDensity,

examples/statistical_performances/overdampedOU.ipynb

@@ -75,12 +75,10 @@
     "fig, axs = plt.subplots(1, len(model_simu.coefficients))\n",
     "\n",
     "for name, transitioncls in zip(\n",
-    "    [\"Exact\", \"Euler\", \"Ozaki\", \"ShojiOzaki\", \"Elerian\", \"Kessler\", \"Drozdov\"],\n",
+    "    [\"Exact\", \"Euler\", \"Elerian\", \"Kessler\", \"Drozdov\"],\n",
     "    [\n",
     "        fl.ExactDensity,\n",
     "        fl.EulerDensity,\n",
-    "        fl.OzakiDensity,\n",
-    "        fl.ShojiOzakiDensity,\n",
     "        fl.ElerianDensity,\n",
     "        fl.KesslerDensity,\n",
     "        fl.DrozdovDensity,\n",

examples/toy_models/plot_1D_Double_Well.py

@@ -13,7 +13,7 @@
 coeff = 0.2 * np.array([0, 0, -4.5, 0, 0.1])
 free_energy = np.polynomial.Polynomial(coeff)
 D = 0.5
-force_coeff = D*np.array([-coeff[1], -2 * coeff[2], -3 * coeff[3], -4 * coeff[4]])
+force_coeff = D * np.array([-coeff[1], -2 * coeff[2], -3 * coeff[3], -4 * coeff[4]])
 
 force_function = fl.functions.Polynomial(deg=3, coefficients=force_coeff)
 diff_function = fl.functions.Polynomial(deg=0, coefficients=np.array(D))
@@ -72,12 +72,10 @@
 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(
-    ["Euler", "Ozaki", "ShojiOzaki", "Elerian", "Kessler", "Drozdov"],
-    ["x", "|", ".", "1", "2", "3"],
+    ["Euler", "Elerian", "Kessler", "Drozdov"],
+    ["x", "1", "2", "3"],
     [
         fl.EulerDensity,
-        fl.OzakiDensity,
-        fl.ShojiOzakiDensity,
         fl.ElerianDensity,
         fl.KesslerDensity,
         fl.DrozdovDensity,

examples/toy_models/plot_biased_1D_Double_Well.py

@@ -13,9 +13,9 @@
 
 coeff = 0.2 * np.array([0, 0, -4.5, 0, 0.1])
 free_energy = np.polynomial.Polynomial(coeff)
-D= np.array([0.5])
+D = np.array([0.5])
 
-force_coeff = D*np.array([-coeff[1], -2 * coeff[2], -3 * coeff[3], -4 * coeff[4]])
+force_coeff = D * np.array([-coeff[1], -2 * coeff[2], -3 * coeff[3], -4 * coeff[4]])
 force_function = fl.functions.Polynomial(deg=3, coefficients=force_coeff)
 diff_function = fl.functions.Polynomial(deg=0, coefficients=D)
 
@@ -73,30 +73,27 @@
 axs[1].plot(xfa, model_simu.diffusion(xfa.reshape(-1, 1)), label="Exact")
 
 domain = fl.MeshedDomain.create_from_range(np.linspace(data.stats.min, data.stats.max, 4).ravel())
-trainmodel = fl.models.Overdamped(force = fl.functions.BSplinesFunction(domain),has_bias=True)
+trainmodel = fl.models.Overdamped(force=fl.functions.BSplinesFunction(domain), has_bias=True)
 
-for name,marker, transitioncls in zip(
-    ["Euler", "Ozaki", "ShojiOzaki", "Elerian", "Kessler", "Drozdov"],
-    ["x", "|",".","1","2","3"],
+for name, marker, transitioncls in zip(
+    ["Euler", "Elerian", "Kessler", "Drozdov"],
+    ["x", "1", "2", "3"],
     [
         fl.EulerDensity,
-        fl.OzakiDensity,
-        fl.ShojiOzakiDensity,
         fl.ElerianDensity,
         fl.KesslerDensity,
         fl.DrozdovDensity,
     ],
 ):
-    trainmodel = fl.models.Overdamped(force = fl.functions.BSplinesFunction(domain),has_bias=True)
-    estimator = fl.LikelihoodEstimator(transitioncls(trainmodel),n_jobs=4) 
-
+    trainmodel = fl.models.Overdamped(force=fl.functions.BSplinesFunction(domain), has_bias=True)
+    estimator = fl.LikelihoodEstimator(transitioncls(trainmodel), n_jobs=4)
 
     res = estimator.fit_fetch(deepcopy(data))
 
     print(name, res.coefficients)
     res.remove_bias()
-    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].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.show()

folie/estimation/overdamped_transitionDensity.py

@@ -165,67 +165,6 @@ def e_step(self, weight, trj, coefficients, mu0, sig0):
         return muh[0, self._model.dim_h :] / weight, Sigh[0, self._model.dim_h :, self._model.dim_h :] / weight  # Return µ0 and sig0
 
 
-class OzakiDensity(TransitionDensity):
-    def __init__(self, model):
-        """
-        Class which represents the Ozaki approximation transition density for a model
-        :param model: the SDE model, referenced during calls to the transition density
-        """
-        super().__init__(model)
-
-    def _logdensity1D(self, x, xt, dt: float, bias=0.0, **kwargs):
-        """
-        The transition density obtained via Ozaki expansion
-        :param x: float or array, the current value
-        :param xt: float or array, the value to transition to  (must be same dimension as x)
-        :param dt: float, the time step between x and xt
-        :return: probability (same dimension as x and xt)
-        """
-        sig = 2 * self._model.diffusion(x, **kwargs).ravel()
-        mu = self._model.meandispl(x, bias, **kwargs).ravel()
-        mu_x = self._model.meandispl.grad_x(x, bias, **kwargs).ravel()
-        temp = mu * (np.exp(mu_x * dt) - 1) / mu_x
-
-        Mt = x.ravel() + temp
-        Kt = (2 / dt) * np.log(1 + temp / x.ravel())
-        Vt = np.sqrt(sig * (np.exp(Kt * dt) - 1) / Kt)
-
-        return gaussian_likelihood_1D(xt, Mt, Vt)
-
-
-class ShojiOzakiDensity(TransitionDensity):
-    def __init__(self, model):
-        """
-        Class which represents the Shoji-Ozaki approximation transition density for a model
-        :param model: the SDE model, referenced during calls to the transition density
-        """
-        super().__init__(model)
-
-    def _logdensity1D(self, x, xt, dt: float, bias=0.0, **kwargs):
-        """
-        The transition density obtained via Shoji-Ozaki expansion
-        :param x: float or array, the current value
-        :param xt: float or array, the value to transition to  (must be same dimension as x)
-        :param dt: float, the time step between x and xt
-        :return: probability (same dimension as x and xt)
-        """
-        sig = np.sqrt(2 * self._model.diffusion(x, **kwargs).ravel())
-        mu = self._model.meandispl(x, bias, **kwargs).ravel()
-
-        Mt = 0.5 * sig**2 * self._model.meandispl.hessian_x(x, bias, **kwargs).ravel()  # + self._model.meandispl_t(x)  #Time homogenous model
-        Lt = self._model.meandispl.grad_x(x, bias, **kwargs).ravel()
-        if (Lt == 0).any():  # TODO: need to fix this
-            B = sig * np.sqrt(dt)
-            A = x.ravel() + mu * dt + Mt * dt**2 / 2
-        else:
-            B = sig * np.sqrt((np.exp(2 * Lt * dt) - 1) / (2 * Lt))
-
-            elt = np.exp(Lt * dt) - 1
-            A = x.ravel() + mu / Lt * elt + Mt / (Lt**2) * (elt - Lt * dt)
-
-        return gaussian_likelihood_1D(xt, A, B)
-
-
 class ElerianDensity(EulerDensity):
     use_jac = False
 

tests/test_likelihood.py

@@ -22,7 +22,7 @@ def data(request):
 
 
 @pytest.mark.parametrize("data", ["numpy", "dask"], indirect=True)
-@pytest.mark.parametrize("transitioncls", [fl.OzakiDensity, fl.ShojiOzakiDensity, fl.ElerianDensity, fl.KesslerDensity, fl.DrozdovDensity])
+@pytest.mark.parametrize("transitioncls", [fl.ElerianDensity, fl.KesslerDensity, fl.DrozdovDensity])
 def test_likelihood(data, request, transitioncls):
     print(data.representative_array())
     fun_lin = fl.functions.Linear().fit(data.representative_array(), np.ones(data.representative_array().shape[0]))