Gauss Algorithm

gauss.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from __future__ import print_function, division

import re
import sys
import numpy as np


class GaussException(Exception):
    def __init__(self, message):
        super(GaussException, self).__init__(message)


def collect_row(row):
    pattern = re.compile(r'[+-]?\d+')
    columns = []
    match = None
    while True:
        index = match.end() if match else 0
        match = pattern.search(row, index)
        if not match:
            break
        columns.append(float(match.group()))

    return columns


def parse(string):
    pattern = re.compile(r'(?:\s*[+-]?\d+\s*)+')
    matrix = []
    match = None
    while True:
        index = match.end() if match else 0
        match = pattern.search(string, index)
        if not match:
            break
        matrix.append(collect_row(match.group()))

    return np.array(matrix)


def first_not_null(row):
    booleans = (row != 0)
    indices = np.argwhere(booleans)

    return indices[0] if len(indices) > 0 else len(row)


def next_row(matrix, row):
    choice, index = row, first_not_null(matrix[row])
    for i in range(row + 1, len(matrix)):
        first = first_not_null(matrix[i])
        if first < index:
            choice, index = i, row

    return choice, index


def check(row, column):
    if column == len(row):
        print("Infinite solutions!")
        return False
    elif column == len(row) - 1:
        print("No solutions!")
        return False

    return True


def swap(matrix, first, second):
    matrix[first], matrix[second] = matrix[second], matrix[first]


def not_strict(matrix):
    for row in range(len(matrix)):
        choice, column = next_row(matrix, row)
        if not check(matrix[choice], column):
            return []
        swap(matrix, choice, row)
        matrix[row] /= matrix[row][column]
        for other_row in matrix[row + 1:]:
            other_row -= other_row[column] * matrix[row]

    return matrix


def strict(matrix):
    reverse = matrix[::-1]
    for i, row in enumerate(reverse[:-1]):
        index = -(i + 2)
        assert row[index] == 1
        for other_row in reverse[i + 1:]:
            other_row -= other_row[index] * row

    return matrix


def gauss(string):
    matrix = parse(string)
    if matrix.shape[0] < matrix.shape[1] - 1:
        raise GaussException("Too few equations to find a solution!")

    matrix = not_strict(matrix)
    if len(matrix) == 0:
        return
    matrix = strict(matrix)

    print("Solutions: {0}".format(matrix[:, -1]))


def main():
    if len(sys.argv) != 2:
        raise RuntimeError("Usage: gauss.py <a b c d, e f g h, ...>")

    gauss(sys.argv[1])


if __name__ == "__main__":
    main()