Python numpy.isclose() Examples

def parameters(params):
        '''
        mdl.parameters is a persistent map of the parameters for the given SchiraModel object mdl.
        '''
        if not pimms.is_pmap(params): params = pyr.pmap(params)
        # do the translations that we need...
        scale = params['scale']
        if pimms.is_number(scale):
            params = params.set('scale', (scale, scale))
        elif not is_tuple(scale):
            params = params.set('scale', tuple(scale))
        shear = params['shear']
        if pimms.is_number(shear) and np.isclose(shear, 0):
            params = params.set('shear', ((1, 0), (0, 1)))
        elif shear[0][0] != 1 or shear[1][1] != 1:
            raise RuntimeError('shear matrix diagonal elements must be 1!')
        elif not is_tuple(shear) or not all(is_tuple(s) for s in shear):
            params.set('shear', tuple([tuple(s) for s in shear]))
        center = params['center']
        if pimms.is_number(center) and np.isclose(center, 0):
            params = params.set('center', (0.0, 0.0))
        return pimms.persist(params, depth=None) 

def test_circuit_generation_and_accuracy():
    for dim in range(2, 10):
        qubits = cirq.LineQubit.range(dim)
        u_generator = numpy.random.random(
            (dim, dim)) + 1j * numpy.random.random((dim, dim))
        u_generator = u_generator - numpy.conj(u_generator).T
        assert numpy.allclose(-1 * u_generator, numpy.conj(u_generator).T)

        unitary = scipy.linalg.expm(u_generator)
        circuit = cirq.Circuit()
        circuit.append(optimal_givens_decomposition(qubits, unitary))

        fermion_generator = QubitOperator(()) * 0.0
        for i, j in product(range(dim), repeat=2):
            fermion_generator += jordan_wigner(
                FermionOperator(((i, 1), (j, 0)), u_generator[i, j]))

        true_unitary = scipy.linalg.expm(
            get_sparse_operator(fermion_generator).toarray())
        assert numpy.allclose(true_unitary.conj().T.dot(true_unitary),
                              numpy.eye(2 ** dim, dtype=complex))

        test_unitary = cirq.unitary(circuit)
        assert numpy.isclose(
            abs(numpy.trace(true_unitary.conj().T.dot(test_unitary))), 2 ** dim) 

def params(self) -> Iterable[sympy.Symbol]:
        """The parameters of the ansatz."""
        for i in range(self.iterations):
            for p in range(len(self.qubits)):
                if (self.include_all_z or not
                        numpy.isclose(self.hamiltonian.one_body[p, p], 0)):
                    yield LetterWithSubscripts('U', p, i)
            for p, q in itertools.combinations(range(len(self.qubits)), 2):
                if (self.include_all_xxyy or not
                        numpy.isclose(self.hamiltonian.one_body[p, q].real, 0)):
                    yield LetterWithSubscripts('T', p, q, i)
                if (self.include_all_yxxy or not
                        numpy.isclose(self.hamiltonian.one_body[p, q].imag, 0)):
                    yield LetterWithSubscripts('W', p, q, i)
                if (self.include_all_cz or not
                        numpy.isclose(self.hamiltonian.two_body[p, q], 0)):
                    yield LetterWithSubscripts('V', p, q, i) 

def calc_anchors(preregistration_map, model, model_hemi,
                 scale=1, sigma=Ellipsis, radius_weight=0, field_sign_weight=0,
                 invert_rh_field_sign=False):
    '''
    calc_anchors is a calculator that creates a set of anchor instructions for a registration.

    Required afferent parameters:
      @ invert_rh_field_sign May be set to True (default is False) to indicate that the right
        hemisphere's field signs will be incorrect relative to the model; this generally should be
        used whenever invert_rh_angle is also set to True.

    '''
    wgts = preregistration_map.prop('weight')
    rads = preregistration_map.prop('radius')
    if np.isclose(radius_weight, 0): radius_weight = 0
    ancs = retinotopy_anchors(preregistration_map, model,
                              polar_angle='polar_angle',
                              eccentricity='eccentricity',
                              radius='radius',
                              weight=wgts, weight_min=0, # taken care of already
                              radius_weight=radius_weight, field_sign_weight=field_sign_weight,
                              scale=scale,
                              invert_field_sign=(model_hemi == 'rh' and invert_rh_field_sign),
                              **({} if sigma is Ellipsis else {'sigma':sigma}))
    return ancs 

def _retinotopic_field_sign_triangles(m, retinotopy):
    t = m.tess if isinstance(m, geo.Mesh) or isinstance(m, geo.Topology) else m
    # get the polar angle and eccen data as a complex number in degrees
    if pimms.is_str(retinotopy):
        (x,y) = as_retinotopy(retinotopy_data(m, retinotopy), 'geographical')
    elif retinotopy is Ellipsis:
        (x,y) = as_retinotopy(retinotopy_data(m, 'any'),      'geographical')
    else:
        (x,y) = as_retinotopy(retinotopy,                     'geographical')
    # Okay, now we want to make some coordinates...
    coords = np.asarray([x, y])
    us = coords[:, t.indexed_faces[1]] - coords[:, t.indexed_faces[0]]
    vs = coords[:, t.indexed_faces[2]] - coords[:, t.indexed_faces[0]]
    (us,vs) = [np.concatenate((xs, np.full((1, t.face_count), 0.0))) for xs in [us,vs]]
    xs = np.cross(us, vs, axis=0)[2]
    xs[np.isclose(xs, 0)] = 0
    return np.sign(xs) 

def test_rhf_func_gen():
    rhf_objective, molecule, parameters, _, _ = make_h6_1_3()
    ansatz, energy, _ = rhf_func_generator(rhf_objective)
    assert np.isclose(molecule.hf_energy, energy(parameters))

    ansatz, energy, _, opdm_func = rhf_func_generator(
        rhf_objective, initial_occ_vec=[1] * 3 + [0] * 3, get_opdm_func=True)
    assert np.isclose(molecule.hf_energy, energy(parameters))
    test_opdm = opdm_func(parameters)
    u = ansatz(parameters)
    initial_opdm = np.diag([1] * 3 + [0] * 3)
    final_odpm = u @ initial_opdm @ u.T
    assert np.allclose(test_opdm, final_odpm)

    result = rhf_minimization(rhf_objective, initial_guess=parameters)
    assert np.allclose(result.x, parameters) 

def point_in_segment(ac, b, atol=1e-8):
    '''
    point_in_segment((a,b), c) yields True if point x is in segment (a,b) and False otherwise. Note
    that this differs from point_on_segment in that a point that if c is equal to a or b it is
    considered 'on' but not 'in' the segment.
    The option atol can be given and is used only to test for difference from 0; by default it is
    1e-8.
    '''
    (a,c) = ac
    abc = [np.asarray(u) for u in (a,b,c)]
    if any(len(u.shape) > 1 for u in abc): (a,b,c) = [np.reshape(u,(len(u),-1)) for u in abc]
    else:                                  (a,b,c) = abc
    vab = b - a
    vbc = c - b
    vac = c - a
    dab = np.sqrt(np.sum(vab**2, axis=0))
    dbc = np.sqrt(np.sum(vbc**2, axis=0))
    dac = np.sqrt(np.sum(vac**2, axis=0))
    return (np.isclose(dab + dbc - dac, 0, atol=atol) &
            ~np.isclose(dac - dab, 0, atol=atol) &
            ~np.isclose(dac - dbc, 0, atol=atol)) 

def test_get_matrix_of_eigs():
    """
    Generate the matrix of [exp(i (li - lj)) - 1] / (i(li - lj)
    :return:
    """
    lam_vals = np.random.randn(4) + 1j * np.random.randn(4)
    mat_eigs = np.zeros((lam_vals.shape[0],
                         lam_vals.shape[0]),
                         dtype=np.complex128)
    for i, j in product(range(lam_vals.shape[0]), repeat=2):
        if np.isclose(abs(lam_vals[i] - lam_vals[j]), 0):
            mat_eigs[i, j] = 1
        else:
            mat_eigs[i, j] = (np.exp(1j * (lam_vals[i] - lam_vals[j])) - 1) / (
                        1j * (lam_vals[i] - lam_vals[j]))

    test_mat_eigs = get_matrix_of_eigs(lam_vals)
    assert np.allclose(test_mat_eigs, mat_eigs) 

def splrep(coordinates, t, order, weights, smoothing, periodic):
        from scipy import interpolate
        (x,y) = coordinates
        # we need to skip anything where t[i] and t[i+1] are too close
        wh = np.where(np.isclose(np.diff(t), 0))[0]
        if len(wh) > 0:
            (t,x,y) = [np.array(u) for u in (t,x,y)]
            ii = np.arange(len(t))
            for i in reversed(wh):
                ii[i+1:-1] = ii[i+2:]
                for u in (t,x,y):
                    u[i] = np.mean(u[i:i+2])
            ii = ii[:-len(wh)]
            (t,x,y) = [u[ii] for u in (t,x,y)]
        xtck = interpolate.splrep(t, x, k=order, s=smoothing, w=weights, per=periodic)
        ytck = interpolate.splrep(t, y, k=order, s=smoothing, w=weights, per=periodic)
        return tuple([tuple([pimms.imm_array(u) for u in tck])
                      for tck in (xtck,ytck)]) 

def _build_normalising_constant(self, re_eval=False):
        balanced_tree = True
        first_constant_value = None
        normalising_constant = {}

        for _base_leaf in self._domain_values:
            constant_value = 0.0

            for _target_leaf in self._domain_values:
                constant_value += self._get_prob(_base_leaf, _target_leaf)

            normalising_constant[_base_leaf] = constant_value

            if first_constant_value is None:
                first_constant_value = constant_value
            elif not np.isclose(constant_value, first_constant_value):
                balanced_tree = False

        # If the tree is balanced, we can eliminate the doubling factor
        if balanced_tree and not re_eval:
            self._balanced_tree = True
            return self._build_normalising_constant(True)

        return normalising_constant 

def color_overlap(color1, *args):
    '''
    color_overlap(color1, color2...) yields the rgba value associated with overlaying color2 on top
      of color1 followed by any additional colors (overlaid left to right). This respects alpha
      values when calculating the results.
    Note that colors may be lists of colors, in which case a matrix of RGBA values is yielded.
    '''
    args = list(args)
    args.insert(0, color1)
    rgba = np.asarray([0.5,0.5,0.5,0])
    for c in args:
        c = to_rgba(c)
        a = c[...,3]
        a0 = rgba[...,3]
        if   np.isclose(a0, 0).all(): rgba = np.ones(rgba.shape) * c
        elif np.isclose(a,  0).all(): continue
        else:                         rgba = times(a, c) + times(1-a, rgba)
    return rgba 

def set_bounds(self, lower, upper):
        """Sets the lower and upper bounds of the mechanism.

        For the folded geometric mechanism, `lower` and `upper` must be integer or half-integer -valued.  Must have
        `lower` <= `upper`.

        Parameters
        ----------
        lower : int or float
            The lower bound of the mechanism.

        upper : int or float
            The upper bound of the mechanism.

        Returns
        -------
        self : class

        """
        if not np.isclose(2 * lower, np.round(2 * lower)) or not np.isclose(2 * upper, np.round(2 * upper)):
            raise ValueError("Bounds must be integer or half-integer floats")

        return super().set_bounds(lower, upper) 

def set_dimension(self, vector_dim):
        """Sets the dimension `vector_dim` of the domain of the mechanism.

        This dimension relates to the size of the input vector of the function being considered by the mechanism.  This
        corresponds to the size of the random vector produced by the mechanism.

        Parameters
        ----------
        vector_dim : int
            Function input dimension.

        Returns
        -------
        self : class

        """
        if not isinstance(vector_dim, Real) or not np.isclose(vector_dim, int(vector_dim)):
            raise TypeError("d must be integer-valued")
        if not vector_dim >= 1:
            raise ValueError("d must be strictly positive")

        self._vector_dim = int(vector_dim)
        return self 

def get_matrix_of_eigs(w: np.ndarray) -> np.ndarray:
    """
    Transform the eigenvalues for getting the gradient

    .. math:
        f(w) \rightarrow \frac{e^{i (\lambda_{i} - \lambda_{j})}{i (\lambda_{i} - \lambda_{j})}

    :param w: eigenvalues of C-matrix
    :return: new array of transformed eigenvalues
    """
    transform_eigs = np.zeros((w.shape[0], w.shape[0]),
                                 dtype=np.complex128)
    for i, j in product(range(w.shape[0]), repeat=2):
        if np.isclose(abs(w[i] - w[j]), 0):
            transform_eigs[i, j] = 1
        else:
            transform_eigs[i, j] = (np.exp(1j * (w[i] - w[j])) - 1) / (
                        1j * (w[i] - w[j]))
    return transform_eigs 

def test_iris(self):
        from sklearn import datasets
        dataset = datasets.load_iris()

        X_train, y_train = dataset.data, dataset.target

        norms = np.linalg.norm(X_train, axis=1)
        clip = (norms[0] + norms[1]) / 2

        X_clipped = clip_to_norm(X_train, clip)
        clipped_norms = np.linalg.norm(X_clipped, axis=1)
        self.assertLessEqual(clipped_norms[0], norms[0])
        self.assertLessEqual(clipped_norms[1], norms[1])
        self.assertTrue(np.isclose(clipped_norms[0], clip) or np.isclose(clipped_norms[1], clip)) 

def test_max_interval_tree():
    tree = MinSegmentTree(4)

    tree[0] = 1.0
    tree[2] = 0.5
    tree[3] = 3.0

    assert np.isclose(tree.min(), 0.5)
    assert np.isclose(tree.min(0, 2), 1.0)
    assert np.isclose(tree.min(0, 3), 0.5)
    assert np.isclose(tree.min(0, -1), 0.5)
    assert np.isclose(tree.min(2, 4), 0.5)
    assert np.isclose(tree.min(3, 4), 3.0)

    tree[2] = 0.7

    assert np.isclose(tree.min(), 0.7)
    assert np.isclose(tree.min(0, 2), 1.0)
    assert np.isclose(tree.min(0, 3), 0.7)
    assert np.isclose(tree.min(0, -1), 0.7)
    assert np.isclose(tree.min(2, 4), 0.7)
    assert np.isclose(tree.min(3, 4), 3.0)

    tree[2] = 4.0

    assert np.isclose(tree.min(), 1.0)
    assert np.isclose(tree.min(0, 2), 1.0)
    assert np.isclose(tree.min(0, 3), 1.0)
    assert np.isclose(tree.min(0, -1), 1.0)
    assert np.isclose(tree.min(2, 4), 3.0)
    assert np.isclose(tree.min(2, 3), 4.0)
    assert np.isclose(tree.min(2, -1), 4.0)
    assert np.isclose(tree.min(3, 4), 3.0) 

def test_tree_set_overlap():
    tree = SumSegmentTree(4)

    tree[2] = 1.0
    tree[2] = 3.0

    assert np.isclose(tree.sum(), 3.0)
    assert np.isclose(tree.sum(2, 3), 3.0)
    assert np.isclose(tree.sum(2, -1), 3.0)
    assert np.isclose(tree.sum(2, 4), 3.0)
    assert np.isclose(tree.sum(1, 2), 0.0) 

def test_simple(self):
        global_seed(3141592653)
        clf = KMeans(5, (0, 1), 3)

        X = np.zeros(1000) + 0.1
        X[:666] = 0.5
        X[:333] = 0.9
        X = X.reshape(-1, 1)

        clf.fit(X)
        centers = clf.cluster_centers_

        self.assertTrue(np.isclose(centers, 0.1, atol=0.05).any())
        self.assertTrue(np.isclose(centers, 0.5, atol=0.05).any())
        self.assertTrue(np.isclose(centers, 0.9, atol=0.05).any()) 

def _print_detection_eval_metrics(self, coco_eval):
    IoU_lo_thresh = 0.5
    IoU_hi_thresh = 0.95

    def _get_thr_ind(coco_eval, thr):
      ind = np.where((coco_eval.params.iouThrs > thr - 1e-5) &
                     (coco_eval.params.iouThrs < thr + 1e-5))[0][0]
      iou_thr = coco_eval.params.iouThrs[ind]
      assert np.isclose(iou_thr, thr)
      return ind

    ind_lo = _get_thr_ind(coco_eval, IoU_lo_thresh)
    ind_hi = _get_thr_ind(coco_eval, IoU_hi_thresh)
    # precision has dims (iou, recall, cls, area range, max dets)
    # area range index 0: all area ranges
    # max dets index 2: 100 per image
    precision = \
      coco_eval.eval['precision'][ind_lo:(ind_hi + 1), :, :, 0, 2]
    ap_default = np.mean(precision[precision > -1])
    print(('~~~~ Mean and per-category AP @ IoU=[{:.2f},{:.2f}] '
           '~~~~').format(IoU_lo_thresh, IoU_hi_thresh))
    print('{:.1f}'.format(100 * ap_default))
    for cls_ind, cls in enumerate(self.classes):
      if cls == '__background__':
        continue
      # minus 1 because of __background__
      precision = coco_eval.eval['precision'][ind_lo:(ind_hi + 1), :, cls_ind - 1, 0, 2]
      ap = np.mean(precision[precision > -1])
      print('{:.1f}'.format(100 * ap))

    print('~~~~ Summary metrics ~~~~')
    coco_eval.summarize() 

def test_tree_set():
    tree = SumSegmentTree(4)

    tree[2] = 1.0
    tree[3] = 3.0

    assert np.isclose(tree.sum(), 4.0)
    assert np.isclose(tree.sum(0, 2), 0.0)
    assert np.isclose(tree.sum(0, 3), 1.0)
    assert np.isclose(tree.sum(2, 3), 1.0)
    assert np.isclose(tree.sum(2, -1), 1.0)
    assert np.isclose(tree.sum(2, 4), 4.0) 

def map_to_matrix(mat):
    if mat.ndim != 4:
        raise TypeError("I only map rank-4 tensors to matices with symmetric support")
    dim = mat.shape[0]
    matform = np.zeros((dim**2, dim**2))
    for p, q, r, s in product(range(dim), repeat=4):
        assert np.isclose(mat[p, q, r, s].imag, 0.0)
        matform[p * dim + q, r * dim + s] = mat[p, q, r, s].real
    return matform 

def higham_root(eigenvalues, target_trace, epsilon=1.0E-15):
    """
    Find the root of f(sigma) = sum_{j}Theta(l_{i} - sigma)(l_{i} - sigma) = T

    :param eigenvalues: ordered list of eigenvalues from least to greatest
    :param target_trace: trace to maintain on new matrix
    :param epsilon: precision on bisection linesearch
    """
    if target_trace < 0.0:
        raise ValueError("Target trace needs to be a non-negative number")

    # when we want the trace to be zero
    if np.isclose(target_trace, 0.0):
        return eigenvalues[-1]

    # find top sigma
    sigma = eigenvalues[-1]
    while higham_polynomial(eigenvalues, sigma) < target_trace:
        sigma -= eigenvalues[-1]

    sigma_low = sigma
    sigma_high = eigenvalues[-1]

    while sigma_high - sigma_low >= epsilon:
        midpoint = sigma_high - (sigma_high - sigma_low) / 2.0
        if higham_polynomial(eigenvalues, midpoint) < target_trace:
            sigma_high = midpoint
        else:
            sigma_low = midpoint

    return sigma_high 

def test_reconstruction():
    dim = 20
    np.random.seed(42)
    mat = np.random.random((dim, dim))
    mat = 0.5 * (mat + mat.T)
    test_mat = np.zeros_like(mat)
    w, v = np.linalg.eigh(mat)
    for i in range(w.shape[0]):
        test_mat += w[i] * v[:, [i]].dot(v[:, [i]].T)
    assert np.allclose(test_mat - mat, 0.0)

    test_mat = fixed_trace_positive_projection(mat, np.trace(mat))
    assert np.isclose(np.trace(test_mat), np.trace(mat))
    w, v = np.linalg.eigh(test_mat)
    assert np.all(w >= -(float(4.0E-15)))

    mat = np.arange(16).reshape((4, 4))
    mat = 0.5 * (mat + mat.T)
    mat_tensor = map_to_tensor(mat)
    trace_mat = np.trace(mat)
    true_mat = fixed_trace_positive_projection(mat, trace_mat)
    test_mat = map_to_matrix(fixed_trace_positive_projection(mat_tensor,
                                                             trace_mat))
    assert np.allclose(true_mat, test_mat)

    w, v = np.linalg.eigh(true_mat)
    assert np.allclose(true_mat, fixed_trace_positive_projection(true_mat,
                                                                 trace_mat)) 

def test_higham_root():
    dim = 20
    np.random.seed(42)
    mat = np.random.random((dim, dim))
    mat = 0.5 * (mat + mat.T)
    w, _ = np.linalg.eigh(mat)
    target_trace = np.round(w[-1]-1)
    sigma = higham_root(w, target_trace)
    assert np.isclose(higham_polynomial(w, shift=sigma), target_trace)

    with pytest.raises(ValueError):
        higham_root(w, target_trace=-1)

    tw = higham_root(w, target_trace=0)
    assert np.isclose(tw, w[-1]) 

def test_highham_polynomial():
    eigs = np.arange(10)
    assert np.isclose(higham_polynomial(eigs, eigs[-1]), 0.0)
    assert np.isclose(higham_polynomial(eigs, 0), sum(eigs))
    assert np.isclose(higham_polynomial(eigs, 5), sum(eigs[5:] - 5))
    assert np.isclose(higham_polynomial(eigs, 8), sum(eigs[8:] - 8)) 

def test_heaviside():
    assert np.isclose(heaviside(0), 1.0)
    assert np.isclose(heaviside(0.5), 1.0)
    assert np.isclose(heaviside(-0.5), 0.0)
    assert np.isclose(heaviside(-0.5, -1), 1.0)
    assert np.isclose(heaviside(-2, -1), 0) 

def test_fidelity():
    parameters = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
    u = sp.linalg.expm(rhf_params_to_matrix(parameters, 6))
    opdm = np.array([[0.766034130, -0.27166330, -0.30936072, -0.08471057, -0.04878244, -0.01285432],
                     [-0.27166330,  0.67657015, -0.37519640, -0.02101843, -0.03568214, -0.05034585],
                     [-0.30936072, -0.37519640,  0.55896791,  0.04267370, -0.02258184, -0.08783738],
                     [-0.08471057, -0.02101843,  0.04267370,  0.05450848,  0.11291253,  0.17131658],
                     [-0.04878244, -0.03568214, -0.02258184,  0.11291253,  0.26821219,  0.42351185],
                     [-0.01285432, -0.05034585, -0.08783738,  0.17131658,  0.42351185,  0.67570713]])

    assert np.isclose(fidelity(u, opdm), 1.0)
    opdm += 0.1
    opdm = 0.5 * (opdm + opdm.T)
    assert np.isclose(fidelity(u, opdm), 0.3532702370138279) 

def test_trace_distance():
    rho = np.arange(16).reshape((4, 4))
    sigma = np.arange(16, 32).reshape((4, 4))
    assert np.isclose(trace_distance(rho, rho), 0.)
    assert np.isclose(trace_distance(rho, sigma), 32.0) 

def test_kdelta():
    assert np.isclose(kdelta(1, 1), 1.)
    assert np.isclose(kdelta(0, 1), 0.) 

def prepare_slater_determinant(qubits: List[cirq.Qid],
                               slater_determinant_matrix: np.ndarray,
                               clean_ryxxy: Optional[Union[bool, int]] = True):
    """
    High level interface to the real basis rotation circuit generator

    :param qubits: List of cirq.Qids denoting logical qubits
    :param slater_determinant_matrix: basis rotation matrix
    :param clean_ryxxy: Optional[True, 1, 2, 3, 4] for indicating an error
                        model of the Givens rotation.
    :return: generator for circuit
    """
    circuit_description = slater_determinant_preparation_circuit(slater_determinant_matrix)
    yield (cirq.X(qubits[j]) for j in range(slater_determinant_matrix.shape[0]))
    for parallel_ops in circuit_description:
        for op in parallel_ops:
            i, j, theta, phi = op
            if not np.isclose(phi, 0):  # testpragma: no cover
                raise ValueError("unitary must be real valued only")
            if clean_ryxxy is True or clean_ryxxy == 1:
                yield ryxxy(qubits[i], qubits[j], theta)
            elif clean_ryxxy == 2:
                yield ryxxy2(qubits[i], qubits[j], theta)
            elif clean_ryxxy == 3:
                yield ryxxy3(qubits[i], qubits[j], theta)
            elif clean_ryxxy == 4:
                yield ryxxy3(qubits[i], qubits[j], theta)
            else:
                raise ValueError("Invalide clean_ryxxy value")