Python numpy.arccos() Examples

The following are 30 code examples of numpy.arccos(). 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 also want to check out all available functions/classes of the module numpy , or try the search function .
Example #1
Source File: util.py    From neuropythy with GNU Affero General Public License v3.0 7 votes vote down vote up
def vector_angle(u, v, direction=None):
    '''
    vector_angle(u, v) yields the angle between the two vectors u and v. The optional argument 
    direction is by default None, which specifies that the smallest possible angle between the
    vectors be reported; if the vectors u and v are 2D vectors and direction parameters True and
    False specify the clockwise or counter-clockwise directions, respectively; if the vectors are
    3D vectors, then direction may be a 3D point that is not in the plane containing u, v, and the
    origin, and it specifies around which direction (u x v or v x u) the the counter-clockwise angle
    from u to v should be reported (the cross product vector that has a positive dot product with
    the direction argument is used as the rotation axis).
    '''
    if direction is None:
        return np.arccos(vector_angle_cos(u, v))
    elif direction is True:
        return np.arctan2(v[1], v[0]) - np.arctan2(u[1], u[0])
    elif direction is False:
        return np.arctan2(u[1], u[0]) - np.arctan2(v[1], v[0])
    else:
        axis1 = normalize(u)
        axis2 = normalize(np.cross(u, v))
        if np.dot(axis2, direction) < 0:
            axis2 = -axis2
        return np.arctan2(np.dot(axis2, v), np.dot(axis1, v)) 
Example #2
Source File: test_optimization_methods.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_both_limit_angle_inactive(self, optimization_variables_avg):
        l, Q, A = optimization_variables_avg
        Q_in = 5e1 * np.eye(len(l))
        max_angle = 45
        m = 2e-3
        m1 = 4e-3
        f = .02
        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 #3
Source File: GFK.py    From transferlearning with MIT License 6 votes vote down vote up
def subspace_disagreement_measure(self, Ps, Pt, Pst):
        """
        Get the best value for the number of subspaces
        For more details, read section 3.4 of the paper.
        **Parameters**
          Ps: Source subspace
          Pt: Target subspace
          Pst: Source + Target subspace
        """

        def compute_angles(A, B):
            _, S, _ = np.linalg.svd(np.dot(A.T, B))
            S[np.where(np.isclose(S, 1, atol=self.eps) == True)[0]] = 1
            return np.arccos(S)

        max_d = min(Ps.shape[1], Pt.shape[1], Pst.shape[1])
        alpha_d = compute_angles(Ps, Pst)
        beta_d = compute_angles(Pt, Pst)
        d = 0.5 * (np.sin(alpha_d) + np.sin(beta_d))
        return np.argmax(d) 
Example #4
Source File: rotation_utils.py    From DOTA_models with Apache License 2.0 6 votes vote down vote up
def rotate_camera_to_point_at(up_from, lookat_from, up_to, lookat_to):
  inputs = [up_from, lookat_from, up_to, lookat_to]
  for i in range(4):
    inputs[i] = normalize(np.array(inputs[i]).reshape((-1,)))
  up_from, lookat_from, up_to, lookat_to = inputs
  r1 = r_between(lookat_from, lookat_to)

  new_x = np.dot(r1, np.array([1, 0, 0]).reshape((-1, 1))).reshape((-1))
  to_x = normalize(np.cross(lookat_to, up_to))
  angle = np.arccos(np.dot(new_x, to_x))
  if angle > ANGLE_EPS:
    if angle < np.pi - ANGLE_EPS:
      ax = normalize(np.cross(new_x, to_x))
      flip = np.dot(lookat_to, ax)
      if flip > 0:
        r2 = get_r_matrix(lookat_to, angle)
      elif flip < 0:
        r2 = get_r_matrix(lookat_to, -1. * angle)
    else:
      # Angle of rotation is too close to 180 degrees, direction of rotation
      # does not matter.
      r2 = get_r_matrix(lookat_to, angle)
  else:
    r2 = np.eye(3)
  return np.dot(r2, r1) 
Example #5
Source File: test_umath.py    From recruit with Apache License 2.0 6 votes vote down vote up
def test_against_cmath(self):
        import cmath

        points = [-1-1j, -1+1j, +1-1j, +1+1j]
        name_map = {'arcsin': 'asin', 'arccos': 'acos', 'arctan': 'atan',
                    'arcsinh': 'asinh', 'arccosh': 'acosh', 'arctanh': 'atanh'}
        atol = 4*np.finfo(complex).eps
        for func in self.funcs:
            fname = func.__name__.split('.')[-1]
            cname = name_map.get(fname, fname)
            try:
                cfunc = getattr(cmath, cname)
            except AttributeError:
                continue
            for p in points:
                a = complex(func(np.complex_(p)))
                b = cfunc(p)
                assert_(abs(a - b) < atol, "%s %s: %s; cmath: %s" % (fname, p, a, b)) 
Example #6
Source File: test_umath.py    From recruit with Apache License 2.0 6 votes vote down vote up
def test_branch_cuts_complex64(self):
        # check branch cuts and continuity on them
        _check_branch_cut(np.log,   -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log2,  -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log10, -0.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.log1p, -1.5, 1j, 1, -1, True, np.complex64)
        _check_branch_cut(np.sqrt,  -0.5, 1j, 1, -1, True, np.complex64)

        _check_branch_cut(np.arcsin, [ -2, 2],   [1j, 1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arccos, [ -2, 2],   [1j, 1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arctan, [0-2j, 2j],  [1,  1], -1, 1, True, np.complex64)

        _check_branch_cut(np.arcsinh, [0-2j,  2j], [1,   1], -1, 1, True, np.complex64)
        _check_branch_cut(np.arccosh, [ -1, 0.5], [1j,  1j], 1, -1, True, np.complex64)
        _check_branch_cut(np.arctanh, [ -2,   2], [1j, 1j], 1, -1, True, np.complex64)

        # check against bogus branch cuts: assert continuity between quadrants
        _check_branch_cut(np.arcsin, [0-2j, 2j], [ 1,  1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arccos, [0-2j, 2j], [ 1,  1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arctan, [ -2,  2], [1j, 1j], 1, 1, False, np.complex64)

        _check_branch_cut(np.arcsinh, [ -2,  2, 0], [1j, 1j, 1], 1, 1, False, np.complex64)
        _check_branch_cut(np.arccosh, [0-2j, 2j, 2], [1,  1,  1j], 1, 1, False, np.complex64)
        _check_branch_cut(np.arctanh, [0-2j, 2j, 0], [1,  1,  1j], 1, 1, False, np.complex64) 
Example #7
Source File: core.py    From neuropythy with GNU Affero General Public License v3.0 6 votes vote down vote up
def arccosine(x, null=(-np.inf, np.inf)):
    '''
    arccosine(x) is equivalent to acos(x) except that it also works on sparse arrays.

    The optional argument null (default, (-numpy.inf, numpy.inf)) may be specified to indicate what
    value(s) should be assigned when x < -1 or x > 1. If only one number is given, then it is used
    for both values; otherwise the first value corresponds to <-1 and the second to >1.  If null is
    None, then an error is raised when invalid values are encountered.
    '''
    if sps.issparse(x): x = x.toarray()
    else:               x = np.asarray(x)
    try:    (nln,nlp) = null
    except Exception: (nln,nlp) = (null,null)
    ii = None if nln is None else np.where(x < -1)
    jj = None if nlp is None else np.where(x > 1)
    if ii: x[ii] = 0
    if jj: x[jj] = 0
    x = np.arccos(x)
    if ii: x[ii] = nln
    if jj: x[jj] = nlp
    return x 
Example #8
Source File: test_umath.py    From recruit with Apache License 2.0 6 votes vote down vote up
def test_branch_cuts(self):
        # check branch cuts and continuity on them
        _check_branch_cut(np.log,   -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log2,  -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log10, -0.5, 1j, 1, -1, True)
        _check_branch_cut(np.log1p, -1.5, 1j, 1, -1, True)
        _check_branch_cut(np.sqrt,  -0.5, 1j, 1, -1, True)

        _check_branch_cut(np.arcsin, [ -2, 2],   [1j, 1j], 1, -1, True)
        _check_branch_cut(np.arccos, [ -2, 2],   [1j, 1j], 1, -1, True)
        _check_branch_cut(np.arctan, [0-2j, 2j],  [1,  1], -1, 1, True)

        _check_branch_cut(np.arcsinh, [0-2j,  2j], [1,   1], -1, 1, True)
        _check_branch_cut(np.arccosh, [ -1, 0.5], [1j,  1j], 1, -1, True)
        _check_branch_cut(np.arctanh, [ -2,   2], [1j, 1j], 1, -1, True)

        # check against bogus branch cuts: assert continuity between quadrants
        _check_branch_cut(np.arcsin, [0-2j, 2j], [ 1,  1], 1, 1)
        _check_branch_cut(np.arccos, [0-2j, 2j], [ 1,  1], 1, 1)
        _check_branch_cut(np.arctan, [ -2,  2], [1j, 1j], 1, 1)

        _check_branch_cut(np.arcsinh, [ -2,  2, 0], [1j, 1j, 1], 1, 1)
        _check_branch_cut(np.arccosh, [0-2j, 2j, 2], [1,  1,  1j], 1, 1)
        _check_branch_cut(np.arctanh, [0-2j, 2j, 0], [1,  1,  1j], 1, 1) 
Example #9
Source File: FakeCatalog.py    From EXOSIMS with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def inverse_method(self,N,d):
        
        t = np.linspace(1e-3,0.999,N)
        f = np.log( t / (1 - t) )
        f = f/f[0]
        
        psi= np.pi*f
        cosPsi = np.cos(psi)
        sinTheta = ( np.abs(cosPsi) + (1-np.abs(cosPsi))*np.random.rand(len(cosPsi)))
        
        theta = np.arcsin(sinTheta)
        theta = np.pi-theta + (2*theta - np.pi)*np.round(np.random.rand(len(t)))
        cosPhi = cosPsi/sinTheta
        phi = np.arccos(cosPhi)*(-1)**np.round(np.random.rand(len(t)))
        
        coords = SkyCoord(phi*u.rad,(np.pi/2-theta)*u.rad,d*np.ones(len(phi))*u.pc)

        return coords 
Example #10
Source File: geometry.py    From connecting_the_dots with MIT License 6 votes vote down vote up
def quat_slerp_space(q0, q1, num=100, endpoint=True):
  q0 = q0.ravel()
  q1 = q1.ravel()
  dot = q0.dot(q1)
  if dot < 0:
    q1 *= -1
    dot *= -1
  t = np.linspace(0, 1, num=num, endpoint=endpoint, dtype=q0.dtype)
  t = t.reshape((-1,1))
  if dot > 0.9995:
    ret = q0 + t * (q1 - q0)
    return ret
  dot = np.clip(dot, -1, 1)
  theta0 = np.arccos(dot)
  theta = theta0 * t
  s0 = np.cos(theta) - dot * np.sin(theta) / np.sin(theta0)
  s1 = np.sin(theta) / np.sin(theta0)
  return (s0 * q0) + (s1 * q1) 
Example #11
Source File: geometry.py    From connecting_the_dots with MIT License 6 votes vote down vote up
def axisangle_from_rotm(R):
  # logarithm of rotation matrix
  # R = R.reshape(-1,3,3)
  # tr = np.trace(R, axis1=1, axis2=2)
  # phi = np.arccos(np.clip((tr - 1) / 2, -1, 1))
  # scale = np.zeros_like(phi)
  # div = 2 * np.sin(phi)
  # np.divide(phi, div, out=scale, where=np.abs(div) > 1e-6)
  # A = (R - R.transpose(0,2,1)) * scale.reshape(-1,1,1)
  # aa = np.stack((A[:,2,1], A[:,0,2], A[:,1,0]), axis=1)
  # return aa.squeeze()
  R = R.reshape(-1,3,3)
  omega = np.empty((R.shape[0], 3), dtype=R.dtype)
  omega[:,0] = R[:,2,1] - R[:,1,2]
  omega[:,1] = R[:,0,2] - R[:,2,0]
  omega[:,2] = R[:,1,0] - R[:,0,1]
  r = np.linalg.norm(omega, axis=1).reshape(-1,1)
  t = np.trace(R, axis1=1, axis2=2).reshape(-1,1)
  omega = np.arctan2(r, t-1) * omega
  aa = np.zeros_like(omega)
  np.divide(omega, r, out=aa, where=r != 0)
  return aa.squeeze() 
Example #12
Source File: reportables.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def model_model_angles_btwn_axes(A, B, mxBasis):  # Note: default 'gm' basis
    """ Angle between the rotation axes of A and B (1-qubit gates)"""
    decomp = _tools.decompose_gate_matrix(A)
    decomp2 = _tools.decompose_gate_matrix(B)
    axisOfRotn = decomp.get('axis of rotation', None)
    rotnAngle = decomp.get('pi rotations', 'X')
    axisOfRotn2 = decomp2.get('axis of rotation', None)
    rotnAngle2 = decomp2.get('pi rotations', 'X')

    if rotnAngle == 'X' or abs(rotnAngle) < 1e-4 or \
       rotnAngle2 == 'X' or abs(rotnAngle2) < 1e-4:
        return _np.nan

    if axisOfRotn is None or axisOfRotn2 is None:
        return _np.nan

    real_dot = _np.clip(_np.real(_np.dot(axisOfRotn, axisOfRotn2)), -1.0, 1.0)
    return _np.arccos(abs(real_dot)) / _np.pi
    #Note: abs() allows axis to be off by 180 degrees -- if showing *angle* as
    #      well, must flip sign of angle of rotation if you allow axis to
    #      "reverse" by 180 degrees. 
Example #13
Source File: utils.py    From iGAN with MIT License 6 votes vote down vote up
def interp_z(z0, z1, ratio, interp='linear'):
    if interp == 'linear':
        z_t = (1 - ratio) * z0 + ratio * z1

    if interp == 'slerp':
        N = len(z0)
        z_t = []
        for i in range(N):
            z0_i = z0[i]
            z1_i = z1[i]
            z0_n = z0_i / np.linalg.norm(z0_i)
            z1_n = z1_i / np.linalg.norm(z1_i)
            omega = np.arccos(np.dot(z0_n, z1_n))
            sin_omega = np.sin(omega)
            if sin_omega == 0:
                z_i = interp_z(z0_i, z1_i, ratio, 'linear')
            else:
                z_i = np.sin((1 - ratio) * omega) / sin_omega * z0_i + np.sin(ratio * omega) / sin_omega * z1_i
            z_t.append(z_i[np.newaxis, ...])
        z_t = np.concatenate(z_t, axis=0)
    return z_t 
Example #14
Source File: test_optimization_methods.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_both_limit_angle_limit(self, optimization_variables_avg):
        l, Q, A = optimization_variables_avg
        Q_in = 5e1 * np.eye(len(l))
        max_angle = 25
        m = 2e-3
        m1 = 4e-3
        f = .02
        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 #15
Source File: test_optimization_methods.py    From simnibs with 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 #16
Source File: test_optimization_methods.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_both_limit_angle_infeasible_field(self, optimization_variables_avg_QCQP):
        l, Q, A, Q_in = optimization_variables_avg_QCQP
        max_angle = 15
        m = 2e-3
        m1 = 4e-3
        f = 2
        x = optimization_methods.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.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle) 
Example #17
Source File: fisheye.py    From DualFisheye with MIT License 6 votes vote down vote up
def get_uv(self, xyz_vec):
        # Extract lens parameters of interest.
        fov_rad = self.lens.fov_deg * pi / 180
        fov_scale = np.float32(2 * self.lens.radius_px / fov_rad)
        # Normalize the input vector and rotate to match lens reference axes.
        xyz_rot = get_rotation_matrix(self.lens.center_qq) * matrix_norm(xyz_vec)
        # Convert to polar coordinates relative to lens boresight.
        # (In lens coordinates, unit vector's X axis gives boresight angle;
        #  normalize Y/Z to get a planar unit vector for the bearing.)
        # Note: Image +Y maps to 3D +Y, and image +X maps to 3D +Z.
        theta_rad = np.arccos(xyz_rot[0,:])
        proj_vec = matrix_norm(np.concatenate((xyz_rot[2,:], xyz_rot[1,:])))
        # Fisheye lens maps 3D angle to focal-plane radius.
        # TODO: Do we need a better model for lens distortion?
        rad_px = theta_rad * fov_scale
        # Convert back to focal-plane rectangular coordinates.
        uv = np.multiply(rad_px, proj_vec) + self.lens.center_px
        return np.asarray(uv + 0.5, dtype=int)

    # Given an 2xN array of UV pixel coordinates, check if each pixel is
    # within the fisheye field of view. Returns N-element boolean mask. 
Example #18
Source File: test_optimization_methods.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_limit_nr_angle(seld, optimization_variables_avg_QCQP):
        l, Q, A, Q_in = optimization_variables_avg_QCQP
        max_angle = 15
        m = 2e-3
        m1 = 4e-3
        f = 2
        n = 4
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle,
                                                   max_active_electrodes=n)

        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.linalg.norm(x, 0) <= n
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle) 
Example #19
Source File: test_optimization_methods.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def test_limit_nr_angle_change_Q(seld, optimization_variables_avg_QCQP):
        l, Q, A, Q_in = optimization_variables_avg_QCQP
        max_angle = 15
        m = 2e-3
        m1 = 4e-3
        f = .01
        n = 4
        x = optimization_methods.optimize_focality(l, Q, f,
                                                   max_el_current=m,
                                                   max_total_current=m1,
                                                   Qin=Q_in,
                                                   max_angle=max_angle,
                                                   max_active_electrodes=n)

        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.linalg.norm(x, 0) <= n
        assert np.arccos(l.dot(x) / np.sqrt(x.dot(Q_in).dot(x))) <= np.deg2rad(max_angle) 
Example #20
Source File: test_geom.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def make12(b):
    theta1 = numpy.arccos(1/numpy.sqrt(5))
    theta2 = (numpy.pi - theta1) * .5
    r = b/2./numpy.sin(theta1/2)
    rot72 = rotmatz(numpy.pi*2/5)
    s1 = numpy.sin(theta1)
    c1 = numpy.cos(theta1)
    p1 = numpy.array(( s1*r,  0,  c1*r))
    p2 = numpy.array((-s1*r,  0, -c1*r))
    coord = [(  0,  0,    r)]
    for i in range(5):
        coord.append(p1)
        p1 = numpy.dot(p1, rot72)
    for i in range(5):
        coord.append(p2)
        p2 = numpy.dot(p2, rot72)
    coord.append((  0,  0,  -r))
    return numpy.array(coord) 
Example #21
Source File: rpetools.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def extract_theta(model, rpeconfig_inst):
    """
    For a given model, obtain the angle between the estimated "loose axis" and
    the target "loose axis".

    WARNING:  This is a gauge-covariant parameter!  (I think!)  Gauge must be
    fixed prior to estimating.

    Parameters
    ----------
    model : Model
        The model whose loose axis misalignment is to be calculated.

    rpeconfig_inst : rpeconfig object
        Declares which model configuration RPE should be trying to fit;
        determines particular functions and values to be used.

    Returns
    -------
    thetaVal : float
        The value of theta for the input model.
    """
    op_label = rpeconfig_inst.loose_axis_gate_label
    decomp = _decompose_gate_matrix(model.operations[op_label])
    target_axis = rpeconfig_inst.loose_axis_target

    decomp = _decompose_gate_matrix(model.operations[op_label])
    thetaVal = _np.real_if_close([_np.arccos(
        _np.dot(decomp['axis of rotation'], target_axis))])[0]
    if thetaVal > _np.pi / 2:
        thetaVal = _np.pi - thetaVal
    elif thetaVal < -_np.pi / 2:
        thetaVal = _np.pi + thetaVal
    return thetaVal 
Example #22
Source File: reportables.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def angles_btwn_rotn_axes(model):
    """
    Array of angles between the rotation axes of the gates of `model`.

    Returns
    -------
    numpy.ndarray
        Of size `(nOperations,nGate)` where `nOperations=len(model.operations)`
    """
    opLabels = list(model.operations.keys())
    angles_btwn_rotn_axes = _np.zeros((len(opLabels), len(opLabels)), 'd')

    for i, gl in enumerate(opLabels):
        decomp = _tools.decompose_gate_matrix(model.operations[gl].todense())
        rotnAngle = decomp.get('pi rotations', 'X')
        axisOfRotn = decomp.get('axis of rotation', None)

        for j, gl_other in enumerate(opLabels[i + 1:], start=i + 1):
            decomp_other = _tools.decompose_gate_matrix(model.operations[gl_other])
            rotnAngle_other = decomp_other.get('pi rotations', 'X')

            if str(rotnAngle) == 'X' or abs(rotnAngle) < 1e-4 or \
               str(rotnAngle_other) == 'X' or abs(rotnAngle_other) < 1e-4:
                angles_btwn_rotn_axes[i, j] = _np.nan
            else:
                axisOfRotn_other = decomp_other.get('axis of rotation', None)
                if axisOfRotn is not None and axisOfRotn_other is not None:
                    real_dot = _np.clip(_np.real(_np.dot(axisOfRotn, axisOfRotn_other)), -1.0, 1.0)
                    angles_btwn_rotn_axes[i, j] = _np.arccos(real_dot) / _np.pi
                else:
                    angles_btwn_rotn_axes[i, j] = _np.nan

            angles_btwn_rotn_axes[j, i] = angles_btwn_rotn_axes[i, j]
    return angles_btwn_rotn_axes 
Example #23
Source File: functions.py    From tangent with Apache License 2.0 5 votes vote down vote up
def numpy_arccos(a):
  return np.arccos(a) 
Example #24
Source File: eval_ecp.py    From pyqmc with MIT License 5 votes vote down vote up
def get_rot(nconf, naip):
    """
    Returns the integration weights (naip), and the positions of the rotated electron e (nconf x naip x 3)
    Parameters: 
      configs[:,e,:]: epos of the electron e to be rotated
    Returns:
      weights: naip array
      epos_rot: positions of the rotated electron, nconf x naip x 3
      
    """
    # t and p are sampled randomly over a sphere around the atom
    t = np.random.uniform(low=0.0, high=np.pi, size=nconf)
    p = np.random.uniform(low=0.0, high=2 * np.pi, size=nconf)

    def sphere(t_, p_):
        s = np.sin(t_)
        return s * np.cos(p_), s * np.sin(p_), np.cos(t_)

    # rotated unit vectors:
    rot = np.zeros([3, 3, nconf])
    rot[0, :, :] = sphere(np.zeros(nconf) + np.pi / 2.0, p - np.pi / 2.0)
    rot[1, :, :] = sphere(t + np.pi / 2.0, p)
    rot[2, :, :] = sphere(t, p)

    if naip == 6:
        d1 = np.array([0.0, 1.0, 0.5, 0.5, 0.5, 0.5]) * np.pi
        d2 = np.array([0.0, 0.0, 0.0, 0.5, 1.0, 1.5]) * np.pi
    elif naip == 12:
        tha = np.arccos(1.0 / np.sqrt(5.0))
        d1 = np.array([0, np.pi] + [tha, np.pi - tha] * 5)
        d2 = np.array([0, 0] + list(range(10))) * np.pi / 5

    rot_vec = np.einsum("ilj,ik->jkl", rot, sphere(d1, d2))
    weights = 1.0 / naip * np.ones(naip)

    return weights, rot_vec 
Example #25
Source File: sasa.py    From pytim with GNU General Public License v3.0 5 votes vote down vote up
def _overlap(self, Ri, Rj, pij, dzi):
        dzj = pij[:, 2] - dzi
        ri2 = Ri**2 - dzi**2
        rj2 = Rj**2 - dzj**2

        cond = np.where(rj2 >= 0.0)[0]
        if len(cond) == 0:
            return [], []
        rj2 = rj2[cond]
        pij = pij[cond]

        ri, rj = np.sqrt(ri2), np.sqrt(rj2)
        pij2 = pij**2
        dij2 = pij2[:, 0] + pij2[:, 1]
        dij = np.sqrt(dij2)

        # slice within neighboring one
        if np.any(dij <= rj - ri):
            return np.array([2. * np.pi]), np.array([0.0])

        # c1 => no overlap; c2 => neighboring slice enclosed
        c1, c2 = dij < ri + rj, dij > ri - rj
        cond = np.where(np.logical_and(c1, c2))[0]
        if len(cond) == 0:
            return [], []

        arg = (ri2 + dij2 - rj2)[cond] / (ri * dij * 2.)[cond]
        alpha = 2. * np.arccos(arg)
        argx, argy = pij[:, 0], pij[:, 1]
        beta = np.arctan2(argx[cond], argy[cond])
        return alpha, beta 
Example #26
Source File: mesh_io.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def triangle_angles(self):
        ''' Calculates minimum angle of the mesh triangles
        
        Returns
        --------
        min_triangle_angles: ElementData
            ElementData field with the angles for each triangle, in degrees
        '''
        tr_indexes = self.elm.triangles
        node_tr = self.nodes[self.elm[tr_indexes, :3]]
        a = np.linalg.norm(node_tr[:, 1] - node_tr[:, 0], axis=1)
        b = np.linalg.norm(node_tr[:, 2] - node_tr[:, 0], axis=1)
        c = np.linalg.norm(node_tr[:, 2] - node_tr[:, 1], axis=1)
        cos_angles = np.vstack([
            (a**2 + b**2 - c**2) / (2*a*b),
            (a**2 + c**2 - b**2) / (2*a*c),
            (b**2 + c**2 - a**2) / (2*b*c),
        ]).T
        angles = np.rad2deg(np.arccos(cos_angles))

        angle_field  = ElementData(
            np.nan*np.zeros((self.elm.nr, 3)),
            'min_triangle_angles'
        )
        angle_field[tr_indexes] = angles
        return angle_field 
Example #27
Source File: atlas.py    From ibllib with MIT License 5 votes vote down vote up
def cart2sph(x, y, z):
    """
    Converts cartesian to spherical Coordinates
    theta: polar angle, phi: azimuth
    """
    r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
    phi = np.arctan2(y, x) * 180 / np.pi
    theta = np.zeros_like(r)
    iok = r != 0
    theta[iok] = np.arccos(z[iok] / r[iok]) * 180 / np.pi
    return r, theta, phi 
Example #28
Source File: grads.py    From tangent with Apache License 2.0 5 votes vote down vote up
def arccos(y, x):
  d[x] = -d[y] / numpy.sqrt(1.0 - x * x) 
Example #29
Source File: sin.py    From cgpm with Apache License 2.0 5 votes vote down vote up
def _valid_x(self, y):
        """Compute valid regions of x for y \in [0, 1], with overflow."""
        assert 0<=y<=1
        x_max = np.arccos(y)
        if y+self.noise < 1:
            x_min = np.arccos(y+self.noise)
        else:
            x_min = 0
        # compute overflow
        overflow = 0
        if y < self.noise:
            overflow = np.arccos(y-self.noise) - np.pi / 2
        return x_max - x_min, overflow 
Example #30
Source File: evaluation.py    From hierarchical_loc with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def is_gt_match_3D(query_poses, ref_poses, distance_thresh, angle_thresh):
    distances = np.linalg.norm(np.expand_dims(query_poses['pos'], axis=1)
                               - np.expand_dims(ref_poses['pos'], axis=0), axis=-1)
    angle_errors = np.arccos(
        (np.trace(
            np.matmul(np.expand_dims(np.linalg.inv(query_poses['rot']), axis=1),
                      np.expand_dims(ref_poses['rot'], axis=0)),
            axis1=2, axis2=3) - 1)/2)
    return np.logical_and(distances < distance_thresh, angle_errors < angle_thresh)