python类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