import Numeric as num
import math, random
from scipy import gplt
A = num.array([[0.1,0.],[0.,0.1]])
omega = num.array([0.04, 0.01])
c = num.array([0.0, 0.0]) # num.matrixmultiply(A,num.array([1,1]))
w = num.array([1., 0.0])
r = num.array([0.0, 0.0])
S_in = num.array([0.0, 0.0])
S_out = 1.
v = 0.
tau = 350.
eta = 0.0001
p = num.array([0.0, 0.0])
d = 0.
a_p = 0.0001
a_d = 0.0001
tau_p = 350.
tau_d = 350.
dp = num.array([0.0, 0.0])
dw = num.array([0.0, 0.0])
its = 350000
data = range(its)
for t in range(its):
r = num.matrixmultiply(A,[math.sin(omega[0]*t),math.sin(omega[1]*t)])+c
r = r.astype('f')
for i in range(2):
S_in[i] = float(random.random() <= (1 - math.exp(-r[i])))
dv = eta*num.matrixmultiply(w,num.transpose(S_in)) - (v/tau)
v += dv
S_out = float(random.random() <= (1 - math.exp(-v-2)))
dd = S_out - d/tau_d
for i in range(2):
dp[i] = S_in[i] - p[i]/tau_p
d += dd
p += dp
if t%1000 == 0:
print w
for i in range(2):
dw[i] = (a_p*S_out*p[i] - a_d*S_in[i]*d)
w += dw
w = w * (1/math.sqrt(w[0]**2+w[1]**2))
data[t] = [w[0],w[1],v+2]
print w
gplt.plot(range(100000,its), data[100000:])