Source code for irfpy.util.nptools

"""In house numpy helpers.

.. codeauthor:: Yoshifumi Futaana

The aim of this module is to extend the numpy functionality.

.. autosummary::

    lingrids
    loggrids
    min_where
    max_where
    shrink
    squeeze_right
    squeeze_left
    isascending
    isdescending
    guess_upper
    guess_lower
    guess_interval
    padlin
    sort_periodic
    unjump_periodic
    span
    broadcast
    broadcast_r
    rtp
    lonlat2xyz
    trim_epsilon
    nearest_index
    nearest_index_right
    nearest_index_left
"""
import bisect as _bisect

from irfpy.util import exception as ex
import numpy as np


[docs]def min_where(arr, cond): """ Conditional minimum. Minimum is calculated but under the specific conditions. For example, it can be used when you want to know the minimum of the given array but above zero. >>> arr = np.array([[0, 0, 7], [3, 2, 4]]) >>> print(min_where(arr, arr > 0)) 2 :param arr: Numpy array. :param cond: Condition. Should be shape shape as ``arr``. :returns: Return the minimum of arr, but where cond is satisfied. If all the cond is ``False``, ``nan`` is returned. >>> print(min_where(arr, arr > 1000)) nan """ if arr.shape != cond.shape: raise ValueError('Shape mismatch.') arr_cond = np.ma.masked_array(arr, mask=(~cond)) minimum_value = arr_cond.min() try: minimum_value = float(minimum_value.filled(np.nan)) # In case the value minimum_value is masked. except AttributeError as e: pass return minimum_value
[docs]def max_where(arr, cond): """ Conditional maximum, similar to :func:`min_where`. :param arr: Numpy array. :param cond: Condition. Should be shape shape as ``arr``. :returns: Return the minimum of arr, but where cond is satisfied. If all the cond is ``False``, ``nan`` is returned. >>> arr = np.array([[0, 0, 7], [3, 2, 4]]) >>> print(max_where(arr, arr < 4)) 3 """ if arr.shape != cond.shape: raise ValueError('Shape mismatch.') arr_cond = np.ma.masked_array(arr, mask=(~cond)) max_value = arr_cond.max() try: max_value = float(max_value.filled(np.nan)) # In case the value max_value is masked. except AttributeError as e: pass # This is nominal case. Usually max_value is float or int, not masked scalar. return max_value
[docs]def shrink(array, newshape, operation='mean'): """ The given array is shrinked in size with given operation. :param array: Numpy array. Multidimensional. :param newshape: New shape. The new shape must have the same number of array dimensions. Number of each new axis should be the divisor of the corresponding array dimension. :keyword operation: 'sum', 'mean', 'max', 'min', 'median' .. note:: Keep in mind, "median" is not multidimensional value, but done dimension-by-dimension from the last dimension. Namely, the operation is like ``array2d.median(axis=1).median(axis=0)``. >>> print(shrink(np.arange(24).reshape(4, 6), (2, 3), operation='mean')) # doctest: +NORMALIZE_WHITESPACE [[ 3.5 5.5 7.5] [15.5 17.5 19.5]] >>> print(shrink(np.arange(24).reshape(4, 6), (2, 3), operation='max')) [[ 7 9 11] [19 21 23]] """ if operation not in ['sum', 'mean', 'max', 'min', 'median']: raise ValueError("Operation ``{}`` is not supported.".format(operation)) if array.ndim != len(newshape): raise ValueError("Dimensions should be kept (Original={}/New={})".format( array.shape, newshape)) shape_num = [(d, c // d) for d, c in zip(newshape, array.shape)] shape_num_index = [shape for num in shape_num for shape in num] try: array = array.reshape(shape_num_index) except ValueError as e: raise ValueError('Given shape is not compatible for shrinking. Original={}/New={}'.format( str(array.shape), str(newshape))) # Operation starts for each dimension. for i in range(len(newshape)): if operation == "sum": array = array.sum(-1 * (i + 1)) elif operation == "mean": array = array.mean(-1 * (i + 1)) elif operation == "max": array = array.max(axis=-1 * (i + 1)) elif operation == "min": array = array.min(axis=-1 * (i + 1)) elif operation == "median": array = np.ma.median(array, axis=-1 * (i + 1)) else: raise ValueError("Operation {} not supported.".format(operation)) return array
[docs]def squeeze_right(arr): """ Similar to :func:`numpy:numpy.squeeze`, the given array is squeezed but only from right axes. >>> arr = np.random.uniform(size=[1, 3, 5, 4, 2, 1]) >>> np.squeeze(arr).shape (3, 5, 4, 2) >>> squeeze_right(arr).shape (1, 3, 5, 4, 2) >>> print((arr[..., 0] == squeeze_right(arr)).all()) True """ shape = arr.shape newshape = [] sq = True for num in shape[::-1]: if sq: if num != 1: sq = False newshape.insert(0, num) else: pass # Ignore the dimension else: newshape.insert(0, num) return arr.reshape(newshape)
[docs]def squeeze_left(arr): """ Similar to :func:`numpy:numpy.squeeze`, the given array is squeezed but only from left axes. >>> arr = np.random.uniform(size=[1, 3, 5, 4, 2, 1]) >>> np.squeeze(arr).shape (3, 5, 4, 2) >>> squeeze_left(arr).shape (3, 5, 4, 2, 1) >>> print((arr[0] == squeeze_left(arr)).all()) True """ shape = arr.shape newshape = [] sq = True for num in shape: if sq: if num != 1: sq = False newshape.append(num) else: pass # Ignore the dimension else: newshape.append(num) return arr.reshape(newshape)
[docs]def padlin(vector, pad_width, iaxis, kwargs): """ A helper function to use :func:`numpy:numpy.pad`. Use :func:`numpy:numpy.pad` function with padding of constant inclination over the whole range. >>> a = [1.0, 3, 7, 5, 2] # 5 element array. Average inclication is 0.25 = (2 -1 ) / 4 >>> b = np.pad(a, (2, 3), mode=padlin) # Specify the function, padlin, as a mode for np.pad. >>> print(b) # doctest: +NORMALIZE_WHITESPACE [0.5 0.75 1. 3. 7. 5. 2. 2.25 2.5 2.75] In this example, the padding before is 2 element, which increments 0.25. At the end, 3 elements are added, using 0.25 increments after the last value of 2. """ incl = (vector[-pad_width[1] - 1] - vector[pad_width[0]]) / (len(vector) - pad_width[0] - pad_width[1] - 1) vector[:pad_width[0]] = -np.arange(pad_width[0], 0, -1) * incl + vector[pad_width[0]] vector[-pad_width[1]:] = (np.arange(pad_width[1]) + 1) * incl + vector[-pad_width[1] - 1] return vector
[docs]def isascending(array1d, allow_equal=True): """ Check if the array is ascending order or not :param array1d: Array to be checked. :keyword allow_equal: Adjacent equal values are treated as ascending. >>> arr = [1, 2, 3, 4, 5] >>> print(isascending(arr)) True >>> arr = [1, 2, -3, 4, 5] >>> print(isascending(arr)) False >>> arr = [1, 2, 2, 5] >>> print(isascending(arr)) True >>> arr = [1, 2, 2, 5] >>> print(isascending(arr, allow_equal=False)) False """ array1d = np.copy(array1d) if len(array1d) == 0: return False elif len(array1d) == 1: return True darr = array1d[1:] - array1d[:-1] if allow_equal: return darr.min() >= 0 else: return darr.min() > 0
[docs]def isdescending(array1d, allow_equal=True): """ Check if the array is ascending order or not :param array1d: Array to be checked. :keyword allow_equal: Adjacent equal values are treated as ascending. >>> arr = [5, 4, 3, 2, 1] >>> print(isdescending(arr)) True >>> arr = [5, 4, -3, 2, 1] >>> print(isdescending(arr)) False >>> arr = [5, 2, 2, 1] >>> print(isdescending(arr)) True >>> arr = [5, 2, 2, 1] >>> print(isdescending(arr, allow_equal=False)) False """ if len(array1d) == 0: return False elif len(array1d) == 1: return True array1d = np.copy(array1d)[::-1] return isascending(array1d, allow_equal=allow_equal)
def _guess_upper_but_last(array1d): """ Return the array of upper bound, except for the last one. """ arr_guess = np.zeros_like(array1d) - 1e20 # Set very strange value as default. arr_guess[:-1] = (array1d[:-1] + array1d[1:]) / 2.0 return arr_guess def _guess_lower_but_first(array1d): arr_guess = np.zeros_like(array1d) - 1e20 # Set very strange value as default. arr_guess[1:] = (array1d[:-1] + array1d[1:]) / 2.0 return arr_guess
[docs]def guess_interval(array1d, log=False, method="median"): """ From the given 1-D array, a guessed interval array is returned. The returned value is :func:`guess_upper` - :func:`guess_lower`. """ return guess_upper(array1d, log=log, method=method) - guess_lower(array1d, log=log, method=method)
[docs]def guess_upper(array1d, log=False, method="median"): """ From the given 1D array, a guessed upper bound array is returned. For the given ascending array, the upper bound is calculate from the mean with the next element. For example, the array (2, 3, 5, 7, 11, 13, 17), the upper bound array is (2.5, 4, 6, 9, 12, 15, *19*). All the elements but the last one is the mean with the next element. The last one is guessed with the given *method*; by default, the median of the difference of the other element. >>> arr = np.array([2., 3, 5, 7, 11, 13, 17]) >>> arr_upper = guess_upper(arr) >>> print(arr_upper) # doctest: +NORMALIZE_WHITESPACE [ 2.5 4. 6. 9. 12. 15. 18. ] The array is assumed to be sorted ascendingly, but no check is done. :param array1d: Ascending np.array. :param log: Set *True* to calculate this in the log space. :return: 1D array. You can get this functionality in log space. >>> arr = [1, 10, 100, 1000, 10000] In this case, the upper value is the "Geometric mean". >>> arr_upper = guess_upper(arr, log=True) >>> print(arr_upper) # doctest: +NORMALIZE_WHITESPACE [3.16227766e+00 3.16227766e+01 3.16227766e+02 3.16227766e+03 3.16227766e+04] """ arr = np.copy(array1d) if log: arr = np.log(arr) ### Get upper bound, but for the last element. arr_guess = _guess_upper_but_last(arr) meth = method.lower() if meth in ('median',): diffarr = arr_guess - arr meandiff = np.median(diffarr[:-1]) arr_guess[-1] = arr[-1] + meandiff ### elif meth in ('average',): # will follow in the future.. else: raise ex.IrfpyException('Only method="median" is supported yet') if log: arr_guess = np.exp(arr_guess) return arr_guess
[docs]def guess_lower(array1d, log=False, method='median'): """ From the given 1D array, a guessed lower bound array is returned. See the similar function :func:`guess_upper` for details. >>> arr = [2, 3, 5, 7, 11, 13, 17] >>> lower = guess_lower(arr) >>> print(lower) # doctest: +NORMALIZE_WHITESPACE [ 1. 2.5 4. 6. 9. 12. 15. ] >>> arr = [1, 10, 100, 1000, 10000] >>> arr_lower = guess_lower(arr, log=True) >>> print(arr_lower) # doctest: +NORMALIZE_WHITESPACE [3.16227766e-01 3.16227766e+00 3.16227766e+01 3.16227766e+02 3.16227766e+03] """ arr = np.copy(array1d) if log: arr = np.log(arr) ### Get lower bound, but for the first element. arr_guess = _guess_lower_but_first(arr) meth = method.lower() if meth in ('median',): diffarr = arr_guess - arr meandiff = np.median(diffarr[1:]) arr_guess[0] = arr[0] + meandiff ### elif meth in ('average',): # will follow in the future.. else: raise ex.IrfpyException('Only method="median" is supported yet') if log: arr_guess = np.exp(arr_guess) return arr_guess
[docs]def guess_bounds(array1d, log=False, method='median'): """ From the given 1D array, a guessed bound array is returned. From the given 1D array, a guessed bound array is returned. Bound array is an array specifying the bounds of the elements of the original array. For example, suppose we have an array ``[0, 1, 2, 3]``. The corresponding bound array is ``[-0.5, 0.5, 1.5, 2.5, 3.5]``. The 1st element (``0``) of the original array's bound is ``-0.5`` and ``0.5``. The returned array must have 1 elements more than that of the original array. :param array1d: Array-like. :param log: If the estimation to be done in ``log`` space, use this option as *True*. :param method: Method to guess. Currently ``median`` is only supported. :return: A new array, with a size of (N+1), where N is the size of the original array (``array1d``). >>> original_array = [0, 1, 2, 3] >>> bound_array = guess_bounds(original_array) >>> print(bound_array) [-0.5, 0.5, 1.5, 2.5, 3.5] """ arr_u = guess_upper(array1d, log=log, method=method) arr_l = guess_lower(array1d, log=log, method=method) arr = arr_l.tolist() + [arr_u[-1]] return arr
[docs]def sort_periodic(array1d, period): """ From an array in an periodic system, ascending array is obtained. :param array1d: An array. :param period: The period. :return: An array. If you have an array [0, 100, 20, 120, 40] in an ascending system with periodicity of 180. In this case, the returned array is [0, 100, *200*, *300*, *400*]. The original element of "20" is smaller than "100", the one element before. This is interpreted as "the value 20 is ascending from 100, but exceed the upper bound so that the value is decrimented by period=180. Thus, the ascending value should be 20+180=200". >>> arr = [0, 100, 20, 120, 40] >>> ascend_arr = sort_periodic(arr, 180) >>> print(ascend_arr) # doctest: +NORMALIZE_WHITESPACE [ 0 100 200 300 400] """ asc_arr = np.copy(array1d) narr = asc_arr.shape[0] for i in range(narr - 1): if array1d[i + 1] < array1d[i]: asc_arr[i + 1:] += period return asc_arr
[docs]def unjump_periodic(array1d, period): """ A given 1-d array's large gap (>period/2) is corrected. *Use case* If you have a series [179, -179, 179, -179] in a periodic system (-180 to 180, say longitude), it could be more natural to consider it as [179, 181, 179, 181] or [-181, -179, -181, -179]. *Sample* >>> arr = [179, -179, 179, -179] >>> arr2 = unjump_periodic(arr, 360) >>> print(arr2) [179 181 179 181] >>> arr = [179, -179, 179, np.nan, -179, 179, -179, np.nan, 179] >>> print(unjump_periodic(arr, 360)) [ 179. 181. 179. nan -179. -181. -179. nan 179.] >>> arr = [0, 180, 40] >>> print(unjump_periodic(arr, 360)) [ 0 -180 -320] """ arr = np.array(array1d, copy=True) diff = arr[1:] - arr[:-1] while np.nanmax(diff) >= period / 2: idx = np.where(diff >= period / 2)[0][-1] # The last one arr[idx+1:] = arr[idx+1:] - period diff = arr[1:] - arr[:-1] while np.nanmin(diff) < -period / 2: idx = np.where(diff < -period / 2)[0][-1] arr[idx+1:] = arr[idx+1:] + period diff = arr[1:] - arr[:-1] if np.nanmax(diff) >= period / 2 or np.nanmin(diff) < -period / 2: # I think this is not needed... return unjump_periodic(arr, period) return arr
[docs]def span(inarray): """ Just return the tuple of (min, max) for the given 1D array. :param inarray: np.ndarray :returns: A tuple, (min, max). Equivalent to (inarray.min(), inarray.max()) >>> arr = np.array([1, 2, 5, -3, 2, 8, 10, -2]) >>> print(span(arr)) (-3, 10) """ return inarray.min(), inarray.max()
[docs]def broadcast(array_in, shape_out): """ Create the broadcasted array with the specific shape from the given array .. note:: You can get the same results by .. code-block:: python broadcasted_array = array_in + np.zeros(shape_out) >>> array_in = np.arange(10).reshape([5, 2]) >>> arr2 = array_in + np.zeros([4, 3, 5, 2], dtype=int) # Equivalent to irfpy.util.nptools.broadcast >>> print(arr2.shape) (4, 3, 5, 2) >>> print(span(arr2[:, :, 0, 0])) (0, 0) >>> print(span(arr2[:, :, 3, 1])) (7, 7) >>> arr2 = array_in + np.zeros([5, 2, 3]) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): File "<stdin>", line 1, in ? ValueError: operands could not be broadcast together with shapes (5,2) (5,2,3) This looks 5 times faster than the broadcast method. .. todo:: Change the algorithm using the above algorithm The array is broadcasted to the given shape_out. :param array_in: Numpy array. Shape of (s0, s1, ...). :param shape_out: The shape_out. The shape_out must be (sX, sY, ..., s0, s1, ...). :returns: The broadcasted array. >>> array_in = np.arange(10).reshape([5, 2]) # (5, 2) shaped array is prepared. >>> arr2 = broadcast(array_in, (4, 3, 5, 2)) # Broadcast to (4, 3, 5, 2) array. >>> print(arr2.shape) (4, 3, 5, 2) >>> print(span(arr2[:, :, 0, 0])) (0, 0) >>> print(span(arr2[:, :, 3, 1])) (7, 7) >>> broadcast(array_in, (5, 2, 3)) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): File "<stdin>", line 1, in ? irfpy.util.exception.IrfpyException: 'Array cannot be broadcasted. Non-broadcastable. In:(5, 2) Out:(5, 2, 3)' >>> broadcast(array_in, (3,)) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): File "<stdin>", line 1, in ? irfpy.util.exception.IrfpyException: 'Array cannot be broadcasted. Lack of dimension. In:(5, 2) Out:(3,)' """ shape_in = array_in.shape dimen_in = array_in.ndim shape_out = tuple(shape_out) dimen_out = len(shape_out) ### Out array dimension is smaller than the in, raise an exception if dimen_in > dimen_out: raise ex.IrfpyException('Array cannot be broadcasted. Lack of dimension. In:{} Out:{}'. format(str(shape_in), str(shape_out))) ### Not broadcastable. if shape_out[-dimen_in:] != shape_in: raise ex.IrfpyException('Array cannot be broadcasted. Non-broadcastable. In:{} Out:{}'. format(str(shape_in), str(shape_out))) ### Tile function shape_tile = np.array(shape_out) shape_tile[-dimen_in:] = 1 # All the "broadcasting indexes" should be 1. array_out = np.tile(array_in, shape_tile) return array_out
[docs]def broadcast_r(array_in, shape_out): """ Broadcast, but in the revered direction .. note:: You can get the same results by .. code-block:: python broadcasted_array = array_in[..., np.newaxis, np.newaxis, np.newaxis] + np.zeros(shape_out) Here ``np.newaxis`` should repeat the proper times to make everythng consistent. >>> array_in = np.arange(10).reshape([5, 2]) # arrayn_in is (5, 2) shape >>> arr2 = array_in[..., np.newaxis, np.newaxis] + np.zeros([5, 2, 10, 3], dtype=int) # Equivalent to irfpy.util.nptools.broadcast_r >>> print(arr2.shape) (5, 2, 10, 3) >>> print(span(arr2[0, 0, :, :])) (0, 0) >>> print(span(arr2[3, 1, :, :])) (7, 7) This looks 10 times faster than the broadcast_r method. .. todo:: Change the algorithm using the above algorithm For example, an array with the shape of (5, 2) will be broadcasted to the array with the shape of (5, 2, 10, 4, 7). New_array[i, j, :, :, :] == Old_array[i, j]. :param array_in: np array with shape of (s0, s1, ...) :param shape_out: The aimed shape. It should be (s0, s1, ..., sX, sY, ...) :return: The broadcasted array from the "top" of the array. :See also: :class:`broadcast` >>> arrin = np.arange(10).reshape([5, 2]) # (5, 2) shaped array is prepared. >>> arrout = broadcast_r(arrin, (5, 2, 10, 3)) # (5, 2, 10, 3) array will be returned. >>> print(arrout.shape) (5, 2, 10, 3) >>> print(arrout[2, 1, :, :].max()) 5 >>> print(arrout[2, 1, :, :].min()) 5 >>> broadcast_r(arrin, (9, 5, 2)) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): File "<stdin>", line 1, in ? irfpy.util.exception.IrfpyException: 'Array cannot be broadcasted. Non-broadcastable. In:(5, 2) Out:(9, 5, 2)' """ shape_in = array_in.shape dimen_in = array_in.ndim shape_out = tuple(shape_out) dimen_out = len(shape_out) ### Out array dimension is smaller than the in, raise an exception if dimen_in > dimen_out: raise ex.IrfpyException('Array cannot be broadcasted. Lack of dimension. In:{} Out:{}'. format(str(shape_in), str(shape_out))) ### Check if broadcast works or not. if shape_out[:dimen_in] != shape_in: raise ex.IrfpyException('Array cannot be broadcasted. Non-broadcastable. In:{} Out:{}'. format(str(shape_in), str(shape_out))) ### Tile function. Shape to give np.tile is prepared. # shape_out is (s0, s1, ..., s(di-1), sX, sY, ...) shape_tile = np.roll(shape_out, -dimen_in) # (sX, sY, ..., s0, s1, ...) shape_tile[-dimen_in:] = 1 # (sX, sY, ..., 1, 1, ...) array_out = np.tile(array_in, shape_tile) # with shape of (sX, sY, ..., s0, s1, ...) ### Transpose. transpose_index = np.arange(dimen_out) # (0, 1, 2, ..., do-1) transpose_index = np.roll(transpose_index, dimen_in) # (..., do-1, 0, 1, 2, ...) array_out = np.transpose(array_out, axes=transpose_index) return array_out
[docs]def rtp(x, y, z): """ Cartesian x, y, and z is converted to r, |theta|, and |phi|. :param x: X :param y: Y :param z: Z :return: (r, |theta|, |phi|). |theta| and |phi| in radians. :rtype: A tuple I don't know why but the numpy or numpy.linalg package does not provide this. (Let me know if there is ofiicial version). It is very simple in math so can be easily implemented, but also good to have a function. >>> print(rtp(0, 0, 1)) (1.0, 0.0, 0.0) >>> r, t, p = rtp(1, 1, 1) >>> print('{:.3f} {:.3f} {:.3f}'.format(r, t, p)) 1.732 0.955 0.785 >>> r, t, p = rtp([1, 0, -1], [1, 1, 1], [0, 1, 1]) # List can be allowed. >>> print(r.shape) (3,) """ xx = np.array(x) yy = np.array(y) zz = np.array(z) r = np.sqrt(np.array(xx ** 2 + yy ** 2 + zz ** 2)) theta = np.arccos(np.clip(zz / r, -1, 1)) phi = np.arctan2(yy, xx) return (r, theta, phi)
[docs]def lonlat2xyz(longitude_deg, latitude_deg): """ Simply convert longitude and latitude to array :param longitude_deg: Longitude in degrees :param latitude_deg: Latitude in degrees :returns: Array with (..., 3) shape >>> lon = 0 >>> lat = 0 >>> print(lonlat2xyz(lon, lat)) [1. 0. 0.] >>> lon = [0, 0, 0, 90, 90, 90] >>> lat = [-90, 0, 90, -90, 0, 90] >>> xyz = lonlat2xyz(lon, lat) >>> print(xyz.shape) (6, 3) >>> print(trim_epsilon(xyz)) [[ 0. 0. -1.] [ 1. 0. 0.] [ 0. 0. 1.] [ 0. 0. -1.] [ 0. 1. 0.] [ 0. 0. 1.]] """ lonr = np.deg2rad(longitude_deg) latr = np.deg2rad(latitude_deg) x = np.cos(lonr) * np.cos(latr) y = np.sin(lonr) * np.cos(latr) z = np.sin(latr) vec = np.array([x, y, z]) vec = np.moveaxis(vec, 0, -1) return vec
[docs]def lingrids(x0, x1, nx): """ Return the grids. :param x0: x0 :param x1: x1 :param nx: Number of bins :return: A tuple, ``(grid_c, grid_b, grid_d)``. grid_c is the center of the grid with (nx) shape, and grid_b is the bound of the grid with (nx+1) shape. grid_d is the width of the bin with (nx) shape. - ``grid_b[0] = x0 and grid_b[nx] = x1``. - ``grid_c = (grid_b[1:] + grid_b[:-1]) / 2.`` - ``grid_b = np.linspace(x0, x1, nx+1)``. - ``grid_d = (grid_b[1:] - grid_b[:-1])``. >>> grid_c, grid_b, grid_d = lingrids(0, 10, 5) >>> print(grid_b[0]) 0.0 >>> print(grid_b[5]) 10.0 >>> print(grid_c[0]) 1.0 >>> print(grid_d[3]) # grid_d is always 2 in this case 2.0 """ grid_b = np.linspace(x0, x1, nx + 1) grid_c = (grid_b[1:] + grid_b[:-1]) / 2. grid_d = grid_b[1:] - grid_b[:-1] return grid_c, grid_b, grid_d
[docs]def loggrids(x0, x1, nx): """ Return the grid. :param x0: x0 :param x1: x1 :param nx: Number of bins :return: A tuple, ``(grid_c, grid_b, grid_d)``. grid_c is the center of the grid with (nx) shape, and grid_b is the bound of the grid with (nx+1) shape. grid_d is the width of the bin with (nx) shape. - ``grid_b[0] = 10^x0 and grid_b[nx] = 10^x1``. - ``grid_c = sqrt(grid_b[1:] * grid_b[:-1])`` - ``grid_b = np.linspace(x0, x1, nx+1)``. - ``grid_d = (grid_b[1:] - grid_b[:-1])``. >>> grid_c, grid_b, grid_d = loggrids(0, 10, 5) >>> print(grid_b[0]) 1.0 >>> print('{:.2e}'.format(grid_b[5])) 1.00e+10 >>> print('{:.2e}'.format(grid_c[0])) 1.00e+01 >>> print('{:.2e}'.format(grid_d[3])) # 10^8 - 10^6 9.90e+07 """ grid_b = np.logspace(x0, x1, nx + 1) grid_c = np.sqrt(grid_b[1:] * grid_b[:-1]) grid_d = grid_b[1:] - grid_b[:-1] return grid_c, grid_b, grid_d
[docs]def trim_epsilon(arr, epsilon=1e-15): r""" Trim the epsilon. :param arr: A numpy array :keyword epsilon: |epsilon| :return: The input ``arr`` but with substituting zero (+0.0) for the components where the absolute value is less than epsilon .. math:: b_i = \left\{ \begin{array}{ll} a_i & (|a_i| \ge \epsilon) \\ 0 & (|a_i| < \epsilon) \end{array}\right. """ return np.where(np.abs(arr) < epsilon, 0, arr)
[docs]def nearest_index(a, x): """ Return the index where the x is nearest to in array a. :param a: A list. It should be sorted in ascending order. :param x: Data :returns: The index of the data in the array ``a`` with the nearest data to ``x``. In other words, the index ``i`` usually satisfies ``a[i] < x <= a[i+1]``. >>> xlist = [100, 200, 400, 800, 1600] >>> print(nearest_index(xlist, 50)) 0 >>> print(nearest_index(xlist, 100)) 0 >>> print(nearest_index(xlist, 149)) 0 >>> print(nearest_index(xlist, 150)) 1 >>> print(nearest_index(xlist, 151)) 1 >>> print(nearest_index(xlist, 200)) 1 >>> print(nearest_index(xlist, 300)) 2 >>> print(nearest_index(xlist, 400)) 2 >>> print(nearest_index(xlist, 1600)) 4 >>> print(nearest_index(xlist, 1800)) 4 """ return nearest_index_right(a, x)
[docs]def nearest_index_right(xlist, x0): """ Same as :func:`nearest_index` """ idxbi = _bisect.bisect(xlist, x0) if idxbi == 0: return idxbi if idxbi == len(xlist): return idxbi - 1 xm = x0 - xlist[idxbi - 1] xp = xlist[idxbi] - x0 if xm < xp: return idxbi - 1 else: return idxbi
[docs]def nearest_index_left(xlist, x0): """ Return the index where the x is nearest to in array a. :param a: A list. It should be sorted in ascending order. :param x: Data :returns: The index of the data in the array ``a`` with the nearest data to ``x``. In other words, the index ``i`` usually satisfies ``a[i] < x <= a[i+1]``. >>> xlist = [100, 200, 400, 800, 1600] >>> print(nearest_index_left(xlist, 50)) 0 >>> print(nearest_index_left(xlist, 100)) 0 >>> print(nearest_index_left(xlist, 149)) 0 >>> print(nearest_index_left(xlist, 150)) 0 >>> print(nearest_index_left(xlist, 151)) 1 >>> print(nearest_index_left(xlist, 200)) 1 >>> print(nearest_index_left(xlist, 300)) 1 >>> print(nearest_index_left(xlist, 400)) 2 >>> print(nearest_index_left(xlist, 1600)) 4 >>> print(nearest_index_left(xlist, 1800)) 4 """ idxbi = _bisect.bisect(xlist, x0) if idxbi == 0: return idxbi if idxbi == len(xlist): return idxbi - 1 xm = x0 - xlist[idxbi - 1] xp = xlist[idxbi] - x0 if xm <= xp: return idxbi - 1 else: return idxbi