add convolve function

This commit is contained in:
wyq 2024-04-14 09:08:16 +08:00
parent 93b1f6b89f
commit b3b0e24d4f
17 changed files with 955 additions and 31 deletions

View File

@ -1,12 +1,6 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<MeteoInfo File="milconfig.xml" Type="configurefile">
<Path OpenPath="D:\Working\MIScript\Jython\mis\plot_types\wind">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\bar"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\contour"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\imshow"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\patch"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\pcolor"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\pie"/>
<Path OpenPath="D:\Working\MIScript\Jython\mis\array">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\plot"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\polar"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\scatter"/>
@ -14,19 +8,25 @@
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\taylor"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\violinplot"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\weather"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\wind"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\fitting"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\linalg"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
</Path>
<File>
<OpenedFiles>
<OpenedFile File="D:\Working\MIScript\Jython\mis\map\geoshow\geoshow_2.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\io\geotiff\geotiff_dem_4.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\wind\windrose_data_func.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\polynomial\polyder_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\array\convolve_1.py"/>
</OpenedFiles>
<RecentFiles>
<RecentFile File="D:\Working\MIScript\Jython\mis\map\geoshow\geoshow_2.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\io\geotiff\geotiff_dem_4.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\wind\windrose_data_func.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\polynomial\polyder_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\array\convolve_1.py"/>
</RecentFiles>
</File>
<Font>
@ -34,5 +34,5 @@
</Font>
<LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
<Figure DoubleBuffering="true"/>
<Startup MainFormLocation="-7,0" MainFormSize="1376,799"/>
<Startup MainFormLocation="-7,0" MainFormSize="1400,825"/>
</MeteoInfo>

View File

@ -15,7 +15,7 @@ import mipylib.numeric.optimize as so
__all__ = [
'equivalent_potential_temperature','exner_function',
'mixing_ratio','mixing_ratio_from_specific_humidity','potential_temperature',
'relative_humidity_from_specific_humidity',
'relative_humidity_from_specific_humidity', 'relative_humidity_from_dewpoint',
'saturation_mixing_ratio','saturation_vapor_pressure','temperature_from_potential_temperature',
'virtual_temperature','dry_static_energy','isentropic_interpolation',
'dewpoint','dewpoint_from_relative_humidity','specific_humidity_from_dewpoint',
@ -162,6 +162,43 @@ def relative_humidity_from_specific_humidity(pressure, temperature, specific_hum
return (mixing_ratio_from_specific_humidity(specific_humidity)
/ saturation_mixing_ratio(pressure, temperature))
def relative_humidity_from_dewpoint(temperature, dewpoint):
r"""Calculate the relative humidity.
Uses temperature and dewpoint to calculate relative humidity as the ratio of vapor
pressure to saturation vapor pressures.
Parameters
----------
temperature : `float`
Air temperature (celsius)
dewpoint : `float`
Dewpoint temperature (celsius)
Returns
-------
`float`
Relative humidity
Examples
--------
>>> relative_humidity_from_dewpoint(25, 12)
<44.2484765>
See Also
--------
saturation_vapor_pressure
Notes
-----
.. math:: RH = \frac{e(T_d)}{e_s(T)}
"""
e = saturation_vapor_pressure(dewpoint + 273.15)
e_s = saturation_vapor_pressure(temperature + 273.15)
return e / e_s
def exner_function(pressure, reference_pressure=constants.P0):
r"""Calculate the Exner function.
@ -208,7 +245,7 @@ def potential_temperature(pressure, temperature):
Returns
-------
array_like
The potential temperature corresponding to the the temperature and
The potential temperature corresponding to the temperature and
pressure.
See Also

View File

@ -214,12 +214,10 @@ def dewpoint2rh(dewpoint, temp):
Relative humidity - %
"""
if isinstance(dewpoint, NDArray):
r = NDArray(MeteoMath.dewpoint2rh(dewpoint.asarray(), temp.asarray()))
if isinstance(dewpoint, DimArray):
r = DimArray(r, dewpoint.dims, dewpoint.fill_value, dewpoint.proj)
return r
r = MeteoMath.dewpoint2rh(dewpoint.asarray(), temp.asarray())
return dewpoint.array_wrap(r)
else:
return MeteoMath.dewpoint2rh(temp, dewpoint)
return MeteoMath.dewpoint2rh(dewpoint, temp)
def rh2dewpoint(rh, temp):
"""

View File

@ -40,7 +40,7 @@ __all__ = [
'argmin','argmax','argsort','array','array_split','amax','amin','asanyarray','asarray',
'arcsin','asin','asmiarray','atleast_1d','atleast_2d','arctan','atan',
'arctan2','atan2','average','histogram','broadcast_to','cdiff','ceil',
'concatenate','conj','conjugate','corrcoef','cos','cosh','cylinder','degrees','delnan','diag','diff',
'concatenate','conj','conjugate','convolve','corrcoef','cos','cosh','cylinder','degrees','delnan','diag','diff',
'datatable','dot','empty','empty_like','exp','eye','flatnonzero','floor',
'fmax','fmin','full','hcurl','hdivg','hstack','hypot','identity','indices','interp2d','interpn','isarray',
'isclose','isfinite','isinf','isnan','isscalar','linspace','log','log10','logical_not','logspace',
@ -2133,6 +2133,102 @@ def hstack(tup):
else:
return concatenate(arrs, 1)
def convolve(a, v, mode='full'):
"""
Returns the discrete, linear convolution of two one-dimensional sequences.
The convolution operator is often seen in signal processing, where it
models the effect of a linear time-invariant system on a signal [1]_. In
probability theory, the sum of two independent random variables is
distributed according to the convolution of their individual
distributions.
If `v` is longer than `a`, the arrays are swapped before computation.
Parameters
----------
a : (N,) array_like
First one-dimensional input array.
v : (M,) array_like
Second one-dimensional input array.
mode : {'full', 'valid', 'same'}, optional
'full':
By default, mode is 'full'. This returns the convolution
at each point of overlap, with an output shape of (N+M-1,). At
the end-points of the convolution, the signals do not overlap
completely, and boundary effects may be seen.
'same':
Mode 'same' returns output of length ``max(M, N)``. Boundary
effects are still visible.
'valid':
Mode 'valid' returns output of length
``max(M, N) - min(M, N) + 1``. The convolution product is only given
for points where the signals overlap completely. Values outside
the signal boundary have no effect.
Returns
-------
out : ndarray
Discrete, linear convolution of `a` and `v`.
Notes
-----
The discrete convolution operation is defined as
.. math:: (a * v)_n = \\sum_{m = -\\infty}^{\\infty} a_m v_{n - m}
It can be shown that a convolution :math:`x(t) * y(t)` in time/space
is equivalent to the multiplication :math:`X(f) Y(f)` in the Fourier
domain, after appropriate padding (padding is necessary to prevent
circular convolution). Since multiplication is more efficient (faster)
than convolution, the function `scipy.signal.fftconvolve` exploits the
FFT to calculate the convolution of large data-sets.
References
----------
.. [1] Wikipedia, "Convolution",
https://en.wikipedia.org/wiki/Convolution
Examples
--------
Note how the convolution operator flips the second array
before "sliding" the two across one another:
>>> np.convolve([1, 2, 3], [0, 1, 0.5])
array([0. , 1. , 2.5, 4. , 1.5])
Only return the middle values of the convolution.
Contains boundary effects, where zeros are taken
into account:
>>> np.convolve([1,2,3],[0,1,0.5], 'same')
array([1. , 2.5, 4. ])
The two arrays are of the same length, so there
is only one position where they completely overlap:
>>> np.convolve([1,2,3],[0,1,0.5], 'valid')
array([2.5])
"""
a, v = asarray(a), asarray(v)
if (len(v) > len(a)):
a, v = v, a
if len(a) == 0:
raise ValueError('a cannot be empty')
if len(v) == 0:
raise ValueError('v cannot be empty')
if mode == 'full':
r = ArrayMath.convolveFull(a._array, v._array)
elif mode == 'same':
r = ArrayMath.convolveSame(a._array, v._array)
else:
r = ArrayMath.convolveValid(a._array, v._array)
return NDArray(r)
def dot(a, b):
"""
Matrix multiplication.

View File

@ -57,7 +57,7 @@ def expfit(x, y, func=False):
def polyfit(x, y, degree, func=False):
"""
Polynomail fitting.
Polynomial fitting.
:param x: (*array_like*) x data array.
:param y: (*array_like*) y data array.
@ -88,7 +88,7 @@ def polyval(p, x):
composite polynomial p(x(t)) is returned.
:param p: (*array_like*) 1D array of polynomial coefficients (including coefficients equal to zero)
from highest degree to the constant term.
from the highest degree to the constant term.
:param x: (*array_like*) A number, an array of numbers, or an instance of poly1d, at which to evaluate
p.

View File

@ -11,6 +11,7 @@ from .type_check import *
from .arraysetops import *
from .npyio import *
from .matio import *
from .polynomial import *
__all__ = []
__all__ += shape_base.__all__
@ -20,4 +21,5 @@ __all__ += stride_tricks.__all__
__all__ += type_check.__all__
__all__ += arraysetops.__all__
__all__ += npyio.__all__
__all__ += matio.__all__
__all__ += matio.__all__
__all__ += polynomial.__all__

View File

@ -0,0 +1,616 @@
"""
Functions to operate on polynomials.
"""
from org.meteoinfo.math.fitting import FittingUtil
from org.meteoinfo.ndarray.math import ArrayMath, ArrayUtil
from ..core import numeric as NX
from ..core import NDArray, atleast_1d
from ..linalg import eigvals
import warnings
__all__ = ['polyfit','polyval','polyder']
def poly(seq_of_zeros):
"""
Find the coefficients of a polynomial with the given sequence of roots.
Returns the coefficients of the polynomial whose leading coefficient
is one for the given sequence of zeros (multiple roots must be included
in the sequence as many times as their multiplicity; see Examples).
A square matrix (or array, which will be treated as a matrix) can also
be given, in which case the coefficients of the characteristic polynomial
of the matrix are returned.
Parameters
----------
seq_of_zeros : array_like, shape (N,) or (N, N)
A sequence of polynomial roots, or a square array or matrix object.
Returns
-------
c : ndarray
1D array of polynomial coefficients from highest to lowest degree:
``c[0] * x**(N) + c[1] * x**(N-1) + ... + c[N-1] * x + c[N]``
where c[0] always equals 1.
Raises
------
ValueError
If input is the wrong shape (the input must be a 1-D or square
2-D array).
See Also
--------
polyval : Compute polynomial values.
roots : Return the roots of a polynomial.
polyfit : Least squares polynomial fit.
poly1d : A one-dimensional polynomial class.
Notes
-----
Specifying the roots of a polynomial still leaves one degree of
freedom, typically represented by an undetermined leading
coefficient. [1]_ In the case of this function, that coefficient -
the first one in the returned array - is always taken as one. (If
for some reason you have one other point, the only automatic way
presently to leverage that information is to use ``polyfit``.)
The characteristic polynomial, :math:`p_a(t)`, of an `n`-by-`n`
matrix **A** is given by
:math:`p_a(t) = \\mathrm{det}(t\\, \\mathbf{I} - \\mathbf{A})`,
where **I** is the `n`-by-`n` identity matrix. [2]_
References
----------
.. [1] M. Sullivan and M. Sullivan, III, "Algebra and Trigonometry,
Enhanced With Graphing Utilities," Prentice-Hall, pg. 318, 1996.
.. [2] G. Strang, "Linear Algebra and Its Applications, 2nd Edition,"
Academic Press, pg. 182, 1980.
Examples
--------
Given a sequence of a polynomial's zeros:
>>> np.poly((0, 0, 0)) # Multiple root example
array([1., 0., 0., 0.])
The line above represents z**3 + 0*z**2 + 0*z + 0.
>>> np.poly((-1./2, 0, 1./2))
array([ 1. , 0. , -0.25, 0. ])
The line above represents z**3 - z/4
>>> np.poly((np.random.random(1)[0], 0, np.random.random(1)[0]))
array([ 1. , -0.77086955, 0.08618131, 0. ]) # random
Given a square array object:
>>> P = np.array([[0, 1./3], [-1./2, 0]])
>>> np.poly(P)
array([1. , 0. , 0.16666667])
Note how in all cases the leading coefficient is always 1.
"""
seq_of_zeros = atleast_1d(seq_of_zeros)
sh = seq_of_zeros.shape
if len(sh) == 2 and sh[0] == sh[1] and sh[0] != 0:
seq_of_zeros = eigvals(seq_of_zeros)
elif len(sh) == 1:
dt = seq_of_zeros.dtype
# Let object arrays slip through, e.g. for arbitrary precision
if dt != object:
seq_of_zeros = seq_of_zeros.astype(mintypecode(dt.char))
else:
raise ValueError("input must be 1d or non-empty square 2d array.")
if len(seq_of_zeros) == 0:
return 1.0
dt = seq_of_zeros.dtype
a = ones((1,), dtype=dt)
for zero in seq_of_zeros:
a = NX.convolve(a, array([1, -zero], dtype=dt), mode='full')
if issubclass(a.dtype.type, NX.complexfloating):
# if complex roots are all complex conjugates, the roots are real.
roots = NX.asarray(seq_of_zeros, complex)
if NX.all(NX.sort(roots) == NX.sort(roots.conjugate())):
a = a.real.copy()
return a
def polyfit(x, y, degree, func=False):
"""
Polynomial fitting.
:param x: (*array_like*) x data array.
:param y: (*array_like*) y data array.
:param degree: (*int*) Degree of the fitting polynomial.
:param func: (*boolean*) Return fit function (for predict function) or not. Default is ``False``.
:returns: Fitting parameters and function (optional).
"""
if isinstance(x, list):
x = NDArray(ArrayUtil.array(x))
if isinstance(y, list):
y = NDArray(ArrayUtil.array(y))
r = FittingUtil.polyFit(x.asarray(), y.asarray(), degree)
if func:
return r[0], r[1], r[2]
else:
return r[0], r[1]
def polyval(p, x):
"""
Evaluate a polynomial at specific values.
If p is of length N, this function returns the value:
p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]
If x is a sequence, then p(x) is returned for each element of x. If x is another polynomial then the
composite polynomial p(x(t)) is returned.
:param p: (*array_like*) 1D array of polynomial coefficients (including coefficients equal to zero)
from the highest degree to the constant term.
:param x: (*array_like*) A number, an array of numbers, or an instance of poly1d, at which to evaluate
p.
:returns: Polynomial value
"""
if isinstance(x, list):
x = NDArray(ArrayUtil.array(x))
return NDArray(ArrayMath.polyVal(p, x.asarray()))
def polyder(p, m=1):
"""
Return the derivative of the specified order of a polynomial.
.. note::
This forms part of the old polynomial API. Since version 1.4, the
new polynomial API defined in `numpy.polynomial` is preferred.
A summary of the differences can be found in the
:doc:`transition guide </reference/routines.polynomials>`.
Parameters
----------
p : poly1d or sequence
Polynomial to differentiate.
A sequence is interpreted as polynomial coefficients, see `poly1d`.
m : int, optional
Order of differentiation (default: 1)
Returns
-------
der : poly1d
A new polynomial representing the derivative.
See Also
--------
polyint : Anti-derivative of a polynomial.
poly1d : Class for one-dimensional polynomials.
Examples
--------
The derivative of the polynomial :math:`x^3 + x^2 + x^1 + 1` is:
>>> p = np.poly1d([1,1,1,1])
>>> p2 = np.polyder(p)
>>> p2
poly1d([3, 2, 1])
which evaluates to:
>>> p2(2.)
17.0
We can verify this, approximating the derivative with
``(f(x + h) - f(x))/h``:
>>> (p(2. + 0.001) - p(2.)) / 0.001
17.007000999997857
The fourth-order derivative of a 3rd-order polynomial is zero:
>>> np.polyder(p, 2)
poly1d([6, 2])
>>> np.polyder(p, 3)
poly1d([6])
>>> np.polyder(p, 4)
poly1d([0])
"""
m = int(m)
if m < 0:
raise ValueError("Order of derivative must be positive (see polyint)")
#truepoly = isinstance(p, poly1d)
p = NX.asarray(p)
n = len(p) - 1
y = p[:-1] * NX.arange(n, 0, -1)
if m == 0:
val = p
else:
val = polyder(y, m - 1)
#if truepoly:
# val = poly1d(val)
return val
class poly1d:
"""
A one-dimensional polynomial class.
A convenience class, used to encapsulate "natural" operations on
polynomials so that said operations may take on their customary
form in code (see Examples).
Parameters
----------
c_or_r : array_like
The polynomial's coefficients, in decreasing powers, or if
the value of the second parameter is True, the polynomial's
roots (values where the polynomial evaluates to 0). For example,
``poly1d([1, 2, 3])`` returns an object that represents
:math:`x^2 + 2x + 3`, whereas ``poly1d([1, 2, 3], True)`` returns
one that represents :math:`(x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x -6`.
r : bool, optional
If True, `c_or_r` specifies the polynomial's roots; the default
is False.
variable : str, optional
Changes the variable used when printing `p` from `x` to `variable`
(see Examples).
Examples
--------
Construct the polynomial :math:`x^2 + 2x + 3`:
>>> p = np.poly1d([1, 2, 3])
>>> print(np.poly1d(p))
2
1 x + 2 x + 3
Evaluate the polynomial at :math:`x = 0.5`:
>>> p(0.5)
4.25
Find the roots:
>>> p.r
array([-1.+1.41421356j, -1.-1.41421356j])
>>> p(p.r)
array([ -4.44089210e-16+0.j, -4.44089210e-16+0.j]) # may vary
These numbers in the previous line represent (0, 0) to machine precision
Show the coefficients:
>>> p.c
array([1, 2, 3])
Display the order (the leading zero-coefficients are removed):
>>> p.order
2
Show the coefficient of the k-th power in the polynomial
(which is equivalent to ``p.c[-(i+1)]``):
>>> p[1]
2
Polynomials can be added, subtracted, multiplied, and divided
(returns quotient and remainder):
>>> p * p
poly1d([ 1, 4, 10, 12, 9])
>>> (p**3 + 4) / p
(poly1d([ 1., 4., 10., 12., 9.]), poly1d([4.]))
``asarray(p)`` gives the coefficient array, so polynomials can be
used in all functions that accept arrays:
>>> p**2 # square of polynomial
poly1d([ 1, 4, 10, 12, 9])
>>> np.square(p) # square of individual coefficients
array([1, 4, 9])
The variable used in the string representation of `p` can be modified,
using the `variable` parameter:
>>> p = np.poly1d([1,2,3], variable='z')
>>> print(p)
2
1 z + 2 z + 3
Construct a polynomial from its roots:
>>> np.poly1d([1, 2], True)
poly1d([ 1., -3., 2.])
This is the same polynomial as obtained by:
>>> np.poly1d([1, -1]) * np.poly1d([1, -2])
poly1d([ 1, -3, 2])
"""
__hash__ = None
@property
def coeffs(self):
""" The polynomial coefficients """
return self._coeffs
@coeffs.setter
def coeffs(self, value):
# allowing this makes p.coeffs *= 2 legal
if value is not self._coeffs:
raise AttributeError("Cannot set attribute")
@property
def variable(self):
""" The name of the polynomial variable """
return self._variable
# calculated attributes
@property
def order(self):
""" The order or degree of the polynomial """
return len(self._coeffs) - 1
@property
def roots(self):
""" The roots of the polynomial, where self(x) == 0 """
return roots(self._coeffs)
# our internal _coeffs property need to be backed by __dict__['coeffs'] for
# scipy to work correctly.
@property
def _coeffs(self):
return self.__dict__['coeffs']
@_coeffs.setter
def _coeffs(self, coeffs):
self.__dict__['coeffs'] = coeffs
# alias attributes
r = roots
c = coef = coefficients = coeffs
o = order
def __init__(self, c_or_r, r=False, variable=None):
if isinstance(c_or_r, poly1d):
self._variable = c_or_r._variable
self._coeffs = c_or_r._coeffs
if set(c_or_r.__dict__) - set(self.__dict__):
msg = ("In the future extra properties will not be copied "
"across when constructing one poly1d from another")
warnings.warn(msg, FutureWarning, stacklevel=2)
self.__dict__.update(c_or_r.__dict__)
if variable is not None:
self._variable = variable
return
if r:
c_or_r = poly(c_or_r)
c_or_r = atleast_1d(c_or_r)
if c_or_r.ndim > 1:
raise ValueError("Polynomial must be 1d only.")
c_or_r = trim_zeros(c_or_r, trim='f')
if len(c_or_r) == 0:
c_or_r = NX.array([0], dtype=c_or_r.dtype)
self._coeffs = c_or_r
if variable is None:
variable = 'x'
self._variable = variable
def __array__(self, t=None):
if t:
return NX.asarray(self.coeffs, t)
else:
return NX.asarray(self.coeffs)
def __repr__(self):
vals = repr(self.coeffs)
vals = vals[6:-1]
return "poly1d(%s)" % vals
def __len__(self):
return self.order
def __str__(self):
thestr = "0"
var = self.variable
# Remove leading zeros
coeffs = self.coeffs[NX.logical_or.accumulate(self.coeffs != 0)]
N = len(coeffs)-1
def fmt_float(q):
s = '%.4g' % q
if s.endswith('.0000'):
s = s[:-5]
return s
for k, coeff in enumerate(coeffs):
if not iscomplex(coeff):
coefstr = fmt_float(real(coeff))
elif real(coeff) == 0:
coefstr = '%sj' % fmt_float(imag(coeff))
else:
coefstr = '(%s + %sj)' % (fmt_float(real(coeff)),
fmt_float(imag(coeff)))
power = (N-k)
if power == 0:
if coefstr != '0':
newstr = '%s' % (coefstr,)
else:
if k == 0:
newstr = '0'
else:
newstr = ''
elif power == 1:
if coefstr == '0':
newstr = ''
elif coefstr == 'b':
newstr = var
else:
newstr = '%s %s' % (coefstr, var)
else:
if coefstr == '0':
newstr = ''
elif coefstr == 'b':
newstr = '%s**%d' % (var, power,)
else:
newstr = '%s %s**%d' % (coefstr, var, power)
if k > 0:
if newstr != '':
if newstr.startswith('-'):
thestr = "%s - %s" % (thestr, newstr[1:])
else:
thestr = "%s + %s" % (thestr, newstr)
else:
thestr = newstr
return _raise_power(thestr)
def __call__(self, val):
return polyval(self.coeffs, val)
def __neg__(self):
return poly1d(-self.coeffs)
def __pos__(self):
return self
def __mul__(self, other):
if isscalar(other):
return poly1d(self.coeffs * other)
else:
other = poly1d(other)
return poly1d(polymul(self.coeffs, other.coeffs))
def __rmul__(self, other):
if isscalar(other):
return poly1d(other * self.coeffs)
else:
other = poly1d(other)
return poly1d(polymul(self.coeffs, other.coeffs))
def __add__(self, other):
other = poly1d(other)
return poly1d(polyadd(self.coeffs, other.coeffs))
def __radd__(self, other):
other = poly1d(other)
return poly1d(polyadd(self.coeffs, other.coeffs))
def __pow__(self, val):
if not isscalar(val) or int(val) != val or val < 0:
raise ValueError("Power to non-negative integers only.")
res = [1]
for _ in range(val):
res = polymul(self.coeffs, res)
return poly1d(res)
def __sub__(self, other):
other = poly1d(other)
return poly1d(polysub(self.coeffs, other.coeffs))
def __rsub__(self, other):
other = poly1d(other)
return poly1d(polysub(other.coeffs, self.coeffs))
def __div__(self, other):
if isscalar(other):
return poly1d(self.coeffs/other)
else:
other = poly1d(other)
return polydiv(self, other)
__truediv__ = __div__
def __rdiv__(self, other):
if isscalar(other):
return poly1d(other/self.coeffs)
else:
other = poly1d(other)
return polydiv(other, self)
__rtruediv__ = __rdiv__
def __eq__(self, other):
if not isinstance(other, poly1d):
return NotImplemented
if self.coeffs.shape != other.coeffs.shape:
return False
return (self.coeffs == other.coeffs).all()
def __ne__(self, other):
if not isinstance(other, poly1d):
return NotImplemented
return not self.__eq__(other)
def __getitem__(self, val):
ind = self.order - val
if val > self.order:
return self.coeffs.dtype.type(0)
if val < 0:
return self.coeffs.dtype.type(0)
return self.coeffs[ind]
def __setitem__(self, key, val):
ind = self.order - key
if key < 0:
raise ValueError("Does not support negative powers.")
if key > self.order:
zr = NX.zeros(key-self.order, self.coeffs.dtype)
self._coeffs = NX.concatenate((zr, self.coeffs))
ind = 0
self._coeffs[ind] = val
return
def __iter__(self):
return iter(self.coeffs)
def integ(self, m=1, k=0):
"""
Return an antiderivative (indefinite integral) of this polynomial.
Refer to `polyint` for full documentation.
See Also
--------
polyint : equivalent function
"""
return poly1d(polyint(self.coeffs, m=m, k=k))
def deriv(self, m=1):
"""
Return a derivative of this polynomial.
Refer to `polyder` for full documentation.
See Also
--------
polyder : equivalent function
"""
return poly1d(polyder(self.coeffs, m=m))

View File

@ -12,8 +12,8 @@ from org.meteoinfo.math.stats import StatsUtil
from .. import core as np
__all__ = ['solve', 'cholesky', 'cond', 'det', 'lu', 'qr', 'svd', 'eig', 'inv', 'lstsq', 'slogdet',
'solve_triangular', 'norm', 'pinv', 'LinAlgError']
__all__ = ['solve', 'cholesky', 'cond', 'det', 'lu', 'qr', 'svd', 'eig', 'eigvals', 'inv',
'lstsq', 'slogdet', 'solve_triangular', 'norm', 'pinv', 'LinAlgError']
class LinAlgError(Exception):
@ -224,7 +224,7 @@ def svd(a, full_matrices=True, compute_uv=True):
def eig(a):
"""
Compute the eigenvalues and right eigenvectors of a square array.
Parameters
----------
a : (M, M) array
@ -252,6 +252,39 @@ def eig(a):
return w, v
def eigvals(a):
"""
Compute the eigenvalues of a general matrix.
Main difference between `eigvals` and `eig`: the eigenvectors aren't
returned.
Parameters
----------
a : (M, M) array
Matrices for which the eigenvalues and right eigenvectors will
be computed
Returns
-------
w : (M) array
The eigenvalues, each repeated according to its multiplicity.
The eigenvalues are not necessarily ordered. The resulting
array will be of complex type, unless the imaginary part is
zero in which case it will be cast to a real type. When `a`
is real the resulting eigenvalues will be real (0 imaginary
part) or occur in conjugate pairs
v : (M, M) array
The normalized (unit "length") eigenvectors, such that the
column ``v[:,i]`` is the eigenvector corresponding to the
eigenvalue ``w[i]``.
"""
r = LinalgUtil.eigen(a.asarray())
# r = LinalgUtil.eigen_EJML(a.asarray())
w = np.NDArray(r[0])
return w
def inv(a):
"""
Compute the (multiplicative) inverse of a matrix.

View File

@ -1,6 +1,6 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<MeteoInfo File="config.xml" Type="configurefile">
<Path OpenPath="D:\Temp\grads"/>
<Path OpenPath="D:\Temp\traj\Sample"/>
<Font>
<TextFont FontName="YaHei Consolas Hybrid" FontSize="14"/>
<LegendFont FontName="宋体" FontSize="12"/>

View File

@ -102,11 +102,11 @@ public class MeteoMath {
/**
* Calculate relative humidity from dewpoint
*
* @param tc Temperature
* @param tdc Dewpoint temperature
* @param tc Temperature
* @return Relative humidity as percent (i.e. 80%)
*/
public static double dewpoint2rh(double tc, double tdc) {
public static double dewpoint2rh(double tdc, double tc) {
double es = cal_Es(tc);
double e = cal_Es(tdc);
return e / es * 100;
@ -125,7 +125,7 @@ public class MeteoMath {
IndexIterator iter_tdc = tdc.getIndexIterator();
IndexIterator iter_tc = tc.getIndexIterator();
while(iter.hasNext()) {
iter.setDoubleNext(MeteoMath.dewpoint2rh(iter_tc.getDoubleNext(), iter_tdc.getDoubleNext()));
iter.setDoubleNext(MeteoMath.dewpoint2rh(iter_tdc.getDoubleNext(), iter_tc.getDoubleNext()));
}
return r;

View File

@ -9276,6 +9276,148 @@ public class ArrayMath {
return r;
}
/**
* Returns the discrete, linear convolution of two one-dimensional sequences with full mode
*
* @param a First one-dimensional input array
* @param v Second one-dimensional input array
* @return Discrete, linear convolution of a and v
*/
public static Array convolveFull(Array a, Array v) {
a = a.copyIfView();
v = v.copyIfView();
int n = (int) a.getSize() + (int) v.getSize() - 1;
DataType dataType = commonType(a.getDataType(), v.getDataType());
Array r = Array.factory(dataType, new int[]{n});
if (dataType == DataType.COMPLEX) {
double real, real1, real2, image, image1, image2;
for (int i = 0; i < a.getSize(); i++) {
real1 = a.getComplex(i).real();
image1 = a.getComplex(i).imag();
for (int j = 0; j < v.getSize(); j++) {
real2 = v.getComplex(j).real();
image2 = v.getComplex(j).imag();
real = real1 * real2 - image1 * image2;
image = image1 * real2 + real1 * image2;
r.setComplex(i + j, r.getComplex(i + j).add(new Complex(real, image)));
}
}
} else {
for (int i = 0; i < a.getSize(); i++) {
for (int j = 0; j < v.getSize(); j++) {
r.setObject(i + j, r.getDouble(i + j) + a.getDouble(i) * v.getDouble(j));
}
}
}
return r;
}
/**
* Returns the discrete, linear convolution of two one-dimensional sequences with same mode
*
* @param a First one-dimensional input array
* @param v Second one-dimensional input array
* @return Discrete, linear convolution of a and v
*/
public static Array convolveSame(Array a, Array v) {
a = a.copyIfView();
v = v.copyIfView();
int n = (int) a.getSize();
DataType dataType = commonType(a.getDataType(), v.getDataType());
Array r = Array.factory(dataType, new int[]{n});
int nv = (int) v.getSize() / 2;
int idx;
if (dataType == DataType.COMPLEX) {
double real, real1, real2, image, image1, image2;
for (int i = 0; i < a.getSize(); i++) {
real1 = a.getComplex(i).real();
image1 = a.getComplex(i).imag();
for (int j = 0; j < v.getSize(); j++) {
idx = i + j - nv;
if (idx < 0) {
continue;
} else if (idx == n) {
break;
}
real2 = v.getComplex(j).real();
image2 = v.getComplex(j).imag();
real = real1 * real2 - image1 * image2;
image = image1 * real2 + real1 * image2;
r.setComplex(idx, r.getComplex(idx).add(new Complex(real, image)));
}
}
} else {
for (int i = 0; i < a.getSize(); i++) {
for (int j = 0; j < v.getSize(); j++) {
idx = i + j - nv;
if (idx < 0) {
continue;
} else if (idx == n) {
break;
}
r.setObject(idx, r.getDouble(idx) + a.getDouble(i) * v.getDouble(j));
}
}
}
return r;
}
/**
* Returns the discrete, linear convolution of two one-dimensional sequences with valid mode
*
* @param a First one-dimensional input array
* @param v Second one-dimensional input array
* @return Discrete, linear convolution of a and v
*/
public static Array convolveValid(Array a, Array v) {
a = a.copyIfView();
v = v.copyIfView();
int n = (int) a.getSize() - (int) v.getSize() + 1;
DataType dataType = commonType(a.getDataType(), v.getDataType());
Array r = Array.factory(dataType, new int[]{n});
int nv = (int) v.getSize();
int idx;
if (dataType == DataType.COMPLEX) {
double real, real1, real2, image, image1, image2;
for (int i = 0; i < a.getSize(); i++) {
real1 = a.getComplex(i).real();
image1 = a.getComplex(i).imag();
for (int j = 0; j < v.getSize(); j++) {
idx = i + j - nv + 1;
if (idx < 0) {
continue;
} else if (idx == n) {
break;
}
real2 = v.getComplex(j).real();
image2 = v.getComplex(j).imag();
real = real1 * real2 - image1 * image2;
image = image1 * real2 + real1 * image2;
r.setComplex(idx, r.getComplex(idx).add(new Complex(real, image)));
}
}
} else {
for (int i = 0; i < a.getSize(); i++) {
for (int j = 0; j < v.getSize(); j++) {
idx = i + j - nv + 1;
if (idx < 0) {
continue;
} else if (idx == n) {
break;
}
r.setObject(idx, r.getDouble(idx) + a.getDouble(i) * v.getDouble(j));
}
}
}
return r;
}
// </editor-fold>
// <editor-fold desc="Meteo">
/**