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:])