Python numpy.cos() Examples

def periodic_hann(window_length):
  """Calculate a "periodic" Hann window.

  The classic Hann window is defined as a raised cosine that starts and
  ends on zero, and where every value appears twice, except the middle
  point for an odd-length window.  Matlab calls this a "symmetric" window
  and np.hanning() returns it.  However, for Fourier analysis, this
  actually represents just over one cycle of a period N-1 cosine, and
  thus is not compactly expressed on a length-N Fourier basis.  Instead,
  it's better to use a raised cosine that ends just before the final
  zero value - i.e. a complete cycle of a period-N cosine.  Matlab
  calls this a "periodic" window. This routine calculates it.

  Args:
    window_length: The number of points in the returned window.

  Returns:
    A 1D np.array containing the periodic hann window.
  """
  return 0.5 - (0.5 * np.cos(2 * np.pi / window_length *
                             np.arange(window_length))) 

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) 

def unit_vec(doa):
    """
    This function takes a 2D (phi) or 3D (phi,theta) polar coordinates
    and returns a unit vector in cartesian coordinates.

    :param doa: (ndarray) An (D-1)-by-N array where D is the dimension and
                N the number of vectors.

    :return: (ndarray) A D-by-N array of unit vectors (each column is a vector)
    """

    if doa.ndim != 1 and doa.ndim != 2:
        raise ValueError("DoA array should be 1D or 2D.")

    doa = np.array(doa)

    if doa.ndim == 0 or doa.ndim == 1:
        return np.array([np.cos(doa), np.sin(doa)])

    elif doa.ndim == 2 and doa.shape[0] == 1:
        return np.array([np.cos(doa[0]), np.sin(doa[0])])

    elif doa.ndim == 2 and doa.shape[0] == 2:
        s = np.sin(doa[1])
        return np.array([s * np.cos(doa[0]), s * np.sin(doa[0]), np.cos(doa[1])]) 

def rotation_matrix_3D(u, th):
    """
    rotation_matrix_3D(u, t) yields a 3D numpy matrix that rotates any vector about the axis u
    t radians counter-clockwise.
    """
    # normalize the axis:
    u = normalize(u)
    # We use the Euler-Rodrigues formula;
    # see https://en.wikipedia.org/wiki/Euler-Rodrigues_formula
    a = math.cos(0.5 * th)
    s = math.sin(0.5 * th)
    (b, c, d) = -s * u
    (a2, b2, c2, d2) = (a*a, b*b, c*c, d*d)
    (bc, ad, ac, ab, bd, cd) = (b*c, a*d, a*c, a*b, b*d, c*d)
    return np.array([[a2 + b2 - c2 - d2, 2*(bc + ad),         2*(bd - ac)],
                     [2*(bc - ad),       a2 + c2 - b2 - d2,   2*(cd + ab)],
                     [2*(bd + ac),       2*(cd - ab),         a2 + d2 - b2 - c2]]) 

def test_cmag(self):
        '''
        test_cmag() ensures that the neuropythy.vision cortical magnification function is working.
        '''
        import neuropythy.vision as vis
        logging.info('neuropythy: Testing areal cortical magnification...')
        dset = ny.data['benson_winawer_2018']
        sub = dset.subjects['S1202']
        hem = [sub.lh, sub.rh][np.random.randint(2)]
        cm = vis.areal_cmag(hem.midgray_surface, 'prf_',
                            mask=('inf-prf_visual_area', 1),
                            weight='prf_variance_explained')
        # cmag should get smaller in general
        ths = np.arange(0, 2*np.pi, np.pi/3)
        es = [0.5, 1, 2, 4]
        x = np.diff([np.mean(cm(e*np.cos(ths), e*np.sin(ths))) for e in es])
        self.assertTrue((x < 0).all()) 

def get_rotation_matrix(rotationVector, angle):
    """
    Calculate the rotation (3X3) matrix about an axis (rotationVector)
    by a rotation angle.

    :Parameters:
        #. rotationVector (list, tuple, numpy.ndarray): Rotation axis
           coordinates.
        #. angle (float): Rotation angle in rad.

    :Returns:
        #. rotationMatrix (numpy.ndarray): Computed (3X3) rotation matrix
    """
    angle = float(angle)
    axis = rotationVector/np.sqrt(np.dot(rotationVector , rotationVector))
    a = np.cos(angle/2)
    b,c,d = -axis*np.sin(angle/2.)
    return np.array( [ [a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
                       [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
                       [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c] ] , dtype = FLOAT_TYPE) 

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 

def sample(self):
    """Samples new points around some existing point.

    Removes the sampling base point and also stores the new jksampled points if
    they are far enough from all existing points.
    """
    active_point = self._active_list.pop()
    for _ in xrange(self._max_sample_size):
      # Generate random points near the current active_point between the radius
      random_radius = np.random.uniform(self._min_radius, 2 * self._min_radius)
      random_angle = np.random.uniform(0, 2 * math.pi)

      # The sampled 2D points near the active point
      sample = random_radius * np.array(
          [np.cos(random_angle), np.sin(random_angle)]) + active_point

      if not self._is_in_grid(sample):
        continue

      if self._is_close_to_existing_points(sample):
        continue

      self._active_list.append(sample)
      self._grid[self._point_to_index_1d(sample)] = sample 

def alive_bonus(self, z, pitch):
		if self.frame%30==0 and self.frame>100 and self.on_ground_frame_counter==0:
			target_xyz  = np.array(self.body_xyz)
			robot_speed = np.array(self.robot_body.speed())
			angle = self.np_random.uniform(low=-3.14, high=3.14)
			from_dist   = 4.0
			attack_speed   = self.np_random.uniform(low=20.0, high=30.0)  # speed 20..30 (* mass in cube.urdf = impulse)
			time_to_travel = from_dist / attack_speed
			target_xyz += robot_speed*time_to_travel  # predict future position at the moment the cube hits the robot
			position = [target_xyz[0] + from_dist*np.cos(angle),
				target_xyz[1] + from_dist*np.sin(angle),
				target_xyz[2] + 1.0]
			attack_speed_vector = target_xyz - np.array(position)
			attack_speed_vector *= attack_speed / np.linalg.norm(attack_speed_vector)
			attack_speed_vector += self.np_random.uniform(low=-1.0, high=+1.0, size=(3,))
			self.aggressive_cube.reset_position(position)
			self.aggressive_cube.reset_velocity(linearVelocity=attack_speed_vector)
		if z < 0.8:
			self.on_ground_frame_counter += 1
		elif self.on_ground_frame_counter > 0:
			self.on_ground_frame_counter -= 1
		# End episode if the robot can't get up in 170 frames, to save computation and decorrelate observations.
		self.frame += 1
		return self.potential_leak() if self.on_ground_frame_counter<170 else -1 

def calc_state(self):
		theta, self.theta_dot = self.central_joint.current_relative_position()
		self.gamma, self.gamma_dot = self.elbow_joint.current_relative_position()
		target_x, _ = self.jdict["target_x"].current_position()
		target_y, _ = self.jdict["target_y"].current_position()
		self.to_target_vec = np.array(self.fingertip.pose().xyz()) - np.array(self.target.pose().xyz())
		return np.array([
			target_x,
			target_y,
			self.to_target_vec[0],
			self.to_target_vec[1],
			np.cos(theta),
			np.sin(theta),
			self.theta_dot,
			self.gamma,
			self.gamma_dot,
		]) 

def random_quat(rand=None):
    """Return uniform random unit quaternion.
    rand: array like or None
        Three independent random variables that are uniformly distributed
        between 0 and 1.
    >>> q = random_quat()
    >>> np.allclose(1.0, vector_norm(q))
    True
    >>> q = random_quat(np.random.random(3))
    >>> q.shape
    (4,)
    """
    if rand is None:
        rand = np.random.rand(3)
    else:
        assert len(rand) == 3
    r1 = np.sqrt(1.0 - rand[0])
    r2 = np.sqrt(rand[0])
    pi2 = math.pi * 2.0
    t1 = pi2 * rand[1]
    t2 = pi2 * rand[2]
    return np.array(
        (np.sin(t1) * r1, np.cos(t1) * r1, np.sin(t2) * r2, np.cos(t2) * r2),
        dtype=np.float32,
    ) 

def test_variable(self, model):
        """ Test that variable updating works. """
        p1 = AnnualHarmonicSeriesParameter(model, 0.5, [0.25], [np.pi/4], is_variable=True)

        assert p1.double_size == 3
        assert p1.integer_size == 0

        new_var = np.array([0.6, 0.1, np.pi/2])
        p1.set_double_variables(new_var)
        np.testing.assert_allclose(p1.get_double_variables(), new_var)

        with pytest.raises(NotImplementedError):
            p1.set_integer_variables(np.arange(3, dtype=np.int32))

        with pytest.raises(NotImplementedError):
            p1.get_integer_variables()

        si = ScenarioIndex(0, np.array([0], dtype=np.int32))

        for ts in model.timestepper:
            doy = (ts.datetime.dayofyear - 1)/365
            np.testing.assert_allclose(p1.value(ts, si), 0.6 + 0.1*np.cos(doy*2*np.pi + np.pi/2)) 

def inverse_method(self,N,d):
        
        t = np.linspace(1e-3,0.999,N)
        f = np.log( t / (1 - t) )
        f = f/f[0]
        
        psi= np.pi*f
        cosPsi = np.cos(psi)
        sinTheta = ( np.abs(cosPsi) + (1-np.abs(cosPsi))*np.random.rand(len(cosPsi)))
        
        theta = np.arcsin(sinTheta)
        theta = np.pi-theta + (2*theta - np.pi)*np.round(np.random.rand(len(t)))
        cosPhi = cosPsi/sinTheta
        phi = np.arccos(cosPhi)*(-1)**np.round(np.random.rand(len(t)))
        
        coords = SkyCoord(phi*u.rad,(np.pi/2-theta)*u.rad,d*np.ones(len(phi))*u.pc)

        return coords 

def Jac(self, b):
        """Calculates determinant of the Jacobian transformation matrix to get
        the joint probability density of dMag and s
        
        Args:
            b (ndarray):
                Phase angles
                
        Returns:
            f (ndarray):
                Determinant of Jacobian transformation matrix
        
        """
        
        f = -2.5/(self.Phi(b)*np.log(10.0))*self.dPhi(b)*np.sin(b) - 5./np.log(10.0)*np.cos(b)
        
        return f 

def psi2c2c3(self, psi0):

        c2 = np.zeros(len(psi0))
        c3 = np.zeros(len(psi0))

        psi12 = np.sqrt(np.abs(psi0))
        pos = psi0 >= 0
        neg = psi0 < 0
        if np.any(pos):
            c2[pos] = (1 - np.cos(psi12[pos]))/psi0[pos]
            c3[pos] = (psi12[pos] - np.sin(psi12[pos]))/psi12[pos]**3.
        if any(neg):
            c2[neg] = (1 - np.cosh(psi12[neg]))/psi0[neg]
            c3[neg] = (np.sinh(psi12[neg]) - psi12[neg])/psi12[neg]**3.

        tmp = c2+c3 == 0
        if any(tmp):
            c2[tmp] = 1./2.
            c3[tmp] = 1./6.

        return c2,c3 

def calc_Phi(self, beta):
        """Calculate the phase function. Prototype method uses the Lambert phase 
        function from Sobolev 1975.
        
        Args:
            beta (astropy Quantity array):
                Planet phase angles at which the phase function is to be calculated,
                in units of rad
                
        Returns:
            Phi (ndarray):
                Planet phase function
        
        """
        
        beta = beta.to('rad').value
        Phi = (np.sin(beta) + (np.pi - beta)*np.cos(beta))/np.pi
        
        return Phi 

def equirect_facetype(h, w):
    '''
    0F 1R 2B 3L 4U 5D
    '''
    tp = np.roll(np.arange(4).repeat(w // 4)[None, :].repeat(h, 0), 3 * w // 8, 1)

    # Prepare ceil mask
    mask = np.zeros((h, w // 4), np.bool)
    idx = np.linspace(-np.pi, np.pi, w // 4) / 4
    idx = h // 2 - np.round(np.arctan(np.cos(idx)) * h / np.pi).astype(int)
    for i, j in enumerate(idx):
        mask[:j, i] = 1
    mask = np.roll(np.concatenate([mask] * 4, 1), 3 * w // 8, 1)

    tp[mask] = 4
    tp[np.flip(mask, 0)] = 5

    return tp.astype(np.int32) 

def rotate_point_cloud_by_angle_y(batch_data, rotation_angle):
	""" Rotate the point cloud along up direction with certain angle.
		Input:
		  BxNx3 array, original batch of point clouds
		Return:
		  BxNx3 array, rotated batch of point clouds
	"""
	rotated_data = np.zeros(batch_data.shape, dtype=np.float32)
	for k in range(batch_data.shape[0]):
		#rotation_angle = np.random.uniform() * 2 * np.pi
		cosval = np.cos(rotation_angle)
		sinval = np.sin(rotation_angle)
		rotation_matrix = np.array([[cosval, 0, sinval],
									[0, 1, 0],
									[-sinval, 0, cosval]])
		shape_pc = batch_data[k, ...]
		# rotated_data[k, ...] = np.dot(shape_pc.reshape((-1, 3)), rotation_matrix)
		rotated_data[k, ...] = np.dot(rotation_matrix, shape_pc.reshape((-1, 3)).T).T 		# Pre-Multiplication (changes done)
	return rotated_data 

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'})) 

def get_random_trajectory(self, seq_length):
    ''' Generate a random sequence of a MNIST digit '''
    canvas_size = self.image_size_ - self.digit_size_
    x = random.random()
    y = random.random()
    theta = random.random() * 2 * np.pi
    v_y = np.sin(theta)
    v_x = np.cos(theta)

    start_y = np.zeros(seq_length)
    start_x = np.zeros(seq_length)
    for i in range(seq_length):
      # Take a step along velocity.
      y += v_y * self.step_length_
      x += v_x * self.step_length_

      # Bounce off edges.
      if x <= 0:
        x = 0
        v_x = -v_x
      if x >= 1.0:
        x = 1.0
        v_x = -v_x
      if y <= 0:
        y = 0
        v_y = -v_y
      if y >= 1.0:
        y = 1.0
        v_y = -v_y
      start_y[i] = y
      start_x[i] = x

    # Scale to the size of the canvas.
    start_y = (canvas_size * start_y).astype(np.int32)
    start_x = (canvas_size * start_x).astype(np.int32)
    return start_y, start_x 

def align(self, marker, axis):
        '''
        Rotate the marker set around the given axis until it is aligned onto the given marker

        Parameters
        ----------
        marker : int or str
            the index or label of the marker onto which to align the set
        axis : int
            the axis around which the rotation happens
        '''

        index = self.key2ind(marker)
        axis = ['x','y','z'].index(axis) if isinstance(marker, (str, unicode)) else axis

        # swap the axis around which to rotate to last position
        Y = self.X
        if self.dim == 3:
            Y[axis,:], Y[2,:] = Y[2,:], Y[axis,:]

        # Rotate around z to align x-axis to second point
        theta = np.arctan2(Y[1,index],Y[0,index])
        c = np.cos(theta)
        s = np.sin(theta)
        H = np.array([[c, s],[-s, c]])
        Y[:2,:] = np.dot(H,Y[:2,:])

        if self.dim == 3:
            Y[axis,:], Y[2,:] = Y[2,:], Y[axis,:] 

def polar2cart(rho, phi):
    """
    convert from polar to cartesian coordinates
    :param rho: radius
    :param phi: azimuth
    :return:
    """
    x = rho * np.cos(phi)
    y = rho * np.sin(phi)
    return x, y 

def polar_distance(x1, x2):
    """
    Given two arrays of numbers x1 and x2, pairs the cells that are the
    closest and provides the pairing matrix index: x1(index(1,:)) should be as
    close as possible to x2(index(2,:)). The function outputs the average of 
    the absolute value of the differences abs(x1(index(1,:))-x2(index(2,:))).
    :param x1: vector 1
    :param x2: vector 2
    :return: d: minimum distance between d
             index: the permutation matrix
    """
    x1 = np.reshape(x1, (1, -1), order='F')
    x2 = np.reshape(x2, (1, -1), order='F')
    N1 = x1.size
    N2 = x2.size
    diffmat = np.arccos(np.cos(x1 - np.reshape(x2, (-1, 1), order='F')))
    min_N1_N2 = np.min([N1, N2])
    index = np.zeros((min_N1_N2, 2), dtype=int)
    if min_N1_N2 > 1:
        for k in range(min_N1_N2):
            d2 = np.min(diffmat, axis=0)
            index2 = np.argmin(diffmat, axis=0)
            index1 = np.argmin(d2)
            index2 = index2[index1]
            index[k, :] = [index1, index2]
            diffmat[index2, :] = float('inf')
            diffmat[:, index1] = float('inf')
        d = np.mean(np.arccos(np.cos(x1[:, index[:, 0]] - x2[:, index[:, 1]])))
    else:
        d = np.min(diffmat)
        index = np.argmin(diffmat)
        if N1 == 1:
            index = np.array([1, index])
        else:
            index = np.array([index, 1])
    return d, index 

def load_RSM(filename):
    om, tt, psd = xu.io.getxrdml_map(filename)
    om = np.deg2rad(om)
    tt = np.deg2rad(tt)
    wavelength = 1.54056

    q_y = (1 / wavelength) * (np.cos(tt) - np.cos(2 * om - tt))
    q_x = (1 / wavelength) * (np.sin(tt) - np.sin(2 * om - tt))

    xi = np.linspace(np.min(q_x), np.max(q_x), 100)
    yi = np.linspace(np.min(q_y), np.max(q_y), 100)
    psd[psd < 1] = 1
    data_grid = griddata(
        (q_x, q_y), psd, (xi[None, :], yi[:, None]), fill_value=1, method="cubic"
    )
    nx, ny = data_grid.shape

    range_values = [np.min(q_x), np.max(q_x), np.min(q_y), np.max(q_y)]
    output_data = (
        Panel(np.log(data_grid).reshape(nx, ny, 1), minor_axis=["RSM"])
        .transpose(2, 0, 1)
        .to_frame()
    )

    return range_values, output_data 

def _lb2RN_northcap(self, l, b):
        R = 100. + (90. - b) * np.sin(np.radians(l)) / 0.3
        N = 100. + (90. - b) * np.cos(np.radians(l)) / 0.3
        return np.round(R).astype('i4'), np.round(N).astype('i4') 

def _lb2RN_southcap(self, l, b):
        R = 100. + (90. + b) * np.sin(np.radians(l)) / 0.3
        N = 100. + (90. + b) * np.cos(np.radians(l)) / 0.3
        return np.round(R).astype('i4'), np.round(N).astype('i4') 

def f(x):
    return x * np.cos(np.pi * x) 

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 

def confounds_file(deriv_dir, preproc_file,
                   deriv_regressor_fname=deriv_regressor_fname):
    confounds_file = deriv_dir.ensure(deriv_regressor_fname)
    confound_dict = {}
    tp = nib.load(str(preproc_file)).shape[-1]
    ix = np.arange(tp)
    # csf
    confound_dict['csf'] = np.cos(2*np.pi*ix*(50/tp)) * 0.1
    # white matter
    confound_dict['white_matter'] = np.sin(2*np.pi*ix*(22/tp)) * 0.1
    # framewise_displacement
    confound_dict['framewise_displacement'] = np.random.random_sample(tp)
    confound_dict['framewise_displacement'][0] = np.nan
    # motion outliers
    for motion_outlier in range(0, 5):
        mo_name = 'motion_outlier0{}'.format(motion_outlier)
        confound_dict[mo_name] = np.zeros(tp)
        confound_dict[mo_name][motion_outlier] = 1
    # derivatives
    derive1 = [
        'csf_derivative1',
        'csf_derivative1_power2',
        'global_signal_derivative1_power2',
        'trans_x_derivative1',
        'trans_y_derivative1',
        'trans_z_derivative1',
        'trans_x_derivative1_power2',
        'trans_y_derivative1_power2',
        'trans_z_derivative1_power2',
    ]
    for d in derive1:
        confound_dict[d] = np.random.random_sample(tp)
        confound_dict[d][0] = np.nan

    # transformations
    for dir in ["trans_x", "trans_y", "trans_z"]:
        confound_dict[dir] = np.random.random_sample(tp)

    confounds_df = pd.DataFrame(confound_dict)
    confounds_df.to_csv(str(confounds_file), index=False, sep='\t', na_rep='n/a')
    return confounds_file 

def spherical_distance(pt0, pt1):
    '''
    spherical_distance(a, b) yields the angular distance between points a and b, both of which
      should be expressed in spherical coordinates as (longitude, latitude).
    If a and/or b are (2 x n) matrices, then the calculation is performed over all columns.
    The spherical_distance function uses the Haversine formula; accordingly it may suffer from
    rounding errors in the case of nearly antipodal points.
    '''
    dtheta = pt1[0] - pt0[0]
    dphi   = pt1[1] - pt0[1]
    a = np.sin(dphi/2)**2 + np.cos(pt0[1]) * np.cos(pt1[1]) * np.sin(dtheta/2)**2
    return 2 * np.arcsin(np.sqrt(a))