spatial.py 文件源码

def calcAvec(new, dQ, W, lambda_, active_set, M, positive):
    # TODO: comment
    r, c = np.nonzero(active_set)
#    [r,c] = find(active_set);
    Mm = -M.take(r, axis=0).take(r, axis=1)

    Mm = (Mm + Mm.T) / 2

    #% verify that there is no numerical instability
    if len(Mm) > 1:
        #        print Mm.shape
        eigMm, _ = scipy.linalg.eig(Mm)
        eigMm = np.real(eigMm)
#        check_here
    else:
        eigMm = Mm

    if any(eigMm < 0):
        np.min(eigMm)
        #%error('The matrix Mm has negative eigenvalues')
        flag = 1

    b = np.sign(W)

    if new >= 0:
        b[new] = np.sign(dQ[new])

    b = b[active_set == 1]

    if len(Mm) > 1:
        avec = np.linalg.solve(Mm, b)
    else:
        avec = b / Mm

    if positive:
        if new >= 0:
            in_ = np.sum(active_set[:new])
            if avec[in_] < 0:
                # new;
                #%error('new component of a is negative')
                flag = 1

    one_vec = np.ones(W.shape)

    dQa = np.zeros(W.shape)
    for j in range(len(r)):
        dQa = dQa + np.expand_dims(avec[j] * M[:, r[j]], axis=1)

    gamma_plus = (lambda_ - dQ) / (one_vec + dQa)
    gamma_minus = (lambda_ + dQ) / (one_vec - dQa)

    return avec, gamma_plus, gamma_minus