In [1]:
%matplotlib inline
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from sys import path
path.append('..')
from oasis.functions import gen_data, gen_sinusoidal_data, deconvolve, estimate_parameters
from oasis.plotting import simpleaxis
from oasis.oasis_methods import oasisAR1, oasisAR2
In [3]:
def plot_trace(groundtruth=False):
plt.figure(figsize=(20,4))
plt.subplot(211)
plt.plot(b+c, lw=2, label='denoised')
if groundtruth:
plt.plot(true_b+true_c, c='r', label='truth', zorder=-11)
plt.plot(y, label='data', zorder=-12, c='y')
plt.legend(ncol=3, frameon=False, loc=(.02,.85))
simpleaxis(plt.gca())
plt.subplot(212)
plt.plot(s, lw=2, label='deconvolved', c='g')
if groundtruth:
for k in np.where(true_s)[0]:
plt.plot([k,k],[-.1,1], c='r', zorder=-11, clip_on=False)
plt.ylim(0,1.3)
plt.legend(ncol=3, frameon=False, loc=(.02,.85));
simpleaxis(plt.gca())
print("Correlation of deconvolved activity with ground truth ('spikes') : %.4f" % np.corrcoef(s,true_s)[0,1])
print("Correlation of denoised fluorescence with ground truth ('calcium'): %.4f" % np.corrcoef(c,true_c)[0,1])
In [4]:
# here we generate some simulated fluorescence data and plot it
true_b = 2
y, true_c, true_s = map(np.squeeze, gen_data(N=1, b=true_b, seed=0))
plt.figure(figsize=(20,4))
plt.plot(y, c='y');
In [5]:
%time c, s, b, g, lam = deconvolve(y, penalty=1)
In [6]:
plot_trace(True)
In [7]:
%time c, s, b, g, lam = deconvolve(y, penalty=0)
plot_trace(True)
$g$ is related to the decay time $\tau_d$ of the exponetial Ca response kernel $e^{-t/\tau_d}$ as $g=e^{-\frac{1}{\tau_d r}}$ with decay time $\tau_d$ in seconds and imaging rate r in Hz.
In [8]:
# Here we provide the ground truth values for b and g.
# The Ca response kernel to 1 AP has maximal amplitude 1 and we pick s_min slightly larger than 1/2.
%time c, s = oasisAR1(y-true_b, g=.95, s_min=.55)
plot_trace(True)
In function deconvolve, we already multiply the autocorrelation estimate by some fudge_factor that is close to but smaller than 1, which increases robustness
In [9]:
# no optimization of g
y, true_c, true_s = map(np.squeeze, gen_sinusoidal_data(N=1, b=true_b, seed=0))
%time c, s, b, g, lam = deconvolve(y, penalty=0)
# the line below uses the direct autocorrelation estimate without fudge_factor, yielding worse results
# %time c, s, b, g, lam = deconvolve(y, g=estimate_parameters(y, 1)[0], penalty=0)
plot_trace(True)
In [10]:
# optimization of g. optimize_g=5 uses 5 large isolated calcium events to update g.
%time c, s, b, g, lam = deconvolve(y, penalty=0, optimize_g=5)
plot_trace(True)
In [ ]:
In [ ]:
In [11]:
# here we generate some simulated fluorescence data and plot it
y, true_c, true_s = map(np.squeeze, gen_data([1.7,-.712], sn=.5, N=1, b=true_b, seed=0))
plt.figure(figsize=(20,4))
plt.plot(y, c='y');
In [12]:
# g=(None,None) because we want to use an AR(2) model, but don't know the parameters
%time c, s, b, g, lam = deconvolve(y, g=(None,None), penalty=1)
plot_trace(True)
In [13]:
%time c, s, b, g, lam = deconvolve(y, g=(None,None), penalty=0)
plot_trace(True)
In [14]:
# this is currently slow, we hope to still improve speed as well as robustness
%time c, s, b, g, lam = deconvolve(y, g=(None,None), penalty=0, optimize_g=5, max_iter=5)
plot_trace(True)
$g=(g_1,g_2)$ is related to the decay time $\tau_d$ and rise time $\tau_r$ (in seconds) of the Ca response kernel $e^{-t/\tau_d}-e^{-t/\tau_r}$ as
$g_1=e^{-\frac{1}{\tau_d r}}+e^{-\frac{1}{\tau_r r}}$ and
$g_2=-e^{-\frac{1}{\tau_d r}}\cdot e^{-\frac{1}{\tau_r r}}$ with imaging rate r in Hz.
In [15]:
# Here we provide the ground truth values for g.
%time c, s, b, g, lam = deconvolve(y, (1.7,-.712))
plot_trace(True)
In [16]:
# Here we provide the ground truth values for b and g and we pick s_min slightly larger than 1/2.
%time c, s = oasisAR2(y-true_b, 1.7,-.712, s_min=.55)
plot_trace(True)
In [ ]: