Python numpy.isclose() Examples

The following are 30 code examples for showing how to use numpy.isclose(). 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: neuropythy   Author: noahbenson   File: models.py    License: GNU Affero General Public License v3.0 7 votes vote down vote up
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) 
Example 2
Project: neuropythy   Author: noahbenson   File: core.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
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 
Example 3
Project: neuropythy   Author: noahbenson   File: util.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
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)) 
Example 4
Project: neuropythy   Author: noahbenson   File: retinotopy.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
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) 
Example 5
Project: neuropythy   Author: noahbenson   File: retinotopy.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
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 
Example 6
Project: neuropythy   Author: noahbenson   File: core.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
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)]) 
Example 7
Project: OpenFermion-Cirq   Author: quantumlib   File: optimal_givens_decomposition_test.py    License: Apache License 2.0 6 votes vote down vote up
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) 
Example 8
Project: OpenFermion-Cirq   Author: quantumlib   File: objective.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example 9
Project: OpenFermion-Cirq   Author: quantumlib   File: objective_test.py    License: Apache License 2.0 6 votes vote down vote up
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) 
Example 10
Project: OpenFermion-Cirq   Author: quantumlib   File: gradient_hf_test.py    License: Apache License 2.0 6 votes vote down vote up
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) 
Example 11
Project: OpenFermion-Cirq   Author: quantumlib   File: swap_network_trotter.py    License: Apache License 2.0 6 votes vote down vote up
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) 
Example 12
Project: differential-privacy-library   Author: IBM   File: exponential.py    License: MIT License 6 votes vote down vote up
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 
Example 13
Project: differential-privacy-library   Author: IBM   File: geometric.py    License: MIT License 6 votes vote down vote up
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) 
Example 14
Project: differential-privacy-library   Author: IBM   File: vector.py    License: MIT License 6 votes vote down vote up
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 
Example 15
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() 
Example 16
Project: gated-graph-transformer-network   Author: hexahedria   File: ggtnn_train.py    License: MIT License 5 votes vote down vote up
def test_accuracy(m, story_buckets, bucket_sizes, num_answer_words, format_spec, batch_size, batch_auto_adjust=None, test_graph=False):
    correct = 0
    out_of = 0
    for bucket, bucket_size in zip(story_buckets, bucket_sizes):
        cur_batch_size = adj_size(m, bucket_size, batch_size, batch_auto_adjust)
        for start_idx in range(0, len(bucket), cur_batch_size):
            stories = bucket[start_idx:start_idx+cur_batch_size]
            batch = assemble_batch(stories, num_answer_words, format_spec)
            answers = batch[2]
            args = batch[:2] + ((answers.shape[1],) if format_spec == model.ModelOutputFormat.sequence else ())

            if test_graph:
                _, batch_close, _ = m.eval(*batch, with_accuracy=True)
            else:
                out_answers, out_strengths, out_ids, out_states, out_edges = m.snap_test_fn(*args)
                close = np.isclose(out_answers, answers)
                batch_close = np.all(close, (1,2))

            print(batch_close)

            batch_correct = np.sum(batch_close).tolist()
            batch_out_of = len(stories)
            correct +=  batch_correct
            out_of += batch_out_of

    return correct/out_of 
Example 17
Project: mmdetection   Author: open-mmlab   File: test_masks.py    License: Apache License 2.0 5 votes vote down vote up
def test_polygon_mask_area():
    # area of empty polygon masks
    raw_masks = dummy_raw_polygon_masks((0, 28, 28))
    polygon_masks = PolygonMasks(raw_masks, 28, 28)
    assert polygon_masks.areas.sum() == 0

    # area of polygon masks contain 1 instance
    # here we hack a case that the gap between the area of bitmap and polygon
    # is minor
    raw_masks = [[np.array([1, 1, 5, 1, 3, 4])]]
    polygon_masks = PolygonMasks(raw_masks, 6, 6)
    polygon_area = polygon_masks.areas
    bitmap_area = polygon_masks.to_bitmap().areas
    assert len(polygon_area) == 1
    assert np.isclose(polygon_area, bitmap_area).all() 
Example 18
Project: neural-fingerprinting   Author: StephanZheng   File: test_attacks.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def test_generate_np_high_confidence_untargeted_examples(self):

        trivial_model = TrivialModel()

        for CONFIDENCE in [0, 2.3]:
            x_val = np.random.rand(10, 1) - .5
            x_val = np.array(x_val, dtype=np.float32)

            orig_labs = np.argmax(self.sess.run(trivial_model.get_logits(x_val)), axis=1)
            attack = CarliniWagnerL2(trivial_model, sess=self.sess)
            x_adv = attack.generate_np(x_val,
                                       max_iterations=100,
                                       binary_search_steps=2,
                                       learning_rate=1e-2,
                                       initial_const=1,
                                       clip_min=-10, clip_max=10,
                                       confidence=CONFIDENCE,
                                       batch_size=10)

            new_labs = self.sess.run(trivial_model.get_logits(x_adv))

            good_labs = new_labs[np.arange(10), 1 - orig_labs]
            bad_labs = new_labs[np.arange(10), orig_labs]

            self.assertTrue(np.mean(np.argmax(new_labs, axis=1) == orig_labs)
                            == 0)
            self.assertTrue(np.isclose(
                0, np.min(good_labs - (bad_labs + CONFIDENCE)), atol=1e-1)) 
Example 19
Project: neural-fingerprinting   Author: StephanZheng   File: test_attacks.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def test_generate_np_high_confidence_targeted_examples(self):

        trivial_model = TrivialModel()

        for CONFIDENCE in [0, 2.3]:
            x_val = np.random.rand(10, 1) - .5
            x_val = np.array(x_val, dtype=np.float32)

            feed_labs = np.zeros((10, 2))
            feed_labs[np.arange(10), np.random.randint(0, 2, 10)] = 1
            attack = CarliniWagnerL2(trivial_model, sess=self.sess)
            x_adv = attack.generate_np(x_val,
                                       max_iterations=100,
                                       binary_search_steps=2,
                                       learning_rate=1e-2,
                                       initial_const=1,
                                       clip_min=-10, clip_max=10,
                                       confidence=CONFIDENCE,
                                       y_target=feed_labs,
                                       batch_size=10)

            new_labs = self.sess.run(trivial_model.get_logits(x_adv))

            good_labs = new_labs[np.arange(10), np.argmax(feed_labs, axis=1)]
            bad_labs = new_labs[np.arange(
                10), 1 - np.argmax(feed_labs, axis=1)]

            self.assertTrue(np.isclose(
                0, np.min(good_labs - (bad_labs + CONFIDENCE)), atol=1e-1))
            self.assertTrue(np.mean(np.argmax(new_labs, axis=1) ==
                                    np.argmax(feed_labs, axis=1)) > .9) 
Example 20
Project: neural-fingerprinting   Author: StephanZheng   File: test_attacks.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def test_generate_np_high_confidence_untargeted_examples(self):

        trivial_model = TrivialModel()

        for CONFIDENCE in [0, 2.3]:
            x_val = np.random.rand(10, 1) - .5
            x_val = np.array(x_val, dtype=np.float32)

            orig_labs = np.argmax(self.sess.run(trivial_model.get_logits(x_val)), axis=1)
            attack = CarliniWagnerL2(trivial_model, sess=self.sess)
            x_adv = attack.generate_np(x_val,
                                       max_iterations=100,
                                       binary_search_steps=2,
                                       learning_rate=1e-2,
                                       initial_const=1,
                                       clip_min=-10, clip_max=10,
                                       confidence=CONFIDENCE,
                                       batch_size=10)

            new_labs = self.sess.run(trivial_model.get_logits(x_adv))

            good_labs = new_labs[np.arange(10), 1 - orig_labs]
            bad_labs = new_labs[np.arange(10), orig_labs]

            self.assertTrue(np.mean(np.argmax(new_labs, axis=1) == orig_labs)
                            == 0)
            self.assertTrue(np.isclose(
                0, np.min(good_labs - (bad_labs + CONFIDENCE)), atol=1e-1)) 
Example 21
Project: models   Author: kipoi   File: test_basset_model.py    License: MIT License 5 votes vote down vote up
def test_ref_seq():
    # Get pure fasta predictions
    model_dir = model_root + "./"
    model = kipoi.get_model(model_dir, source="dir")
    # The preprocessor
    Dataloader = kipoi.get_dataloader_factory(model_dir, source="dir")
    dataloader_arguments = {
        "fasta_file": "/nfs/research1/stegle/users/rkreuzhu/opt/manuscript_code/data/raw/dataloader_files/shared/hg19.fa",
        "intervals_file": "test_files/test_encode_roadmap.bed"
    }
    # predict using results
    preds = model.pipeline.predict(dataloader_arguments)
    #
    res_orig = pd.read_csv("/nfs/research1/stegle/users/rkreuzhu/deeplearning/Basset/data/test_encode_roadmap_short_pred.txt", "\t", header=None)
    assert np.isclose(preds, res_orig.values, atol=1e-3).all() 
Example 22
Project: neuropythy   Author: noahbenson   File: core.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def white_to_pial_vectors(white_surface, pial_surface):
        '''
        cortex.white_to_pial_vectors is a (3 x n) matrix of the unit direction vectors that point
          from the n vertices in the cortex's white surface to their equivalent positions in the
          pial surface.
        '''
        u = pial_surface.coordinates - white_surface.coordinates
        d = np.sqrt(np.sum(u**2, axis=0))
        z = np.isclose(d, 0)
        return pimms.imm_array(np.logical_not(z) / (d + z)) 
Example 23
Project: neuropythy   Author: noahbenson   File: util.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def normalize(u):
    '''
    normalize(u) yields a vetor with the same direction as u but unit length, or, if u has zero
    length, yields u.
    '''
    u = np.asarray(u)
    unorm = np.sqrt(np.sum(u**2, axis=0))
    z = np.isclose(unorm, 0)
    c = np.logical_not(z) / (unorm + z)
    return u * c 
Example 24
Project: neuropythy   Author: noahbenson   File: util.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def point_on_line(ab, c, atol=1e-8):
    '''
    point_on_line((a,b), c) yields True if point x is on line (a,b) and False otherwise.
    '''
    (a,b) = ab
    abc = [np.asarray(u) for u in (a,b,c)]
    if any(len(u.shape) == 2 for u in abc): (a,b,c) = [np.reshape(u,(len(u),-1)) for u in abc]
    else:                                   (a,b,c) = abc
    vca = a - c
    vcb = b - c
    uba = czdivide(vba, np.sqrt(np.sum(vba**2, axis=0)))
    uca = czdivide(vca, np.sqrt(np.sum(vca**2, axis=0)))
    return (np.isclose(np.sqrt(np.sum(vca**2, axis=0)), 0, atol=atol) |
            np.isclose(np.sqrt(np.sum(vcb**2, axis=0)), 0, atol=atol) |
            np.isclose(np.abs(np.sum(uba*uca, axis=0)), 1, atol=atol)) 
Example 25
Project: neuropythy   Author: noahbenson   File: util.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def points_close(a,b, atol=1e-8):
    '''
    points_close(a,b) yields True if points a and b are close to each other and False otherwise.
    '''
    (a,b) = [np.asarray(u) for u in (a,b)]
    if len(a.shape) == 2 or len(b.shape) == 2: (a,b) = [np.reshape(u,(len(u),-1)) for u in (a,b)]
    return np.isclose(np.sqrt(np.sum((a - b)**2, axis=0)), 0, atol=atol) 
Example 26
Project: neuropythy   Author: noahbenson   File: util.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def line_intersection_2D(abarg, cdarg, atol=1e-8):
    '''
    line_intersection((a, b), (c, d)) yields the intersection point between the lines that pass
    through the given pairs of points. If any lines are parallel, (numpy.nan, numpy.nan) is
    returned; note that a, b, c, and d can all be 2 x n matrices of x and y coordinate row-vectors.
    '''
    ((x1,y1),(x2,y2)) = abarg
    ((x3,y3),(x4,y4)) = cdarg
    dx12 = (x1 - x2)
    dx34 = (x3 - x4)
    dy12 = (y1 - y2)
    dy34 = (y3 - y4)
    denom = dx12*dy34 - dy12*dx34
    unit = np.isclose(denom, 0, atol=atol)
    if unit is True: return (np.nan, np.nan)
    denom = unit + denom
    q12 = (x1*y2 - y1*x2) / denom
    q34 = (x3*y4 - y3*x4) / denom
    xi = q12*dx34 - q34*dx12
    yi = q12*dy34 - q34*dy12
    if   unit is False: return (xi, yi)
    elif unit is True:  return (np.nan, np.nan)
    else:
        xi = np.asarray(xi)
        yi = np.asarray(yi)
        xi[unit] = np.nan
        yi[unit] = np.nan
        return (xi, yi) 
Example 27
Project: neuropythy   Author: noahbenson   File: util.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def point_in_tetrahedron(tetra, pt):
    '''
    point_in_tetrahedron(tetrahedron, point) yields True if the given point is in the given 
      tetrahedron. If either tetrahedron or point (or both) are lists of shapes/points, then this
      calculation is automatically threaded over all the given arguments.
    '''
    bcs = tetrahedral_barycentric_coordinates(tetra, pt)
    return np.logical_not(np.all(np.isclose(bcs, 0), axis=0)) 
Example 28
Project: neuropythy   Author: noahbenson   File: util.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def point_in_prism(tri1, tri2, pt):
    '''
    point_in_prism(tri1, tri2, pt) yields True if the given point is inside the prism that stretches
      between triangle 1 and triangle 2. Will automatically thread over extended dimensions. If
      multiple triangles are given, then the vertices must be an earlier dimension than the
      coordinates; e.g., a 3 x 3 x n array will be assumed to organized such that element [0,1,k] is
      the y coordinate of the first vertex of the k'th triangle.
    '''
    bcs = prism_barycentric_coordinates(tri1, tri2, pt)
    return np.logical_not(np.isclose(np.sum(bcs, axis=0), 0)) 
Example 29
Project: neuropythy   Author: noahbenson   File: core.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def __add__(self, x):
        x = to_potential(x)
        if is_const_potential(x):
            if   np.isclose(x.c, 0).all(): return self
            elif is_const_potential(self): return const_potential(self.c + x.c)
            else:                          return PotentialPlusConstant(self, x.c)
        elif is_const_potential(self):     return PotentialPlusConstant(x.c, self)
        else:                              return PotentialPlusPotential(self, x) 
Example 30
Project: neuropythy   Author: noahbenson   File: core.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def __mul__(self, x):
        x = to_potential(x)
        if is_const_potential(x):
            if   np.isclose(x.c, 1).all(): return self
            elif np.isclose(x.c, 0).all(): return x
            elif is_const_potential(self): return const_potential(self.c * x.c)
            else:                          return PotentialTimesConstant(self, x.c)
        elif is_const_potential(self):     return PotentialTimesConstant(x, self.c)
        else:                              return PotentialTimesPotential(self, x)