But when you see things, you get more than one input at a time
What happens if I see a dot at $<h_0, v_0>$ and another dot at $<h_1, v_1>$ at the same time?
Is that the same as seeing one dot at $<h_0 + h_1, v_0 + v_1>$?
We see them as two distinct points
Same thing with auditory (multiple frequencies at once, plus other sound features)
Maybe even the same thing with the superior colliculus and eye movements (although that's less clear)
So what should this representation be?
import numpy
def gabor(width, height, lambd, theta, psi, sigma, gamma, x_offset, y_offset):
x = numpy.linspace(-1, 1, width)
y = numpy.linspace(-1, 1, height)
X, Y = numpy.meshgrid(x, y)
X = X - x_offset
Y = Y + y_offset
cosTheta = numpy.cos(theta)
sinTheta = numpy.sin(theta)
xTheta = X * cosTheta + Y * sinTheta
yTheta = -X * sinTheta + Y * cosTheta
e = numpy.exp( -(xTheta**2 + yTheta**2 * gamma**2) / (2 * sigma**2) )
cos = numpy.cos(2 * numpy.pi * xTheta / lambd + psi)
return e * cos
g = gabor(500, 500, lambd=0.6,
theta=numpy.pi/4,
psi=numpy.pi/2,
sigma=0.2,
gamma=0.7,
x_offset=0,
y_offset=0)
import pylab
pylab.imshow(g, extent=(-1, 1, -1, 1), interpolation='none', vmin=-1, vmax=1, cmap='gray')
pylab.show()
width = 50
height = 50
N = 100
def make_random_gabor(width, height):
return gabor(width, height,
lambd=random.uniform(0.3, 0.8),
theta=random.uniform(0, 2*numpy.pi),
psi=random.uniform(0, 2*numpy.pi),
sigma=random.uniform(0.2, 0.5),
gamma=random.uniform(0.4, 0.8),
x_offset=random.uniform(-1, 1),
y_offset=random.uniform(-1, 1))
encoders = [make_random_gabor(width, height) for i in range(N)]
import pylab
pylab.figure(figsize=(10,8))
for i in range(N):
w = i%12
h = i/12
pylab.imshow(encoders[i], extent=(w, w+0.95, h, h+0.95), interpolation='none',
vmin=-1, vmax=1, cmap='gray')
pylab.xticks([])
pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, N/12))
pylab.show()
S = 100
width = 50
height = 50
x = numpy.random.normal(size=(width*height, S))
import pylab
pylab.figure(figsize=(10,8))
for i in range(S):
w = i%12
h = i/12
img = x[:,i]
img.shape = width, height
pylab.imshow(img, extent=(w, w+0.95, h, h+0.95), interpolation='none',
vmin=-1, vmax=1, cmap='gray')
pylab.xticks([])
pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/12))
pylab.show()
S = 100
width = 50
height = 50
def normalize(v):
return v/numpy.linalg.norm(v)
C = 4
x = [normalize(numpy.sum([make_random_gabor(width, height) for i in range(C)],axis=0).flatten()) for j in range(S)]
x = numpy.array(x).T
import pylab
pylab.figure(figsize=(10,8))
for i in range(S):
w = i%12
h = i/12
img = x[:,i]
img.shape = height, width
pylab.imshow(img, extent=(w, w+0.95, h, h+0.95), interpolation='none',
vmin=-0.1, vmax=0.1, cmap='gray')
pylab.xticks([])
pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/12))
pylab.show()
import syde556
S = 100
x = []
for i in range(S):
sx, A = syde556.generate_signal(T=width, dt=1.0, rms=1, limit=0.05, seed=i)
sy, A = syde556.generate_signal(T=height, dt=1.0, rms=1, limit=0.05, seed=i+2000)
SX, SY = numpy.meshgrid(sx, sy)
img = SX*SY
x.append(normalize(img.flatten()))
x = numpy.array(x).T
import pylab
pylab.figure(figsize=(10,8))
for i in range(S):
w = i%12
h = i/12
img = x[:,i]
img.shape = height, width
pylab.imshow(img, extent=(w, w+0.95, h, h+0.95), interpolation='none',
vmin=-0.1, vmax=0.1, cmap='gray')
pylab.xticks([])
pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/12))
pylab.show()
import syde556
reload(syde556)
N = 100
width = 50
height = 50
encoders = [make_random_gabor(width, height).flatten() for i in range(N)]
vision = syde556.Ensemble(neurons=N, dimensions=width*height, encoders=encoders)
decoder = vision.compute_decoder(x, x, noise=0.1)
A, xhat = vision.simulate_rate(x, decoder)
import pylab
pylab.figure(figsize=(10,12))
for i in range(S):
w = i%6
h = i/6
img = x[:,i]
img.shape = height, width
pylab.imshow(img, extent=(2*w, 2*w+0.9, h, h+0.9), interpolation='none',
vmin=-0.1, vmax=0.1, cmap='gray')
img = xhat[:,i]
img.shape = height, width
pylab.imshow(img, extent=(2*w+0.95, 2*w+1.85, h, h+0.9), interpolation='none',
vmin=-0.1, vmax=0.1, cmap='gray')
pylab.xticks([])
pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/6))
pylab.show()
width = 50
height = 50
# make samples
S = 100
x = []
for i in range(S):
sx, A = syde556.generate_signal(T=width, dt=1.0, rms=1, limit=2.5/width, seed=i)
sy, A = syde556.generate_signal(T=height, dt=1.0, rms=1, limit=2.5/height, seed=i+2000)
SX, SY = numpy.meshgrid(sx, sy)
img = SX*SY
x.append(normalize(img.flatten()))
x = numpy.array(x).T
print x.shape
#make neurons
N = 500
encoders = [make_random_gabor(width, height).flatten() for i in range(N)]
vision = syde556.Ensemble(neurons=N, dimensions=width*height, encoders=encoders)
decoder = vision.compute_decoder(x, x, noise=0.1)
A, xhat = vision.simulate_rate(x, decoder)
scale = 5/np.sqrt(width*height)
#plot results
import pylab
pylab.figure(figsize=(10,12))
for i in range(S):
w = i%6
h = i/6
img = x[:,i]
img.shape = height, width
pylab.imshow(img, extent=(2*w, 2*w+0.9, h, h+0.9), interpolation='none',
vmin=-scale, vmax=scale, cmap='gray')
img = xhat[:,i]
img.shape = height, width
pylab.imshow(img, extent=(2*w+0.95, 2*w+1.85, h, h+0.9), interpolation='none',
vmin=-scale, vmax=scale, cmap='gray')
pylab.xticks([])
pylab.yticks([])
pylab.xlim((0, 12))
pylab.ylim((0, S/6))
pylab.show()
(2500L, 100L)
Is that all we need?
Still have to be able to compute $\int x(v) e(v) dv$
$\int x(v) e(v) dv$
$\int (\sum_i c_{x,i} b_i(v)) (\sum_i c_{e,i} b_i(v)) dv$
$\int (\sum_i c_{x,i} c_{e,i} b_i(v)^2)dv + (\mbox{ ... all the cross terms ... })$
so those cross terms are ugly
Let's assume $b_i(v)$ are orthogonal
$\int (\sum_i c_{x,i} c_{e,i} b_i(v)^2)dv$
now rearrange, pulling the sum out
$\sum c_{x,i} c_{e,i} \int b_i(v)^2 dv$
what if we make the basis functions have unit length?
so $\int x(v) e(v) dv$ becomes $\sum c_{x,i} c_{e,i}$, which is just $c_x \cdot c_e$
the only works assuming the basis vectors $b_i(v)$ are:
this is called an orthonormal basis
import numpy
v = numpy.linspace(-1, 1, 1000)
b_i = v**1
b_j = v**2
import pylab
pylab.plot(v, b_i, label='$b_i$')
pylab.plot(v, b_j, label='$b_j$')
pylab.plot(v, b_i*b_j, label='$b_i b_j$')
pylab.legend(loc='best')
xlabel('$v$')
pylab.show()
dv = 2.0/len(v)
print 'sum:', numpy.sum(b_i*b_j*dv)
sum: 3.90312782095e-18
for i in range(5):
plot(v, numpy.polynomial.legendre.legval(v, numpy.eye(5)[i]))
xlabel('$v$')
show()
import numpy
v = numpy.linspace(-1, 1, 1000)
i = 1
j = 2
b_i = numpy.polynomial.legendre.legval(v, numpy.eye(i+1)[i])
b_j = numpy.polynomial.legendre.legval(v, numpy.eye(j+1)[j])
import pylab
pylab.plot(v, b_i, label='$b_i$', linewidth=3)
pylab.plot(v, b_j, label='$b_j$', linewidth=3)
pylab.plot(v, b_i*b_j, label='$b_i b_j$')
pylab.legend(loc='best')
xlabel('$v$')
pylab.show()
dv = 2.0/len(v)
print 'sum:', numpy.sum(b_i*b_j*dv)
sum: 9.28077059648e-17
v = numpy.linspace(-1, 1, 1000)
dv = 2.0/len(v)
length = []
for i in range(10):
b_i = numpy.polynomial.legendre.legval(v, numpy.eye(i+1)[i])
length.append(sum(b_i*b_i)*dv)
plot(length)
xlabel('$i$')
show()
v = numpy.linspace(-1, 1, 1000)
dv = 2.0/len(v)
length = []
for i in range(10):
b_i = numpy.polynomial.legendre.legval(v, numpy.eye(i+1)[i])
b_i = b_i * numpy.sqrt((2*i + 1)/2.0)
length.append(sum(b_i*b_i)*dv)
plot(length)
xlabel('$i$')
ylim(0, 2)
show()
for i in range(5):
plot(v, numpy.polynomial.legendre.legval(v, numpy.eye(i+1)[i])* numpy.sqrt((2*i + 1)/2.0))
xlabel('$v$')
show()
D = 5
S = 20
x_poly = numpy.random.randn(D, S)/numpy.sqrt(D)
v = numpy.linspace(-1, 1, 1000)
for i in range(S):
plot(v, numpy.polynomial.polynomial.polyval(v, x_poly[:,i]))
xlabel('$v$')
show()
x = []
for i in range(S):
x.append(numpy.polynomial.legendre.poly2leg(x_poly[:,i])/numpy.sqrt((2*i + 1)/2.0))
x = numpy.array(x).T
v = numpy.linspace(-1, 1, 1000)
for i in range(S):
plot(v, numpy.polynomial.legendre.legval(v, x[:,i])*numpy.sqrt((2*i + 1)/2.0))
xlabel('$v$')
show()
N = 50
ensemble = syde556.Ensemble(neurons=N, dimensions=5)
decoder = ensemble.compute_decoder(x, x, noise=0.1)
A, xhat = ensemble.simulate_rate(x, decoder)
figure(figsize=(10,4))
subplot(1, 2, 1)
for i in range(S):
plot(v, numpy.polynomial.legendre.legval(v, x[:,i])*numpy.sqrt((2*i + 1)/2.0))
ylim(-3, 3)
xlabel('$v$')
title('original')
subplot(1, 2, 2)
for i in range(S):
plot(v, numpy.polynomial.legendre.legval(v, xhat[:,i])*numpy.sqrt((2*i + 1)/2.0))
ylim(-3, 3)
xlabel('$v$')
title('decoded')
show()
def area(x):
fv = numpy.polynomial.legendre.legval(v, x)*numpy.sqrt((2*i + 1)/2.0)
a = sum(fv)*dv
return [a]
fx = numpy.array([area(xx) for xx in x.T]).T
decoder = ensemble.compute_decoder(x, fx, noise=0.1)
A, fxhat = ensemble.simulate_rate(x, decoder)
for i in range(10):
print x[:,i], fx[:,i], fxhat[:,i]
[-0.45249468 0.06689621 0.65862321 -0.22781363 0.04283414] [-3.99012265] [-3.692147] [-0.13508849 -0.45085403 -0.29704591 -0.14351425 -0.04635415] [-1.19610606] [-0.99389806] [ 0.03477385 0.59657275 -0.18919625 0.00530379 -0.06779654] [ 0.30484099] [ 0.29829337] [ 0.15631005 -0.07939436 -0.10217138 -0.09528619 -0.02525089] [ 1.37936599] [ 1.31180309] [-0.45944895 0.02760009 -0.27678184 0.08308695 -0.0547291 ] [-4.06067516] [-3.88563593] [ 0.25147024 0.0522454 0.02486886 0.11510253 0.03265593] [ 2.22143421] [ 1.61950481] [ 0.03977521 -0.06903214 0.03796262 0.09403689 0.02458351] [ 0.35183856] [ 0.24418011] [-0.34478934 0.26520947 -0.15590749 0.04113701 -0.00300367] [-3.04650196] [-2.60464404] [-0.08437981 0.28442028 0.19314713 0.05949096 -0.01813703] [-0.74367552] [-0.55534047] [-0.01543183 -0.11761898 0.06478489 0.07507399 -0.00886165] [-0.13579599] [-0.09870697]
def maximum(x):
fv = numpy.polynomial.legendre.legval(v, x)*numpy.sqrt((2*i + 1)/2.0)
index = numpy.argmax(fv)
return [v[index]]
fx = numpy.array([maximum(xx) for xx in x.T]).T
decoder = ensemble.compute_decoder(x, fx, noise=0.01)
A, fxhat = ensemble.simulate_rate(x, decoder)
for i in range(20):
print x[:,i], fx[:,i], fxhat[:,i]
[-0.45249468 0.06689621 0.65862321 -0.22781363 0.04283414] [-1.] [-1.03426409] [-0.13508849 -0.45085403 -0.29704591 -0.14351425 -0.04635415] [-0.94394394] [-0.94790649] [ 0.03477385 0.59657275 -0.18919625 0.00530379 -0.06779654] [ 0.78178178] [ 0.80341805] [ 0.15631005 -0.07939436 -0.10217138 -0.09528619 -0.02525089] [ 0.21721722] [ 0.03715061] [-0.45944895 0.02760009 -0.27678184 0.08308695 -0.0547291 ] [-0.17317317] [-0.14223553] [ 0.25147024 0.0522454 0.02486886 0.11510253 0.03265593] [ 1.] [ 0.98869827] [ 0.03977521 -0.06903214 0.03796262 0.09403689 0.02458351] [ 1.] [ 0.47082208] [-0.34478934 0.26520947 -0.15590749 0.04113701 -0.00300367] [ 1.] [ 1.01930901] [-0.08437981 0.28442028 0.19314713 0.05949096 -0.01813703] [ 1.] [ 0.98959533] [-0.01543183 -0.11761898 0.06478489 0.07507399 -0.00886165] [-0.7997998] [-0.36437875] [-0.147665 0.049168 -0.06538418 -0.00584206 -0.02632643] [ 0.47347347] [ 0.32563004] [ 0.02441503 -0.1699608 0.05092713 -0.02245458 -0.00107661] [-1.] [-0.90374256] [ 0.14030728 0.06652304 0.04019052 0.00865063 0.0332155 ] [ 1.] [ 0.9777674] [-0.04854667 0.19268431 0.12686373 -0.00477267 -0.00366835] [ 1.] [ 0.95878839] [ 0.14449928 0.04155529 0.23186591 0.01521891 0.03536285] [ 1.] [ 0.94144684] [-0.04972338 -0.09333049 0.195611 -0.03597476 0.02794989] [-1.] [-0.93327958] [ 0.0649419 -0.02901023 0.03064615 -0.03862747 -0.00013708] [-1.] [-0.87645674] [ 0.11407232 -0.05903258 0.02590225 -0.0141109 0.01488021] [-1.] [-0.69223594] [ 0.07999296 0.02683351 0.0958427 0.03251879 0.02420397] [ 1.] [ 0.71611585] [ 0.13212116 0.19791227 0.16534859 0.02255834 0.03551369] [ 1.] [ 1.22114539]
def negative(x):
return -x
fx = numpy.array([negative(xx) for xx in x.T]).T
decoder = ensemble.compute_decoder(x, fx, noise=0.01)
A, fxhat = ensemble.simulate_rate(x, decoder)
figure(figsize=(10,4))
subplot(1, 2, 1)
for i in range(S):
plot(v, numpy.polynomial.legendre.legval(v, x[:,i])*numpy.sqrt((2*i + 1)/2.0))
xlabel('$v$')
ylim(-3, 3)
title('original')
subplot(1, 2, 2)
for i in range(S):
plot(v, numpy.polynomial.legendre.legval(v, fxhat[:,i])*numpy.sqrt((2*i + 1)/2.0))
xlabel('$v$')
ylim(-3, 3)
title('decoded')
show()
N = 200
v = numpy.linspace(-1, 1, 1000)
def make_gaussian(centre, sigma):
return numpy.exp(-(v-centre)**2/(2*sigma**2))
encoders = numpy.array([make_gaussian(centre=random.uniform(-1, 1),
sigma=0.1) for i in range(N)])
plot(v, encoders[:20].T)
xlabel('$v$')
show()
U, S, V = numpy.linalg.svd(encoders.T)
plot(v, U[:,:10])
xlabel('$v$')
show()
import numpy
v = numpy.linspace(-1, 1, 1000)
i = 0
j = 1
b_i = U[:,i]
b_j = U[:,j]
import pylab
pylab.plot(v, b_i, label='$b_i$', linewidth=3)
pylab.plot(v, b_j, label='$b_j$', linewidth=3)
pylab.plot(v, b_i*b_j, label='$b_i b_j$')
pylab.legend(loc='best')
xlabel('$v$')
pylab.show()
dv = 2.0/len(v)
print 'sum:', numpy.sum(b_i*b_j*dv)
sum: -2.50635984098e-20
v = numpy.linspace(-1, 1, 1000)
dv = 2.0/len(v)
length = []
for i in range(10):
b_i = U[:,i]
length.append(sum(b_i*b_i)*dv)
plot(length)
xlabel('$i$')
ylim(0, 0.004)
show()
v = numpy.linspace(-1, 1, 1000)
dv = 2.0/len(v)
basis = U/(math.sqrt(dv))
length = []
for i in range(10):
b_i = basis[:,i]
length.append(sum(b_i*b_i)*dv)
plot(length)
xlabel('$i$')
ylim(0, 2)
show()
plot(v, basis[:, :10])
xlabel('$v$')
show()
plot(S)
show()
bases = 20
figure(figsize=(10,4))
subplot(1, 2, 1)
plot(v, encoders.T[:,:20])
ylim(0,1)
title('original')
subplot(1, 2, 2)
basis = U[:,:bases]/(math.sqrt(dv))
e = numpy.dot(encoders, basis)*dv
plot(v, numpy.dot(basis, e.T)[:,:20])
ylim(0,1)
title('bases=%d'%bases)
show()
S = 100
import syde556
x = numpy.array([syde556.generate_signal(2.0, 2.0/1000,
rms=0.5, limit=2, seed=i)[0] for i in range(S)])
plot(v, x[:5].T)
xlabel('$v$')
show()
figure(figsize=(10,4))
subplot(1, 2, 1)
plot(v, x.T[:,:20])
title('original')
subplot(1, 2, 2)
basis = U[:,:bases]/(math.sqrt(dv))
x_vector = numpy.dot(x, basis)*dv
plot(v, numpy.dot(basis, x_vector.T)[:,:20])
title('bases=%d'%bases)
show()
N = 200
bases = 20
S = 200
encoders = numpy.array([make_gaussian(centre=random.uniform(-1, 1),
sigma=0.1) for i in range(N)])
U, Sing, V = numpy.linalg.svd(encoders.T)
basis = U[:,:bases]/(math.sqrt(dv))
e = numpy.dot(encoders, basis)*dv
x_func = numpy.array([syde556.generate_signal(2.0, 2.0/1000,
rms=0.5, limit=2, seed=i)[0] for i in range(S)])
x = numpy.dot(x_func, basis).T*dv
ensemble = syde556.Ensemble(neurons=N, dimensions=bases, encoders=e)
decoder = ensemble.compute_decoder(x, x, noise=0.1)
A, xhat = ensemble.simulate_rate(x, decoder)
figure(figsize=(10,4))
subplot(1, 2, 1)
plot(v, x_func.T[:,:20])
ylim(-4, 4)
title('original')
subplot(1, 2, 2)
basis = U[:,:bases]/(math.sqrt(dv))
xhat_func = numpy.dot(basis, xhat)
plot(v, xhat_func[:,:20])
title('bases=%d'%bases)
ylim(-4, 4)
show()
N = 200
bases = 20
S = 200
encoders = numpy.array([make_gaussian(centre=random.uniform(-1, 1),
sigma=0.1) for i in range(N)])
U, Sing, V = numpy.linalg.svd(encoders.T)
basis = U[:,:bases]/(math.sqrt(dv))
e = numpy.dot(encoders, basis)*dv
x_func = numpy.array([syde556.generate_signal(2.0, 2.0/1000,
rms=0.5, limit=2, seed=i)[0] for i in range(S)])
x = numpy.dot(x_func, basis).T*dv
ensemble = syde556.Ensemble(neurons=N, dimensions=bases, encoders=e)
def my_function(x):
value = numpy.dot(x, basis.T)
value = -value
#value = numpy.hstack([value[100:], value[:100]])
return numpy.dot(value, basis)*dv
fx = numpy.array([my_function(xx) for xx in x.T]).T
decoder = ensemble.compute_decoder(x, fx, noise=0.1)
A, fxhat = ensemble.simulate_rate(x, decoder)
figure(figsize=(10,4))
subplot(1, 2, 1)
plot(v, x_func.T[:,:20])
ylim(-4, 4)
title('original')
subplot(1, 2, 2)
basis = U[:,:bases]/(math.sqrt(dv))
fxhat_func = numpy.dot(basis, fxhat)
plot(v, fxhat_func[:,:20])
title('bases=%d'%bases)
ylim(-4, 4)
show()