update on machine learning

doc/pub/How2ReadData/html/src/Hudson_Bay.csv

@@ -0,0 +1,22 @@
+Year,Hares (x1000),Lynx (x1000)
+1900,30.0,4.0
+1901,47.2,6.1
+1902,70.2,9.8
+1903,77.4,35.2
+1904,36.3,59.4
+1905,20.6,41.7
+1906,18.1,19.0
+1907,21.4,13.0
+1908,22.0,8.3
+1909,25.4,9.1
+1910,27.1,7.4
+1911,40.3,8.0
+1912,57,12.3
+1913,76.6,19.5
+1914,52.3,45.7
+1915,19.5,51.1
+1916,11.2,29.7
+1917,7.6,15.8
+1918,14.6,9.7
+1919,16.2,10.1
+1920,24.7,8.6

doc/pub/How2ReadData/html/src/Hudson_Bay.py

@@ -0,0 +1,43 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def solver(m, H0, L0, dt, a, b, c, d, t0):
+    """Solve the difference equations for H and L over m years
+    with time step dt (measured in years."""
+
+    num_intervals = int(m/float(dt))
+    t = np.linspace(t0, t0 + m, num_intervals+1)
+    H = np.zeros(t.size)
+    L = np.zeros(t.size)
+
+    print('Init:', H0, L0, dt)
+    H[0] = H0
+    L[0] = L0
+
+    for n in range(0, len(t)-1):
+        H[n+1] = H[n] + a*dt*H[n] - b*dt*H[n]*L[n]
+        L[n+1] = L[n] + d*dt*H[n]*L[n] - c*dt*L[n]
+    return H, L, t
+
+# Load in data file
+data = np.loadtxt('src/Hudson_Bay.csv', delimiter=',', skiprows=1)
+# Make arrays containing x-axis and hares and lynx populations
+t_e = data[:,0]
+H_e = data[:,1]
+L_e = data[:,2]
+
+# Simulate using the model
+H, L, t = solver(m=20, H0=34.91, L0=3.857, dt=0.1,
+                 a=0.4807, b=0.02482, c=0.9272, d=0.02756,
+                 t0=1900)
+
+# Visualize simulations and data
+plt.plot(t_e, H_e, 'b-+', t_e, L_e, 'r-o', t, H, 'm--', t, L, 'k--')
+plt.xlabel('Year')
+plt.ylabel('Numbers of hares and lynx')
+plt.axis([1900, 1920, 0, 140])
+plt.title(r'Population of hares and lynx 1900-1920 (x1000)')
+plt.legend(('H_e', 'L_e', 'H', 'L'), loc='upper left')
+plt.savefig('Hudson_Bay_sim.pdf')
+plt.savefig('Hudson_Bay_sim.png')
+plt.show()

doc/pub/How2ReadData/html/src/diffeq.py

@@ -0,0 +1,12 @@
+import numpy as np
+
+t = np.linspace(0, 10, 21)  # 20 intervals in [0, 10]
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 1
+r = 0.5
+
+for n in range(0, N.size-1, 1):
+    N[n+1] = N[n] + r*dt*N[n]
+    print 'N[%d]=%.1f' % (n+1, N[n+1])

doc/pub/How2ReadData/html/src/ecoli.csv

@@ -0,0 +1,11 @@
+0,100
+600,140
+1200,250
+1800,360
+2400,480
+3000,820
+3600,1300
+4200,1700
+4800,2900
+5400,3900
+6000,7000

doc/pub/How2ReadData/html/src/ecoli.py

@@ -0,0 +1,27 @@
+import numpy as np
+
+# Estimate r
+data = np.loadtxt('ecoli.csv', delimiter=',')
+t_e = data[:,0]
+N_e = data[:,1]
+i = 2  # Data point (i,i+1) used to estimate r
+r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
+print 'Estimated r=%.5f' % r
+# Can experiment with r values and see if the model can
+# match the data better
+
+T = 1200     # cell can divide after T sec
+t_max = 5*T  # 5 generations in experiment
+t = np.linspace(0, t_max, 1000)
+dt = t[1] - t[0]
+N = np.zeros(t.size)
+
+N[0] = 100
+for n in range(0, len(t)-1, 1):
+    N[n+1] = N[n] + r*dt*N[n]
+
+import matplotlib.pyplot as plt
+plt.plot(t, N, 'r-', t_e, N_e, 'bo')
+plt.xlabel('time [s]');  plt.ylabel('N')
+plt.legend(['model', 'experiment'], loc='upper left')
+plt.show()

doc/pub/How2ReadData/html/src/ecoli_inverse.py

@@ -0,0 +1,27 @@
+import numpy as np
+
+data = np.loadtxt('ecoli.csv', delimiter=',')
+t_experiment = data[:,0]
+N_experiment = data[:,1]
+
+def error(p):
+    r = p[0]
+    T = 1200     # cell can divide after T sec
+    t_max = 5*T  # 5 generations in experiment
+    t = np.linspace(0, t_max, len(t_experiment))
+    dt = (t[1] - t[0])
+    N = np.zeros(t.size)
+
+    N[0] = 100
+    for n in range(0, len(t)-1, 1):
+        N[n+1] = N[n] + r*dt*N[n]
+
+    e = np.sqrt(np.sum((N - N_experiment)**2))/N[0]  # error measure
+    e = abs(N[-1] - N_experiment[-1])/N[0]
+    print 'r=', r, 'e=',e
+    return e
+
+from scipy.optimize import minimize
+
+p = minimize(error, [0.0006], tol=1E-5)
+print p

doc/pub/How2ReadData/html/src/plot_Hudson.py

@@ -0,0 +1,19 @@
+import numpy as np
+from  matplotlib import pyplot as plt
+
+# Load in data file
+data = np.loadtxt('src/Hudson_Bay.csv', delimiter=',', skiprows=1)
+# Make arrays containing x-axis and hares and lynx populations
+year = data[:,0]
+hares = data[:,1]
+lynx = data[:,2]
+
+plt.plot(year, hares ,'b-+', year, lynx, 'r-o')
+plt.axis([1900,1920,0, 100.0])
+plt.xlabel(r'Year')
+plt.ylabel(r'Numbers of hares and lynx ')
+plt.legend(('Hares','Lynx'), loc='upper right')
+plt.title(r'Population of hares and lynx from 1900-1920 (x1000)}')
+plt.savefig('Hudson_Bay_data.pdf')
+plt.savefig('Hudson_Bay_data.png')
+plt.show()