Python numpy.linalg.solve() Examples

The following are code examples for showing how to use numpy.linalg.solve(). 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: LaserTOF   Author: kyleuckert   File: test_linalg.py    MIT License 6 votes vote down vote up
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 2
Project: laplacian-meshes   Author: bmershon   File: Utilities2D.py    GNU General Public License v3.0 6 votes vote down vote up
def getBarycentricCoords(A, B, C, X, checkValidity = True):
    T = np.array( [ [A.x - C.x, B.x - C.x ], [A.y - C.y, B.y - C.y] ] )
    y = np.array( [ [X.x - C.x], [X.y - C.y] ] )
    lambdas = linalg.solve(T, y)
    lambdas = lambdas.flatten()
    lambdas = np.append(lambdas, 1 - (lambdas[0] + lambdas[1]))
    if checkValidity:
        if (lambdas[0] < 0 or lambdas[1] < 0 or lambdas[2] < 0):
            print "ERROR: Not a convex combination; lambda = %s"%lambdas
            print "pointInsideConvexPolygon2D = %s"%pointInsideConvexPolygon2D([A, B, C], X, 0)
            plt.plot([A.x, B.x, C.x, A.x], [A.y, B.y, C.y, A.y], 'r')
            plt.hold(True)
            plt.plot([X.x], [X.y], 'b.')
            plt.show()
        assert (lambdas[0] >= 0 and lambdas[1] >= 0 and lambdas[2] >= 0)
    else:
        lambdas[0] = max(lambdas[0], 0)
        lambdas[1] = max(lambdas[1], 0)
        lambdas[2] = max(lambdas[2], 0)
    return lambdas 
Example 3
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_linalg.py    MIT License 6 votes vote down vote up
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 4
Project: recruit   Author: Frank-qlu   File: test_linalg.py    Apache License 2.0 6 votes vote down vote up
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: GenEML   Author: nirbhayjm   File: model.py    MIT License 6 votes vote down vote up
def update_V(m_opts, m_vars):
    P, N = E_x_omega_col(m_opts, m_vars) # expectation of xi_{nl} for n = i, expecation of omega_{nl} for n = i
    Kappa = PG_col(m_opts, m_vars) # polyagamma kappa_{nl} for n = i
    PN = P*N
    PK = P*Kappa

    for i in range(m_vars['n_labels']):
        PN_i = PN[i][:,np.newaxis]
        PK_i = PK[i]

        sigma = m_vars['U_batch'].T.dot(PN_i*m_vars['U_batch'])# + m_opts['lam_v']*np.eye(m_opts['n_components'])
        x = m_vars['U_batch'].T.dot(PK_i)

        m_vars['sigma_v'][i] = (1-m_vars['gamma'])*m_vars['sigma_v'][i] + m_vars['gamma']*sigma
        m_vars['x_v'][i] = (1-m_vars['gamma'])*m_vars['x_v'][i] + m_vars['gamma']*x
        m_vars['V'][i] = linalg.solve(m_vars['sigma_v'][i], m_vars['x_v'][i]) 
Example 6
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_linalg.py    GNU General Public License v3.0 6 votes vote down vote up
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 7
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 6 votes vote down vote up
def iterSolver(self, X, F, FF_init, nWay, r, orderWays):
        # solve NNLS problems for each factor
        for k in range(nWay):
            curWay = orderWays[k]
            ways = range(nWay)
            ways.remove(curWay)
            XF = X.uttkrp(F, curWay)
            # Compute the inner-product matrix
            FF = ones((r, r))
            for i in ways:
                FF = FF * FF_init[i]  # (F[i].T.dot(F[i]))
            ow = 0
            Fthis, temp = nnlsm_activeset(FF, XF.T, ow, 1, F[curWay].T)
            F[curWay] = Fthis.T
            FF_init[curWay] = (F[curWay].T.dot(F[curWay]))
        return F, FF_init 
Example 8
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_linalg.py    MIT License 6 votes vote down vote up
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 9
Project: MARRtino-2.0   Author: DaniAffCH   File: test_linalg.py    GNU General Public License v3.0 6 votes vote down vote up
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 10
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_linalg.py    MIT License 6 votes vote down vote up
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 11
Project: vnpy_crypto   Author: birforce   File: test_linalg.py    MIT License 6 votes vote down vote up
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 12
Project: ble5-nrf52-mac   Author: tomasero   File: test_linalg.py    MIT License 6 votes vote down vote up
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 13
Project: Computable   Author: ktraunmueller   File: test_linalg.py    MIT License 6 votes vote down vote up
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 14
Project: poker   Author: surgebiswas   File: test_linalg.py    MIT License 6 votes vote down vote up
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 15
Project: P3_image_processing   Author: latedude2   File: test_linalg.py    MIT License 6 votes vote down vote up
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 16
Project: LaserTOF   Author: kyleuckert   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) 
Example 17
Project: LaserTOF   Author: kyleuckert   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype 
Example 18
Project: LaserTOF   Author: kyleuckert   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 19
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b, tags):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(consistent_subclass(x, b)) 
Example 20
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self, dtype):
        x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
        assert_equal(linalg.solve(x, x).dtype, dtype) 
Example 21
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
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 22
Project: reconstruction   Author: microelly2   File: projectiontools.py    GNU Lesser General Public License v3.0 5 votes vote down vote up
def mittelpunkt(A,B,C,D):
	''' schnitt der Diagonalen AC und BD '''
	M=np.array([A-C,D-B])
	R=D-C
	x=npl.solve(M.T,R)
	(np.dot(M,x)-R)
	S=x[0]*A+(1-x[0])*C
	return S 
Example 23
Project: reconstruction   Author: microelly2   File: projectiontools.py    GNU Lesser General Public License v3.0 5 votes vote down vote up
def mittelpunkt(A,B,C,D):
	''' schnitt der Diagonalen AC und BD '''
	M=np.array([A-C,D-B])
	R=D-C
	x=npl.solve(M.T,R)
	(np.dot(M,x)-R)
	S=x[0]*A+(1-x[0])*C
	return S 
Example 24
Project: recruit   Author: Frank-qlu   File: test_linalg.py    Apache License 2.0 5 votes vote down vote up
def do(self, a, b, tags):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(consistent_subclass(x, b)) 
Example 25
Project: recruit   Author: Frank-qlu   File: test_linalg.py    Apache License 2.0 5 votes vote down vote up
def test_types(self, dtype):
        x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
        assert_equal(linalg.solve(x, x).dtype, dtype) 
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
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 27
Project: GenEML   Author: nirbhayjm   File: model.py    MIT License 5 votes vote down vote up
def update_U(m_opts, m_vars):
    for it in range(m_opts['PG_iters']):
        P, N = E_x_omega_row(m_opts, m_vars) # expectation of xi_{nl} for n = i, expecation of omega_{nl} for n = i
        Kappa = PG_row(m_opts, m_vars) # polyagamma kappa_{nl} for n = i
        PN = P*N
        PK = P*Kappa

        for i in range(m_opts['batch_size']):
            PN_i = PN[i][:,np.newaxis]
            PK_i = PK[i]
        
            sigma = m_vars['V'].T.dot(PN_i*m_vars['V']) + m_opts['lam_u']*np.eye(m_opts['n_components'])
            x = m_vars['V'].T.dot(PK_i) + np.asarray((m_opts['lam_u']*m_vars['W']).dot(m_vars['X_batch'][i].todense().T)).reshape(-1)
            m_vars['U_batch'][i] = linalg.solve(sigma, x) 
Example 28
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_linalg.py    GNU General Public License v3.0 5 votes vote down vote up
def do(self, a, b):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) 
Example 29
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_linalg.py    GNU General Public License v3.0 5 votes vote down vote up
def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype 
Example 30
Project: att   Author: Centre-Alt-Rendiment-Esportiu   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
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0,:])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2) # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 31
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: bench_basic.py    GNU General Public License v3.0 5 votes vote down vote up
def bench_random(self):
        numpy_solve = nl.solve
        scipy_solve = sl.solve
        print()
        print('      Solving system of linear equations')
        print('      ==================================')

        print('      |    contiguous     |   non-contiguous ')
        print('----------------------------------------------')
        print(' size |  scipy  | numpy   |  scipy  | numpy ')

        for size,repeat in [(20,1000),(100,150),(500,2),(1000,1)][:-1]:
            repeat *= 2
            print('%5s' % size, end=' ')
            sys.stdout.flush()

            a = random([size,size])
            # larger diagonal ensures non-singularity:
            for i in range(size):
                a[i,i] = 10*(.1+a[i,i])
            b = random([size])

            print('| %6.2f ' % measure('scipy_solve(a,b)',repeat), end=' ')
            sys.stdout.flush()

            print('| %6.2f ' % measure('numpy_solve(a,b)',repeat), end=' ')
            sys.stdout.flush()

            a = a[-1::-1,-1::-1]  # turn into a non-contiguous array
            assert_(not a.flags['CONTIGUOUS'])

            print('| %6.2f ' % measure('scipy_solve(a,b)',repeat), end=' ')
            sys.stdout.flush()

            print('| %6.2f ' % measure('numpy_solve(a,b)',repeat), end=' ')
            sys.stdout.flush()

            print('   (secs for %s calls)' % (repeat)) 
Example 32
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 5 votes vote down vote up
def poly_to_jacobi(x):
    """x is a poly1d object"""
    xc = x.coeffs
    N = x.order+1
    matrix = zeros(shape=(N,N), dtype=float)
    for i in range(N):
        matrix[N-i-1:N, i] = jacobi(i,a,b).coeffs
    return solve(matrix, xc) 
Example 33
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 5 votes vote down vote up
def RobustContinuousClustering(X, W, offset_mu=100, delta=0.05, eps=1e-4,
                               max_iter=100):
    n_samples = X.shape[0]
    d = X.shape[1]
    chi = norm(X, ord=2)

    U = X.copy()

    L = np.ones((n_samples, n_samples))
    dists = cdist(X, X, metric='euclidean')
    lower_bound_mu = np.max(dists**2)
    mu = offset_mu + lower_bound_mu
    A = computeA(W, L)
    landa = chi / norm(A, ord=2)

    conv_diff = 100
    old_C = 100
    i = 0
    while conv_diff > eps and i < max_iter:
        # Update L
        L = update_L(U, mu)
        # Update A
        A = computeA(W, L)
        M = np.eye(n_samples) + landa * A
        # Update U
        # U = np.dot(X, inv(M))
        U = solve(M, X)
        # Evaluate objective
        C = objective(X, U, W, L, landa, mu)
        print(C)
        conv_diff = np.abs(C - old_C)
        # Update landa and mu
        if i % 4 == 0:
            landa = chi / norm(A, ord=2)
            mu = max(mu / 2, delta / 2)

        # Keep iterating
        old_c = C
        i += 1
    return U 
Example 34
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 35
Project: scikit-tensor-py3   Author: evertrol   File: rescal.py    GNU General Public License v3.0 5 votes vote down vote up
def _updateA(X, A, R, P, Z, lmbdaA, orthogonalize):
    """Update step for A"""
    n, rank = A.shape
    F = np.zeros((n, rank), dtype=A.dtype)
    E = np.zeros((rank, rank), dtype=A.dtype)

    AtA = np.dot(A.T, A)

    for i in range(len(X)):
        F += X[i].dot(np.dot(A, R[i].T)) + X[i].T.dot(np.dot(A, R[i]))
        E += np.dot(R[i], np.dot(AtA, R[i].T)) + np.dot(R[i].T, np.dot(AtA, R[i]))

    # regularization
    I = lmbdaA * np.eye(rank, dtype=A.dtype)

    # attributes
    for i in range(len(Z)):
        F += P[i].dot(Z[i].T)
        E += np.dot(Z[i], Z[i].T)

    # finally compute update for A
    A = solve(I + E.T, F.T).T
    return orth(A) if orthogonalize else A


# ------------------ Update R ------------------------------------------------ 
Example 36
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b, tags):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(consistent_subclass(x, b)) 
Example 37
Project: FUTU_Stop_Loss   Author: BigtoC   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            check(dtype) 
Example 38
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
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 39
Project: SparseSC   Author: microsoft   File: fit_ct.py    MIT License 5 votes vote down vote up
def ct_weights(
    X, V, w_pen, treated_units=None, control_units=None, custom_donor_pool=None
):
    """ 
    Fit the weights using the cross-train gradient approach
    """
    if treated_units is None:
        if control_units is None:
            raise ValueError(
                "At least on of treated_units or control_units is required"
            )
        # Set the treated units to the not-control units
        treated_units = list(set(range(X.shape[0])) - set(control_units))
    if control_units is None:
        control_units = list(set(range(X.shape[0])) - set(treated_units))

    def _calc_W_ct(X_treated, X_control, V, w_pen):
        A = X_control.dot(2 * V).dot(X_control.T) + 2 * w_pen * diag(
            ones(X_control.shape[0])
        )  # 5
        B = (
            X_treated.dot(2 * V).dot(X_control.T).T + 2 * w_pen / X_control.shape[0]
        )  # 6
        try:
            weights = linalg.solve(A, B).T
        except linalg.LinAlgError as exc:
            print("Unique weights not possible.")
            if w_pen == 0:
                print("Try specifying a very small w_pen rather than 0.")
            raise exc
        return weights

    if custom_donor_pool is None:
        weights = _calc_W_ct(X[treated_units, :], X[control_units, :], V, w_pen)
    else:
        weights = np.zeros((len(treated_units), len(control_units)))
        for i, treated_unit in enumerate(treated_units):
            donors = np.where(custom_donor_pool[treated_unit, :])
            weights[i, donors] = _calc_W_ct(X[treated_unit, :], X[donors, :], V, w_pen)

    return weights 
Example 40
Project: MARRtino-2.0   Author: DaniAffCH   File: test_linalg.py    GNU General Public License v3.0 5 votes vote down vote up
def do(self, a, b, tags):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(consistent_subclass(x, b)) 
Example 41
Project: MARRtino-2.0   Author: DaniAffCH   File: test_linalg.py    GNU General Public License v3.0 5 votes vote down vote up
def test_types(self, dtype):
        x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
        assert_equal(linalg.solve(x, x).dtype, dtype) 
Example 42
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
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 43
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) 
Example 44
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype 
Example 45
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 46
Project: vnpy_crypto   Author: birforce   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b, tags):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) 
Example 47
Project: vnpy_crypto   Author: birforce   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype 
Example 48
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
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 49
Project: miRge   Author: mhalushka   File: vector.py    GNU General Public License v3.0 5 votes vote down vote up
def change_basis(coords, new_basis, old_basis):
    '''
    Change the basis of coordinates to a new basis. In a regular structure
    we have the coordinates in the regular cartesian coordinate system. For helix-helix
    orientations, however, we want to express the coordinates in a coordinate system
    defined by the first helix.

    The new basis will consist of the axis of the first helix, one of its twist elements
    (the one closest to the second vector) and a third vector orthogonal to the previous
    two.

    # http://tutorial.math.lamar.edu/Classes/LinAlg/ChangeOfBasis.aspx

    @param coords: The coordinates to transform (array of n elements).
    @param new_basis: The new basis vectors (n x n matrix)
    @param old_basis: The old basis for the coordinates(n x n matrix)
    @return: The new coordinates according to the new basis
    '''
    #assert(len(coords) == len(new_basis))
    #assert(len(new_basis) == len(old_basis))

    dim = len(coords)
    #print "coords:", coords
    standard_coords = np.dot(old_basis.transpose(), coords)
    '''
    #print "standard_coords:", standard_coords
    standard_to_new = inv(new_basis.transpose())
    #print "standard_to_new:", standard_to_new
    new_coords = np.dot(standard_to_new, standard_coords)
    print "new_coords:", new_coords
    '''

    new_coords = nl.solve(new_basis.transpose(), standard_coords)
    #print "new_coords1:", new_coords1

    return new_coords 
Example 50
Project: miRge   Author: mhalushka   File: vector.py    GNU General Public License v3.0 5 votes vote down vote up
def change_basis2(coords, new_basis, old_basis):
    '''
    '''

    dim = len(coords)
    standard_coords = np.dot(old_basis.T, coords)
    new_coords = nl.solve(new_basis.T, standard_coords)

    return new_coords 
Example 51
Project: ble5-nrf52-mac   Author: tomasero   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b, tags):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(consistent_subclass(x, b)) 
Example 52
Project: ble5-nrf52-mac   Author: tomasero   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self, dtype):
        x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
        assert_equal(linalg.solve(x, x).dtype, dtype) 
Example 53
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
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 54
Project: Computable   Author: ktraunmueller   File: linear_algebra.py    MIT License 5 votes vote down vote up
def solve_linear_equations(a, b):
    return linalg.solve(a, b)

# Matrix inversion 
Example 55
Project: Computable   Author: ktraunmueller   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) 
Example 56
Project: Computable   Author: ktraunmueller   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype 
Example 57
Project: Computable   Author: ktraunmueller   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0,:])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2) # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 58
Project: Computable   Author: ktraunmueller   File: bench_basic.py    MIT License 5 votes vote down vote up
def bench_random(self):
        numpy_solve = nl.solve
        scipy_solve = sl.solve
        print()
        print('      Solving system of linear equations')
        print('      ==================================')

        print('      |    contiguous     |   non-contiguous ')
        print('----------------------------------------------')
        print(' size |  scipy  | numpy   |  scipy  | numpy ')

        for size,repeat in [(20,1000),(100,150),(500,2),(1000,1)][:-1]:
            repeat *= 2
            print('%5s' % size, end=' ')
            sys.stdout.flush()

            a = random([size,size])
            # larger diagonal ensures non-singularity:
            for i in range(size):
                a[i,i] = 10*(.1+a[i,i])
            b = random([size])

            print('| %6.2f ' % measure('scipy_solve(a,b)',repeat), end=' ')
            sys.stdout.flush()

            print('| %6.2f ' % measure('numpy_solve(a,b)',repeat), end=' ')
            sys.stdout.flush()

            a = a[-1::-1,-1::-1]  # turn into a non-contiguous array
            assert_(not a.flags['CONTIGUOUS'])

            print('| %6.2f ' % measure('scipy_solve(a,b)',repeat), end=' ')
            sys.stdout.flush()

            print('| %6.2f ' % measure('numpy_solve(a,b)',repeat), end=' ')
            sys.stdout.flush()

            print('   (secs for %s calls)' % (repeat)) 
Example 59
Project: poker   Author: surgebiswas   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(imply(isinstance(b, matrix), isinstance(x, matrix))) 
Example 60
Project: poker   Author: surgebiswas   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self):
        def check(dtype):
            x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
            assert_equal(linalg.solve(x, x).dtype, dtype)
        for dtype in [single, double, csingle, cdouble]:
            yield check, dtype 
Example 61
Project: poker   Author: surgebiswas   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_0_size(self):
        class ArraySubclass(np.ndarray):
            pass
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 62
Project: morphMan   Author: KVSlab   File: common.py    GNU General Public License v3.0 5 votes vote down vote up
def compute_least_square_plane(cl_points, region_points):
    """
    Find the least squares plane through
    the points in P and approximating the points
    in Z.

    Args:
        cl_points (ndarray): Array of points to approximate.
        region_points (ndarray): Array of points used as constraints.

    Returns:
        n (ndarray): Normal vector to plane.
    """
    # Defined matrices
    b = np.ones(len(cl_points))
    d = np.ones(len(region_points))

    # Create complete matrix
    ata = np.transpose(cl_points).dot(cl_points)
    m0 = np.c_[ata, np.transpose(region_points)]
    m1 = np.c_[region_points, np.zeros((len(region_points), len(region_points)))]
    m = np.r_[m0, m1]
    y = np.r_[np.transpose(cl_points).dot(b), d]

    # Solve system
    x = la.solve(m, y)
    a, b, c = x[0], x[1], x[2]
    n = np.array([a, b, c])
    n = n / la.norm(n)

    return n 
Example 63
Project: pygeom   Author: Xero64   File: cubicspline2d.py    MIT License 5 votes vote down vote up
def calc_d2r_closed(self):
        u"""This function calculates the curvature of a closed spline."""
        n = self.npnts
        inda = [i-1 for i in range(n)]
        indb = [i for i in range(n)]
        inda[0] = n-1
        pnl_dx = [pnl.vec.x/pnl.length for pnl in self.pnls]
        pnl_dy = [pnl.vec.y/pnl.length for pnl in self.pnls]
        del_dx = [0.]*n
        del_dy = [0.]*n
        for i in range(n):
            del_dx[i] = pnl_dx[indb[i]]-pnl_dx[inda[i]]
            del_dy[i] = pnl_dy[indb[i]]-pnl_dy[inda[i]]
        A = zeros((n, n))
        B = zeros((n, 2))
        for i in range(n):
            sA = self.pnls[inda[i]].length
            sB = self.pnls[indb[i]].length
            if i-1 < 0:
                A[i, n-1] = sA/6
            else:
                A[i, i-1] = sA/6
            A[i, i] = (sA+sB)/3
            if i+1 > n-1:
                A[i, 0] = sB/6
            else:
                A[i,i+1] = sB/6
            B[i, 0] = del_dx[i]
            B[i, 1] = del_dy[i]
        X = solve(A, B)
        d2x = [X[i, 0] for i in range(n)]
        d2y = [X[i, 1] for i in range(n)]
        d2r = [Vector2D(d2x[i], d2y[i]) for i in range(self.npnts)]
        return d2r 
Example 64
Project: pygeom   Author: Xero64   File: cubicspline.py    MIT License 5 votes vote down vote up
def calc_d2r_closed(self):
        u"""This function calculates the curvature of a closed spline."""
        n = self.npnts
        inda = [i-1 for i in range(n)]
        indb = [i for i in range(n)]
        inda[0] = n-1
        pnl_dx = [pnl.vec.x/pnl.length for pnl in self.pnls]
        pnl_dy = [pnl.vec.y/pnl.length for pnl in self.pnls]
        pnl_dz = [pnl.vec.z/pnl.length for pnl in self.pnls]
        del_dx = [0.]*n
        del_dy = [0.]*n
        del_dz = [0.]*n
        for i in range(n):
            del_dx[i] = pnl_dx[indb[i]]-pnl_dx[inda[i]]
            del_dy[i] = pnl_dy[indb[i]]-pnl_dy[inda[i]]
            del_dz[i] = pnl_dz[indb[i]]-pnl_dz[inda[i]]
        A = zeros((n, n))
        B = zeros((n, 3))
        for i in range(n):
            sA = self.pnls[inda[i]].length
            sB = self.pnls[indb[i]].length
            if i-1 < 0:
                A[i, n-1] = sA/6
            else:
                A[i, i-1] = sA/6
            A[i, i] = (sA+sB)/3
            if i+1 > n-1:
                A[i, 0] = sB/6
            else:
                A[i, i+1] = sB/6
            B[i, 0] = del_dx[i]
            B[i, 1] = del_dy[i]
            B[i, 2] = del_dz[i]
        X = solve(A, B)
        d2x = [X[i, 0] for i in range(n)]
        d2y = [X[i, 1] for i in range(n)]
        d2z = [X[i, 2] for i in range(n)]
        d2r = [Vector(d2x[i], d2y[i], d2z[i]) for i in range(self.npnts)]
        return d2r 
Example 65
Project: pygeom   Author: Xero64   File: matrixvector.py    MIT License 5 votes vote down vote up
def solve_matrix_vector(a: matrix, b: MatrixVector):
    from numpy.linalg import solve
    newb = zeros((b.shape[0], b.shape[1]*3))
    for i in range(b.shape[1]):
        newb[:, 3*i+0] = b[:, i].x
        newb[:, 3*i+1] = b[:, i].y
        newb[:, 3*i+2] = b[:, i].z
    newc = solve(a, newb)
    c = zero_matrix_vector(b.shape)
    for i in range(b.shape[1]):
        c[:, i] = MatrixVector(newc[:, 3*i+0], newc[:, 3*i+1], newc[:, 3*i+2])
    return c 
Example 66
Project: egk   Author: mlds-lab   File: gp.py    MIT License 5 votes vote down vote up
def posterior_mean(t_train, y_train, t_test, parms, mean_shift=MEAN_SHIFT):
    a, b, c = parms

    K_train = reg_cov_mat(t_train, a, b, c)
    K_test_train = cov_mat(t_test, t_train, a, b)

    if mean_shift:
        mu = np.mean(y_train)
        post_mean = mu + K_test_train.dot(solve(K_train, y_train - mu))
    else:
        post_mean = K_test_train.dot(solve(K_train, y_train))

    return post_mean 
Example 67
Project: egk   Author: mlds-lab   File: gp.py    MIT License 5 votes vote down vote up
def _dfunc(dk, cov, Kinv_y):
    return (Kinv_y.dot(dk).dot(Kinv_y) - np.trace(solve(cov, dk))) * 0.5 
Example 68
Project: P3_image_processing   Author: latedude2   File: test_linalg.py    MIT License 5 votes vote down vote up
def do(self, a, b, tags):
        x = linalg.solve(a, b)
        assert_almost_equal(b, dot_generalized(a, x))
        assert_(consistent_subclass(x, b)) 
Example 69
Project: P3_image_processing   Author: latedude2   File: test_linalg.py    MIT License 5 votes vote down vote up
def test_types(self, dtype):
        x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
        assert_equal(linalg.solve(x, x).dtype, dtype) 
Example 70
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
        # Test system of 0x0 matrices
        a = np.arange(8).reshape(2, 2, 2)
        b = np.arange(6).reshape(1, 2, 3).view(ArraySubclass)

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

        # Test errors for non-square and only b's dimension being 0
        assert_raises(linalg.LinAlgError, linalg.solve, a[:, 0:0, 0:1], b)
        assert_raises(ValueError, linalg.solve, a, b[:, 0:0, :])

        # Test broadcasting error
        b = np.arange(6).reshape(1, 3, 2)  # broadcasting error
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])

        # Test zero "single equations" with 0x0 matrices.
        b = np.arange(2).reshape(1, 2).view(ArraySubclass)
        expected = linalg.solve(a, b)[:, 0:0]
        result = linalg.solve(a[:, 0:0, 0:0], b[:, 0:0])
        assert_array_equal(result, expected)
        assert_(isinstance(result, ArraySubclass))

        b = np.arange(3).reshape(1, 3)
        assert_raises(ValueError, linalg.solve, a, b)
        assert_raises(ValueError, linalg.solve, a[0:0], b[0:0])
        assert_raises(ValueError, linalg.solve, a[:, 0:0, 0:0], b) 
Example 71
Project: QuantEcon.lectures.code   Author: QuantEcon   File: asset_pricing.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def tree_price(ap):
    """
    Computes the price-dividend ratio of the Lucas tree.

    Parameters
    ----------
    ap: AssetPriceModel
        An instance of AssetPriceModel containing primitives

    Returns
    -------
    v : array_like(float)
        Lucas tree price-dividend ratio

    """
    # == Simplify names, set up matrices  == #
    β, γ, P, y = ap.β, ap.γ, ap.mc.P, ap.mc.state_values
    J = P * ap.g(y)**(1 - γ)

    # == Make sure that a unique solution exists == #
    ap.test_stability(J)

    # == Compute v == #
    I = np.identity(ap.n)
    Ones = np.ones(ap.n)
    v = solve(I - β * J, β * J @ Ones)

    return v 
Example 72
Project: QuantEcon.lectures.code   Author: QuantEcon   File: asset_pricing.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def consol_price(ap, ζ):
    """
    Computes price of a consol bond with payoff ζ

    Parameters
    ----------
    ap: AssetPriceModel
        An instance of AssetPriceModel containing primitives

    ζ : scalar(float)
        Coupon of the console

    Returns
    -------
    p : array_like(float)
        Console bond prices

    """
    # == Simplify names, set up matrices  == #
    β, γ, P, y = ap.β, ap.γ, ap.mc.P, ap.mc.state_values
    M = P * ap.g(y)**(- γ)

    # == Make sure that a unique solution exists == #
    ap.test_stability(M)

    # == Compute price == #
    I = np.identity(ap.n)
    Ones = np.ones(ap.n)
    p = solve(I - β * M, β * ζ * M @ Ones)

    return p 
Example 73
Project: pypiv   Author: jr7   File: peak_detection.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def gaussian2D(window):
    """Real 2D Gaussian interpolation for sub pixel shift"""
    #ref on paper
    w = np.ones((3, 3))*(1./9)
    rhs = np.zeros(6)
    M = np.zeros((6,6))
    for i in [-1, 0, 1]:
        for j in [-1, 0, 1]:
            rhs = rhs     +     np.array([i*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                          j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        i*j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        i*i*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        j*j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                            w[i+1, j+1]*np.log(np.abs(window[i+1, j+1]))], dtype='float')

            M = M + w[i+1, j+1]*np.array([[  i*i,   i*j,   i*i*j,   i*i*i,   i*j*j,   i],
                                          [  i*j,   j*j,   i*j*j,   i*i*j,   j*j*j,   j],
                                          [i*i*j, i*j*j, i*i*j*j, i*i*i*j, i*j*j*j, i*j],
                                          [i*i*i, i*i*j, i*i*i*j, i*i*i*i, i*i*j*j, i*i],
                                          [i*j*j, j*j*j, i*j*j*j, i*i*j*j, j*j*j*j, j*j],
                                          [    i,     j,     i*j,     i*i,     j*j,   1]], dtype='float')
    solution = nl.solve(M, rhs)

    dx = (    solution[2]*solution[1] - 2.0*solution[0]*solution[4])/ \
         (4.0*solution[3]*solution[4] -     solution[2]*solution[2])

    dy = (    solution[2]*solution[0] - 2.0*solution[1]*solution[3])/ \
         (4.0*solution[3]*solution[4] -     solution[2]*solution[2])

    return dx, dy 
Example 74
Project: Old-school-processing   Author: cianfrocco-lab   File: imagefilter.py    MIT License 4 votes vote down vote up
def planeRegression(imgarray, msg=True):
	"""
	performs a two-dimensional linear regression and subtracts it from an image
	essentially a fast high pass filter
	"""


	### create index arrays, e.g., [1, 2, 3, 4, 5, ..., N]
	def retx(y,x):
		return x
	def rety(y,x):
		return y
	xarray = numpy.fromfunction(retx, imgarray.shape, dtype=numpy.float32)
	yarray = numpy.fromfunction(rety, imgarray.shape, dtype=numpy.float32)
	xsize = imgarray.shape[0]
	ysize = imgarray.shape[1]
	xarray = xarray/(ysize-1.0) - 0.5
	yarray = yarray/(xsize-1.0) - 0.5

	### get running sums
	count =  float(xsize*ysize)
	xsum =   xarray.sum()
	xsumsq = (xarray*xarray).sum()
	ysum =   yarray.sum()
	ysumsq = (yarray*yarray).sum()
	xysum =  (xarray*yarray).sum()
	xzsum =  (xarray*imgarray).sum()
	yzsum =  (yarray*imgarray).sum()
	zsum =   imgarray.sum()
	zsumsq = (imgarray*imgarray).sum()

	### create linear algebra matrices
	leftmat = numpy.array( [[xsumsq, xysum, xsum], [xysum, ysumsq, ysum], [xsum, ysum, count]], dtype=numpy.float64)
	rightmat = numpy.array( [xzsum, yzsum, zsum], dtype=numpy.float64)

	### solve eigen vectors
	resvec = linalg.solve(leftmat,rightmat)

	### show solution
	if msg is True:
		apDisplay.printMsg("plane_regress: x-slope: %.3f, y-slope: %.3f, xy-intercept: %.3f"
			%(resvec[0], resvec[1], resvec[2]))

	### subtract plane from array
	newarray = imgarray - xarray*resvec[0] - yarray*resvec[1] - resvec[2]
	return newarray

#========================= 
Example 75
Project: pyqmc   Author: WagnerGroup   File: test_obdm.py    MIT License 4 votes vote down vote up
def test():

    mol = gto.M(
        atom="Li 0. 0. 0.; Li 0. 0. 1.5", basis="sto-3g", unit="bohr", verbose=0
    )
    mf = scf.RHF(mol).run()

    # Lowdin orthogonalized AO basis.
    lowdin = lo.orth_ao(mol, "lowdin")

    # MOs in the Lowdin basis.
    mo = solve(lowdin, mf.mo_coeff)

    # make AO to localized orbital coefficients.
    mfobdm = mf.make_rdm1(mo, mf.mo_occ)

    ### Test OBDM calculation.
    nconf = 500
    nsteps = 400
    obdm_steps = 4
    warmup = 15
    wf = PySCFSlaterUHF(mol, mf)
    configs = initial_guess(mol, nconf)
    energy = EnergyAccumulator(mol)
    obdm = OBDMAccumulator(mol=mol, orb_coeff=lowdin, nsweeps=1)
    obdm_up = OBDMAccumulator(mol=mol, orb_coeff=lowdin, nsweeps=1, spin=0)
    obdm_down = OBDMAccumulator(mol=mol, orb_coeff=lowdin, nsweeps=1, spin=1)

    df, coords = vmc(
        wf,
        configs,
        nsteps=nsteps,
        accumulators={
            "energy": energy,
            "obdm": obdm,
            "obdm_up": obdm_up,
            "obdm_down": obdm_down,
        },
    )
    df = DataFrame(df)

    obdm_est = {}
    for k in ["obdm", "obdm_up", "obdm_down"]:
        avg_norm = np.array(df.loc[warmup:, k + "norm"].values.tolist()).mean(axis=0)
        avg_obdm = np.array(df.loc[warmup:, k + "value"].values.tolist()).mean(axis=0)
        obdm_est[k] = normalize_obdm(avg_obdm, avg_norm)

    print("Average OBDM(orb,orb)", obdm_est["obdm"].diagonal().round(3))
    print("mf obdm", mfobdm.diagonal().round(3))
    assert np.max(np.abs(obdm_est["obdm"] - mfobdm)) < 0.05
    print(obdm_est["obdm_up"].diagonal().round(3))
    print(obdm_est["obdm_down"].diagonal().round(3))
    assert np.mean(np.abs(obdm_est["obdm_up"] + obdm_est["obdm_down"] - mfobdm)) < 0.05 
Example 76
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 4 votes vote down vote up
def normalEqComb(AtA, AtB, PassSet=None):
    """ Solve many systems of linear equations using combinatorial grouping.

    M. H. Van Benthem and M. R. Keenan, J. Chemometrics 2004; 18: 441-450

    Parameters
    ----------
    AtA : numpy.array, shape (n,n)
    AtB : numpy.array, shape (n,k)

    Returns
    -------
    Z : numpy.array, shape (n,k) - solution
    """
    if AtB.size == 0:
        Z = np.zeros([])
    elif (PassSet == None) or np.all(PassSet):
        Z = nla.solve(AtA, AtB)
    else:
        Z = np.zeros(AtB.shape)
        if PassSet.shape[1] == 1:
            if np.any(PassSet):
                cols = PassSet.nonzero()[0]
                Z[cols] = nla.solve(AtA[np.ix_(cols, cols)], AtB[cols])
        else:
            #
            # Both _column_group_loop() and _column_group_recursive() work well.
            # Based on preliminary testing,
            # _column_group_loop() is slightly faster for tiny k(<10), but
            # _column_group_recursive() is faster for large k's.
            #
            grps = _column_group_recursive(PassSet)
            for gr in grps:
                cols = PassSet[:, gr[0]].nonzero()[0]
                if cols.size > 0:
                    ix1 = np.ix_(cols, gr)
                    ix2 = np.ix_(cols, cols)
                    #
                    # scipy.linalg.cho_solve can be used instead of numpy.linalg.solve.
                    # For small n(<200), numpy.linalg.solve appears faster, whereas
                    # for large n(>500), scipy.linalg.cho_solve appears faster.
                    # Usage example of scipy.linalg.cho_solve:
                    # Z[ix1] = sla.cho_solve(sla.cho_factor(AtA[ix2]),AtB[ix1])
                    #
                    Z[ix1] = nla.solve(AtA[ix2], AtB[ix1])
    return Z 
Example 77
Project: Echobase   Author: akhambhati   File: nnls.py    GNU General Public License v3.0 4 votes vote down vote up
def normal_eq_comb(AtA, AtB, PassSet=None):
    """ Solve many systems of linear equations using combinatorial grouping.

    M. H. Van Benthem and M. R. Keenan, J. Chemometrics 2004; 18: 441-450

    Parameters
    ----------
    AtA : numpy.array, shape (n,n)
    AtB : numpy.array, shape (n,k)

    Returns
    -------
    (Z,num_cholesky,num_eq)
    Z : numpy.array, shape (n,k) - solution
    num_cholesky : int - the number of unique cholesky decompositions done
    num_eq: int - the number of systems of linear equations solved
    """
    num_cholesky = 0
    num_eq = 0
    if AtB.size == 0:
        Z = np.zeros([])
    elif (PassSet is None) or np.all(PassSet):
        Z = nla.solve(AtA, AtB)
        num_cholesky = 1
        num_eq = AtB.shape[1]
    else:
        Z = np.zeros(AtB.shape)
        if PassSet.shape[1] == 1:
            if np.any(PassSet):
                cols = PassSet.nonzero()[0]
                Z[cols] = nla.solve(AtA[np.ix_(cols, cols)], AtB[cols])
                num_cholesky = 1
                num_eq = 1
        else:
            #
            # Both _column_group_loop() and _column_group_recursive() work well.
            # Based on preliminary testing,
            # _column_group_loop() is slightly faster for tiny k(<10), but
            # _column_group_recursive() is faster for large k's.
            #
            grps = _column_group_recursive(PassSet)
            for gr in grps:
                cols = PassSet[:, gr[0]].nonzero()[0]
                if cols.size > 0:
                    ix1 = np.ix_(cols, gr)
                    ix2 = np.ix_(cols, cols)
                    #
                    # scipy.linalg.cho_solve can be used instead of numpy.linalg.solve.
                    # For small n(<200), numpy.linalg.solve appears faster, whereas
                    # for large n(>500), scipy.linalg.cho_solve appears faster.
                    # Usage example of scipy.linalg.cho_solve:
                    # Z[ix1] = sla.cho_solve(sla.cho_factor(AtA[ix2]),AtB[ix1])
                    #
                    Z[ix1] = nla.solve(AtA[ix2], AtB[ix1])
                    num_cholesky += 1
                    num_eq += len(gr)
                    num_eq += len(gr)
    return Z, num_cholesky, num_eq 
Example 78
Project: egk   Author: mlds-lab   File: full_marginal.py    MIT License 4 votes vote down vote up
def compute_means_covs(ts, t_ref, gp_parms, winsize=0, mean_shift=True):
    """
    Compute the posterior GP means and covariance matrices.

    ts: time series
    t_ref: reference time points the posterior GP is marginalized over
    gp_parms: GP hyperparameters and the noise term
    winsize: window size, 0 for using the full Gaussian (over t_ref)
    """
    a, b, c = gp_parms
    K_test = cov_mat(t_ref, t_ref, a, b)

    n_ts = len(ts)
    n_sample = len(t_ref)
    if winsize == 0:
        post_means = np.empty((n_ts, n_sample))
        post_covs = np.empty((n_ts, n_sample, n_sample))
    else:
        n_kernel = n_sample - winsize + 1
        post_means = np.empty((n_ts, n_kernel, winsize))
        post_covs = np.empty((n_ts, n_kernel, winsize, winsize))

    for idx, (t, y) in enumerate(ts):
        mean_list, cov_list = [], []
        K_train = reg_cov_mat(t, a, b, c)
        K_train_test = cov_mat(t, t_ref, a, b)
        Ktr_inv_Ktt = solve(K_train, K_train_test)
        if mean_shift:
            mu = np.mean(y)
            mean_test = mu + Ktr_inv_Ktt.T.dot(y - mu)
        else:
            mean_test = Ktr_inv_Ktt.T.dot(y)

        full_cov = K_test - K_train_test.T.dot(Ktr_inv_Ktt)
        if winsize == 0:
            post_means[idx] = mean_test
            post_covs[idx] = full_cov
        else:
            for i in xrange(n_sample - winsize + 1):
                post_means[idx, i] = mean_test[i:(i + winsize)]
                post_covs[idx, i] = full_cov[i:(i + winsize), i:(i + winsize)]

    return post_means, post_covs 
Example 79
Project: egk   Author: mlds-lab   File: gp.py    MIT License 4 votes vote down vote up
def posterior_mean_cov(t_train, y_train, t_test, parms, mean_shift=MEAN_SHIFT):
    """Prediction for Gaussian process regression

    Returns the predictive mean and covariance matrix

    Parameters
    ----------
    t_train : array_like, shape (n_training_samples, )
              time points of training samples

    y_train : array_like, shape (n_training_samples, )
              values of corresponding time points of training samples

    t_test : array_like, shape (n_test_samples, )
             time points of test samples

    parms : tuple, length 3
            hyperparameters (a, b, c) for covariance function of Gaussian
            process.
            [K(t)]_ij = a * exp(-b * (t_i - t_j)^2) + c * I(i == j)

    Returns
    -------
    mean_test : array, shape (n_test_samples, )
                predictive mean

    cov_test : array, shape (n_test_samples, n_test_samples)
               predictive covariance matrix
    """

    #assert t_train.shape[0] == y_train.shape[0]
    a, b, c = parms

    K_train = reg_cov_mat(t_train, a, b, c)
    K_train_test = cov_mat(t_train, t_test, a, b)

    Ktr_inv_Ktt = solve(K_train, K_train_test)

    if mean_shift:
        mu = np.mean(y_train)
        mean_test = mu + Ktr_inv_Ktt.T.dot(y_train - mu)
    else:
        mean_test = Ktr_inv_Ktt.T.dot(y_train)

    cov_test = cov_mat(t_test, t_test, a, b) - K_train_test.T.dot(Ktr_inv_Ktt)

    return mean_test, cov_test 
Example 80
Project: egk   Author: mlds-lab   File: gp.py    MIT License 4 votes vote down vote up
def pointwise_posterior_mean_var(t_train, y_train, t_test, parms,
                                 mean_shift=MEAN_SHIFT):
    """Prediction for Gaussian process regression

    Returns the pointwise predictive mean and variance

    Parameters
    ----------
    t_train : array_like, shape (n_training_samples, )
              time points of training samples

    y_train : array_like, shape (n_training_samples, )
              values of corresponding time points of training samples

    t_test : array_like, shape (n_test_samples, )
             time points of test samples

    parms : tuple, length 3
            hyperparameters (a, b, c) for covariance function of Gaussian
            process.
            [K(t)]_ij = a * exp(-b * (t_i - t_j)^2) + c * I(i == j)

    Returns
    -------
    mean_test : array, shape (n_test_samples, )
                predictive mean

    var_test : array, shape (n_test_samples, )
               predictive variance
    """

    #assert t_train.shape[0] == y_train.shape[0]
    a, b, c = parms

    if len(y_train) == 0:
        return np.zeros_like(t_test), np.ones_like(t_test) * a

    K_train = reg_cov_mat(t_train, a, b, c)
    K_train_test = cov_mat(t_train, t_test, a, b)

    Ktr_inv_Ktt = solve(K_train, K_train_test)

    if mean_shift:
        mu = np.mean(y_train)
        mean_test = mu + Ktr_inv_Ktt.T.dot(y_train - mu)
    else:
        mean_test = Ktr_inv_Ktt.T.dot(y_train)

    var_test = a * np.ones_like(t_test) - (K_train_test * Ktr_inv_Ktt).sum(axis=0)

    return mean_test, var_test