Python numpy.linalg() Examples

def rotation_matrix_axis_and_angle_2(R, debug=False, errorThreshold=10**-7, msg=None):
    w, v = numpy.linalg.eig(R) #this method is not used at the primary method as numpy.linalg.eig does not return answers in high enough precision
    angle, axis = None, None
    eigErrors = abs(w -1) #errors from 1+0j
    i = (eigErrors == min(eigErrors)).tolist().index(True)
    axis = numpy.real(v[:,i])
    if i != 1:
        angle = arccos2(  numpy.real( w[1] ) )
    else:
        angle = arccos2(  numpy.real( w[0] ) )
    error  = norm(axis_rotation_matrix(angle, *axis) - R)
    if debug: print('rotation_matrix_axis_and_angle error %1.1e' % error)
    if error > errorThreshold:
        angle = -angle
        error = norm(axis_rotation_matrix(angle, *axis) - R)
        if error > errorThreshold:
            R_pickle_str = pickle.dumps(R)
            #R_abs_minus_identity = abs(R) - numpy.eye(3)
            print(R*R.transpose())
            raise ValueError( 'rotation_matrix_axis_and_angle_2: no solution found! locals %s' % str(locals()))
    return axis, angle 

def eig(h, s):
    from pyscf import symm
    nirrep = len(mol.symm_orb)
    h = symm.symmetrize_matrix(h, mol.symm_orb)
    s = symm.symmetrize_matrix(s, mol.symm_orb)
    cs = []
    es = []
#
# Linear dependency are removed by looping over different symmetry irreps.
#
    for ir in range(nirrep):
        d, t = numpy.linalg.eigh(s[ir])
        x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8])
        xhx = reduce(numpy.dot, (x.T, h[ir], x))
        e, c = numpy.linalg.eigh(xhx)
        cs.append(reduce(numpy.dot, (mol.symm_orb[ir], x, c)))
        es.append(e)
    e = numpy.hstack(es)
    c = numpy.hstack(cs)
    return e, c 

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 

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) 

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) 

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) 

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())) 

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,
        ]) 

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) 

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) 

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 

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 

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) 

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) 

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) 

def xavier_vector(word, D=300):
    """
    Returns a D-dimensional vector for the word. 

    We hash the word to always get the same vector for the given word. 
    """

    seed_value = hash_string(word)
    numpy.random.seed(seed_value)

    neg_value = - math.sqrt(6)/math.sqrt(D)
    pos_value = math.sqrt(6)/math.sqrt(D)

    rsample = numpy.random.uniform(low=neg_value, high=pos_value, size=(D,))
    norm = numpy.linalg.norm(rsample)
    rsample_normed = rsample/norm

    return rsample_normed 

def calculateNormalsFromIndexedVertices(vertices: numpy.ndarray, indices: numpy.ndarray, face_count: int) -> numpy.ndarray:
    """Calculate the normals of this mesh of triagles using indexes.

    :param vertices: :type{narray} list of vertices as a 1D list of float triples
    :param indices: :type{narray} list of indices as a 1D list of integers
    :param face_count: :type{integer} the number of triangles defined by the indices array
    :return: :type{narray} list normals as a 1D array of floats, each group of 3 floats is a vector
    """

    start_time = time()
    # Numpy magic!
    # First, reset the normals
    normals = numpy.zeros((face_count*3, 3), dtype=numpy.float32)

    for face in indices[0:face_count]:
        normals[face[0]] = numpy.cross(vertices[face[0]] - vertices[face[1]], vertices[face[0]] - vertices[face[2]])
        length = numpy.linalg.norm(normals[face[0]])
        normals[face[0]] /= length
        normals[face[1]] = normals[face[0]]
        normals[face[2]] = normals[face[0]]
    end_time = time()
    Logger.log("d", "Calculating normals took %s seconds", end_time - start_time)
    return normals 

def eig(h, s):
    from pyscf import symm
    nirrep = len(mol.symm_orb)
    h = symm.symmetrize_matrix(h, mol.symm_orb)
    s = symm.symmetrize_matrix(s, mol.symm_orb)
    cs = []
    es = []
#
# Linear dependency are removed by looping over different symmetry irreps.
#
    for ir in range(nirrep):
        d, t = numpy.linalg.eigh(s[ir])
        x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8])
        xhx = reduce(numpy.dot, (x.T, h[ir], x))
        e, c = numpy.linalg.eigh(xhx)
        cs.append(reduce(numpy.dot, (mol.symm_orb[ir], x, c)))
        es.append(e)
    e = numpy.hstack(es)
    c = numpy.hstack(cs)
    return e, c 

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 

def distance_between_axes_fmin( p1, u1, p2, u2):
    from scipy.optimize import fmin_bfgs
    def distance(T):
        t1, t2 = T
        return numpy.linalg.norm( p1 + u1*t1 - (p2 + u2*t2) )
    T_opt = fmin_bfgs( distance, [0 , 0], disp=False)
    return distance(T_opt) 

def test_numpy_linalg():
    bad_results = check_dir(np.linalg)
    assert bad_results == {} 

def closest_perp_distance_between_sample_points(NullcA, NullcB, xa, x0B, x1B,
                                                newt_tol=1e-6):
    """Closest perpendicular distance from (xa, ya) on Nullcline A to Nullcline B,
    given that it's known to be between sample points x0B and x1B on Nullcline B.

    Uses a Newton-like step, solving up to an angular tolerance given by newt_tol.
    """
    global plotter
    ya = NullcA(xa)
    a = np.array([xa,ya])
    tgt_vec_at_a = NullcA.tgt_vec(xa, rescale=True)
    # normal has length 1
    normal_at_a = make_vec_at_A_face_B(get_orthonormal(tgt_vec_at_a),
                                ya, NullcB(xa))
    #plotter.plot_line_from_points(a, a+tgt_vec_at_a, 'k:')
    phi = angle_to_vertical(normal_at_a)

    A = Point2D(a)
    B = Point2D(a + 0.2*normal_at_a)
    C = Point2D(x0B, NullcB(x0B))
    D = Point2D(x1B, NullcB(x1B))

    P0 = line_intersection(A, B, C, D)
    Q0 = Point2D(P0.x, NullcB(P0.x))  # project onto spline
    #theta0 = angle_to_vertical(Q0-a)

    #plotter.plot_line_from_points(A, B, 'k-')
    #plotter.plot_line_from_points(C, D, 'r--')
    #plotter.plot_point(P0, 'gs')
    #plotter.plot_point(Q0, 'ys')

    try:
        Q = _newton_step(Q0, NullcB, A, B, phi, newt_tol)
    except ValueError:
        return np.Inf
    else:
        vec = Q-A
        #vec[0] = vec[0] / 100.  # hack to investigate differing x and y scales
        return np.linalg.norm(vec) 

def angle_to_vertical(v):
    """Return an angle between 0 and 2*pi measured clockwise from vertical."""
    up = np.array((0.,1.))
    theta = np.arccos(np.dot(up, v)/np.linalg.norm(v))
    if np.cross(up, v) > 0:
        # angle will be greater than pi
        return 2*np.pi - theta
    else:
        return theta


# not currently used, but helpful if lines are not arranged correctly for intersection test 

def test_numpy_linalg():
    bad_results = check_dir(np.linalg)
    assert bad_results == {} 

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) 

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) 

def eig(h, s):
    d, t = numpy.linalg.eigh(s)
# Removing the eigenvectors assoicated to the smallest eigenvalue, the new
# basis defined by x matrix has 139 vectors.
    x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8])
    xhx = reduce(numpy.dot, (x.T, h, x))
    e, c = numpy.linalg.eigh(xhx)
    c = numpy.dot(x, c)
# Return 139 eigenvalues and 139 eigenvectors.
    return e, c
#
# Replacing the default eig function with the above one,  the HF solver
# generate only 139 canonical orbitals
# 

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 

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)) 

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)