Python numpy.outer() Examples

The following are 30 code examples for showing how to use numpy.outer(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: discomll   Author: romanorac   File: linear_svm.py    License: Apache License 2.0 6 votes vote down vote up
def map_fit(interface, state, label, inp):
    """
    Function calculates matrices ete and etde for every sample, aggregates and output them.
    """
    import numpy as np
    ete, etde = 0, 0
    out = interface.output(0)

    for row in inp:
        row = row.strip().split(state["delimiter"])  # split row
        if len(row) > 1:  # check if row is empty
            # intercept term is added to every sample
            x = np.array([(0 if v in state["missing_vals"] else float(v)) for i, v in enumerate(row) if
                          i in state["X_indices"]] + [-1])
            # map label value to 1 or -1. If label does not match set error
            y = 1 if state["y_map"][0] == row[state["y_index"]] else -1 if state["y_map"][1] == row[
                state["y_index"]] else "Error"
            ete += np.outer(x, x)
            etde += x * y
    out.add("etde", etde)
    for i, row in enumerate(ete):
        out.add(i, row) 
Example 2
Project: discomll   Author: romanorac   File: locally_weighted_linear_regression.py    License: Apache License 2.0 6 votes vote down vote up
def map_fit(interface, state, label, inp):
    import numpy as np
    combiner = dict([(k, [0, 0]) for k in state["samples"].keys()])
    out = interface.output(0)

    for row in inp:
        row = row.strip().split(state["delimiter"])
        if len(row) > 1:
            x = np.array([1] + [(0 if v in state["missing_vals"] else float(v)) for i, v in enumerate(row) if
                                i in state["X_indices"]])
            y = float(row[state["y_index"]])

            for test_id, x1 in state["samples"].iteritems():
                w = np.exp(np.dot(-(x1 - x).transpose(), x1 - x) / state["tau"])
                combiner[test_id][0] += w * np.outer(x, x)
                combiner[test_id][1] += w * x * y

    for k, v in combiner.iteritems():
        out.add(k + state["delimiter"] + "b", v[1])
        for i in range(len(v[0])):
            out.add(k + state["delimiter"] + "A" + state["delimiter"] + str(i), v[0][i]) 
Example 3
Project: discomll   Author: romanorac   File: linear_regression.py    License: Apache License 2.0 6 votes vote down vote up
def map_fit(interface, state, label, inp):
    import numpy as np
    a, b = 0, 0
    out = interface.output(0)

    for row in inp:
        row = row.strip().split(state["delimiter"])
        if len(row) > 1:
            x = np.array([1] + [(0 if v in state["missing_vals"] else float(v)) for i, v in enumerate(row) if
                                i in state["X_indices"]])
            y = float(row[state["y_index"]])
            a += np.outer(x, x)
            b += x * y
    out.add("b", b)
    for i, row in enumerate(a):
        out.add(i, row) 
Example 4
Project: OpenFermion-Cirq   Author: quantumlib   File: fermionic_simulation.py    License: Apache License 2.0 6 votes vote down vote up
def _eigen_components(self):
        components = [(0, np.diag([1, 1, 1, 0, 1, 0, 0, 1]))]
        nontrivial_part = np.zeros((3, 3), dtype=np.complex128)
        for ij, w in zip([(1, 2), (0, 2), (0, 1)], self.weights):
            nontrivial_part[ij] = w
            nontrivial_part[ij[::-1]] = w.conjugate()
        assert np.allclose(nontrivial_part, nontrivial_part.conj().T)
        eig_vals, eig_vecs = np.linalg.eigh(nontrivial_part)
        for eig_val, eig_vec in zip(eig_vals, eig_vecs.T):
            exp_factor = -eig_val / np.pi
            proj = np.zeros((8, 8), dtype=np.complex128)
            nontrivial_indices = np.array([3, 5, 6], dtype=np.intp)
            proj[nontrivial_indices[:, np.newaxis], nontrivial_indices] = (
                np.outer(eig_vec.conjugate(), eig_vec))
            components.append((exp_factor, proj))
        return components 
Example 5
Project: robosuite   Author: StanfordVL   File: transform_utils.py    License: MIT License 6 votes vote down vote up
def quat2mat(quaternion):
    """
    Converts given quaternion (x, y, z, w) to matrix.

    Args:
        quaternion: vec4 float angles

    Returns:
        3x3 rotation matrix
    """
    q = np.array(quaternion, dtype=np.float32, copy=True)[[3, 0, 1, 2]]
    n = np.dot(q, q)
    if n < EPS:
        return np.identity(3)
    q *= math.sqrt(2.0 / n)
    q = np.outer(q, q)
    return np.array(
        [
            [1.0 - q[2, 2] - q[3, 3], q[1, 2] - q[3, 0], q[1, 3] + q[2, 0]],
            [q[1, 2] + q[3, 0], 1.0 - q[1, 1] - q[3, 3], q[2, 3] - q[1, 0]],
            [q[1, 3] - q[2, 0], q[2, 3] + q[1, 0], 1.0 - q[1, 1] - q[2, 2]],
        ]
    ) 
Example 6
Project: scarlet   Author: pmelchior   File: interpolation.py    License: MIT License 6 votes vote down vote up
def quintic_spline(dx, dtype=np.float64):
    def inner(x):
        return 1 + x ** 3 / 12 * (-95 + 138 * x - 55 * x ** 2)

    def middle(x):
        return (x - 1) * (x - 2) / 24 * (-138 + 348 * x - 249 * x ** 2 + 55 * x ** 3)

    def outer(x):
        return (x - 2) * (x - 3) ** 2 / 24 * (-54 + 50 * x - 11 * x ** 2)

    window = np.arange(-3, 4)
    x = np.abs(dx - window)
    result = np.piecewise(
        x,
        [x <= 1, (x > 1) & (x <= 2), (x > 2) & (x <= 3)],
        [lambda x: inner(x), lambda x: middle(x), lambda x: outer(x)],
    )
    return result, window 
Example 7
Project: pyscf   Author: pyscf   File: test_common.py    License: Apache License 2.0 6 votes vote down vote up
def tdhf_frozen_mask(eri, kind="ov"):
    if isinstance(eri.nocc, int):
        nocc = int(eri.model.mo_occ.sum() // 2)
        mask = eri.space
    else:
        nocc = numpy.array(tuple(int(i.sum() // 2) for i in eri.model.mo_occ))
        assert numpy.all(nocc == nocc[0])
        assert numpy.all(eri.space == eri.space[0, numpy.newaxis, :])
        nocc = nocc[0]
        mask = eri.space[0]
    mask_o = mask[:nocc]
    mask_v = mask[nocc:]
    if kind == "ov":
        mask_ov = numpy.outer(mask_o, mask_v).reshape(-1)
        return numpy.tile(mask_ov, 2)
    elif kind == "1ov":
        return numpy.outer(mask_o, mask_v).reshape(-1)
    elif kind == "sov":
        mask_ov = numpy.outer(mask_o, mask_v).reshape(-1)
        nk = len(eri.model.mo_occ)
        return numpy.tile(mask_ov, 2 * nk ** 2)
    elif kind == "o,v":
        return mask_o, mask_v 
Example 8
Project: pyscf   Author: pyscf   File: test_common.py    License: Apache License 2.0 6 votes vote down vote up
def tdhf_frozen_mask(eri, kind="ov"):
    if isinstance(eri.nocc, int):
        nocc = int(eri.model.mo_occ.sum() // 2)
        mask = eri.space
    else:
        nocc = numpy.array(tuple(int(i.sum() // 2) for i in eri.model.mo_occ))
        assert numpy.all(nocc == nocc[0])
        assert numpy.all(eri.space == eri.space[0, numpy.newaxis, :])
        nocc = nocc[0]
        mask = eri.space[0]
    mask_o = mask[:nocc]
    mask_v = mask[nocc:]
    if kind == "ov":
        mask_ov = numpy.outer(mask_o, mask_v).reshape(-1)
        return numpy.tile(mask_ov, 2)
    elif kind == "1ov":
        return numpy.outer(mask_o, mask_v).reshape(-1)
    elif kind == "sov":
        mask_ov = numpy.outer(mask_o, mask_v).reshape(-1)
        nk = len(eri.model.mo_occ)
        return numpy.tile(mask_ov, 2 * nk ** 2)
    elif kind == "o,v":
        return mask_o, mask_v 
Example 9
Project: pyscf   Author: pyscf   File: test_common.py    License: Apache License 2.0 6 votes vote down vote up
def tdhf_frozen_mask(eri, kind="ov"):
    if isinstance(eri.nocc, int):
        nocc = int(eri.model.mo_occ.sum() // 2)
        mask = eri.space
    else:
        nocc = numpy.array(tuple(int(i.sum() // 2) for i in eri.model.mo_occ))
        assert numpy.all(nocc == nocc[0])
        assert numpy.all(eri.space == eri.space[0, numpy.newaxis, :])
        nocc = nocc[0]
        mask = eri.space[0]
    mask_o = mask[:nocc]
    mask_v = mask[nocc:]
    if kind == "ov":
        mask_ov = numpy.outer(mask_o, mask_v).reshape(-1)
        return numpy.tile(mask_ov, 2)
    elif kind == "1ov":
        return numpy.outer(mask_o, mask_v).reshape(-1)
    elif kind == "sov":
        mask_ov = numpy.outer(mask_o, mask_v).reshape(-1)
        nk = len(eri.model.mo_occ)
        return numpy.tile(mask_ov, 2 * nk ** 2)
    elif kind == "o,v":
        return mask_o, mask_v 
Example 10
Project: me-ica   Author: ME-ICA   File: test_dwiparams.py    License: GNU Lesser General Public License v2.1 6 votes vote down vote up
def test_b2q():
    # conversion of b matrix to q
    q = np.array([1,2,3])
    s = np.sqrt(np.sum(q * q)) # vector norm
    B = np.outer(q, q)
    assert_array_almost_equal(q*s, B2q(B))
    q = np.array([1,2,3])
    # check that the sign of the vector as positive x convention
    B = np.outer(-q, -q)
    assert_array_almost_equal(q*s, B2q(B))
    q = np.array([-1, 2, 3])
    B = np.outer(q, q)
    assert_array_almost_equal(-q*s, B2q(B))
    B = np.eye(3) * -1
    assert_raises(ValueError, B2q, B)
    # no error if we up the tolerance
    q = B2q(B, tol=1) 
Example 11
Project: simnibs   Author: simnibs   File: test_optimization_methods.py    License: GNU General Public License v3.0 6 votes vote down vote up
def optimization_variables_avg_QCQP():
    l = np.array([4, 1, 1, -5, 3], dtype=float)
    A = np.array([[1, 2, 3, 4],
                  [3, 1, 2, 2],
                  [2, 3, 1, 1],
                  [1, 3, 5, 1]], dtype=float)
    A = np.vstack([-np.mean(A, axis=1), A - np.mean(A, axis=1)])
    Q = A.dot(A.T)
    A_in = np.array([[1, 5, 2, 4],
                     [5, 2, 1, 2],
                     [0, 1, 1, 1],
                     [2, 4, 3, 5]], dtype=float)
    A_in = np.vstack([-np.mean(A_in, axis=1), A_in - np.mean(A_in, axis=1)])
    Q_in = A_in.dot(A_in.T)
    Q_in += np.outer(l, l)
    return l, Q, np.ones(5), Q_in 
Example 12
Project: simnibs   Author: simnibs   File: test_optimization_methods.py    License: GNU General Public License v3.0 6 votes vote down vote up
def test_both_limit_angle_Q_iteration(self, optimization_variables_avg_QCQP):
        l, Q, A, Q_in = optimization_variables_avg_QCQP
        # l, Q, A = optimization_variables_avg
        # Q_in = 6 * np.eye(len(l)) + np.outer(l, l)
        max_angle = 20
        m = 2e-3
        m1 = 4e-3
        f = .01
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle)
        x_sp = optimize_focality(
            l, Q, f, max_el_current=m, max_total_current=m1,
            Qin=Q_in, max_angle=max_angle)
        assert np.linalg.norm(x, 1) <= 2 * m1 + 1e-4
        assert np.all(np.abs(x) <= m + 1e-4)
        assert np.isclose(np.sum(x), 0)
        assert np.isclose(l.dot(x), f)
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle)
        assert np.allclose(x.dot(Q.dot(x)), x_sp.dot(Q.dot(x_sp)), rtol=1e-4, atol=1e-5) 
Example 13
Project: simnibs   Author: simnibs   File: misc.py    License: GNU General Public License v3.0 6 votes vote down vote up
def mutcoh(A):
    """ mutcoh calculates the mutual coherence of a matrix A, i.e. the cosine
        of the smallest angle between two columns
      
        mutual_coherence = mutcoh(A)
 
        Input:
            A ... m x n matrix
 
        Output:
            mutual_coherence """

    T = np.dot(A.conj().T,A)
    s = np.sqrt(np.diag(T))    
    S = np.diag(s)
    mutual_coherence = np.max(1.0*(T-S)/np.outer(s,s))
    
    return mutual_coherence 
Example 14
Project: rai-python   Author: MarcToussaint   File: transformations.py    License: MIT License 6 votes vote down vote up
def reflection_matrix(point, normal):
    """Return matrix to mirror at plane defined by point and normal vector.
    >>> v0 = numpy.random.random(4) - 0.5
    >>> v0[3] = 1.0
    >>> v1 = numpy.random.random(3) - 0.5
    >>> R = reflection_matrix(v0, v1)
    >>> numpy.allclose(2., numpy.trace(R))
    True
    >>> numpy.allclose(v0, numpy.dot(R, v0))
    True
    >>> v2 = v0.copy()
    >>> v2[:3] += v1
    >>> v3 = v0.copy()
    >>> v2[:3] -= v1
    >>> numpy.allclose(v2, numpy.dot(R, v3))
    True
    """
    normal = unit_vector(normal[:3])
    M = numpy.identity(4)
    M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
    M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
    return M 
Example 15
Project: rai-python   Author: MarcToussaint   File: transformations.py    License: MIT License 6 votes vote down vote up
def shear_matrix(angle, direction, point, normal):
    """Return matrix to shear by angle along direction vector on shear plane.
    The shear plane is defined by a point and normal vector. The direction
    vector must be orthogonal to the plane's normal vector.
    A point P is transformed by the shear matrix into P" such that
    the vector P-P" is parallel to the direction vector and its extent is
    given by the angle of P-P'-P", where P' is the orthogonal projection
    of P onto the shear plane.
    >>> angle = (random.random() - 0.5) * 4*math.pi
    >>> direct = numpy.random.random(3) - 0.5
    >>> point = numpy.random.random(3) - 0.5
    >>> normal = numpy.cross(direct, numpy.random.random(3))
    >>> S = shear_matrix(angle, direct, point, normal)
    >>> numpy.allclose(1.0, numpy.linalg.det(S))
    True
    """
    normal = unit_vector(normal[:3])
    direction = unit_vector(direction[:3])
    if abs(numpy.dot(normal, direction)) > 1e-6:
        raise ValueError("direction and normal vectors are not orthogonal")
    angle = math.tan(angle)
    M = numpy.identity(4)
    M[:3, :3] += angle * numpy.outer(direction, normal)
    M[:3, 3] = -angle * numpy.dot(point[:3], normal) * direction
    return M 
Example 16
Project: recruit   Author: Frank-qlu   File: timer_comparison.py    License: Apache License 2.0 6 votes vote down vote up
def test_4(self):
        """
        Test of take, transpose, inner, outer products.

        """
        x = self.arange(24)
        y = np.arange(24)
        x[5:6] = self.masked
        x = x.reshape(2, 3, 4)
        y = y.reshape(2, 3, 4)
        assert self.allequal(np.transpose(y, (2, 0, 1)), self.transpose(x, (2, 0, 1)))
        assert self.allequal(np.take(y, (2, 0, 1), 1), self.take(x, (2, 0, 1), 1))
        assert self.allequal(np.inner(self.filled(x, 0), self.filled(y, 0)),
                            self.inner(x, y))
        assert self.allequal(np.outer(self.filled(x, 0), self.filled(y, 0)),
                            self.outer(x, y))
        y = self.array(['abc', 1, 'def', 2, 3], object)
        y[2] = self.masked
        t = self.take(y, [0, 3, 4])
        assert t[0] == 'abc'
        assert t[1] == 2
        assert t[2] == 3 
Example 17
Project: recruit   Author: Frank-qlu   File: test_old_ma.py    License: Apache License 2.0 6 votes vote down vote up
def test_testTakeTransposeInnerOuter(self):
        # Test of take, transpose, inner, outer products
        x = arange(24)
        y = np.arange(24)
        x[5:6] = masked
        x = x.reshape(2, 3, 4)
        y = y.reshape(2, 3, 4)
        assert_(eq(np.transpose(y, (2, 0, 1)), transpose(x, (2, 0, 1))))
        assert_(eq(np.take(y, (2, 0, 1), 1), take(x, (2, 0, 1), 1)))
        assert_(eq(np.inner(filled(x, 0), filled(y, 0)),
                   inner(x, y)))
        assert_(eq(np.outer(filled(x, 0), filled(y, 0)),
                   outer(x, y)))
        y = array(['abc', 1, 'def', 2, 3], object)
        y[2] = masked
        t = take(y, [0, 3, 4])
        assert_(t[0] == 'abc')
        assert_(t[1] == 2)
        assert_(t[2] == 3) 
Example 18
Project: recruit   Author: Frank-qlu   File: test_core.py    License: Apache License 2.0 6 votes vote down vote up
def test_TakeTransposeInnerOuter(self):
        # Test of take, transpose, inner, outer products
        x = arange(24)
        y = np.arange(24)
        x[5:6] = masked
        x = x.reshape(2, 3, 4)
        y = y.reshape(2, 3, 4)
        assert_equal(np.transpose(y, (2, 0, 1)), transpose(x, (2, 0, 1)))
        assert_equal(np.take(y, (2, 0, 1), 1), take(x, (2, 0, 1), 1))
        assert_equal(np.inner(filled(x, 0), filled(y, 0)),
                     inner(x, y))
        assert_equal(np.outer(filled(x, 0), filled(y, 0)),
                     outer(x, y))
        y = array(['abc', 1, 'def', 2, 3], object)
        y[2] = masked
        t = take(y, [0, 3, 4])
        assert_(t[0] == 'abc')
        assert_(t[1] == 2)
        assert_(t[2] == 3) 
Example 19
Project: ibllib   Author: int-brain-lab   File: _knockoff.py    License: MIT License 6 votes vote down vote up
def _get_knmat(exog, xcov, sl):
    # Utility function, see equation 2.2 of Barber & Candes.

    nobs, nvar = exog.shape

    ash = np.linalg.inv(xcov)
    ash *= -np.outer(sl, sl)
    i, j = np.diag_indices(nvar)
    ash[i, j] += 2 * sl

    umat = np.random.normal(size=(nobs, nvar))
    u, _ = np.linalg.qr(exog)
    umat -= np.dot(u, np.dot(u.T, umat))
    umat, _ = np.linalg.qr(umat)

    ashr, xc, _ = np.linalg.svd(ash, 0)
    ashr *= np.sqrt(xc)
    ashr = ashr.T

    ex = (sl[:, None] * np.linalg.solve(xcov, exog.T)).T
    exogn = exog - ex + np.dot(umat, ashr)

    return exogn 
Example 20
Project: cozmo_driver   Author: OTL   File: transformations.py    License: Apache License 2.0 6 votes vote down vote up
def reflection_matrix(point, normal):
    """Return matrix to mirror at plane defined by point and normal vector.

    >>> v0 = numpy.random.random(4) - 0.5
    >>> v0[3] = 1.0
    >>> v1 = numpy.random.random(3) - 0.5
    >>> R = reflection_matrix(v0, v1)
    >>> numpy.allclose(2., numpy.trace(R))
    True
    >>> numpy.allclose(v0, numpy.dot(R, v0))
    True
    >>> v2 = v0.copy()
    >>> v2[:3] += v1
    >>> v3 = v0.copy()
    >>> v2[:3] -= v1
    >>> numpy.allclose(v2, numpy.dot(R, v3))
    True

    """
    normal = unit_vector(normal[:3])
    M = numpy.identity(4)
    M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
    M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
    return M 
Example 21
Project: pylogit   Author: timothyb0912   File: choice_calcs.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def robust_outer_product(vec_1, vec_2):
    """
    Calculates a 'robust' outer product of two vectors that may or may not
    contain very small values.

    Parameters
    ----------
    vec_1 : 1D ndarray
    vec_2 : 1D ndarray

    Returns
    -------
    outer_prod : 2D ndarray. The outer product of vec_1 and vec_2
    """
    mantissa_1, exponents_1 = np.frexp(vec_1)
    mantissa_2, exponents_2 = np.frexp(vec_2)
    new_mantissas = mantissa_1[None, :] * mantissa_2[:, None]
    new_exponents = exponents_1[None, :] + exponents_2[:, None]
    return new_mantissas * np.exp2(new_exponents) 
Example 22
Project: pylogit   Author: timothyb0912   File: choice_calcs.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def calc_asymptotic_covariance(hessian, fisher_info_matrix):
    """
    Parameters
    ----------
    hessian : 2D ndarray.
        It should have shape `(num_vars, num_vars)`. It is the matrix of second
        derivatives of the total loss across the dataset, with respect to each
        pair of coefficients being estimated.
    fisher_info_matrix : 2D ndarray.
        It should have a shape of `(num_vars, num_vars)`.  It is the
        approximation of the negative of the expected hessian formed by taking
        the outer product of (each observation's gradient of the loss function)
        with itself, and then summing across all observations.

    Returns
    -------
    huber_white_matrix : 2D ndarray.
        Will have shape `(num_vars, num_vars)`. The entries in the returned
        matrix are calculated by the following formula:
        `hess_inverse * fisher_info_matrix * hess_inverse`.
    """
    # Calculate the inverse of the hessian
    hess_inv = scipy.linalg.inv(hessian)

    return np.dot(hess_inv, np.dot(fisher_info_matrix, hess_inv)) 
Example 23
Project: pylogit   Author: timothyb0912   File: test_choice_calcs.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_robust_outer_product(self):
        """
        Ensure that robust_outer_product returns the expected results.
        Unfortunately, I cannot find a good case now where using the regular
        outer product gives incorrect results. However without a compelling
        reason to remove the function, I'll trust my earlier judgement in
        creating it in the first place.
        """
        # Define a vector whose outer product we want to take
        x = np.array([1e-100, 0.01])
        outer_product = np.outer(x, x)
        robust_outer_product = cc.robust_outer_product(x, x)

        # Perform the desired tests
        self.assertIsInstance(robust_outer_product, np.ndarray)
        self.assertEqual(robust_outer_product.shape, outer_product.shape)
        npt.assert_allclose(outer_product, robust_outer_product)

        return None 
Example 24
Project: KAIR   Author: cszn   File: utils_image.py    License: MIT License 6 votes vote down vote up
def ssim(img1, img2):
    C1 = (0.01 * 255)**2
    C2 = (0.03 * 255)**2

    img1 = img1.astype(np.float64)
    img2 = img2.astype(np.float64)
    kernel = cv2.getGaussianKernel(11, 1.5)
    window = np.outer(kernel, kernel.transpose())

    mu1 = cv2.filter2D(img1, -1, window)[5:-5, 5:-5]  # valid
    mu2 = cv2.filter2D(img2, -1, window)[5:-5, 5:-5]
    mu1_sq = mu1**2
    mu2_sq = mu2**2
    mu1_mu2 = mu1 * mu2
    sigma1_sq = cv2.filter2D(img1**2, -1, window)[5:-5, 5:-5] - mu1_sq
    sigma2_sq = cv2.filter2D(img2**2, -1, window)[5:-5, 5:-5] - mu2_sq
    sigma12 = cv2.filter2D(img1 * img2, -1, window)[5:-5, 5:-5] - mu1_mu2

    ssim_map = ((2 * mu1_mu2 + C1) * (2 * sigma12 + C2)) / ((mu1_sq + mu2_sq + C1) *
                                                            (sigma1_sq + sigma2_sq + C2))
    return ssim_map.mean() 
Example 25
Project: bayeslite   Author: probcomp   File: stats.py    License: Apache License 2.0 6 votes vote down vote up
def chi2_contingency(contingency):
    """Pearson chi^2 test of independence statistic on contingency table.

    https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test#Test_of_independence
    """
    contingency = numpy.array(contingency, dtype=int, ndmin=2)
    assert contingency.ndim == 2
    n = float(numpy.sum(contingency))
    n0 = contingency.shape[0]
    n1 = contingency.shape[1]
    assert 0 < n0
    assert 0 < n1
    p0 = numpy.sum(contingency, axis=1)/n
    p1 = numpy.sum(contingency, axis=0)/n
    expected = n * numpy.outer(p0, p1)
    return numpy.sum(((contingency - expected)**2)/expected) 
Example 26
Project: tenpy   Author: tenpy   File: test_mps.py    License: GNU General Public License v3.0 6 votes vote down vote up
def test_expectation_value_multisite():
    s = spin_half
    psi = mps.MPS.from_singlets(s, 6, [(0, 1), (2, 3), (4, 5)], lonely=[], bc='finite')
    SpSm = npc.outer(s.Sp.replace_labels(['p', 'p*'], ['p0', 'p0*']),
                     s.Sm.replace_labels(['p', 'p*'], ['p1', 'p1*']))
    psi1 = psi.copy()
    ev = psi.expectation_value(SpSm)
    npt.assert_almost_equal(ev, [-0.5, 0., -0.5, 0., -0.5])
    env1 = mps.MPSEnvironment(psi1, psi)
    ev = env1.expectation_value(SpSm)
    npt.assert_almost_equal(ev, [-0.5, 0., -0.5, 0., -0.5])

    psi1.apply_local_op(2, SpSm)  # multi-site operator
    ev = psi1.expectation_value(SpSm)  # normalized!
    npt.assert_almost_equal(ev, [-0.5, 0., 0.0, 0., -0.5])
    env1 = mps.MPSEnvironment(psi1, psi)
    ev = env1.expectation_value(SpSm) / psi1.overlap(psi)  # normalize
    npt.assert_almost_equal(ev, [-0.5, 0., -1., 0., -0.5]) 
Example 27
Project: tenpy   Author: tenpy   File: test_purification.py    License: GNU General Public License v3.0 6 votes vote down vote up
def gen_disentangler_psi_prod(psiP, psiQ):
    """generate a PurificationMPS as tensorproduct (psi_P x psi_Q).

    psiQ should have the same `sites` as psiP.
    """
    L = psiP.L
    Bs = []
    for i in range(L):
        BP = psiP.get_B(i)
        BQ = psiQ.get_B(i)
        B2 = npc.outer(BP, BQ.conj())
        B2 = B2.combine_legs([['vL', 'vL*'], ['vR', 'vR*']], qconj=[+1, -1])
        B2.ireplace_labels(['(vL.vL*)', '(vR.vR*)', 'p*'], ['vL', 'vR', 'q'])
        Bs.append(B2)
    Ss = [np.outer(S, S2).flatten() for S, S2 in zip(psiP._S, psiQ._S)]
    return purification_mps.PurificationMPS(psiP.sites, Bs, Ss) 
Example 28
Project: pylops   Author: equinor   File: test_basicoperators.py    License: GNU Lesser General Public License v3.0 6 votes vote down vote up
def test_Flip2D(par):
    """Dot-test, forward and adjoint for Flip operator on 2d signal
    """
    np.random.seed(10)
    x = {}
    x['0'] = np.outer(np.arange(par['ny']), np.ones(par['nx'])) + \
             par['imag'] * np.outer(np.arange(par['ny']), np.ones(par['nx']))
    x['1'] = np.outer(np.ones(par['ny']), np.arange(par['nx'])) + \
             par['imag'] * np.outer(np.ones(par['ny']), np.arange(par['nx']))

    for dir in [0, 1]:
        Fop = Flip(par['ny']*par['nx'], dims=(par['ny'], par['nx']),
                   dir=dir, dtype=par['dtype'])
        assert dottest(Fop, par['ny']*par['nx'], par['ny']*par['nx'])

        y = Fop * x[str(dir)].flatten()
        xadj = Fop.H * y.flatten()
        xadj = xadj.reshape(par['ny'], par['nx'])
        assert_array_equal(x[str(dir)], xadj) 
Example 29
Project: pylops   Author: equinor   File: test_basicoperators.py    License: GNU Lesser General Public License v3.0 6 votes vote down vote up
def test_Symmetrize2D(par):
    """Dot-test, forward and inverse for Symmetrize operator on 2d signal
    """
    np.random.seed(10)
    x = {}
    x['0'] = np.outer(np.arange(par['ny']), np.ones(par['nx'])) + \
             par['imag'] * np.outer(np.arange(par['ny']), np.ones(par['nx']))
    x['1'] = np.outer(np.ones(par['ny']), np.arange(par['nx'])) + \
             par['imag'] * np.outer(np.ones(par['ny']), np.arange(par['nx']))

    for dir in [0, 1]:
        Sop = Symmetrize(par['ny']*par['nx'],
                         dims=(par['ny'], par['nx']),
                         dir=dir, dtype=par['dtype'])
        y = Sop * x[str(dir)].flatten()
        assert dottest(Sop, y.size, par['ny']*par['nx'])

        xinv = Sop / y
        assert_array_almost_equal(x[str(dir)].ravel(), xinv, decimal=3) 
Example 30
Project: FRIDA   Author: LCAV   File: point_cloud.py    License: MIT License 5 votes vote down vote up
def EDM(self):
        ''' Computes the EDM corresponding to the marker set '''
        if self.X is None:
            raise ValueError('No marker set')

        G = np.dot(self.X.T, self.X)
        return np.outer(np.ones(self.m), np.diag(G)) \
            - 2*G + np.outer(np.diag(G), np.ones(self.m))