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 8 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: 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 4
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 5
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 6
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 7
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 8
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 9
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 10
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 11
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 12
Project: skylab   Author: coenders   File: ps_model.py    GNU General Public License v3.0 6 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 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]