add python examples

Author: Nico-Curti

python/bernoulli2D.py

@@ -0,0 +1,38 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np # numerical library
+import matplotlib.pylab as plt # plot library
+import matplotlib.animation as animation # animation plot
+
+#%% Bernoulli2D
+bernoulli2D = lambda x : np.mod(2*x, 1) # bernoulli formula
+
+# initial condition
+N = 100 # number of points
+x = np.linspace(0, 1, N) # x steps
+y = np.linspace(0, 1, N) # y steps
+meanx, stdx = np.mean(x), np.std(x)
+meany, stdy = np.mean(y), np.std(y)
+x, y = np.meshgrid(x, y)
+G = np.exp( - ( .5*(x-meanx)**2 / stdx**2 + .5*(y-meany)**2 / stdy**2 ) )
+
+
+fig = plt.figure(figsize=(8,8))
+time = 100 # number of iterations
+ims = np.empty(time - 1, dtype=np.object)
+ev0 = G
+for i in range(1, time):
+    ev = bernoulli2D(ev0)
+    ims[i-1] = [plt.imshow( ev, animated=True )]
+    ev0 = ev
+movie = animation.ArtistAnimation(fig,
+        								  ims,
+                                   interval=50, 
+                                   blit=True,
+                                   repeat_delay=100
+                                  )

python/bifurcation.py

@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np # numerical library
+import matplotlib.pylab as plt # plot library
+
+# logistic map
+logistic = lambda x, mu : 4*mu*x*(1-x)
+
+# Bifurcation logistic map x_{n+1} = \mu x_n (1 - x_n)
+it = 100000
+mu = np.linspace(.5, 1, it) # coefficient values
+x = np.empty(it)
+x[0] = .5
+
+for t in range(0, it-1):
+	x[t+1] = logistic(x[t], mu[i])
+
+critical_y = x[np.where(np.diff(x) > 0)[0][0]]
+critical_x = mu[critical_y]
+
+fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8,8))
+ax.scatter(mu, x, marker=',', color='k')
+ax.vlines(critical_x, 0, 1, color='r', linestyle="dashed")
+ax.set_xlabel("$\mu$", fontsize=14)
+ax.set_ylabel("$x_k$", fontsize=14)
+ax.text(critical_y, critical_x, "$\mu = %.3f$"critical[0])
+ax.set_xlim(.5, 1.01)

python/bisection_root.py

@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np # numerical library
+import matplotlib.pylab as plt # plot library
+
+#%% Bisection method
+
+f = lambda x : x**2 - 2
+bisection = lambda x0, x1 : (x0, x1) if f(x0)*f(x1) < 0 else (x1, x0)
+
+n = 10000
+toll = 1e-4
+x0 = 5
+x1 = 1
+x = np.linspace(start = 0, stop = 10, num = n)
+y = f(x)
+step = [x0]
+tst = 1
+while tst >= toll:
+	c = (x0 - x1) * .5 # middle point
+	x0, x1 = bisection(c, x0) # bisection update
+	step.append(c) # array of steps
+	tst = abs(step[-1] - step[-2]) # check convergence
+
+fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8,8))
+ax.plot(x, y, 'b-', label="f(x)") # plot function
+ax.plot(step, f(np.asarray(step)), 'r--', label="bisection step") # plot steps
+ax.plot(step[-1], f(step[-1]), 'go') # plot root
+ax.set_xlabel("x", fontsize=14)
+ax.set_ylabel("f(x)", fontsize=14)
+ax.set_title("Bisection method", fontsize=14)

python/duffin_van_der_pol.py

@@ -0,0 +1,30 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np # numerical library
+import matplotlib.pylab as plt # plot library
+from mpl_toolkits.mplot3d import Axes3D # 3D plot
+
+# Duffin-Van Der Pol oscillator
+dv_der_pol = lambda x, y, z, mu, f, o, dt: (dt*y, dt*(mu*1-x*x)*y - x*x*x + f*np.cos(z), z + dt*o)
+
+mu    = .2
+f     = 1.
+omega = .94
+dt    = 1e-3
+it    = 10000
+
+x, y, z = np.empty(it), np.empty(it), np.empty(it)
+x[0], y[0], z[0] = 1., 0., 1. # initial conditions
+time = np.arange(0, dt*it, it) # time
+
+for t in range(0, it-1):
+	x[t+1], y[t+1], z[t+1] = dv_der_pol(x[t], y[t], z[t], mu, f, o, dt)
+
+fig = plt.figure(figsize=(8,8))
+ax = Axes3D(fig)
+ax.plot(x, y, np.sin(z))

python/false_nearest_neighbours.py

@@ -0,0 +1,63 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np
+import matplotlib.pylab as plt
+import scipy
+
+dt = .01
+it = 10000
+sigma = 16.
+b = 4.
+r = 45.92
+x, y, z = np.empty(shape=(it)), np.empty(shape=(it)), np.empty(shape=(it))
+x[0], y[0], z[0] = 10, 1, 1
+
+Vx = lambda x, y, sigma : sigma*(y-x)
+Vy = lambda x, y, z, r : -x*z + r*x - y
+Vz = lambda x, y, z, b : -b*z + x*y
+
+for i in range(0, it-1):
+	x[i+1] = x[i] + dt*Vx(x[i], y[i], sigma)
+	y[i+1] = y[i] + dt*Vy(x[i], y[i], z[i], r)
+	z[i+1] = z[i] + dt*Vz(x[i], y[i], z[i], b)
+
+# False Nearest Neighbours
+
+RT = 15
+AT = 2
+sigmay = np.std(x)
+
+maxEmbDim = 10
+delay= 1
+rEEM = it - (maxEmbDim*delay - delay)
+EEM = np.concatenate([y[delay*maxEmbDim-delay:].reshape(1, len(y[delay*maxEmbDim-delay:])), 
+                      [y[delay*maxEmbDim-(i+1)*delay:-i*delay] for i in range(1, maxEmbDim)]
+                      ], axis=0).T
+Ind1 = np.empty(maxEmbDim)
+Ind2 = np.empty(maxEmbDim)
+embedm = 0 # only for plot
+for k in range(1, maxEmbDim+1):
+    D = scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(EEM[:, :k], "euclidean"))        
+    np.fill_diagonal(a=D, val=np.inf)
+    l = np.argmin(D[:rEEM - maxEmbDim - k, :], axis=1)
+    fnn1 = np.asarray([abs(y[i + maxEmbDim + k - 1]-y[li + maxEmbDim + k - 1])/D[i, li] for i,li in enumerate(l) if D[i, li] > 0 and li + maxEmbDim + k - 1 < it])
+    fnn2 = np.asarray([abs(y[i + maxEmbDim + k - 1]-y[li + maxEmbDim + k - 1])/sigmay for i,li in enumerate(l) if D[i, li] > 0 and li + maxEmbDim + k - 1 < it])
+    Ind1[k-1] = len(np.where(np.asarray(fnn1) > RT)[0])
+    Ind2[k-1] = len(np.where(np.asarray(fnn2) > AT)[0])
+    
+    if embedm == 0: # only for plot
+	    if Ind1[k-1] / len(fnn1) < .1 and Ind2[-1] / len(fnn1) < .1 and Ind1[k-1] != 0:
+	        embedm = k
+	        #break # uncomment for true algorithm
+fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8,8))
+ax.plot(np.arange(0, maxEmbDim), Ind1)
+ax.set_xlabel('Embedding dimension', fontsize=14)
+ax.set_ylabel('% FNN', fontsize=14)
+ax.set_title('Optimal Embedding Dimension with FNN', fontsize=16)
+ax.plot(embedm, Ind1[embedm], 'r.') 
+plt.text(embedm, Ind1[embedm] + 100, "EmbDim = $%d$"%(embedm))

python/grassberger_procaccia.py

@@ -0,0 +1,53 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np
+import matplotlib.pylab as plt
+import scipy
+import pandas as pd
+import seaborn as sns
+
+dt = .01
+it = 10000
+sigma = 16.
+b = 4.
+r = 45.92
+x, y, z = np.empty(shape=(it)), np.empty(shape=(it)), np.empty(shape=(it))
+x[0], y[0], z[0] = 10, 1, 1
+
+Vx = lambda x, y, sigma : sigma*(y-x)
+Vy = lambda x, y, z, r : -x*z + r*x - y
+Vz = lambda x, y, z, b : -b*z + x*y
+
+for i in range(0, it-1):
+	x[i+1] = x[i] + dt*Vx(x[i], y[i], sigma)
+	y[i+1] = y[i] + dt*Vy(x[i], y[i], z[i], r)
+	z[i+1] = z[i] + dt*Vz(x[i], y[i], z[i], b)
+
+# Grassberger Procaccia
+
+omit_pts = 3
+ED = scipy.spatial.distance.pdist(np.concatenate([x.reshape(len(x), 1), x.reshape(len(x),1)]), "euclidean")
+min_eps = ED[ED!=0].min()
+max_eps = 2**np.ceil(np.log(ED.max())/np.log(2))
+eps_vec_size = int(np.floor( (np.log(max_eps/min_eps) / np.log(2)))) + 1
+eps_vec = max_eps*2.**(-np.arange(0, eps_vec_size))
+
+Npairs = it*(it-1)*.5
+C_eps = np.asarray([len(ED[(ED < eps) & (ED > 0)])/Npairs for eps in eps_vec])
+
+k1 = omit_pts + 1
+k2 = eps_vec_size - omit_pts
+xp = np.log(eps_vec[np.arange(k1, k2+1)]) / np.log(2)
+yp = np.log(C_eps[np.arange(k1, k2+1)]) / np.log(2)
+
+data = pd.DataFrame(data=[xp, yp], index=["xp", "yp"]).T
+fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8,8))
+plot = sns.regplot(x="xp",
+                   y="yp",
+                   data=data,
+                   ax=ax).legend()

python/harmonic.py

@@ -0,0 +1,40 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np # numerical library
+import matplotlib.pylab as plt # plot library
+
+# Harmonic oscillator
+Vx = lambda q, p : p
+Vy = lambda q, p : -q
+
+# Integrators
+euler = lambda q, p, dt, k : ( q + dt*Vx(q, p), q + dt*k*Vy(q, p) )
+simplettic   = lambda q, p, dt, k : ( q + dt+Vx(q, p + dt*k*Vy(q, p)), p + dt*k*Vy(q, p))
+
+
+k = .5
+it = 1000
+dt = .1
+q_eu, p_eu = np.empty(it), np.empty(it)
+q_eu[0], p_eu[0] = .5, 0 # initial conditions
+q_si, p_si = np.empty(it), np.empty(it)
+q_si[0], p_si[0] = .5, 0 # initial conditions
+
+for t in range(0, it-1):
+	q_eu[t+1], p_eu[t+1] = euler(q_eu[t], p_eu[t], dt, k)
+	q_si[t+1], p_si[t+1] = rk2(q_si[t], p_si[t], dt, k)
+
+time = np.arange(0, it*dt, dt)
+fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, figsize=(8,8))
+ax1.plot(time, q_eu, 'b-', label="q euler")
+ax1.plot(time, p_eu, 'r-', label="p euler")
+ax2.plot(q_eu, p_eu, 'b-', label="phase space")
+
+ax3.plot(time, q_si, 'b-', label="q simplettic")
+ax3.plot(time, p_si, 'r-', label="p simplettic")
+ax4.plot(q_si, p_si, 'b-', label="phase space")

python/integration.py

@@ -0,0 +1,22 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np # numerical library
+
+#%% Integration
+
+f = lambda x : np.sin(x)
+rettangoli = lambda dx, N : sum( [dx * f(i * dx) for i in range(N)] )
+trapezi = lambda dx, N : sum( [dx * (f(i*dx) + f((i+1)*dx)) *.5 for i in range(N)] )
+simpson = lambda dx, N : sum( [dx/6*(f(i*dx) + 4*f((i*dx + (i+1)*dx)*.5) + f((i+1)*dx)) for i in range(N)] )
+
+truth = 1. - np.cos(1.) # sin primitive in [0, 1]
+dx = 1.
+for n in range(1, 100):
+	N = n*n
+	print( N, abs(rettangoli(dx / N, N) - truth), abs(trapezi(dx / N, N) - truth), abs(simpson(dx / N, N) - truth) )
+

python/logistic_map.py

@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np # numerical library
+import matplotlib.pylab as plt # plot library
+from scipy.integrate import odeint # python integrator
+
+# Logistic map
+logistic = lambda x, t, mu : mu * (np.mod(x, 1) - np.mod(x**2, 1)) # logistic formula
+
+it = 60
+dt = .1
+x0 = np.mod(.1, 1)
+mu = 1.3
+t = np.arange(0, dt*it, dt)
+x = np.empty(shape=(it))
+x[0] = x0
+
+for i in range(0, it-1):
+	x[i+1] = np.mod( np.mod(x[i], 1) + dt*logistic(x[i], t[i], mu) , 1)
+
+xt = odeint(logistic, x0, t, args=(mu, ))
+
+fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8,8))
+ax.scatter(t, x, marker='o', facecolors="none", edgecolor='k', s=20, label="euler")
+ax.plot(t, xt, 'r-', lw=2, alpha=.5, label="scipy")
+ax.set_ylim(0, 1)
+ax.legend(loc='best', fontsize=14)

python/lorentz_attractor.py

@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 22:30:58 2018
+
+@author: NICO
+"""
+
+import numpy as np # numerical library
+import matplotlib.pylab as plt # plot library
+from mpl_toolkits.mplot3d import Axes3D # 3D plot
+
+# Lorentz attractor
+
+dt = .01
+it = 10000
+sigma = 16.
+b = 4.
+r = 45.92
+x, y, z = np.empty(shape=(it)), np.empty(shape=(it)), np.empty(shape=(it))
+x[0], y[0], z[0] = 10, 1, 1 # initial conditions
+
+# lorentz formula
+Vx = lambda x, y, sigma : sigma*(y-x)
+Vy = lambda x, y, z, r : -x*z + r*x - y
+Vz = lambda x, y, z, b : -b*z + x*y
+
+# euler integration
+for i in range(0, it-1):
+	x[i+1] = x[i] + dt*Vx(x[i], y[i], sigma)
+	y[i+1] = y[i] + dt*Vy(x[i], y[i], z[i], r)
+	z[i+1] = z[i] + dt*Vz(x[i], y[i], z[i], b)
+
+fig = plt.figure(figsize=(8,8))
+ax = Axes3D(fig)
+ax.plot(x, y, z)