In [1]:
# Load data
import scipy.io as sci
import numpy as np
import wget
import ipyparallel as ipp
import time
from tqdm import tqdm
import matplotlib.pyplot as plt
In [2]:
# Load data
# Use wget or manually download to local directory
# wget.download('https://people.duke.edu/~jmp33/quantitative-neurobio/data/week3/DirectionTuning_V1_dec.mat')
vars = sci.loadmat("DirectionTuning_V1_dec.mat")
data= vars['data']
dirTuningExp = vars['dirTuningExp']
## shape of the data
Nvert, Nhorz, Nframes = data.shape
In [3]:
# first clean data
labels = -1 * np.ones([Nframes, 1])
codes = dirTuningExp['tGratingDirectionDeg'][0][0][0]
stim_names = np.unique(codes)
Ntrials = dirTuningExp['nTrials'][0][0][0][0]
Noff = dirTuningExp['stimOffFrames'][0][0][0][0]
Non = dirTuningExp['stimOnFrames'][0][0][0][0]
In [4]:
# get a data structure with one difference image per trial
def make_trial_stack(data, Noff, Non, Ntrials):
# given the data (images x times) in data and the trial labels in codes,
# return a 3d stack of images (images x trials) representing the
# *difference* in fluorescence from baseline
# Noff is the number of baseline frames per trial and Non the number of
# stimulus presentation frames
# preallocate
Nv, Nh, Nf = data.shape
imstack = np.full((Nv, Nh, Ntrials), np.nan)
# compute difference images
offset = 0
for i in range(Ntrials):
baseline = data[:, :, offset + np.array(range(Noff))].mean(axis = 2)
offset = offset + Noff
stim = data[:, :, offset + np.array(range(Non))].mean(axis = 2)
offset = offset + Non
imstack[:, :, i] = stim - baseline
return imstack
tstack = make_trial_stack(data, Noff, Non, Ntrials)
In [5]:
# new tuning function
def find_pixel_tuning(pixvals, codes, min_idx):
# given a set of pixel values and a vector of stimulus
# identities for each trial, find the preferred direction of each pixel
# based on the following:
# 1) Significant anova for pixel value vs stimulus
# 2) direction index 1 - (response_orthogonal/response_preferred) > min_idx
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
import numpy as np
alpha = 0.05
dat = {'pixels': pixvals, 'codes': np.uint16(codes).ravel()}
dat = pd.DataFrame(data=dat)
dat['codes'] = dat['codes'].astype('category')
# first ANOVA
mod = ols('pixels ~ codes', data = dat).fit()
anv = sm.stats.anova_lm(mod)
pval = anv['PR(>F)'].codes
# if passed, compute tuning index
if pval < alpha:
pixmeans = dat.groupby(by='codes', observed=False).mean()
rvals = pixmeans.values
cvals = np.uint16(pixmeans.index.values)
rpref = np.max(rvals)
pref = np.argmax(rvals)
cpref = cvals[pref]
corthogonal = (cpref + 180) % 360
rorthogonal = rvals[cvals==corthogonal]
R = 1 - rorthogonal / rpref
if R > min_idx:
preferred = pref
else:
preferred = np.nan
else:
preferred = np.nan
return preferred
In [6]:
# loop through all pixels in one order
# can also use %timeit but it performs multiple runs
# tqdm provides a progress bar that estimates the time
t0 = time.time()
preferred_img = np.full((Nvert, Nhorz), np.nan)
for hh in tqdm(range(Nhorz)):
for vv in range(Nvert):
preferred_img[vv, hh] = find_pixel_tuning(tstack[vv, hh, :], codes, 0.33)
t1 = time.time()
print(t1 - t0)
100%|█████████████████████████████████████████| 350/350 [03:31<00:00, 1.65it/s]
211.98852825164795
In [7]:
# visualize
stims = np.uint16(stim_names).astype(str)
fig, ax = plt.subplots()
cmap = plt.get_cmap('jet', len(stim_names))
cmap.set_bad('black')
cax = ax.imshow(preferred_img, cmap=cmap, aspect='auto')
cbar = fig.colorbar(cax)
cbar.ax.set_yticks(np.arange(len(stim_names)), labels=stims)
plt.show()
In [8]:
# loop through all pixels in the other order; no difference
t0 = time.time()
preferred_img = np.full((Nvert, Nhorz), np.nan)
for vv in tqdm(range(Nvert)):
for hh in range(Nhorz):
preferred_img[vv, hh] = find_pixel_tuning(tstack[vv, hh, :], codes, 0.33)
t1 = time.time()
print(t1 - t0)
100%|█████████████████████████████████████████| 125/125 [03:34<00:00, 1.72s/it]
214.91616463661194
In [9]:
# profile the codes
# -s cumulative means sorting by the cumulative time of each operation
# -l 20 means limiting to the top 20
# We see that anova is the most time consuming operation,
# except for find_pixel_tuning as a whole.
%prun -l 10 -s cumulative find_pixel_tuning(tstack[0, 0, :], codes, 0.33)
13767 function calls (13438 primitive calls) in 0.015 seconds
Ordered by: cumulative time
List reduced from 964 to 10 due to restriction <10>
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.015 0.015 {built-in method builtins.exec}
1 0.000 0.000 0.015 0.015 <string>:1(<module>)
1 0.000 0.000 0.015 0.015 3589185578.py:2(find_pixel_tuning)
1 0.000 0.000 0.007 0.007 anova.py:277(anova_lm)
1 0.000 0.000 0.007 0.007 anova.py:35(anova_single)
1 0.000 0.000 0.007 0.007 anova.py:95(anova1_lm_single)
1 0.000 0.000 0.006 0.006 model.py:147(from_formula)
1 0.000 0.000 0.005 0.005 formulatools.py:24(handle_formula_data)
1 0.000 0.000 0.005 0.005 highlevel.py:297(dmatrices)
1 0.000 0.000 0.005 0.005 highlevel.py:154(_do_highlevel_design)
In [10]:
# parallelize the codes
# another way of doing it: ipcluster start -n <desired number of engines>
n_clusters = 6
rc = ipp.Cluster(n=n_clusters).start_and_connect_sync()
# convert tstack to iterable
tstack_par = tstack.reshape((Nvert * Nhorz, Ntrials), order='F').tolist()
def func_par(x): # x is an entry of tstack_par
return main_func(x, codes, min_idx)
# run parallel
t0 = time.time()
rc[:].push(dict(main_func=find_pixel_tuning, codes=codes, min_idx=0.33))
preferred_img_par = rc[:].map_sync(func_par, tstack_par)
preferred_img_par = np.array(preferred_img_par).reshape((Nvert, Nhorz), order='F')
t1 = time.time()
print(t1 - t0)
rc.shutdown() # when you don't need it anymore, shut it down
Starting 6 engines with <class 'ipyparallel.cluster.launcher.LocalEngineSetLauncher'>
0%| | 0/6 [00:00<?, ?engine/s]
55.16378688812256
In [11]:
# visualize
fig, ax = plt.subplots()
cmap = plt.get_cmap('jet', stims.shape[0])
cmap.set_bad('black')
cax = ax.imshow(preferred_img_par, cmap=cmap, aspect='auto')
cbar = fig.colorbar(cax)
cbar.ax.set_yticklabels(stims)
plt.show()
/tmp/ipykernel_64088/1240356369.py:7: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator. cbar.ax.set_yticklabels(stims)
In [ ]: