# 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
sort_periodic
unjump_periodic
span
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.')

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

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 == 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)

""" 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
>>> 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.
"""
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 = arr + 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

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)[-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)[-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()

""" Create the broadcasted array with the specific shape from the given array

.. note::

You can get the same results by

.. code-block:: python

>>> 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, ...).

>>> 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 ?

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

if shape_out[-dimen_in:] != shape_in:
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

""" 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 ?

"""
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:
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.]]
"""

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 = 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
>>> print(grid_b)
10.0
>>> print(grid_c)
1.0
>>> print(grid_d)   # 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 = 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)
1.0
>>> print('{:.2e}'.format(grid_b))
1.00e+10
>>> print('{:.2e}'.format(grid_c))
1.00e+01
>>> print('{:.2e}'.format(grid_d))   # 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