import collections
import numpy as np
from numpy.linalg import matrix_rank
import scipy.signal as sig
from scipy.linalg import block_diag, lu, matrix_balance, solve, lstsq
from ._aux_linalg import haroldsvd, e_i
__all__ = ['haroldlcm', 'haroldgcd', 'haroldcompanion', 'haroldpoly',
'haroldpolyadd', 'haroldpolymul', 'haroldpolydiv']
def haroldlcm(*args, compute_multipliers=True, cleanup_threshold=1e-9):
"""
Takes n-many 1D numpy arrays and computes the numerical
least common multiple polynomial. The polynomials are
assumed to be in decreasing powers, e.g. s^2 + 5 should
be given as ``[1,0,5]``
Returns a numpy array holding the polynomial coefficients
of LCM and a list, of which entries are the polynomial
multipliers to arrive at the LCM of each input element.
For the multiplier computation, a variant of [1]_ is used.
Parameters
----------
args : iterable
Input arrays. 1-D arrays or array_like sequences of polynomial
coefficients
compute_multipliers : bool, optional
After the computation of the LCM, this switch decides whether the
multipliers of the given arguments should be computed or skipped.
A multiplier in this context is ``[1,3]`` for the argument ``[1,2]``
if the LCM turns out to be ``[1,5,6]``.
cleanup_threshold : float
The computed polynomials might contain some numerical noise and after
finishing everything this value is used to clean up the tiny entries.
Set this value to zero to turn off this behavior. The default value
is :math:`10^{-9}`.
Returns
--------
lcmpoly : ndarray
Resulting 1D polynomial coefficient array for the LCM.
mults : list
The multipliers given as a list of 1D arrays, for each given argument.
Notes
-----
If complex-valued arrays are given, only real parts are taken into account.
Examples
--------
>>> a , b = haroldlcm([1,3,0,-4], [1,-4,-3,18], [1,-4,3], [1,-2,-8])
>>> a
array([ 1., -7., 3., 59., -68., -132., 144.]
>>> b
[array([ 1., -10., 33., -36.]),
array([ 1., -3., -6., 8.]),
array([ 1., -3., -12., 20., 48.]),
array([ 1., -5., 1., 21., -18.])]
>>> np.convolve([1, 3, 0, -4], b[0]) # or haroldpolymul() for poly mult
(array([ 1., -7., 3., 59., -68., -132., 144.]),
References
----------
.. [1] Karcanias, Mitrouli, "System theoretic based characterisation and
computation of the least common multiple of a set of polynomials",
2004, :doi:`10.1016/j.laa.2003.11.009`
"""
# Regularize the arguments
args = [np.array(a).squeeze().real for a in args]
# Add dimension if any scalar arrays such as np.array(1)
args = [a if a.ndim > 0 else np.atleast_1d(a) for a in args]
if not all([x.ndim == 1 for x in args]):
raise ValueError('Input arrays must be 1D.')
if not all([x.size > 0 for x in args]):
raise ValueError('Empty arrays are not allowed.')
# All scalars
if all([x.size == 1 for x in args]):
if compute_multipliers:
return np.array([1.]), [np.array([1.]) for _ in range(len(args))]
else:
return np.array([1.])
# Remove if there are constant polynomials but return their multiplier!
poppedargs = [x for x in args if x.size > 1]
# Get the index number of the ones that are popped
p_ind, l_ind = [], []
[p_ind.append(ind) if x.size == 1 else l_ind.append(ind)
for ind, x in enumerate(args)]
# If there are more than one nonconstant polynomial to consider
if len(poppedargs) > 1:
a = block_diag(*(map(haroldcompanion, poppedargs)))
b = np.concatenate([e_i(x.size-1, -1) for x in poppedargs])
n = a.shape[0]
# Balance A
As, (sca, _) = matrix_balance(a, permute=False, separate=True)
Bs = b*np.reciprocal(sca)[:, None]
# Computing full c'bility matrix is redundant we just need to see where
# the rank drop is (if any!). Due to matrix power, things grow quickly!
C = Bs
for _ in range(n-1):
C = np.hstack([C, As @ C[:, [-1]]])
if matrix_rank(C) != C.shape[1]:
break
else:
# No break
C = np.hstack([C, As @ C[:, [-1]]])
cols = C.shape[1]
_, s, v = haroldsvd(C)
temp = s @ v
lcmpoly = solve(temp[:cols-1, :-1], -temp[:cols-1, -1])
# Add monic coefficient and flip
lcmpoly = np.append(lcmpoly, 1)[::-1]
else:
lcmpoly = np.trim_zeros(poppedargs[0], 'f')
lcmpoly = lcmpoly/lcmpoly[0]
if compute_multipliers:
n_lcm = lcmpoly.size - 1
if len(poppedargs) > 1:
c = block_diag(*[e_i(x.size-1, 0).T for x in poppedargs]) * sca
b_lcm, _, _, _ = lstsq(C[:c.shape[1], :-1], Bs)
c_lcm = c @ C[:c.shape[1], :-1]
# adj(sI-A) formulas with A being a companion matrix. Use a 3D
# array where x,y,z = adj(sI-A)[x,y] and z is the coefficient array
adjA = np.zeros([n_lcm, n_lcm, n_lcm])
# fill in the adjoint
for x in range(n_lcm):
# Diagonal terms
adjA[x, x, :n_lcm-x] = lcmpoly[:n_lcm-x]
for y in range(n_lcm):
if y < x: # Upper Triangular terms
adjA[x, y, x-y:] = adjA[x, x, :n_lcm-(x-y)]
elif y > x: # Lower Triangular terms
adjA[x, y, n_lcm-y:n_lcm+1-y+x] = \
-lcmpoly[-x-1:n_lcm+1]
# C*adj(sI-A)*B
mults = c_lcm @ np.sum(adjA * b_lcm, axis=1)
else:
mults = np.zeros((1, n_lcm))
mults[0, -1] = 1.
if len(p_ind) > 0:
temp = np.zeros((len(args), lcmpoly.size), dtype=float)
temp[p_ind] = lcmpoly
temp[l_ind, 1:] = mults
mults = temp
lcmpoly[abs(lcmpoly) < cleanup_threshold] = 0.
mults[abs(mults) < cleanup_threshold] = 0.
mults = [np.trim_zeros(z, 'f') for z in mults]
return lcmpoly, mults
else:
return lcmpoly
def haroldgcd(*args):
"""
Takes 1D numpy arrays and computes the numerical greatest common
divisor polynomial. The polynomials are assumed to be in decreasing
powers, e.g. :math:`s^2 + 5` should be given as ``numpy.array([1,0,5])``.
Returns a numpy array holding the polynomial coefficients
of GCD. The GCD does not cancel scalars but returns only monic roots.
In other words, the GCD of polynomials :math:`2` and :math:`2s+4` is
still computed as :math:`1`.
Parameters
----------
args : iterable
A collection of 1D array_likes.
Returns
--------
gcdpoly : ndarray
Computed GCD of args.
Examples
--------
>>> a = haroldgcd(*map(haroldpoly,([-1,-1,-2,-1j,1j],
[-2,-3,-4,-5],
[-2]*10)))
>>> a
array([ 1., 2.])
.. warning:: It uses the LU factorization of the Sylvester matrix.
Use responsibly. It does not check any certificate of
success by any means (maybe it will in the future).
I have played around with ERES method but probably due
to my implementation, couldn't get satisfactory results.
I am still interested in better methods.
"""
raw_arr_args = [np.atleast_1d(np.squeeze(x)) for x in args]
arr_args = [np.trim_zeros(x, 'f') for x in raw_arr_args if x.size > 0]
dimension_list = [x.ndim for x in arr_args]
# do we have 2d elements?
if max(dimension_list) > 1:
raise ValueError('Input arrays must be 1D arrays, rows, or columns')
degree_list = np.array([x.size-1 for x in arr_args])
max_degree = np.max(degree_list)
max_degree_index = np.argmax(degree_list)
try:
# There are polynomials of lesser degree
second_max_degree = np.max(degree_list[degree_list < max_degree])
except ValueError:
# all degrees are the same
second_max_degree = max_degree
n, p, h = max_degree, second_max_degree, len(arr_args) - 1
# If a single item is passed then return it back
if h == 0:
return arr_args[0]
if n == 0:
return np.array([1])
if n > 0 and p == 0:
return arr_args.pop(max_degree_index)
# pop out the max degree polynomial and zero pad
# such that we have n+m columns
S = np.array([np.hstack((
arr_args.pop(max_degree_index),
np.zeros((1, p-1)).squeeze()
))]*p)
# Shift rows to the left
for rows in range(S.shape[0]):
S[rows] = np.roll(S[rows], rows)
# do the same to the remaining ones inside the regular_args
for item in arr_args:
_ = np.array([np.hstack((item, [0]*(n+p-item.size)))]*(
n+p-item.size+1))
for rows in range(_.shape[0]):
_[rows] = np.roll(_[rows], rows)
S = np.r_[S, _]
rank_of_sylmat = np.linalg.matrix_rank(S)
if rank_of_sylmat == min(S.shape):
return np.array([1])
else:
p, l, u = lu(S)
u[abs(u) < 1e-8] = 0
for rows in range(u.shape[0]-1, 0, -1):
if not any(u[rows, :]):
u = np.delete(u, rows, 0)
else:
break
gcdpoly = np.real(np.trim_zeros(u[-1, :], 'f'))
# make it monic
gcdpoly /= gcdpoly[0]
return gcdpoly
def haroldcompanion(somearray):
"""
Takes a 1D numpy array or list and returns the companion matrix
of the monic polynomial of somearray. Hence ``[0.5,1,2]`` will be first
converted to ``[1,2,4]``.
Examples
--------
>>> haroldcompanion([2,4,6])
array([[ 0., 1.],
[-3., -2.]])
>>> haroldcompanion([1,3])
array([[-3.]])
>>> haroldcompanion([1])
array([], dtype=float64)
"""
if not isinstance(somearray, (list, np.ndarray)):
raise TypeError('Companion matrices are meant only for '
'1D lists or 1D Numpy arrays. I found '
'a \"{0}\"'.format(type(somearray).__name__))
if len(somearray) == 0:
return np.array([])
# regularize to flat 1D np.array
somearray = np.array(somearray).flatten()
ta = np.trim_zeros(somearray, 'f')
# convert to monic polynomial.
# Note: ta *=... syntax doesn't autoconvert to float
ta = np.array(1/ta[0])*ta
ta = -ta[-1:0:-1]
n = ta.size
if n == 0: # Constant polynomial
return np.array([])
elif n == 1: # First-order --> companion matrix is a scalar
return np.atleast_2d(np.array(ta))
else: # Other stuff
return np.vstack((np.hstack((np.zeros((n-1, 1)), np.eye(n-1))), ta))
def haroldpoly(rootlist):
"""
Takes a 1D array-like numerical elements as roots and forms the polynomial
"""
if isinstance(rootlist, collections.abc.Iterable):
r = np.array([x for x in rootlist], dtype=complex)
else:
raise TypeError('The argument must be something iterable,\nsuch as '
'list, numpy array, tuple etc. I don\'t know\nwhat '
'to do with a \"{0}\" object.'
''.format(type(rootlist).__name__))
n = r.size
if n == 0:
return np.ones(1)
else:
p = np.array([0.+0j for x in range(n+1)], dtype=complex)
p[0] = 1 # Monic polynomial
p[1] = -rootlist[0]
for x in range(1, n):
p[x+1] = -p[x]*r[x]
for y in range(x, 0, -1):
p[y] -= p[y-1] * r[x]
return p
def haroldpolyadd(*args, trim_zeros=True):
"""
A wrapper around NumPy's :func:`numpy.polyadd` but allows for multiple
args and offers a trimming option.
Parameters
----------
args : iterable
An iterable with 1D array-like elements.
trim_zeros : bool, optional
If True, the zeros at the front of the input and output arrays are
truncated. Default is True.
Returns
-------
p : ndarray
The polynomial coefficients of the sum.
Examples
--------
>>> a = np.array([2, 3, 5, 8])
>>> b = np.array([1, 3, 4])
>>> c = np.array([6, 9, 10, -8, 6])
>>> haroldpolyadd(a, b, c)
array([ 6., 11., 14., 0., 18.])
>>> d = np.array([-2, -4 ,3, -1])
>>> haroldpolyadd(a, b, d)
array([ 0., 0., 11., 11.])
>>> haroldpolyadd(a, b, d, trim_zeros=False)
array([ 0., 0., 11., 11.])
"""
if trim_zeros:
trimmedargs = [np.trim_zeros(x, 'f') for x in args]
else:
trimmedargs = args
degs = [len(m) for m in trimmedargs] # Get the max len of args
s = np.zeros((1, max(degs)))
for ind, x in enumerate(trimmedargs):
s[0, max(degs)-degs[ind]:] += np.real(x)
return s[0]
def haroldpolymul(*args, trim_zeros=True):
"""
Simple wrapper around the :func:`numpy.convolve` function for polynomial
multiplication with multiple args. The arguments are passed through
the left zero trimming function first.
See Also
--------
haroldpolydiv, :func:`numpy.convolve`, :func:`scipy.signal.convolve`
Parameters
----------
args : iterable
An iterable with 1D array-like elements.
trim_zeros : bool, optional
If True, the zeros at the front of the arrays are truncated.
Default is True.
Returns
-------
p : ndarray
The polynomial coefficients of the product.
Examples
--------
>>> haroldpolymul([0,2,0], [0,0,0,1,3,3,1], [0,0.5,0.5])
array([ 1., 4., 6., 4., 1., 0.])
"""
if trim_zeros:
trimmedargs = [np.trim_zeros(x.flatten(), 'f') for x in args]
else:
trimmedargs = args
p = trimmedargs[0]
for x in trimmedargs[1:]:
if x.size == 0: # it was all zeros
p = []
break
p = np.convolve(p, x)
return p if np.any(p) else np.array([0.])
def haroldpolydiv(dividend, divisor):
"""
Polynomial division wrapped around :func:`scipy.signal.deconvolve`
function. Takes two arguments and divides the first
by the second.
Parameters
----------
dividend : (n,) array_like
The polynomial to be divided
divisor : (m,) array_like
The polynomial that divides
Returns
-------
factor : ndarray
The resulting polynomial coeffients of the factor
remainder : ndarray
The resulting polynomial coefficients of the remainder
Examples
--------
>>> a = np.array([2, 3, 4 ,6])
>>> b = np.array([1, 3, 6])
>>> haroldpolydiv(a, b)
(array([ 2., -3.]), array([ 1., 24.]))
>>> c = np.array([1, 3, 3, 1])
>>> d = np.array([1, 2, 1])
>>> haroldpolydiv(c, d)
(array([1., 1.]), array([], dtype=float64))
"""
h_factor, h_remainder = (np.trim_zeros(x, 'f') for x
in sig.deconvolve(dividend, divisor))
return h_factor, h_remainder
