In [1]:
import sys
sys.path.append('..')

In [2]:
import numpy as np
import electrode2currentmap as e2cm
import effectivecurrent2brightness as ec2b
from scipy import interpolate
from utils import TimeSeries
import matplotlib.pyplot as plt
%matplotlib inline

In [3]:
r = e2cm.Retina(axon_map='../retina_1700by2800.npz',
                sampling=25, ylo=-1700, yhi=1700, xlo=-2800, xhi=2800)

In [4]:
xlist=[]
ylist=[]
rlist=[]
e_spacing=525

for x in np.arange(-2362, 2364, e_spacing):
    for y in np.arange(-1312, 1314, e_spacing):
        xlist.append(x)
        ylist.append(y)
        rlist.append(100)

e_all = e2cm.ElectrodeArray(rlist,xlist,ylist)

e_rf=[]
for e in e_all.electrodes:
    e_rf.append(e2cm.receptive_field(e, r.gridx, r.gridy,e_spacing))

In [5]:
# create movie
# original screen was [52.74, 63.32]  visual angle
# res=[768 ,1024] # resolution of screen
#pixperdeg=degscreen/res
# no need to simulate the whole movie, just match it to the retina
# xhi+xlo/294 (microns per degree)

degscreen=[11.6, 19.1] # array visual angle,
res=[e_rf[0].shape[0],e_rf[1].shape[1]] # resolution of screen

fps=30
bar_width=6.7
[X,Y]=np.meshgrid(np.linspace(-degscreen[1]/2, degscreen[1]/2, res[1]),
np.linspace(-degscreen[0]/2, degscreen[0]/2, res[0]));

In [6]:
orientations = np.arange(0, 2*np.pi, 2*np.pi/4)

In [7]:
#for o in orientations: # each orientation
o = orientations[0]
M=np.cos(o)*X +np.sin(o)*Y
#   for sp in range (32:32): # DEBUG each speed, eventually 8:32
sp=8

In [8]:
movie=np.zeros((res[0],res[1], int(np.ceil((70/5)*30))))
st=np.min(M)
fm_ct=1
while (st<np.max(M)):
    img=np.zeros(M.shape)
    ind=np.where((M>st) & (M<st+bar_width))
    img[ind]=1
    movie[:,:, fm_ct]=img
    fm_ct=fm_ct+1
    st=st+(sp/fps)

movie=movie[:,:, 0:fm_ct-1]
moviedur=movie.shape[2]/fps

In [9]:
pt=[]
for rf in e_rf:
    rflum = e2cm.retinalmovie2electrodtimeseries(rf, movie)
    #plt.plot(rflum)
    ptrain = e2cm.Movie2Pulsetrain(rflum)
    #plt.plot(ptrain.data)
    pt.append(ptrain)

[ecs_list, cs_list]  = r.electrode_ecs(e_all)

In [10]:
tm1 = ec2b.TemporalModel()
#fr=np.zeros([e_rf[0].shape[0],e_rf[0].shape[1], len(pt[0].data)])

In [11]:
def brightness(r, xx, yy, ecs_list, pt, tm):
    ecm = r.ecm(xx, yy, ecs_list, pt)
    fr = tm.fast_response(ecm, dojit=True)    
    ca = tm.charge_accumulation(fr, ecm)
    sn = tm.stationary_nonlinearity(ca)
    sr = tm.slow_response(sn)
    sr.resample(25)
    return sr.data

In [12]:
sr_0_0 = brightness(r, 0, 0, ecs_list, pt, tm1)

In [13]:
sr_0_0.shape


Out[13]:
(23776,)

In [52]:
brightness_movie = np.zeros((4, 4, sr_0_0.shape[0]))

In [54]:
from itertools import product
idx = list(product(*(range(s) for s in brightness_movie.shape[:-1])))

In [55]:
def calc_pixel(arr, idx, r, ecs_list, pt, tm):
    ecm = r.ecm(*idx, ecs_list, pt)
    fr = tm.fast_response(ecm, dojit=True)    
    ca = tm.charge_accumulation(fr, ecm)
    sn = tm.stationary_nonlinearity(ca)
    sr = tm.slow_response(sn)
    sr.resample(25)
    return sr.data

In [56]:
%%timeit
ff = np.array([calc_pixel(brightness_movie, i, r, ecs_list, pt, tm1) for i in idx]).reshape(brightness_movie.shape)


1 loops, best of 3: 5.41 s per loop

In [57]:
import imp

In [58]:
import utils
imp.reload(utils)


Out[58]:
<module 'utils' from '../utils.py'>

In [59]:
%%timeit
answer = utils.parfor(brightness_movie, calc_pixel, r, ecs_list, pt, tm1, n_jobs=10, axis=-1)


1 loops, best of 3: 3.39 s per loop

In [ ]: