<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Estimating noise and stabilizing variance in calcium imaging data
The purpose of this notebook is to show how to use the VST transform for estimating the noise profile of calcium imaging data and apply a generalized Anscombe transform that aims to transform the noise into white Gaussian.
Step1:  Below is a function that will compute and apply the transformation and its inverse. The underlying noise model is scaled Poisson plus Gaussian, i.e., the underlying fluorescence value $x$ is related to the observed value $y$ by the equation
$$y = \alpha*\text{Poisson}(x) + \varepsilon$$
where $\alpha$ is non-negative scalar, and $\varepsilon \sim \mathcal{N}(\mu,\sigma^2)$ is distributed according to a Gaussian distribution.
Step2:  The transformed movie should have more uniform dynamic range (press q to exit)
Step3:  The movie might appear more noisy but information is preserved as seen from the correlation image
Step4:  The noise estimates in space should also be more uniform
Step5:  If we apply the inverse transform we approximately get back the original movie (press q to exit)
<ASSISTANT_TASK:>
Python Code:
get_ipython().magic('load_ext autoreload')
get_ipython().magic('autoreload 2')
import logging
import matplotlib.pyplot as plt
import numpy as np
import os
import timeit
logging.basicConfig(format=
                          "%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s] [%(process)d] %(message)s",
                    # filename="/tmp/caiman.log",
                    level=logging.DEBUG)
import caiman.external.houghvst.estimation as est
from caiman.external.houghvst.gat import compute_gat, compute_inverse_gat
import caiman as cm
from caiman.paths import caiman_datadir
Explanation: Estimating noise and stabilizing variance in calcium imaging data
The purpose of this notebook is to show how to use the VST transform for estimating the noise profile of calcium imaging data and apply a generalized Anscombe transform that aims to transform the noise into white Gaussian.
End of explanation
def main():
    fnames = [os.path.join(caiman_datadir(), 'example_movies', 'demoMovie.tif')]
    movie = cm.load(fnames)
    movie = movie.astype(np.float)
    # makes estimation numerically better:
    movie -= movie.mean()
    # use one every 200 frames
    temporal_stride = 100
    # use one every 8 patches (patches are 8x8 by default)
    spatial_stride = 6
    movie_train = movie[::temporal_stride]
    t = timeit.default_timer()
    estimation_res = est.estimate_vst_movie(movie_train, stride=spatial_stride)
    print('\tTime', timeit.default_timer() - t)
    alpha = estimation_res.alpha
    sigma_sq = estimation_res.sigma_sq
    movie_gat = compute_gat(movie, sigma_sq, alpha=alpha)
    # save movie_gat here
    movie_gat_inv = compute_inverse_gat(movie_gat, sigma_sq, alpha=alpha,
                                        method='asym')
    # save movie_gat_inv here
    return movie, movie_gat, movie_gat_inv
movie, movie_gat, movie_gat_inv = main()
Explanation: Below is a function that will compute and apply the transformation and its inverse. The underlying noise model is scaled Poisson plus Gaussian, i.e., the underlying fluorescence value $x$ is related to the observed value $y$ by the equation
$$y = \alpha*\text{Poisson}(x) + \varepsilon$$
where $\alpha$ is non-negative scalar, and $\varepsilon \sim \mathcal{N}(\mu,\sigma^2)$ is distributed according to a Gaussian distribution.
End of explanation
movie_gat.play(magnification=4, q_max=99.8)
Explanation: The transformed movie should have more uniform dynamic range (press q to exit):
End of explanation
CI = movie.local_correlations(swap_dim=False)
CI_gat = movie_gat.local_correlations(swap_dim=False)
plt.figure(figsize=(15,5))
plt.subplot(1,2,1); plt.imshow(CI); plt.colorbar(); plt.title('Correlation Image (original)')
plt.subplot(1,2,2); plt.imshow(CI_gat); plt.colorbar(); plt.title('Correlation Image (transformed)')
Explanation: The movie might appear more noisy but information is preserved as seen from the correlation image:
End of explanation
sn = cm.source_extraction.cnmf.pre_processing.get_noise_fft(movie.transpose(1,2,0), noise_method='mean')[0]
sn_gat = cm.source_extraction.cnmf.pre_processing.get_noise_fft(movie_gat.transpose(1,2,0), noise_method='mean')[0]
# sn = np.std(movie.transpose(1,2,0), axis=-1)
# sn_gat = np.std(movie_gat.transpose(1,2,0), axis=-1)
plt.figure(figsize=(15,5))
plt.subplot(1,2,1); plt.imshow(sn); plt.colorbar(); plt.title('Noise standard deviation (original)')
plt.subplot(1,2,2); plt.imshow(sn_gat); plt.colorbar(); plt.title('Noise standard deviation (transformed)')
Explanation: The noise estimates in space should also be more uniform
End of explanation
cm.concatenate([movie,movie_gat_inv],axis=2).play(magnification=5, q_max=99.5)
Explanation: If we apply the inverse transform we approximately get back the original movie (press q to exit):
End of explanation
<END_TASK>