def log_likelihood(self, p):
        r"""Returns the log of the likelihood for the stored data given
        parameters ``p``.  The likelihood is computed in time
        proportional to :math:`\mathcal{O}(N)`, where :math:`N` is the
        number of data samples.

        """
        p = self.to_params(p)

        alphas, betas = self._alphas_betas(p)
        bc = self._banded_covariance(p)

        ys = self.ys.copy() - p['mu']
        ys[1:] = ys[1:] - alphas*ys[0:-1]

        dts = self.ts.reshape((-1, 1)) - self.ts.reshape((1, -1))
        tau = np.exp(p['lntau'])
        nu = self._inv_logit(p['logitnu'])
        sigma = np.exp(p['lnsigma'])
        full_cov = sigma*sigma*((1-nu)*np.exp(-np.abs(dts)/tau) + nu)

        llow = sl.cholesky_banded(bc, lower=True)

        logdet = np.sum(np.log(llow[0, :]))

        return -0.5*self.ts.shape[0]*np.log(2.0*np.pi) - logdet - 0.5*np.dot(ys, sl.cho_solve_banded((llow, True), ys))

 def __init__(self, Q, cholesky=None, banded=False):
     self.primal_shape = Q.shape[0]
     self.dual_shape = Q.shape[0]
     self.affine_offset = None
     self._Q = Q
     self.banded = banded
     if cholesky is None:
         if not self.banded:
             self._cholesky = cho_factor(Q)
         else:
             self._cholesky = cholesky_banded(Q)
     else:
         self._cholesky = cholesky

    def test_upper_real(self):
        # Symmetric positive definite banded matrix `a`
        a = array([[4.0, 1.0, 0.0, 0.0], [1.0, 4.0, 0.5, 0.0], [0.0, 0.5, 4.0, 0.2], [0.0, 0.0, 0.2, 4.0]])
        # Banded storage form of `a`.
        ab = array([[-1.0, 1.0, 0.5, 0.2], [4.0, 4.0, 4.0, 4.0]])
        c = cholesky_banded(ab, lower=False)
        ufac = zeros_like(a)
        ufac[list(range(4)), list(range(4))] = c[-1]
        ufac[(0, 1, 2), (1, 2, 3)] = c[0, 1:]
        assert_array_almost_equal(a, dot(ufac.T, ufac))

        b = array([0.0, 0.5, 4.2, 4.2])
        x = cho_solve_banded((c, False), b)
        assert_array_almost_equal(x, [0.0, 0.0, 1.0, 1.0])

    def test_lower_complex(self):
        # Hermitian positive definite banded matrix `a`
        a = array([[4.0, 1.0, 0.0, 0.0], [1.0, 4.0, 0.5, 0.0], [0.0, 0.5, 4.0, -0.2j], [0.0, 0.0, 0.2j, 4.0]])
        # Banded storage form of `a`.
        ab = array([[4.0, 4.0, 4.0, 4.0], [1.0, 0.5, 0.2j, -1.0]])
        c = cholesky_banded(ab, lower=True)
        lfac = zeros_like(a)
        lfac[list(range(4)), list(range(4))] = c[0]
        lfac[(1, 2, 3), (0, 1, 2)] = c[1, :3]
        assert_array_almost_equal(a, dot(lfac, lfac.conj().T))

        b = array([0.0, 0.5j, 3.8j, 3.8])
        x = cho_solve_banded((c, True), b)
        assert_array_almost_equal(x, [0.0, 0.0, 1.0j, 1.0])

    def test_lower_real(self):
        # Symmetric positive definite banded matrix `a`
        a = array([[4.0, 1.0, 0.0, 0.0], [1.0, 4.0, 0.5, 0.0], [0.0, 0.5, 4.0, 0.2], [0.0, 0.0, 0.2, 4.0]])
        # Banded storage form of `a`.
        ab = array([[4.0, 4.0, 4.0, 4.0], [1.0, 0.5, 0.2, -1.0]])
        c = cholesky_banded(ab, lower=True)
        lfac = zeros_like(a)
        lfac[list(range(4)), list(range(4))] = c[0]
        lfac[(1, 2, 3), (0, 1, 2)] = c[1, :3]
        assert_array_almost_equal(a, dot(lfac, lfac.T))

        b = array([0.0, 0.5, 4.2, 4.2])
        x = cho_solve_banded((c, True), b)
        assert_array_almost_equal(x, [0.0, 0.0, 1.0, 1.0])

    def test_upper_complex(self):
        # Hermitian positive definite banded matrix `a`
        a = array([[4.0, 1.0,  0.0,  0.0],
                    [1.0, 4.0,  0.5,  0.0],
                    [0.0, 0.5,  4.0, -0.2j],
                    [0.0, 0.0,  0.2j, 4.0]])
        # Banded storage form of `a`.
        ab = array([[-1.0, 1.0, 0.5, -0.2j],
                     [4.0, 4.0, 4.0,  4.0]])
        c = cholesky_banded(ab, lower=False)
        ufac = zeros_like(a)
        ufac[range(4),range(4)] = c[-1]
        ufac[(0,1,2),(1,2,3)] = c[0,1:]
        assert_array_almost_equal(a, dot(ufac.conj().T, ufac))

        b = array([0.0, 0.5, 4.0-0.2j, 0.2j + 4.0])
        x = cho_solve_banded((c, False), b)
        assert_array_almost_equal(x, [0.0, 0.0, 1.0, 1.0])

    def cholesky(self, overwrite_ab=False):
        """
        Chlesky decomposition.
        Operator needs to be positive-definite.

        Uses scipy.linalg.cholesky_banded.

        Returns a matrix in ab form
        """
        from scipy.linalg import cholesky_banded

        ab_chol = cholesky_banded(self.ab,
                               overwrite_ab=overwrite_ab,
                               lower=self.lower)
        if self.lower:
            out = LowerTriangularOperator(self.shape, ab_chol, **self.kwargs)
        else:
            out = UpperTriangularOperator(self.shape, ab_chol, **self.kwargs)
        return out

def influence_matrix_diag_chol(p, Q, R, y, w, N):
    ## Following Hutchinson & de Hoog 1985:
    #LDL factorization of Bp = 6*(1-p)*QT*D2*QT' + p*R
    D2 = sparse.diags(1./np.sqrt(w), 0, shape=(N,N))
    Bp = 6*(1-p) * Q.T*D2*Q + p*R
    Bp = Bp.todia()
    u = la.cholesky_banded(Bp.data[Bp.offsets >= 0][::-1])
    U = sparse.dia_matrix((u, Bp.offsets[Bp.offsets >= 0][::-1]), Bp.shape)
    ## may fail as the input is sparse (if so use todense), need testing
    #U = la.cholesky(A , lower=False)
    ## To get from ut u -> Ut D U:
    ## U you are looking for is the above matrix with diagonal elements
    ## replaced by 1, and D is the diagonal matrix whose diagonal elements are 
    ## the squares of the diagonal elements in the above matrix.
    d = 1./u[-1]
    U = sparse.diags(d, 0, shape=U.shape)*U
    d = d**2
    ## TODO: this should probably change...
    #5 central bands in the inverse of 6*(1-p)*QT*diag(1 ./ w)*QT' + p*R
    Binv = banded_matrix_inverse(d, U, 2)
    Hd = np.diag( sparse.eye(N,N) - (6*(1-p))*D2*Q*Binv*Q.T )
    return Hd

    def whitened_residuals(self, p):
        r"""Returns an array of the same shape as the input data, but whose
        values are indepenent :math:`N(0,1)` variables under the model
        with parameters ``p``.

        """
        p = self.to_params(p)

        alpha = self._alpha_matrix(p)
        beta = self._beta_matrix(p)

        xs = np.zeros(self.n)

        al.log_likelihood_xs_loop(self.n, self.p, alpha, self.ys - p['mu'], xs)

        beta12 = sl.cholesky_banded(beta, lower=False)
        beta12T = np.zeros((self.p, self.n))

        for i in range(self.p):
            beta12T[-(i+1), :self.n-self.p+1+i] = beta12[i, self.p-1-i:]

        return sl.solve_banded((self.p-1, 0), beta12T, xs)

 def test_cho_solve_banded(self):
     x = array([[0, -1, -1], [2, 2, 2]])
     xcho = cholesky_banded(x)
     assert_no_overwrite(lambda b: cho_solve_banded((xcho, False), b),
                         [(3,)])

#.........这里部分代码省略.........
    Here we make the knot vector (k+1)-regular by adding boundary knots:

    >>> from scipy.interpolate import make_lsq_spline, BSpline
    >>> t = [-1, 0, 1]
    >>> k = 3
    >>> t = np.r_[(x[0],)*(k+1),
    ...           t,
    ...           (x[-1],)*(k+1)]
    >>> spl = make_lsq_spline(x, y, t, k)

    For comparison, we also construct an interpolating spline for the same
    set of data:

    >>> from scipy.interpolate import make_interp_spline
    >>> spl_i = make_interp_spline(x, y)

    Plot both:

    >>> import matplotlib.pyplot as plt
    >>> xs = np.linspace(-3, 3, 100)
    >>> plt.plot(x, y, 'ro', ms=5)
    >>> plt.plot(xs, spl(xs), 'g-', lw=3, label='LSQ spline')
    >>> plt.plot(xs, spl_i(xs), 'b-', lw=3, alpha=0.7, label='interp spline')
    >>> plt.legend(loc='best')
    >>> plt.show()

    **NaN handling**: If the input arrays contain ``nan`` values, the result is
    not useful since the underlying spline fitting routines cannot deal with
    ``nan``. A workaround is to use zero weights for not-a-number data points:

    >>> y[8] = np.nan
    >>> w = np.isnan(y)
    >>> y[w] = 0.
    >>> tck = make_lsq_spline(x, y, t, w=~w)

    Notice the need to replace a ``nan`` by a numerical value (precise value
    does not matter as long as the corresponding weight is zero.)

    See Also
    --------
    BSpline : base class representing the B-spline objects
    make_interp_spline : a similar factory function for interpolating splines
    LSQUnivariateSpline : a FITPACK-based spline fitting routine
    splrep : a FITPACK-based fitting routine

    """
    x = _as_float_array(x, check_finite)
    y = _as_float_array(y, check_finite)
    t = _as_float_array(t, check_finite)
    if w is not None:
        w = _as_float_array(w, check_finite)
    else:
        w = np.ones_like(x)
    k = operator.index(k)

    if not -y.ndim <= axis < y.ndim:
        raise ValueError("axis {} is out of bounds".format(axis))
    if axis < 0:
        axis += y.ndim

    y = np.rollaxis(y, axis)    # now internally interp axis is zero

    if x.ndim != 1 or np.any(x[1:] - x[:-1] <= 0):
        raise ValueError("Expect x to be a 1-D sorted array_like.")
    if x.shape[0] < k+1:
        raise ValueError("Need more x points.")
    if k < 0:
        raise ValueError("Expect non-negative k.")
    if t.ndim != 1 or np.any(t[1:] - t[:-1] < 0):
        raise ValueError("Expect t to be a 1-D sorted array_like.")
    if x.size != y.shape[0]:
        raise ValueError('x & y are incompatible.')
    if k > 0 and np.any((x < t[k]) | (x > t[-k])):
        raise ValueError('Out of bounds w/ x = %s.' % x)
    if x.size != w.size:
        raise ValueError('Incompatible weights.')

    # number of coefficients
    n = t.size - k - 1

    # construct A.T @ A and rhs with A the collocation matrix, and
    # rhs = A.T @ y for solving the LSQ problem  ``A.T @ A @ c = A.T @ y``
    lower = True
    extradim = prod(y.shape[1:])
    ab = np.zeros((k+1, n), dtype=np.float_, order='F')
    rhs = np.zeros((n, extradim), dtype=y.dtype, order='F')
    _bspl._norm_eq_lsq(x, t, k,
                      y.reshape(-1, extradim),
                      w,
                      ab, rhs)
    rhs = rhs.reshape((n,) + y.shape[1:])

    # have observation matrix & rhs, can solve the LSQ problem
    cho_decomp = cholesky_banded(ab, overwrite_ab=True, lower=lower,
                                 check_finite=check_finite)
    c = cho_solve_banded((cho_decomp, lower), rhs, overwrite_b=True,
                         check_finite=check_finite)

    c = np.ascontiguousarray(c)
    return BSpline.construct_fast(t, c, k, axis=axis)

def qso_engine(time,data,error,ltau=3.,lvar=-1.7,sys_err=0.,return_model=False):
    """Calculates the fit quality of a damped random walk to a qso lightcurve.
    The formalism is from Rybicki & Press (1994; arXiv:comp-gas/9405004)

    Data are modelled with a covariance function
        Lij = 0.5*var*tau*exp(-|time_i-time_j|/tau) .

    Input:
        time - measurement times, typically days
        data - measured magnitudes
        error - uncertainty in measured magnitudes

    Output (dictionary):

        chi2/nu - classical variability measure
        chi2_qso/nu - for goodness of fit given fixed parameters
        chi2_qso/nu_extra - for parameter fitting, add to chi2/nu
        chi^2/nu_NULL - expected chi2/nu for non-qso variable

        signif_qso - significance chi^2/nu<chi^2/nu_NULL (rule out false alarm)
        signif_not_qso - significance chi^2/nu>1 (rule out qso)
        signif_vary - significance that source is variable
        class - resulting source type (ambiguous, not_qso, qso)

        model - time series prediction for each datum given all others (iff return_model==True)
        dmodel - model uncertainty, including uncertainty in data

    Notes:
        T = L^(-1)
        Data variance is D
        Full covariance C^(-1) = (L+D)^(-1) = T [T+D^(-1)]^(-1) D^(-1)
        Code takes advantage of the tridiagonality of T and T+D^(-1).
    """

    out_dict = {}
    out_dict['chi2_qso/nu']=999; out_dict['chi2_qso/nu_extra']=0.;
    out_dict['signif_qso']=0.; out_dict['signif_not_qso']=0.;  out_dict['signif_vary']=0.
    out_dict['chi2_qso/nu_NULL']=0.; out_dict['chi2/nu']=0.; out_dict['nu']=0
    out_dict['model']=[]; out_dict['dmodel']=[];
    out_dict['class']='ambiguous'

    lvar0 = np.log10(0.5) + lvar + ltau

    ln = len(data)
    dt = abs(time[1:]-time[:-1])

    # first make sure all dt>0
    g = np.where(dt>0.)[0]; lg = len(g)
    # must have at least 2 data points
    if lg <= 0:
        return out_dict

    if return_model:
        model = 1.*data; dmodel = -1.*error

    if lg < ln:
      dt = dt[g]
      gg = np.zeros(lg+1,dtype='int64'); gg[1:] = g+1
      dat = data[gg]; wt = 1./(sys_err**2+error[gg]**2)
      ln = lg+1
    else:
      dat = 1.*data
      wt = 1./(sys_err**2+error**2)

    out_dict['nu'] = ln-1.
    varx = np.var(dat)
    dat0 = (dat * wt).sum() / wt.sum()
    out_dict['chi2/nu'] = ((dat - dat0)**2 * wt).sum() / out_dict['nu']

    # define tridiagonal matrix T = L^(-1)
    # sparse matrix form: ab[u + i - j, j] == a[i,j]   i<=j, (here u=1)
    T = np.zeros((2,ln),dtype='float64')
    arg = dt*np.exp(-np.log(10)*ltau); ri = np.exp(-arg); ei = 1./(1./ri-ri)
    T[0,1:] = -ei; T[1,:-1] = 1.+ri*ei; T[1,1:] += ri*ei; T[1,ln-1] += 1.
    T0 = np.median(T[1,:]); T /= T0

    # equation for chi2_qso is [ (dat-x0) T Tp^(-1) D^(-1) (dat-x0) ]  , where Tp=T+D^(-1) and D^(-1)=wt
    fac = np.exp(np.log(10)*lvar0)/T0
    Tp = 1.*T
    Tp[1,:] += wt*fac
    # solve Tp*z=y for z (y=wt*dat)
    # This works for scipy __version__>='0.9.0' on anathem (20120809)
    b1 = (wt*dat).reshape((1,ln))
    b2 = b1.T
    #(Tpc,z) = solveh_banded(Tp,b2)
    z = solveh_banded(Tp,b2)
    Tpc = cholesky_banded(Tp) # the solveh_banded() function used to return the cholesky matrix, now we get seperately
    z = z.T
    z = z[0,:]
    c1 = wt.reshape((1,ln))
    c2 = c1.T
    #(Tpc,z0) = solveh_banded(Tp,c2)
    z0 = solveh_banded(Tp,c2)
    #HAS NOT CHANGED#Tpc2 = cholesky_banded(Tp)
    z0 = z0.T
    z0 = z0[0,:]


    #finally, get u=T*z
    u = T[1,:]*z; u[1:] += T[0,1:]*z[:-1]; u[:-1] += T[0,1:]*z[1:]
#.........这里部分代码省略.........