Python scipy.linalg 模块，solve() 实例源码

def lpfls(N,wp,ws,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
    b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
    b[0] = wp/np.pi
    q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    a = ln.solve(Q,b)
    h = list(nq)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def bpfls(N,ws1,wp1,wp2,ws2,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    q = W*np.sinc(nq) - (W*ws2/np.pi) * np.sinc(nq* (ws2/np.pi)) + (wp2/np.pi) * np.sinc(nq*(wp2/np.pi)) - (wp1/np.pi) * np.sinc(nq*(wp1/np.pi)) + (W*ws1/np.pi) * np.sinc(nq*(ws1/np.pi))
    b = (wp2/np.pi)*np.sinc((wp2/np.pi)*nb) - (wp1/np.pi)*np.sinc((wp1/np.pi)*nb)
    b[0] = wp2/np.pi - wp1/np.pi
    q[0] = W - W*ws2/np.pi + wp2/np.pi - wp1/np.pi + W*ws1/np.pi # since sin(pi*n)/pi*n = 1, not 0
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    a = ln.solve(Q,b)
    h = list(nq)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def hpfls(N,ws,wp,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    b = 1 - (wp/np.pi)* np.sinc(nb * wp/np.pi)
    b[0] = 1- wp/np.pi
    q = 1 - (wp/np.pi)* np.sinc(nq * wp/np.pi) + W * (ws/np.pi) * np.sinc(nq * ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
    q[0] = b[0] + W* ws/np.pi
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    a = ln.solve(Q,b)
    h = list(nq)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def solve_linear2(model:Model.fem_model):
    K_bar,F_bar,index=model.K_,model.F_,model.index
    Dvec=model.D
    Logger.info('Solving linear model with %d DOFs...'%model.DOF)
    n_nodes=model.node_count
    #sparse matrix solution
    delta_bar = sl.spsolve(sp.csc_matrix(K_bar),F_bar)
    #delta_bar=linalg.solve(K_bar,F_bar,sym_pos=True)
    delta = delta_bar
    #fill original displacement vector
    prev = 0
    for idx in index:
        gap=idx-prev
        if gap>0:
            delta=np.insert(delta,prev,[0]*gap)
        prev = idx + 1               
        if idx==index[-1] and idx!=n_nodes-1:
            delta = np.insert(delta,prev, [0]*(n_nodes*6-prev))
    delta += Dvec

    model.is_solved=True
    return delta

def demand_mass_balance_h(odemand, close_prob, recapture_rate):
    """Solve Demand Mass Balance equation for host-level

    Parameters
    ----------
    odemand: int
        Observerd demand
    close_prob: float
        Probability of selecting a closed element from host set H
    recapture_rate: float
        Estimated recapture rate

    Returns
    -------
    tuple
        Estimated demand, spill and recapture
    """

    A = np.array([[1, -1, 1], [-close_prob, 1, 0], [0, -recapture_rate, 1]])
    B = np.array([odemand, 0, 0])

    demand, spill, recapture = solve(A, B)

    return demand, spill, recapture

def fit(self, Y):
        """
        Generates the RBF coefficients to fit a set of given data values Y for centers self.centers
        :param Y: A set of dependent data values corresponding to self.centers
        :return: Void, sets the self.coefs values
        """
        kernel_matrix = self.EvaluateCentersKernel()
        kernel_matrix[np.isinf(kernel_matrix)] = 0 # TODO: Is there a better way to avoid the diagonal?
        monomial_basis = poly.GetMonomialBasis(self.dimension, self.poly_degree)
        poly_matrix = poly.BuildPolynomialMatrix(monomial_basis, self.centers.transpose()) # TODO: Probably remove transpose requirement
        poly_shape = np.shape(poly_matrix)
        # Get the number of columns, as we need to make an np.zeros((num_cols,num_cols))
        num_cols = poly_shape[1]
        num_rbf_coefs = len(self.centers)
        zero_mat = np.zeros((num_cols,num_cols))
        upper_matrix = np.hstack((kernel_matrix, poly_matrix))
        lower_matrix = np.hstack((poly_matrix.transpose(),zero_mat))
        rbf_matrix = np.vstack((upper_matrix,lower_matrix))
        Y = np.concatenate((Y,np.zeros((num_cols)))) # Extend with zeros for the polynomial annihilation
        self.coefs = sl.solve(rbf_matrix, Y, sym_pos=False)

def coeffs(u1):
    """ Calculate coefficients of basis polynomials and weights
    """
    wL  = solve(ML, u1[:N+1])
    wR  = solve(MR, u1[N:])
    oL  = weights(wL, ?s)
    oR  = weights(wR, ?s)
    if N==1:
        return (mult(wL,oL) + mult(wR,oR)) / (oL + oR)

    wCL = solve(MCL, u1[fhN:fhN2])
    oCL = weights(wCL, ?c)
    if nStencils==3:
        return (mult(wL,oL) + mult(wCL,oCL) + mult(wR,oR)) / (oL + oCL + oR)

    oCR = weights(wCR, ?c)
    wCR = solve(MCR, u1[chN:chN2])
    return (mult(wL,oL) + mult(wCL,oCL) + mult(wCR,oCR) + mult(wR,oR)) / (oL + oCL + oCR + oR)

def logistic_regression(x, t, w, eps=1e-2, max_iter=int(1e3)):
    N = x.shape[1]
    Phi = np.vstack([np.ones(N), phi(x)]).T

    for k in range(max_iter):
        y = expit(Phi.dot(w))
        R = np.diag(np.ones(N) * (y * (1 - y)))
        H = Phi.T.dot(R).dot(Phi)
        g = Phi.T.dot(y - t)

        w_new = w - linalg.solve(H, g)

        diff = linalg.norm(w_new - w) / linalg.norm(w)
        if (diff < eps):
            break

        w = w_new
        print('{0:5d} {1:10.6f}'.format(k, diff))

    return w

def compute_mtx_obj(GtG_lst, Tbeta_lst, Rc0, num_bands, K):
    """
    compute the matrix (M) in the objective function:
        min   c^H M c
        s.t.  c0^H c = 1

    :param GtG_lst: list of G^H * G
    :param Tbeta_lst: list of Teoplitz matrices for beta-s
    :param Rc0: right dual matrix for the annihilating filter (same of each block -> not a list)
    :return:
    """
    mtx = np.zeros((K + 1, K + 1), dtype=float)  # <= assume G, Tbeta and Rc0 are real-valued
    for loop in range(num_bands):
        Tbeta_loop = Tbeta_lst[loop]
        GtG_loop = GtG_lst[loop]
        mtx += np.dot(Tbeta_loop.T,
                      linalg.solve(np.dot(Rc0, linalg.solve(GtG_loop, Rc0.T)),
                                   Tbeta_loop)
                      )
    return mtx

def mldivide(a, b):
    dimensions = a.shape
    if dimensions[0] == dimensions[1]:
        return la.solve(a, b)
    else:
        return la.lstsq(a, b)[0]

def qr_ls(A, b):
    '''
    least square using QR (A must be full column rank)
    '''
    m = A.shape[0]
    n = A.shape[1]
    if rank(A) < n:
        raise Exception('Rank deficient')

    A = qr_householder(A)
    for j in range(n):
        v = np.hstack((1, A[j+1:, j]))
        A[j+1:, j] = 0
        b[j:] = (np.eye(m - j + 1) - 2 * np.outer(v, v) / la.norm(v, 2)).dot(b[j:])

    x_ls = la.solve(A[:n, :n], b[:n])

    return x_ls

def ls_qr(A, b):
    '''
    least square using QR (A must be full column rank)
    '''
    m = A.shape[0]
    n = A.shape[1]
    if rank(A) < n:
        raise Exception('Rank deficient')

    A = qr.qr_householder(A)
    for j in range(n):
        v = np.hstack((1, A[j+1:, j]))
        A[j+1:, j] = 0
        b[j:] = (np.eye(m - j + 1) - 2 * np.outer(v, v) / la.norm(v, 2)).dot(b[j:])

    x_ls = la.solve(A, b)

    return x_ls0

def single_estimate(sigin, f, fs, hmax):
    """Estimate the contribution of electrical network
    signal in sigin, if ENF is f.

    return the estimated signal
    """
    X = np.empty((len(sigin), hmax))
    fact = 2.0 * np.pi * f / fs * np.arange(1, hmax + 1)[:, None]
    p = np.arange(len(sigin))[:, None]
    X = np.dot(fact, p.T)
    X = np.concatenate([X, X + np.pi / 2]).T
    X = np.cos(X)
    XX = np.dot(X.T, X)
    Xy = np.dot(X.T, sigin)
    theta = linalg.solve(XX, Xy)
    return np.dot(X, theta)

def make_var(self):

        # self.predict: distinction between prediction and estimation
        self.Dnew.make_A(s2=(self.par.sigma**2), predict=self.predict) 

        invA_H = linalg.solve( self.Dold.A , self.Dold.H )

        temp1 = self.Dnew.H - ( self.covar.T ).dot( invA_H )
        temp2 = ( self.Dold.H.T ).dot( invA_H )
        temp3 = self.Dnew.A \
          - (self.covar.T).dot( linalg.solve( self.Dold.A , self.covar ) )

        self.var = (self.par.sigma**2) \
          * ( temp3 + temp1.dot( linalg.solve( temp2 , temp1.T ) ) )


    # create vectors of lower and upper 95% confidence intervals

def mahalanobis_distance(self):
        retrain=False

        # calculate expected value Mahalanobis distance
        MDtheo = self.Dnew.outputs.size
        try:
            MDtheovar = 2*self.Dnew.outputs.size*\
                (self.Dnew.outputs.size + self.Dold.outputs.size\
                    -self.par.beta.size - 2.0)/\
                (self.Dold.outputs.size-self.par.beta.size-4.0)
            print("theoretical Mahalanobis_distance (mean, var):" \
                  "(", MDtheo, "," , MDtheovar, ")")
        except ZeroDivisionError as e:
            print("theoretical Mahalanobis_distance mean:", MDtheo, \
                  "(too few data for variance)")

        # calculate actual Mahalanobis distance from data
        MD = ( (self.Dnew.outputs-self.mean).T ).dot\
             ( linalg.solve( self.var , (self.Dnew.outputs-self.mean) ) )
        print("calculated Mahalanobis_distance:", MD)
        retrain=True
        return retrain

def _int_var_rbf(self, X, hyp, jitter=1e-8):
        """
        Posterior integral variance of the Gaussian Process quadrature.
        X - vector (1, 2*xdim**2+xdim)
        hyp - kernel hyperparameters [s2, el_1, ... el_d]
        """
        # reshape X to SP matrix
        X = np.reshape(X, (self.n, self.d))
        # set kernel hyper-parameters
        s2, el = hyp[0], hyp[1:]
        self.kern.param_array[0] = s2  # variance
        self.kern.param_array[1:] = el  # lengthscale
        K = self.kern.K(X)
        L = np.diag(el ** 2)
        # posterior variance of the integral
        ks = s2 * np.sqrt(det(L + np.eye(self.d))) * multivariate_normal(mean=np.zeros(self.d), cov=L).pdf(X)
        postvar = -ks.dot(solve(K + jitter * np.eye(self.n), ks.T))
        return postvar

def _int_var_rbf_hyp(self, hyp, X, jitter=1e-8):
        """
        Posterior integral variance as a function of hyper-parameters
        :param hyp: RBF kernel hyper-parameters [s2, el_1, ..., el_d]
        :param X: sigma-points
        :param jitter: numerical jitter (for stabilizing computations)
        :return: posterior integral variance
        """
        # reshape X to SP matrix
        X = np.reshape(X, (self.n, self.d))
        # set kernel hyper-parameters
        s2, el = 1, hyp  # sig_var hyper always set to 1
        self.kern.param_array[0] = s2  # variance
        self.kern.param_array[1:] = el  # lengthscale
        K = self.kern.K(X)
        L = np.diag(el ** 2)
        # posterior variance of the integral
        ks = s2 * np.sqrt(det(L + np.eye(self.d))) * multivariate_normal(mean=np.zeros(self.d), cov=L).pdf(X)
        postvar = s2 * np.sqrt(det(2 * inv(L) + np.eye(self.d))) ** -1 - ks.dot(
            solve(K + jitter * np.eye(self.n), ks.T))
        return postvar

def get_ht(basis,stateA,traj_edges,test_set=None,delay=1,dt_eff=1.):
    """
    Calculates the hitting time using a galerkin method.
    """
    if test_set is None:
        test_set = basis
    # Check if any of the basis functions are nonzero on target state.
    A_locs = np.where(stateA)[0]

    if np.any(basis[A_locs]):
        raise RuntimeWarning("Some of the basis vectors are nonzero in state A.")

    L = get_generator(basis,traj_edges,test_set=test_set,delay=delay,dt_eff=dt_eff)
    beta = get_beta(stateA-1.,test_set,traj_edges,delay=delay)
    coeffs = spl.solve(L,beta)
    ht = np.dot(basis,coeffs)
    return ht,coeffs

def get_stationary_distrib(basis,traj_edges,test_set=None,delay=1,fix=0):
    """

    """
    # Initialize params. 
    N, k = np.shape(basis)
    if test_set is None:
        test_set = basis
    # Calculate Generator, Stiffness matrix
    L = get_generator(basis,traj_edges,test_set=test_set,delay=delay,dt_eff=1)
    # Get stationary distribution
    nfi = range(0,fix)+range(fix+1,len(L)) #not fixed indices.
    b = -1*L[fix,nfi]
    L_submat_transpose = (L[nfi,:][:,nfi]).T
    rho_notfixed = spl.solve(L_submat_transpose,b)
    pi = np.ones(k)
    pi[nfi] = rho_notfixed
    # Convert to values in realspace.
    pi_realspace = np.dot(basis,pi)
    pi_realspace *= np.sign(np.median(pi_realspace))
    return pi_realspace

def _solve_cholesky(X, y, alpha):
    # w = inv(X^t X + alpha*Id) * X.T y
    n_samples, n_features = X.shape
    n_targets = y.shape[1]

    A = safe_sparse_dot(X.T, X, dense_output=True)
    Xy = safe_sparse_dot(X.T, y, dense_output=True)

    one_alpha = np.array_equal(alpha, len(alpha) * [alpha[0]])

    if one_alpha:
        A.flat[::n_features + 1] += alpha[0]
        return linalg.solve(A, Xy, sym_pos=True,
                            overwrite_a=True).T
    else:
        coefs = np.empty([n_targets, n_features])
        for coef, target, current_alpha in zip(coefs, Xy.T, alpha):
            A.flat[::n_features + 1] += current_alpha
            coef[:] = linalg.solve(A, target, sym_pos=True,
                                   overwrite_a=False).ravel()
            A.flat[::n_features + 1] -= current_alpha
        return coefs

def test_axis(ST, quad, axis):
    ST = ST(N, quad=quad, plan=True)
    points, weights = ST.points_and_weights(N)
    f_hat = np.random.random(N)

    B = inner_product((ST, 0), (ST, 0))

    c = np.zeros_like(f_hat)
    c = B.solve(f_hat, c)

    # Multidimensional version
    bc = [np.newaxis,]*3
    bc[axis] = slice(None)
    fk = np.broadcast_to(f_hat[bc], (N,)*3).copy()
    ST.plan((N,)*3, axis, np.float, {})
    if hasattr(ST, 'bc'):
        ST.bc.set_tensor_bcs(ST) # To set Dirichlet boundary conditions on multidimensional array
    ck = np.zeros_like(fk)
    ck = B.solve(fk, ck, axis=axis)
    cc = [0,]*3
    cc[axis] = slice(None)
    assert np.allclose(ck[cc], c)

#test_axis(cbases.ShenDirichletBasis, "GC", 1)

def theta_hat(H_hat=None, h_hat=None, lambdaRegularizer=0.0, kernelBasis=None) :
        """
        Calculates theta_hat given H_hat, h_hat, lambda, and the kernel basis function
        Treat as a system of lienar equations and find the exact, optimal
        solution
        """
        theta_hat = linalg.solve(H_hat + (lambdaRegularizer * numpy.eye(kernelBasis)), h_hat)

        return theta_hat

def fit(self, X, y, L):
        """Fit the model according to the given training data.

        Prameters
        ---------
        X : array-like, shpae = [n_samples, n_features]
            Training data.

        y : array-like, shpae = [n_samples]
            Target values (unlabeled points are marked as 0).

        L : array-like, shpae = [n_samples, n_samples]
            Graph Laplacian.
        """

        labeled               = y != 0
        X_labeled             = X[labeled]
        y_labeled             = y[labeled]
        n_samples, n_features = X.shape
        n_labeled_samples     = y_labeled.size
        I                     = sp.eye(n_features)
        M                     = X_labeled.T @ X_labeled \
            + self.gamma_a * n_labeled_samples * I \
            + self.gamma_i * n_labeled_samples / n_samples**2 * X.T @ L**self.p @ X

        # Train a classifer
        self.coef_            = LA.solve(M, X_labeled.T @ y_labeled)

        return self

def fit(self, X, y, L):
        """Fit the model according to the given training data.

        Prameters
        ---------
        X : array-like, shpae = [n_samples, n_features]
            Training data.

        y : array-like, shpae = [n_samples]
            Target values (unlabeled points are marked as 0).

        L : array-like, shpae = [n_samples, n_samples]
            Graph Laplacian.
        """

        labeled               = y != 0
        y_labeled             = y[labeled]
        n_samples, n_features = X.shape
        n_labeled_samples     = y_labeled.size
        I                     = sp.eye(n_samples)
        J                     = sp.diags(labeled.astype(np.float64))
        K                     = rbf_kernel(X, gamma=self.gamma_k)
        M                     = J @ K \
            + self.gamma_a * n_labeled_samples * I \
            + self.gamma_i * n_labeled_samples / n_samples**2 * L**self.p @ K

        # Train a classifer
        self.dual_coef_       = LA.solve(M, y)

        return self

def _gaus_condition(self, xi):

        if np.ma.count_masked(xi) == 0:
            return xi

        a = xi.mask
        b = ~xi.mask

        xb = xi[b].data
        Laa = self.prec[np.ix_(a, a)]
        Lab = self.prec[np.ix_(a, b)]

        xfill = np.empty_like(xi)
        xfill[b] = xb
        xfill[a] = self.mean[a] - solve(Laa, Lab.dot(xb - self.mean[b]))
        return xfill

def _lobatto(alpha, beta, xl1, xl2):
    '''Compute the Lobatto nodes and weights with the preassigned node xl1, xl2.
    Based on the section 7 of the paper

        Some modified matrix eigenvalue problems,
        Gene Golub,
        SIAM Review Vol 15, No. 2, April 1973, pp.318--334,

    and

        http://www.scientificpython.net/pyblog/radau-quadrature
    '''
    from scipy.linalg import solve_banded, solve
    n = len(alpha)-1
    en = numpy.zeros(n)
    en[-1] = 1
    A1 = numpy.vstack((numpy.sqrt(beta), alpha-xl1))
    J1 = numpy.vstack((A1[:, 0:-1], A1[0, 1:]))
    A2 = numpy.vstack((numpy.sqrt(beta), alpha-xl2))
    J2 = numpy.vstack((A2[:, 0:-1], A2[0, 1:]))
    g1 = solve_banded((1, 1), J1, en)
    g2 = solve_banded((1, 1), J2, en)
    C = numpy.array(((1, -g1[-1]), (1, -g2[-1])))
    xl = numpy.array((xl1, xl2))
    ab = solve(C, xl)

    alphal = alpha
    alphal[-1] = ab[0]
    betal = beta
    betal[-1] = ab[1]
    x, w = orthopy.line.schemes.custom(alphal, betal, mode='numpy')
    return x, w

def lpfls2notch(N,wp,ws,wn1,wn2,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
    b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
    q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
    b = np.asmatrix(b)
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    G1 = np.cos(wn1*nb)
    G2 = np.cos(wn2*nb)
    G = np.matrix([G1,G2])

    d = np.array([0,0])
    d = np.asmatrix(d)
    d = d.transpose()

    c = np.asmatrix(ln.solve(Q,b))

    mu = ln.solve(G*ln.inv(Q)*G.transpose(),G*c - d)

    a = c - ln.solve(Q,G.transpose()*mu)
    h = np.zeros(N)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def lpfls1notch(N,wp,ws,wn1,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
    b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
    q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
    b = np.asmatrix(b)
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    G1 = np.cos(wn1*nb)
    G = np.matrix([G1])

    d = np.array([0])
    d = np.asmatrix(d)

    c = np.asmatrix(ln.solve(Q,b))

    mu = ln.solve(G*ln.inv(Q)*G.transpose(),G*c - d)

    a = c - ln.solve(Q,G.transpose()*mu)
    h = np.zeros(N)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def lowess(x, y, f=1. / 3., iter=5):
    """lowess(x, y, f=2./3., iter=3) -> yest
    Lowess smoother: Robust locally weighted regression.
    The lowess function fits a nonparametric regression curve to a scatterplot.
    The arrays x and y contain an equal number of elements; each pair
    (x[i], y[i]) defines a data point in the scatterplot. The function returns
    the estimated (smooth) values of y.
    The smoothing span is given by f. A larger value for f will result in a
    smoother curve. The number of robustifying iterations is given by iter. The
    function will run faster with a smaller number of iterations.
    """
    n = len(x)
    r = int(np.ceil(f * n))
    h = [np.sort(np.abs(x - x[i]))[r] for i in range(n)]
    w = np.clip(np.abs((x[:, None] - x[None, :]) / h), 0.0, 1.0)
    w = (1 - w ** 3) ** 3
    yest = np.zeros(n)
    delta = np.ones(n)
    for iteration in range(iter):
        for i in range(n):
            weights = delta * w[:, i]
            b = np.array([np.sum(weights * y), np.sum(weights * y * x)])
            A = np.array([[np.sum(weights), np.sum(weights * x)],
                          [np.sum(weights * x), np.sum(weights * x * x)]])
            beta = linalg.solve(A, b)
            yest[i] = beta[0] + beta[1] * x[i]

        residuals = y - yest
        s = np.median(np.abs(residuals))
        delta = np.clip(residuals / (6.0 * s), -1, 1)
        delta = (1 - delta ** 2) ** 2
    return yest

def control_systems(request):
    ct_sys, ref = request.param
    Ac, Bc, Cc = ct_sys.data
    Dc = np.zeros((Cc.shape[0], 1))

    Q = np.eye(Ac.shape[0])
    R = np.eye(Bc.shape[1] if len(Bc.shape) > 1 else 1)

    Sc = linalg.solve_continuous_are(Ac, Bc.reshape(-1, 1), Q, R,)
    Kc = linalg.solve(R, Bc.T @ Sc).reshape(1, -1)
    ct_ctr = LTISystem(Kc)

    evals = np.sort(np.abs(
        linalg.eig(Ac, left=False, right=False, check_finite=False)
    ))
    dT = 1/(2*evals[-1])

    Tsim = (8/np.min(evals[~np.isclose(evals, 0)])
            if np.sum(np.isclose(evals[np.nonzero(evals)], 0)) > 0
            else 8
            )

    dt_data = signal.cont2discrete((Ac, Bc.reshape(-1, 1), Cc, Dc), dT)
    Ad, Bd, Cd, Dd = dt_data[:-1]
    Sd = linalg.solve_discrete_are(Ad, Bd.reshape(-1, 1), Q, R,)
    Kd = linalg.solve(Bd.T @ Sd @ Bd + R, Bd.T @ Sd @ Ad)

    dt_sys = LTISystem(Ad, Bd, dt=dT)
    dt_sys.initial_condition = ct_sys.initial_condition
    dt_ctr = LTISystem(Kd, dt=dT)

    yield ct_sys, ct_ctr, dt_sys, dt_ctr, ref, Tsim

def test_events():
    # use bouncing ball to test events work

    # simulate in block diagram
    int_opts = block_diagram.DEFAULT_INTEGRATOR_OPTIONS.copy()
    int_opts['rtol'] = 1E-12
    int_opts['atol'] = 1E-15
    int_opts['nsteps'] = 1000
    int_opts['max_step'] = 2**-3
    x = x1, x2 = Array(dynamicsymbols('x_1:3'))
    mu, g = sp.symbols('mu g')
    constants = {mu: 0.8, g: 9.81}
    ic = np.r_[10, 15]
    sys = SwitchedSystem(
        x1, Array([0]),
        state_equations=r_[x2, -g],
        state_update_equation=r_[sp.Abs(x1), -mu*x2],
        state=x,
        constants_values=constants,
        initial_condition=ic
    )
    bd = BlockDiagram(sys)
    res = bd.simulate(5, integrator_options=int_opts)

    # compute actual impact time
    tvar = dynamicsymbols._t
    impact_eq = (x2*tvar - g*tvar**2/2 + x1).subs(
        {x1: ic[0], x2: ic[1], g: 9.81}
    )
    t_impact = sp.solve(impact_eq, tvar)[-1]

    # make sure simulation actually changes velocity sign around impact
    abs_diff_impact = np.abs(res.t - t_impact)
    impact_idx = np.where(abs_diff_impact == np.min(abs_diff_impact))[0]
    assert np.sign(res.x[impact_idx-1, 1]) != np.sign(res.x[impact_idx+1, 1])

def create_markov_operator(self, base_operator, dt, method="CrankNicolson"):
        """ Let u be the Markov operator and L be the base generator. The model dynamic
            is given by
                            u ' = L \cdot u
            Following the Method Of Line (MOL) discretization,
            let u(n) = u(t_n) where t_n is the discretized time domain.

            The following schemes can be used:
            - Explicit Euler --> u(n + 1) = (I + L * dt) \cdot u(n)
            - Implicit Euler --> (I - L * dt) \cdot u(n + 1) = u(n)
            - Crank Nicolson --> (I - L * dt / 2) \cdot u(n + 1) = (I + L * dt / 2) \cdot u(n)
        :param base_operator:
        :param dt:
        :param method:
        :return:
        """
        if method == "ExplicitEuler":
            u = np.eye(self.d) + dt * base_operator
        elif method == "ImplicitEuler":
            u = linalg.solve(np.eye(self.d) - dt * base_operator, np.eye(self.d))
        elif method == "CrankNicolson":
            u = linalg.solve(np.eye(self.d) - .5 * dt * base_operator,
                             np.eye(self.d) + .5 * dt * base_operator)
        else:
            raise NotImplementedError("Unsupported SDE resolution method")
        return u

def build_base_operator(self, t):
        """
        :param t: Not used as mu and sigma are constant
        :return:
        """
        # Update drift and volatility
        self.build_moment_vectors(t)

        base_operator = np.zeros((self.d, self.d))

        nabla = linalg.block_diag(*[self.build_gradient_matrix(x) for x in range(1, self.d - 1)])

        moments = np.zeros(2 * (self.d - 2))
        for i in range(0, self.d - 2):
            moments[2 * i] = self.drift[i + 1]
            moments[2 * i + 1] = self.volatility[i + 1]

        generator_elements = linalg.solve(nabla, moments)

        r_idx, c_idx = np.diag_indices_from(base_operator)
        base_operator[r_idx[1:-1], c_idx[:-2]] = generator_elements[::2]
        base_operator[r_idx[1:-1], c_idx[2:]] = generator_elements[1::2]
        np.fill_diagonal(base_operator, -np.sum(base_operator, axis=1))

        # -- Boundary Condition: Volatility Matching --
        nabla_0 = self.grid[1] - self.grid[0]
        base_operator[0, 0] = - 1. * self.volatility[0] / (nabla_0 * nabla_0)
        base_operator[0, 1] = - base_operator[0, 0]

        nabla_d = self.grid[self.d - 1] - self.grid[self.d - 2]
        base_operator[self.d - 1, self.d - 1] = - 1. * self.volatility[self.d - 1] / (nabla_d * nabla_d)
        base_operator[self.d - 1, self.d - 2] = - base_operator[self.d - 1, self.d - 1]
        # ----------------------------------------------

        self.sanity_check_base_operator(base_operator)

        return base_operator

def demand_mass_balance_c(host_odemand, class_odemand, avail, host_recapture):
    """Solve Demand Mass Balance equation for class-level

    Parameters
    ----------
    host_odemand: int
        Observerd host demand
    class_odemand: int
        Observed class demand
    avail: dict
        Availability of demand open during period considered
    host_recapture: float
        Estimated host level recapture

    Returns
    -------
    tuple
        Estimated demand, spill and recapture
    """

    # if observed demand of a class is 0 demand mass balance can't
    # estimate demand and spill alone without additioanl information
    demand = spill = recapture = 0
    if class_odemand:
        recapture = host_recapture * class_odemand / host_odemand

        # availability of demand closed during period considered
        k = 1 - avail
        A = np.array([[1, -1], [-k, 1]])
        B = np.array([class_odemand - recapture, 0])
        demand, spill = solve(A, B)

    return demand, spill, recapture

def generate_lin_system_from_regression_problem(n, d, sigma, filepath=None):
    (X, y, beta, e) = generate_lin_regression(n, d, sigma)
    # lambda_ = 6. * sigma**2. / n
    lambda_ = sigma**2. / (n * numpy.linalg.norm(beta) ** 2)
    A = 1. / (d * n) * X.T.dot(X) + numpy.identity(d) * lambda_
    b = 1. / (d * n) * X.T.dot(y)
    x = linalg.solve(A, b)
    cn = numpy.linalg.cond(A)
    obj = objective(X, y, x, lambda_, n)
    if filepath:
        write_system(A, b, x, filepath)
    return (A, b, X, y, lambda_, x, cn, obj)

def fit(self, X, y):
        """
        Fits a Gaussian Process regressor

        Parameters
        ----------
        X: np.ndarray, shape=(nsamples, nfeatures)
            Training instances to fit the GP.
        y: np.ndarray, shape=(nsamples,)
            Corresponding continuous target values to X.

        """
        self.X = X
        self.y = y
        self.nsamples = self.X.shape[0]
        if self.optimize:
            grads = None
            if self.usegrads:
                grads = self._grad
            self.optHyp(param_key=self.covfunc.parameters, param_bounds=self.covfunc.bounds, grads=grads)

        self.K = self.covfunc.K(self.X, self.X)
        self.L = cholesky(self.K).T
        self.alpha = solve(self.L.T, solve(self.L, y - self.mprior))
        self.logp = -.5 * np.dot(self.y, self.alpha) - np.sum(np.log(np.diag(self.L))) - self.nsamples / 2 * np.log(
            2 * np.pi)

def param_grad(self, k_param):
        """
        Returns gradient over hyperparameters. It is recommended to use `self._grad` instead.

        Parameters
        ----------
        k_param: dict
            Dictionary with keys being hyperparameters and values their queried values.

        Returns
        -------
        np.ndarray
            Gradient corresponding to each hyperparameters. Order given by `k_param.keys()`
        """
        k_param_key = list(k_param.keys())
        covfunc = self.covfunc.__class__(**k_param)
        n = self.X.shape[0]
        K = covfunc.K(self.X, self.X)
        L = cholesky(K).T
        alpha = solve(L.T, solve(L, self.y))
        inner = np.dot(np.atleast_2d(alpha).T, np.atleast_2d(alpha)) - np.linalg.inv(K)
        grads = []
        for param in k_param_key:
            gradK = covfunc.gradK(self.X, self.X, param=param)
            gradK = .5 * np.trace(np.dot(inner, gradK))
            grads.append(gradK)
        return np.array(grads)

def predict(self, Xstar, return_std=False):
        """
        Returns mean and covariances for the posterior Gaussian Process.

        Parameters
        ----------
        Xstar: np.ndarray, shape=((nsamples, nfeatures))
            Testing instances to predict.
        return_std: bool
            Whether to return the standard deviation of the posterior process. Otherwise,
            it returns the whole covariance matrix of the posterior process.

        Returns
        -------
        np.ndarray
            Mean of the posterior process for testing instances.
        np.ndarray
            Covariance of the posterior process for testing instances.
        """
        Xstar = np.atleast_2d(Xstar)
        kstar = self.covfunc.K(self.X, Xstar).T
        fmean = self.mprior + np.dot(kstar, self.alpha)
        v = solve(self.L, kstar.T)
        fcov = self.covfunc.K(Xstar, Xstar) - np.dot(v.T, v)
        if return_std:
            fcov = np.diag(fcov)
        return fmean, fcov

def fit(self, Y):
        """
        Generates the RBF coefficients to fit a set of given data values Y for centers self.centers
        :param Y: A set of dependent data values corresponding to self.centers
        :return: Void, sets the self.coefs values
        """
        kernel_matrix = self.EvaluateCentersKernel()
        self.coefs = sl.solve(kernel_matrix, Y, sym_pos=True)

def Aint(qL, qR, d):
    """ Returns the Osher-Solomon jump matrix for A, in the dth direction
        NB: eig function should be replaced with analytical function, if known
    """
    ret = zeros(n, dtype=complex128)
    ?q = qR - qL
    for i in range(N+1):
        q = qL + nodes[i] * ?q
        J = jacobian(q, d)
        ?, R = eig(J, overwrite_a=1, check_finite=0)
        b = solve(R, ?q, check_finite=0)
        ret += weights[i] * dot(R, abs(?)*b)
    return ret.real

def additive_decomp(self):
        """
        Return values for the martingale decomposition 
            - ?         : unconditional mean difference in Y
            - H         : coefficient for the (linear) martingale component (kappa_a)
            - g         : coefficient for the stationary component g(x)
            - Y_0       : it should be the function of X_0 (for now set it to 0.0)
        """
        I = np.identity(self.nx)
        A_res = la.solve(I - self.A, I)
        g = self.D @ A_res
        H = self.F + self.D @ A_res @ self.B

        return self.?, H, g

def maximization_tv(self, LU, RU):
    # ML re-estimation of the total subspace matrix or the factor loading
    # matrix
    for mix in range(self.nmix):
      idx = np.arange(self.ndim) + mix * self.ndim
      Lu = np.zeros((self.tv_dim, self.tv_dim))
      Lu[self.itril] = LU[mix, :]
      Lu += np.tril(Lu, -1).T
      self.Tm[:, idx] = solve(Lu, RU[:, idx])

def fit(self, X, T):
        Phi = self.nonlinear_transformation(X)

        # MLE for w
        self.w = linalg.solve(
            Phi.T.dot(Phi) + np.eye(self.n_basis) * self.lamb, Phi.T.dot(T))

        # MLE for beta
        n_output = 1 if T.ndim < 2 else T.shape[1]
        self.beta = n_output / np.mean(linalg.norm(T - Phi.dot(self.w))**2)

        return self

def compute_b(G_lst, GtG_lst, beta_lst, Rc0, num_bands, a_ri):
    """
    compute the uniform sinusoidal samples b from the updated annihilating
    filter coeffiients.
    :param GtG_lst: list of G^H G for different subbands
    :param beta_lst: list of beta-s for different subbands
    :param Rc0: right-dual matrix, here it is the convolution matrix associated with c
    :param num_bands: number of bands
    :param L: size of b: L by 1
    :param a_ri: a 2D numpy array. each column corresponds to the measurements within a subband
    :return:
    """
    b_lst = []
    a_Gb_lst = []
    for loop in range(num_bands):
        GtG_loop = GtG_lst[loop]
        beta_loop = beta_lst[loop]
        b_loop = beta_loop - \
                 linalg.solve(GtG_loop,
                              np.dot(Rc0.T,
                                     linalg.solve(np.dot(Rc0, linalg.solve(GtG_loop, Rc0.T)),
                                                  np.dot(Rc0, beta_loop)))
                              )

        b_lst.append(b_loop)
        a_Gb_lst.append(a_ri[:, loop] - np.dot(G_lst[loop], b_loop))

    return np.column_stack(b_lst), linalg.norm(np.concatenate(a_Gb_lst))