Python numpy.linalg.cholesky() Examples

The following are code examples for showing how to use numpy.linalg.cholesky(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_linalg.py    MIT License 7 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 2
Project: recruit   Author: Frank-qlu   File: test_linalg.py    Apache License 2.0 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 3
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 4
Project: MARRtino-2.0   Author: DaniAffCH   File: test_linalg.py    GNU General Public License v3.0 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 5
Project: vnpy_crypto   Author: birforce   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 6
Project: ble5-nrf52-mac   Author: tomasero   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 7
Project: P3_image_processing   Author: latedude2   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 8
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 9
Project: predictive-maintenance-using-machine-learning   Author: awslabs   File: test_linalg.py    Apache License 2.0 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 10
Project: fund   Author: Frank-qlu   File: test_linalg.py    Apache License 2.0 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 11
Project: pySINDy   Author: luckystarufo   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 12
Project: pymanopt   Author: pymanopt   File: psd.py    BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def exp(self, x, u):
        # TODO: Check which method is faster depending on n, k.
        x_inv_u = la.solve(x, u)
        if self._k > 1:
            e = np.zeros(np.shape(x))
            for i in range(self._k):
                e[i] = sl.expm(x_inv_u[i])
        else:
            e = sl.expm(x_inv_u)
        return multiprod(x, e)
        # This alternative implementation is sometimes faster though less
        # stable. It can return a matrix with small negative determinant.
        #    c = la.cholesky(x)
        #    c_inv = la.inv(c)
        #    e = multiexp(multiprod(multiprod(c_inv, u), multitransp(c_inv)),
        #                 sym=True)
        #    return multiprod(multiprod(c, e), multitransp(c)) 
Example 13
Project: Programming-for-Non-Technical-Roles-   Author: PacktPublishing   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 14
Project: nispat   Author: amarquand   File: gp.py    GNU General Public License v3.0 6 votes vote down vote up
def post(self, hyp, covfunc, X, y):
        """ Generic function to compute posterior distribution.
        """
        
        if len(hyp.shape) > 1: # force 1d hyperparameter array
            hyp = hyp.flatten()

        if len(X.shape) == 1:
            X = X[:, np.newaxis]
        self.N, self.D = X.shape

        # hyperparameters
        sn2 = np.exp(2*hyp[0])       # noise variance
        theta = hyp[1:]            # (generic) covariance hyperparameters

        if self.verbose:
            print("estimating posterior ... | hyp=", hyp)

        self.K = covfunc.cov(theta, X)
        self.L = chol(self.K + sn2*np.eye(self.N))
        self.alpha = solve(self.L.T, solve(self.L, y))
        self.hyp = hyp
        self.covfunc = covfunc 
Example 15
Project: inference-tools   Author: C-bowman   File: gp_tools.py    MIT License 6 votes vote down vote up
def loo_likelihood(self, theta):
        """
        Calculates the 'leave-one out' (LOO) log-likelihood.

        This implementation is based on equations (5.10, 5.11, 5.12) from
        Rasmussen & Williams.
        """
        try:
            K_xx = self.cov.build_covariance(theta) + self.sig
            L = cholesky(K_xx)

            # Use the Cholesky decomposition of the covariance to find its inverse
            I = eye(len(self.x))
            iK = solve_triangular(L.T, solve_triangular(L, I, lower = True))
            alpha = solve_triangular(L.T, solve_triangular(L, self.y, lower = True))
            var = 1. / diag(iK)
            return -0.5*(var*alpha**2 + log(var)).sum()
        except:
            warn('Cholesky decomposition failure in loo_likelihood')
            return -1e50 
Example 16
Project: inference-tools   Author: C-bowman   File: gp_tools.py    MIT License 6 votes vote down vote up
def loo_likelihood_gradient(self, theta):
        """
        Calculates the 'leave-one out' (LOO) log-likelihood, as well as its
        gradient with respect to the hyperparameters.

        This implementation is based on equations (5.10, 5.11, 5.12, 5.13, 5.14)
        from Rasmussen & Williams.
        """
        K_xx, grad_K = self.cov.covariance_and_gradients(theta)
        K_xx += self.sig
        L = cholesky(K_xx)

        # Use the Cholesky decomposition of the covariance to find its inverse
        I = eye(len(self.x))
        iK = solve_triangular(L.T, solve_triangular(L, I, lower = True))
        alpha = solve_triangular(L.T, solve_triangular(L, self.y, lower = True))
        var = 1. / diag(iK)
        LOO = -0.5*(var*alpha**2 + log(var)).sum()

        grad = zeros(len(grad_K))
        for i,dK in enumerate(grad_K):
            Z = iK.dot(dK)
            grad[i] = ((alpha*Z.dot(alpha) - 0.5*(1 + var*alpha**2)*diag(Z.dot(iK))) * var).sum()

        return LOO, grad*exp(theta) 
Example 17
Project: inference-tools   Author: C-bowman   File: gp_tools.py    MIT License 6 votes vote down vote up
def marginal_likelihood_gradient(self, theta):
        """
        returns the log-marginal likelihood and its gradient with respect
        to the hyperparameters.

        This implementation is based on equations (5.8, 5.9) from Rasmussen & Williams.
        """
        K_xx, grad_K = self.cov.covariance_and_gradients(theta)
        K_xx += self.sig

        # get the cholesky decomposition
        L = cholesky(K_xx)
        # calculate the log-marginal likelihood
        alpha = solve_triangular(L.T, solve_triangular(L, self.y, lower = True))
        LML = -0.5*dot( self.y.T, alpha ) - log(diagonal(L)).sum()
        # calculate the gradients
        iK = solve_triangular(L.T, solve_triangular(L, eye(L.shape[0]), lower = True))
        Q = alpha[:,None]*alpha[None,:] - iK
        grad = array([ 0.5*(Q*dK.T).sum() for dK in grad_K])*exp(theta)
        return LML, grad 
Example 18
Project: facethin   Author: ParkerGod   File: test_linalg.py    GNU General Public License v3.0 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 19
Project: islam-buddy   Author: hamir   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 20
Project: MCT   Author: robbert-harms   File: rCovSoS.py    GNU Lesser General Public License v3.0 6 votes vote down vote up
def __init__(self, channels, covariance_noise_matrix, **kwargs):
        """Instantiate the rCovSos method.

        This will do a cholesky decomposition on the input covariance noise matrix, inverts it and takes the
        transpose. We then take per voxel the dot product of the signals with this resulting matrix.

        Args:
            channels (list): the list of input nifti files, one for each channel element. Every nifti file
                    should be a 4d matrix with on the 4th dimension all the time series. The length of this list
                    should equal the number of input channels.
            covariance_noise_matrix (str or ndarray): the corresponding noise matrix to use. If a string is given it is
                supposed to be a nifti file path.
        """
        super().__init__(channels, **kwargs)
        if isinstance(covariance_noise_matrix, str):
            covariance_noise_matrix = load_nifti(covariance_noise_matrix).get_data()
        self._inverse_covar = cholesky(inv(covariance_noise_matrix)) 
Example 21
Project: Deribit_funding_rate_indicator   Author: Dimasik007   File: test_linalg.py    MIT License 6 votes vote down vote up
def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 
Example 22
Project: LaserTOF   Author: kyleuckert   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 23
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 24
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_unsupported_commontype():
    # linalg gracefully handles unsupported type
    arr = np.array([[1, -2], [2, 5]], dtype='float16')
    with assert_raises_regex(TypeError, "unsupported in linalg"):
        linalg.cholesky(arr) 
Example 25
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 26
Project: recruit   Author: Frank-qlu   File: test_linalg.py    Apache License 2.0 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 27
Project: recruit   Author: Frank-qlu   File: test_regression.py    Apache License 2.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 28
Project: Efficient_Augmentation   Author: mkuchnik   File: PSD_util.py    MIT License 5 votes vote down vote up
def isPD(B):
    """Returns true when input is positive-definite, via Cholesky"""
    try:
        _ = la.cholesky(B)
        return True
    except la.LinAlgError:
        return False 
Example 29
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_regression.py    GNU General Public License v3.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367 ],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 30
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 5 votes vote down vote up
def wishartrand(nu, phi):
    dim = phi.shape[0]
    chol = cholesky(phi)
    #nu = nu+dim - 1
    #nu = nu + 1 - np.arange(1,dim+1)
    foo = np.zeros((dim,dim))
    
    for i in range(dim):
        for j in range(i+1):
            if i == j:
                foo[i,j] = np.sqrt(chi2.rvs(nu-(i+1)+1))
            else:
                foo[i,j]  = npr.normal(0,1)
    return np.dot(chol, np.dot(foo, np.dot(foo.T, chol.T))) 
Example 31
Project: NiaPy   Author: NiaOrg   File: es.py    MIT License 5 votes vote down vote up
def CovarianceMaatrixAdaptionEvolutionStrategyF(task, epsilon=1e-20, rnd=rand):
	lam, alpha_mu, hs, sigma0 = (4 + round(3 * log(task.D))) * 10, 2, 0, 0.3 * task.bcRange()
	mu = int(round(lam / 2))
	w = log(mu + 0.5) - log(range(1, mu + 1))
	w = w / sum(w)
	mueff = 1 / sum(w ** 2)
	cs = (mueff + 2) / (task.D + mueff + 5)
	ds = 1 + cs + 2 * max(sqrt((mueff - 1) / (task.D + 1)) - 1, 0)
	ENN = sqrt(task.D) * (1 - 1 / (4 * task.D) + 1 / (21 * task.D ** 2))
	cc, c1 = (4 + mueff / task.D) / (4 + task.D + 2 * mueff / task.D), 2 / ((task.D + 1.3) ** 2 + mueff)
	cmu, hth = min(1 - c1, alpha_mu * (mueff - 2 + 1 / mueff) / ((task.D + 2) ** 2 + alpha_mu * mueff / 2)), (1.4 + 2 / (task.D + 1)) * ENN
	ps, pc, C, sigma, M = full(task.D, 0.0), full(task.D, 0.0), eye(task.D), sigma0, full(task.D, 0.0)
	x = rnd.uniform(task.bcLower(), task.bcUpper())
	x_f = task.eval(x)
	while not task.stopCondI():
		pop_step = asarray([rnd.multivariate_normal(full(task.D, 0.0), C) for _ in range(int(lam))])
		pop = asarray([task.repair(x + sigma * ps, rnd) for ps in pop_step])
		pop_f = apply_along_axis(task.eval, 1, pop)
		isort = argsort(pop_f)
		pop, pop_f, pop_step = pop[isort[:mu]], pop_f[isort[:mu]], pop_step[isort[:mu]]
		if pop_f[0] < x_f: x, x_f = pop[0], pop_f[0]
		M = sum(w * pop_step.T, axis=1)
		ps = solve(chol(C).conj() + epsilon, ((1 - cs) * ps + sqrt(cs * (2 - cs) * mueff) * M + epsilon).T)[0].T
		sigma *= exp(cs / ds * (norm(ps) / ENN - 1)) ** 0.3
		ifix = where(sigma == inf)
		if any(ifix): sigma[ifix] = sigma0
		if norm(ps) / sqrt(1 - (1 - cs) ** (2 * (task.Iters + 1))) < hth: hs = 1
		else: hs = 0
		delta = (1 - hs) * cc * (2 - cc)
		pc = (1 - cc) * pc + hs * sqrt(cc * (2 - cc) * mueff) * M
		C = (1 - c1 - cmu) * C + c1 * (tile(pc, [len(pc), 1]) * tile(pc.reshape([len(pc), 1]), [1, len(pc)]) + delta * C)
		for i in range(mu): C += cmu * w[i] * tile(pop_step[i], [len(pop_step[i]), 1]) * tile(pop_step[i].reshape([len(pop_step[i]), 1]), [1, len(pop_step[i])])
		E, V = eig(C)
		if any(E < epsilon):
			E = fmax(E, 0)
			C = lstsq(V.T, dot(V, diag(E)).T, rcond=None)[0].T
	return x, x_f 
Example 32
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 33
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 34
Project: inception_loop2019   Author: cajal   File: multi_mei.py    MIT License 5 votes vote down vote up
def gauss2d(vx, vy, mu, cov):
    input_shape = vx.shape
    mu_x, mu_y = mu
    v = np.stack([vx.ravel() - mu_x, vy.ravel() - mu_y])
    cinv = inv(cholesky(cov))
    y = cinv @ v
    g = np.exp(-0.5 * (y * y).sum(axis=0))
    return g.reshape(input_shape) 
Example 35
Project: MARRtino-2.0   Author: DaniAffCH   File: test_linalg.py    GNU General Public License v3.0 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 36
Project: MARRtino-2.0   Author: DaniAffCH   File: test_linalg.py    GNU General Public License v3.0 5 votes vote down vote up
def test_unsupported_commontype():
    # linalg gracefully handles unsupported type
    arr = np.array([[1, -2], [2, 5]], dtype='float16')
    with assert_raises_regex(TypeError, "unsupported in linalg"):
        linalg.cholesky(arr) 
Example 37
Project: MARRtino-2.0   Author: DaniAffCH   File: test_regression.py    GNU General Public License v3.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 38
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 39
Project: vnpy_crypto   Author: birforce   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 40
Project: vnpy_crypto   Author: birforce   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 41
Project: GPt   Author: ialong   File: test_mvn_logp.py    Apache License 2.0 5 votes vote down vote up
def prepare_L(self):
        L = mvn(self.D, self.D)
        L = cholesky(L @ L.T)
        return L 
Example 42
Project: GPt   Author: ialong   File: test_KL.py    Apache License 2.0 5 votes vote down vote up
def choleskify(mats):
    ret = []
    for m in mats:
        if m.ndim == 1:
            m = np.abs(m)
        else:
            m = cholesky(m @ (m.T if m.ndim == 2 else np.transpose(m, [0, 2, 1])))
        ret.append(m)
    return ret 
Example 43
Project: ble5-nrf52-mac   Author: tomasero   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 44
Project: ble5-nrf52-mac   Author: tomasero   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 45
Project: Computable   Author: ktraunmueller   File: linear_algebra.py    MIT License 5 votes vote down vote up
def cholesky_decomposition(a):
    return linalg.cholesky(a)

# Eigenvalues 
Example 46
Project: Computable   Author: ktraunmueller   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367 ],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 47
Project: poker   Author: surgebiswas   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 48
Project: P3_image_processing   Author: latedude2   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 49
Project: P3_image_processing   Author: latedude2   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_unsupported_commontype():
    # linalg gracefully handles unsupported type
    arr = np.array([[1, -2], [2, 5]], dtype='float16')
    with assert_raises_regex(TypeError, "unsupported in linalg"):
        linalg.cholesky(arr) 
Example 50
Project: P3_image_processing   Author: latedude2   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 51
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 52
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 53
Project: 3dprinteros-client   Author: panasevychol   File: test_regression.py    GNU Affero General Public License v3.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367 ],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 54
Project: 3dprinteros-client   Author: panasevychol   File: test_regression.py    GNU Affero General Public License v3.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367 ],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 55
Project: predictive-maintenance-using-machine-learning   Author: awslabs   File: test_linalg.py    Apache License 2.0 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 56
Project: predictive-maintenance-using-machine-learning   Author: awslabs   File: test_regression.py    Apache License 2.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 57
Project: fund   Author: Frank-qlu   File: test_linalg.py    Apache License 2.0 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 58
Project: fund   Author: Frank-qlu   File: test_regression.py    Apache License 2.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 59
Project: PyMDNet   Author: zhyj3038   File: bbox_regressor.py    GNU General Public License v3.0 5 votes vote down vote up
def solve(A, y, delta, method):
  if method == 'ridge_reg_chol':
    R = cholesky(dot(A.T, A) + delta*np.identity(A.shape[1]))
    z = lstsq(R.T, dot(A.T, y))[0]
    x = lstsq(R, z)[0]
  elif method == 'ridge_reg_inv':
    x = dot(dot(inv(dot(A.T, A) + delta*np.identity(A.shape[1])), A.T), y)
  elif method == 'ls_mldivide':
    if delta > 0:
      print('ignoring lambda; no regularization used')
    x = lstsq(A, y)[0]
  loss = 0.5 * (dot(A, x) - y) **2
  return x.reshape(-1, 1) 
Example 60
Project: copulae   Author: DanielBok   File: _multivariate_t.py    MIT License 5 votes vote down vote up
def cdf(cls, x: Numeric, mean: Numeric = None, cov: Numeric = 1, df: float = None):
        # TODO implement CDF
        # dim, mean, cov, df = cls._process_parameters(None, mean, cov, df)
        # x = cls._process_input(x, dim)
        #
        # C = la.cholesky(cov)
        # err, eps, = np.inf, 1e-6
        #
        # _g = 3  # standard deviations for monte carlo method
        # N = len(x)
        # int_vals, var_sums = np.zeros(N), np.zeros(N)
        #
        # for xi in range(N):
        #     int_val, var_sum, n = 0, 0, 0
        #     while err > eps:
        #         ws = np.random.uniform(size=dim)
        #
        #         ys = np.zeros(dim)
        #         n += 1
        #         for i in range(dim):
        #             dfi = df + i
        #             ys[i] = np.sqrt((df + (ys ** 2).sum()) / dfi) / t.pdf(ws[i], dfi)
        #             if (C[i, :] * ys).sum() > x[xi, i]:
        #                 break
        #         else:
        #             var_sum += (n - 1) / n * (1 - int_val) ** 2
        #             int_val += (1 - int_val) / n
        #             err = _g * np.sqrt(var_sum / (n * (n - 1)))
        #             print(err)
        #     int_vals[xi] = int_val
        #     var_sums[xi] = var_sum
        #
        # return int_vals
        raise NotImplementedError 
Example 61
Project: pySINDy   Author: luckystarufo   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 62
Project: pySINDy   Author: luckystarufo   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 63
Project: safe-exploration   Author: befelix   File: utils_visualization.py    MIT License 5 votes vote down vote up
def plot_ellipsoid_2D(p, q, ax, n_points=100, color="r"):
    """ Plot an ellipsoid in 2D

    TODO: Untested!

    Parameters
    ----------
    p: 3x1 array[float]
        Center of the ellipsoid
    q: 3x3 array[float]
        Shape matrix of the ellipsoid
    ax: matplotlib.Axes object
        Ax on which to plot the ellipsoid

    Returns
    -------
    ax: matplotlib.Axes object
        The Ax containing the ellipsoid
    """
    plt.sca(ax)
    r = nLa.cholesky(q).T;  # checks spd inside the function
    t = np.linspace(0, 2 * np.pi, n_points);
    z = [np.cos(t), np.sin(t)];
    ellipse = np.dot(r, z) + p;
    handle, = ax.plot(ellipse[0, :], ellipse[1, :], color)

    return ax, handle 
Example 64
Project: safe-exploration   Author: befelix   File: utils_visualization.py    MIT License 5 votes vote down vote up
def plot_ellipsoid_2D(p, q, ax, n_points = 100, color = "r"):
    """ Plot an ellipsoid in 2D

    TODO: Untested!

    Parameters
    ----------
    p: 3x1 array[float]
        Center of the ellipsoid
    q: 3x3 array[float]
        Shape matrix of the ellipsoid
    ax: matplotlib.Axes object
        Ax on which to plot the ellipsoid

    Returns
    -------
    ax: matplotlib.Axes object
        The Ax containing the ellipsoid
    """
    plt.sca(ax)
    r = nLa.cholesky(q).T; #checks spd inside the function
    t = np.linspace(0, 2*np.pi, n_points);
    z = [np.cos(t), np.sin(t)];
    ellipse = np.dot(r,z) + p;
    handle, = ax.plot(ellipse[0,:], ellipse[1,:],color)

    return ax, handle 
Example 65
Project: pymanopt   Author: pymanopt   File: psd.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def dist(self, x, y):
        # Adapted from equation 6.13 of "Positive definite matrices". Chol
        # decomp gives the same result as matrix sqrt. There may be a more
        # efficient way to compute this!
        c = la.cholesky(x)
        c_inv = la.inv(c)
        logm = multilog(multiprod(multiprod(c_inv, y), multitransp(c_inv)),
                        pos_def=True)
        return la.norm(logm) 
Example 66
Project: pymanopt   Author: pymanopt   File: psd.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def norm(self, x, u):
        # This implementation is as fast as np.linalg.solve_triangular and is
        # more stable, as the above solver tends to output non positive
        # definite results.
        c = la.cholesky(x)
        c_inv = la.inv(c)
        return la.norm(multiprod(multiprod(c_inv, u), multitransp(c_inv))) 
Example 67
Project: pymanopt   Author: pymanopt   File: psd.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def log(self, x, y):
        c = la.cholesky(x)
        c_inv = la.inv(c)
        logm = multilog(multiprod(multiprod(c_inv, y), multitransp(c_inv)),
                        pos_def=True)
        return multiprod(multiprod(c, logm), multitransp(c)) 
Example 68
Project: Programming-for-Non-Technical-Roles-   Author: PacktPublishing   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 69
Project: Programming-for-Non-Technical-Roles-   Author: PacktPublishing   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 70
Project: klustakwik2   Author: kwikteam   File: linear_algebra.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def cholesky(self):
        M_chol = self.new_with_same_masks()
        lower = None
        if self.block.size:
            block = cholesky(self.block, lower=True)
        else:
            block = zeros_like(self.block)
        if (self.diagonal<=0).any():
            raise LinAlgError
        diagonal = sqrt(self.diagonal)
        return self.new_with_same_masks(block, diagonal) 
Example 71
Project: inference-tools   Author: C-bowman   File: gp_tools.py    MIT License 5 votes vote down vote up
def set_hyperparameters(self, theta):
        self.hyperpars = theta
        self.K_xx = self.cov.build_covariance(theta) + self.sig
        self.L = cholesky(self.K_xx)
        self.alpha = solve_triangular(self.L.T, solve_triangular(self.L, self.y, lower = True)) 
Example 72
Project: inference-tools   Author: C-bowman   File: gp_tools.py    MIT License 5 votes vote down vote up
def marginal_likelihood(self, theta):
        """
        returns the log-marginal likelihood for the supplied hyperparameter values.

        This implementation is based on equation (5.8) from Rasmussen & Williams.
        """
        K_xx = self.cov.build_covariance(theta) + self.sig

        try: # protection against singular matrix error crash
            L = cholesky(K_xx)
            alpha = solve_triangular(L.T, solve_triangular(L, self.y, lower = True))
            return -0.5*dot( self.y.T, alpha ) - log(diagonal(L)).sum()
        except:
            warn('Cholesky decomposition failure in marginal_likelihood')
            return -1e50 
Example 73
Project: linear_neuron   Author: uglyboxer   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 74
Project: facethin   Author: ParkerGod   File: test_linalg.py    GNU General Public License v3.0 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 75
Project: facethin   Author: ParkerGod   File: test_regression.py    GNU General Public License v3.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 76
Project: islam-buddy   Author: hamir   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 77
Project: islam-buddy   Author: hamir   File: test_regression.py    MIT License 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 78
Project: mxnet-lambda   Author: awslabs   File: test_linalg.py    Apache License 2.0 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        a = np.zeros((0, 1, 1), dtype=np.int_).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.float64)
        # for documentation purpose:
        assert_(isinstance(res, np.ndarray))

        a = np.zeros((1, 0, 0), dtype=np.complex64).view(ArraySubclass)
        res = linalg.cholesky(a)
        assert_equal(a.shape, res.shape)
        assert_(res.dtype.type is np.complex64)
        assert_(isinstance(res, np.ndarray)) 
Example 79
Project: mxnet-lambda   Author: awslabs   File: test_regression.py    Apache License 2.0 5 votes vote down vote up
def test_lapack_endian(self):
        # For bug #1482
        a = array([[5.7998084,  -2.1825367],
                   [-2.1825367,   9.85910595]], dtype='>f8')
        b = array(a, dtype='<f8')

        ap = linalg.cholesky(a)
        bp = linalg.cholesky(b)
        assert_array_equal(ap, bp) 
Example 80
Project: vnpy_crypto   Author: birforce   File: irf.py    MIT License 4 votes vote down vote up
def __init__(self, model, P=None, periods=10, order=None, svar=False,
                 vecm=False):
        self.model = model
        self.periods = periods
        self.neqs, self.lags, self.T = model.neqs, model.k_ar, model.nobs

        self.order = order

        if P is None:
            sigma = model.sigma_u

            # TODO, may be difficult at the moment
            # if order is not None:
            #     indexer = [model.get_eq_index(name) for name in order]
            #     sigma = sigma[:, indexer][indexer, :]

            #     if sigma.shape != model.sigma_u.shape:
            #         raise ValueError('variable order is wrong length')

            P = la.cholesky(sigma)

        self.P = P

        self.svar = svar

        self.irfs = model.ma_rep(periods)
        if svar:
            self.svar_irfs = model.svar_ma_rep(periods, P=P)
        else:
            self.orth_irfs = model.orth_ma_rep(periods, P=P)

        self.cum_effects = self.irfs.cumsum(axis=0)
        if svar:
            self.svar_cum_effects = self.svar_irfs.cumsum(axis=0)
        else:
            self.orth_cum_effects = self.orth_irfs.cumsum(axis=0)

        # long-run effects may be infinite for VECMs.
        if not vecm:
            self.lr_effects = model.long_run_effects()
            if svar:
                self.svar_lr_effects = np.dot(model.long_run_effects(), P)
            else:
                self.orth_lr_effects = np.dot(model.long_run_effects(), P)


        # auxiliary stuff
        if vecm:
            self._A = util.comp_matrix(model.var_rep)
        else:
            self._A = util.comp_matrix(model.coefs)