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
def _get_MAP_spikes ( F , c_hat , theta , dt , tol = 1E-6 , maxiter = 100 , verbosity = 0 ) : 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
Used internally by deconvolve to compute the maximum a posteriori spike train for a given set of fluorescence traces and model parameters .