Changing meandispl and force into drift

docs/notebooks/estimation.ipynb

@@ -96,7 +96,7 @@
     "ax.set_xlabel(\"$x$\")\n",
     "ax.set_ylabel(\"$F(x)$\")\n",
     "ax.grid()\n",
-    "ax.plot(xfa, model.force(xfa.reshape(-1, 1)))\n"
+    "ax.plot(xfa, model.drift(xfa.reshape(-1, 1)))\n"
    ]
   },
   {

examples/example_em.py

@@ -33,11 +33,11 @@
 
 fig, axs = plt.subplots(1, 3)
 # Force plot
-axs[0].set_title("Force")
+axs[0].set_title("Drift")
 axs[0].set_xlabel("$x$")
 axs[0].set_ylabel("$F(x)$")
 axs[0].grid()
-axs[0].plot(xfa, model.force(xfa.reshape(-1, 1)))
+axs[0].plot(xfa, model.pos_drift(xfa.reshape(-1, 1)))
 
 
 # Friction plot

examples/plot_biasedOU.py

@@ -25,8 +25,13 @@
     axs[0].plot(trj["x"])
     axs[1].plot(xmax[:, n])
 
+axs[0].set_title("Trajectories")
+axs[0].set_xlabel("step")
+axs[1].set_title("Bias")
+axs[1].set_xlabel("step")
+
 fig, axs = plt.subplots(1, 2)
-axs[0].set_title("Force")
+axs[0].set_title("Drift")
 axs[0].set_xlabel("$x$")
 axs[0].set_ylabel("$F(x)$")
 axs[0].grid()
@@ -39,15 +44,15 @@
 
 xfa = np.linspace(-7.0, 7.0, 75)
 model_simu.remove_bias()
-axs[0].plot(xfa, model_simu.force(xfa.reshape(-1, 1)), label="Exact")
+axs[0].plot(xfa, model_simu.drift(xfa.reshape(-1, 1)), label="Exact")
 axs[1].plot(xfa, model_simu.diffusion(xfa.reshape(-1, 1)), label="Exact")
 
 name = "KramersMoyal"
 estimator = fl.KramersMoyalEstimator(fl.models.OrnsteinUhlenbeck(has_bias=True))
 res = estimator.fit_fetch(data)
 print(name, res.coefficients, res.is_biased)
 res.remove_bias()
-axs[0].plot(xfa, res.force(xfa.reshape(-1, 1)), "--", label=name)
+axs[0].plot(xfa, res.drift(xfa.reshape(-1, 1)), "--", label=name)
 axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), "--", label=name)
 
 for name, marker, transitioncls in zip(
@@ -64,7 +69,7 @@
     res = estimator.fit_fetch(data)
     print(name, res.coefficients, res.is_biased)
     res.remove_bias()
-    axs[0].plot(xfa, res.force(xfa.reshape(-1, 1)), marker, label=name)
+    axs[0].plot(xfa, res.drift(xfa.reshape(-1, 1)), marker, label=name)
     axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), marker, label=name)
 axs[0].legend()
 axs[1].legend()

examples/plot_example.py

@@ -32,11 +32,11 @@
 
 fig, axs = plt.subplots(1, 2)
 # Force plot
-axs[0].set_title("Force")
+axs[0].set_title("Drift")
 axs[0].set_xlabel("$x$")
 axs[0].set_ylabel("$F(x)$")
 axs[0].grid()
-axs[0].plot(xfa, model.force(xfa.reshape(-1, 1)))
+axs[0].plot(xfa, model.drift(xfa.reshape(-1, 1)))
 
 # Diffusion plot
 axs[1].set_title("Diffusion")

examples/plot_fem.py

@@ -23,7 +23,7 @@
 
 fig, axs = plt.subplots(1, 2)
 # Force plot
-axs[0].set_title("Force")
+axs[0].set_title("Drift")
 axs[0].set_xlabel("$x$")
 axs[0].set_ylabel("$F(x)$")
 axs[0].grid()
@@ -49,7 +49,7 @@
 
     model = estimator.fit_fetch(data)
 
-    axs[0].plot(xfa, model.force(xfa.reshape(-1, 1)))
+    axs[0].plot(xfa, model.drift(xfa.reshape(-1, 1)))
     axs[1].plot(xfa, model.diffusion(xfa.reshape(-1, 1)))
 
 plt.show()

examples/plot_likelihood.py

@@ -22,7 +22,7 @@
 model = fl.models.BrownianMotion()
 
 fig, axs = plt.subplots(1, 2)
-axs[0].set_title("Force")
+axs[0].set_title("Drift")
 axs[0].set_xlabel("$f$")
 axs[0].set_ylabel("$L(f,1.0)$")
 axs[0].grid()
@@ -34,7 +34,7 @@
 axs[1].set_ylabel("$L(1.0,\\sigma)$")
 
 
-force_range = np.linspace(-1, 1, 25)
+drift_range = np.linspace(-1, 1, 25)
 diff_range = np.linspace(0.5, 2, 25)
 
 
@@ -49,10 +49,10 @@
 ):
     likelihood = transitioncls(model)
     likelihood.preprocess_traj(data[0])
-    likelihood_vals_force = np.zeros_like(force_range)
-    for n, f in enumerate(force_range):
-        likelihood_vals_force[n] = likelihood(1.0, data[0], np.array([f, 1.0]))[0]
-    axs[0].plot(force_range, likelihood_vals_force, label=name)
+    likelihood_vals_drift = np.zeros_like(drift_range)
+    for n, f in enumerate(drift_range):
+        likelihood_vals_drift[n] = likelihood(1.0, data[0], np.array([f, 1.0]))[0]
+    axs[0].plot(drift_range, likelihood_vals_drift, label=name)
     likelihood_vals_diff = np.zeros_like(diff_range)
     for n, d in enumerate(diff_range):
         likelihood_vals_diff[n] = likelihood(1.0, data[0], np.array([1.0, d]))[0]

examples/toy_models/plot_1D_Double_Well.py

@@ -13,16 +13,16 @@
 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]])
+drift_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)
+drift_function = fl.functions.Polynomial(deg=3, coefficients=drift_coeff)
 diff_function = fl.functions.Polynomial(deg=0, coefficients=np.array(D))
 
 # 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[1].plot(x_values, drift_function(x_values.reshape(len(x_values), 1)))
 # axs[0].set_title("Potential")
 # axs[0].set_xlabel("$x$")
 # axs[0].set_ylabel("$V(x)$")
@@ -34,7 +34,7 @@
 
 # Define model to simulate and type of simulator to use
 dt = 1e-3
-model_simu = fl.models.overdamped.Overdamped(force_function, diffusion=diff_function)
+model_simu = fl.models.overdamped.Overdamped(drift_function, diffusion=diff_function)
 simulator = fl.simulations.Simulator(fl.simulations.EulerStepper(model_simu), dt)
 
 
@@ -66,10 +66,17 @@
 axs[1].grid()
 
 xfa = np.linspace(-7.0, 7.0, 75)
-axs[0].plot(xfa, model_simu.force(xfa.reshape(-1, 1)), label="Exact")
+axs[0].plot(xfa, model_simu.pos_drift(xfa.reshape(-1, 1)), label="Exact")
 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)
+traindrift = fl.functions.Polynomial(deg=3, coefficients=np.array([0, 0, 0, 0]))
+trainmodel = fl.models.Overdamped(traindrift, diffusion=fl.functions.Polynomial(deg=0, coefficients=np.asarray([0.9])), has_bias=False)
+
+name = "KramersMoyal"
+estimator = fl.KramersMoyalEstimator(trainmodel)
+res = estimator.fit_fetch(data)
+axs[0].plot(xfa, res.drift(xfa.reshape(-1, 1)), "--", label=name)
+axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), "--", label=name)
+
 
 for name, marker, transitioncls in zip(
     ["Euler", "Elerian", "Kessler", "Drozdov"],
@@ -84,13 +91,13 @@
     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[0].plot(xfa, res.pos_drift(xfa.reshape(-1, 1)), marker=marker, label=name)
     axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), marker=marker, label=name)
 
 
 name = "Kramers Moyal"
 res = fl.KramersMoyalEstimator(trainmodel).fit_fetch(data)
-axs[0].plot(xfa, res.force(xfa.reshape(-1, 1)), "--", label=name)
+axs[0].plot(xfa, res.pos_drift(xfa.reshape(-1, 1)), "--", label=name)
 axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), "--", label=name)
 
 

examples/toy_models/plot_2D_Double_Well.py

@@ -17,11 +17,11 @@
 y = np.linspace(-1.8, 1.8, 36)
 input = np.transpose(np.array([x, y]))
 
-D=0.5
-diff_function= fl.functions.Polynomial(deg=0,coefficients=D * np.eye(2,2))
-a,b = 5, 10
-drift_quartic2d= fl.functions.Quartic2D(a=D*a,b=D*b)  # simple way to multiply D*Potential here force is the SDE force (meandispl)  ## use this when you need the drift ###
-quartic2d= fl.functions.Quartic2D(a=a,b=b)            # Real potential , here force is just -grad pot ## use this when you need the potential energy ###
+D = 0.5
+diff_function = fl.functions.Polynomial(deg=0, coefficients=D * np.eye(2, 2))
+a, b = 5, 10
+drift_quartic2d = fl.functions.Quartic2D(a=D * a, b=D * b)  # simple way to multiply D*Potential here force is the SDE force (meandispl)  ## use this when you need the drift ###
+quartic2d = fl.functions.Quartic2D(a=a, b=b)  # Real potential , here force is just -grad pot ## use this when you need the potential energy ###
 
 X, Y = np.meshgrid(x, y)
 
@@ -39,8 +39,8 @@
 ax.quiver(x, y, U, V)
 ax.set_title("Force")
 
-dt= 5e-4
-model_simu = fl.models.overdamped.Overdamped(force=drift_quartic2d, diffusion=diff_function)
+dt = 5e-4
+model_simu = fl.models.overdamped.Overdamped(drift_quartic2d, diffusion=diff_function)
 simulator = fl.simulations.Simulator(fl.simulations.EulerStepper(model_simu), dt)
 
 # initialize positions

examples/toy_models/plot_biased_1D_Double_Well.py

@@ -15,26 +15,26 @@
 free_energy = np.polynomial.Polynomial(coeff)
 D = np.array([0.5])
 
-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)
+drift_coeff = D * np.array([-coeff[1], -2 * coeff[2], -3 * coeff[3], -4 * coeff[4]])
+drift_function = fl.functions.Polynomial(deg=3, coefficients=drift_coeff)
 diff_function = fl.functions.Polynomial(deg=0, coefficients=D)
 
 # 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[1].plot(x_values, drift_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_title("Drift")
 axs[1].set_xlabel("$x$")
 axs[1].set_ylabel("$F(x)$")
 axs[1].grid()
 
 # Define model to simulate and type of simulator to use
-model_simu = fl.models.overdamped.Overdamped(force_function, diffusion=diff_function)
+model_simu = fl.models.overdamped.Overdamped(drift_function, diffusion=diff_function)
 simulator = fl.simulations.ABMD_Simulator(fl.simulations.EulerStepper(model_simu), 1e-3, k=10.0, xstop=6.0)
 
 # initialize positions
@@ -57,7 +57,7 @@
     axs[1].grid()
 
 fig, axs = plt.subplots(1, 2)
-axs[0].set_title("Force")
+axs[0].set_title("Drift")
 axs[0].set_xlabel("$x$")
 axs[0].set_ylabel("$F(x)$")
 axs[0].grid()
@@ -69,11 +69,19 @@
 
 xfa = np.linspace(-7.0, 7.0, 75)
 model_simu.remove_bias()
-axs[0].plot(xfa, model_simu.force(xfa.reshape(-1, 1)), label="Exact")
+axs[0].plot(xfa, model_simu.drift(xfa.reshape(-1, 1)), label="Exact")
 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(fl.functions.BSplinesFunction(domain), has_bias=True)
+
+name = "KramersMoyal"
+estimator = fl.KramersMoyalEstimator(trainmodel)
+res = estimator.fit_fetch(data)
+res.remove_bias()
+axs[0].plot(xfa, res.drift(xfa.reshape(-1, 1)), "--", label=name)
+axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), "--", label=name)
+
 
 for name, marker, transitioncls in zip(
     ["Euler", "Elerian", "Kessler", "Drozdov"],
@@ -85,14 +93,14 @@
         fl.DrozdovDensity,
     ],
 ):
-    trainmodel = fl.models.Overdamped(force=fl.functions.BSplinesFunction(domain), has_bias=True)
+    trainmodel = fl.models.Overdamped(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[0].plot(xfa, res.drift(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()

examples/toy_models/plot_biased_2D_Double_Well.py

@@ -18,11 +18,11 @@
 y = np.linspace(-1.8, 1.8, 36)
 input = np.transpose(np.array([x, y]))
 
-D=0.5
-diff_function= fl.functions.Polynomial(deg=0,coefficients=D * np.eye(2,2))
-a,b = 5, 10
-drift_quartic2d= fl.functions.Quartic2D(a=D*a,b=D*b)  # simple way to multiply D*Potential here force is the SDE force (meandispl)  ## use this when you need the drift ###
-quartic2d= fl.functions.Quartic2D(a=a,b=b)            # Real potential , here force is just -grad pot ## use this when you need the potential energy ###
+D = 0.5
+diff_function = fl.functions.Polynomial(deg=0, coefficients=D * np.eye(2, 2))
+a, b = 5, 10
+drift_quartic2d = fl.functions.Quartic2D(a=D * a, b=D * b)  # simple way to multiply D*Potential here force is the SDE force (meandispl)  ## use this when you need the drift ###
+quartic2d = fl.functions.Quartic2D(a=a, b=b)  # Real potential , here force is just -grad pot ## use this when you need the potential energy ###
 
 X, Y = np.meshgrid(x, y)
 
@@ -47,7 +47,7 @@ def colvar(x, y):
 
 
 dt = 1e-3
-model_simu = fl.models.overdamped.Overdamped(force=drift_quartic2d, diffusion=diff_function)
+model_simu = fl.models.overdamped.Overdamped(drift_quartic2d, diffusion=diff_function)
 simulator = fl.simulations.ABMD_2D_to_1DColvar_Simulator(fl.simulations.EulerStepper(model_simu), dt, colvar=colvar, k=25.0, qstop=1.2)
 
 # Choose number of trajectories and initialize positions
@@ -139,10 +139,10 @@ def colvar(x, y):
 trainmodel = fl.models.OverdampedSplines1D(domain=domain)
 
 xfa = np.linspace(proj_data.stats.min, proj_data.stats.max, 75)
-force_exact = (xfa**2 - 1.0) ** 2
+drift_exact = (xfa**2 - 1.0) ** 2
 
 fig, axs = plt.subplots(1, 2)
-axs[0].set_title("Force")
+axs[0].set_title("Drift")
 axs[0].set_xlabel("$x$")
 axs[0].set_ylabel("$F(x)$")
 axs[0].grid()
@@ -155,12 +155,12 @@ def colvar(x, y):
 KM_Estimator = fl.KramersMoyalEstimator(deepcopy(trainmodel))
 res_KM = KM_Estimator.fit_fetch(proj_data)
 
-axs[0].plot(xfa, res_KM.force(xfa.reshape(-1, 1)),  marker="x",label="KramersMoyal")
-axs[1].plot(xfa, res_KM.diffusion(xfa.reshape(-1, 1)), marker="x",label="KramersMoyal")
+axs[0].plot(xfa, res_KM.drift(xfa.reshape(-1, 1)), marker="x", label="KramersMoyal")
+axs[1].plot(xfa, res_KM.diffusion(xfa.reshape(-1, 1)), marker="x", label="KramersMoyal")
 print("KramersMoyal ", res_KM.coefficients)
-for name,marker, transitioncls in zip(
+for name, marker, transitioncls in zip(
     ["Euler", "Elerian", "Kessler", "Drozdov"],
-        ["|","1","2","3"],
+    ["|", "1", "2", "3"],
     [
         fl.EulerDensity,
         fl.ElerianDensity,
@@ -171,8 +171,8 @@ def colvar(x, y):
     estimator = fl.LikelihoodEstimator(transitioncls(deepcopy(trainmodel)))
     res = estimator.fit_fetch(deepcopy(proj_data))
     print(name, 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].plot(xfa, res.drift(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()