Source code for nutils.matrix

# Copyright (c) 2014 Evalf
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.

"""
The matrix module defines a number of 2D matrix objects, notably the
:func:`ScipyMatrix` and :func:`NumpyMatrix`. Matrix objects support basic
addition and subtraction operations and provide a consistent insterface for
solving linear systems. Matrices can be converted to numpy arrays via
``toarray`` or scipy matrices via ``toscipy``.
"""

from . import util, numpy, log, numeric, warnings
import functools


class MyCallback:

  def __init__ (self, matrix, rhs, tol, callback):
    self.matrix = matrix
    self.rhs = rhs
    self.norm = max(numpy.linalg.norm(rhs), 1) # scipy terminates on minimum of relative and absolute residual
    self.niter = 0
    self.tol = tol
    self.callback = callback

  def __call__(self, arg):
    self.niter += 1
    # some solvers provide the residual, others the left hand side vector
    res = float(numpy.linalg.norm(self.rhs - self.matrix * arg) / self.norm if numpy.ndim(arg) == 1 else arg)
    if self.callback:
      self.callback(res)
    with log.context('residual {:.2e} ({:.0f}%)'.format(res, 100. * numpy.log10(res) / numpy.log10(self.tol) if res > 0 else 0)):
      pass

[docs]class Matrix: 'matrix base class' def __init__(self, shape): 'constructor' assert len(shape) == 2 self.shape = shape @property def size(self): return numpy.prod(self.shape)
[docs] def cond(self, constrain=None, lconstrain=None, rconstrain=None): 'condition number' x, I, J = parsecons(constrain, lconstrain, rconstrain, self.shape) matrix = self.toarray()[numpy.ix_(I,J)] return numpy.linalg.cond(matrix)
[docs] def res(self, x, b=0, constrain=None, lconstrain=None, rconstrain=None, scaled=True): 'residual' x0, I, J = parsecons(constrain, lconstrain, rconstrain, self.shape) res = numpy.linalg.norm((self.matvec(x)-b)[I]) if scaled: res /= numpy.linalg.norm((self.matvec(x0)-b)[I]) return res
def clone(self): warnings.deprecation('warning: arrays are immutable; clone returns self for backwards compatibility') return self
[docs]class ScipyMatrix(Matrix): '''matrix based on any of scipy's sparse matrices''' def __init__(self, core): self.core = core super().__init__(core.shape) matvec = lambda self, vec: self.core.dot(vec) toarray = lambda self: self.core.toarray() toscipy = lambda self: self.core __add__ = lambda self, other: ScipyMatrix(self.core + (other.toscipy() if isinstance(other,Matrix) else other)) __sub__ = lambda self, other: ScipyMatrix(self.core - (other.toscipy() if isinstance(other,Matrix) else other)) __mul__ = lambda self, other: ScipyMatrix(self.core * (other.toscipy() if isinstance(other,Matrix) else other)) __radd__ = __add__ __rmul__ = __mul__ __div__ = lambda self, other: ScipyMatrix(self.core / other) T = property(lambda self: ScipyMatrix(self.core.transpose()))
[docs] def rowsupp(self, tol=0): 'return row indices with nonzero/non-small entries' supp = numpy.empty(self.shape[0], dtype=bool) for irow in range(self.shape[0]): a, b = self.core.indptr[irow:irow+2] supp[irow] = a != b and numpy.any(numpy.abs(self.core.data[a:b]) > tol) return supp
[docs] @log.title def solve(self, rhs=None, constrain=None, lconstrain=None, rconstrain=None, tol=0, lhs0=None, solver=None, symmetric=False, callback=None, precon=None, **solverargs): 'solve' import scipy.sparse.linalg lhs, I, J = parsecons(constrain, lconstrain, rconstrain, self.shape) A = self.core if not I.all(): A = A[I,:] b = (rhs[I] if rhs is not None else 0) - A * lhs if not J.all(): A = A[:,J] if not b.any(): log.info('right hand side is zero') return lhs if solver is None: solver = 'spsolve' if tol == 0 else 'cg' if symmetric else 'gmres' x0 = lhs0[J] if lhs0 is not None else numpy.zeros(A.shape[0]) res0 = numpy.linalg.norm(b - A * x0) / numpy.linalg.norm(b) if res0 < tol: log.info('initial residual is below tolerance') x = x0 elif solver == 'spsolve': log.info('solving system using sparse direct solver') x = scipy.sparse.linalg.spsolve(A, b) else: assert tol, 'tolerance must be specified for iterative solver' log.info('solving system using {} iterative solver'.format(solver)) solverfun = getattr(scipy.sparse.linalg, solver) # keep scipy from making things circular by shielding the nature of A A = scipy.sparse.linalg.LinearOperator(A.shape, A.__mul__, dtype=float) if isinstance(precon, str): M = self.getprecon(precon, constrain, lconstrain, rconstrain) elif callable(precon): M = precon(A) elif precon is None: # create identity operator, because scipy's native identity operator has circular references M = scipy.sparse.linalg.LinearOperator(A.shape, matvec=lambda x:x, rmatvec=lambda x:x, matmat=lambda x:x, dtype=float) else: M = precon assert isinstance(M, scipy.sparse.linalg.LinearOperator) mycallback = MyCallback(matrix=A, rhs=b, tol=tol, callback=callback) x, status = solverfun(A, b, M=M, tol=tol, x0=x0, callback=mycallback, **solverargs) assert status == 0, '{} solver failed with status {}'.format(solver, status) res = numpy.linalg.norm(b - A * x) / mycallback.norm log.info('solver converged in {} iterations to residual {:.1e}'.format(mycallback.niter, res)) lhs[J] = x return lhs
def getprecon(self, name='SPLU', constrain=None, lconstrain=None, rconstrain=None): import scipy.sparse.linalg name = name.lower() x, I, J = parsecons(constrain, lconstrain, rconstrain, self.shape) A = self.core[I,:][:,J] assert A.shape[0] == A.shape[1], 'constrained matrix must be square' log.info('building {} preconditioner'.format(name)) if name == 'splu': precon = scipy.sparse.linalg.splu(A.tocsc()).solve elif name == 'spilu': precon = scipy.sparse.linalg.spilu(A.tocsc(), drop_tol=1e-5, fill_factor=None, drop_rule=None, permc_spec=None, diag_pivot_thresh=None, relax=None, panel_size=None, options=None).solve elif name == 'diag': precon = numpy.reciprocal(A.diagonal()).__mul__ else: raise Exception('invalid preconditioner {!r}'.format(name)) return scipy.sparse.linalg.LinearOperator(A.shape, precon, dtype=float)
[docs]class NumpyMatrix(Matrix): '''matrix based on numpy array''' def __init__(self, core): assert numeric.isarray(core) self.core = core super().__init__(core.shape) matvec = lambda self, vec: numpy.dot(self.core, vec) toarray = lambda self: self.core __add__ = lambda self, other: NumpyMatrix(self.core + other.toarray()) __sub__ = lambda self, other: NumpyMatrix(self.core - other.toarray()) __mul__ = lambda self, other: NumpyMatrix(self.core * other) __rmul__ = __mul__ __div__ = lambda self, other: NumpyMatrix(self.core / other) T = property(lambda self: NumpyMatrix(self.core.T))
[docs] @log.title def solve(self, b=None, constrain=None, lconstrain=None, rconstrain=None, tol=0): 'solve' if b is None: b = numpy.zeros(self.shape[0]) else: b = numpy.asarray(b, dtype=float) assert b.ndim == 1, 'right-hand-side has shape {}, expected a vector'.format(b.shape) assert b.shape == self.shape[:1] x, I, J = parsecons(constrain, lconstrain, rconstrain, self.shape) if I.all() and J.all(): return numpy.linalg.solve(self.core, b) data = self.core[I] x[J] = numpy.linalg.solve(data[:,J], b[I] - numpy.dot(data[:,~J], x[~J])) return x
# UTILITY FUNCTIONS
[docs]def assemble(data, index, shape, force_dense=False): 'create data from values and indices' if len(shape) == 0: retval = data.sum() elif len(shape) == 2 and not force_dense: import scipy.sparse.linalg csr = scipy.sparse.csr_matrix((data,index), shape) retval = ScipyMatrix(csr) else: flatindex = numpy.dot(numpy.cumprod((1,)+shape[:0:-1])[::-1], index) retval = numpy.bincount(flatindex, data, numpy.prod(shape)).reshape(shape).astype(data.dtype, copy=False) if retval.ndim == 2: retval = NumpyMatrix(retval) assert retval.shape == shape log.debug('assembled {}({})'.format(retval.__class__.__name__, ','.join(str(n) for n in shape))) return retval
[docs]def parsecons(constrain, lconstrain, rconstrain, shape): 'parse constraints' I = numpy.ones(shape[0], dtype=bool) x = numpy.empty(shape[1]) x[:] = numpy.nan if constrain is not None: assert lconstrain is None assert rconstrain is None assert numeric.isarray(constrain) I[:] = numpy.isnan(constrain) x[:] = constrain if lconstrain is not None: assert numeric.isarray(lconstrain) x[:] = lconstrain if rconstrain is not None: assert numeric.isarray(rconstrain) I[:] = rconstrain J = numpy.isnan(x) assert numpy.sum(I) == numpy.sum(J), 'constrained matrix is not square: {}x{}'.format(numpy.sum(I), numpy.sum(J)) x[J] = 0 return x, I, J
# vim:shiftwidth=2:softtabstop=2:expandtab:foldmethod=indent:foldnestmax=2