diff --git a/numpy/lib/__init__.py b/numpy/lib/__init__.py index 17557605def5..b84b8f61e319 100644 --- a/numpy/lib/__init__.py +++ b/numpy/lib/__init__.py @@ -19,6 +19,7 @@ import math from arrayterator import * from arraypad import * +from neighbor import * __all__ = ['emath','math'] __all__ += type_check.__all__ @@ -34,6 +35,7 @@ __all__ += arraysetops.__all__ __all__ += npyio.__all__ __all__ += financial.__all__ +__all__ += ['neighbor'] from numpy.testing import Tester test = Tester().test diff --git a/numpy/lib/neighbor.py b/numpy/lib/neighbor.py new file mode 100644 index 000000000000..2951b0e9756d --- /dev/null +++ b/numpy/lib/neighbor.py @@ -0,0 +1,305 @@ + +""" +The neighbor module contains a function to calculate a new value +based on neighboring values. +""" + +import numpy as np + +__all__ = ['neighbor'] + + + +################################################################################ +# Public functions + + +def shapeintersect(intersect_arrays): + """ + Returns the maximally overlaping intersection of input arrays. + + After alignment the intersection is determined by expanding outward from + the origin/alignment elements to the common maximal extent in all + arrays in each dimension. + + Parameters + ---------- + intersect_arrays : sequence of two item sequences + A sequence of (array, origin) pairs. The origin index represent the + alignment element between all of the arrays. + + Returns + ------- + shapeintersect : list of ndarry + List of length equal to the number of input arrays. Each returned + array is the maximal intesected shape from the input arrays, in turn. + + Examples + -------- + >>> a = (((7, 3, 5, 1, 5), (1,)), ((4, 5, 6), (1,))) + >>> shapeintersect(a) + [array([7, 3, 5]), array([4, 5, 6])] + + >>> a = (((7, 3, 5, 1, 5), (1,)), ((4, 5, 6), (2,))) + >>> shapeintersect(a) + [array([7, 3]), array([5, 6])] + """ + + # Bring in just one array to test rank. Since all arrays must have the + # same rank, it doesn't matter which one I use. + tarr = np.array(intersect_arrays[0][0]) + tori = np.array(intersect_arrays[0][1], dtype='int') + # Might be a better way. Move through the arrays to collect the minimum + # lengths between the origins and edges of the arrays. + for narr, norigin in intersect_arrays: + oarr = np.array(narr) + oori = np.array(norigin, dtype='int') + if tarr.ndim != oarr.ndim: + raise ValueError(( + 'All input arrays should have the same number ' + 'of dimensions.')) + if len(tori) != len(oori): + raise ValueError(( + 'All origins should have the same length.')) + if len(oori) != oarr.ndim: + raise ValueError(( + 'The length of the origin should be equal to the number ' + 'of dimensions in array.')) + # Test to make sure origin is within array by taking advantage of index + # error. Error checking is all the following statement is used for. + a = oarr[oori] + try: + lindex = np.minimum(lindex, oori) + except NameError: + lindex = np.copy(oori) + try: + uindex = np.minimum(uindex, np.array(oarr.shape) - oori) + except NameError: + uindex = np.array(oarr.shape) - oori + + # Now lets make the real calculation of the slice indices + accum_arrays = [] + for narr, norigin in intersect_arrays: + oarr = np.array(narr) + oori = np.array(norigin, dtype='int') + lslice = oori - lindex + uslice = oori + uindex + slc = [slice(l, u) for l, u in zip(lslice, uslice)] + accum_arrays.append(oarr[slc]) + return accum_arrays + + +def neighbor(array, weight, func, mode='same', weight_origin=None, + pad='reflect', constant_values=0, end_values=0, stat_length=None, + reflect_type='even', **kwargs): + """ + Calculates a new value from a collection of neighboring elements. + + Parameters + ---------- + array : array_like + Input array + weight : array_like + Array that selects and weights values to be part of the neighborhood. + The center of the `weight` array is aligned with each element in + `array`. The center is exactly determined if the axis is odd, and + approximated as the nearest lower indexed element if even. The center + can be specified with the `weight_origin` keyword. + + Identical rank to `array`. + func : function, optional + Function `func` must accept a single rank 1 ndarray as an argument (and + optionally `kwargs`) and return a scalar. + mode : {'same', 'valid'}, optional + 'same' Mode same is the default and returns the same shape as + `array`. Boundary effects are still visible, though padding + can be used to minimize. + 'valid' Mode valid returns will return an array of + shape = `array.shape` - `weight.shape` + 1 + The element is only returned where `array` and + `weight` overlap completely. + weight_origin : array_like + The origin is the element that is aligned with each element in `array`. + The rank of `weight_origin` must equal the rank of the `weight` array, + and in each dimension must be less than the length of the corresponding + dimension in `array`. + pad : {str, function, None}, optional + Default is 'reflect'. + + See np.pad documentation for details. + + One of the following string values or a user supplied function. + + 'constant' Pads with a constant value. + 'edge' Pads with the edge values of array. + 'linear_ramp' Pads with the linear ramp between end_value and the + array edge value. + 'maximum' Pads with the maximum value of all or part of the + vector along each axis. + 'mean' Pads with the mean value of all or part of the + vector along each axis. + 'median' Pads with the median value of all or part of the + vector along each axis. + 'minimum' Pads with the minimum value of all or part of the + vector along each axis. + 'reflect' Pads with the reflection of the vector mirrored on + the first and last values of the vector along each + axis. + 'symmetric' Pads with the reflection of the vector mirrored + along the edge of the array. + 'wrap' Pads with the wrap of the vector along the axis. + The first values are used to pad the end and the + end values are used to pad the beginning. + Padding function, see Notes in the np.pad function. + constant_values : {sequence, int} + Optional for `pad` of 'constant'. Defaults to 0. + See np.pad documentation for details. + end_values : {sequence, int} + Optional for `pad` of 'linear_ramp'. Defaults to 0. + See np.pad documentation for details. + stat_length : {sequence, int} + Optional for `pad` in ['minimum', 'maximum', 'mean', 'median']. + Defaults to using the entire vector along each axis. + See np.pad documentation for details. + reflect_type : {sequence, int} + Optional for `pad` in ['reflect', 'symmetric']. Defaults to 'even'. + See np.pad documentation for details. + kwargs : varies + Any additional keyword arguments that should be passed to `func`. + + Returns + ------- + neighbor : ndarray + Array with values calculated from neighbors. Rank will be equivalent + to rank of `array`, but shape will depend on `mode` keyword. + + Examples + -------- + >>> a = [7, 3, 5, 1, 5] + >>> np.lib.neighbor(a, [1]*3, np.mean, pad=None) + array([5, 5, 3, 3, 3]) + + >>> np.lib.neighbor(a, [1]*3, np.min) + array([3, 3, 1, 1, 1]) + + >>> np.lib.neighbor(a, [1]*3, np.max, mode='valid') + array([7, 5, 5]) + + >>> np.lib.neighbor(a, [1]*3, np.max) + array([7, 7, 5, 5, 5]) + + The 'neighbor' and 'convolution' concepts are related. + See the next two examples. + + >>> np.convolve(a, [2]*3, mode='same') + array([20, 30, 18, 22, 12]) + + To emulate convolution with 'neighbor' use the np.sum function. + + >>> np.lib.neighbor(a, [2]*3, np.sum, pad='constant') + array([20, 30, 18, 22, 12]) + + Example of 'neighbor' support for n-dimensional arrays: + + >>> a = np.arange(25).reshape((5, 5)) + >>> kernel = np.ones((3, 3)) + >>> np.lib.neighbor(a, kernel, np.min) + array([[ 0, 0, 1, 2, 3], + [ 0, 0, 1, 2, 3], + [ 5, 5, 6, 7, 8], + [10, 10, 11, 12, 13], + [15, 15, 16, 17, 18]]) + + Create your own function to use in the neighborhood. + >>> def counter(arr): + ... # here 'arr' will always be a 1-D array. + ... return len(arr) + + >>> np.lib.neighbor(a, kernel, counter, pad=None) + array([[4, 6, 6, 6, 4], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [6, 9, 9, 9, 6], + [4, 6, 6, 6, 4]]) + """ + + narray = np.array(array) + weight_array = np.array(weight) + + # Error checks + if narray.ndim != weight_array.ndim: + raise ValueError(( + 'The "array" (rank %d) and "weight" (rank %d) arguments ' + 'must have the same rank.') % + (narray.ndim, weight_array.ndim)) + + if mode not in ['same', 'valid']: + raise ValueError( + 'The "mode" keyword must be one of ["same", "valid"]') + + # Find the origin of the weight array + if weight_origin: + weight_center = np.array(weight_origin) + if np.any(np.array(weight_center.shape) > np.array(weight_array.shape)): + raise ValueError(( + 'The "weight_origin" keyword must less than "weight" array ' + 'in each dimension. Instead %s !< %s') % + (weight_center.shape, weight_array.shape)) + if np.any(weight_center < 0): + raise ValueError(( + 'All values in array_like "weight_origin" keyword must be' + 'positive. Instead have %s') % + (weight_center)) + else: + weight_center = (np.array(weight_array.shape) - 1)/2 + + if mode == 'valid': + pad = None + + if pad: + # Pad width dependent on weight. + pad_width = zip(weight_center, + np.array(weight_array.shape) - weight_center - 1) + # Different modes only allow certain keywords. + padkw = {} + if pad == 'constant': + padkw['constant_values'] = constant_values + if pad == 'linear_ramp': + padkw['linear_ramp'] = end_values + if pad in ['mean', 'minimum', 'maximum', 'median']: + padkw['stat_length'] = stat_length + if pad in ['reflect', 'symmetric']: + padkw['reflect_type'] = reflect_type + narray = np.pad(narray, pad_width, mode=pad, **padkw) + + if pad or mode == 'valid': + array_offset = weight_center + else: + array_offset = np.zeros(weight_center.shape) + + if mode == 'valid': + # Need to convert element numbering to slice notation, by subtracting 1 + # for the lower index and adding 1 to the upper index. Then need to + # make sure the the lower index is greater than zero, and the upper + # index is less than shape of the input array. + lindex = np.array(weight_center) - 1 + lindex[lindex < 0] = 0 + uindex = np.array(weight_array.shape) - lindex + uindex = np.where(uindex > weight_array.shape, + np.array(weight_array.shape) + 1, uindex) + valid_slices = tuple(slice(l, u) for l, u in zip(lindex, uindex)) + else: + valid_slices = (slice(None)) + + out_array = np.empty_like(np.array(array)[valid_slices]) + for mindex in np.ndindex(*out_array.shape): + offset_index = array_offset + np.array(mindex) + tarray, tweight = shapeintersect(((narray, offset_index), + (weight_array, weight_center))) + tarray = tarray*tweight + if kwargs: + out_array[mindex] = func(np.extract(np.isfinite(tweight), tarray), + kwargs) + else: + out_array[mindex] = func(np.extract(np.isfinite(tweight), tarray)) + return out_array diff --git a/numpy/lib/tests/test_neighbor.py b/numpy/lib/tests/test_neighbor.py new file mode 100644 index 000000000000..43a865f52a12 --- /dev/null +++ b/numpy/lib/tests/test_neighbor.py @@ -0,0 +1,278 @@ +''' +Tests for the pad functions. +''' + +from numpy.testing import TestCase, run_module_suite, assert_array_equal +from numpy.testing import assert_raises, assert_array_almost_equal +import numpy as np +from numpy.lib import neighbor + + +class TestStatistic(TestCase): + def test_check_two_d(self): + a = np.arange(100).astype('f').reshape((10,10)) + kern = np.ones((3,3)) + a = neighbor(a, kern, np.mean, pad=None) + b = np.array( + [[ 5.5, 6., 7., 8., 9., 10., 11., 12., 13., 13.5], + [ 10.5, 11., 12., 13., 14., 15., 16., 17., 18., 18.5], + [ 20.5, 21., 22., 23., 24., 25., 26., 27., 28., 28.5], + [ 30.5, 31., 32., 33., 34., 35., 36., 37., 38., 38.5], + [ 40.5, 41., 42., 43., 44., 45., 46., 47., 48., 48.5], + [ 50.5, 51., 52., 53., 54., 55., 56., 57., 58., 58.5], + [ 60.5, 61., 62., 63., 64., 65., 66., 67., 68., 68.5], + [ 70.5, 71., 72., 73., 74., 75., 76., 77., 78., 78.5], + [ 80.5, 81., 82., 83., 84., 85., 86., 87., 88., 88.5], + [ 85.5, 86., 87., 88., 89., 90., 91., 92., 93., 93.5]]) + assert_array_equal(a, b) + + def test_check_one_d(self): + a = np.arange(10).astype('f') + kern = [1]*3 + a = neighbor(a, kern, np.max, pad=None) + b = np.array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 9.]) + assert_array_equal(a, b) + + def test_check_collapsed_axis_kern(self): + a = np.arange(25).astype('f').reshape((5,5)) + kern = [1]*3 + kern = np.array(kern).reshape((3,1)) + a = neighbor(a, kern, np.mean, pad=None) + b = np.array( + [[ 2.5, 3.5, 4.5, 5.5, 6.5], + [ 5. , 6. , 7. , 8. , 9. ], + [ 10. , 11. , 12. , 13. , 14. ], + [ 15. , 16. , 17. , 18. , 19. ], + [ 17.5, 18.5, 19.5, 20.5, 21.5]]) + assert_array_equal(a, b) + + def test_check_1_d_convolve(self): + a = np.arange(25).astype('f') + wei = [0.1, 0.2, 0.5, 0.2, 0.1] + b = neighbor(a, wei, np.sum, pad=None) + c = np.convolve(a, wei, mode='same') + assert_array_almost_equal(b, c) + + def test_check_2_d_convolve(self): + a = np.arange(25).astype('f') + a = np.reshape(a, (5,5)) + wei = [ + [ 2.5, 5.5, 6.5], + [ 5. , 8. , 9. ], + [ 17.5, 20.5, 21.5]] + flwei = np.fliplr(np.flipud(wei)) + b = neighbor(a, flwei, np.sum, pad='constant') + # scipy.signal.convolve2d(a, wei) + c = np.array([[ 47.5, 101. , 137.5, 174. , 160.5], + [ 170. , 339. , 435. , 531. , 452. ], + [ 465. , 819. , 915. , 1011. , 807. ], + [ 760. , 1299. , 1395. , 1491. , 1162. ], + [ 852.5, 1406. , 1487.5, 1569. , 1175.5]]) + assert_array_almost_equal(b, c) + + def test_ndimage_convolve_example_1(self): + a = np.array([[1, 2, 0, 0], + [5, 3, 0, 4], + [0, 0, 0, 7], + [9, 3, 0, 0]]) + b = np.array([[1,1,1],[1,1,0],[1,0,0]]) + # to get the same answer as ndimage.convolve b must be flipped + b = np.fliplr(np.flipud(b)) + # ndimage.convolve(a, b, mode='constant', cval=0.0) + # scipy.signal.convolve2d(a, b, mode='valid') + neigh = neighbor(a, b, np.sum, pad='constant') + c = np.array([[11, 10, 7, 4], + [10, 3, 11, 11], + [15, 12, 14, 7], + [12, 3, 7, 0]]) + assert_array_almost_equal(neigh, c) + + def test_ndimage_convolve_example_2(self): + a = np.array([[1, 2, 0, 0], + [5, 3, 0, 4], + [0, 0, 0, 7], + [9, 3, 0, 0]]) + b = np.array([[1,1,1],[1,1,0],[1,0,0]]) + # to get the same answer as ndimage.convolve b must be flipped + b = np.fliplr(np.flipud(b)) + # ndimage.convolve(a, b, mode='constant', cval=1.0) + neigh = neighbor(a, b, np.sum, pad='constant', constant_values=1) + c = np.array([[13, 11, 8, 7], + [11, 3, 11, 14], + [16, 12, 14, 10], + [15, 6, 10, 5]]) + assert_array_almost_equal(neigh, c) + + def test_ndimage_convolve_example_reflect_1(self): + a = np.array([[2, 0, 0], + [1, 0, 0], + [0, 0, 0]]) + b = np.array([[0,1,0],[0,1,0],[0,1,0]]) + # to get the same answer as ndimage.convolve b must be flipped + b = np.fliplr(np.flipud(b)) + # 'symmetric' in np.pad/np.neighbor it what ndimage.convolve calls + # 'reflect' + neigh = neighbor(a, b, np.sum, pad='symmetric') + c = np.array([[5, 0, 0], + [3, 0, 0], + [1, 0, 0]]) + assert_array_almost_equal(neigh, c) + + def test_ndimage_convolve_example_reflect_2(self): + a = np.array([[2, 0, 1], + [1, 0, 0], + [0, 0, 0]]) + b = np.array([[0,1,0], + [0,1,0], + [0,1,0], + [0,1,0], + [0,1,0]]) + # to get the same answer as ndimage.convolve b must be flipped + b = np.fliplr(np.flipud(b)) + # 'edge' in np.pad/np.neighbor it what ndimage.convolve calls + # 'nearest' + neigh = neighbor(a, b, np.sum, pad='edge') + c = np.array([[7, 0, 3], + [5, 0, 2], + [3, 0, 1]]) + assert_array_almost_equal(neigh, c) + + + def test_check_3_d(self): + a = np.arange(360) + a = np.reshape(a, (12, 5, 6)) + kern = np.ones((3,3,3)) + a = neighbor(a, kern, np.sum, pad=None) + b = np.array( + [[[ 148, 228, 240, 252, 264, 180], + [ 258, 396, 414, 432, 450, 306], + [ 330, 504, 522, 540, 558, 378], + [ 402, 612, 630, 648, 666, 450], + [ 292, 444, 456, 468, 480, 324]], + + [[ 402, 612, 630, 648, 666, 450], + [ 657, 999, 1026, 1053, 1080, 729], + [ 765, 1161, 1188, 1215, 1242, 837], + [ 873, 1323, 1350, 1377, 1404, 945], + [ 618, 936, 954, 972, 990, 666]], + + [[ 762, 1152, 1170, 1188, 1206, 810], + [1197, 1809, 1836, 1863, 1890, 1269], + [1305, 1971, 1998, 2025, 2052, 1377], + [1413, 2133, 2160, 2187, 2214, 1485], + [ 978, 1476, 1494, 1512, 1530, 1026]], + + [[1122, 1692, 1710, 1728, 1746, 1170], + [1737, 2619, 2646, 2673, 2700, 1809], + [1845, 2781, 2808, 2835, 2862, 1917], + [1953, 2943, 2970, 2997, 3024, 2025], + [1338, 2016, 2034, 2052, 2070, 1386]], + + [[1482, 2232, 2250, 2268, 2286, 1530], + [2277, 3429, 3456, 3483, 3510, 2349], + [2385, 3591, 3618, 3645, 3672, 2457], + [2493, 3753, 3780, 3807, 3834, 2565], + [1698, 2556, 2574, 2592, 2610, 1746]], + + [[1842, 2772, 2790, 2808, 2826, 1890], + [2817, 4239, 4266, 4293, 4320, 2889], + [2925, 4401, 4428, 4455, 4482, 2997], + [3033, 4563, 4590, 4617, 4644, 3105], + [2058, 3096, 3114, 3132, 3150, 2106]], + + [[2202, 3312, 3330, 3348, 3366, 2250], + [3357, 5049, 5076, 5103, 5130, 3429], + [3465, 5211, 5238, 5265, 5292, 3537], + [3573, 5373, 5400, 5427, 5454, 3645], + [2418, 3636, 3654, 3672, 3690, 2466]], + + [[2562, 3852, 3870, 3888, 3906, 2610], + [3897, 5859, 5886, 5913, 5940, 3969], + [4005, 6021, 6048, 6075, 6102, 4077], + [4113, 6183, 6210, 6237, 6264, 4185], + [2778, 4176, 4194, 4212, 4230, 2826]], + + [[2922, 4392, 4410, 4428, 4446, 2970], + [4437, 6669, 6696, 6723, 6750, 4509], + [4545, 6831, 6858, 6885, 6912, 4617], + [4653, 6993, 7020, 7047, 7074, 4725], + [3138, 4716, 4734, 4752, 4770, 3186]], + + [[3282, 4932, 4950, 4968, 4986, 3330], + [4977, 7479, 7506, 7533, 7560, 5049], + [5085, 7641, 7668, 7695, 7722, 5157], + [5193, 7803, 7830, 7857, 7884, 5265], + [3498, 5256, 5274, 5292, 5310, 3546]], + + [[3642, 5472, 5490, 5508, 5526, 3690], + [5517, 8289, 8316, 8343, 8370, 5589], + [5625, 8451, 8478, 8505, 8532, 5697], + [5733, 8613, 8640, 8667, 8694, 5805], + [3858, 5796, 5814, 5832, 5850, 3906]], + + [[2548, 3828, 3840, 3852, 3864, 2580], + [3858, 5796, 5814, 5832, 5850, 3906], + [3930, 5904, 5922, 5940, 5958, 3978], + [4002, 6012, 6030, 6048, 6066, 4050], + [2692, 4044, 4056, 4068, 4080, 2724]]]) + assert_array_equal(a, b) + + def test_lower_edge_origin(self): + a = np.array([[1, 2, 3], + [4, 5, 6], + [7, 8, 9]]) + b = np.array([[1,2,3],[4,5,6],[7,8,9]]) + neigh = neighbor(a, b, np.sum, pad='symmetric', weight_origin=[0,0]) + c = np.array([[285, 312, 306], + [348, 375, 369], + [294, 321, 315]]) + assert_array_almost_equal(neigh, c) + + def test_upper_edge_origin(self): + a = np.array([[1, 2, 3], + [4, 5, 6], + [7, 8, 9]]) + b = np.array([[1,2,3],[4,5,6],[7,8,9]]) + neigh = neighbor(a, b, np.sum, pad=None, weight_origin=[2,2]) + c = np.array([[ 9, 26, 50], + [ 42, 94, 154], + [ 90, 186, 285]]) + assert_array_almost_equal(neigh, c) + + def test_valid(self): + a = np.array([[1, 2, 3], + [4, 5, 6], + [7, 8, 9]]) + b = np.array([[1,2,3],[4,5,6],[7,8,9]]) + neigh = neighbor(a, b, np.sum, pad=None, weight_origin=[2,2], mode='valid') + c = np.array([[285]]) + assert_array_almost_equal(neigh, c) + + def test_valid_01(self): + a = np.arange(10) + neigh = neighbor(a, [1]*3, np.max, mode='valid') + c = np.array([2,3,4]) + assert_array_almost_equal(neigh, c) + + def test_failing_docstring(self): + a = np.array([7, 3, 5, 1, 5]) + neigh = np.lib.neighbor(a, [1]*3, np.max, mode='valid') + c = np.array([7, 5, 5]) + assert_array_almost_equal(neigh, c) + +class ValueError1(TestCase): + def test_check_same_rank(self): + arr = np.arange(30) + arr = np.reshape(arr, (6, 5)) + kern = [True]*3 + assert_raises(ValueError, neighbor, arr, kern, np.sum) + + def test_check_kern_odd(self): + arr = np.arange(30) + arr = np.reshape(arr, (6, 5)) + kern = [True]*4 + assert_raises(ValueError, neighbor, arr, kern, np.sum) + + +if __name__ == "__main__": + run_module_suite()