def deconvolve(F, c0=None, theta0=((None,) * 5), dt=0.02, rate=0.5, tau=1.,
               learn_theta=((0,) * 5), norm_alpha=True, params_tol=1E-3,
               spikes_tol=1E-3, spikes_maxiter=100, params_maxiter=20,
               decimate=0, verbosity=0):
    """

    Fast Non-Negative Deconvolution
    ---------------------------------------------------------------------------
    This function uses an interior point method to solve the following
    optimization problem:

        n_hat = argmax_{n >= 0} P(n | F)

    where n_hat_best is a maximum a posteriori estimate for the most likely
    spike train, given the fluorescence signal F, and the model:

    c_{t} = gamma * c_{t-1} + n_{t},            n_{t} ~ Poisson(lambda * dt)
    F_{t} = alpha * c_{t} + beta + epsilon,     epsilon ~ N(0, sigma)

    It is also possible to estimate the model parameters sigma, alpha, beta and
    lambda from the data using pseudo-EM updates.

    Arguments:
    ---------------------------------------------------------------------------
    F: ndarray, [nt] or [npix, nt]
        measured fluorescence values

    c0: ndarray, [nt]
        initial estimate of the calcium concentration for each time bin

    theta0: len(5) sequence
        initial estimates of the model parameters (sigma, alpha, beta, lambda,
        gamma). any parameters in theta0 that are None will be estimated from
        F.

    dt: float scalar
        duration of each time bin (s)

    rate: float scalar
        estimate of mean firing rate (Hz), ignored if theta0[3] is not None

    tau: float scalar
        estimate of calcium decay time constant (s), ignored if theta0[4] is
        not None

    learn_theta: len(5) bool sequence
        specifies which of the model parameters to attempt learn via pseudo-EM
        iterations. currently gamma cannot be optimised.

    norm_alpha: bool
        if True (default), we impose that sum(alpha) == 1 by scaling n_hat and
        c_hat

    spikes_tol: float scalar
        termination condition for interior point spike train estimation:
            params_tol > abs((LL_prev - LL) / LL)

    params_tol: float scalar
        as above, but for the model parameter estimation

    spikes_maxiter: int scalar
        maximum number of interior point iterations to estimate MAP spike train

    params_maxiter: int scalar
        maximum number of pseudo-EM iterations to estimate model parameters

    decimate: int scalar
        option to decimate the input array over time when initalizing and
        updating the theta parameters, i.e. Fd = F[:, ::decimate]. this can
        be useful for very large input arrays. the default value (0) means no
        decimation is performed.

    verbosity: int scalar
        0: no convergence messages (default)
        1: convergence messages for model parameters
        2: convergence messages for model parameters & MAP spike train

    Returns:
    ---------------------------------------------------------------------------
    n_hat_best: ndarray, [nt]
        MAP estimate of the most likely spike train

    c_hat_best: ndarray, [nt]
        estimated intracellular calcium concentration (A.U.)

    LL_best: float scalar
        posterior log-likelihood of F given n_hat_best and theta_best

    theta_best: len(5) tuple
        model parameters, updated according to learn_theta

    Reference:
    ---------------------------------------------------------------------------
    Vogelstein, J. T., Packer, A. M., Machado, T. A., Sippy, T., Babadi, B.,
    Yuste, R., & Paninski, L. (2010). Fast nonnegative deconvolution for spike
    train inference from population calcium imaging. Journal of
    Neurophysiology, 104(6), 3691-704. doi:10.1152/jn.01073.2009

    """

    tstart = time.time()

    F = np.atleast_2d(F.astype(DTYPE))

    npix, nt = F.shape

    # ensure that F is non-negative
    offset = F.min() - EPS
    F = F - offset

    theta = _init_theta(F, theta0, offset, dt=dt, rate=rate, tau=tau,
                        decimate=decimate)

    sigma, alpha, beta, lamb, gamma = theta

    if c0 is None:

        # let n0 be a uniform vector of estimated mean spike probability, then
        # push this through the forward model to get c0. this way c0 is
        # guaranteed not to result in negative spike probabities in n_hat on
        # the first iteration.
        n0 = lamb * dt * np.ones(nt, dtype=DTYPE)
        c0 = signal.lfilter(np.r_[1.], np.r_[1., -gamma], n0, axis=0)

    else:
        # force a copy because c0 will be modified in place.
        c0 = c0.copy()

    # if we're not learning the parameters, this step is all we need to do
    n_hat, c_hat, LL = _get_MAP_spikes(F, c0, theta, dt, spikes_tol,
                                       spikes_maxiter, verbosity)

    # pseudo-EM iterations to optimize the model parameters
    if any(learn_theta):

        if verbosity >= 1:
            print('params: iter=%3i; LL=%-10.4g; delta_LL= N/A' % (0, LL))

        nloop1 = 1
        done = False

        while not done:

            # a 'full' parameter update, as used in the Vogelstein paper/code
            theta1 = _update_theta(n_hat, c_hat, F, theta, dt, learn_theta,
                                     decimate)

            # get the new n_hat, c_hat, and LL
            n1, c_hat1, LL1 = _get_MAP_spikes(F, c_hat, theta1, dt, spikes_tol,
                                              spikes_maxiter, verbosity)

            # test for convergence
            delta_LL = -((LL1 - LL) / LL)

            if verbosity >= 1:
                print('params: iter=%3i; LL=%-10.4g; delta_LL=%-10.4g'
                      % (nloop1, LL1, delta_LL))

            if delta_LL > 0:

                # keep the new parameters
                n_hat, c_hat, LL, theta = n1, c_hat1, LL1, theta1

                # if the LL is not improving significantly, time to terminate
                if delta_LL < params_tol:
                    if verbosity >= 1:
                        print("Parameter optimization converged after %i "
                              "iterations" % nloop1)
                    done = True

                elif nloop1 == params_maxiter:
                    if verbosity >= 1:
                        print('Parameter optimization failed to converge '
                              'before maxiter was reached (%i)' % nloop1)
                    done = True

            else:
                if verbosity >= 1:
                    print('Terminating parameter optimization on %i '
                          'iterations: LL is decreasing' % nloop1)
                done = True

            if done:
                if verbosity >= 1:
                    print("Last delta log-likelihood:\t%-10.4g" % delta_LL)
                    print("Best posterior log-likelihood:\t%10.3f" % LL)

            # increment the loop counter
            nloop1 += 1

    if verbosity >= 1:
        time_taken = time.time() - tstart
        print("Completed: %s" % s2h(time_taken))

    sigma, alpha, beta, lamb, gamma = theta

    # correct for the offset we originally applied to F
    beta += offset

    # we can impose that sum(alpha) == 1 by scaling c_hat and n_hat to
    # compensate
    if norm_alpha:

        alpha_sum = np.sum(alpha)
        if alpha_sum < 0:
            warnings.warn('sum(alpha) is negative - n_hat may therefore '
                          'contain negative values')
        alpha /= alpha_sum
        c_hat *= alpha_sum

        n_hat = c_hat[1:] - gamma * c_hat[:-1]

    # since we can't use FNND to estimate the spike probabilities in the 0th
    # timebin, for convenience we just concatenate 0 to the start of
    # n_hat so that it has the same shape as F and c_hat
    n_hat = np.r_[0, n_hat]

    theta = sigma, alpha, beta, lamb, gamma

    return n_hat, c_hat, LL, theta
def deconvolve ( F , c0 = None , theta0 = ( ( None , ) * 5 ) , dt = 0.02 , rate = 0.5 , tau = 1. , learn_theta = ( ( 0 , ) * 5 ) , norm_alpha = True , params_tol = 1E-3 , spikes_tol = 1E-3 , spikes_maxiter = 100 , params_maxiter = 20 , decimate = 0 , verbosity = 0 ) : tstart = time . time ( ) F = np . atleast_2d ( F . astype ( DTYPE ) ) npix , nt = F . shape offset = F . min ( ) - EPS F = F - offset theta = _init_theta ( F , theta0 , offset , dt = dt , rate = rate , tau = tau , decimate = decimate ) sigma , alpha , beta , lamb , gamma = theta if c0 is None : n0 = lamb * dt * np . ones ( nt , dtype = DTYPE ) c0 = signal . lfilter ( np . r_ [ 1. ] , np . r_ [ 1. , - gamma ] , n0 , axis = 0 ) else : c0 = c0 . copy ( ) n_hat , c_hat , LL = _get_MAP_spikes ( F , c0 , theta , dt , spikes_tol , spikes_maxiter , verbosity ) if any ( learn_theta ) : if verbosity >= 1 : print ( 'params: iter=%3i; LL=%-10.4g; delta_LL= N/A' % ( 0 , LL ) ) nloop1 = 1 done = False while not done : theta1 = _update_theta ( n_hat , c_hat , F , theta , dt , learn_theta , decimate ) n1 , c_hat1 , LL1 = _get_MAP_spikes ( F , c_hat , theta1 , dt , spikes_tol , spikes_maxiter , verbosity ) delta_LL = - ( ( LL1 - LL ) / LL ) if verbosity >= 1 : print ( 'params: iter=%3i; LL=%-10.4g; delta_LL=%-10.4g' % ( nloop1 , LL1 , delta_LL ) ) if delta_LL > 0 : n_hat , c_hat , LL , theta = n1 , c_hat1 , LL1 , theta1 if delta_LL < params_tol : if verbosity >= 1 : print ( "Parameter optimization converged after %i " "iterations" % nloop1 ) done = True elif nloop1 == params_maxiter : if verbosity >= 1 : print ( 'Parameter optimization failed to converge ' 'before maxiter was reached (%i)' % nloop1 ) done = True else : if verbosity >= 1 : print ( 'Terminating parameter optimization on %i ' 'iterations: LL is decreasing' % nloop1 ) done = True if done : if verbosity >= 1 : print ( "Last delta log-likelihood:\t%-10.4g" % delta_LL ) print ( "Best posterior log-likelihood:\t%10.3f" % LL ) nloop1 += 1 if verbosity >= 1 : time_taken = time . time ( ) - tstart print ( "Completed: %s" % s2h ( time_taken ) ) sigma , alpha , beta , lamb , gamma = theta beta += offset if norm_alpha : alpha_sum = np . sum ( alpha ) if alpha_sum < 0 : warnings . warn ( 'sum(alpha) is negative - n_hat may therefore ' 'contain negative values' ) alpha /= alpha_sum c_hat *= alpha_sum n_hat = c_hat [ 1 : ] - gamma * c_hat [ : - 1 ] n_hat = np . r_ [ 0 , n_hat ] theta = sigma , alpha , beta , lamb , gamma return n_hat , c_hat , LL , theta
Fast Non-Negative Deconvolution
    ---------------------------------------------------------------------------
    This function uses an interior point method to solve the following
    optimization problem:

        n_hat = argmax_{n >= 0} P(n | F)

    where n_hat_best is a maximum a posteriori estimate for the most likely
    spike train, given the fluorescence signal F, and the model:

    c_{t} = gamma * c_{t-1} + n_{t},            n_{t} ~ Poisson(lambda * dt)
    F_{t} = alpha * c_{t} + beta + epsilon,     epsilon ~ N(0, sigma)

    It is also possible to estimate the model parameters sigma, alpha, beta and
    lambda from the data using pseudo-EM updates.

    Arguments:
    ---------------------------------------------------------------------------
    F: ndarray, [nt] or [npix, nt]
        measured fluorescence values

    c0: ndarray, [nt]
        initial estimate of the calcium concentration for each time bin

    theta0: len(5) sequence
        initial estimates of the model parameters (sigma, alpha, beta, lambda,
        gamma). any parameters in theta0 that are None will be estimated from
        F.

    dt: float scalar
        duration of each time bin (s)

    rate: float scalar
        estimate of mean firing rate (Hz), ignored if theta0[3] is not None

    tau: float scalar
        estimate of calcium decay time constant (s), ignored if theta0[4] is
        not None

    learn_theta: len(5) bool sequence
        specifies which of the model parameters to attempt learn via pseudo-EM
        iterations. currently gamma cannot be optimised.

    norm_alpha: bool
        if True (default), we impose that sum(alpha) == 1 by scaling n_hat and
        c_hat

    spikes_tol: float scalar
        termination condition for interior point spike train estimation:
            params_tol > abs((LL_prev - LL) / LL)

    params_tol: float scalar
        as above, but for the model parameter estimation

    spikes_maxiter: int scalar
        maximum number of interior point iterations to estimate MAP spike train

    params_maxiter: int scalar
        maximum number of pseudo-EM iterations to estimate model parameters

    decimate: int scalar
        option to decimate the input array over time when initalizing and
        updating the theta parameters, i.e. Fd = F[:, ::decimate]. this can
        be useful for very large input arrays. the default value (0) means no
        decimation is performed.

    verbosity: int scalar
        0: no convergence messages (default)
        1: convergence messages for model parameters
        2: convergence messages for model parameters & MAP spike train

    Returns:
    ---------------------------------------------------------------------------
    n_hat_best: ndarray, [nt]
        MAP estimate of the most likely spike train

    c_hat_best: ndarray, [nt]
        estimated intracellular calcium concentration (A.U.)

    LL_best: float scalar
        posterior log-likelihood of F given n_hat_best and theta_best

    theta_best: len(5) tuple
        model parameters, updated according to learn_theta

    Reference:
    ---------------------------------------------------------------------------
    Vogelstein, J. T., Packer, A. M., Machado, T. A., Sippy, T., Babadi, B.,
    Yuste, R., & Paninski, L. (2010). Fast nonnegative deconvolution for spike
    train inference from population calcium imaging. Journal of
    Neurophysiology, 104(6), 3691-704. doi:10.1152/jn.01073.2009