Merge pull request #9 from Dani-dB/toy-models

toy-models branch pull request

examples/plot_biasedOU.py

@@ -14,7 +14,7 @@
 
 model_simu = fl.models.OrnsteinUhlenbeck(0.0, 1.2, 2.0)
 simulator = fl.simulations.ABMD_Simulator(fl.simulations.EulerStepper(model_simu), 1e-3, k=10.0, xstop=6.0)
-data = simulator.run(5000, np.zeros((25,)), 25)
+data = simulator.run(5000, np.zeros((25,)), 1)
 xmax = np.concatenate(simulator.xmax_hist, axis=1).T
 
 # Plot the resulting trajectories
@@ -30,7 +30,7 @@
 axs[0].set_ylabel("$F(x)$")
 axs[0].grid()
 
-axs[1].set_title("Force")
+axs[1].set_title("Diffusion")
 axs[1].set_xlabel("$x$")
 axs[1].set_ylabel("$D(x)$")
 axs[1].grid()

examples/toy_models/1D_simulations/plot_1D_Double_Well.py

@@ -0,0 +1,81 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import folie as fl
+
+coeff=0.1*np.array([0,0,-4.5,0,0.1])
+free_energy = np.polynomial.Polynomial(coeff)
+force_coeff=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.asarray([0.5]))
+
+# 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()
+
+# Define model to simulate and type of simulator to use
+dt=1e-3
+model_simu = fl.models.overdamped.Overdamped(force_function,diffusion=diff_function)
+simulator = fl.simulations.Simulator(fl.simulations.EulerStepper(model_simu), dt) 
+
+
+# initialize positions 
+ntraj=30
+q0= np.empty(ntraj)
+for i in range(len(q0)):
+    q0[i]=0
+# Calculate Trajectory
+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")
+
+
+fig, axs = plt.subplots(1, 2)
+axs[0].set_title("Force")
+axs[0].set_xlabel("$x$")
+axs[0].set_ylabel("$F(x)$")
+axs[0].grid()
+
+axs[1].set_title("Diffusion Coefficient")
+axs[1].set_xlabel("$x$")
+axs[1].set_ylabel("$D(x)$") 
+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[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, transitioncls in zip(
+   ["Euler", "Ozaki", "ShojiOzaki", "Elerian", "Kessler", "Drozdov"],
+   [
+       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)), label=name)
+   axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), label=name)
+axs[0].legend()
+axs[1].legend()
+plt.show() 

examples/toy_models/1D_simulations/plot_biased_1D_Double_Well.py

@@ -0,0 +1,84 @@
+
+import numpy as np
+import matplotlib.pyplot as plt
+import folie as fl
+
+coeff=0.1*np.array([0,0,-4.5,0,0.1])
+free_energy = np.polynomial.Polynomial(coeff)
+
+force_coeff=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.asarray([0.5]))
+
+# 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()
+
+# Define model to simulate and type of simulator to use
+model_simu = fl.models.overdamped.Overdamped(force_function,diffusion=diff_function)
+simulator = fl.simulations.ABMD_Simulator(fl.simulations.EulerStepper(model_simu), 1e-3, k=1.0, xstop=6.0)
+
+# initialize positions 
+ntraj=30
+q0= np.empty(ntraj)
+for i in range(len(q0)):
+    q0[i]=0
+# Calculate Trajectory
+time_steps=10000
+data = simulator.run(time_steps, q0, save_every=1)
+xmax = np.concatenate(simulator.xmax_hist, axis=1).T
+
+# Plot the resulting trajectories
+fig, axs = plt.subplots(1,2)
+for n, trj in enumerate(data):
+    axs[0].plot(trj["x"])
+    axs[1].plot(xmax[:, n])
+    axs[1].set_xlabel("$timestep$")
+    axs[1].set_ylabel("$x(t)$")
+    axs[1].grid()
+
+fig, axs = plt.subplots(1, 2)
+axs[0].set_title("Force")
+axs[0].set_xlabel("$x$")
+axs[0].set_ylabel("$F(x)$")
+axs[0].grid()
+
+axs[1].set_title("Diffusion Coefficient") # i think should be diffusion coefficient
+axs[1].set_xlabel("$x$")
+axs[1].set_ylabel("$D(x)$") 
+axs[1].grid()
+
+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[1].plot(xfa, model_simu.diffusion(xfa.reshape(-1, 1)), label="Exact")
+for name, transitioncls in zip(
+   ["Euler", "Ozaki", "ShojiOzaki", "Elerian", "Kessler", "Drozdov"],
+   [
+       fl.EulerDensity,
+       fl.OzakiDensity,
+       fl.ShojiOzakiDensity,
+       fl.ElerianDensity,
+       fl.KesslerDensity,
+       fl.DrozdovDensity,
+   ],
+):
+   estimator = fl.LikelihoodEstimator(transitioncls(fl.models.Overdamped(force_function,has_bias=True))) #diffusion= diff_function,
+   res = estimator.fit_fetch(data)
+   print(name, res.coefficients)
+   res.remove_bias()
+   axs[0].plot(xfa, res.force(xfa.reshape(-1, 1)), label=name)
+   axs[1].plot(xfa, res.diffusion(xfa.reshape(-1, 1)), label=name)
+axs[0].legend()
+axs[1].legend()
+plt.show() 

examples/toy_models/2D_simulations/plot_2D_Double_Well.py

@@ -0,0 +1,110 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import folie as fl
+from mpl_toolkits.mplot3d import Axes3D
+from copy import deepcopy
+
+""" Script for simulation of 2D double well and projection along user provided direction, No fitting is carried out   """
+x = np.linspace(-1.8,1.8,36)
+y = np.linspace(-1.8,1.8,36)
+input=np.transpose(np.array([x,y]))
+
+diff_function= fl.functions.Polynomial(deg=0,coefficients=np.asarray([0.5]) * np.eye(2,2))
+a,b = 5.0, 10.0
+quartic2d= fl.functions.Quartic2D(a=a,b=b)
+
+X,Y =np.meshgrid(x,y)
+print(X.shape,Y.shape)
+
+# Plot potential surface 
+pot = quartic2d.potential_plot(X,Y)
+print(pot.shape)
+fig = plt.figure()
+ax = plt.axes(projection='3d')
+ax.plot_surface(X,Y,pot, rstride=1, cstride=1,cmap='jet', edgecolor = 'none')
+
+# Plot Force function
+ff=quartic2d.force(input) # returns x and y components of the force : x_comp =ff[:,0] , y_comp =ff[:,1]
+U,V = np.meshgrid(ff[:,0],ff[:,1])
+fig, ax =plt.subplots()
+ax.quiver(x,y,U,V)
+ax.set_title('Force')
+
+model_simu=fl.models.overdamped.Overdamped(force=quartic2d,diffusion=diff_function)
+simulator=fl.simulations.Simulator(fl.simulations.EulerStepper(model_simu), 1e-3)
+
+# initialize positions 
+ntraj=30
+q0= np.empty(shape=[ntraj,2])
+for i in range(ntraj):
+    for j in range(2):
+        q0[i][j]=0.000
+
+# Calculate Trajectory
+time_steps=1000
+data = simulator.run(time_steps, q0,save_every=1)  
+
+# Plot the resulting trajectories
+fig, axs = plt.subplots()
+for n, trj in enumerate(data):
+    axs.plot(trj["x"][:,0],trj["x"][:,1])
+    axs.spines['left'].set_position('center')
+    axs.spines['right'].set_color('none')
+    axs.spines['bottom'].set_position('center')
+    axs.spines['top'].set_color('none')
+    axs.xaxis.set_ticks_position('bottom')
+    axs.yaxis.set_ticks_position('left')
+    axs.set_xlabel("$X(t)$")
+    axs.set_ylabel("$Y(t)$")
+    axs.set_title("X-Y Trajectory")
+    axs.grid()
+
+# plot Trajectories 
+fig,bb =  plt.subplots(1,2)
+for n, trj in enumerate(data):
+    bb[0].plot(trj["x"][:,0])
+    bb[1].plot(trj["x"][:,1])
+
+# Set visible  axis
+    bb[0].spines['right'].set_color('none')
+    bb[0].spines['bottom'].set_position('center')
+    bb[0].spines['top'].set_color('none')
+    bb[0].xaxis.set_ticks_position('bottom')
+    bb[0].yaxis.set_ticks_position('left')
+    bb[0].set_xlabel("$timestep$")
+    bb[0].set_ylabel("$X(t)$")
+
+# Set visible axis
+    bb[1].spines['right'].set_color('none')
+    bb[1].spines['bottom'].set_position('center')
+    bb[1].spines['top'].set_color('none')
+    bb[1].xaxis.set_ticks_position('bottom')
+    bb[1].yaxis.set_ticks_position('left')
+    bb[1].set_xlabel("$timestep$")
+    bb[1].set_ylabel("$Y(t)$")
+
+    bb[0].set_title("X Dynamics")
+    bb[1].set_title("Y Dynamics")
+
+
+#########################################
+#  PROJECTION ALONG CHOSEN COORDINATE  #
+#########################################
+
+# Choose unit versor of direction 
+u = np.array([1,1])
+u_norm= (1/np.linalg.norm(u,2))*u
+w = np.empty_like(trj["x"][:,0])
+proj_data = fl.Trajectories(dt= 1e-3)
+fig, axs =plt.subplots()
+for n, trj in enumerate(data):
+    for i in range(len(trj["x"])):
+        w[i]=np.dot(trj["x"][i],u_norm)
+    proj_data.append(fl.Trajectory(1e-3,deepcopy(w.reshape(len(trj["x"][:,0]),1))))
+    axs.plot(proj_data[n]["x"])
+    axs.set_xlabel("$timesteps$")
+    axs.set_ylabel("$w(t)$")
+    axs.set_title("trajectory projected along $u =$"  + str(u) + " direction")
+    axs.grid()
+
+plt.show()