mantis/grandalf/utils/linalg.py

from math import sqrt
from array import array as _array

constants = (int, float)


def coerce_(types):
    if types is None:
        types = []
    if str in types:
        raise TypeError
    if complex in types:
        raise TypeError
    dtype = ("i", int)
    if float in types:
        dtype = ("d", float)
    return dtype


def _mkslice(i, n):
    if not isinstance(i, slice):
        i = slice(i, i + 1, 1)
    start, stop, stride = i.indices(n)
    return slice(start, stop, stride)


def make_ij_slices(f):
    def wrapper(self, ij, *args):
        # I = slice(0,self.n,1)
        J = slice(0, self.p, 1)
        if isinstance(ij, tuple):
            I = _mkslice(ij[0], self.n)
            J = _mkslice(ij[1], self.p)
        else:
            I = _mkslice(ij, self.n)
        return f(self, (I, J), *args)

    return wrapper


# minimalistic numpy.array replacement class used as fallback
# when numpy is not found in geometry module
class array(object):
    def __init__(self, data, dtype=None, copy=True):
        self.dim = len(data)
        types = None
        if self.dim > 0:
            types = set([type(x) for x in data])
        if dtype is not None:
            types = (dtype,)
        tc, self.dtype = coerce_(types)
        data = [self.dtype(x) for x in data]
        if copy is True:
            self.data = _array(tc, data)
        else:
            raise NotImplementedError

    def coerce(self, dtype):
        data = [dtype(x) for x in self.data]
        tc, dtype = coerce_((dtype,))
        self.data = _array(tc, data)
        self.dtype = dtype

    @property
    def typecode(self):
        return self.data.typecode

    def __len__(self):
        return self.dim

    def __str__(self):
        s = " ".join(("%.12s" % x).ljust(12) for x in self)
        return "[%s]" % s.strip()

    def copy(self):
        return array(self.data, self.dtype)

    def __add__(self, v):
        if isinstance(v, constants):
            v = array([v] * self.dim)
        assert v.dim == self.dim
        return array([x + y for (x, y) in zip(self.data, v.data)])

    def __sub__(self, v):
        if isinstance(v, constants):
            v = array([v] * self.dim)
        assert v.dim == self.dim
        return array([x - y for (x, y) in zip(self.data, v.data)])

    def __neg__(self):
        return array([-x for x in self.data], dtype=self.dtype)

    def __radd__(self, v):
        return self + v

    def __rsub__(self, v):
        return (-self) + v

    def dot(self, v):
        assert v.dim == self.dim
        return sum([x * y for (x, y) in zip(self.data, v.data)])

    def __rmul__(self, k):
        return array([k * x for x in self.data])

    def __mul__(self, v):
        if isinstance(v, constants):
            v = array([v] * self.dim)
        assert v.dim == self.dim
        return array([x * y for (x, y) in zip(self.data, v.data)])

    def __truediv__(self, v):
        if isinstance(v, constants):
            v = array([v] * self.dim)
        assert v.dim == self.dim
        return array([x / y for (x, y) in zip(self.data, v.data)])

    __div__ = __truediv__

    def __rtruediv__(self, v):
        if isinstance(v, constants):
            v = array([v] * self.dim)
        assert v.dim == self.dim
        return array([x / y for (x, y) in zip(v.data, self.data)])

    __rdiv__ = __rtruediv__

    def __floordiv__(self, v):
        if isinstance(v, constants):
            v = array([v] * self.dim)
        assert v.dim == self.dim
        return array([x // y for (x, y) in zip(self.data, v.data)])

    def __rfloordiv__(self, v):
        if isinstance(v, constants):
            v = array([v] * self.dim)
        assert v.dim == self.dim
        return array([x // y for (x, y) in zip(v.data, self.data)])

    def norm(self):
        return sqrt(self.dot(self))

    def max(self):
        return max(self.data)

    def min(self):
        return min(self.data)

    def __iter__(self):
        for x in self.data:
            yield x

    def __setitem__(self, i, v):
        assert isinstance(i, int)
        self.data[i] = self.dtype(v)

    def __getitem__(self, i):
        i = _mkslice(i, self.dim)
        res = self.data[i]
        if len(res) == 1:
            return res[0]
        return array(res)

    def transpose(self):
        return matrix(self.data, self.dtype)

    def __float__(self):
        assert self.dim == 1
        return float(self.data[0])


# ------------------------------------------------------------------------------
# minimalistic numpy.matrix replacement class used as fallback
# when numpy is not found in geometry module
class matrix(object):
    def __init__(self, data, dtype=None, copy=True, transpose=False):
        # check input data types:
        types = set([type(v) for v in data])
        if len(types) > 1:
            raise TypeError
        t = types.pop()
        # import data:
        if t in constants:
            self.data = [array(data, dtype, copy)]
        else:
            if transpose:
                data = zip(*data)
            self.data = [array(v, dtype, copy) for v in data]
        # define matrix sizes:
        self.n = len(self.data)
        sizes = set([len(v) for v in self.data])
        if len(sizes) > 1:
            raise ValueError
        self.p = sizes.pop()
        if dtype is None:
            # coerce types of arrays of matrix:
            types = set([v.dtype for v in self.data])
            tc, dtype = coerce_(types)
            for v in self.data:
                v.coerce(dtype)
        self.dtype = dtype

    def __len__(self):
        return self.n * self.p

    def __str__(self):
        s = "\n ".join([str(v) for v in self.data])
        return "[%s]" % s.strip()

    @property
    def shape(self):
        return (self.n, self.p)

    def lvecs(self):
        return self.data

    def cvecs(self):
        return [array(v, self.dtype) for v in zip(*self.data)]

    def copy(self):
        return matrix(self.data, self.dtype)

    def transpose(self):
        return matrix(self.data, dtype=self.dtype, transpose=True)

    def sum(self):
        return sum([sum(v) for v in self.data])

    @make_ij_slices
    def __getitem__(self, ij):
        I, J = ij
        l = self.lvecs()[I]
        m = matrix([v[J] for v in l])
        if m.n == 1:
            v = m.data[0]
            if v.dim == 1:
                return v[0]
            if len(l) > 1:
                return v
        return m

    @make_ij_slices
    def __setitem__(self, ij, v):
        I, J = ij
        Ri = range(I.start, I.stop, I.step)
        Rj = range(J.start, J.stop, J.step)
        if type(v) in constants:
            v = (v,)
        value = (x for x in v)
        for i in Ri:
            for j in Rj:
                self.data[i][j] = next(value)

    def __add__(self, m):
        if isinstance(m, constants):
            return matrix([u + m for u in self.data])
        else:
            assert self.shape == m.shape
            return matrix([u + v for (u, v) in zip(self.data, m.data)])

    def __sub__(self, m):
        if isinstance(m, constants):
            return matrix([u - m for u in self.data])
        else:
            assert self.shape == m.shape
            return matrix([u - v for (u, v) in zip(self.data, m.data)])

    def __neg__(self):
        return matrix([-x for x in self.data], dtype=self.dtype)

    def __float__(self):
        assert self.n == 1 and self.p == 1
        return self[0, 0]

    def __radd__(self, v):
        return self + v

    def __rsub__(self, v):
        return (-self) + v

    def __rmul__(self, k):
        if not isinstance(k, constants):
            raise TypeError
        return matrix([k * v for v in self.data])

    def __mul__(self, X):
        if isinstance(X, constants):
            return X * self
        if isinstance(X, array):
            assert X.dim == self.p
            return array([v.dot(X) for v in self.data])
        if isinstance(X, matrix):
            assert X.n == self.p
            return matrix([self * v for v in X.cvecs()])

    def __pow__(self, v):
        S = [self] * v
        assert len(S) > 0
        return reduce(lambda x, y: x * y, S)

    def __iter__(self):
        for l in self.data:
            for v in l:
                yield v


class SimplexMin(object):
    def __init__(self, A, b, c):
        self.A = A
        self.b = b
        self.c = c
        self.tableau()

    def tableau(self):
        self.T = []

    def setup(self):
        self.enter = []
        delf.outer = []