numpy.linalg.solve

Here are the examples of the python api numpy.linalg.solve taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

93 Examples 7

Page 1 Selected Page 2

Example 1

Project: lfd Source File: solver.py
Function: solve
    def solve(self, wt_n, y_nd, bend_coef, f_res):
        if y_nd.shape[0] != self.n or y_nd.shape[1] != self.d:
            raise RuntimeError("The dimensions of y_nd doesn't match the dimensions of x_nd")
        WQN = wt_n[:, None] * self.QN
        lhs = self.QN.T.dot(WQN) + bend_coef * self.NKN + self.NRN
        rhs = self.NR + WQN.T.dot(y_nd)
        z = np.linalg.solve(lhs, rhs)
        theta = self.N.dot(z)
        f_res.update(self.x_nd, y_nd, bend_coef, self.rot_coef, wt_n, theta, N=self.N, z=z)

Example 2

Project: pydy Source File: ode_function_generators.py
    @linear_sys_solver.setter
    def linear_sys_solver(self, v):

        if isinstance(v, type(lambda x: x)):
            self._solve_linear_system = v
        elif v == 'numpy':
            self._solve_linear_system = numpy.linalg.solve
        elif v == 'scipy':
            self._solve_linear_system = scipy.linalg.solve
        else:
            msg = '{} is not a valid solver.'
            raise ValueError(msg.format(self.linear_sys_solver))

Example 3

Project: hedge Source File: local.py
        @memoize_method
        def stiffness_t_matrices(self):
            return [numpy.asarray(
                la.solve(
                    self.ldis.vandermonde().T,
                    numpy.dot(
                        diff_vdm.T,
                        numpy.diag(self.volume_weights))),
                    order="C")
                    for diff_vdm in self.diff_vandermonde_matrices()]

Example 4

Project: AWS-Lambda-ML-Microservice-Skeleton Source File: test_linalg.py
    def test_0_size_k(self):
        # test zero multiple equation (K=0) case.
        class ArraySubclass(np.ndarray):
            pass
        a = np.arange(4).reshape(1, 2, 2)
        b = np.arange(6).reshape(3, 2, 1).view(ArraySubclass)

        expected = linalg.solve(a, b)[:, :, 0:0]
        result = linalg.solve(a, b[:, :, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        # test both zero.
        expected = linalg.solve(a, b)[:, 0:0, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

Example 5

Project: FreeCAD_drawing_dimensioning Source File: svgLib_dd.py
    def t_of_position( self, pos ):
        pos_element = numpy.linalg.solve( self.T_element, pos - self.c_element )
        A = numpy.linalg.solve(self.T_ellipse, pos_element)
        B = self.center_circular
        C = B + array([cos(self.theta_start), sin(self.theta_start)])
        AB = (A - B) / norm(A-B)
        BC = C - B
        theta = arccos2(dot(AB,BC))
        D = array([ A[0] - B[0], A[1] - B[1], 0.0 ])
        E = array([ A[0] - C[0], A[1] - C[1], 0.0 ])
        F = cross(D,E)
        if F[2] < 0:
            theta = -theta
        #print(theta, self.dtheta, theta / self.dtheta )
        return theta / self.dtheta

Example 6

Project: dynts Source File: ols.py
Function: beta
    def beta(self):
        '''\
The linear estimation of the parameter vector :math:`\beta` given by

.. math::

    \beta = (X^T X)^-1 X^T y
        
'''
        t = self.X.transpose()
        XX = dot(t,self.X)
        XY = dot(t,self.y)
        return linalg.solve(XX,XY)

Example 7

Project: pymatgen Source File: analyzer.py
    def get_facet_chempots(self, facet):
        """
        Calculates the chemical potentials for each element within a facet.

        Args:
            facet: Facet of the phase diagram.

        Returns:
            { element: chempot } for all elements in the phase diagram.
        """
        complist = [self._pd.qhull_entries[i].composition for i in facet]
        energylist = [self._pd.qhull_entries[i].energy_per_atom for i in facet]
        m = self._make_comp_matrix(complist)
        chempots = np.linalg.solve(m, energylist)
        return dict(zip(self._pd.elements, chempots))

Example 8

Project: hedge Source File: local.py
        @memoize_method
        def face_mass_matrix(self):
            return numpy.asarray(
                    la.solve(
                        self.ldis.face_vandermonde().T,
                        numpy.dot(
                            self.face_vandermonde().T,
                            numpy.diag(self.face_weights))),
                    order="C")

Example 9

Project: pyNastran Source File: beam.py
Function: do_problem
def do_problem(elements):
    Kg = build_global_stiffness()     # global  stiffness (Kg)
    Kr = apply_boundary_conditions()  # reduced stiffness (Kr)

    #F = K*x
    F = get_forces()
    x = solve(Kr, F)  # deflections

    (stress, strain) = recover_stress_Strain(elements, x)

Example 10

Project: PyOP2 Source File: test_matrices.py
Function: test_solve
    def test_solve(self, mat, b, x, f):
        """Solve a linear system where the solution is equal to the right-hand
        side and check the result."""
        mat.assemble()
        x = np.linalg.solve(mat.values, b.data)
        eps = 1.e-8
        assert_allclose(x, f.data, eps)

Example 11

Project: refinery Source File: WishartDistr.py
Function: e_dist_mahalanobis
  def E_dist_mahalanobis(self, X ):
    '''Calculate Mahalanobis distance to x
             dist(x) = dX'*E[Lambda]*dX
       If X is a matrix, we compute a vector of distances to each row vector
             dist(X)[n] = dX[n]'*E[Lambda]*dX[n]
    '''
    xWx = np.linalg.solve(self.cholinvW(), X.T)
    xWx *= xWx
    return self.v * np.sum(xWx, axis=0)

Example 12

Project: GPflow Source File: test_variational.py
def referenceMultivariatePriorKL(meanA, covA, meanB, covB):
    # KL[ qA | qB ] = E_{qA} \log [qA / qB] where qA and aB are K dimensional multivariate normal distributions.
    # Analytically tractable and equal to...
    # 0.5 * (Tr(covB^{-1} covA) + (meanB - meanA)^T covB^{-1} (meanB - meanA) - K + log(det(covB)) - log (det(covA)))
    K = covA.shape[0]
    traceTerm = 0.5 * np.trace(np.linalg.solve(covB, covA))
    delta = meanB - meanA
    mahalanobisTerm = 0.5 * np.dot(delta.T, np.linalg.solve(covB, delta))
    constantTerm = -0.5 * K
    priorLogDeterminantTerm = 0.5*np.linalg.slogdet(covB)[1]
    variationalLogDeterminantTerm = -0.5 * np.linalg.slogdet(covA)[1]
    return traceTerm + mahalanobisTerm + constantTerm + priorLogDeterminantTerm + variationalLogDeterminantTerm

Example 13

Project: pymanopt Source File: test_manifold_positive_definite.py
Function: test_exp
    def test_exp(self):
        man = self.man
        x = man.rand()
        u = man.randvec(x)
        e = sp.linalg.expm(la.solve(x, u))

        np_testing.assert_allclose(multiprod(x, e), man.exp(x, u))
        u = u * 1e-6
        np_testing.assert_allclose(man.exp(x, u), x + u)

Example 14

Project: ignition Source File: test_iterative_methods.py
def cg_alg_driver(alg, n=100):
    A, b = get_fd_poisson(n)
    iters, x = alg(A, b, n)
    x_numpy = solve(A, b)
    diff = x - x_numpy
    l_inf = max(diff)
    l_2 = sqrt(dot(diff, diff))
#    assert(l_inf < 1e-10)
#    assert(l_2 < 1e-10)
    return iters, l_2, l_inf

Example 15

Project: pymanopt Source File: grassmann.py
    def log(self, X, Y):
        ytx = multiprod(multitransp(Y), X)
        At = multitransp(Y) - multiprod(ytx, multitransp(X))
        Bt = np.linalg.solve(ytx, At)
        u, s, vt = svd(multitransp(Bt), full_matrices=False)
        arctan_s = np.expand_dims(np.arctan(s), -2)

        U = multiprod(u * arctan_s, vt)
        return U

Example 16

Project: hedge Source File: gmsh.py
    def __init__(self, nodes, ldis):
        self.nodes = nodes
        self.ldis  = ldis

        node_src_indices = np.array(
                ldis.get_lexicographic_gmsh_node_indices(),
                dtype=np.intp)

        nodes = np.array(nodes, dtype=np.float64)
        reordered_nodes = nodes[node_src_indices, :]

        self.modal_coeff = la.solve(
                ldis.equidistant_vandermonde(), reordered_nodes)
        # axis 0: node number, axis 1: xyz axis

        if False:
            for i, c in zip(ldis.generate_mode_identifiers(), self.modal_coeff):
                print i, c

Example 17

Project: pymatgen Source File: analyzer.py
    def get_decomposition(self, comp):
        """
        Provides the decomposition at a particular composition.

        Args:
            comp: A composition

        Returns:
            Decomposition as a dict of {Entry: amount}
        """
        facet = self._get_facet(comp)
        comp_list = [self._pd.qhull_entries[i].composition for i in facet]
        m = self._make_comp_matrix(comp_list)
        compm = self._make_comp_matrix([comp])
        decomp_amts = np.linalg.solve(m.T, compm.T)
        return {self._pd.qhull_entries[f]: amt[0]
                for f, amt in zip(facet, decomp_amts)
                if abs(amt[0]) > PDAnalyzer.numerical_tol}

Example 18

Project: pyhsmm Source File: general.py
def solve_psd(A,b,chol=None,overwrite_b=False,overwrite_A=False):
    if A.shape[0] < 5000 and chol is None:
        return np.linalg.solve(A,b)
    else:
        if chol is None:
            chol = np.linalg.cholesky(A)
        return scipy.linalg.solve_triangular(
                chol.T,
                scipy.linalg.solve_triangular(chol,b,lower=True,overwrite_b=overwrite_b),
                lower=False,overwrite_b=True)

Example 19

Project: hedge Source File: local.py
Function: mass_matrix
        @memoize_method
        def mass_matrix(self):
            return numpy.asarray(
                    la.solve(
                        self.ldis.vandermonde().T,
                        numpy.dot(
                            self.vandermonde().T,
                            numpy.diag(self.volume_weights))),
                    order="C")

Example 20

Project: elephant Source File: icsd.py
    def get_csd(self, ):
        '''
        Perform the CSD estimate from the LFP and forward matrix F, i.e as
        CSD=F**-1*LFP

        Arguments
        ---------

        Returns
        -------
        csd : np.ndarray * quantity.Quantity
            Array with the csd estimate
        '''
        csd = np.linalg.solve(self.f_matrix, self.lfp)

        return csd * (self.f_matrix.units**-1*self.lfp.units).simplified

Example 21

Project: AWS-Lambda-ML-Microservice-Skeleton Source File: math.py
Function: solve
def solve(a, b):
    """Returns the solution of A X = B."""
    try:
        return linalg.solve(a, b)
    except linalg.LinAlgError:
        return np.dot(linalg.pinv(a), b)

Example 22

Project: refinery Source File: GaussWishDistr.py
Function: dist_mahalanobis
  def dist_mahalanobis(self, X):
    '''Calculate Mahalanobis distance to x
             dist(x) = (x-m)'*W*(x-m)
       If X is a matrix, we compute a vector of distances to each row vector
             dist(X)[n] = (x_n-m)'*W*(x_n-m)
    '''           
    dX = (X-self.m).T
    Q  = np.linalg.solve( self.cholinvW(), dX )
    Q *= Q
    return np.sum( Q, axis=0)

Example 23

Project: PRML Source File: curve_fitting.py
Function: estimate
def estimate(xlist, tlist):
    """訓練データからパラメータを推定"""
    # M次多項式のときはM+1個のパラメータがある
    A = []
    for i in range(M + 1):
        for j in range(M + 1):
            temp = (xlist ** (i + j)).sum()
            A.append(temp)
    A = array(A).reshape(M + 1, M + 1)

    T = []
    for i in range(M + 1):
        T.append(((xlist ** i) * tlist).sum())
    T = array(T)

    # パラメータwは線形連立方程式の解
    wlist = np.linalg.solve(A, T)

    return wlist

Example 24

Project: refinery Source File: ZMGaussDistr.py
Function: dist_mahalanobis
  def dist_mahalanobis(self, X):
    '''  Given NxD matrix X, compute  Nx1 vector Dist
            Dist[n] = ( X[n]-m )' L (X[n]-m)
    '''
    if self.doSigma:
      Q = np.linalg.solve(self.cholSigma(), X.T)
    else:
      Q = dotABT(self.cholL(), X)
    Q *= Q
    return np.sum(Q, axis=0)

Example 25

Project: python-control Source File: statesp.py
Function: horner
    def horner(self, s):
        '''Evaluate the systems's transfer function for a complex variable

        Returns a matrix of values evaluated at complex variable s.
        '''
        resp = self.C * solve(s * eye(self.states) - self.A,
                              self.B) + self.D
        return array(resp)

Example 26

Project: dolo Source File: decision_rules.py
Function: risky_ss
    @memoized
    def risky_ss(self):
        oo = self.order
        if oo == 1:
            return self['ys']
        elif oo in (2,3):
            #TODO: RSS for order 3 should be computed differently
            import numpy.linalg      
            A = self['g_a']
            I = np.eye(A.shape[0])
            D = self['g_ss']/2
            dx = numpy.linalg.solve( I - A, D)
            return self['ys'] + dx

Example 27

Project: PRML Source File: maximum_likelihood_fitting.py
Function: estimate
def estimate(xlist, tlist):
    A = []
    for i in range(M + 1):
        for j in range(M + 1):
            temp = (xlist ** (i + j)).sum()
            A.append(temp)
    A = array(A).reshape(M + 1, M + 1)

    T = []
    for i in range(M + 1):
        T.append(((xlist ** i) * tlist).sum())
    T = array(T)

    # パラメータwは線形連立方程式の解
    wlist = np.linalg.solve(A, T)

    return wlist

Example 28

Project: pybasicbayes Source File: regression.py
    def resample(self,data=[],stats=None):
        if len(data) > 0 or stats is not None:
            stats = self._get_statistics(data) if stats is None else stats
            for itr in range(self.niter):
                self.A, self.sigma = \
                    sample_mniw(*self._natural_to_standard(
                        self.natural_hypparam + stats))

                mat = self.M_0 - self.A
                self.resample_K(1./2*np.einsum(
                    'ij,ij->j',mat,np.linalg.solve(self.sigma,mat)))
        else:
            self.resample_K()
            super(ARDRegression,self).resample()

Example 29

Project: FreeCAD_assembly2 Source File: degreesOfFreedom.py
    def determine_R_about_axis(self, R_effective, checkAnswer=True, tol=10**-12): #not used anymore
        'determine R_about_axis so that R_effective = R_about_axis * R_to_align_axis'
        A = self.R_to_align_axis.transpose()
        X = numpy.array([
                numpy.linalg.solve(A, R_effective[row,:]) for row in range(3)
                ])
        #prettyPrintArray(X)
        if checkAnswer:
            print('  determine_R_about_axis: diff between R_effective and R_about_axis * R_to_align_axis (should be all close to zero):')
            error = R_effective - dotProduct(X, self.R_to_align_axis)
            assert norm(error) <= tol
        return X

Example 30

Project: pymanopt Source File: test_manifold_positive_definite.py
Function: test_exp
    def test_exp(self):
        # Test against manopt implementation, test that for small vectors
        # exp(x, u) = x + u.
        man = self.man
        x = man.rand()
        u = man.randvec(x)
        e = np.zeros((self.k, self.n, self.n))
        for i in range(self.k):
            e[i] = sp.linalg.expm(la.solve(x[i], u[i]))
        np_testing.assert_allclose(multiprod(x, e), man.exp(x, u))
        u = u * 1e-6
        np_testing.assert_allclose(man.exp(x, u), x + u)

Example 31

Project: pgmult Source File: distributions.py
    def conditional_psi(self, x):
        """
        Compute the conditional distribution over psi given observation x and omega
        :param x:
        :return:
        """
        assert x.ndim == 2
        Omega = np.diag(self.omega)
        Sigma_cond = inv(Omega + inv(self.Sigma))

        # kappa is the mean dot precision, i.e. the sufficient statistic of a Gaussian
        # therefore we can sum over datapoints
        kappa = kappa_vec(x).sum(0)
        mu_cond = Sigma_cond.dot(kappa +
                                 solve(self.Sigma, self.mu))

        return mu_cond, Sigma_cond

Example 32

Project: pymanopt Source File: test_manifold_positive_definite.py
Function: test_dist
    def test_dist(self):
        k = self.k
        # n = self.n
        man = self.man
        x = man.rand()
        y = man.rand()

        # Compare to the implementation in manopt
        d = la.solve(x, y)
        for i in range(k):
            d[i] = sp.linalg.logm(d[i])

        np_testing.assert_almost_equal(man.dist(x, y), la.norm(d))

Example 33

Project: menpofit Source File: algorithm.py
Function: pre_compute
    def _precompute(self):
        # call super method
        super(InverseCompositional, self)._precompute()
        # compute appearance model mean gradient
        nabla_t = self.interface.gradient(self.template)
        # compute masked inverse Jacobian
        self.J_m = self.interface.steepest_descent_images(-nabla_t, self.dW_dp)
        # compute masked inverse Hessian
        self.JJ_m = self.J_m.T.dot(self.J_m)
        # compute masked Jacobian pseudo-inverse
        self.pinv_J_m = np.linalg.solve(self.JJ_m, self.J_m.T)

Example 34

Project: pyBAST Source File: classes.py
    def llik(self,P=np.array( [0., 0., 0., 0., 0., 1., 1.] ) ):
        """ Compute log-likelihood of parameter set P within 
        likelihood distribution bgmap object.
        """

        delta = self.mu - P
        sigma = self.sigma.copy()

        # Interpret infs in covariance matrix as contributing
        #  zero to the chi^2 (set delta[i] = 0).
        I = np.nonzero(np.isinf(np.diag(sigma)))[0]
        delta[I] = 0
        sigma[I,:] = 0
        sigma[:,I] = 0
        sigma[I,I] = 1

        return -0.5 * delta.dot( solve( sigma, delta ) )

Example 35

Project: pymanopt Source File: psd.py
Function: exp
    def exp(self, x, u):
        # TODO: Check which method is faster depending on n, k.
        x_inv_u = np.linalg.solve(x, u)
        if self._k > 1:
            e = np.zeros(np.shape(x))
            for i in range(self._k):
                e[i] = sp.linalg.expm(x_inv_u[i])
        else:
            e = sp.linalg.expm(x_inv_u)
        return multiprod(x, e)

Example 36

Project: dolo Source File: decision_rules.py
Function: risky_ss
    @memoized
    def risky_ss(self):
        oo = self.order
        if oo == 1:
            return self['ys']
        elif oo in (2,3):
            #TODO: RSS for order 3 should be computed differently
            import numpy.linalg
            A = self['g_a']
            I = np.eye(A.shape[0])
            D = self['g_ss']/2
            dx = numpy.linalg.solve( I - A, D)
            return self['ys'] + dx

Example 37

Project: tfdeploy Source File: tfdeploy.py
@Operation.factory(types=("MatrixSolve", "BatchMatrixSolve"))
def MatrixSolve(a, b):
    """
    Matrix solve op.
    """
    return np.linalg.solve(a, b),

Example 38

Project: python-control Source File: statesp.py
Function: dcgain
    def dcgain(self):
        """Return the zero-frequency gain

        The zero-frequency gain of a state-space system is given by:

        .. math: G(0) = - C A^{-1} B + D

        Returns
        -------
        gain : ndarray
            The zero-frequency gain, or np.nan if the system has a pole
            at the origin
        """
        try:
            gain = np.asarray(self.D -
                              self.C.dot(np.linalg.solve(self.A, self.B)))
        except LinAlgError:
            # zero eigenvalue: singular matrix
            return np.nan
        return np.squeeze(gain)

Example 39

Project: python-control Source File: statesp.py
def _rss_generate(states, inputs, outputs, type):
    """Generate a random state space.

    This does the actual random state space generation expected from rss and
    drss.  type is 'c' for continuous systems and 'd' for discrete systems.

    """

    # Probability of repeating a previous root.
    pRepeat = 0.05
    # Probability of choosing a real root.  Note that when choosing a complex
    # root, the conjugate gets chosen as well.  So the expected proportion of
    # real roots is pReal / (pReal + 2 * (1 - pReal)).
    pReal = 0.6
    # Probability that an element in B or C will not be masked out.
    pBCmask = 0.8
    # Probability that an element in D will not be masked out.
    pDmask = 0.3
    # Probability that D = 0.
    pDzero = 0.5

    # Check for valid input arguments.
    if states < 1 or states % 1:
        raise ValueError("states must be a positive integer.  states = %g." %
            states)
    if inputs < 1 or inputs % 1:
        raise ValueError("inputs must be a positive integer.  inputs = %g." %
            inputs)
    if outputs < 1 or outputs % 1:
        raise ValueError("outputs must be a positive integer.  outputs = %g." %
            outputs)

    # Make some poles for A.  Preallocate a complex array.
    poles = zeros(states) + zeros(states) * 0.j
    i = 0

    while i < states:
        if rand() < pRepeat and i != 0 and i != states - 1:
            # Small chance of copying poles, if we're not at the first or last
            # element.
            if poles[i-1].imag == 0:
                # Copy previous real pole.
                poles[i] = poles[i-1]
                i += 1
            else:
                # Copy previous complex conjugate pair of poles.
                poles[i:i+2] = poles[i-2:i]
                i += 2
        elif rand() < pReal or i == states - 1:
            # No-oscillation pole.
            if type == 'c':
                poles[i] = -exp(randn()) + 0.j
            elif type == 'd':
                poles[i] = 2. * rand() - 1.
            i += 1
        else:
            # Complex conjugate pair of oscillating poles.
            if type == 'c':
                poles[i] = complex(-exp(randn()), 3. * exp(randn()))
            elif type == 'd':
                mag = rand()
                phase = 2. * pi * rand()
                poles[i] = complex(mag * cos(phase),
                    mag * sin(phase))
            poles[i+1] = complex(poles[i].real, -poles[i].imag)
            i += 2

    # Now put the poles in A as real blocks on the diagonal.
    A = zeros((states, states))
    i = 0
    while i < states:
        if poles[i].imag == 0:
            A[i, i] = poles[i].real
            i += 1
        else:
            A[i, i] = A[i+1, i+1] = poles[i].real
            A[i, i+1] = poles[i].imag
            A[i+1, i] = -poles[i].imag
            i += 2
    # Finally, apply a transformation so that A is not block-diagonal.
    while True:
        T = randn(states, states)
        try:
            A = dot(solve(T, A), T) # A = T \ A * T
            break
        except LinAlgError:
            # In the unlikely event that T is rank-deficient, iterate again.
            pass

    # Make the remaining matrices.
    B = randn(states, inputs)
    C = randn(outputs, states)
    D = randn(outputs, inputs)

    # Make masks to zero out some of the elements.
    while True:
        Bmask = rand(states, inputs) < pBCmask
        if any(Bmask): # Retry if we get all zeros.
            break
    while True:
        Cmask = rand(outputs, states) < pBCmask
        if any(Cmask): # Retry if we get all zeros.
            break
    if rand() < pDzero:
        Dmask = zeros((outputs, inputs))
    else:
        Dmask = rand(outputs, inputs) < pDmask

    # Apply masks.
    B = B * Bmask
    C = C * Cmask
    D = D * Dmask

    return StateSpace(A, B, C, D)

Example 40

Project: QuantEcon.py Source File: matrix_eqn.py
def solve_discrete_riccati(A, B, Q, R, N=None, tolerance=1e-10, max_iter=500):
    """
    Solves the discrete-time algebraic Riccati equation

        X = A'XA - (N + B'XA)'(B'XB + R)^{-1}(N + B'XA) + Q

    via a modified structured doubling algorithm.  An explanation of the
    algorithm can be found in the reference below.

    Note that SciPy also has a discrete riccati equation solver.  However it
    cannot handle the case where R is not invertible, or when N is nonzero.
    Both of these cases can be handled in the algorithm implemented below.

    Parameters
    ----------
    A : array_like(float, ndim=2)
        k x k array.
    B : array_like(float, ndim=2)
        k x n array
    Q : array_like(float, ndim=2)
        k x k, should be symmetric and non-negative definite
    R : array_like(float, ndim=2)
        n x n, should be symmetric and positive definite
    N : array_like(float, ndim=2)
        n x k array
    tolerance : scalar(float), optional(default=1e-10)
        The tolerance level for convergence
    max_iter : scalar(int), optional(default=500)
        The maximum number of iterations allowed

    Returns
    -------
    X : array_like(float, ndim=2)
        The fixed point of the Riccati equation; a  k x k array
        representing the approximate solution

    References
    ----------
    Chiang, Chun-Yueh, Hung-Yuan Fan, and Wen-Wei Lin. "STRUCTURED DOUBLING
    ALGORITHM FOR DISCRETE-TIME ALGEBRAIC RICCATI EQUATIONS WITH SINGULAR
    CONTROL WEIGHTING MATRICES." Taiwanese Journal of Mathematics 14, no. 3A
    (2010): pp-935.

    """
    # == Set up == #
    error = tolerance + 1
    fail_msg = "Convergence failed after {} iterations."

    # == Make sure that all array_likes are np arrays, two-dimensional == #
    A, B, Q, R = np.atleast_2d(A, B, Q, R)
    n, k = R.shape[0], Q.shape[0]
    I = np.identity(k)
    if N is None:
        N = np.zeros((n, k))
    else:
        N = np.atleast_2d(N)

    # == Choose optimal value of gamma in R_hat = R + gamma B'B == #
    current_min = np.inf
    candidates = (0.0, 0.01, 0.1, 0.25, 0.5, 1.0, 2.0, 10.0, 100.0, 10e5)
    BB = dot(B.T, B)
    BTA = dot(B.T, A)
    for gamma in candidates:
        Z = R + gamma * BB
        cn = np.linalg.cond(Z)
        if np.isfinite(cn):
            Q_tilde = - Q + dot(N.T, solve(Z, N + gamma * BTA)) + gamma * I
            G0 = dot(B, solve(Z, B.T))
            A0 = dot(I - gamma * G0, A) - dot(B, solve(Z, N))
            H0 = gamma * dot(A.T, A0) - Q_tilde
            f1 = np.linalg.cond(Z, np.inf)
            f2 = gamma * f1
            f3 = np.linalg.cond(I + dot(G0, H0))
            f_gamma = max(f1, f2, f3)
            if f_gamma < current_min:
                best_gamma = gamma
                current_min = f_gamma

    # == If no candidate successful then fail == #
    if current_min == np.inf:
        msg = "Unable to initialize routine due to ill conditioned arguments"
        raise ValueError(msg)

    gamma = best_gamma
    R_hat = R + gamma * BB

    # == Initial conditions == #
    Q_tilde = - Q + dot(N.T, solve(R_hat, N + gamma * BTA)) + gamma * I
    G0 = dot(B, solve(R_hat, B.T))
    A0 = dot(I - gamma * G0, A) - dot(B, solve(R_hat, N))
    H0 = gamma * dot(A.T, A0) - Q_tilde
    i = 1

    # == Main loop == #
    while error > tolerance:

        if i > max_iter:
            raise ValueError(fail_msg.format(i))

        else:
            A1 = dot(A0, solve(I + dot(G0, H0), A0))
            G1 = G0 + dot(dot(A0, G0), solve(I + dot(H0, G0), A0.T))
            H1 = H0 + dot(A0.T, solve(I + dot(H0, G0), dot(H0, A0)))

            error = np.max(np.abs(H1 - H0))
            A0 = A1
            G0 = G1
            H0 = H1
            i += 1

    return H1 + gamma * I  # Return X

Example 41

Project: ray Source File: linalg.py
@ray.remote
def solve(a, b):
  return np.linalg.solve(a, b)

Example 42

Project: pymanopt Source File: psd.py
    def inner(self, x, u, v):
        return np.tensordot(la.solve(x, u), la.solve(x, v), axes=x.ndim)

Example 43

Project: pymc3 Source File: quadpotential.py
Function: energy
    def energy(self, x):
        L1x = solve(self.L, x)
        return .5 * dot(L1x.T, L1x)

Example 44

Project: structural_engineering Source File: system.py
Function: solve
    def solve(self):
        assert(self.system_force_vector is not None), "There are no forces on the structure"
        self.__assemble_system_matrix()
        self.__process_conditions()

        # solution of the reduced system (reduced due to support conditions)
        reduced_displacement_vector = np.linalg.solve(self.reduced_system_matrix, self.reduced_force_vector)

        # add the solution of the reduced system in the complete system displacement vector
        self.system_displacement_vector = np.zeros(self.shape_system_matrix)
        count = 0
        for i in self.remainder_indexes:
            self.system_displacement_vector[i] = reduced_displacement_vector[count]
            count += 1

        # determine the displacement vector of the elements
        for el in self.elements:
            index_node_1 = (el.node_1.ID - 1) * 3
            index_node_2 = (el.node_2.ID - 1) * 3

            for i in range(3):  # node 1 ux, uz, phi
                el.element_displacement_vector[i] = self.system_displacement_vector[index_node_1 + i]

            for i in range(3):  # node 2 ux, uz, phi
                el.element_displacement_vector[3 + i] = self.system_displacement_vector[index_node_2 + i]

            el.determine_force_vector()

        # determining the node results in post processing class
        self.post_processor.node_results()
        self.post_processor.reaction_forces()
        self.post_processor.element_results()

        # check the values in the displacement vector for extreme values, indicating a flawed calculation
        for value in np.nditer(self.system_displacement_vector):
            assert(value < 1e6), "The displacements of the structure exceed 1e6. Check your support conditions," \
                                 "or your elements Young's modulus"
        return self.system_displacement_vector

Example 45

Project: pymc3 Source File: quadpotential.py
Function: random
    def random(self):
        n = normal(size=self.L.shape[0])
        return solve(self.L.T, n)

Example 46

Project: pymanopt Source File: psd.py
Function: norm
    def norm(self, x, u):
        return la.norm(la.solve(x, u))

Example 47

Project: python-control Source File: frdata.py
Function: feedback
    def feedback(self, other=1, sign=-1):
        """Feedback interconnection between two FRD objects."""

        other = _convertToFRD(other, omega=self.omega)

        if (self.outputs != other.inputs or
            self.inputs != other.outputs):
            raise ValueError(
                "FRD.feedback, inputs/outputs mismatch")
        fresp = empty((self.outputs, self.inputs, len(other.omega)),
                      dtype=complex)
        # TODO: vectorize this
        # TODO: handle omega re-mapping
        for k, w in enumerate(other.omega):
            fresp[:, :, k] = self.fresp[:, :, k].view(type=matrix)* \
                linalg.solve(
                eye(self.inputs) +
                other.fresp[:, :, k].view(type=matrix) *
                self.fresp[:, :, k].view(type=matrix),
                eye(self.inputs))

        return FRD(fresp, other.omega, smooth=(self.ifunc is not None))

Example 48

Project: python-whiteboard Source File: wiimote.py
    def calibrate(self, p_screen, p_wii):
        l = []
        for i in range(0,4):
            l.append( [p_wii[i][0], p_wii[i][1], 1, 0, 0, 0, 
                (-p_screen[i][0] * p_wii[i][0]), 
                (-p_screen[i][0] * p_wii[i][1])] )
            l.append( [0, 0, 0, p_wii[i][0], p_wii[i][1], 1, 
                (-p_screen[i][1] * p_wii[i][0]), 
                (-p_screen[i][1] * p_wii[i][1])] )

        A = matrix(l)

        x = matrix( [
            [p_screen[0][0]],
            [p_screen[0][1]],
            [p_screen[1][0]],
            [p_screen[1][1]],
            [p_screen[2][0]],
            [p_screen[2][1]],
            [p_screen[3][0]],
            [p_screen[3][1]],
        ])

        self.hCoefs = linalg.solve(A, x)
        self.h11 = self.hCoefs[0]
        self.h12 = self.hCoefs[1]
        self.h13 = self.hCoefs[2]
        self.h21 = self.hCoefs[3]
        self.h22 = self.hCoefs[4]
        self.h23 = self.hCoefs[5]
        self.h31 = self.hCoefs[6]
        self.h32 = self.hCoefs[7]
        
        self.calibrationPoints = list(p_wii)
        self.screenPoints = list(p_screen)
        self.state = Wiimote.CALIBRATED
        
        area_inside = calculateArea(self.calibrationPoints)
        total_area = Wiimote.MAX_X * Wiimote.MAX_Y
        self.utilization = float(area_inside)/float(total_area)

Example 49

Project: pydy Source File: test_ode_function_generator.py
    def test_set_linear_system_solver(self):

        g = ODEFunctionGenerator(self.sys.eom_method.forcing_full,
                                 self.sys.coordinates,
                                 self.sys.speeds,
                                 self.sys.constants_symbols,
                                 mass_matrix=self.sys.eom_method.mass_matrix_full)

        assert g._solve_linear_system == np.linalg.solve

        g = ODEFunctionGenerator(self.sys.eom_method.forcing_full,
                                 self.sys.coordinates,
                                 self.sys.speeds,
                                 self.sys.constants_symbols,
                                 mass_matrix=self.sys.eom_method.mass_matrix_full,
                                 linear_sys_solver='numpy')

        assert g._solve_linear_system == np.linalg.solve

        g = ODEFunctionGenerator(self.sys.eom_method.forcing_full,
                                 self.sys.coordinates,
                                 self.sys.speeds,
                                 self.sys.constants_symbols,
                                 mass_matrix=self.sys.eom_method.mass_matrix_full,
                                 linear_sys_solver='scipy')

        assert g._solve_linear_system == sp.linalg.solve

        solver = lambda A, b: np.dot(np.inv(A), b)

        g = ODEFunctionGenerator(self.sys.eom_method.forcing_full,
                                 self.sys.coordinates,
                                 self.sys.speeds,
                                 self.sys.constants_symbols,
                                 mass_matrix=self.sys.eom_method.mass_matrix_full,
                                 linear_sys_solver=solver)

        assert g._solve_linear_system == solver

Example 50

Project: python-control Source File: statesp.py
    def feedback(self, other=1, sign=-1):
        """Feedback interconnection between two LTI systems."""

        other = _convertToStateSpace(other)

        # Check to make sure the dimensions are OK
        if ((self.inputs != other.outputs) or (self.outputs != other.inputs)):
                raise ValueError("State space systems don't have compatible \
inputs/outputs for feedback.")

        # Figure out the sampling time to use
        if (self.dt == None and other.dt != None):
            dt = other.dt       # use dt from second argument
        elif (other.dt == None and self.dt != None) or \
                timebaseEqual(self, other):
            dt = self.dt        # use dt from first argument
        else:
            raise ValueError("Systems have different sampling times")

        A1 = self.A
        B1 = self.B
        C1 = self.C
        D1 = self.D
        A2 = other.A
        B2 = other.B
        C2 = other.C
        D2 = other.D

        F = eye(self.inputs) - sign * D2 * D1
        if matrix_rank(F) != self.inputs:
            raise ValueError("I - sign * D2 * D1 is singular.")

        # Precompute F\D2 and F\C2 (E = inv(F))
        E_D2 = solve(F, D2)
        E_C2 = solve(F, C2)

        T1 = eye(self.outputs) + sign * D1 * E_D2
        T2 = eye(self.inputs) + sign * E_D2 * D1

        A = concatenate(
            (concatenate(
                (A1 + sign * B1 * E_D2 * C1, sign * B1 * E_C2), axis=1),
            concatenate(
                (B2 * T1 * C1, A2 + sign * B2 * D1 * E_C2), axis=1)),
            axis=0)
        B = concatenate((B1 * T2, B2 * D1 * T2), axis=0)
        C = concatenate((T1 * C1, sign * D1 * E_C2), axis=1)
        D = D1 * T2

        return StateSpace(A, B, C, D, dt)
See More Examples - Go to Next Page
Page 1 Selected Page 2