diff --git a/meteoinfo-lab/milconfig.xml b/meteoinfo-lab/milconfig.xml index 42bdf7f0..063f604e 100644 --- a/meteoinfo-lab/milconfig.xml +++ b/meteoinfo-lab/milconfig.xml @@ -1,14 +1,10 @@ - - - - + - @@ -16,19 +12,23 @@ + + + + - + - + diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/_ndarray.py b/meteoinfo-lab/pylib/mipylib/numeric/core/_ndarray.py index b7a38fa0..63e9671a 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/core/_ndarray.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/core/_ndarray.py @@ -20,8 +20,8 @@ class NDArray(object): if not isinstance(array, Array): array = ArrayUtil.array(array, None) - if array.getRank() == 0 and array.getDataType() != DataType.STRUCTURE: - array = ArrayUtil.array([array.getIndexIterator().getObjectNext()]) + # if array.getRank() == 0 and array.getDataType() != DataType.STRUCTURE: + # array = ArrayUtil.array([array.getIndexIterator().getObjectNext()]) self._array = array self.ndim = array.getRank() @@ -153,7 +153,7 @@ class NDArray(object): indices = nindices if len(indices) != self.ndim: - print 'indices must be ' + str(self.ndim) + ' dimensions!' + print('indices must be ' + str(self.ndim) + ' dimensions!') raise IndexError() ranges = [] @@ -582,6 +582,24 @@ class NDArray(object): r = NDArray(ArrayMath.inValues(self._array, other)) return r + def all(self, axis=None): + """ + Test whether all array element along a given axis evaluates to True. + + :param x: (*array_like or list*) Input values. + :param axis: (*int*) Axis along which a logical OR reduction is performed. + The default (axis = None) is to perform a logical OR over all the + dimensions of the input array. + + :returns: (*array_like*) All result + """ + if axis is None: + return ArrayMath.all(self._array) + else: + if axis < 0: + axis += self.ndim + return NDArray(ArrayMath.all(self._array, axis)) + def any(self, axis=None): """ Test whether any array element along a given axis evaluates to True. @@ -942,7 +960,7 @@ class NDArray(object): :param axis: (*int*) Axis along which the standard deviation is computed. The default is to compute the standard deviation of the flattened array. :param ddof: (*int*) Delta Degrees of Freedom: the divisor used in the calculation is - N - ddof, where N represents the number of elements. By default ddof is zero. + N - ddof, where N represents the number of elements. By default, ddof is zero. returns: (*array_like*) Standart deviation result. """ @@ -1222,7 +1240,7 @@ class NDArray(object): :param indices: (*array_like*) The indices of the values to extract. :param axis: (*int*) The axis over which to select values. - :returns: (*array*) The returned array has the same type as a. + :returns: (*array*) The returned array has the same type as this array. """ if isinstance(indices, (list, tuple)): indices = NDArray(indices) diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/fromnumeric$py.class b/meteoinfo-lab/pylib/mipylib/numeric/core/fromnumeric$py.class index 842e057d..401055b6 100644 Binary files a/meteoinfo-lab/pylib/mipylib/numeric/core/fromnumeric$py.class and b/meteoinfo-lab/pylib/mipylib/numeric/core/fromnumeric$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/fromnumeric.py b/meteoinfo-lab/pylib/mipylib/numeric/core/fromnumeric.py index 3f5dc02f..3af84f48 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/core/fromnumeric.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/core/fromnumeric.py @@ -1,6 +1,6 @@ # coding=utf-8 -from .numeric import asarray, array +from .numeric import asarray, array, isscalar from ._ndarray import NDArray from .dimarray import DimArray from org.meteoinfo.ndarray.math import ArrayUtil, ArrayMath @@ -39,6 +39,9 @@ def ndim(a): >>> np.ndim(1) 0 """ + if isscalar(a): + return 0 + try: return a.ndim except AttributeError: diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class index 8ddb3749..5c4d87ea 100644 Binary files a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class and b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py index b64ffa4e..cfbd8bc8 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py @@ -110,14 +110,15 @@ def array(object, dtype=None, copy=True, order='K', subok=False, ndmin=0): elif miutil.iscomplex(object): a = NDArray(JythonUtil.toComplexArray(object)) - if isinstance(dtype, basestring): - dtype = _dtype.DataType(dtype) - - if not dtype is None: - dtype = dtype._dtype - if a is None: - a = NDArray(ArrayUtil.array(object, dtype)) + if isinstance(dtype, basestring): + dtype = _dtype.DataType(dtype) + + if not dtype is None: + dtype = dtype._dtype + a = NDArray(ArrayUtil.array(object, dtype)) + else: + a = NDArray(object) if a.ndim < ndmin: shape = [] @@ -1233,12 +1234,7 @@ def all(x, axis=None): if isinstance(x, list): x = array(x) - if axis is None: - return ArrayMath.all(x._array) - else: - if axis < 0: - axis += x.ndim - return NDArray(ArrayMath.all(x._array, axis)) + return x.all(axis) def isclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False): """ diff --git a/meteoinfo-lab/pylib/mipylib/numeric/lib/function_base.py b/meteoinfo-lab/pylib/mipylib/numeric/lib/function_base.py index f48b9383..e8a17e20 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/lib/function_base.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/lib/function_base.py @@ -1,10 +1,11 @@ +from mipylib import numeric as np from ..core import numeric as _nx from ..core import dtype from ..core._ndarray import NDArray from ..core.fromnumeric import (ravel, nonzero) from org.meteoinfo.ndarray.math import ArrayMath -__all__ = ['angle','extract', 'place', 'grid_edge'] +__all__ = ['angle','extract', 'place', 'grid_edge', 'gradient'] def extract(condition, arr): @@ -161,3 +162,207 @@ def angle(z, deg=False): if deg: a *= 180 / _nx.pi return a + +def gradient(f, *varargs, **kwargs): + """ + Return the gradient of an N-dimensional array. + + The gradient is computed using second order accurate central differences + in the interior points and either first or second order accurate one-sides + (forward or backwards) differences at the boundaries. + The returned gradient hence has the same shape as the input array. + + Parameters + ---------- + f : array_like + An N-dimensional array containing samples of a scalar function. + varargs : list of scalar or array, optional + Spacing between f values. Default unitary spacing for all dimensions. + Spacing can be specified using: + + 1. single scalar to specify a sample distance for all dimensions. + 2. N scalars to specify a constant sample distance for each dimension. + i.e. `dx`, `dy`, `dz`, ... + 3. N arrays to specify the coordinates of the values along each + dimension of F. The length of the array must match the size of + the corresponding dimension + 4. Any combination of N scalars/arrays with the meaning of 2. and 3. + + If `axis` is given, the number of varargs must equal the number of axes. + Default: 1. + + edge_order : {1, 2}, optional + Gradient is calculated using N-th order accurate differences + at the boundaries. Default: 1. + + axis : None or int or tuple of ints, optional + Gradient is calculated only along the given axis or axes + The default (axis = None) is to calculate the gradient for all the axes + of the input array. axis may be negative, in which case it counts from + the last to the first axis. + + Returns + ------- + gradient : ndarray or list of ndarray + A list of ndarrays (or a single ndarray if there is only one dimension) + corresponding to the derivatives of f with respect to each dimension. + Each derivative has the same shape as f. + """ + f = _nx.asanyarray(f) + N = f.ndim # number of dimensions + + axis = kwargs.pop('axis', None) + edge_order = kwargs.pop('edge_order', 1) + + if axis is None: + axes = tuple(range(N)) + else: + axes = _nx.normalize_axis_tuple(axis, N) + + len_axes = len(axes) + n = len(varargs) + if n == 0: + # no spacing argument - use 1 in all axes + dx = [1.0] * len_axes + elif n == 1 and np.ndim(varargs[0]) == 0: + # single scalar for all axes + dx = varargs * len_axes + elif n == len_axes: + # scalar or 1d array for each axis + dx = list(varargs) + for i, distances in enumerate(dx): + distances = _nx.asanyarray(distances) + if distances.ndim == 0: + continue + elif distances.ndim != 1: + raise ValueError("distances must be either scalars or 1d") + if len(distances) != f.shape[axes[i]]: + raise ValueError("when 1d, distances must match " + "the length of the corresponding dimension") + if distances.dtype == dtype.int: + # Convert numpy integer types to float to avoid modular + # arithmetic in np.diff(distances). + distances = distances.astype(dtype.float) + diffx = _nx.diff(distances) + # if distances are constant reduce to the scalar case + # since it brings a consistent speedup + if (diffx == diffx[0]).all(): + diffx = diffx[0] + dx[i] = diffx + else: + raise TypeError("invalid number of arguments") + + if edge_order > 2: + raise ValueError("'edge_order' greater than 2 not supported") + + # use central differences on interior and one-sided differences on the + # endpoints. This preserves second order-accuracy over the full domain. + + outvals = [] + + # create slice objects --- initially all are [:, :, ..., :] + slice1 = [slice(None)]*N + slice2 = [slice(None)]*N + slice3 = [slice(None)]*N + slice4 = [slice(None)]*N + + otype = f.dtype + if otype == dtype.int: + otype = dtype.float + + for axis, ax_dx in zip(axes, dx): + if f.shape[axis] < edge_order + 1: + raise ValueError( + "Shape of array too small to calculate a numerical gradient, " + "at least (edge_order + 1) elements are required.") + # result allocation + out = np.empty_like(f, dtype=otype) + + # spacing for the current axis + uniform_spacing = np.ndim(ax_dx) == 0 + + # Numerical differentiation: 2nd order interior + slice1[axis] = slice(1, -1) + slice2[axis] = slice(None, -2) + slice3[axis] = slice(1, -1) + slice4[axis] = slice(2, None) + + if uniform_spacing: + out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx) + else: + dx1 = ax_dx[0:-1] + dx2 = ax_dx[1:] + a = -(dx2)/(dx1 * (dx1 + dx2)) + b = (dx2 - dx1) / (dx1 * dx2) + c = dx1 / (dx2 * (dx1 + dx2)) + # fix the shape for broadcasting + shape = np.ones(N, dtype=dtype.int) + shape[axis] = -1 + a.shape = b.shape = c.shape = shape + # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:] + out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)] + + # Numerical differentiation: 1st order edges + if edge_order == 1: + slice1[axis] = 0 + slice2[axis] = 1 + slice3[axis] = 0 + dx_0 = ax_dx if uniform_spacing else ax_dx[0] + # 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0]) + out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0 + + slice1[axis] = -1 + slice2[axis] = -1 + slice3[axis] = -2 + dx_n = ax_dx if uniform_spacing else ax_dx[-1] + # 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2]) + out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n + + # Numerical differentiation: 2nd order edges + else: + slice1[axis] = 0 + slice2[axis] = 0 + slice3[axis] = 1 + slice4[axis] = 2 + if uniform_spacing: + a = -1.5 / ax_dx + b = 2. / ax_dx + c = -0.5 / ax_dx + else: + dx1 = ax_dx[0] + dx2 = ax_dx[1] + a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2)) + b = (dx1 + dx2) / (dx1 * dx2) + c = - dx1 / (dx2 * (dx1 + dx2)) + # 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2] + out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)] + + slice1[axis] = -1 + slice2[axis] = -3 + slice3[axis] = -2 + slice4[axis] = -1 + if uniform_spacing: + a = 0.5 / ax_dx + b = -2. / ax_dx + c = 1.5 / ax_dx + else: + dx1 = ax_dx[-2] + dx2 = ax_dx[-1] + a = (dx2) / (dx1 * (dx1 + dx2)) + b = - (dx2 + dx1) / (dx1 * dx2) + c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2)) + # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1] + out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)] + + outvals.append(out) + + # reset the slice object in this dimension to ":" + slice1[axis] = slice(None) + slice2[axis] = slice(None) + slice3[axis] = slice(None) + slice4[axis] = slice(None) + + if len_axes == 1: + return outvals[0] + else: + return outvals \ No newline at end of file diff --git a/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayMath.java b/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayMath.java index 206aaeae..eaceb98d 100644 --- a/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayMath.java +++ b/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayMath.java @@ -6022,8 +6022,10 @@ public class ArrayMath { public static Array setSection(Array a, List ranges, Array v) throws InvalidRangeException { Array r = a.section(ranges); IndexIterator iter = r.getIndexIterator(); - if (r.getShape() != v.getShape()) { - v = ArrayUtil.broadcast(v, r.getShape()); + if (r.getSize() != v.getSize()) { + if (r.getShape() != v.getShape()) { + v = ArrayUtil.broadcast(v, r.getShape()); + } } Index index = v.getIndex(); while (iter.hasNext()) { diff --git a/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayUtil.java b/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayUtil.java index ce1c2591..91541f41 100644 --- a/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayUtil.java +++ b/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayUtil.java @@ -652,7 +652,7 @@ public class ArrayUtil { } else { if (dt == null) dt = ArrayMath.getDataType(data); - Array a = Array.factory(dt, new int[]{1}); + Array a = Array.factory(dt, new int[0]); if (data instanceof String) { a.setString(0, (String)data); } else {