"""Solving linear equations in integer coefficients.
This is desirable for unit-related calculations, when inferring how a set of
values, of known kinds, imply units for values of assorted other kinds. For
this use, the addition and integer scaling of our linear system correspond to
multiplication and taking powers of the kinds of value. Each available vector
describes the kind of one of the given values by the powers of the 'base kinds'
that make up the value's kind. The standard base units (i.e. canonical basis)
may then be inferred from the values by using the linear system's inverse.
$Id: reduce.py,v 1.12 2007/06/03 16:43:10 eddy Exp $
"""
import natural, permute
from study.snake.lazy import Lazy
class System (Lazy):
"""Analyzer for integer-valued linear systems.
The system is described by a canonical basis and a set of available vectors
expressed in terms of this basis. The object representing the system does
what it can to express the canonical basis in terms of the available
vectors. Where the available vectors' span coincides with that of some
subset of the canonical basis, the object expresses this subset of the
canonical basis in terms of it, via the .inverse attribute. Otherwise,
attempting to access .inverse shall get you a value error and the object
does the best it can, via other attributes (to be invented).
Each matrix attribute is encoded as a tuple of rows; each row is itself a
tuple giving co-ordinates of the row with respect to either the canonical
basis or the available vectors. In matrix attributes marked as
'(rational)', each row's last entry is an implicit denominator for all
earlier entries; others are simple integer sequences, with no extra entry on
each row. Aside from .problem, all attributes are lazily-evaluated.
Public attributes:
problem -- cleaned-up form of the constructor's input (rational)
inverse -- 'inverse' of .problem (rational), where available; each row,
.inverse[i], specifies a linear combination of available
vectors which yields the canonical basis member with index i,
where that is accessible; otherwise, inaccessibility is
signalled by a denominator of 0 and the rest of the row is
taken from kernel or, if those run out, is all zero. Where two
canonical basis members are inseparable, you'll get a
ValueError (there's only so much I'm willing to kludge).
available -- each row is a vector, in terms of the canonical basis,
which can be expressed in terms of the available vectors
recipe -- each row specifies the linear comibination of the available
vectors that yields the matching row of .available
kernel -- each row is a linear dependency among the available vectors
solution -- tuple, (.available, .recipe, .kernel), describing what can be
achieved.
The rows of .available and of .recipe are linearly independent, as are those
of .inverse with non-zero denominator, when available. Note that .kernel is
not guaranteed to be exhaustive, or linearly independent, but it's quite
likely to be both.
Although dealing notionally with integer-valued data, this class actually
(via the implicit denominator mechanism) supports rational values; general
integer-valued problems can perfectly readilly imply solutions involving
rationals, so it makes sense to support the converse case. The crucial thing
is that all arithmetic is approached from an exact integer perspective,
rather than using floating-point approximations.\n"""
def __init__(self, n, *rows):
"""Prepare to analyze an integer-valued linear system.
First input, n, is the dimension of a linear space; each subsequent
input is a sequence of integers, denoting a vector in that linear space.
Any [n] entry of such a sequence denotes a denominator by which the [:n]
entries are implicitly divided. Any sequence shorter than n is
implicitly padded with zero entries to length n, and given an [n] entry
of 1. The sequences are not modified. The new object retains copies of
them, so subsequent modifications to them shall not affect the new
object.\n"""
self.__ca = n, len(rows) # numbers of vectors: in canonical basis; and available
self.problem = self.__ingest(n, rows)
def contract(self, row):
"""Weighted sum of available values.
Single argument is a sequence of coefficients for the available vectors
supplied to the constructor; each entry is multiplied by the matching
available vector, the results are summed and a rational vector
(i.e. final entry is an implied denominator for the earlier ones)
relative to the canonical basis is returned. If the given sequence has
fewer entries than there are available values, (a copy of) it is padded
with zeros; otherwise, it may have one extra value, indicating an
implied denominator for all the others.\n"""
n = self.__ca[1]
if len(row) > 1 + n:
raise ValueError('too many entries', row, n)
elif len(row) < n:
row = list(row)
while len(row) < n: row.append(0)
if len(row) > n: row, scale = row[:n], row[n]
else: scale = 1
return self.__lincomb(row, scale)
def obtain(self, row, scale=1):
"""Find how to combine availables to yield a given member of our linear space.
Singe argument, row, is a member of the linear space; returns a sequence
of coefficients which, if fed to .contract(), would yield this row as
output.\n"""
n = self.__ca[0]
if len(row) > 1 + n:
raise ValueError('too many entries', row, n)
elif len(row) < n:
row = list(row)
while len(row) < n: row.append(0)
if len(row) > n: row, scale = row[:n], row[n]
else: scale = 1
return self.__obtain(row, scale)
def check(self):
"""Compute solution and test for correctness.
The test is done via assertions, so only works if debug is enabled !\n"""
dim, ava = self.__ca
try: inv = self.inverse
except ValueError: pass
else:
assert ava >= dim == len(inv)
i = dim
while i > 0:
i -= 1
self.__check(i, inv[i])
av, re, ker = self.solution # re * problem == av
i = len(re)
while i > 0:
i -= 1
self.__confirm(re[i], av[i])
i, nul = len(ker), ( 0, ) * dim
while i > 0:
i -= 1
self.__confirm(ker[i], nul)
if len(re) + len(ker) != len(self.problem):
print "It surprises me that kernel", len(ker), "and rank", len(re), \
"don't add up to system size", len(self.problem)
# but this isn't an error
def _lazy_get_available_(self, ig): return self.solution[0]
def _lazy_get_recipe_(self, ig): return self.solution[1]
def _lazy_get_kernel_(self, ig): return self.solution[2]
def _lazy_get_solution_(self, ig):
self.__setup()
self.__upper()
self.__diag()
return self.__tidy()
def freeze(row): return tuple(map(tuple, row)) # local function, del'd later
def _lazy_get_inverse_(self, ig, gcd=natural.hcf, safe=freeze):
can, how, ker = self.solution
how = map(list, how) # deep copy; it'll be our result
# Conveniently, .__tidy() is sure to have left this diagonal if invertible.
dim, ava = self.__ca
i, j, k, co = 0, 0, 0, []
while i < dim:
row = can[j]
if row[i]:
j += 1
how[i].append(row[i])
f = apply(gcd, how[i])
if f and how[i][-1] < 0: f = -f
if f > 1 or f < 0:
k = len(how[i])
while k > 0:
k -= 1
how[i][k] /= f
if filter(None, row[:i] + row[1+i:]):
co.append(j)
elif k < len(ker):
how.insert(i, ker[k] + (0,))
k += 1
else:
how.insert(i, (0,) * (ava + 1))
i += 1
if co:
raise ValueError('Irreducible composites', co)
return safe(how)
def __ingest(self, n, rows, gcd=natural.hcf, safe=freeze):
ans = []
for row in rows:
if len(row) > 1 + n:
raise ValueError('over-long row', row)
r = list(row)
while len(r) < n: r.append(0)
if len(r) == n: r.append(1)
elif not r[n]:
raise ValueError('zero denominator on row', row)
else:
# The two prior cases make f be 0 or 1, so don't need this:
f = apply(gcd, r)
if r[n] < 0: f = -f
if f > 1 or f < 0:
i = len(r)
while i > 0:
i -= 1
r[i] /= f
ans.append(r)
return safe(ans)
def __confirm(self, left, out):
# left times column j of problem should yield out[j]
j = self.__ca[0]
assert j == len(out)
while j > 0:
j -= 1
tot, den = self.__contract(j, left)
assert tot == out[j] * den
def __check(self, i, row):
# row i of inverse times column j of problem should yield 0, or 1 when i == j
# Final entry in each row of each matrix is a denominator for that row.
row, scale = row[:-1], row[-1]
j = self.__ca[0]
while j > 0:
j -= 1
tot, den = self.__contract(j, row) # sum of product entry, implicitly as tot/den
if j == i:
assert tot == den * scale, ("Non-unit diagonal entry", i, tot, den, scale)
else:
assert tot == 0, ("Non-zero off-diagonal entry", i, j, tot, den, scale)
def dragwith(matrix, j, row, gcd=natural.gcd):
k, tot, den = len(row), 0, 1
assert k == len(matrix)
assert j < len(matrix[0])
while k > 0:
k -= 1
right = matrix[k]
v, s = row[k] * right[j], right[-1]
tot = tot * s + v * den
den *= s
f = gcd(tot, den)
if f > 1: tot, den = tot / f, den / f
return tot, den
def denominate(pairs, scale lcm=natural.lcm, gcd=natural.hcf):
den = apply(lcm, map(lambda (n, d): d, pairs))
ans = map(lambda (n, d), s=den: n * s / d, pairs) + [ den * scale ]
f = apply(gcd, ans)
if ans[-1] < 0: f = -f
if f > 1 or f < 0:
i = len(ans)
while i > 0:
i -= 1
ans[i] /= f
return tuple(ans)
def __contract(self, j, row, dot=dragwith):
"""Contract column j of matrix with given row."""
return dot(self.problem, j, row)
def __lincomb(self, row, scale, comb=denominate):
"""As for .contract(), but with row of correct length and scale separate."""
return comb(map(lambda j, c=self.__contract, r=row: c(j, r),
range(self.__ca[0])), scale)
def __obtain(self, row, scale=1, comb=denominate, dot=dragwith):
"""As for .obtain(), but with row of correct length and scale separate."""
return comb(map(lambda j, r=row, d=dot, m=self.inverse: d(m, j, r),
range(len(self.inverse[0]))), scale)
del dragwith, denominate
# The actual analysis of the problem:
def __tidy(self, order=permute.order, shuffle=permute.permute, safe=freeze):
"""Tidy-up after attempted diagonalization.
We have .__matrix as near as we can hope for to diagonal form, albeit as
seen through the distorting lens of .__column; we need to extract from
it, and from the accompanying .__result, the information we were looking
for.\n"""
avail, recip, degen = self.__matrix, self.__result, []
del self.__matrix, self.__result, self.__column
# Purge any zero rows:
i = len(avail)
assert i == len(recip)
while i > 0:
i -= 1
if filter(None, avail[i]):
assert filter(None, recip[i])
else:
if filter(None, recip[i]):
degen.append(recip[i])
del avail[i], recip[i]
self.kernel = safe(degen)
# Order rows by increasing length of initial sequence of zeros:
indent, i = [], 0
while i < len(avail):
j, row = 0, avail[i]
while not row[j]: j += 1 # we know row has at least one non-zero entry.
indent.append(j)
i += 1
perm = order(indent)
assert len(avail) == len(perm) == len(recip)
self.available = safe(shuffle(avail, perm))
self.recipe = safe(shuffle(recip, perm))
return self.available, self.recipe, self.kernel
del freeze
def scaling(row):
nul, ans, i = ( 0, ) * len(row), [], 0
while i < len(row):
it = list(nul)
it[i] = row[i]
ans.append(it)
i += 1
return ans
def __setup(self, diag=scaling):
"""Prepare for analysis of our linear system.
Uses .problem's denominators as diagonal entries for (otherwise zero)
.__result and the rest of .problem as .__matrix; sets up __column as the
identity permutation.
When computation is trying to make it be upper triangular and, later,
diagonal, .__matrix only appears so (in so far as we succeed) when
viewed via the permutation of its columns encoded in .__column; this
view of .__matrix is what .__entry() implements.
If we multiply .__matrix on its right by a column whose entries are the
members of our canonical basis, the result is a column of vectors. If
we multiply .__result on its right by a column whose entries are our
initial available vectors, we get the same column of vectors. This
column is initially just the available vectors, scaled if necessary to
have integer components. Note that it is this column of vectors that is
permuted by .__swap(), not the basis or availables.\n"""
self.__matrix = map(lambda r: list(r[:-1]), self.problem)
self.__result = diag(map(lambda r: r[-1], self.problem))
self.__column = range(self.__ca[0])
del scaling
def __entry(self, i, j): return self.__matrix[i][self.__column[j]]
def __diag(self):
"""Reduce .__matrix from upper triangular form to diagonal form."""
i, ava = self.__ca
if i > ava: i = ava
while i > 0:
i -= 1
num = self.__entry(i, i)
if not num: continue # can't change anything usefully
j = ava
while j > 0:
j -= 1
if j == i: continue # leave that row itself alone !
n = self.__entry(j, i)
if n:
assert j > self.__ca[0] or j < i, ('expected zeros', i)
self.__take(i, n, j, num)
def __upper(self):
"""Reduce .__matrix to upper triangular form."""
dim, ava = self.__ca
if dim > ava: dim = ava
i = 0
while i < dim:
self.__select(i)
key = self.__entry(i, i)
if not key:
break # Remainder of .__matrix is zero, so upper triangular already
# Use .__matrix[i] to clear .__entry(j, i) for all j > i:
j = i + 1
while j < ava:
q = self.__entry(j, i)
if q: # else: nothing to do
self.__take(i, q, j, key)
assert not self.__entry(j, i), ('Bad arithmetic', i, self.__entry(j, i))
j += 1
i += 1
# TODO: study choices for peak; natural.lcm, number of non-zero entries ...
def __select(self, i, peak=lambda r: max(map(abs, r))):
"""Find a good candidate for self.__entry(i, i).
This prefers a row with low, but non-zero, peak entry; and selects a
smallest non-zero entry in this row to be the nominal i-th entry on the
diagonal of __matrix.\n"""
dim, ava = self.__ca
# Find row with smallest peak, swap it into row i:
j, k, m = i+1, i, peak(self.__matrix[i])
while j < ava and m != 1:
n = peak(self.__matrix[j])
if n and (not m or n < m):
k, m = j, n
j += 1
if not m: return # all remaining entries are zero !
if k != i: self.__swap(i, k)
# Find smallest non-zero entry in row, swap it into virtual column i
j, k, m = i+1, i, abs(self.__entry(i, i))
while j < dim and m != 1:
n = abs(self.__entry(i, j))
if n and (not m or n < m):
k, m = j, n
j += 1
assert m, 'How did we find a non-zero lcm of this row, to like ?'
if k != i:
col = self.__column
col[k], col[i] = col[i], col[k]
def __swap(self, i, j):
"""Swap rows i and j of the system."""
for seq in self.__matrix, self.__result:
seq[i], seq[j] = seq[j], seq[i]
def __take(self, i, ni, j, nj, gcd=natural.hcf):
"""Replace row j of system by a linear combination with row i."""
for seq in self.__matrix, self.__result:
top, row = seq[i], seq[j]
k = len(row)
while k > 0:
k -= 1
row[k] = row[k] * nj - top[k] * ni
f = apply(gcd, self.__matrix[j] + self.__result[j])
r = filter(None, self.__matrix[j])
if r and r[0] < 0: f = -f
if f > 1 or f < 0:
mat, row = self.__matrix[j], self.__result[j]
k = len(row)
while k > 0:
k -= 1
row[k] /= f
mat[k] /= f
del Lazy, permute, natural