pyfnnd/_fnndeconv.py
def _get_MAP_spikes(F, c_hat, theta, dt, tol=1E-6, maxiter=100, verbosity=0):
    """
    Used internally by deconvolve to compute the maximum a posteriori
    spike train for a given set of fluorescence traces and model parameters.

    See the documentation for deconvolve for the meaning of the
    arguments

    Returns:    n_hat_best, c_hat_best, LL_best

    """
    npix, nt = F.shape

    sigma, alpha, beta, lamb, gamma = theta

    # we project everything onto the alpha mask so that we only ever have to
    # deal with 1D vector norms
    alpha_ss = np.dot(alpha, alpha)
    c = np.dot(alpha, F) - np.dot(alpha, beta)

    # used for computing the LL and gradient
    scale_var = 1. / (2 * sigma * sigma)
    lD = lamb * dt

    # used for computing the gradient (M.T.dot(lamb * dt))
    grad_lnprior = np.zeros(nt, dtype=DTYPE)
    grad_lnprior[1:] = lD
    grad_lnprior[:-1] -= lD * gamma

    # initialize the weight of the barrier term to 1
    z = 1.

    # initial estimate of spike probabilities
    n_hat = c_hat[1:] - gamma * c_hat[:-1]
    # assert not np.any(n_hat < 0), "spike probabilities < 0"

    # (actual - predicted) fluorescence
    res = c - alpha_ss * c_hat

    # best overall posterior log-likelihood of the fluorescence
    LL_best = _post_LL(n_hat, res, scale_var, lD, z)
    LL_barrier = LL_best

    nloop1 = 0
    terminate_interior = False

    # in the outer loop we'll progressively reduce the weight of the barrier
    # term and check the interior point termination criteria
    while not terminate_interior:

        nloop2 = 0
        terminate_barrier = False

        # converge for this barrier weight
        while not terminate_barrier:

            # by projecting everything onto alpha, we reduce this to a 1D
            # vector norm
            res = c - alpha_ss * c_hat

            # compute direction of newton step
            d = _direction(n_hat, res, alpha_ss, gamma, scale_var,
                           grad_lnprior, z)

            terminate_linesearch = False

            # find the largest step we can take in direction d without
            # violating the non-negativity constraint on n_hat
            s_upper_bnd = -n_hat / (d[1:] - gamma * d[:-1])

            # we are only interested in positive step sizes
            feasible = (s_upper_bnd > 0)

            if np.any(feasible):
                # largest allowable step size is 1.
                s = min(1., 0.999 * np.min(s_upper_bnd[feasible]))
            else:
                # if there is no step size that will keep n_hat >= 0, just
                # reduce the barrier weight and try again
                terminate_linesearch = True
                terminate_barrier = True

                if verbosity >= 2:
                    print("skipping linesearch: no positive step size will "
                          "keep n_hat >= 0")

            nloop3 = 0

            # backtracking line search for the largest step size that increases
            # the posterior log-likelihood of the fluorescence
            while not terminate_linesearch:

                # update estimated calcium
                c_hat_line = c_hat + (s * d)

                # update spike probabilities
                n_hat_line = c_hat_line[1:] - gamma * c_hat_line[:-1]
                # assert not np.any(n_hat_line < 0), "spike probabilities < 0"

                # (actual - predicted) fluorescence
                res = c - alpha_ss * c_hat_line

                # compute the new posterior log-likelihood
                LL_line = _post_LL(n_hat_line, res, scale_var, lD, z)
                # assert not np.any(np.isnan(LL1)), "nan LL"

                if verbosity >= 2:
                    print('spikes: iter=%3i, %3i, %3i; z=%-10.4g; s=%-10.4g;'
                          ' LL=%-10.4g'
                          % (nloop1, nloop2, nloop3, z, s, LL_line))

                # if the step size gets too small without making any progress,
                # we terminate the linesearch and reduce the barrier weight
                if s < S_TOL:
                    if verbosity >= 2:
                        print('--> terminated linesearch: s < %.3g on %i '
                              'iterations' % (S_TOL, nloop3))
                    terminate_linesearch = True
                    terminate_barrier = True

                # only update c_hat & LL if LL improved
                if LL_line > LL_barrier:
                    LL_barrier, n_hat, c_hat = LL_line, n_hat_line, c_hat_line
                    terminate_linesearch = True

                # reduce the step size
                else:
                    s /= S_FAC

                nloop3 += 1

            # if d gets too small, reduce the barrier weight
            if (np.linalg.norm(d) < D_TOL):
                terminate_barrier = True

            nloop2 += 1

        # only test for convergence if we were actually able to enter the
        # linesearch
        if nloop3:
            delta_LL = -(LL_barrier - LL_best) / LL_best
            LL_best = LL_barrier

            if (delta_LL < tol):
                terminate_interior = True

        elif z < Z_TOL:
            if verbosity >= 2:
                print('MAP spike train failed to converge before z < %.3g'
                      % Z_TOL)
            terminate_interior = True

        elif nloop1 > maxiter:
            if verbosity >= 2:
                print('MAP spike train failed to converge within maxiter (%i)'
                       % maxiter)
            terminate_interior = True

        # increment the outer loop counter, reduce the barrier weight
        nloop1 += 1
        z /= Z_FAC

    return n_hat, c_hat, LL_best
def _get_MAP_spikes(F, c_hat, theta, dt, tol=1E-6, maxiter=100, verbosity=0):
    """
    Used internally by deconvolve to compute the maximum a posteriori
    spike train for a given set of fluorescence traces and model parameters.

    See the documentation for deconvolve for the meaning of the
    arguments

    Returns:    n_hat_best, c_hat_best, LL_best

    """
    npix, nt = F.shape

    sigma, alpha, beta, lamb, gamma = theta

    # we project everything onto the alpha mask so that we only ever have to
    # deal with 1D vector norms
    alpha_ss = np.dot(alpha, alpha)
    c = np.dot(alpha, F) - np.dot(alpha, beta)

    # used for computing the LL and gradient
    scale_var = 1. / (2 * sigma * sigma)
    lD = lamb * dt

    # used for computing the gradient (M.T.dot(lamb * dt))
    grad_lnprior = np.zeros(nt, dtype=DTYPE)
    grad_lnprior[1:] = lD
    grad_lnprior[:-1] -= lD * gamma

    # initialize the weight of the barrier term to 1
    z = 1.

    # initial estimate of spike probabilities
    n_hat = c_hat[1:] - gamma * c_hat[:-1]
    # assert not np.any(n_hat < 0), "spike probabilities < 0"

    # (actual - predicted) fluorescence
    res = c - alpha_ss * c_hat

    # best overall posterior log-likelihood of the fluorescence
    LL_best = _post_LL(n_hat, res, scale_var, lD, z)
    LL_barrier = LL_best

    nloop1 = 0
    terminate_interior = False

    # in the outer loop we'll progressively reduce the weight of the barrier
    # term and check the interior point termination criteria
    while not terminate_interior:

        nloop2 = 0
        terminate_barrier = False

        # converge for this barrier weight
        while not terminate_barrier:

            # by projecting everything onto alpha, we reduce this to a 1D
            # vector norm
            res = c - alpha_ss * c_hat

            # compute direction of newton step
            d = _direction(n_hat, res, alpha_ss, gamma, scale_var,
                           grad_lnprior, z)

            terminate_linesearch = False

            # find the largest step we can take in direction d without
            # violating the non-negativity constraint on n_hat
            s_upper_bnd = -n_hat / (d[1:] - gamma * d[:-1])

            # we are only interested in positive step sizes
            feasible = (s_upper_bnd > 0)

            if np.any(feasible):
                # largest allowable step size is 1.
                s = min(1., 0.999 * np.min(s_upper_bnd[feasible]))
            else:
                # if there is no step size that will keep n_hat >= 0, just
                # reduce the barrier weight and try again
                terminate_linesearch = True
                terminate_barrier = True

                if verbosity >= 2:
                    print("skipping linesearch: no positive step size will "
                          "keep n_hat >= 0")

            nloop3 = 0

            # backtracking line search for the largest step size that increases
            # the posterior log-likelihood of the fluorescence
            while not terminate_linesearch:

                # update estimated calcium
                c_hat_line = c_hat + (s * d)

                # update spike probabilities
                n_hat_line = c_hat_line[1:] - gamma * c_hat_line[:-1]
                # assert not np.any(n_hat_line < 0), "spike probabilities < 0"

                # (actual - predicted) fluorescence
                res = c - alpha_ss * c_hat_line

                # compute the new posterior log-likelihood
                LL_line = _post_LL(n_hat_line, res, scale_var, lD, z)
                # assert not np.any(np.isnan(LL1)), "nan LL"

                if verbosity >= 2:
                    print('spikes: iter=%3i, %3i, %3i; z=%-10.4g; s=%-10.4g;'
                          ' LL=%-10.4g'
                          % (nloop1, nloop2, nloop3, z, s, LL_line))

                # if the step size gets too small without making any progress,
                # we terminate the linesearch and reduce the barrier weight
                if s < S_TOL:
                    if verbosity >= 2:
                        print('--> terminated linesearch: s < %.3g on %i '
                              'iterations' % (S_TOL, nloop3))
                    terminate_linesearch = True
                    terminate_barrier = True

                # only update c_hat & LL if LL improved
                if LL_line > LL_barrier:
                    LL_barrier, n_hat, c_hat = LL_line, n_hat_line, c_hat_line
                    terminate_linesearch = True

                # reduce the step size
                else:
                    s /= S_FAC

                nloop3 += 1

            # if d gets too small, reduce the barrier weight
            if (np.linalg.norm(d) < D_TOL):
                terminate_barrier = True

            nloop2 += 1

        # only test for convergence if we were actually able to enter the
        # linesearch
        if nloop3:
            delta_LL = -(LL_barrier - LL_best) / LL_best
            LL_best = LL_barrier

            if (delta_LL < tol):
                terminate_interior = True

        elif z < Z_TOL:
            if verbosity >= 2:
                print('MAP spike train failed to converge before z < %.3g'
                      % Z_TOL)
            terminate_interior = True

        elif nloop1 > maxiter:
            if verbosity >= 2:
                print('MAP spike train failed to converge within maxiter (%i)'
                       % maxiter)
            terminate_interior = True

        # increment the outer loop counter, reduce the barrier weight
        nloop1 += 1
        z /= Z_FAC

    return n_hat, c_hat, LL_best
Used internally by deconvolve to compute the maximum a posteriori
    spike train for a given set of fluorescence traces and model parameters.

    See the documentation for deconvolve for the meaning of the
    arguments

    Returns:    n_hat_best, c_hat_best, LL_best