From b42059183591dcd873bd663a76f667f4041ef55c Mon Sep 17 00:00:00 2001 From: Kiyoshi Masui Date: Thu, 26 Jan 2012 11:51:53 -0500 Subject: [PATCH] Moved core/utils.py to utils/misc.py as it was screwing up some namespace stuff. Moved cubic interpolation module to utils. --- core/algebra.py | 6 +- core/base_data.py | 2 +- core/cubic_conv_interpolation.py | 403 ------------------------------- core/data_block.py | 2 +- core/test_algebra.py | 10 +- core/test_utils.py | 158 ------------ core/utils.py | 286 ---------------------- correlate/freq_slices.py | 3 +- noise/measure_noise.py | 4 +- plotting/call_summary_plot.py | 2 +- time_stream/rebin_time.py | 2 +- time_stream/split_bands.py | 2 +- 12 files changed, 17 insertions(+), 863 deletions(-) delete mode 100644 core/cubic_conv_interpolation.py delete mode 100644 core/test_utils.py delete mode 100644 core/utils.py diff --git a/core/algebra.py b/core/algebra.py index 04ca1587..52c7a41f 100644 --- a/core/algebra.py +++ b/core/algebra.py @@ -59,14 +59,14 @@ import sys import warnings from functools import wraps +import operator import scipy as sp import numpy as np import numpy.lib.format as npfor from numpy.lib.utils import safe_eval -import operator -import cubic_conv_interpolation as cci +import utils.cubic_conv_interpolation as cci # TODO: # when submitted as batch on Sunnyvale, the PYTHONPATH seems to get clobbered @@ -850,7 +850,7 @@ def slice_interpolate(self, axes, coord, kind='linear') : # Get the nodes needed and their weights. out = cci.interpolate(axes, self, pnt, x0, step_sizes) else: - # Get the contributing points and their wiehgts. + # Get the contributing points and their weights. points, weights = self.slice_interpolate_weights(axes, coord, kind) # Sum up the points q = points.shape[0] diff --git a/core/base_data.py b/core/base_data.py index e859744f..e5e3ec72 100644 --- a/core/base_data.py +++ b/core/base_data.py @@ -4,7 +4,7 @@ import numpy.ma as ma import kiyopy.custom_exceptions as ce -import utils +import utils.misc as utils from hist import History, merge_histories class BaseData(object) : diff --git a/core/cubic_conv_interpolation.py b/core/cubic_conv_interpolation.py deleted file mode 100644 index 4de2bccb..00000000 --- a/core/cubic_conv_interpolation.py +++ /dev/null @@ -1,403 +0,0 @@ -import numpy as np -#import matplotlib.pyplot as plt - -## If python complains of indentation errors, use tabs after def, if, for, -## etc... not spaces. - -def cross(set_list): - '''Given a list of sets, return the cross product.''' - - # By associativity of cross product, cross the first two sets together - # then cross that with the rest. The big conditional in the list - # comprehension is just to make sure that there are no nested lists - # in the final answer. - if len(set_list) == 1: - ans = [] - for elem in set_list[0]: - # In the 1D case, these elements are not lists. - if type(elem) == list: - ans.append(np.array(elem)) - else: - ans.append(np.array([elem])) - return ans - else: - A = set_list[0] - B = set_list[1] - cross_2 = [a+b if ((type(a) == list) and (type(b) == list)) else \ - (a+[b] if type(a) == list else ([a]+b if type(b) == list else \ - [a]+[b])) for a in A for b in B] - remaining = set_list[2:] - remaining.insert(0,cross_2) - return cross(remaining) - - -def point_list(node): - '''Return the list of lists containing the numbers, in order, of the - indeces that must be looked at in the interpolation. node is a - coordinate.''' - - ans = [] - for coord in node: - ans.append([coord-1,coord,coord+1,coord+2]) - return ans - - -def left_interpolation_node(pnt, x0, step_sizes): - '''Suppose pnt (list of coordinates in each dimension) is where the - sampled data is to be interpolated. Then pnt must be between two - consecutive interpolation nodes (in each dimension), x_j and x_(j+1). - Return x_j. The first point x0 and the step size in each dimension - is needed for x_j to be found. - Note: all parameters must have same dimensions. - Note: returns the node in indeces, not actual units.''' - - pnt = np.array(pnt, dtype=float) - left_node = np.floor((pnt - x0)/step_sizes) - return np.array(left_node,dtype=int) - - -def get_value(arr,node): - '''Get the value of n-dimensional array arr at n-coordinate index node. - Note node is a 1D array in pixels, not actual units.''' - - # Copy node since we needed it later and shouldn't mess with it. - node_copy = node.copy() - # The max index possible for each dimension of arr. - max_indeces = np.array(arr.shape) - 1 - # Initialize the bad_index. If it's stll -99 after looking through - # all the coordinates of pnt, then pnt is in arr. - bad_index = -99 - # True if one coordinate of pnt is under the min possible (0) and - # false if one coordinate of pnt is over the max possible. - under = False - # Find any index out of bounds (lower). - for i in range(len(node)): - if node[i] < 0: - bad_index = i - under = True - elif node[i] > max_indeces[i]: - bad_index = i - under = False - - if bad_index == -99: - # No bad indeces so just return the value of array at that point. - return arr[tuple(node)] - elif under: - # There is a bad index and it is -1. Use method similar to set of - # equations on page 1158 of article to find value. - node_copy[bad_index] += 1 - fst_node = node_copy.copy() - node_copy[bad_index] += 1 - snd_node = node_copy.copy() - node_copy[bad_index] += 1 - thd_node = node_copy.copy() - return 3.*get_value(arr,fst_node) - 3.*get_value(arr,snd_node) + \ - get_value(arr,thd_node) - else: - # There is a bad index and it is one more than the max possible. - # Note there cannot be any other cases. - node_copy[bad_index] -= 1 - fst_node = node_copy.copy() - node_copy[bad_index] -= 1 - snd_node = node_copy.copy() - node_copy[bad_index] -= 1 - thd_node = node_copy.copy() - return 3.*get_value(arr,fst_node) - 3.*get_value(arr,snd_node) + \ - get_value(arr,thd_node) - - -def u(s): - '''Return the cubic convolution interpolation kernel given the pixel - distance. See equation 15 from Robert G. Keys 1981 article.''' - - # Note: the s==0 and else cases are put in for correctness but - # these should never be reached since this function does not need - # to be called in those cases. - s = abs(s) - if 0 < s < 1: - return 3./2.*s**3 - 5./2.*s**2 + 1. - elif 1 < s < 2: - return -1./2.*s**3 + 5./2.*s**2 - 4.*s + 2. - elif s == 0: - return 1. - else: - return 0. - - -def u_product(pnt, node, x0, step_sizes): - '''Return the product of the "u" in each dimension. Want to get the - equivalent of the u part of the equation on pg 1157. pnt and node - are arrays containing the point we wish to interpolate and the points - with known values in the array, respectively. pnt is in actual units - while node is in indeces. x0 is the initial point in n-dimensions. - step_sizes is the distance between pixels in each dimension.''' - - step_sizes = np.array(step_sizes, dtype='float') - node_pixel = x0 + node*step_sizes - return np.product(map(u,(pnt-node_pixel)/step_sizes)) - - -def interpolate_weights(axes, pnt, x0, step_sizes): - '''Return the points and weights of all non-zero weighted points - needed in the cubic interpolation. 'axes' is the list of axes to be - interpolated over. Only the dimensions in axes will be looked at to - get the appropriate nodes and all other dimensions will be ignored. - See interpolate for description of other parameters.''' - - # Get only the needed axes out of the info we have to be able to - # call the code already written. - new_pnt = [] - new_x0 = [] - new_step_sizes = [] - - for i in range(len(axes)): - new_pnt.append(pnt[i]) - new_x0.append(x0[i]) - new_step_sizes.append(step_sizes[i]) - - new_pnt = np.array(new_pnt) - new_x0 = np.array(new_x0) - new_step_sizes = np.array(new_step_sizes) - - # Get the closest node to the "left". - left_node = left_interpolation_node(new_pnt, new_x0, new_step_sizes) - # Get the list of nodes that we need. - needed_nodes = cross(point_list(left_node)) - - # Get the corresponding list of weights. - weights = [] - for node in needed_nodes: - weights.append(u_product(new_pnt,node,new_x0,new_step_sizes)) - weights = np.array(weights,dtype=float) - - return np.array(needed_nodes),weights - - -def interpolate(axes, arr, pnt, x0, step_sizes): - '''Return cubic convolution interpolation at m-coordinate point pnt. - Note pnt is in the actual units of the axes. arr is the n-dimentional - array that contains the known discrete values to be interpolated. - x0 is a 1D m-element array that contains the first value of every - dimension in actual units. step_sizes is a 1D m-element array that - contains the step size between pixels in each dimension. - axes is the list of axes (length m) to interpolate on. If m==n, then - the estimated value of arr at point pnt is returned else, an interpolated - array with shape of the axes not in 'axes' is returned.''' - - # Check if axes is a subset of all the axes. - if len(axes) == arr.ndim: - return interpolate_one_step(arr,pnt,x0,step_sizes) - else: - # Get the points we want to look at and their associated weights. - points, weights = interpolate_weights(axes, pnt, x0, step_sizes) - - all_slices = [] - # Since a point no longer represent one node, we have to make - # the array it represents, assign the weights to each of these - # arrays, then sum up these arrays for the final answer. - for node in points: - # Get lists of the indices of each axis not in axes. - indices_lists = [] - sub_arr_shape = [] - for dim in range(arr.ndim): - if dim not in axes: - indices = range(arr.shape[dim]) - indices_lists.append(indices) - sub_arr_shape.append(arr.shape[dim]) - # Now put in the indices of node in the proper place - # based on its axis. - for j in range(len(axes)): - indices_lists.insert(axes[j], [node[j]]) - # Get the slice of the array we want to return. This contains - # the indices we need to look at, not the values yet. This also - # is flattened (1D) but can be made the right shape - # with sub_arr_shape. - slice = cross(indices_lists) - # Get the values at each index. Note we must use the "c" - # business since some of our node may be out of bounds. - values_slice = [] - for index in slice: - values_slice.append(get_value(arr,index)) - # Put out slice with the proper values into the proper shape. - values_slice = np.array(values_slice) - values_slice.shape = tuple(sub_arr_shape) - all_slices.append(values_slice) - - # Now that we have all the slices (and their weights from before), - # Multiply by their weights and add. - interpol_sum = 0. - for i in range(len(all_slices)): - interpol_sum += all_slices[i] * weights[i] - return interpol_sum - - -def interpolate_one_step(arr, pnt, x0, step_sizes): - '''Return cubic convolution interpolation at n-coordinate point pnt. - Note pnt is in the actual units of the axes. arr is the n-dimentional - array that contains the known discrete values to be interpolated. - x0 is a 1D n-element array that contains the first value of every - dimension in actual units. step_sizes is a 1D n-element array that - contains the step size between pixels in each dimension.''' - - # Get one of the nodes that surround the point we wish to interpolate. - left_node = left_interpolation_node(pnt, x0, step_sizes) - # Get the list of the 4 closest nodes (in each dimension) around the - # interpolation point. - needed_nodes = cross(point_list(left_node)) - # Do the n-dimensional equivalent to the formula on page 1157. - total = 0. - for node in needed_nodes: - c_node = get_value(arr,node) - u_node = u_product(pnt, node, x0, step_sizes) - total += c_node * u_node - return total - - -if __name__ == '__main__': - - # Make sure the kernel is right. - x = np.arange(-3,3,0.001); - ux = []; - for i in x: - ux.append(u(i)); -# plt.plot(x,ux); -# plt.plot(x,np.zeros(len(x)),'k'); -# plt.xlabel('s (pixel distance)'); plt.ylabel('kernel'); -# plt.title('Interpolation Kernel. Figure 1.(d)'); -# plt.show(); - - ## Some tests. - print '\nx**2 with x from -5 to 5' - print 'Point ..... Expected ..... Interpolation' - - # - # 1D quadratic function. - arr = np.array([25,16,9,4,1,0,1,4,9,16,25]); - x0 = np.array([-5]); - step_sizes = np.array([1]); - - # Some obvious points. - pnt = np.array([0]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '0 ..... 0 ..... ' + str(ans) - - pnt = np.array([-1]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '-1 ..... 1 ..... ' + str(ans) - - pnt = np.array([3]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '3 ..... 9 ..... ' + str(ans) - - # Some real tests. - pnt = np.array([0.5]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '0.5 ..... 0.25 ..... ' + str(ans) - - pnt = np.array([2.5]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '2.5 ..... 6.25 ..... ' + str(ans) - - pnt = np.array([-3.456]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '-3.456 ..... 11.943936 ..... ' + str(ans) - - pnt = np.array([4.34]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '4.34 ..... 18.8356 ..... ' + str(ans) - - pnt = np.array([-4.78]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '-4.78 ..... 22.8484 ..... ' + str(ans) - - - print '\nx**3 with x from -5 to 5' - print 'Point ..... Expected ..... Interpolation' - - # - # 1D cubic function. - arr = np.array([125,64,27,16,1,0,1,16,27,64,125]); - x0 = np.array([-5]); - step_sizes = np.array([1]); - - # Some obvious points. - pnt = np.array([0]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '0 ..... 0 ..... ' + str(ans) - - pnt = np.array([-3]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '-3 ..... 27 ..... ' + str(ans) - - # Real tests. Watch how the approximation is not as good for non quadratics - pnt = np.array([0.9]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '0.9 ..... 0.729 ..... ' + str(ans) - - pnt = np.array([1.17]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '1.17 ..... 1.601613 ..... ' + str(ans) - - pnt = np.array([3.65]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '3.65 ..... 48.627125 ..... ' + str(ans) - - pnt = np.array([-4.5]); - ans = interpolate([0], arr, pnt, x0, step_sizes); - print '-4.5 ..... 91.125 ..... ' + str(ans) - - - def f(x,y): - return x**2 + 4.*x*y + 6.*y**2 - - print '\nx**2 + 4xy + 6y**2 with x from -5 to 5, y from -5 to 5' - print 'Point ..... Expected ..... Interpolation' - - # - # 2D quadratic function. - # mesgrid does what I want backwards. - y,x = np.meshgrid(np.arange(-5,5),np.arange(-5,5)); - arr = f(x,y); - x0 = np.array([-5,-5]); - step_sizes = np.array([1,1]); - - # Simple cases - pnt = np.array([2,-1]); - ans = interpolate([0,1], arr, pnt, x0, step_sizes); - print '(2,-1) ..... 2 ..... ' + str(ans) - - pnt = np.array([-3,4]); - ans = interpolate([0,1], arr, pnt, x0, step_sizes); - print '(-3,4) ..... 57 ..... ' + str(ans) - - # Real tests. - pnt = np.array([3.7,-2.3]); - ans = interpolate([0,1], arr, pnt, x0, step_sizes); - print '(3.7,-2.3) ..... 11.39 ..... ' + str(ans) - - pnt = np.array([4.001,3.999]); - ans = interpolate([0,1], arr, pnt, x0, step_sizes); - print '(4.001,3.999) ..... 175.960003 ..... ' + str(ans) - - pnt = np.array([4.75,-4.23]); - ans = interpolate([0,1], arr, pnt, x0, step_sizes); - print '(4.75,-4.23) ..... 49.5499 ..... ' + str(ans) - - # - # Interpolation of sub axes. - print '\nSame setup as before but doing interpolation on one axes.' - print 'Using first axis and point 3.7' - print interpolate([0],arr,np.array([3.7]),np.array([-5]),np.array([1])); - - print '\nCheck by doing full interpolation on points (3.7,x)' - print 'where x are all the values from the second axis.' - for i in range(-5,5): - print interpolate([0,1],arr,np.array([3.7,i]),x0,step_sizes); - - # x_large = np.arange(0,1,0.01); - # f_large = x_large; - - # f = np.arange(0,1,0.05); - # plt.plot(f_large,f_large,'rx'); - # plt.plot(f,f,'g+'); - # plt.show(); diff --git a/core/data_block.py b/core/data_block.py index 25f3a1a5..2a13fc79 100644 --- a/core/data_block.py +++ b/core/data_block.py @@ -4,7 +4,7 @@ import numpy.ma as ma -import utils +import utils.misc as utils import base_data diff --git a/core/test_algebra.py b/core/test_algebra.py index ee64a524..dab13219 100755 --- a/core/test_algebra.py +++ b/core/test_algebra.py @@ -668,8 +668,8 @@ def test_slice_interpolate_nearest(self): def test_slice_interpolate_cubic(self) : # Construct a 3D array that is a quadratic function. - data = sp.arange(60, dtype=float) - data.shape = (5, 4, 3) + data = sp.arange(140, dtype=float) + data.shape = (5, 4, 7) vect = algebra.info_array(data) vect = algebra.vect_array(vect,axis_names=('freq', 'a', 'b')) @@ -678,12 +678,12 @@ def test_slice_interpolate_cubic(self) : a.shape = (5, 1, 1) b = sp.arange(-1,3)**2 b.shape = (1, 4, 1) - c = sp.arange(-1,2)**2 - c.shape = (1, 1, 3) + c = sp.arange(-1,6)**2 + c.shape = (1, 1, 7) v[:,:,:] = a + b + c v.set_axis_info('freq', 0, 1) v.set_axis_info('a', 1, 1) - v.set_axis_info('b', 0, 1) + v.set_axis_info('b', 2, 1) #### First test the weights. diff --git a/core/test_utils.py b/core/test_utils.py deleted file mode 100644 index be8c2de8..00000000 --- a/core/test_utils.py +++ /dev/null @@ -1,158 +0,0 @@ -"""Unit tests for utils module.""" - -import unittest - -import scipy as sp -from numpy import random -import scipy.linalg as linalg - -import utils - -class TestElAz2RaDec_LST(unittest.TestCase) : - - def test_known_cases(self) : - el = [90, 21, 60] - az = [12, 180, 90] - lst = [0, 26.*86400./360., 0] - lat = [13, 0, 0] - - ra, dec = utils.elaz2radec_lst(el, az, lst, lat) - - self.assertAlmostEqual(ra[0], 0) - self.assertAlmostEqual(dec[0], lat[0]) - - self.assertAlmostEqual(ra[1], 26, 4) - self.assertAlmostEqual(dec[1], -90+el[1]) - - self.assertAlmostEqual(ra[2], 360-30) - self.assertAlmostEqual(dec[2], 0) - -class TestElAz2RaDec(unittest.TestCase) : - - def test_works(self) : - UT = '2000-01-01T12:00:00.00' - el = 90. - az = 50.5 - ra, dec = utils.elaz2radecGBT(el, az, UT) - self.assertAlmostEqual(dec, 38.43312, 1) #GBT Latitude - el = 38.43312 - az = 0. - ra, dec = utils.elaz2radecGBT(el, az, UT) - self.assertAlmostEqual(dec, 90, 1) - el = 90 - 38.43312 - az = 180. - ra, dec = utils.elaz2radecGBT(el, az, UT) - self.assertAlmostEqual(dec, 0, 1) - -class TestMakeMapGrid(unittest.TestCase) : - - def test_basic(self) : - centre = [80., 20.] - shape = (11, 11) - spacing = 1. - grid_ra, grid_dec = utils.mk_map_grid(centre, shape, spacing) - self.assertEqual(sp.shape(grid_ra), shape) - self.assertEqual(sp.shape(grid_dec), shape) - self.assertTrue(sp.allclose(grid_ra[:,5], centre[0])) - self.assertTrue(sp.allclose(grid_dec[5,:], centre[1])) - self.assertAlmostEqual(sp.amax(grid_dec), 25) - self.assertAlmostEqual(sp.amin(grid_dec), 15) - self.assertTrue(sp.allclose(grid_ra[:,7], centre[0] + - 2/sp.cos(centre[1]*sp.pi/180.))) - -class TestPolInt2String(unittest.TestCase): - - def test_all_functionality(self): - self.assertEqual(utils.polint2str(1), 'I') - self.assertEqual(utils.polint2str(2), 'Q') - self.assertEqual(utils.polint2str(3), 'U') - self.assertEqual(utils.polint2str(4), 'V') - self.assertEqual(utils.polint2str(-5), 'XX') - self.assertEqual(utils.polint2str(-6), 'YY') - self.assertEqual(utils.polint2str(-7), 'XY') - self.assertEqual(utils.polint2str(-8), 'YX') - self.assertEqual(utils.polint2str(-1), 'RR') - self.assertEqual(utils.polint2str(-2), 'LL') - self.assertEqual(utils.polint2str(-3), 'RL') - self.assertEqual(utils.polint2str(-4), 'LR') - self.assertRaises(ValueError, utils.polint2str, 7) - -class TestAmpFit(unittest.TestCase): - - def test_uncorrelated_noscatter(self): - data = sp.arange(10, dtype=float) - theory = data/2.0 - C = sp.identity(10) - a, s = utils.ampfit(data, C, theory) - self.assertAlmostEqual(a, 2) - - def test_uncorrelated_noscatter_error(self): - data = sp.arange(10, dtype=float) + 1.0 - amp = 5.0 - theory = data/amp - C = sp.diag((sp.arange(10, dtype=float) + 1.0)**2) - a, s = utils.ampfit(data, C, theory) - self.assertAlmostEqual(a, amp) - self.assertAlmostEqual(s/a, 1.0/sp.sqrt(10)) - - def test_uncorrelated_scatter_error(self): - n = 1000 - data = sp.arange(n, dtype=float) + 1.0 - amp = 5.0 - theory = data/amp - C = sp.diag((sp.arange(n, dtype=float) + 1.0)**2) - data += random.normal(size=n)*data - a, s = utils.ampfit(data, C, theory) - self.assertTrue(sp.allclose(a, amp, rtol=5.0/sp.sqrt(n), atol=0)) - # Expect the next line to fail 1/100 trials. - self.assertFalse(sp.allclose(a, amp, rtol=0.01/sp.sqrt(n), atol=0)) - self.assertAlmostEqual(s, amp/sp.sqrt(1000)) - - def test_correlated_scatter(self) : - n = 50 - r = (sp.arange(n, dtype=float) + 10.0*n)/10.0*n - data = sp.sin(sp.arange(n)) * r - amp = 25.0 - theory = data/amp - # Generate correlated matrix. - C = random.rand(n, n) # [0, 1) - # Raise to high power to make values near 1 rare. - C = (C**10) * 0.2 - C = (C + C.T)/2.0 - C += sp.identity(n) - C *= r[:, None]/2.0 - C *= r[None, :]/2.0 - # Generate random numbers in diagonal frame. - h, R = linalg.eigh(C) - self.assertTrue(sp.alltrue(h>0)) - rand_vals = random.normal(size=n)*sp.sqrt(h) - # Rotate back. - data += sp.dot(R.T, rand_vals) - a, s = utils.ampfit(data, C, theory) - self.assertTrue(sp.allclose(a, amp, atol=5.0*s, rtol=0)) - # Expect the next line to fail 1/100 trials. - self.assertFalse(sp.allclose(a, amp, atol=0.01*s, rtol=0)) - -class TestFloatTime(unittest.TestCase): - - def test_near_epoch(self): - UT = "2000-01-01T00:00:54.51" - self.assertAlmostEqual(54.51, utils.time2float(UT)) - self.assertEqual(utils.float2time(54.51), UT) - - def test_full_circle(self): - seconds = 451732642.56 - UT = utils.float2time(seconds) - seconds2 = utils.time2float(UT) - self.assertEqual(UT, utils.float2time(seconds2)) - self.assertAlmostEqual(seconds, seconds2) - - def test_full_circle_vectorized(self): - seconds = sp.array([[34252672.12, 623421.65], [23464.1, 644656784.56]]) - UT = utils.float2time(seconds) - seconds2 = utils.time2float(UT) - self.assertTrue(sp.all(UT == utils.float2time(seconds2))) - self.assertTrue(sp.allclose(seconds, seconds2)) - -if __name__ == '__main__' : - unittest.main() diff --git a/core/utils.py b/core/utils.py deleted file mode 100644 index 6e6cb9e7..00000000 --- a/core/utils.py +++ /dev/null @@ -1,286 +0,0 @@ -"""Some general utilities that are useful.""" - -import datetime -import time - -import scipy as sp -from scipy.interpolate import interp1d -from scipy import linalg -import ephem - -def elaz2radec_lst(el, az, lst, lat = 38.43312) : - """DO NOT USE THIS ROUTINE FOR ANTHING THAT NEEDS TO BE RIGHT. IT DOES NOT - CORRECT FOR PRECESSION. - - Calculates the Ra and Dec from elavation, aximuth, LST and Latitude. - - This function is vectorized with numpy so should be fast. Standart numpy - broadcasting should also work. - - All angles in degrees, lst in seconds. Latitude defaults to GBT. - """ - - # Convert everything to radians. - el = sp.radians(el) - az = sp.radians(az) - lst = sp.array(lst, dtype = float)*2*sp.pi/86400 - lat = sp.radians(lat) - # Calculate dec. - dec = sp.arcsin(sp.sin(el)*sp.sin(lat) + - sp.cos(el)*sp.cos(lat)*sp.cos(az)) - # Calculate the hour angle - ha = sp.arccos((sp.sin(el) - sp.sin(lat)*sp.sin(dec)) / - (sp.cos(lat)*sp.cos(dec))) - ra = sp.degrees(lst - ha) % 360 - - return ra, sp.degrees(dec) - -def elaz2radecGBT(el, az, UT) : - """Calculates the Ra and Dec from the elevation, azimuth and UT for an - observer at GBT. - - All input should be formated to correspond to the data in a GBT fits file. - El, and Az in degrees and UT a string like in the GBT DATE-OBS field. - - Largely copied from Kevin's code. - """ - - GBT = ephem.Observer() - GBT.long = '-79:50:23.4' - GBT.lat = '38:25:59.23' - GBT.pressure = 0 # no refraction correction. - GBT.temp = 0 - - UT_wholesec, partial_sec = UT.split('.', 1) - time_obj = time.strptime(UT_wholesec, "%Y-%m-%dT%H:%M:%S") - UT_reformated = time.strftime("%Y/%m/%d %H:%M:%S", time_obj) - GBT.date = UT_reformated + "." + partial_sec - - el_r = el*sp.pi/180.0 - az_r = az*sp.pi/180.0 - ra, dec = GBT.radec_of(az_r,el_r) - - return ra*180.0/sp.pi, dec*180.0/sp.pi - -def LSTatGBT(UT) : - """Calculates the LST from the UT of an observer at GBT. - - All input should be formated to correspond to the data in a GBT fits file. - UT a string like in the GBT DATE-OBS field. - - Largely copied from Kevin's code. - """ - - GBT = ephem.Observer() - GBT.long = '-79:50:23.4' - GBT.lat = '38:25:59.23' - GBT.pressure = 0 # no refraction correction. - GBT.temp = 0 - - UT_wholesec, partial_sec = UT.split('.', 1) - time_obj = time.strptime(UT_wholesec, "%Y-%m-%dT%H:%M:%S") - UT_reformated = time.strftime("%Y/%m/%d %H:%M:%S", time_obj) - GBT.date = UT_reformated + "." + partial_sec - - LST = GBT.sidereal_time() #IN format xx:xx:xx.xx ? - - return LST*180.0/sp.pi - -def time2float(UT) : - """Calculates float seconds from a time string. - - Convert a time string in format %Y-%m-%dT%H:%M:%S.partial to a float number - of seconds ignoring many corrections (such as leap seconds). However it - should be internally consistant, and it should be safe to use this method - to convert to a float, do some operations, then use `float2time` to convert - back. Works for array inputs.""" - - if not isinstance(UT, sp.ndarray) : - UT = sp.array([UT]) - single_time = True - else: - single_time = False - - time_array = sp.empty(UT.shape, dtype=float) - - for ii in xrange(UT.size) : - UT_wholesec, partial_sec = UT.flat[ii].split('.', 1) - to = datetime.datetime(*time.strptime(UT_wholesec, - "%Y-%m-%dT%H:%M:%S")[:6]) - epoch_start = datetime.datetime(2000, 1, 1) - td = to - epoch_start - td_seconds = ((td.microseconds + (td.seconds + td.days * 24 * 3600) * - 10**6) / 10**6) - time_array.flat[ii] = td_seconds + float("0." + partial_sec) - if single_time : - return time_array[0] - else : - return time_array - -def float2time(t) : - """Does the reverse operation of `time2float`.""" - - if not isinstance(t, sp.ndarray) : - t = sp.array([t]) - single_time = True - else: - single_time = False - - time_str_array = sp.empty(t.shape, dtype='S22') - - for ii in range(t.size) : - td = datetime.timedelta(seconds=t.flat[ii]) - epoch_start = datetime.datetime(2000, 1, 1) - time_obj = epoch_start + td - - time_str = time_obj.strftime("%Y-%m-%dT%H:%M:%S.%f") - # Truncate the fraction seconds to hundredths. - time_str = time_str[:22] - time_str_array.flat[ii] = time_str - if single_time : - return time_str_array[0] - else : - return time_str_array - -def mk_map_grid(centre, shape, spacing) : - """Make a grid of coordinates in Ra and Dec. - - This function accepts a field centre (tuple 2 floats), map-shape (tuple 2 - ints) and pixel spacing (float) for a map. It returns two arrays with - the specified map shape, and meshed like sp.meshgrid. However, the spacing - in Ra is converted to real degrees (devided by cos(dec) at field centre). - - All units in degrees, values are pixel centres. - """ - - dec = centre[1] + spacing*sp.arange(-(shape[1]-1.)/2., shape[1]/2.) - ra = centre[0] + (spacing/sp.cos(centre[1]*sp.pi/180.) * - sp.arange(-(shape[0]-1.)/2., shape[0]/2.)) - - grid_ra, grid_dec = sp.meshgrid(ra, dec) - - return grid_ra, grid_dec - -def get_beam(freq) : - """Get the GBT beam width at a frequency (or an array of frequencies). - - This is currently pretty rough and only uses on scans worth of data. - """ - - # This data is pretty rough. Just cut and pasted from one scan, not - # averaged. - beam_data = [0.316148488246, 0.306805630985, 0.293729620792, - 0.281176247549, 0.270856788455, 0.26745856078, - 0.258910010848, 0.249188429031] - freq_data = sp.array([695, 725, 755, 785, 815, 845, 875, 905], dtype=float) - freq_data *= 1.0e6 - f = interp1d(freq_data, beam_data, bounds_error=False, fill_value = -1) - b = f(freq) - b[b<0] = 0.316148488246 - return b - -def polint2str(pol_int) : - """Convert an interger representing a polarization to a representing the - polarization. - - This is based on the SDfits convention that I pulled from: - https://safe.nrao.edu/wiki/bin/view/Main/SdfitsDetails - - Here are the return values based on the passed integer. - - RR -1 - LL -2 - RL -3 - LR -4 - XX -5 - YY -6 - XY -7 - YX -8 - I 1 - Q 2 - U 3 - V 4 - Otherwise raises a ValueError. - """ - - if pol_int == -1 : - return 'RR' - elif pol_int == -2 : - return 'LL' - elif pol_int == -3 : - return 'RL' - elif pol_int == -4 : - return 'LR' - elif pol_int == -5 : - return 'XX' - elif pol_int == -6 : - return 'YY' - elif pol_int == -7 : - return 'XY' - elif pol_int == -8 : - return 'YX' - elif pol_int == 1 : - return 'I' - elif pol_int == 2 : - return 'Q' - elif pol_int == 3 : - return 'U' - elif pol_int == 4 : - return 'V' - else : - raise ValueError("Polarization integer must be in range(-8, 5) and " - "nonzero") - -def ampfit(data, covariance, theory): - """Fits the amplitude of the theory curve to the data. - - Finds `amp` such that `amp`*`theory` is the best fit to `data`. - - Returns - ------- - amp : float - Fitted amplitude. - errir : float - Error on fitted amplitude. - """ - - data = sp.asarray(data) - covariance = sp.asarray(covariance) - theory = sp.asarray(theory) - # Sanity check inputs. - if len(data.shape) != 1: - raise ValueError("`data` must be a 1D vector.") - n = len(data) - if data.shape != theory.shape: - raise ValueError("`theory` must be the same shape as `data`.") - if covariance.shape != (n,n): - msg = "`covariance` must be a square matrix compatible with data." - raise ValueError(msg) - # Linear fit for the amplitude. Formulas July 24, 2011 of Kiyo's notebook. - covariance_inverse = linalg.inv(covariance) - weighted_data = sp.dot(covariance_inverse, data) - amp = sp.dot(theory, weighted_data) - normalization = sp.dot(covariance_inverse, theory) - normalization = sp.dot(theory, normalization) - amp /= normalization - # Calculate the Error. - error = sp.sqrt(1/normalization) - - return amp, error - -def rebin_1D(array, reduce=4, axis=-1): - shape = array.shape - if axis < 0: - axis = array.ndim + axis - array.shape = shape[:axis] + (shape[axis]//reduce, reduce) + shape[axis+1:] - out = sp.zeros(shape[:axis] + (shape[axis]//reduce,) + shape[axis+1:], - dtype=array.dtype) - for ii in range(reduce): - index = [slice(None),]*array.ndim - index[axis+1] = ii - index = tuple(index) - out += array[index] - out /= reduce - array.shape = shape - return out - diff --git a/correlate/freq_slices.py b/correlate/freq_slices.py index 8256c798..6a650f61 100644 --- a/correlate/freq_slices.py +++ b/correlate/freq_slices.py @@ -26,7 +26,8 @@ from kiyopy import parse_ini import kiyopy.utils import kiyopy.custom_exceptions as ce -from core import utils, algebra +from core import algebra +import utils.misc as utils from core import handythread as ht import itertools import map.tools diff --git a/noise/measure_noise.py b/noise/measure_noise.py index 852c0f80..94d773c5 100644 --- a/noise/measure_noise.py +++ b/noise/measure_noise.py @@ -15,7 +15,7 @@ import kiyopy.utils import kiyopy.custom_exceptions as ce import core.fitsGBT -import core.utils +import utils.misc from map.constants import T_infinity, T_small # XXX @@ -140,7 +140,7 @@ def process_file(self, file_middle, Pipe=None) : ii * n_bands_proc + jj) a.set_xlabel("frequency (Hz)") a.set_ylabel("power (K^2)") - pol_str = core.utils.polint2str(pols[jj]) + pol_str = utils.misc.polint2str(pols[jj]) band = band_centres[ii] a.set_title("Band %dMHz, polarization " % band + pol_str) diff --git a/plotting/call_summary_plot.py b/plotting/call_summary_plot.py index 5dfd2fe6..9661b9dd 100644 --- a/plotting/call_summary_plot.py +++ b/plotting/call_summary_plot.py @@ -7,7 +7,7 @@ from correlate.freq_slices import MapPair import multiprocessing import shelve -from core import utils +import utils.misc as utils import numpy as np # note that notes are purely human-readable and the keys do not mean anything diff --git a/time_stream/rebin_time.py b/time_stream/rebin_time.py index 0b7653bf..615e2e13 100644 --- a/time_stream/rebin_time.py +++ b/time_stream/rebin_time.py @@ -9,7 +9,7 @@ import base_single import kiyopy.custom_exceptions as ce -from core import utils +import utils.misc as utils class RebinTime(base_single.BaseSingle) : diff --git a/time_stream/split_bands.py b/time_stream/split_bands.py index f0c90bb1..c7aacad5 100644 --- a/time_stream/split_bands.py +++ b/time_stream/split_bands.py @@ -11,7 +11,7 @@ import base_single import kiyopy.custom_exceptions as ce -from core import utils +import utils.misc as utils from core import data_block