Python numpy.pi() Examples

The following are code examples for showing how to use numpy.pi(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: b2ac   Author: hbldh   File: ellipse.py    MIT License 9 votes vote down vote up
def polygonize(self, n=73):
        """Gets a approximate polygon array representing the ellipse.

        Note that the last point is the same as the first point, creating a closed
        polygon.

        :param n: The number of points to generate. Default is 73 (one vertex every 5 degrees).
        :type n: int
        :return: An [n  x 2] numpy array, describing the boundary vertices of
                 the polygonized ellipse.
        :rtype: :py:class:`numpy.ndarray`

        """
        t = np.linspace(0, 2 * np.pi, num=n, endpoint=True)
        out = np.zeros((len(t), 2), dtype='float')
        out[:, 0] = (self.center_point[0] +
                     self.radii[0] * np.cos(t) * np.cos(self.rotation_angle) -
                     self.radii[1] * np.sin(t) * np.sin(self.rotation_angle))
        out[:, 1] = (self.center_point[1] +
                     self.radii[0] * np.cos(t) * np.sin(self.rotation_angle) +
                     self.radii[1] * np.sin(t) * np.cos(self.rotation_angle))
        return out 
Example 2
Project: FRIDA   Author: LCAV   File: doa.py    MIT License 7 votes vote down vote up
def compute_mode(self):
        """
        Pre-compute mode vectors from candidate locations (in spherical 
        coordinates).
        """
        if self.num_loc is None:
            raise ValueError('Lookup table appears to be empty. \
                Run build_lookup().')
        self.mode_vec = np.zeros((self.max_bin,self.M,self.num_loc), 
            dtype='complex64')
        if (self.nfft % 2 == 1):
            raise ValueError('Signal length must be even.')
        f = 1.0 / self.nfft * np.linspace(0, self.nfft / 2, self.max_bin) \
            * 1j * 2 * np.pi
        for i in range(self.num_loc):
            p_s = self.loc[:, i]
            for m in range(self.M):
                p_m = self.L[:, m]
                if (self.mode == 'near'):
                    dist = np.linalg.norm(p_m - p_s, axis=1)
                if (self.mode == 'far'):
                    dist = np.dot(p_s, p_m)
                # tau = np.round(self.fs*dist/self.c) # discrete - jagged
                tau = self.fs * dist / self.c  # "continuous" - smoother
                self.mode_vec[:, m, i] = np.exp(f * tau) 
Example 3
Project: skylab   Author: coenders   File: ps_model.py    GNU General Public License v3.0 7 votes vote down vote up
def background(self, ev):
        r"""Spatial background distribution.

        For IceCube is only declination dependent, in a more general scenario,
        it is dependent on zenith and
        azimuth, e.g. in ANTARES, KM3NET, or using time dependent information.

        Parameters
        -----------
        ev : structured array
            Event array, importand information *sinDec* for this calculation

        Returns
        --------
        P : array-like
            spatial background probability for each event to be found
            at *sinDec*

        """
        return 1. / 2. / np.pi * np.exp(self.bckg_spline(ev["sinDec"])) 
Example 4
Project: UR5_Controller   Author: tsinghua-rll   File: quaternion.py    MIT License 6 votes vote down vote up
def to_euler(q):
    # rpy
    sinr = 2.0 * (q[0] * q[1] + q[2] * q[3])
    cosr = 1.0 - 2.0 * (q[1] * q[1] + q[2] * q[2])
    roll = math.atan2(sinr, cosr)

    sinp = 2.0 * (q[0] * q[2] - q[3] * q[1])
    if math.fabs(sinp) >= 1.:
        pitch = math.copysign(np.pi / 2., sinp)
    else:
        pitch = math.asin(sinp)

    siny = 2.0 * (q[0] * q[3] + q[1] * q[2])
    cosy = 1.0 - 2.0 * (q[2] * q[2] + q[3] * q[3])
    yaw = math.atan2(siny, cosy)

    return np.asarray((roll, pitch, yaw), np.float32) 
Example 5
Project: cs207-FinalProject   Author: PYNE-AD   File: elemFunctions_Dual_test.py    MIT License 6 votes vote down vote up
def test_tan_ad_results():
	# Defined Realue and Dualivative when cos(Real)!=0
	# Positive reals
	x = Dual(0.5, 2.0)
	f = ef.tan(x)
	assert f.Real == np.array([[np.tan(0.5)]])
	assert f.Dual == np.array([[2.0 / (np.cos(0.5)**2)]])
	# Negative reals
	y = Dual(-0.5, 2.0)
	f = ef.tan(y)
	assert f.Real == np.array([[np.tan(-0.5)]])
	assert f.Dual == np.array([[2.0 / (np.cos(-0.5)**2)]])
	# Zero
	z = Dual(0.0, 2.0)
	f = ef.tan(z)
	assert f.Real == np.array([[np.tan(0)]])
	assert f.Dual == np.array([[2.0]])

	# Undefined Value and Derivative when cos(Real)==0
	with pytest.warns(RuntimeWarning):
		h = Dual(np.pi/2, 1.0)
		f = ef.tan(h)
		print("HERE",f.Real)
		assert np.isnan(f.Real)
		assert np.isnan(f.Dual) 
Example 6
Project: aospy   Author: spencerahill   File: vertcoord.py    Apache License 2.0 6 votes vote down vote up
def to_radians(arr, is_delta=False):
    """Force data with units either degrees or radians to be radians."""
    # Infer the units from embedded metadata, if it's there.
    try:
        units = arr.units
    except AttributeError:
        pass
    else:
        if units.lower().startswith('degrees'):
            warn_msg = ("Conversion applied: degrees -> radians to array: "
                        "{}".format(arr))
            logging.debug(warn_msg)
            return np.deg2rad(arr)
    # Otherwise, assume degrees if the values are sufficiently large.
    threshold = 0.1*np.pi if is_delta else 4*np.pi
    if np.max(np.abs(arr)) > threshold:
        warn_msg = ("Conversion applied: degrees -> radians to array: "
                    "{}".format(arr))
        logging.debug(warn_msg)
        return np.deg2rad(arr)
    return arr 
Example 7
Project: xrft   Author: xgcm   File: test_xrft.py    MIT License 6 votes vote down vote up
def test_cross_phase_2d(self, dask):
        Ny, Nx = (32, 16)
        x = np.linspace(0, 1, num=Nx, endpoint=False)
        y = np.ones(Ny)
        f = 6
        phase_offset = np.pi/2
        signal1 = np.cos(2*np.pi*f*x)  # frequency = 1/(2*pi)
        signal2 = np.cos(2*np.pi*f*x - phase_offset)
        da1 = xr.DataArray(data=signal1*y[:,np.newaxis], name='a',
                          dims=['y','x'], coords={'y':y, 'x':x})
        da2 = xr.DataArray(data=signal2*y[:,np.newaxis], name='b',
                          dims=['y','x'], coords={'y':y, 'x':x})
        with pytest.raises(ValueError):
            xrft.cross_phase(da1, da2, dim=['y','x'])

        if dask:
            da1 = da1.chunk({'x': 16})
            da2 = da2.chunk({'x': 16})
        cp = xrft.cross_phase(da1, da2, dim=['x'])
        actual_phase_offset = cp.sel(freq_x=f).values
        npt.assert_almost_equal(actual_phase_offset, phase_offset) 
Example 8
Project: FRIDA   Author: LCAV   File: point_cloud.py    MIT License 6 votes vote down vote up
def doa(self, receiver, source):
        ''' Computes the direction of arrival wrt a source and receiver '''

        s_ind = self.key2ind(source)
        r_ind = self.key2ind(receiver)

        # vector from receiver to source
        v = self.X[:,s_ind] - self.X[:,r_ind]

        azimuth = np.arctan2(v[1], v[0])
        elevation = np.arctan2(v[2], la.norm(v[:2]))

        azimuth = azimuth + 2*np.pi if azimuth < 0. else azimuth
        elevation = elevation + 2*np.pi if elevation < 0. else elevation

        return np.array([azimuth, elevation]) 
Example 9
Project: DataHack2018   Author: InnovizTech   File: vis.py    BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def draw_cuboids(self, boxes):
        for box_idx, box in enumerate(boxes):
            color = box['color'] if hasattr(box, 'color') else np.ones(4)
            size = box['size']
            translation = box['translation']
            rotation = box['rotation'] / np.pi * 180.
            try:
                text = box['text']
            except:
                text = ''

            if box_idx < len(self._cuboids):
                self._cuboids[box_idx].show()
                self._cuboids[box_idx].update_values(size, translation, rotation, color, text)
            else:
                self._cuboids.append(Cuboid(size, translation, rotation, color, text))
        for c in self._cuboids[len(boxes):]:
            c.hide() 
Example 10
Project: DataHack2018   Author: InnovizTech   File: vis_utils.py    BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def _zoom(self, value):
        if type(value) is str:
            if value == '-z':
                value = -base.camera.getPos()[2]
            if value == '+z':
                value = base.camera.getPos()[2]
        scale = 0.3
        hpr = base.camera.getHpr()
        rotation_angles = np.array((-hpr[1], 0., -hpr[0])) / 180. * np.pi
        rotation = euler_matrix(rotation_angles[2], rotation_angles[1], rotation_angles[0])
        # inverse
        R = rotation.T
        offset = R.dot(np.array([[0., 1., 0.]]).T)

        new_pos = np.array(base.camera.getPos()) + offset.ravel() * scale * value
        new_pos[2] = max(new_pos[2], 0.)
        base.camera.setPos(new_pos[0], new_pos[1], new_pos[2])
        self._update_location_text() 
Example 11
Project: ddd-utils   Author: inconvergent   File: random.py    MIT License 6 votes vote down vote up
def random_points_in_circle(n,xx,yy,rr):
  """
  get n random points in a circle.
  """


  rnd = random(size=(n,3))
  t = 2.*PI*rnd[:,0]
  u = rnd[:,1:].sum(axis=1)
  r = zeros(n,'float')
  mask = u>1.
  xmask = logical_not(mask)
  r[mask] = 2.-u[mask]
  r[xmask] = u[xmask]
  xyp = reshape(rr*r,(n,1))*column_stack( (cos(t),sin(t)) )
  dartsxy  = xyp + array([xx,yy])
  return dartsxy 
Example 12
Project: skylab   Author: coenders   File: ps_injector.py    GNU General Public License v3.0 6 votes vote down vote up
def _setup(self):
        r"""Reset solid angle and declination band.

        """
        A, B = self.sinDec_range

        m = (A - B + 2. * self.sinDec_bandwidth) / (A - B)
        b = self.sinDec_bandwidth * (A + B) / (B - A)

        sinDec = m * np.sin(self.src_dec) + b

        min_sinDec = max(A, sinDec - self.sinDec_bandwidth)
        max_sinDec = min(B, sinDec + self.sinDec_bandwidth)

        self._min_dec = np.arcsin(min_sinDec)
        self._max_dec = np.arcsin(max_sinDec)

        # Solid angle of selected events
        self._omega = 2. * np.pi * (max_sinDec - min_sinDec) 
Example 13
Project: skylab   Author: coenders   File: data.py    GNU General Public License v3.0 6 votes vote down vote up
def exp(N=100):
    r"""Create uniformly distributed data on sphere. """
    g = 3.7

    arr = np.empty((N, ), dtype=[("ra", np.float), ("sinDec", np.float),
                                 ("sigma", np.float), ("logE", np.float)])

    arr["ra"] = np.random.uniform(0., 2.*np.pi, N)
    arr["sinDec"] = np.random.uniform(-1., 1., N)

    E = np.log10(np.random.pareto(g, size=N) + 1)
    arr["sigma"] = np.random.lognormal(mean=np.log((mrs - mrs_min) * np.exp(-np.log(10)*E) + mrs_min),
                                       sigma=log_sig)
    arr["logE"] = E + logE_res * np.random.normal(size=N)

    return arr 
Example 14
Project: skylab   Author: coenders   File: data.py    GNU General Public License v3.0 6 votes vote down vote up
def MC(N=1000):
    r"""Create uniformly distributed MC data on sphere. """
    g = 2.

    arr = np.empty((N, ), dtype=[("ra", np.float), ("sinDec", np.float),
                                 ("sigma", np.float), ("logE", np.float),
                                 ("trueRa", np.float), ("trueDec", np.float),
                                 ("trueE", np.float), ("ow", np.float)])

    # true information

    arr["trueRa"] = np.random.uniform(0., 2.*np.pi, N)
    arr["trueDec"] = np.arcsin(np.random.uniform(-1., 1., N))
    arr["trueE"] = np.random.pareto(g, size=N) + 1
    arr["ow"] = arr["trueE"]**(g)
    arr["ow"] /= arr["ow"].sum()

    eta = np.random.uniform(0., 2.*np.pi, len(arr))
    arr["sigma"] = np.random.lognormal(mean=np.log((mrs - mrs_min) * np.exp(-np.log(10)*np.log10(arr["trueE"])) + mrs_min),
                                       sigma=log_sig)
    arr["ra"] = arr["trueRa"] + np.cos(eta) * arr["sigma"] / np.cos(arr["trueDec"])
    arr["sinDec"] = np.sin(arr["trueDec"] + np.sin(eta) * arr["sigma"])
    arr["logE"] = np.log10(arr["trueE"]) + logE_res * np.random.normal(size=len(arr))

    return arr 
Example 15
Project: spqrel_tools   Author: LCAS   File: urdf.py    MIT License 6 votes vote down vote up
def to_openrave_xml(self, doc):
        xml = doc.createElement("translation")
        set_content(doc, xml, self.position)
        elements = [xml]
        if not self.rotation[0] == 0:
            rotr = doc.createElement("rotationaxis")
            set_content(doc, rotr, [1, 0, 0, self.rotation[0] * 180.0/pi])
            elements.append(rotr)
        if not self.rotation[1] == 0:
            rotp = doc.createElement("rotationaxis")
            set_content(doc, rotp, [0, 1, 0, self.rotation[1] * 180.0/pi])
            elements.append(rotp)
        if not self.rotation[2] == 0:
            roty = doc.createElement("rotationaxis")
            set_content(doc, roty, [0, 0, 1, self.rotation[2] * 180.0/pi])
            elements.append(roty)

        return elements 
Example 16
Project: DOTA_models   Author: ringringyi   File: rotation_utils.py    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 17
Project: DOTA_models   Author: ringringyi   File: nav_env.py    Apache License 2.0 6 votes vote down vote up
def get_loc_axis(self, node, delta_theta, perturb=None):
    """Based on the node orientation returns X, and Y axis. Used to sample the
    map in egocentric coordinate frame.
    """
    if type(node) == tuple:
      node = np.array([node])
    if perturb is None:
      perturb = np.zeros((node.shape[0], 4))
    xyt = self.to_actual_xyt_vec(node)
    x = xyt[:,[0]] + perturb[:,[0]]
    y = xyt[:,[1]] + perturb[:,[1]]
    t = xyt[:,[2]] + perturb[:,[2]]
    theta = t*delta_theta
    loc = np.concatenate((x,y), axis=1)
    x_axis = np.concatenate((np.cos(theta), np.sin(theta)), axis=1)
    y_axis = np.concatenate((np.cos(theta+np.pi/2.), np.sin(theta+np.pi/2.)),
                            axis=1)
    # Flip the sampled map where need be.
    y_axis[np.where(perturb[:,3] > 0)[0], :] *= -1.
    return loc, x_axis, y_axis, theta 
Example 18
Project: DOTA_models   Author: ringringyi   File: distributions.py    Apache License 2.0 6 votes vote down vote up
def diag_gaussian_log_likelihood(z, mu=0.0, logvar=0.0):
  """Log-likelihood under a Gaussian distribution with diagonal covariance.
    Returns the log-likelihood for each dimension.  One should sum the
    results for the log-likelihood under the full multidimensional model.

  Args:
    z: The value to compute the log-likelihood.
    mu: The mean of the Gaussian
    logvar: The log variance of the Gaussian.

  Returns:
    The log-likelihood under the Gaussian model.
  """

  return -0.5 * (logvar + np.log(2*np.pi) + \
                 tf.square((z-mu)/tf.exp(0.5*logvar))) 
Example 19
Project: DOTA_models   Author: ringringyi   File: distributions.py    Apache License 2.0 6 votes vote down vote up
def gaussian_pos_log_likelihood(unused_mean, logvar, noise):
  """Gaussian log-likelihood function for a posterior in VAE

  Note: This function is specialized for a posterior distribution, that has the
  form of z = mean + sigma * noise.

  Args:
    unused_mean: ignore
    logvar: The log variance of the distribution
    noise: The noise used in the sampling of the posterior.

  Returns:
    The log-likelihood under the Gaussian model.
  """
  # ln N(z; mean, sigma) = - ln(sigma) - 0.5 ln 2pi - noise^2 / 2
  return - 0.5 * (logvar + np.log(2 * np.pi) + tf.square(noise)) 
Example 20
Project: speed_estimation   Author: NeilNie   File: helper.py    MIT License 5 votes vote down vote up
def optical_flow_rgb(previous, current):

    ''' perform optical flow on two consequtive frames
    and then return the RGB results.

    :param previous: the previous frame (rgb)
    :param current: the current frame (rgb)
    :return: rgb image after optical flow
    '''

    previous_gray = cv2.cvtColor(previous, cv2.COLOR_RGB2GRAY)
    gray = cv2.cvtColor(current, cv2.COLOR_RGB2GRAY)

    flow = cv2.calcOpticalFlowFarneback(previous_gray, gray, None, 0.5, 3, 15, 3, 5, 1.2, 0)

    hsvImg = np.zeros_like(previous)
    hsvImg[..., 1] = 255
    # Obtain the flow magnitude and direction angle
    mag, ang = cv2.cartToPolar(flow[..., 0], flow[..., 1])

    # Update the color image
    hsvImg[..., 0] = 0.5 * ang * 180 / np.pi
    hsvImg[..., 2] = cv2.normalize(mag, None, 0, 255, cv2.NORM_MINMAX)
    rgbImg = cv2.cvtColor(hsvImg, cv2.COLOR_HSV2BGR)
    rgbImg = cv2.resize(rgbImg, (configs.IMG_WIDTH, configs.IMG_HEIGHT))

    return rgbImg 
Example 21
Project: b2ac   Author: hbldh   File: conversion.py    MIT License 5 votes vote down vote up
def conic_to_general_reference(conic_coeffs):
    """Transform from conic section format to general format.

    :param conic_coeffs: The six coefficients defining the ellipse as a conic shape.
    :type conic_coeffs: :py:class:`numpy.ndarray` or tuple
    :param verbose: If debug printout is desired.
    :type verbose: bool
    :return:
    :rtype: tuple

    """
    a, b, c, d, e, f = conic_coeffs

    angle = np.arctan2(b, a - c) / 2

    cos_theta = np.cos(angle)  # np.sqrt((1 + np.cos(2*angle)) / 2)
    sin_theta = np.sin(angle)  # np.sqrt((1 - np.cos(2*angle)) / 2)

    a_prime = a * (cos_theta ** 2) + (b * cos_theta * sin_theta) + c * (sin_theta ** 2)
    c_prime = a * (sin_theta ** 2) - (b * cos_theta * sin_theta) + c * (cos_theta ** 2)
    d_prime = (d * cos_theta) + (e * sin_theta)
    e_prime = (-(d * sin_theta)) + (e * cos_theta)
    f_prime = f

    x_prime = (-d_prime) / (2 * a_prime)
    y_prime = (-e_prime) / (2 * c_prime)
    major_axis = np.sqrt(
        ((-4 * f_prime * a_prime * c_prime) + (c_prime * (d_prime ** 2)) + (a_prime * (e_prime ** 2))) /
        (4 * a_prime * (c_prime ** 2)))
    minor_axis = np.sqrt(
        ((-4 * f_prime * a_prime * c_prime) + (c_prime * (d_prime ** 2)) + (a_prime * (e_prime ** 2))) /
        (4 * (a_prime ** 2) * c_prime))

    if a_prime > c_prime:
        angle += np.pi / 2

    x = x_prime * cos_theta - y_prime * sin_theta
    y = x_prime * sin_theta + y_prime * cos_theta

    return [x, y], [major_axis, minor_axis], angle 
Example 22
Project: b2ac   Author: hbldh   File: ellipse.py    MIT License 5 votes vote down vote up
def get_area(self):
        """Returns the area covered of the B2ACEllipse."""
        return np.pi * self.radii[0] * self.radii[1] 
Example 23
Project: cs207-FinalProject   Author: PYNE-AD   File: elemFunctions_test.py    MIT License 5 votes vote down vote up
def test_tan_ad_results():
	# Defined value and derivative when cos(val)!=0
	# Positive reals
	x = AutoDiff(0.5, 2.0)
	f = ef.tan(x)
	assert f.val == np.array([[np.tan(0.5)]])
	assert f.der == np.array([[2.0 / (np.cos(0.5)**2)]])
	assert f.jacobian == np.array([[1.0 / (np.cos(0.5)**2)]])
	# Negative reals
	y = AutoDiff(-0.5, 2.0)
	f = ef.tan(y)
	assert f.val == np.array([[np.tan(-0.5)]])
	assert f.der == np.array([[2.0 / (np.cos(-0.5)**2)]])
	assert f.jacobian == np.array([[1.0 / (np.cos(-0.5)**2)]])
	# Zero
	z = AutoDiff(0.0, 2.0)
	f = ef.tan(z)
	assert f.val == np.array([[np.tan(0)]])
	assert f.der == np.array([[2.0]])
	assert f.jacobian == np.array([[1.0]])

	# Undefined value and derivative when cos(val)==0
	with pytest.warns(RuntimeWarning):
		h = AutoDiff(np.pi/2, 1.0)
		f = ef.tan(h)
		assert np.isnan(f.val)
		assert np.isnan(f.der)
		assert np.isnan(f.jacobian) 
Example 24
Project: cs207-FinalProject   Author: PYNE-AD   File: elemFunctions_test.py    MIT License 5 votes vote down vote up
def test_tan_constant_results():
	a = ef.tan(5)
	assert a == np.tan(5)
	b = ef.tan(-5)
	assert b == np.tan(-5)
	c = ef.tan(0)
	assert c == np.tan(0)
	# Value undefined when cos(val)==0
	with pytest.warns(RuntimeWarning):
		d = ef.tan(np.pi/2)
		assert np.isnan(d) 
Example 25
Project: cs207-FinalProject   Author: PYNE-AD   File: elemFunctions_Dual_test.py    MIT License 5 votes vote down vote up
def test_tan_constant_results():
	a = ef.tan(5)
	assert a == np.tan(5)
	b = ef.tan(-5)
	assert b == np.tan(-5)
	c = ef.tan(0)
	assert c == np.tan(0)
	# Realue undefined when cos(Real)==0
	with pytest.warns(RuntimeWarning):
		d = ef.tan(np.pi/2)
		assert np.isnan(d) 
Example 26
Project: xrft   Author: xgcm   File: test_xrft.py    MIT License 5 votes vote down vote up
def test_cross_phase_1d(self, dask):
        N = 32
        x = np.linspace(0, 1, num=N, endpoint=False)
        f = 6
        phase_offset = np.pi/2
        signal1 = np.cos(2*np.pi*f*x)  # frequency = 1/(2*pi)
        signal2 = np.cos(2*np.pi*f*x - phase_offset)
        da1 = xr.DataArray(data=signal1, name='a', dims=['x'], coords={'x': x})
        da2 = xr.DataArray(data=signal2, name='b', dims=['x'], coords={'x': x})

        if dask:
            da1 = da1.chunk({'x': 32})
            da2 = da2.chunk({'x': 32})
        cp = xrft.cross_phase(da1, da2, dim=['x'])

        actual_phase_offset = cp.sel(freq_x=f).values
        npt.assert_almost_equal(actual_phase_offset, phase_offset)
        assert cp.name == 'a_b_phase'

        xrt.assert_equal(xrft.cross_phase(da1, da2), cp)

        with pytest.raises(ValueError):
            xrft.cross_phase(da1, da2.isel(x=0).drop('x'))

        with pytest.raises(ValueError):
            xrft.cross_phase(da1, da2.rename({'x':'y'})) 
Example 27
Project: FRIDA   Author: LCAV   File: utils.py    MIT License 5 votes vote down vote up
def polar_error(x1, x2):

    tp = 2*np.pi
    e = np.minimum(np.mod(x1-x2, tp), np.mod(x2-x1, tp))

    return e 
Example 28
Project: FRIDA   Author: LCAV   File: generators.py    MIT License 5 votes vote down vote up
def gen_sig_at_mic(sigmak2_k, phi_k, pos_mic_x,
                   pos_mic_y, omega_band, sound_speed,
                   SNR, Ns=256):
    """
    generate complex base-band signal received at microphones
    :param sigmak2_k: the variance of the circulant complex Gaussian signal
                emitted by the K sources
    :param phi_k: source locations (azimuths)
    :param pos_mic_x: a vector that contains microphones' x coordinates
    :param pos_mic_y: a vector that contains microphones' y coordinates
    :param omega_band: mid-band (ANGULAR) frequency [radian/sec]
    :param sound_speed: speed of sound
    :param SNR: SNR for the received signal at microphones
    :param Ns: number of snapshots used to estimate the covariance matrix
    :return: y_mic: received (complex) signal at microphones
    """
    num_mic = pos_mic_x.size
    xk, yk = polar2cart(1, phi_k)  # source locations in cartesian coordinates
    # reshape to use broadcasting
    xk = np.reshape(xk, (1, -1), order='F')
    yk = np.reshape(yk, (1, -1), order='F')
    pos_mic_x = np.reshape(pos_mic_x, (-1, 1), order='F')
    pos_mic_y = np.reshape(pos_mic_y, (-1, 1), order='F')

    t = np.reshape(np.linspace(0, 10 * np.pi, num=Ns), (1, -1), order='F')
    K = sigmak2_k.size
    sigmak2_k = np.reshape(sigmak2_k, (-1, 1), order='F')

    # x_tilde_k size: K x length_of_t
    # circular complex Gaussian process
    x_tilde_k = np.sqrt(sigmak2_k / 2.) * (np.random.randn(K, Ns) + 1j *
                                           np.random.randn(K, Ns))
    y_mic = np.dot(np.exp(-1j * (xk * pos_mic_x + yk * pos_mic_y) / (sound_speed / omega_band)),
                   x_tilde_k * np.exp(1j * omega_band * t))
    signal_energy = linalg.norm(y_mic, 'fro') ** 2
    noise_energy = signal_energy / 10 ** (SNR * 0.1)
    sigma2_noise = noise_energy / (Ns * num_mic)
    noise = np.sqrt(sigma2_noise / 2.) * (np.random.randn(*y_mic.shape) + 1j *
                                          np.random.randn(*y_mic.shape))
    y_mic_noisy = y_mic + noise
    return y_mic_noisy, y_mic 
Example 29
Project: DataHack2018   Author: InnovizTech   File: vis_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def euler_matrix(yaw, pitch, roll):
    # angles in radians
    # negating pitch for easier computation
    pi=np.pi+1e-7
    assert -pi <= yaw <= pi and -np.pi <= roll <= np.pi and -np.pi/2 <= pitch <= np.pi/2, \
        "Erroneous yaw, pitch, roll={},{},{}".format(yaw, pitch, roll)
    rotX = np.eye(3)
    rotY = np.eye(3)
    rotZ = np.eye(3)
    rotX[1:, 1:] = rot_mat_2d(roll)
    rotY[::2, ::2] = rot_mat_2d(-pitch)
    rotZ[:2, :2] = rot_mat_2d(yaw)

    return rotX.dot(rotY.dot(rotZ)) 
Example 30
Project: DataHack2018   Author: InnovizTech   File: vis_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def _pan_camera(self, movement):
        camera_pos = base.camera.getPos()
        hpr = base.camera.getHpr()
        rot_mat = lambda angle: np.array([[np.cos(angle), -np.sin(angle)],
                                          [np.sin(angle), np.cos(angle)]])
        movement_vec_xy = np.dot(rot_mat(hpr[0] / 180. * np.pi), np.array([movement[0], movement[1]]))
        base.camera.setPos(camera_pos[0] + movement_vec_xy[0],
                           camera_pos[1] + movement_vec_xy[1],
                           camera_pos[2] + movement[2])
        self._update_location_text() 
Example 31
Project: DataHack2018   Author: InnovizTech   File: math_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def yaw(points):
    # calculate yaw of points given in nx3 format. yaw in [-pi, pi]
    assert isinstance(points, np.ndarray)
    assert points.ndim == 2
    assert points.shape[1] == 3
    return np.arctan2(points[:, 1], points[:, 0]) 
Example 32
Project: DataHack2018   Author: InnovizTech   File: math_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def pitch(points):
    # calculate pitch of points given in nx3 format. pitch in [-pi/2, pi/2]
    assert isinstance(points, np.ndarray)
    assert points.ndim == 2
    assert points.shape[1] == 3
    return np.arctan2(points[:, 2], np.linalg.norm(points[:, :2], axis=1)) 
Example 33
Project: DataHack2018   Author: InnovizTech   File: math_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def pose_from_rt(rotation, translation):
    # rotation is received as roll, pitch, yaw
    # Convert to rotation translation. Roll and yaw angles are in [-pi, pi] while pitch is [-pi/2, pi/2]
    return rotation_translation_to_4d(rotation_angles_to_rotation(rotation), translation) 
Example 34
Project: DataHack2018   Author: InnovizTech   File: math_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def axis_angle_to_rotation(axis, angle):
    # converts axis angle notation to a rotation matrix. angle is assumed in radians and in [0, 2pi],
    # axis should be normalized.
    # formula from https://en.wikipedia.org/wiki/Rotation_matrix#Conversion_from_and_to_axis-angle
    assert isinstance(axis, np.ndarray) and axis.shape == (3,)
    assert np.abs(np.linalg.norm(axis) - 1.) < 1e-6
    assert 0 <= angle <= np.pi * 2
    rotation_matrix = np.cos(angle) * np.eye(3) + np.sin(angle) * cross_product_matrix(axis) + \
                      (1 - np.cos(angle)) * np.tensordot(axis, axis, axes=0)
    return rotation_matrix 
Example 35
Project: DataHack2018   Author: InnovizTech   File: math_utils.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def extract_axis_angle(rot_mat):
    # Convert from rotation matrix to axis angle. This conversion is good for angles in [0, pi], angles in [-pi, 0] will
    # be mapped to [pi, 0] (pi downto 0) with the negative axis. To handle this issue we can use quaternions.
    assert isinstance(rot_mat, np.ndarray) and rot_mat.shape == (3, 3,)
    u = np.array([rot_mat[2, 1] - rot_mat[1, 2],
                  rot_mat[0, 2] - rot_mat[2, 0],
                  rot_mat[1, 0] - rot_mat[0, 1]])
    angle = np.arccos(np.trace(rot_mat[:3, :3]) / 2 - 0.5)
    if np.linalg.norm(u) == 0.:
        return np.array([0., 0., 0., 1.])
    else:
        u = u / np.linalg.norm(u)
        return np.array([angle, u[0], u[1], u[2]]) 
Example 36
Project: StructEngPy   Author: zhuoju36   File: frame_section.py    MIT License 5 votes vote down vote up
def __init__(self,mat,d,name=None):
        """
        d - diameter
        """
        self.d=d
        A=np.pi*d**2/4
        J=np.pi*d**4/32
        I33=np.pi*d**4/64
        I22=I33
        W33=I33/d*2
        W22=W33
        super(Circle,self).__init__(mat,A,J,I33,I22,W33,W22,name)
#        self.gamma33=1.15
#        self.gamma22=1.15 
Example 37
Project: StructEngPy   Author: zhuoju36   File: frame_section.py    MIT License 5 votes vote down vote up
def __init__(self,mat,d,t,name=None):
        """
        d - diameter\n
        t - thickness of wall\n
        fab - fabrication\n
            'r' - rolled\n
            'w' - welded\n
        """
        self._d=d
        self._t=t
        A=np.pi*d**2/4-np.pi*(d-2*t)**2/4
        J=np.pi*(d-t)/t*2*A
        I33=np.pi*d**4/64*(1-((d-2*t)/d)**4)
        I22=I33
        W33=I33/d*2
        W22=W33
        super(Pipe,self).__init__(mat,A,J,I33,I22,W33,W22,name)
        
#        self.gamma33=1.15
#        self.gamma22=1.15
#        if fab=='r':
#            self.cls33='b'
#            self.cls22='b'
#        elif fab=='w':
#            self.cls33='c'
#            self.cls22='c'
#        else:
#            raise ValueError('wrong fabrication!') 
Example 38
Project: StructEngPy   Author: zhuoju36   File: __init__.py    MIT License 5 votes vote down vote up
def period(self):
        return 2*np.pi/(self.omega_) 
Example 39
Project: Griffin_lim   Author: candlewill   File: audio.py    MIT License 5 votes vote down vote up
def _griffin_lim(S):
    angles = np.exp(2j * np.pi * np.random.rand(*S.shape))
    S_complex = np.abs(S).astype(np.complex)
    for i in range(hparams.griffin_lim_iters):
        if i > 0:
            angles = np.exp(1j * np.angle(_stft(y)))
        y = _istft(S_complex * angles)
    return y 
Example 40
Project: cplot   Author: sunchaoatmo   File: taylorDiagram.py    GNU General Public License v3.0 5 votes vote down vote up
def add_contours(self, levels=5, **kwargs):
        """Add constant centered RMS difference contours."""

        rs,ts = NP.meshgrid(NP.linspace(self.smin,self.smax),
                            NP.linspace(0,NP.pi/2.0))
        # Compute centered RMS difference
        rms = NP.sqrt(self.refstd**2 + rs**2 - 2*self.refstd*rs*NP.cos(ts))
        
        contours = self.ax.contour(ts, rs, rms, levels,linewidths=0.5,ls='--',zorder=0, **kwargs)
        #contours = self.ax.contour(ts, rs, rms, levels,linewidths=0.01, **kwargs)

        return contours 
Example 41
Project: skylab   Author: coenders   File: psLLH.py    GNU General Public License v3.0 5 votes vote down vote up
def __init__(self, exp, mc, livetime, llh_model, scramble=True, mode="box",
                 delta_ang=np.deg2rad(10.), thresh_S=0., **kwargs):
        super(PointSourceLLH, self).__init__(**kwargs)

        # Add sine declination to experimental data if not available.
        if "sinDec" not in exp.dtype.fields:
            exp = numpy.lib.recfunctions.append_fields(
                exp, names="sinDec", data=np.sin(exp["dec"]),
                dtypes=np.float, usemask=False)

        # Scramble experimental data in right ascension if data was not
        # unblinded yet.
        if scramble:
            exp["ra"] = self.random.uniform(0., 2.*np.pi, exp.size)
        else:
            logger = logging.getLogger(self._logname)
            logger.warn("Working on >> UNBLINDED << data.")

        self.exp = exp
        self.livetime = livetime
        self.set_llh_model(llh_model, mc=mc)

        # Calculate background probability here because it does not change.
        self.exp = numpy.lib.recfunctions.append_fields(
            self.exp, names="B", data=llh_model.background(exp), usemask=False)

        self.mode = mode
        self.delta_ang = delta_ang
        self.thresh_S = thresh_S

        self._events = None
        self._signal = None 
Example 42
Project: skylab   Author: coenders   File: ps_model.py    GNU General Public License v3.0 5 votes vote down vote up
def background(self, ev):
        return np.full(len(ev), 1. / 4. / np.pi) 
Example 43
Project: skylab   Author: coenders   File: utils.py    GNU General Public License v3.0 5 votes vote down vote up
def skymap(plt, vals, **kwargs):
    fig, ax = plt.subplots(subplot_kw=dict(projection="aitoff"))

    gridsize = 1000

    x = np.linspace(np.pi, -np.pi, 2 * gridsize)
    y = np.linspace(np.pi, 0., gridsize)

    X, Y = np.meshgrid(x, y)

    r = hp.rotator.Rotator(rot=(-180., 0., 0.))

    YY, XX = r(Y.ravel(), X.ravel())

    pix = hp.ang2pix(hp.npix2nside(len(vals)), YY, XX)

    Z = np.reshape(vals[pix], X.shape)

    lon = x[::-1]
    lat = np.pi /2.  - y

    cb = kwargs.pop("colorbar", dict())
    cb.setdefault("orientation", "horizontal")
    cb.setdefault("fraction", 0.075)

    title = cb.pop("title", None)

    p = ax.pcolormesh(lon, lat, Z, **kwargs)

    cbar = fig.colorbar(p, **cb)

    cbar.solids.set_edgecolor("face")
    cbar.update_ticks()
    if title is not None:
        cbar.set_label(title)

    ax.xaxis.set_ticks([])

    return fig, ax 
Example 44
Project: skylab   Author: coenders   File: data.py    GNU General Public License v3.0 5 votes vote down vote up
def init(Nexp, NMC, energy=True, **kwargs):
    Nsrc = kwargs.pop("Nsrc", 0)

    arr_exp = exp(Nexp - Nsrc)
    arr_mc = MC(NMC)

    if Nsrc > 0:
        inj = PointSourceInjector(2, sinDec_bandwidth=1, seed=0)
        inj.fill(0., arr_mc, 333.)

        source = inj.sample(np.pi, Nsrc, poisson=False).next()[1]

        arr_exp = np.append(arr_exp, source)

    if energy:
        llh_model = PowerLawLLH(["logE"], min(50, Nexp // 50),
                                range=[[0.9 * arr_mc["logE"].min(),
                                        1.1 * arr_mc["logE"].max()]],
                                sinDec_bins=min(50, Nexp // 50),
                                sinDec_range=[-1., 1.],
                                bounds=(0, 5))
    else:
        llh_model = UniformLLH(sinDec_bins=max(3, Nexp // 200),
                               sinDec_range=[-1., 1.])

    llh = PointSourceLLH(arr_exp, arr_mc, 365., llh_model=llh_model,
                         mode="all", nsource=25, scramble=False,
                         nsource_bounds=(-Nexp / 2., Nexp / 2.)
                                        if not energy else (0., Nexp / 2.),
                         seed=np.random.randint(2**32),
                         **kwargs)

    return llh 
Example 45
Project: deep-learning-note   Author: wdxtub   File: 4_simulate_sin.py    MIT License 5 votes vote down vote up
def draw_correct_line():
    x = np.arange(0, 2 * np.pi, 0.01)
    x = x.reshape((len(x), 1))
    y = np.sin(x)

    pylab.plot(x, y, label='标准 sin 曲线')
    plt.axhline(linewidth=1, color='r')


# 返回训练样本 
Example 46
Project: deep-learning-note   Author: wdxtub   File: 4_simulate_sin.py    MIT License 5 votes vote down vote up
def get_train_data():
    train_x = np.random.uniform(0.0, 2*np.pi, (1))
    train_y = np.sin(train_x)
    return train_x, train_y


# 定义前向结构 
Example 47
Project: deep-learning-note   Author: wdxtub   File: 4_simulate_sin.py    MIT License 5 votes vote down vote up
def train():
    # 学习率
    lr = 0.01

    x = tf.placeholder(tf.float32)
    y = tf.placeholder(tf.float32)

    y_ = inference(x)
    # 损失函数
    loss = tf.square(y_ - y)

    # 随机梯度下降
    opt = tf.train.GradientDescentOptimizer(lr)
    train_op = opt.minimize(loss)

    init = tf.global_variables_initializer()

    with tf.Session() as sess:
        sess.run(init)
        print("start training")
        for i in range(1000000):
            train_x, train_y = get_train_data()
            sess.run(train_op, feed_dict={x: train_x, y: train_y})

            if i % 10000 == 0:
                times = int(i / 10000)
                test_x_ndarray = np.arange(0, 2 * np.pi, 0.01)
                test_y_ndarray = np.zeros([len(test_x_ndarray)])
                ind = 0
                for test_x in test_x_ndarray:
                    test_y = sess.run(y_, feed_dict={x: test_x, y: 1})
                    np.put(test_y_ndarray, ind, test_y)
                    ind += 1
                draw_correct_line()
                pylab.plot(test_x_ndarray, test_y_ndarray, '--', label= str(times)+'times')
                pylab.show() 
Example 48
Project: NiBetaSeries   Author: HBClab   File: conftest.py    MIT License 5 votes vote down vote up
def preproc_file(deriv_dir, sub_metadata, deriv_bold_fname=deriv_bold_fname):
    deriv_bold = deriv_dir.ensure(deriv_bold_fname)
    with open(str(sub_metadata), 'r') as md:
        bold_metadata = json.load(md)
    tr = bold_metadata["RepetitionTime"]
    # time_points
    tp = 200
    ix = np.arange(tp)
    # create voxel timeseries
    task_onsets = np.zeros(tp)
    # add activations at every 40 time points
    # waffles
    task_onsets[0::40] = 1
    # fries
    task_onsets[3::40] = 1.5
    # milkshakes
    task_onsets[6::40] = 2
    signal = np.convolve(task_onsets, spm_hrf(tr))[0:len(task_onsets)]
    # csf
    csf = np.cos(2*np.pi*ix*(50/tp)) * 0.1
    # white matter
    wm = np.sin(2*np.pi*ix*(22/tp)) * 0.1
    # voxel time series (signal and noise)
    voxel_ts = signal + csf + wm
    # a 4d matrix with 2 identical timeseries
    img_data = np.array([[[voxel_ts, voxel_ts]]])
    # make a nifti image
    img = nib.Nifti1Image(img_data, np.eye(4))
    # save the nifti image
    img.to_filename(str(deriv_bold))

    return deriv_bold 
Example 49
Project: NiBetaSeries   Author: HBClab   File: conftest.py    MIT License 5 votes vote down vote up
def confounds_file(deriv_dir, preproc_file,
                   deriv_regressor_fname=deriv_regressor_fname):
    confounds_file = deriv_dir.ensure(deriv_regressor_fname)
    tp = nib.load(str(preproc_file)).shape[-1]
    ix = np.arange(tp)
    # csf
    csf = np.cos(2*np.pi*ix*(50/tp)) * 0.1
    # white matter
    wm = np.sin(2*np.pi*ix*(22/tp)) * 0.1
    confounds_df = pd.DataFrame({'WhiteMatter': wm, 'CSF': csf})
    confounds_df.to_csv(str(confounds_file), index=False, sep='\t')
    return confounds_file 
Example 50
Project: spqrel_tools   Author: LCAS   File: urdf.py    MIT License 5 votes vote down vote up
def to_openrave_xml(self, doc):
        limit = create_element(doc, 'limitsdeg', [round(self.lower*180/pi, 1), round(self.upper*180/pi, 1)])
        maxvel = create_element(doc, 'maxvel', self.velocity)
        maxtrq = create_element(doc, 'maxtorque', self.effort)
        return [limit, maxvel, maxtrq] 
Example 51
Project: dynamic-training-with-apache-mxnet-on-aws   Author: awslabs   File: loss.py    Apache License 2.0 5 votes vote down vote up
def hybrid_forward(self, F, pred, target, sample_weight=None, epsilon=1e-08):
        target = _reshape_like(F, target, pred)
        if self._from_logits:
            loss = F.exp(pred) - target * pred
        else:
            loss = pred - target * F.log(pred + epsilon)
        if self._compute_full:
            # Using numpy's pi value
            stirling_factor = target * F.log(target)- target + 0.5 * F.log(2 * target * np.pi)
            target_gt_1 = target > 1
            stirling_factor *= target_gt_1
            loss += stirling_factor
        loss = _apply_weighting(F, loss, self._weight, sample_weight)
        return F.mean(loss) 
Example 52
Project: dynamic-training-with-apache-mxnet-on-aws   Author: awslabs   File: image.py    Apache License 2.0 5 votes vote down vote up
def __call__(self, src):
        """Augmenter body.
        Using approximate linear transfomation described in:
        https://beesbuzz.biz/code/hsv_color_transforms.php
        """
        alpha = random.uniform(-self.hue, self.hue)
        u = np.cos(alpha * np.pi)
        w = np.sin(alpha * np.pi)
        bt = np.array([[1.0, 0.0, 0.0],
                       [0.0, u, -w],
                       [0.0, w, u]])
        t = np.dot(np.dot(self.ityiq, bt), self.tyiq).T
        src = nd.dot(src, nd.array(t))
        return src 
Example 53
Project: dynamic-training-with-apache-mxnet-on-aws   Author: awslabs   File: test_loss.py    Apache License 2.0 5 votes vote down vote up
def test_poisson_nllloss():
    pred = mx.nd.random.normal(shape=(3, 4))
    min_pred = mx.nd.min(pred)
    #This is necessary to ensure only positive random values are generated for prediction,
    # to avoid ivalid log calculation
    pred[:] = pred + mx.nd.abs(min_pred)
    target = mx.nd.random.normal(shape=(3, 4))
    min_target = mx.nd.min(target)
    #This is necessary to ensure only positive random values are generated for prediction,
    # to avoid ivalid log calculation
    target[:] += mx.nd.abs(min_target)

    Loss = gluon.loss.PoissonNLLLoss(from_logits=True)
    Loss_no_logits = gluon.loss.PoissonNLLLoss(from_logits=False)
    #Calculating by brute formula for default value of from_logits = True

    # 1) Testing for flag logits = True
    brute_loss = np.mean(np.exp(pred.asnumpy()) - target.asnumpy() * pred.asnumpy())
    loss_withlogits = Loss(pred, target)
    assert_almost_equal(brute_loss, loss_withlogits.asscalar())

    #2) Testing for flag logits = False
    loss_no_logits = Loss_no_logits(pred, target)
    np_loss_no_logits = np.mean(pred.asnumpy() - target.asnumpy() * np.log(pred.asnumpy() + 1e-08))
    if np.isnan(loss_no_logits.asscalar()):
        assert_almost_equal(np.isnan(np_loss_no_logits), np.isnan(loss_no_logits.asscalar()))
    else:
        assert_almost_equal(np_loss_no_logits, loss_no_logits.asscalar())

    #3) Testing for Sterling approximation
    np_pred = np.random.uniform(1, 5, (2, 3))
    np_target = np.random.uniform(1, 5, (2, 3))
    np_compute_full = np.mean((np_pred - np_target * np.log(np_pred + 1e-08)) + ((np_target * np.log(np_target)-\
     np_target + 0.5 * np.log(2 * np_target * np.pi))*(np_target > 1)))
    Loss_compute_full = gluon.loss.PoissonNLLLoss(from_logits=False, compute_full=True)
    loss_compute_full = Loss_compute_full(mx.nd.array(np_pred), mx.nd.array(np_target))
    assert_almost_equal(np_compute_full, loss_compute_full.asscalar()) 
Example 54
Project: helloworld   Author: pip-uninstaller-python   File: matplotlibTest.py    GNU General Public License v2.0 5 votes vote down vote up
def main():
    # line
    x = np.linspace(-np.pi, np.pi, 256, endpoint=True)
    c, s = np.cos(x), np.sin(x)
    plt.figure(1)
    plt.plot(x, c, color="blue", linewidth=1.0, linestyle="-", label="COS", alpha=0.5)  # 自变量, 因变量
    plt.plot(x, s, "r.", label="SIN")  # 正弦  "-"/"r-"/"r."
    plt.title("COS & SIN")
    ax = plt.gca()
    ax.spines["right"].set_color("none")
    ax.spines["top"].set_color("none")
    ax.spines["left"].set_position(("data", 0))  # 横轴位置
    ax.spines["bottom"].set_position(("data", 0))  # 纵轴位置
    ax.xaxis.set_ticks_position("bottom")
    ax.yaxis.set_ticks_position("left")
    plt.xticks([-np.pi, -np.pi / 2.0, np.pi / 2, np.pi],
               [r'$-\pi/2$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$-\pi$'])
    plt.yticks(np.linspace(-1, 1, 5, endpoint=True))
    for label in ax.get_xticklabels() + ax.get_yticklabels():
        label.set_fontsize(16)
        label.set_bbox(dict(facecolor="white", edgecolor="None", alpha=0.2))
    plt.legend(loc="upper left")  # 左上角的显示图标
    plt.grid()  # 网格线
    # plt.axis([-1, 1, -0.5, 1])  # 显示范围
    plt.fill_between(x, np.abs(x) < 0.5, c, c < 0.5, color="green", alpha=0.25)
    t = 1
    plt.plot([t, t], [0, np.cos(t)], "y", linewidth=3, linestyle="--")
    # 注释
    plt.annotate("cos(1)", xy=(t, np.cos(1)), xycoords="data", xytext=(+10, +30),
                 textcoords="offset points", arrowprops=dict(arrowstyle="->", connectionstyle="arc3, rad=.2"))
    plt.show()


# Scatter --> 散点图 
Example 55
Project: FCOS_GluonCV   Author: DetectionTeamUCAS   File: pose.py    Apache License 2.0 5 votes vote down vote up
def get_affine_transform(center,
                         scale,
                         rot,
                         output_size,
                         shift=np.array([0, 0], dtype=np.float32),
                         inv=0):
    cv2 = try_import_cv2()
    if not isinstance(scale, np.ndarray) and not isinstance(scale, list):
        scale = np.array([scale, scale])

    scale_tmp = scale * 200.0
    src_w = scale_tmp[0]
    dst_w = output_size[0]
    dst_h = output_size[1]

    rot_rad = np.pi * rot / 180
    src_dir = get_dir([0, src_w * -0.5], rot_rad)
    dst_dir = np.array([0, dst_w * -0.5], np.float32)

    src = np.zeros((3, 2), dtype=np.float32)
    dst = np.zeros((3, 2), dtype=np.float32)
    src[0, :] = center + scale_tmp * shift
    src[1, :] = center + src_dir + scale_tmp * shift
    dst[0, :] = [dst_w * 0.5, dst_h * 0.5]
    dst[1, :] = np.array([dst_w * 0.5, dst_h * 0.5]) + dst_dir

    src[2:, :] = get_3rd_point(src[0, :], src[1, :])
    dst[2:, :] = get_3rd_point(dst[0, :], dst[1, :])

    if inv:
        trans = cv2.getAffineTransform(np.float32(dst), np.float32(src))
    else:
        trans = cv2.getAffineTransform(np.float32(src), np.float32(dst))

    return trans 
Example 56
Project: estim2bapi   Author: fredhatt   File: motion.py    MIT License 5 votes vote down vote up
def calc_angles(self, pos=None):
        if pos is None:
            xyz = np.array(self.hist)[:, 1:4]
        else:
            xyz = np.array(self.hist)[pos, 1:4]
            xyz = np.reshape(xyz, (-1, 3)) # single entry
        pitch = np.arctan(xyz[:, 0] / np.sqrt(xyz[:, 1]**2 + xyz[:, 2]**2)) * 180. / np.pi
        roll = np.arctan(xyz[:, 1] / np.sqrt(xyz[:, 0]**2 + xyz[:, 2]**2)) * 180. / np.pi
        return pitch, roll 
Example 57
Project: estim2bapi   Author: fredhatt   File: motion.py    MIT License 5 votes vote down vote up
def calc_angles(self, pos=None):
        if pos is None:
            xyz = np.array(self.hist)[:, 1:4]
        else:
            xyz = np.array(self.hist)[pos, 1:4]
            xyz = np.reshape(xyz, (-1, 3)) # single entry
        pitch = np.arctan(xyz[:, 0] / np.sqrt(xyz[:, 1]**2 + xyz[:, 2]**2)) * 180. / np.pi
        roll = np.arctan(xyz[:, 1] / np.sqrt(xyz[:, 0]**2 + xyz[:, 2]**2)) * 180. / np.pi
        return pitch, roll 
Example 58
Project: UBQC   Author: BlindQlouder   File: simulation.py    GNU General Public License v3.0 5 votes vote down vote up
def client_protocol(client):

    def compute_delta(angle, theta, m, r1, r2, r3, s1, s2): # this is the angle the server will have to measure at
        return ((-1)**(1+s2+r2)*angle + theta + (m+s1+r1+r3)*pi)%(2*pi)

    angle = np.random.randint(0, 3, 4)*pi/4     # these are the actual angles only known to the client; the set {0, pi/4, pi/2} is sufficient for universal computation
    angle[3] = 0    # this will just do a Hadamard at the end
    theta = np.random.randint(0, 8, 5)*pi/8     # the hiding angles
    r = np.random.randint(0, 2, 4)              # for hiding the measurement outcome
    s1 = 0
    s2 = 0
    r1 = 0
    r2 = 0
    r3 = r[0]
    m = client.RS(theta[0])
    delta = compute_delta(angle[0], theta[0], m, r1, r2, r3, s1, s2)

    for i in range(1,5):
        client.send_angle(delta)
        m = client.RS(theta[i])
        s1 = s2
        s2 = client.get_measurement()
        if i < 4:
            r1 = r2
            r2 = r3
            r3 = r[i]
            delta = compute_delta(angle[i], theta[i], m, r1, r2, r3, s1, s2)



    clean_state = qt.rz(angle[2]) * qt.rx(angle[1]) * qt.rz(angle[0]) * 1./np.sqrt(2)*qt.Qobj([[1],[1]])
    expected_state = qt.rz(theta[4]) * qt.rx(pi*(s2+r[3])) * qt.rz(pi*(m+s1+r[2])) *  qt.rz(angle[2]) * qt.rx(angle[1]) * qt.rz(angle[0]) * 1./np.sqrt(2)*qt.Qobj([[1],[1]])
    server_state = client.conn.recv()

    print("--------------------")
    print("C: final state info: \nprep_state correction m = {}, hidden angle : theta_5 =  {:.4f}, dependancies : s4 = {} s3 = {},  hidden signaling : r4 = {} r3 = {}".format(m,theta[4],s2,s1,r[3],r[2]))
    #print("|psi_out> = R_Z(theta_5) X**(s_4+r_4) Z**(m+s_3+r_3) R_Z(alpha) R_X(beta) R_Z(gamma) |psi_in>")
    print("non-hidden state:\n", clean_state.full())
    print("state on the server side:\n", server_state.full())
    print("state should be:\n", expected_state.full())
    print("overlap: ", abs(expected_state.overlap(server_state))) 
Example 59
Project: UBQC   Author: BlindQlouder   File: cliserv.py    GNU General Public License v3.0 5 votes vote down vote up
def CPHASE(self):
        print("S: applying CPHASE")
        self.system = qt.cphase(pi)*self.system
        return 1 
Example 60
Project: UBQC   Author: BlindQlouder   File: cliserv.py    GNU General Public License v3.0 5 votes vote down vote up
def RS(self, theta):    # projection onto a basis { |0> + exp(i*theta)|1>, |0> - exp(i*theta)|1> }
        print("C: applying RS")
        m = random.randint(0,1) # measurement outcome
        self.conn.send(theta+m*pi)
        return m 
Example 61
Project: UBQC   Author: BlindQlouder   File: DNMEnc.py    GNU General Public License v3.0 5 votes vote down vote up
def plot_statistics(states):
    print('plotting histogram...')
    data = []
    for state in states:
        #rho = state * state.dag()
        w = np.abs(state[0])**2 - np.abs(state[1])**2
        u = 2. * (state[0] * state[1].conj()).real
        v = 2. * (state[1] * state[0].conj()).imag
        theta = np.arccos(w)
        if v>=0 :
            phi = np.arccos(u/np.sin(theta))
        else:
            phi = 2*np.pi - np.arccos(u/np.sin(theta))


        data.append([theta, phi])

    data = np.array(data)
    #data.reshape(data.shape[0],data.shape[1])
    h1, b1 = np.histogram(data[:,0], bins=200, range=(0,np.pi))
    h2, b2 = np.histogram(data[:,1], bins=200, range=(0,2*np.pi))
    plt.figure()
    plt.plot(b1[1:]*2, h1, label='2*theta')
    plt.plot(b2[1:], h2, label='phi')
    plt.legend()
    plt.show() 
Example 62
Project: DOTA_models   Author: ringringyi   File: real_nvp_utils.py    Apache License 2.0 5 votes vote down vote up
def standard_normal_ll(input_):
    """Log-likelihood of standard Gaussian distribution."""
    res = -.5 * (tf.square(input_) + numpy.log(2. * numpy.pi))

    return res 
Example 63
Project: DOTA_models   Author: ringringyi   File: map_utils.py    Apache License 2.0 5 votes vote down vote up
def get_map_to_predict(src_locs, src_x_axiss, src_y_axiss, map, map_size,
                       interpolation=cv2.INTER_LINEAR):
  fss = []
  valids = []

  center = (map_size-1.0)/2.0
  dst_theta = np.pi/2.0
  dst_loc = np.array([center, center])
  dst_x_axis = np.array([np.cos(dst_theta), np.sin(dst_theta)])
  dst_y_axis = np.array([np.cos(dst_theta+np.pi/2), np.sin(dst_theta+np.pi/2)])

  def compute_points(center, x_axis, y_axis):
    points = np.zeros((3,2),dtype=np.float32)
    points[0,:] = center
    points[1,:] = center + x_axis
    points[2,:] = center + y_axis
    return points

  dst_points = compute_points(dst_loc, dst_x_axis, dst_y_axis)
  for i in range(src_locs.shape[0]):
    src_loc = src_locs[i,:]
    src_x_axis = src_x_axiss[i,:]
    src_y_axis = src_y_axiss[i,:]
    src_points = compute_points(src_loc, src_x_axis, src_y_axis)
    M = cv2.getAffineTransform(src_points, dst_points)

    fs = cv2.warpAffine(map, M, (map_size, map_size), None, flags=interpolation,
                        borderValue=np.NaN)
    valid = np.invert(np.isnan(fs))
    valids.append(valid)
    fss.append(fs)
  return fss, valids 
Example 64
Project: DOTA_models   Author: ringringyi   File: nav_env.py    Apache License 2.0 5 votes vote down vote up
def _get_relative_goal_loc(goal_loc, loc, theta):
  r = np.sqrt(np.sum(np.square(goal_loc - loc), axis=1))
  t = np.arctan2(goal_loc[:,1] - loc[:,1], goal_loc[:,0] - loc[:,0])
  t = t-theta[:,0] + np.pi/2
  return np.expand_dims(r,axis=1), np.expand_dims(t, axis=1) 
Example 65
Project: DOTA_models   Author: ringringyi   File: nav_env.py    Apache License 2.0 5 votes vote down vote up
def pre(self, inputs):
    inputs['imgs'] = image_pre(inputs['imgs'], self.task_params.modalities)
    if inputs['loc_on_map'] is not None:
      inputs['loc_on_map'] = inputs['loc_on_map'] - inputs['loc_on_map'][:,[0],:]
    if inputs['theta_on_map'] is not None:
      inputs['theta_on_map'] = np.pi/2. - inputs['theta_on_map']
    return inputs 
Example 66
Project: DOTA_models   Author: ringringyi   File: policy.py    Apache License 2.0 5 votes vote down vote up
def entropy(self, logits,
              sampling_dim, act_dim, act_type):
    """Calculate entropy of distribution."""
    if self.env_spec.is_discrete(act_type):
      entropy = tf.reduce_sum(
          -tf.nn.softmax(logits) * tf.nn.log_softmax(logits), -1)
    elif self.env_spec.is_box(act_type):
      means = logits[:, :sampling_dim / 2]
      std = logits[:, sampling_dim / 2:]
      entropy = tf.reduce_sum(
          0.5 * (1 + tf.log(2 * np.pi * tf.square(std))), -1)
    else:
      assert False

    return entropy 
Example 67
Project: DOTA_models   Author: ringringyi   File: dsn_eval.py    Apache License 2.0 5 votes vote down vote up
def angle_diff(true_q, pred_q):
  angles = 2 * (
      180.0 /
      np.pi) * np.arccos(np.abs(np.sum(np.multiply(pred_q, true_q), axis=1)))
  return angles 
Example 68
Project: soccer-matlab   Author: utra-robosoccer   File: kuka_diverse_object_gym_env.py    BSD 2-Clause "Simplified" License 5 votes vote down vote up
def _randomly_place_objects(self, urdfList):
    """Randomly places the objects in the bin.

    Args:
      urdfList: The list of urdf files to place in the bin.

    Returns:
      The list of object unique ID's.
    """


    # Randomize positions of each object urdf.
    objectUids = []
    for urdf_name in urdfList:
      xpos = 0.4 +self._blockRandom*random.random()
      ypos = self._blockRandom*(random.random()-.5)
      angle = np.pi/2 + self._blockRandom * np.pi * random.random()
      orn = p.getQuaternionFromEuler([0, 0, angle])
      urdf_path = os.path.join(self._urdfRoot, urdf_name)
      uid = p.loadURDF(urdf_path, [xpos, ypos, .15],
        [orn[0], orn[1], orn[2], orn[3]])
      objectUids.append(uid)
      # Let each object fall to the tray individual, to prevent object
      # intersection.
      for _ in range(500):
        p.stepSimulation()
    return objectUids 
Example 69
Project: soccer-matlab   Author: utra-robosoccer   File: robot_locomotors.py    BSD 2-Clause "Simplified" License 5 votes vote down vote up
def robot_specific_reset(self, bullet_client):
		WalkerBase.robot_specific_reset(self, bullet_client)
		self.motor_names  = ["abdomen_z", "abdomen_y", "abdomen_x"]
		self.motor_power  = [100, 100, 100]
		self.motor_names += ["right_hip_x", "right_hip_z", "right_hip_y", "right_knee"]
		self.motor_power += [100, 100, 300, 200]
		self.motor_names += ["left_hip_x", "left_hip_z", "left_hip_y", "left_knee"]
		self.motor_power += [100, 100, 300, 200]
		self.motor_names += ["right_shoulder1", "right_shoulder2", "right_elbow"]
		self.motor_power += [75, 75, 75]
		self.motor_names += ["left_shoulder1", "left_shoulder2", "left_elbow"]
		self.motor_power += [75, 75, 75]
		self.motors = [self.jdict[n] for n in self.motor_names]
		if self.random_yaw:
			position = [0,0,0]
			orientation = [0,0,0]
			yaw = self.np_random.uniform(low=-3.14, high=3.14)
			if self.random_lean and self.np_random.randint(2)==0:
				cpose.set_xyz(0, 0, 1.4)
				if self.np_random.randint(2)==0:
					pitch = np.pi/2
					position = [0, 0, 0.45]
				else:
					pitch = np.pi*3/2
					position = [0, 0, 0.25]
				roll = 0
				orientation = [roll, pitch, yaw]
			else:
				position = [0, 0, 1.4]
				orientation = [0, 0, yaw]  # just face random direction, but stay straight otherwise
			self.robot_body.reset_position(position)
			self.robot_body.reset_orientation(orientation)
		self.initial_z = 0.8 
Example 70
Project: oceanmapper   Author: vtamsitt   File: oceanmapper.py    Apache License 2.0 4 votes vote down vote up
def topography3d(mode,topo_x=None,topo_y=None,topo_z=None,topo_limits=None,zscale=500.,topo_vmin=None,topo_vmax=None,topo_cmap='bone',topo_cmap_reverse=False,land_constant=False,land_color=(0.7,0.7,0.7),set_view=None):
    """
    mode: string; coordinate system of 3D projection. Options are 'rectangle' (default), 'sphere' or 'cylinder'
    topo: array_like, optional; input topography file, default is etopo 30 
##TODO: need to define sign of topography 
    topo_limits: array_like, optional; longitude and latitude limits for 3d topography plot [lon min, lon max, lat min, lat max], longitudes range -180 to 180, latitude -90 to 90, default is entire globe
    zscale: scalar, optional; change vertical scaling for plotting, such that the vertical axis is scaled as topo_z/zscale (assumes topo_z units are m); default zscale is 500
    topo_cmap: string, optional; set colormap for topography, default is bone 
    topo_cmap_reverse: string, optional; reverse topography colormap, default is false
    land_constant: string optional; if True, land is set to one colour, default is False
    land_color: color, optional; RGB triplet specifying land colour, defauly is gret
    set_view: array_like, optional; set the mayavi camera angle with input [azimuth, elevation, distance, focal point], default is None 
    """
        
    #load topo data
    if topo_x is not None and topo_y is not None and topo_z is not None:
        xraw = topo_x
        yraw = topo_y
        zraw = topo_z

    else:
        tfile = np.load('etopo1_30min.npz')
        xraw = tfile['x']
        yraw = tfile['y']
        zraw = np.swapaxes(tfile['z'][:,:],0,1)
    

    #create coordinate variables
    phi = (yraw[:]*np.pi*2)/360.+np.pi/2.
    theta = (xraw[:]*np.pi*2)/360.
    c = zraw
    theta=np.append(theta,theta[0])
    c = np.concatenate((c,np.expand_dims(c[0,:],axis=0)),axis=0)

    if topo_vmin is None:
        tvmin = 0
    else:
        tvmin = topo_vmin
    if topo_vmax is None:
        tvmax = 7000
    else:
        tvmax = topo_vmax
   
    if topo_limits is not None:
        phi_1 = topo_limits[2]
        phi_2 = topo_limits[3]
        theta_1 = topo_limits[0]
        theta_2 = topo_limits[1]

        phi_ind1 = np.argmin(np.abs(yraw-phi_1))
        phi_ind2 = np.argmin(np.abs(yraw-phi_2))
        theta_ind1 = np.argmin(np.abs(xraw-theta_1))
        theta_ind2 = np.argmin(np.abs(xraw-theta_2))

        #restrict topo extent
        phi=phi[phi_ind1:phi_ind2]
        theta=theta[theta_ind1:theta_ind2]
        c = c[theta_ind1:theta_ind2:,phi_ind1:phi_ind2] 
Example 71
Project: oceanmapper   Author: vtamsitt   File: oceanmapper.py    Apache License 2.0 4 votes vote down vote up
def surface3d(mode,xdata,ydata,zdata,fig=None,scalardata=None,vmin=None,vmax=None,data_cmap='blue-red',data_color=(0.5,0.5,0.5),data_alpha=1,zscale=500.,set_view=None):
    """
    fig: integer or string, optional. Figure key will plot data on corresponding mlab figure, if it exists, or create a new one
    mode: string; coordinate system of 3D projection. Options are 'rectangle' (default), 'sphere' or 'cylinder'
    xdata: 1D numpy array; longitude values for data array
    ydata: 1D numpy array; latitude values for data array
    zdata: 1D numpy array; depth values for data array
    scalardata: 2D numpy array, optional; 2D scalar field to plot colors on surface
    vmin: float, optional; colorbar minimum for data
    vmax: float, optional; colorbar maximum for data
    zscale: scalar, optional; change vertical scaling for plotting, such that the vertical axis is scaled as topo_z/zscale (assumes topo_z units are m); default zscale is 500 
    data_cmap: string, optional; colormap for data surface, default is blue-red
    data_color: triplet of floats ranging from 0 to 1, optional; sets color of surface, overrides colormap when scalardata is not included
    data_alpha: float or int, optional; opacity for data surface from 0 to 1, default is 1
    set_view: array_like, optional; set the mayavi camera angle with input [azimuth, elevation, distance, focal point], default is 
    """
        
    
    #make figure
    if fig is None: 
        mlab.figure(size = (1024,768),bgcolor = (1,1,1), fgcolor = (0.5, 0.5, 0.5))
        mlab.clf()
    else:
        mlab.figure(figure=fig)
    
    #do coordinate transformation
    if xdata is not None and ydata is not None and zdata is not None:
        #TODO add an error message if not all data fields are provided
        #prep data grid
        phi_iso, theta_iso = np.meshgrid(((ydata*np.pi*2)/360.)+np.pi/2.,(xdata*np.pi*2)/360.)

        if mode is 'sphere':
            x_iso = np.sin(phi_iso) * np.cos(theta_iso[::-1]) * (1 -zdata/zscale)
            y_iso = np.sin(phi_iso) * np.sin(theta_iso[::-1]) * (1 -zdata/zscale)
            z_iso = np.cos(phi_iso) * (1 -zdata/zscale)
        elif mode is 'cylinder':
            x_iso = np.sin(phi_iso) * np.cos(theta_iso[::-1])
            y_iso = np.sin(phi_iso) * np.sin(theta_iso[::-1])
            z_iso = zdata/zscale

        elif mode is 'rectangle':
            y_iso,z_iso = np.meshgrid(ydata,zdata)
            x_iso,z_iso = np.meshgrid(xdata,zdata)
            z_iso =-z_iso/zscale
    else:
        #raise error if all three fields are not provided
        print('ERROR: not all data fields are provided. Must provide 1D data x, y and z data points')  
    
    #map data surface
    if scalardata is not None:
        m = mlab.mesh(x_iso, y_iso, z_iso,scalars=scalardata,colormap=data_cmap,vmin =vmin,vmax=vmax,opacity=data_alpha)
        m.module_manager.scalar_lut_manager.lut.nan_color = [0,0,0,0]

    else:
        m = mlab.mesh(x_iso, y_iso, z_iso,color=data_color,vmin =vmin,vmax=vmax,opacity=data_alpha)
        m.module_manager.scalar_lut_manager.lut.nan_color = [0,0,0,0]
    
    return m 
Example 72
Project: oceanmapper   Author: vtamsitt   File: oceanmapper.py    Apache License 2.0 4 votes vote down vote up
def vector3d(mode,xdata,ydata,zdata,udata,vdata,wdata,scalardata=None,fig=None,zscale=500.,vector_color=(0,0,0),vector_cmap=None,alpha=1.0,vector_mode='2darrow', scale=1, spacing=8., set_view=None):
    """
    fig: integer or string, optional. Figure key will plot data on corresponding mlab figure, if it exists, or create a new one
    mode: string; coordinate system of 3D projection. Options are 'rectangle' (default), 'sphere' or 'cylinder'
    xdata: 1D array; longitude values for data array
    ydata: 1D array; latitude values for data array
    zdata: 1D array; depth values for data array
    udata: 2D or 3D array; u vector component
    vdata: 2D or 3D array; v vector component
    wdata: 2D or 3D array; w vector component
    zscale: scalar, optional; change vertical scaling for plotting, such that the vertical axis is scaled as topo_z/zscale (assumes topo_z units are m); default zscale is 500 
    vector_mode: string, optional; style of vector plot
    color: colormap or rgb triplet,optional; color of quiver plot default is black (0,0,0). 
    alpha: float or int, optional; opacity for data surface from 0 to 1, default is 1
    scale: float or int, optional; scaling for length of vectors, default is 1. 
    spacing: int, optional; If supplied, only one out of 'spacing' data points is displayed. This option is useful to reduce the number of points displayed on large datasets Must be an integer (int or long) or None
    set_view: array_like, optional; set the mayavi camera angle with input [azimuth, elevation, distance, focal point], default is 
    """
        
    
    #make figure
    if fig is None: 
        mlab.figure(size = (1024,768),bgcolor = (1,1,1))
        mlab.clf()
    else:
        mlab.figure(figure=fig,bgcolor = (1,1,1))
    
    #do coordinate transformation
    if xdata is not None and ydata is not None and zdata is not None:
        #TODO add an error message if not all data fields are provided
        #prep data grid
        phi_iso, theta_iso = np.meshgrid(((ydata*np.pi*2)/360.)+np.pi/2.,(xdata*np.pi*2)/360.)

        if mode is 'sphere':
            x_iso = np.sin(phi_iso) * np.cos(theta_iso[::-1]) * (1 -zdata/zscale)
            y_iso = np.sin(phi_iso) * np.sin(theta_iso[::-1]) * (1 -zdata/zscale)
            z_iso = np.cos(phi_iso) * (1 -zdata/zscale)
        elif mode is 'cylinder':
            x_iso = np.sin(phi_iso) * np.cos(theta_iso[::-1])
            y_iso = np.sin(phi_iso) * np.sin(theta_iso[::-1])
            z_iso = zdata/zscale

        elif mode is 'rectangle':
            y_iso,z_iso = np.meshgrid(ydata,zdata)
            x_iso,z_iso = np.meshgrid(xdata,zdata)
            z_iso =-z_iso/zscale
    else:
        #raise error if all three fields are not provided
        print('ERROR: not all data fields are provided. Must provide 1D data x, y and z data points')  
    
    #do quiver plot 
    if scalardata is not None:
        m = mlab.quiver3d(x_iso, y_iso, z_iso, udata, vdata, wdata, scalars=scalardata, scale_mode=None,colormap=vector_cmap,mode=vector_mode,opacity=alpha,scale_factor=scale,mask_points=spacing)   
    elif vector_cmap is not None:
        m = mlab.quiver3d(x_iso, y_iso, z_iso, udata, vdata, wdata, colormap=vector_cmap,mode=vector_mode,opacity=alpha,scale_factor=scale,mask_points=spacing)   
    else:
        m = mlab.quiver3d(x_iso, y_iso, z_iso, udata, vdata, wdata, color=vector_color,mode=vector_mode,opacity=alpha,scale_factor=scale,mask_points=spacing)   
     
  
    #optional: change mayavi camera settings
    return m 
Example 73
Project: oceanmapper   Author: vtamsitt   File: oceanmapper.py    Apache License 2.0 4 votes vote down vote up
def trajectory3d(mode,xdata,ydata,zdata,fig=None,scalardata=None,vmin=None,vmax=None,color=(0,0,0),data_cmap=None,data_alpha=1,zscale=500.,tube_radius=0.01,tube_sides=15,set_view=None):
    """
    fig: integer or string, optional. Figure key will plot data on corresponding mlab figure, if it exists, or create a new one
    mode: string; coordinate system of 3D projection. Options are 'rectangle' (default), 'sphere' or 'cylinder'
    xdata: 1D array; longitude values for data array
    ydata: 1D array; latitude values for data array
    zdata: 1D array; depth values for data array
    scalardata: 1D array, optional; 1D scalar field to plot colors along trajectoy
    zscale: scalar, optional; change vertical scaling for plotting, such that the vertical axis is scaled as topo_z/zscale (assumes topo_z units are m); default zscale is 500 
    vmin: float, optional; colorbar minimum for data
    vmax: float, optional; colorbar maximum for data
    data_cmap: string, optional; colormap for data surface, default is blue-red
    data_alpha: float or int, optional; opacity for data surface from 0 to 1, default is 1
    tube_radius: float, optional; radius of tube
    tube_sides: int, optional; number of sides of tube
    set_view: array_like, optional; set the mayavi camera angle with input [azimuth, elevation, distance, focal point], default is 
    """
        
    
    #make figure
    if fig is None: 
        mlab.figure(size = (1024,768),bgcolor = (1,1,1), fgcolor = (0.5, 0.5, 0.5))
        mlab.clf()
    else:
        mlab.figure(figure=fig,bgcolor = (1,1,1))
    
    #do coordinate transformation

    if xdata is not None and ydata is not None and zdata is not None:
        phi_iso, theta_iso = np.meshgrid(((ydata*np.pi*2)/360.)+np.pi/2.,(xdata*np.pi*2)/360.)

        # Create variable dimensions
        if mode == 'sphere':
            x_iso = np.sin(phi_iso) * np.cos(theta_iso[::-1]) * (1 -zdata/zscale)
            y_iso = np.sin(phi_iso) * np.sin(theta_iso[::-1]) * (1 -zdata/zscale)
            z_iso = np.cos(phi_iso) * (1 -zdata/zscale)
        elif mode == 'cylinder':
            x_iso = np.sin(phi_iso) * np.cos(theta_iso[::-1])
            y_iso = np.sin(phi_iso) * np.sin(theta_iso[::-1])
            z_iso = -zdata/zscale
        elif mode == 'rectangle':
            x_iso = xdata
            y_iso = ydata
            z_iso =-zdata/zscale
    else:
        #raise error if all three fields are not provided
        print('ERROR: not all data fields are provided. Must provide 1D data x, y and z data points')  
    

    #map data surface
    if scalardata is not None:
        m = mlab.plot3d(x_iso,y_iso,z_iso, scalardata,opacity=data_alpha,tube_radius=tube_radius,tube_sides=tube_sides,color=color,vmin=vmin,vmax=vmax)
    
    else:
        m = mlab.plot3d(x_iso,y_iso,z_iso, opacity=data_alpha,tube_radius=tube_radius,tube_sides=tube_sides,color=color,vmin=vmin,vmax=vmax)
    
    return m 
Example 74
Project: cs207-FinalProject   Author: PYNE-AD   File: elemFunctions.py    MIT License 4 votes vote down vote up
def tan(x):
	''' Compute the tangent of an AutoDiff object and its derivative.

	INPUTS
	======
	x: an AutoDiff object

	'''
	try:
		# Value and derivative undefined when divisible by pi/2 but not pi
		# To make sure the asymptotes are undefined:
		if x.val%(np.pi/2)==0 and x.val%np.pi!=0:
			new_val = np.nan
			new_der = np.nan
			new_jacobian = np.nan
			warnings.warn('Undefined at value', RuntimeWarning)
		else:
			new_val = np.tan(x.val)
			new_der = x.der / (np.cos(x.val)**2.0)
			new_jacobian = x.jacobian / (np.cos(x.val)**2.0)
		return AutoDiff(new_val, new_der, x.n, 0, new_jacobian)
	except AttributeError:
		try:
			if x.Real%(np.pi/2)==0 and x.Real%np.pi!=0:
				ans = Dual(np.nan,np.nan)
				warnings.warn('Undefined at value', RuntimeWarning)
				return ans
			else:
				return Dual(np.tan(x.Real), x.Dual / (np.cos(x.Real))**2)
		except AttributeError:
			try:
				if x.Real%(np.pi/2)==0 and x.Real%np.pi!=0:
					ans = Dual(np.nan,np.nan)
					warnings.warn('Undefined at value', RuntimeWarning)
					return ans
				else:
					# return Dual(tan(x.Real), x.Dual / (cos(x.Real))**2)
					return sin(x)/cos(x)
			except AttributeError:
				if x%(np.pi/2)==0 and x%np.pi!=0:
					warnings.warn('Undefined at value', RuntimeWarning)
					return np.nan
				else:
					return np.tan(x)

#-------------------INVERSE TRIG FUNCTIONS-------------------#
# arc sin 
Example 75
Project: FRIDA   Author: LCAV   File: doa.py    MIT License 4 votes vote down vote up
def __init__(self, L, fs, nfft, c=343.0, num_src=1, mode='far', r=None, 
        theta=None, phi=None):

        self.L = L              # locations of mics
        self.fs = fs            # sampling frequency
        self.c = c              # speed of sound
        self.M = L.shape[1]     # number of microphones
        self.D = L.shape[0]     # number of dimensions (x,y,z)
        self.num_snap = None    # number of snapshots

        self.nfft = nfft
        self.max_bin = int(self.nfft/2) + 1
        self.freq_bins = None
        self.freq_hz = None
        self.num_freq = None

        self.num_src = self._check_num_src(num_src)
        self.sources = np.zeros([self.D, self.num_src])
        self.src_idx = np.zeros(self.num_src, dtype=np.int)
        self.phi_recon = None

        self.mode = mode
        if self.mode is 'far':
            self.r = np.ones(1)
        elif r is None:
            self.r = np.ones(1)
            self.mode = 'far'
        else:
            self.r = r
            if r == np.ones(1):
                mode = 'far'
        if theta is None:
            self.theta = np.linspace(-180., 180., 30) * np.pi / 180
        else:
            self.theta = theta
        if phi is None:
            self.phi = np.pi / 2 * np.ones(1)
        else:
            self.phi = phi

        # spatial spectrum / dirty image (FRI)
        self.P = None

        # build lookup table to candidate locations from r, theta, phi 
        from fri import FRI
        if not isinstance(self, FRI):
            self.loc = None
            self.num_loc = None
            self.build_lookup()
            self.mode_vec = None
            self.compute_mode()
        else:   # no grid search for FRI
            self.num_loc = len(self.theta) 
Example 76
Project: DataHack2018   Author: InnovizTech   File: vis_utils.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def __init__(self, max_r):
        self.max_r = max_r
        format = GeomVertexFormat.getV3c4()
        vdata = GeomVertexData('point', format, Geom.UHDynamic)
        self._pos_writer = GeomVertexWriter(vdata, 'vertex')
        self._color_writer = GeomVertexWriter(vdata, 'color')

        line_num = 60
        vdata.setNumRows(line_num)

        angles = np.linspace(0, np.pi * 2 - np.pi * 2 / line_num , line_num)

        other_rgba = (0., 0., 0.3, 0.1)
        other2_rgba = (0.1, 0.1, 0.4, 0.4)
        axis_rgba = (0.2, 0.2, 0.5, 1.0)
        max_r = 250
        for indx, angle in enumerate(angles):
            if indx % 5 == 0:
                rgba = axis_rgba
            else:
                rgba = other_rgba
            self._pos_writer.addData3d(0, 0, 0.)
            self._color_writer.addData4f(rgba[0], rgba[1], rgba[2], rgba[3])
            self._pos_writer.addData3d(max_r * np.sin(angle), max_r * np.cos(angle), 0.)
            self._color_writer.addData4f(rgba[0], rgba[1], rgba[2], rgba[3])

        grnd_prmtv = GeomLines(Geom.UHStatic)
        grnd_prmtv.addConsecutiveVertices(0, 2 * line_num)
        grnd_prmtv.closePrimitive()
        ground_geom = Geom(vdata)
        ground_geom.addPrimitive(grnd_prmtv)
        snode = GeomNode('ground_lines')
        snode.addGeom(ground_geom)

        self.points_node = base.render.attachNewNode(snode)
        self.points_node.setTwoSided(True)

        for rad in range(int(max_r)):
            color = axis_rgba
            pp = makeArc(angleDegrees=360, numSteps=160, scale=rad, color=color)
            tn = TextNode('dd')
            tn.setText(str(rad))
            tn.setTextScale(0.2)
            tn.setTextColor(color)
            text_geom = GeomNode('text')
            text_geom.addChild(tn)
            tp = NodePath(text_geom)
            tp.setPos((0, rad-0.2, 0))
            tp.setHpr((0, -90, 0))
            tp.reparentTo(self.points_node)
            pp.reparentTo(self.points_node) 
Example 77
Project: skylab   Author: coenders   File: utils.py    GNU General Public License v3.0 4 votes vote down vote up
def rotate(ra1, dec1, ra2, dec2, ra3, dec3):
    r"""Rotation matrix for rotation of (ra1, dec1) onto (ra2, dec2).
    The rotation is performed on (ra3, dec3).

    """
    def cross_matrix(x):
        r"""Calculate cross product matrix

        A[ij] = x_i * y_j - y_i * x_j

        """
        skv = np.roll(np.roll(np.diag(x.ravel()), 1, 1), -1, 0)
        return skv - skv.T

    ra1 = np.atleast_1d(ra1)
    dec1 = np.atleast_1d(dec1)
    ra2 = np.atleast_1d(ra2)
    dec2 = np.atleast_1d(dec2)
    ra3 = np.atleast_1d(ra3)
    dec3 = np.atleast_1d(dec3)

    assert(
        len(ra1) == len(dec1) == len(ra2) == len(dec2) == len(ra3) == len(dec3)
        )

    alpha = np.arccos(np.cos(ra2 - ra1) * np.cos(dec1) * np.cos(dec2)
                      + np.sin(dec1) * np.sin(dec2))
    vec1 = np.vstack([np.cos(ra1) * np.cos(dec1),
                      np.sin(ra1) * np.cos(dec1),
                      np.sin(dec1)]).T
    vec2 = np.vstack([np.cos(ra2) * np.cos(dec2),
                      np.sin(ra2) * np.cos(dec2),
                      np.sin(dec2)]).T
    vec3 = np.vstack([np.cos(ra3) * np.cos(dec3),
                      np.sin(ra3) * np.cos(dec3),
                      np.sin(dec3)]).T
    nvec = np.cross(vec1, vec2)
    norm = np.sqrt(np.sum(nvec**2, axis=1))
    nvec[norm > 0] /= norm[np.newaxis, norm > 0].T

    one = np.diagflat(np.ones(3))
    nTn = np.array([np.outer(nv, nv) for nv in nvec])
    nx = np.array([cross_matrix(nv) for nv in nvec])

    R = np.array([(1.-np.cos(a)) * nTn_i + np.cos(a) * one + np.sin(a) * nx_i
                  for a, nTn_i, nx_i in zip(alpha, nTn, nx)])
    vec = np.array([np.dot(R_i, vec_i.T) for R_i, vec_i in zip(R, vec3)])

    ra = np.arctan2(vec[:, 1], vec[:, 0])
    dec = np.arcsin(vec[:, 2])

    ra += np.where(ra < 0., 2. * np.pi, 0.)

    return ra, dec 
Example 78
Project: skylab   Author: coenders   File: grbllh.py    GNU General Public License v3.0 4 votes vote down vote up
def _select_events(self, src_ra, src_dec, scramble=False, inject=None):
        r"""Select events for log-likelihood evaluation.

        If `scramble` is `True`, `nbackground` (plus Poisson
        fluctuations) events are selected from the off-source time
        range. Otherwise, the on-source events ``data["on"]`` are
        selected.

        Note
        ----
        In the current implementation, the selection depends only on the
        on-source time range. Hence, `src_ra` and `src_dec` are ignored.

        """
        # We will chose new events, so it is time to clean the likelihood
        # model's cache.
        self.llh_model.reset()

        if scramble:
            N = self.random.poisson(self.nbackground)

            if N > 0:
                self._events = self.random.choice(self.data["off"], N)
                self._events["ra"] = self.random.uniform(0., 2.*np.pi, N)
            else:
                self._events = np.empty(0, dtype=self.data["off"].dtype)
        else:
            self._events = self.data["on"]

        if inject is not None:
            remove = np.logical_or(
                inject["sinDec"] < self.llh_model.sinDec_range[0],
                inject["sinDec"] > self.llh_model.sinDec_range[-1])

            if np.any(remove):
                inject = inject[np.logical_not(remove)]

            inject = numpy.lib.recfunctions.append_fields(
                inject, names="B", data=self.llh_model.background(inject),
                usemask=False)

            self._events = np.append(
                self._events, inject[list(self._events.dtype.names)])

        self._signal = self.llh_model.signal(src_ra, src_dec, self._events)

        # Method has to set number of events and number of selected
        # events. Here, both numbers are equal.
        self._nevents = self._events.size
        self._nselected = self._nevents 
Example 79
Project: LHMP   Author: hydrogo   File: wfdei_to_lumped_dataframe.py    GNU General Public License v3.0 4 votes vote down vote up
def PET(data, path_to_scheme):
    # Reference: http://www.fao.org/docrep/x0490e/x0490e07.htm
    # use with caution for latitudes out of range 0-67 degrees

    # Part 1. Avarage latitude calculation
    # read watershed scheme
    schema = pd.read_csv(path_to_scheme, usecols=[0, 1])
    # calculate mean watershed latitude
    # and convert it from degrees to radians
    lat = np.deg2rad(schema.Y.values.mean())

    # Part 2. Extraterrrestrial radiation calculation
    # set solar constant (in W m-2)
    Rsc = 1367
    # calculate day of the year array
    doy = np.array([i for i in range(1, 367)])
    # calculate solar declination dt (in radians)
    dt = 0.409 * np.sin(2 * np.pi / 365 * doy - 1.39)
    # calculate sunset hour angle (in radians)
    ws = np.arccos(-np.tan(lat) * np.tan(dt))
    # Calculate sunshine duration N (in hours)
    N = 24 / np.pi * ws
    # Calculate day angle j (in radians)
    j = 2 * np.pi / 365.25 * doy
    # Calculate relative distance to sun
    dr = 1.0 + 0.03344 * np.cos(j - 0.048869)
    # Calculate extraterrestrial radiation (J m-2 day-1)
    Re = Rsc * 86400 / np.pi * dr * (ws * np.sin(lat) * np.sin(dt)\
           + np.sin(ws) * np.cos(lat) * np.cos(dt))
    # convert from J m-2 day-1 to MJ m-2 day-1
    Re = Re/10**6

    # Part 3. Avearge daily temperatures calculation derived from long-term observations
    Ta = np.array([data.ix[data.index.dayofyear == x, 'Temp'].mean() for x in range(1, 367)])

    # Part 4. PET main equation by (Oudin et al., 2005)
    # lambda (latent heat flux const) = 2.45 MJ kg-1
    # ro (density of water const) = 1000 kg m-3
    # PE im m day -1 should be converted to mm/day (* 10**3)
    # PE = ( Re / (2.45*1000) ) * ( (Ta+5) /100 ) * 10**3
    # threshhold condition
    # if Ta+5>0 - use Oudin formula, else set to zero
    PE = np.where(Ta+5 > 0, ( Re / (2.45*1000) ) * ( (Ta+5) /100 )*10**3, 0)

    return PE


# source: my code from sandbox/WFDEI_to_dataframe.ipynb

# define files and variables been used 
Example 80
Project: DOTA_models   Author: ringringyi   File: nav_env.py    Apache License 2.0 4 votes vote down vote up
def _preprocess_for_task(self, seed):
    if self.task is None or self.task.seed != seed:
      rng = np.random.RandomState(seed)
      origin_loc = get_graph_origin_loc(rng, self.traversible)
      self.task = utils.Foo(seed=seed, origin_loc=origin_loc,
                            n_ori=self.task_params.n_ori)
      G = generate_graph(self.valid_fn_vec,
                                  self.task_params.step_size, self.task.n_ori,
                                  (0, 0, 0))
      gtG, nodes, nodes_to_id = convert_to_graph_tool(G)
      self.task.gtG = gtG
      self.task.nodes = nodes
      self.task.delta_theta = 2.0*np.pi/(self.task.n_ori*1.)
      self.task.nodes_to_id = nodes_to_id
      logging.info('Building %s, #V=%d, #E=%d', self.building_name,
                   self.task.nodes.shape[0], self.task.gtG.num_edges())

      if self.logdir is not None:
        write_traversible = cv2.applyColorMap(self.traversible.astype(np.uint8)*255, cv2.COLORMAP_JET)
        img_path = os.path.join(self.logdir,
                                '{:s}_{:d}_graph.png'.format(self.building_name,
                                                             seed))
        node_xyt = self.to_actual_xyt_vec(self.task.nodes)
        plt.set_cmap('jet');
        fig, ax = utils.subplot(plt, (1,1), (12,12))
        ax.plot(node_xyt[:,0], node_xyt[:,1], 'm.')
        ax.imshow(self.traversible, origin='lower');
        ax.set_axis_off(); ax.axis('equal');
        ax.set_title('{:s}, {:d}, {:d}'.format(self.building_name,
                                               self.task.nodes.shape[0],
                                               self.task.gtG.num_edges()))
        if self.room_dims is not None:
          for i, r in enumerate(self.room_dims['dims']*1):
            min_ = r[:3]*1
            max_ = r[3:]*1
            xmin, ymin, zmin = min_
            xmax, ymax, zmax = max_

            ax.plot([xmin, xmax, xmax, xmin, xmin],
                    [ymin, ymin, ymax, ymax, ymin], 'g')
        with fu.fopen(img_path, 'w') as f:
          fig.savefig(f, bbox_inches='tight', transparent=True, pad_inches=0)
        plt.close(fig)