Python numpy.linalg() Examples

The following are code examples for showing how to use numpy.linalg(). 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: Parallel.GAMIT   Author: demiangomez   File: pyVoronoi.py    GNU General Public License v3.0 6 votes vote down vote up
def produce_array_Voronoi_vertices_on_sphere_surface(facet_coordinate_array_Delaunay_triangulation,sphere_radius,sphere_centroid):
    '''Return shape (N,3) array of coordinates for the vertices of the Voronoi diagram on the sphere surface given a shape (N,3,3) array of Delaunay triangulation vertices.'''
    assert facet_coordinate_array_Delaunay_triangulation.shape[1:] == (3,3), "facet_coordinate_array_Delaunay_triangulation should have shape (N,3,3)."
    #draft numpy vectorized workflow to avoid Python for loop
    facet_normals_array = numpy.cross(facet_coordinate_array_Delaunay_triangulation[...,1,...] - facet_coordinate_array_Delaunay_triangulation[...,0,...],facet_coordinate_array_Delaunay_triangulation[...,2,...] - facet_coordinate_array_Delaunay_triangulation[...,0,...])
    facet_normal_magnitudes = numpy.linalg.norm(facet_normals_array,axis=1)
    facet_normal_unit_vector_array = facet_normals_array / numpy.column_stack((facet_normal_magnitudes,facet_normal_magnitudes,facet_normal_magnitudes))
    #try to ensure that facet normal faces the correct direction (i.e., out of sphere)
    triangle_centroid_array = numpy.average(facet_coordinate_array_Delaunay_triangulation,axis=1)
    #normalize the triangle_centroid to unit sphere distance for the purposes of the following directionality check
    array_triangle_centroid_spherical_coords = convert_cartesian_array_to_spherical_array(triangle_centroid_array)
    array_triangle_centroid_spherical_coords[...,0] = 1.0
    triangle_centroid_array = convert_spherical_array_to_cartesian_array(array_triangle_centroid_spherical_coords)
    #the Euclidean distance between the triangle centroid and the facet normal should be smaller than the sphere centroid to facet normal distance, otherwise, need to invert the vector
    triangle_to_normal_distance_array = numpy.linalg.norm(triangle_centroid_array - facet_normal_unit_vector_array,axis=1)
    sphere_centroid_to_normal_distance_array = numpy.linalg.norm(sphere_centroid-facet_normal_unit_vector_array,axis=1)
    delta_value_array = sphere_centroid_to_normal_distance_array - triangle_to_normal_distance_array
    facet_normal_unit_vector_array[delta_value_array < -0.1] *= -1.0 #need to rotate the vector so that it faces out of the circle
    facet_normal_unit_vector_array *= sphere_radius #adjust for radius of sphere
    array_Voronoi_vertices = facet_normal_unit_vector_array
    assert array_Voronoi_vertices.shape[1] == 3, "The array of Voronoi vertices on the sphere should have shape (N,3)."
    return array_Voronoi_vertices 
Example 2
Project: SparseSC   Author: microsoft   File: fit_fast.py    MIT License 6 votes vote down vote up
def _weights(V , X_treated, X_control, w_pen):
    V = np.diag(V) #make square
    #weights = np.zeros((X_control.shape[0], X_treated.shape[0]))
    w_pen_mat = 2 * w_pen * np.diag(np.ones(X_control.shape[0]))
    A = X_control.dot(2 * V).dot(X_control.T) + w_pen_mat  # 5
    B = (
        X_treated.dot(2 * V).dot(X_control.T).T + 2 * w_pen / X_control.shape[0]
    )  # 6
    try:
        b = scipy.linalg.solve(A, B)
    except scipy.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 b 
Example 3
Project: maze-builder   Author: kcsaff   File: mesh.py    MIT License 6 votes vote down vote up
def __init__(self, mesh, transform_out=None, transform_in=None):
        transform_out = numpy.matrix(transform_out) if transform_out is not None else None
        transform_in = numpy.matrix(transform_in) if transform_in is not None else None

        if transform_in is None and transform_out is None:
            transform_in = numpy.identity(3)
            transform_out = numpy.identity(3)
        elif transform_in is None:
            try:
                transform_in = numpy.linalg.inv(transform_out)
            except:
                transform_in = None
        elif transform_out is None:
            try:
                transform_out = numpy.linalg.inv(transform_in)
            except:
                transform_out = None

        self.transform_out, self.transform_in = transform_out, transform_in

        super().__init__(
            mesh,
            warp_in=lambda vertex: self.transform_in.dot(vertex).tolist()[0][:3] if self.transform_in else None,
            warp_out=lambda vertex: self.transform_out.dot(vertex).tolist()[0][:3] if self.transform_out else None,
        ) 
Example 4
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    MIT License 6 votes vote down vote up
def test_eigvalsh():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the geigvalsh op.")
    import scipy.linalg

    A = theano.tensor.dmatrix('a')
    B = theano.tensor.dmatrix('b')
    f = function([A, B], eigvalsh(A, B))

    rng = numpy.random.RandomState(utt.fetch_seed())
    a = rng.randn(5, 5)
    a = a + a.T
    for b in [10 * numpy.eye(5, 5) + rng.randn(5, 5)]:
        w = f(a, b)
        refw = scipy.linalg.eigvalsh(a, b)
        numpy.testing.assert_array_almost_equal(w, refw)

    # We need to test None separatly, as otherwise DebugMode will
    # complain, as this isn't a valid ndarray.
    b = None
    B = theano.tensor.NoneConst
    f = function([A], eigvalsh(A, B))
    w = f(a)
    refw = scipy.linalg.eigvalsh(a, b)
    numpy.testing.assert_array_almost_equal(w, refw) 
Example 5
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    MIT License 6 votes vote down vote up
def test_perform(self):
        if not imported_scipy:
            raise SkipTest('kron tests need the scipy package to be installed')

        for shp0 in [(2,), (2, 3), (2, 3, 4), (2, 3, 4, 5)]:
            x = tensor.tensor(dtype='floatX',
                              broadcastable=(False,) * len(shp0))
            a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX)
            for shp1 in [(6,), (6, 7), (6, 7, 8), (6, 7, 8, 9)]:
                if len(shp0) + len(shp1) == 2:
                    continue
                y = tensor.tensor(dtype='floatX',
                                  broadcastable=(False,) * len(shp1))
                f = function([x, y], kron(x, y))
                b = self.rng.rand(*shp1).astype(config.floatX)
                out = f(a, b)
                # Newer versions of scipy want 4 dimensions at least,
                # so we have to add a dimension to a and flatten the result.
                if len(shp0) + len(shp1) == 3:
                    scipy_val = scipy.linalg.kron(
                        a[numpy.newaxis, :], b).flatten()
                else:
                    scipy_val = scipy.linalg.kron(a, b)
                utt.assert_allclose(out, scipy_val) 
Example 6
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    MIT License 6 votes vote down vote up
def test_numpy_compare(self):
        rng = numpy.random.RandomState(utt.fetch_seed())

        M = tensor.matrix("A", dtype=theano.config.floatX)
        V = tensor.vector("V", dtype=theano.config.floatX)

        a = rng.rand(4, 4).astype(theano.config.floatX)
        b = rng.rand(4).astype(theano.config.floatX)

        A = (   [None, 'fro', 'inf', '-inf', 1, -1, None, 'inf', '-inf', 0, 1, -1, 2, -2],
                [M, M, M, M, M, M, V, V, V, V, V, V, V, V],
                [a, a, a, a, a, a, b, b, b, b, b, b, b, b],
                [None, 'fro', inf, -inf, 1, -1, None, inf, -inf, 0, 1, -1, 2, -2])

        for i in range(0, 14):
            f = function([A[1][i]], norm(A[1][i], A[0][i]))
            t_n = f(A[2][i])
            n_n = numpy.linalg.norm(A[2][i], A[3][i])
            assert _allclose(n_n, t_n) 
Example 7
Project: vnpy_crypto   Author: birforce   File: linalg.py    MIT License 6 votes vote down vote up
def logdet_symm(m, check_symm=False):
    """
    Return log(det(m)) asserting positive definiteness of m.

    Parameters
    ----------
    m : array-like
        2d array that is positive-definite (and symmetric)

    Returns
    -------
    logdet : float
        The log-determinant of m.
    """
    from scipy import linalg
    if check_symm:
        if not np.all(m == m.T):  # would be nice to short-circuit check
            raise ValueError("m is not symmetric.")
    c, _ = linalg.cho_factor(m, lower=True)
    return 2*np.sum(np.log(c.diagonal())) 
Example 8
Project: OpenFermion   Author: quantumlib   File: _davidson_test.py    Apache License 2.0 6 votes vote down vote up
def setUp(self):
        """Sets up all variables needed for Davidson class."""
        dimension = 10
        matrix = generate_matrix(dimension)

        def mat_vec(vec):
            """Trivial matvec with a numpy matrix."""
            return numpy.dot(matrix, vec)

        self.linear_operator = scipy.sparse.linalg.LinearOperator(
            (dimension, dimension), matvec=mat_vec)
        self.diagonal = numpy.diag(matrix)

        self.davidson = Davidson(linear_operator=self.linear_operator,
                                 linear_operator_diagonal=self.diagonal)

        self.matrix = matrix
        self.initial_guess = numpy.eye(self.matrix.shape[0], 10)

        self.eigen_values = numpy.array([
            1.15675714, 1.59132505, 2.62268014, 4.44533793, 5.3722743,
            5.54393114, 7.73652405, 8.50089897, 9.4229309, 15.54405993,
        ]) 
Example 9
Project: OpenFermion   Author: quantumlib   File: _davidson_test.py    Apache License 2.0 6 votes vote down vote up
def setUp(self):
        """Sets up all variables needed for SparseDavidson class."""
        logging.basicConfig(level=logging.INFO)
        self.dimension = 1000
        self.sparse_matrix = generate_sparse_matrix(self.dimension)
        self.davidson_options = DavidsonOptions(max_subspace=100,
                                                max_iterations=50,
                                                real_only=True)

        # Checks for built-in eigh() function.
        self.eigen_values, self.eigen_vectors = numpy.linalg.eigh(
            self.sparse_matrix)
        self.assertAlmostEqual(get_difference(
            self.sparse_matrix, self.eigen_values, self.eigen_vectors), 0)

        # Makes sure eigenvalues are sorted.
        self.eigen_values = sorted(self.eigen_values) 
Example 10
Project: OpenFermion   Author: quantumlib   File: _sparse_tools.py    Apache License 2.0 6 votes vote down vote up
def variance(operator, state):
    """Compute variance of operator with a state.

    Args:
        operator(scipy.sparse.spmatrix or scipy.sparse.linalg.LinearOperator):
            The operator whose expectation value is desired.
        state(numpy.ndarray or scipy.sparse.spmatrix): A numpy array
            representing a pure state or a sparse matrix representing a density
            matrix.

    Returns:
        A complex number giving the variance.

    Raises:
        ValueError: Input state has invalid format.
    """
    return (expectation(operator ** 2, state) -
            expectation(operator, state) ** 2) 
Example 11
Project: OpenFermion   Author: quantumlib   File: _sparse_tools.py    Apache License 2.0 6 votes vote down vote up
def get_gap(sparse_operator, initial_guess=None):
    """Compute gap between lowest eigenvalue and first excited state.

    Args:
        sparse_operator (LinearOperator): Operator to find the ground state of.
        initial_guess (ndarray): Initial guess for eigenspace.  A good
            guess dramatically reduces the cost required to converge.
    Returns: A real float giving eigenvalue gap.
    """
    if not is_hermitian(sparse_operator):
        raise ValueError('sparse_operator must be Hermitian.')

    values, _ = scipy.sparse.linalg.eigsh(
        sparse_operator, k=2, v0=initial_guess, which='SA', maxiter=1e7)

    gap = abs(values[1] - values[0])
    return gap 
Example 12
Project: OpenFermion   Author: quantumlib   File: _linear_qubit_operator.py    Apache License 2.0 6 votes vote down vote up
def generate_linear_qubit_operator(qubit_operator, n_qubits=None, options=None):
    """ Generates a LinearOperator from a QubitOperator.

    Args:
        qubit_operator(QubitOperator): A qubit operator to be applied on
            vectors.
        n_qubits(int): The total number of qubits
        options(LinearQubitOperatorOptions): Options for the
            ParallelLinearQubitOperator.
    Returns:
        linear_operator(scipy.sparse.linalg.LinearOperator): A linear operator.
    """
    if options is None:
        linear_operator = LinearQubitOperator(qubit_operator, n_qubits)
    else:
        linear_operator = ParallelLinearQubitOperator(
            qubit_operator, n_qubits, options)
    return linear_operator 
Example 13
Project: OpenFermion   Author: quantumlib   File: _davidson.py    Apache License 2.0 6 votes vote down vote up
def generate_random_vectors(row, col, real_only=False):
    """Generates orthonormal random vectors with col columns.

    Args:
        row(int): Number of rows for the vectors.
        col(int): Number of columns for the vectors.
        real_only(bool): Real vectors or complex ones.

    Returns:
        random_vectors(numpy.ndarray(complex)): Orthonormal random vectors.
    """
    random_vectors = numpy.random.rand(row, col)
    if not real_only:
        random_vectors = random_vectors + numpy.random.rand(row, col) * 1.0j
    random_vectors = scipy.linalg.orth(random_vectors)
    return random_vectors 
Example 14
Project: poker   Author: surgebiswas   File: test_slinalg.py    MIT License 6 votes vote down vote up
def test_eigvalsh():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the geigvalsh op.")
    import scipy.linalg

    A = theano.tensor.dmatrix('a')
    B = theano.tensor.dmatrix('b')
    f = function([A, B], eigvalsh(A, B))

    rng = numpy.random.RandomState(utt.fetch_seed())
    a = rng.randn(5, 5)
    a = a + a.T
    for b in [10 * numpy.eye(5, 5) + rng.randn(5, 5)]:
        w = f(a, b)
        refw = scipy.linalg.eigvalsh(a, b)
        numpy.testing.assert_array_almost_equal(w, refw)

    # We need to test None separatly, as otherwise DebugMode will
    # complain, as this isn't a valid ndarray.
    b = None
    B = theano.tensor.NoneConst
    f = function([A], eigvalsh(A, B))
    w = f(a)
    refw = scipy.linalg.eigvalsh(a, b)
    numpy.testing.assert_array_almost_equal(w, refw) 
Example 15
Project: poker   Author: surgebiswas   File: test_slinalg.py    MIT License 6 votes vote down vote up
def test_perform(self):
        if not imported_scipy:
            raise SkipTest('kron tests need the scipy package to be installed')

        for shp0 in [(2,), (2, 3), (2, 3, 4), (2, 3, 4, 5)]:
            x = tensor.tensor(dtype='floatX',
                              broadcastable=(False,) * len(shp0))
            a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX)
            for shp1 in [(6,), (6, 7), (6, 7, 8), (6, 7, 8, 9)]:
                if len(shp0) + len(shp1) == 2:
                    continue
                y = tensor.tensor(dtype='floatX',
                                  broadcastable=(False,) * len(shp1))
                f = function([x, y], kron(x, y))
                b = self.rng.rand(*shp1).astype(config.floatX)
                out = f(a, b)
                # Newer versions of scipy want 4 dimensions at least,
                # so we have to add a dimension to a and flatten the result.
                if len(shp0) + len(shp1) == 3:
                    scipy_val = scipy.linalg.kron(
                        a[numpy.newaxis, :], b).flatten()
                else:
                    scipy_val = scipy.linalg.kron(a, b)
                utt.assert_allclose(out, scipy_val) 
Example 16
Project: poker   Author: surgebiswas   File: test_nlinalg.py    MIT License 6 votes vote down vote up
def test_numpy_compare(self):
        rng = numpy.random.RandomState(utt.fetch_seed())

        M = tensor.matrix("A", dtype=theano.config.floatX)
        V = tensor.vector("V", dtype=theano.config.floatX)

        a = rng.rand(4, 4).astype(theano.config.floatX)
        b = rng.rand(4).astype(theano.config.floatX)

        A = (   [None, 'fro', 'inf', '-inf', 1, -1, None, 'inf', '-inf', 0, 1, -1, 2, -2],
                [M, M, M, M, M, M, V, V, V, V, V, V, V, V],
                [a, a, a, a, a, a, b, b, b, b, b, b, b, b],
                [None, 'fro', inf, -inf, 1, -1, None, inf, -inf, 0, 1, -1, 2, -2])

        for i in range(0, 14):
            f = function([A[1][i]], norm(A[1][i], A[0][i]))
            t_n = f(A[2][i])
            n_n = numpy.linalg.norm(A[2][i], A[3][i])
            assert _allclose(n_n, t_n) 
Example 17
Project: POVME   Author: POVME   File: pymolecule.py    MIT License 6 votes vote down vote up
def get_dihedral_angle(self, pt1, pt2, pt3, pt4):
        '''Calculates the dihedral angle formed by four points (numpy.array objects).
                
            Arguments:
            pt1 -- A numpy.array (x, y, z) representing the first 3D point.
            pt2 -- A numpy.array (x, y, z) representing the second 3D point.
            pt3 -- A numpy.array (x, y, z) representing the third 3D point.
            pt4 -- A numpy.array (x, y, z) representing the fourth 3D point.
            
            Returns:
            A float containing the dihedral angle between the four points, in radians.
            
            '''
        
        b1 = pt2 - pt1
        b2 = pt3 - pt2
        b3 = pt4 - pt3
        
        b2Xb3 = numpy.cross(b2, b3)
        b1Xb2 = numpy.cross(b1, b2)
        
        b1XMagb2 = numpy.linalg.norm(b2) * b1
        
        return numpy.arctan2(numpy.dot(b1XMagb2, b2Xb3), numpy.dot(b1Xb2, b2Xb3)) 
Example 18
Project: attention-lvcsr   Author: rizar   File: test_slinalg.py    MIT License 6 votes vote down vote up
def test_eigvalsh():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the geigvalsh op.")
    import scipy.linalg

    A = theano.tensor.dmatrix('a')
    B = theano.tensor.dmatrix('b')
    f = function([A, B], eigvalsh(A, B))

    rng = numpy.random.RandomState(utt.fetch_seed())
    a = rng.randn(5, 5)
    a = a + a.T
    for b in [10 * numpy.eye(5, 5) + rng.randn(5, 5)]:
        w = f(a, b)
        refw = scipy.linalg.eigvalsh(a, b)
        numpy.testing.assert_array_almost_equal(w, refw)

    # We need to test None separatly, as otherwise DebugMode will
    # complain, as this isn't a valid ndarray.
    b = None
    B = theano.tensor.NoneConst
    f = function([A], eigvalsh(A, B))
    w = f(a)
    refw = scipy.linalg.eigvalsh(a, b)
    numpy.testing.assert_array_almost_equal(w, refw) 
Example 19
Project: attention-lvcsr   Author: rizar   File: test_slinalg.py    MIT License 6 votes vote down vote up
def test_perform(self):
        if not imported_scipy:
            raise SkipTest('kron tests need the scipy package to be installed')

        for shp0 in [(2,), (2, 3), (2, 3, 4), (2, 3, 4, 5)]:
            x = tensor.tensor(dtype='floatX',
                              broadcastable=(False,) * len(shp0))
            a = numpy.asarray(self.rng.rand(*shp0)).astype(config.floatX)
            for shp1 in [(6,), (6, 7), (6, 7, 8), (6, 7, 8, 9)]:
                if len(shp0) + len(shp1) == 2:
                    continue
                y = tensor.tensor(dtype='floatX',
                                  broadcastable=(False,) * len(shp1))
                f = function([x, y], kron(x, y))
                b = self.rng.rand(*shp1).astype(config.floatX)
                out = f(a, b)
                # Newer versions of scipy want 4 dimensions at least,
                # so we have to add a dimension to a and flatten the result.
                if len(shp0) + len(shp1) == 3:
                    scipy_val = scipy.linalg.kron(
                        a[numpy.newaxis, :], b).flatten()
                else:
                    scipy_val = scipy.linalg.kron(a, b)
                utt.assert_allclose(out, scipy_val) 
Example 20
Project: Parallel.GAMIT   Author: demiangomez   File: pyVoronoi.py    GNU General Public License v3.0 5 votes vote down vote up
def calculate_surface_area_of_planar_polygon_in_3D_space(array_ordered_Voronoi_polygon_vertices):
    '''Based largely on: http://stackoverflow.com/a/12653810
    Use this function when spherical polygon surface area calculation fails (i.e., lots of nearly-coplanar vertices and negative surface area).'''
    #unit normal vector of plane defined by points a, b, and c
    def unit_normal(a, b, c):
        x = numpy.linalg.det([[1,a[1],a[2]],
             [1,b[1],b[2]],
             [1,c[1],c[2]]])
        y = numpy.linalg.det([[a[0],1,a[2]],
             [b[0],1,b[2]],
             [c[0],1,c[2]]])
        z = numpy.linalg.det([[a[0],a[1],1],
             [b[0],b[1],1],
             [c[0],c[1],1]])
        magnitude = (x**2 + y**2 + z**2)**.5
        return (x/magnitude, y/magnitude, z/magnitude)

    #area of polygon poly
    def poly_area(poly):
        '''Accepts a list of xyz tuples.'''
        assert len(poly) >= 3, "Not a polygon (< 3 vertices)."
        total = [0, 0, 0]
        N = len(poly)
        for i in range(N):
            vi1 = poly[i]
            vi2 = poly[(i+1) % N]
            prod = numpy.cross(vi1, vi2)
            total[0] += prod[0]
            total[1] += prod[1]
            total[2] += prod[2]
        result = numpy.dot(total, unit_normal(poly[0], poly[1], poly[2]))
        return abs(result/2)

    list_vertices = [] #need a list of tuples for above function
    for coord in array_ordered_Voronoi_polygon_vertices:
        list_vertices.append(tuple(coord))
    planar_polygon_surface_area = poly_area(list_vertices)
    return planar_polygon_surface_area 
Example 21
Project: tangent   Author: google   File: grads.py    Apache License 2.0 5 votes vote down vote up
def adet(z, x):
  """d|A|/dA = adj(A).T

  See  Jacobi's formula: https://en.wikipedia.org/wiki/Jacobi%27s_formula
  """
  adjugate = numpy.linalg.det(x) * numpy.linalg.pinv(x)
  d[x] = d[z] * numpy.transpose(adjugate)


#
# Tangent adjoints
# 
Example 22
Project: Autoenv   Author: intelligent-control-lab   File: cma_es_lib.py    MIT License 5 votes vote down vote up
def expms(A, eig=np.linalg.eigh):
            """matrix exponential for a symmetric matrix"""
            # TODO: check that this works reliably for low rank matrices
            # first: symmetrize A
            D, B = eig(A)
            return np.dot(B, (np.exp(D) * B).T) 
Example 23
Project: Autoenv   Author: intelligent-control-lab   File: cma_es_lib.py    MIT License 5 votes vote down vote up
def likelihood(x, m=None, Cinv=None, sigma=1, detC=None):
        """return likelihood of x for the normal density N(m, sigma**2 * Cinv**-1)"""
        # testing: MC integrate must be one: mean(p(x_i)) * volume(where x_i are uniformely sampled)
        # for i in xrange(3): print mean([cma.likelihood(20*r-10, dim * [0], None, 3) for r in rand(10000,dim)]) * 20**dim
        if m is None:
            dx = x
        else:
            dx = x - m  # array(x) - array(m)
        n = len(x)
        s2pi = (2 * np.pi)**(n / 2.)
        if Cinv is None:
            return exp(-sum(dx**2) / sigma**2 / 2) / s2pi / sigma**n
        if detC is None:
            detC = 1. / np.linalg.linalg.det(Cinv)
        return  exp(-np.dot(dx, np.dot(Cinv, dx)) / sigma**2 / 2) / s2pi / abs(detC)**0.5 / sigma**n 
Example 24
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_public_api.py    MIT License 5 votes vote down vote up
def test_numpy_linalg():
    bad_results = check_dir(np.linalg)
    assert bad_results == {} 
Example 25
Project: recruit   Author: Frank-qlu   File: test_public_api.py    Apache License 2.0 5 votes vote down vote up
def test_numpy_linalg():
    bad_results = check_dir(np.linalg)
    assert bad_results == {} 
Example 26
Project: PyRiemann   Author: DaRochaRomain   File: CovMat.py    GNU General Public License v3.0 5 votes vote down vote up
def determinant(self):
        if self.__determinant is not None:
            return self.__determinant

        self.__determinant = numpy.linalg.det(self._numpy_array)
        return self.__determinant 
Example 27
Project: hoggorm   Author: olivertomic   File: statTools.py    BSD 2-Clause "Simplified" License 5 votes vote down vote up
def ortho(arr1, arr2):
    """
    This function orthogonalises arr1 with respect to arr2. The function then
    returns orthogonalised array arr1_orth.

    PARAMETERS
    ----------
    arr1 : numpy array
        A numpy array containing some data

    arr2 : numpy array
        A numpy array containing some data

    RETURNS
    -------
    numpy array
        A numpy array holding orthogonalised numpy array ``arr1``.

    Examples
    --------
    some examples

    """

    # Find number of rows, such that identity matrix I can be created
    numberRows = numpy.shape(arr1)[0]
    I = numpy.identity(numberRows, float)

    # Compute transpose of arr1
    arr2_T = numpy.transpose(arr2)

    term1 = numpy.linalg.inv(numpy.dot(arr2_T, arr2))
    term2 = numpy.dot(arr2, term1)
    term3 = numpy.dot(term2, arr2_T)
    arr1_orth = numpy.dot((I - term3), arr1)

    return arr1_orth 
Example 28
Project: hoggorm   Author: olivertomic   File: statTools.py    BSD 2-Clause "Simplified" License 5 votes vote down vote up
def matrixRank(arr, tol=1e-8):
    """
    Computes the rank of an array/matrix, i.e. number of linearly independent
    variables. This is not the same as numpy.rank() which only returns the
    number of ways (2-way, 3-way, etc) an array/matrix has.

    PARAMETERS
    ----------
    arrX : numpy array
        A numpy array containing the data

    RETURNS
    -------
    scalar
        Rank of matrix.

    Examples
    --------
    >>> import hoggorm as ho
    >>>
    >>> # Get the rank of the data
    >>> ho.matrixRank(myData)
    >>> 8

    """
    if len(arr.shape) != 2:
        raise ValueError('Input must be a 2-d array or Matrix object')

    s = numpy.linalg.svd(arr, compute_uv=0)
    return numpy.sum(numpy.where(s > tol, 1, 0)) 
Example 29
Project: MARRtino-2.0   Author: DaniAffCH   File: test_public_api.py    GNU General Public License v3.0 5 votes vote down vote up
def test_numpy_linalg():
    bad_results = check_dir(np.linalg)
    assert bad_results == {} 
Example 30
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    MIT License 5 votes vote down vote up
def test_eigvalsh_grad():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the geigvalsh op.")
    import scipy.linalg

    rng = numpy.random.RandomState(utt.fetch_seed())
    a = rng.randn(5, 5)
    a = a + a.T
    b = 10 * numpy.eye(5, 5) + rng.randn(5, 5)
    tensor.verify_grad(lambda a, b: eigvalsh(a, b).dot([1, 2, 3, 4, 5]),
                       [a, b], rng=numpy.random) 
Example 31
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    MIT License 5 votes vote down vote up
def test_solve_correctness(self):
        if not imported_scipy:
            raise SkipTest("Scipy needed for the Cholesky op.")
        rng = numpy.random.RandomState(utt.fetch_seed())
        A = theano.tensor.matrix()
        b = theano.tensor.matrix()
        y = self.op(A, b)
        gen_solve_func = theano.function([A, b], y)

        cholesky_lower = Cholesky(lower=True)
        L = cholesky_lower(A)
        y_lower = self.op(L, b)
        lower_solve_func = theano.function([L, b], y_lower)

        cholesky_upper = Cholesky(lower=False)
        U = cholesky_upper(A)
        y_upper = self.op(U, b)
        upper_solve_func = theano.function([U, b], y_upper)

        b_val = numpy.asarray(rng.rand(5, 1), dtype=config.floatX)
        
        # 1-test general case
        A_val = numpy.asarray(rng.rand(5, 5), dtype=config.floatX)
        # positive definite matrix:
        A_val = numpy.dot(A_val.transpose(), A_val)
        assert numpy.allclose(scipy.linalg.solve(A_val, b_val),
                              gen_solve_func(A_val, b_val))

        # 2-test lower traingular case
        L_val = scipy.linalg.cholesky(A_val, lower=True)
        assert numpy.allclose(scipy.linalg.solve_triangular(L_val, b_val, lower=True),
                              lower_solve_func(L_val, b_val))

        # 3-test upper traingular case
        U_val = scipy.linalg.cholesky(A_val, lower=False)
        assert numpy.allclose(scipy.linalg.solve_triangular(U_val, b_val, lower=False),
                              upper_solve_func(U_val, b_val)) 
Example 32
Project: D-VAE   Author: muhanzhang   File: test_slinalg.py    MIT License 5 votes vote down vote up
def test_expm():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the expm op.")
    rng = numpy.random.RandomState(utt.fetch_seed())
    A = rng.randn(5, 5).astype(config.floatX)

    ref = scipy.linalg.expm(A)

    x = tensor.matrix()
    m = expm(x)
    expm_f = function([x], m)

    val = expm_f(A)
    numpy.testing.assert_array_almost_equal(val, ref) 
Example 33
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    MIT License 5 votes vote down vote up
def test_qr_modes():
    rng = numpy.random.RandomState(utt.fetch_seed())

    A = tensor.matrix("A", dtype=theano.config.floatX)
    a = rng.rand(4, 4).astype(theano.config.floatX)

    f = function([A], qr(A))
    t_qr = f(a)
    n_qr = numpy.linalg.qr(a)
    assert _allclose(n_qr, t_qr)

    for mode in ["reduced", "r", "raw", "full", "economic"]:
        f = function([A], qr(A, mode))
        t_qr = f(a)
        n_qr = numpy.linalg.qr(a, mode)
        if isinstance(n_qr, (list, tuple)):
            assert _allclose(n_qr[0], t_qr[0])
            assert _allclose(n_qr[1], t_qr[1])
        else:
            assert _allclose(n_qr, t_qr)

    try:
        n_qr = numpy.linalg.qr(a, "complete")
        f = function([A], qr(A, "complete"))
        t_qr = f(a)
        assert _allclose(n_qr, t_qr)
    except TypeError as e:
        assert "name 'complete' is not defined" in str(e) 
Example 34
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    MIT License 5 votes vote down vote up
def test_svd():
    rng = numpy.random.RandomState(utt.fetch_seed())
    A = tensor.matrix("A", dtype=theano.config.floatX)
    U, V, T = svd(A)
    fn = function([A], [U, V, T])
    a = rng.rand(4, 4).astype(theano.config.floatX)
    n_u, n_v, n_t = numpy.linalg.svd(a)
    t_u, t_v, t_t = fn(a)

    assert _allclose(n_u, t_u)
    assert _allclose(n_v, t_v)
    assert _allclose(n_t, t_t) 
Example 35
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    MIT License 5 votes vote down vote up
def test_inverse_singular():
    singular = numpy.array([[1, 0, 0]] + [[0, 1, 0]] * 2,
                           dtype=theano.config.floatX)
    a = tensor.matrix()
    f = function([a], matrix_inverse(a))
    try:
        f(singular)
    except numpy.linalg.LinAlgError:
        return
    assert False 
Example 36
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    MIT License 5 votes vote down vote up
def test_det():
    rng = numpy.random.RandomState(utt.fetch_seed())

    r = rng.randn(5, 5).astype(config.floatX)
    x = tensor.matrix()
    f = theano.function([x], det(x))
    assert numpy.allclose(numpy.linalg.det(r), f(r)) 
Example 37
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    MIT License 5 votes vote down vote up
def test_wrong_coefficient_matrix(self):
        x = tensor.vector()
        y = tensor.vector()
        z = tensor.scalar()
        b = theano.tensor.nlinalg.lstsq()(x, y, z)
        f = function([x, y, z], b)
        self.assertRaises(numpy.linalg.linalg.LinAlgError, f, [2, 1], [2, 1], 1) 
Example 38
Project: D-VAE   Author: muhanzhang   File: test_nlinalg.py    MIT License 5 votes vote down vote up
def test_numpy_compare(self):
        rng = numpy.random.RandomState(utt.fetch_seed())
        A = tensor.matrix("A", dtype=theano.config.floatX)
        Q = matrix_power(A, 3)
        fn = function([A], [Q])
        a = rng.rand(4, 4).astype(theano.config.floatX)

        n_p = numpy.linalg.matrix_power(a, 3)
        t_p = fn(a)
        assert numpy.allclose(n_p, t_p) 
Example 39
Project: LEAR   Author: nmrksic   File: lear.py    Apache License 2.0 5 votes vote down vote up
def asymmetric_distance(v1, v2, distance_metric, order):
    """
    NOTE: this must be changed whenever cost function is reoriented.
    """
    #return distance(v1, v2) + norm(v1) - norm(v2) 
    
    cosine_similarity = ( dot(v1, v2) / ( norm(v1) * norm(v2) )) 

    norm1 = scipy.linalg.norm(v1, ord=order)
    norm2 = scipy.linalg.norm(v2, ord=order)

    if distance_metric == "metric_1":
        # |x| - |y|
        return (1-cosine_similarity) + (norm1 - norm2)

    elif distance_metric == "metric_2":
        # (|x| - |y|) / (|x| + |y|)
        
        norm_difference = norm1 - norm2
        norm_sum = norm1 + norm2

        return (1-cosine_similarity) + (norm_difference / norm_sum)

    elif distance_metric == "metric_3":

        max_norm = numpy.maximum(norm1, norm2)
        norm_difference = norm1 - norm2

        return (1-cosine_similarity) + (norm_difference / max_norm) 
Example 40
Project: vnpy_crypto   Author: birforce   File: linalg.py    MIT License 5 votes vote down vote up
def pinv(a, cond=None, rcond=None):
    """Compute the (Moore-Penrose) pseudo-inverse of a matrix.

    Calculate a generalized inverse of a matrix using a least-squares
    solver.

    Parameters
    ----------
    a : array, shape (M, N)
        Matrix to be pseudo-inverted
    cond, rcond : float
        Cutoff for 'small' singular values in the least-squares solver.
        Singular values smaller than rcond*largest_singular_value are
        considered zero.

    Returns
    -------
    B : array, shape (N, M)

    Raises LinAlgError if computation does not converge

    Examples
    --------
    >>> from numpy import *
    >>> a = random.randn(9, 6)
    >>> B = linalg.pinv(a)
    >>> allclose(a, dot(a, dot(B, a)))
    True
    >>> allclose(B, dot(B, dot(a, B)))
    True

    """
    a = asarray_chkfinite(a)
    b = numpy.identity(a.shape[0], dtype=a.dtype)
    if rcond is not None:
        cond = rcond
    return lstsq(a, b, cond=cond)[0] 
Example 41
Project: SimpleSVClustering   Author: josiahw   File: SimpleSVC.py    BSD 2-Clause "Simplified" License 5 votes vote down vote up
def rbfKernel(a,b,gamma):
    return numpy.exp(-gamma * numpy.linalg.norm(a - b)) 
Example 42
Project: ble5-nrf52-mac   Author: tomasero   File: test_public_api.py    MIT License 5 votes vote down vote up
def test_numpy_linalg():
    bad_results = check_dir(np.linalg)
    assert bad_results == {} 
Example 43
Project: OpenFermion   Author: quantumlib   File: _davidson_test.py    Apache License 2.0 5 votes vote down vote up
def test_with_built_in(self):
        """Compare with eigenvalues from built-in functions."""
        eigen_values, _ = numpy.linalg.eig(self.matrix)
        eigen_values = sorted(eigen_values)
        self.assertTrue(numpy.allclose(eigen_values, self.eigen_values))

        # Checks for eigh() function.
        eigen_values, eigen_vectors = numpy.linalg.eigh(self.matrix)
        self.assertAlmostEqual(get_difference(self.davidson.linear_operator,
                                              eigen_values, eigen_vectors), 0) 
Example 44
Project: OpenFermion   Author: quantumlib   File: _sparse_tools.py    Apache License 2.0 5 votes vote down vote up
def sparse_eigenspectrum(sparse_operator):
    """Perform a dense diagonalization.

    Returns:
        eigenspectrum: The lowest eigenvalues in a numpy array.
    """
    dense_operator = sparse_operator.todense()
    if is_hermitian(sparse_operator):
        eigenspectrum = numpy.linalg.eigvalsh(dense_operator)
    else:
        eigenspectrum = numpy.linalg.eigvals(dense_operator)
    return numpy.sort(eigenspectrum) 
Example 45
Project: OpenFermion   Author: quantumlib   File: _davidson.py    Apache License 2.0 5 votes vote down vote up
def orthonormalize(vectors, num_orthonormals=1, eps=1e-6):
    """Orthonormalize vectors, so that they're all normalized and orthogoal.

    The first vector is the same to that of vectors, while vector_i is
    orthogonal to vector_j, where j < i.

    Args:
        vectors(numpy.ndarray(complex)): Input vectors to be
            orthonormalized.
        num_orthonormals(int): First `num_orthonormals` columns are already
            orthonormal, so that one doesn't need to make any changes.
        eps(float): criterion of elements' max absolute value for zero vectors.

    Returns:
        ortho_normals(numpy.ndarray(complex)): Output orthonormal vectors.
    """
    ortho_normals = vectors
    count_orthonormals = num_orthonormals
    # Skip unchanged ones.
    for i in range(num_orthonormals, vectors.shape[1]):
        vector_i = vectors[:, i]
        # Makes sure vector_i is orthogonal to all processed vectors.
        for j in range(i):
            vector_i -= ortho_normals[:, j] * numpy.dot(
                ortho_normals[:, j].conj(), vector_i)

        # Makes sure vector_i is normalized.
        if numpy.max(numpy.abs(vector_i)) < eps:
            continue
        ortho_normals[:, count_orthonormals] = (vector_i /
                                                numpy.linalg.norm(vector_i))
        count_orthonormals += 1
    return ortho_normals[:, :count_orthonormals] 
Example 46
Project: fragalysis   Author: xchem   File: vector_utils.py    Apache License 2.0 5 votes vote down vote up
def get_best_fit_plane(pts, weights=None):
    """

    :param pts:
    :param weights:
    :return:
    """
    linalg = np.linalg
    if weights is None:
        wSum = len(pts)
        origin = np.sum(pts, 0)
    origin /= wSum
    sums = np.zeros((3, 3), np.double)
    for pt in pts:
        dp = pt - origin
        for i in range(3):
            sums[i, i] += dp[i] * dp[i]
            for j in range(i + 1, 3):
                sums[i, j] += dp[i] * dp[j]
                sums[j, i] += dp[i] * dp[j]
    sums /= wSum
    vals, vects = linalg.eigh(sums)
    order = np.argsort(vals)
    normal = vects[:, order[0]]
    plane = np.zeros((4,), np.double)
    plane[:3] = normal
    plane[3] = -1 * normal.dot(origin)
    return plane 
Example 47
Project: poker   Author: surgebiswas   File: test_slinalg.py    MIT License 5 votes vote down vote up
def test_eigvalsh_grad():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the geigvalsh op.")
    import scipy.linalg

    rng = numpy.random.RandomState(utt.fetch_seed())
    a = rng.randn(5, 5)
    a = a + a.T
    b = 10 * numpy.eye(5, 5) + rng.randn(5, 5)
    tensor.verify_grad(lambda a, b: eigvalsh(a, b).dot([1, 2, 3, 4, 5]),
                       [a, b], rng=numpy.random) 
Example 48
Project: poker   Author: surgebiswas   File: test_slinalg.py    MIT License 5 votes vote down vote up
def test_solve_correctness(self):
        if not imported_scipy:
            raise SkipTest("Scipy needed for the Cholesky op.")
        rng = numpy.random.RandomState(utt.fetch_seed())
        A = theano.tensor.matrix()
        b = theano.tensor.matrix()
        y = self.op(A, b)
        gen_solve_func = theano.function([A, b], y)

        cholesky_lower = Cholesky(lower=True)
        L = cholesky_lower(A)
        y_lower = self.op(L, b)
        lower_solve_func = theano.function([L, b], y_lower)

        cholesky_upper = Cholesky(lower=False)
        U = cholesky_upper(A)
        y_upper = self.op(U, b)
        upper_solve_func = theano.function([U, b], y_upper)

        b_val = numpy.asarray(rng.rand(5, 1), dtype=config.floatX)
        
        # 1-test general case
        A_val = numpy.asarray(rng.rand(5, 5), dtype=config.floatX)
        # positive definite matrix:
        A_val = numpy.dot(A_val.transpose(), A_val)
        assert numpy.allclose(scipy.linalg.solve(A_val, b_val),
                              gen_solve_func(A_val, b_val))

        # 2-test lower traingular case
        L_val = scipy.linalg.cholesky(A_val, lower=True)
        assert numpy.allclose(scipy.linalg.solve_triangular(L_val, b_val, lower=True),
                              lower_solve_func(L_val, b_val))

        # 3-test upper traingular case
        U_val = scipy.linalg.cholesky(A_val, lower=False)
        assert numpy.allclose(scipy.linalg.solve_triangular(U_val, b_val, lower=False),
                              upper_solve_func(U_val, b_val)) 
Example 49
Project: poker   Author: surgebiswas   File: test_slinalg.py    MIT License 5 votes vote down vote up
def test_expm():
    if not imported_scipy:
        raise SkipTest("Scipy needed for the expm op.")
    rng = numpy.random.RandomState(utt.fetch_seed())
    A = rng.randn(5, 5).astype(config.floatX)

    ref = scipy.linalg.expm(A)

    x = tensor.matrix()
    m = expm(x)
    expm_f = function([x], m)

    val = expm_f(A)
    numpy.testing.assert_array_almost_equal(val, ref) 
Example 50
Project: poker   Author: surgebiswas   File: test_nlinalg.py    MIT License 5 votes vote down vote up
def test_qr_modes():
    rng = numpy.random.RandomState(utt.fetch_seed())

    A = tensor.matrix("A", dtype=theano.config.floatX)
    a = rng.rand(4, 4).astype(theano.config.floatX)

    f = function([A], qr(A))
    t_qr = f(a)
    n_qr = numpy.linalg.qr(a)
    assert _allclose(n_qr, t_qr)

    for mode in ["reduced", "r", "raw", "full", "economic"]:
        f = function([A], qr(A, mode))
        t_qr = f(a)
        n_qr = numpy.linalg.qr(a, mode)
        if isinstance(n_qr, (list, tuple)):
            assert _allclose(n_qr[0], t_qr[0])
            assert _allclose(n_qr[1], t_qr[1])
        else:
            assert _allclose(n_qr, t_qr)

    try:
        n_qr = numpy.linalg.qr(a, "complete")
        f = function([A], qr(A, "complete"))
        t_qr = f(a)
        assert _allclose(n_qr, t_qr)
    except TypeError as e:
        assert "name 'complete' is not defined" in str(e)