import syde556
import numpy

x = numpy.zeros(100)
for i in range(10):
    x[i*10:(i+1)*10]=numpy.random.uniform(-1,1)

figure()
title('$x$ value changing over time')
plot(x)
ylabel('$x$')
ylim(-1,1)
xlabel('time (ms)')

figure()
title('firing rate changing over time')
plot(syde556.lif(1.5*x+1))
ylabel('firing rate (Hz)')
xlabel('time (ms)')
show()


import syde556
reload(syde556)

import numpy

x = numpy.zeros(100)

v, spikes = syde556.lif_spikes(1.5*x+2)    
    
figure()
title('voltage and spikes over time')
plot(v+spikes)
vlines(numpy.where(spikes>0)[0], 0, 2)
ylabel('firing rate (Hz)')
xlabel('time (ms)')
show()


import syde556
reload(syde556)

import numpy

x = numpy.zeros(100)
for i in range(10):
    x[i*10:(i+1)*10]=numpy.random.uniform(-1,1)

v, spikes = syde556.lif_spikes(1.5*x+1)    
    
figure()
title('$x$ value changing over time')
plot(x)
ylabel('$x$')
ylim(-1,1)
xlabel('time (ms)')

figure()
title('firing rate changing over time')
plot(syde556.lif(1.5*x+1))
ylabel('firing rate (Hz)')
xlabel('time (ms)')
show()


figure()
title('voltage and spikes over time')
plot(v+spikes)
vlines(numpy.where(spikes>0)[0], 0, 2)
ylabel('voltage')
xlabel('time (ms)')
show()



import syde556
reload(syde556)

import numpy

dt = 0.001
t = numpy.arange(600)*dt
x = numpy.sin(10*t)

alpha = 1.5
bias = 2
J = alpha*x+bias

v, spikes = syde556.lif_spikes(J)    
    
figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, v+spikes, label='v')
vlines(dt*numpy.where(spikes>0)[0], 0, 2)
ylabel('$x$')
ylim(-1,2)
xlabel('time (s)')
legend(loc='lower left')

show()


import syde556

import numpy

dt = 0.001
t = numpy.arange(600)*dt
x = numpy.sin(10*t)

alpha = 1.5
bias = 2
J1 = alpha*x+bias
J2 = -alpha*x+bias

v1, spikes1 = syde556.lif_spikes(J1)    
v2, spikes2 = syde556.lif_spikes(J2)    
    
figure(figsize=(8,4))
plot(t, x, label='x')
vlines(dt*numpy.where(spikes1>0)[0], 0, 1)
vlines(dt*numpy.where(spikes2>0)[0], -1, 0)
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')

show()

import syde556

import numpy

dt = 0.001
t = numpy.arange(600)*dt
x = numpy.sin(10*t)

alpha = 1.5
bias = 2
J1 = alpha*x+bias
J2 = -alpha*x+bias

v1, spikes1 = syde556.lif_spikes(J1)    
v2, spikes2 = syde556.lif_spikes(J2)

A = numpy.array([spikes1, spikes2]).T
d = syde556.decoders(A, x)

xhat = numpy.dot(A, d)

    
figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, xhat, label='x')
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')

show()

import syde556

import numpy

N = 50
dt = 0.001
t = numpy.arange(600)*dt
x = numpy.sin(10*t)

encoders = numpy.random.choice([-1, 1], N)
intercepts = numpy.random.uniform(-0.9, 0.9, N)
max_rates = numpy.random.uniform(100, 200, N)
alpha, bias = syde556.find_alpha_and_bias(intercepts, max_rates)

A = syde556.activity_spikes(x, encoders, alpha, bias)

d = syde556.decoders(A, x)

xhat = numpy.dot(A, d)

    
figure(figsize=(8,4))
imshow(A.T, extent=(0,0.6,-1,1), aspect='auto', interpolation='none', cmap='binary')
plot(t, x, label='x')
plot(t, xhat, label='x')
ylabel('$x$')
ylim(-1,1);
xlabel('time (s)')

show()

import syde556
import numpy

dt = 0.001
sigma = 0.007
t_h = numpy.arange(200)*dt-0.1
h = numpy.exp(-t_h**2/(2*sigma*sigma))

figure()
plot(t_h, h)
xlabel('t')
ylabel('h(t)')
show()

t = numpy.arange(600)*dt
x = numpy.sin(10*t)

alpha = 1.5
bias = 2
J1 = alpha*x+bias
J2 = -alpha*x+bias

v1, spikes1 = syde556.lif_spikes(J1)    
v2, spikes2 = syde556.lif_spikes(J2)

fspikes1 = numpy.convolve(spikes1, h, mode='same')
fspikes2 = numpy.convolve(spikes2, h, mode='same')


A = numpy.array([fspikes1, fspikes2]).T


d = syde556.decoders(A, x)

xhat = numpy.dot(A, d)

    
figure(figsize=(8,4))
plot(t, x)
vlines(dt*numpy.where(spikes1>0)[0], 0, 1)
vlines(dt*numpy.where(spikes2>0)[0], -1, 0)
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')

figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, fspikes1*d[0])
plot(t, fspikes2*d[1])
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')

figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, xhat, label='x')
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')

show()


import numpy
x = numpy.random.uniform(-1,1,500)
plot(numpy.arange(500)*0.001, x)
show()

import numpy
import syde556
x, X = syde556.generate_signal(T=1.0, dt=0.001, rms=0.5, limit=10, seed=4)
plot(numpy.arange(1000)*0.001, x)
show()

import numpy
import syde556
T = 1.0
dt = 0.001
x, X = syde556.generate_signal(T, dt, rms=0.5, limit=10, seed=3)
freq = numpy.arange(1000)/T
freq -= (T/dt)/2

figure()
plot(freq, np.abs(X))
xlim(-500, 500)
xlabel('freq (Hz)')
ylabel('$|X(\omega)|$')

figure()
plot(freq, np.abs(X))
xlim(-50, 50)
xlabel('freq (Hz)')
ylabel('$|X(\omega)|$')

figure()
plot(numpy.arange(1000)*0.001, x)
xlabel('time (s)')
ylabel('$x(t)$')
show()



import numpy
import syde556
T = 1.0
dt = 0.001
x, X = syde556.generate_signal(T, dt, rms=0.3, limit=10, seed=3)
t = numpy.arange(int(T/dt))*dt
freq = numpy.arange(int(T/dt))/T - (T/dt)/2
omega = freq*2*numpy.pi

plot(t,x)
xlabel('t')
ylabel('$x$')
show()

alpha = 1.5
bias = 2
J1 = alpha*x+bias
J2 = -alpha*x+bias

v1, spikes1 = syde556.lif_spikes(J1)    
v2, spikes2 = syde556.lif_spikes(J2)

r = spikes1 - spikes2

plot(t, r)
plot(t,x)
xlabel('t')
ylabel('$r$')
show()

R = np.fft.fftshift(np.fft.fft(r))

figure()
plot(omega, np.abs(X))
xlabel('$\omega$')
ylabel('$|X(\omega)|$')

figure()
plot(omega, np.abs(R))
xlabel('$\omega$')
ylabel('$|R(\omega)|$')


show()

H = (X*R.conjugate()) / (R*R.conjugate())
h = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(H))).real

figure()
plot(omega, np.abs(H))
xlabel('$\omega$')
ylabel('$|H(\omega)|$')
xlim(-500, 500)

figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')
show()


fspikes1 = numpy.convolve(spikes1, h, mode='same')
fspikes2 = numpy.convolve(spikes2, h, mode='same')

A = numpy.array([fspikes1, fspikes2]).T

d = syde556.decoders(A, x)

xhat = numpy.dot(A, d)

figure()
plot(t, x)
plot(t, xhat)
show()
print 'RMSE:', numpy.average((x-xhat)**2)

y, Y = syde556.generate_signal(T, dt, rms=0.3, limit=10, seed=1)

J1y = alpha*y+bias
J2y = -alpha*y+bias

v1y, spikes1y = syde556.lif_spikes(J1y)    
v2y, spikes2y = syde556.lif_spikes(J2y)

fspikes1y = numpy.convolve(spikes1y, h, mode='same')
fspikes2y = numpy.convolve(spikes2y, h, mode='same')

Ay = numpy.array([fspikes1y, fspikes2y]).T

yhat = numpy.dot(Ay, d)

figure()
plot(t, y)
plot(t, yhat)
show()
print 'RMSE:', numpy.average((y-yhat)**2)

x2 = numpy.zeros_like(x)
r2 = numpy.zeros_like(r)

x2[400:600] = x[100:300]
r2[400:600] = r[100:300]

figure()
plot(t, x)
plot(t, r)
vlines(0.1, -1.5, 1.5, linewidth=3)
vlines(0.3, -1.5, 1.5, linewidth=3)

figure()
plot(t, x2)
plot(t, r2)
show()



sigma = 0.1
window = np.exp(-(t-0.5)**2/(sigma**2))

x2w = numpy.roll(x, 300)*window
r2w = numpy.roll(r, 300)*window

figure()
plot(t, x)
plot(t, r)
vlines(0.1, -1.5, 1.5, linewidth=3)
vlines(0.3, -1.5, 1.5, linewidth=3)

figure()
plot(t, x2w)
plot(t, r2w)
show()

# original code
H = (X*R.conjugate()) / (R*R.conjugate())

# new code
sigma_t = 0.025
W2 = np.exp(-omega**2*sigma_t**2)
H = (np.convolve(X*R.conjugate(),W2, 'same')) / (np.convolve(R*R.conjugate(),W2, 'same'))

import numpy
import syde556
T = 1.0
dt = 0.001

x, X = syde556.generate_signal(T, dt, rms=0.3, limit=10, seed=3)
y, Y = syde556.generate_signal(T, dt, rms=0.3, limit=10, seed=1)

t = numpy.arange(int(T/dt))*dt
freq = numpy.arange(int(T/dt))/T - (T/dt)/2
omega = freq*2*numpy.pi

alpha = 1.5
bias = 2
J1 = alpha*x+bias
J2 = -alpha*x+bias
J1y = alpha*y+bias
J2y = -alpha*y+bias


v1, spikes1 = syde556.lif_spikes(J1)    
v2, spikes2 = syde556.lif_spikes(J2)
v1y, spikes1y = syde556.lif_spikes(J1y)    
v2y, spikes2y = syde556.lif_spikes(J2y)


r = spikes1 - spikes2
R = np.fft.fftshift(np.fft.fft(r))

sigma_t = 0.025
W2 = np.exp(-omega**2*sigma_t**2)
H = (np.convolve(X*R.conjugate(),W2, 'same')) / (np.convolve(R*R.conjugate(),W2, 'same'))
h = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(H))).real


fspikes1 = numpy.convolve(spikes1, h, mode='same')
fspikes2 = numpy.convolve(spikes2, h, mode='same')
fspikes1y = numpy.convolve(spikes1y, h, mode='same')
fspikes2y = numpy.convolve(spikes2y, h, mode='same')

A = numpy.array([fspikes1, fspikes2]).T
Ay = numpy.array([fspikes1y, fspikes2y]).T

d = syde556.decoders(A, x)

xhat = numpy.dot(A, d)
yhat = numpy.dot(Ay, d)


figure()
plot(t, x)
plot(t, xhat)
title('original signal')
show()
print 'RMSE:', numpy.average((x-xhat)**2)

figure()
title('new testing signal')
plot(t, y)
plot(t, yhat)
show()
print 'RMSE:', numpy.average((y-yhat)**2)

z, Z = syde556.generate_signal(T, dt, rms=0.3, limit=2, seed=1)

J1z = alpha*z+bias
J2z = -alpha*z+bias

v1z, spikes1z = syde556.lif_spikes(J1z)    
v2z, spikes2z = syde556.lif_spikes(J2z)

fspikes1z = numpy.convolve(spikes1z, h, mode='same')
fspikes2z = numpy.convolve(spikes2z, h, mode='same')

Az = numpy.array([fspikes1z, fspikes2z]).T

zhat = numpy.dot(Az, d)

figure()
plot(t, z)
plot(t, zhat)
show()
print 'RMSE:', numpy.average((z-zhat)**2)

figure()
plot(omega, np.abs(H))
xlabel('$\omega$')
ylabel('$|H(\omega)|$')
xlim(-500, 500)

figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')

figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')
xlim(-0.1, 0.1)

show()

h_optimal = h
H_optimal = H
xhat_optimal = xhat

figure(figsize=(8,4))
vlines(dt*numpy.where(spikes1>0)[0], -1, 1)
plot(t, fspikes1*d[0], label='filtered spikes')
ylabel('$x$')
ylim(-1,1)
xlim(0.2, 0.4)
xlabel('time (s)')
legend()
show()


import syde556
import numpy

dt = 0.001
tau = 0.05
t_h = numpy.arange(1000)*dt-0.5
h = numpy.exp(-t_h/tau)
h[numpy.where(t_h<0)]=0

figure()
plot(t_h, h)
xlabel('t')
ylabel('h(t)')
show()

t = numpy.arange(1000)*dt
x, X = syde556.generate_signal(1.0, dt, rms=0.3, limit=10, seed=3)


alpha = 1.5
bias = 2
J1 = alpha*x+bias
J2 = -alpha*x+bias

v1, spikes1 = syde556.lif_spikes(J1)    
v2, spikes2 = syde556.lif_spikes(J2)

fspikes1 = numpy.convolve(spikes1, h, mode='same')
fspikes2 = numpy.convolve(spikes2, h, mode='same')


A = numpy.array([fspikes1, fspikes2]).T


d = syde556.decoders(A, x)

xhat = numpy.dot(A, d)

    
figure(figsize=(8,4))
plot(t, x)
vlines(dt*numpy.where(spikes1>0)[0], 0, 1)
vlines(dt*numpy.where(spikes2>0)[0], -1, 0)
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')

figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, fspikes1*d[0])
plot(t, fspikes2*d[1])
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')

figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, xhat, label='x')
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')

show()


import syde556

import numpy

N = 50
tau = 0.01
dt = 0.001
T = 1.0
t = numpy.arange(T/dt)*dt
x, X = syde556.generate_signal(T, dt, rms=0.3, limit=10, seed=3)

encoders = numpy.random.choice([-1, 1], N)
intercepts = numpy.random.uniform(-0.9, 0.9, N)
max_rates = numpy.random.uniform(100, 200, N)
alpha, bias = syde556.find_alpha_and_bias(intercepts, max_rates)

A = syde556.activity_spikes(x, encoders, alpha, bias)
h = syde556.psc_filter(t-T/2, tau)

Af = syde556.apply_filter(A, h)
d = syde556.decoders(Af, x, dx=1.0/1000, noise=0.1)

xhat = numpy.dot(Af, d)

    
figure(figsize=(8,4))
#imshow(A.T, extent=(0,T,-1,1), aspect='auto', interpolation='none', cmap='binary')
plot(t, x, label='$x$')
plot(t, xhat, label='$\hat{x}$')
ylabel('$x$')
ylim(-1,1)
legend()
xlabel('time (s)')
show()

import syde556
reload(syde556)

import numpy

N = 50
tau = 0.01
dt = 0.001
T = 1.0
t = numpy.arange(T/dt)*dt
x, X = syde556.generate_signal(T, dt, rms=0.3, limit=10, seed=3)

encoders = numpy.random.choice([-1, 1], N)
intercepts = numpy.random.uniform(-0.9, 0.9, N)
max_rates = numpy.random.uniform(100, 200, N)
alpha, bias = syde556.find_alpha_and_bias(intercepts, max_rates)

x_rate = numpy.linspace(-1,1, 1000)
A_rate = syde556.activity(x_rate, encoders, alpha, bias)
d = syde556.decoders(A_rate, x_rate, dx=1.0/1000, noise=0.1)

A = syde556.activity_spikes(x, encoders, alpha, bias)

h = syde556.psc_filter(t-T/2, tau)

Af = syde556.apply_filter(A, h)

xhat = numpy.dot(Af, d)

    
figure(figsize=(8,4))
#imshow(A.T, extent=(0,T,-1,1), aspect='auto', interpolation='none', cmap='binary')
plot(t, x, label='$x$')
plot(t, xhat, label='$\hat{x}$')
ylabel('$x$')
ylim(-1,1)
legend()
xlabel('time (s)')
show()