boikwd

from numpy import array, log10, ones, eye
from numpy.linalg import solve
import functools
import operator

# data
A = array([
    [1, 2 / 3, 2, 5 / 2, 5 / 3, 5], [
        3 / 2, 1, 3, 10 / 3, 3, 9], [
        1 / 2, 1 / 3, 1, 4 / 3, 7 / 8, 5 / 2], [
        2 / 5, 3 / 10, 3 / 4, 1, 5 / 6, 12 / 5], [
        3 / 5, 1 / 3, 8 / 7, 6 / 5, 1, 3], [
        1 / 5, 1 / 9, 2 / 5, 5 / 12, 1 / 3, 1]])
vA = array([3, 1])
positionsA = [5, 6]
B = array([
    [1, 2 / 5, 3, 7 / 3, 1 / 2, 1], [
        5 / 2, 1, 4 / 7, 5 / 8, 1 / 3, 3], [
        1 / 3, 7 / 4, 1, 1 / 2, 2, 1 / 2], [
        3 / 7, 8 / 5, 2, 1, 4, 2], [
        2, 3, 1 / 2, 1 / 4, 1, 1 / 2], [
        1, 1 / 3, 2, 1 / 2, 2, 1]])
vB = array([2, 1 / 2, 1])
positionsB = [4, 5, 6]
C = array([
    [1, 17 / 4, 17 / 20, 8 / 5, 23 / 6, 8 / 3], [
        4 / 17, 1, 1 / 5, 2 / 5, 9 / 10, 2 / 3], [
        20 / 17, 5, 1, 21 / 10, 51 / 10, 10 / 3], [
        5 / 8, 5 / 2, 10 / 21, 1, 5 / 2, 11 / 6], [
        6 / 23, 10 / 9, 10 / 51, 2 / 5, 1, 19 / 30], [
        3 / 8, 3 / 2, 3 / 10, 6 / 11, 30 / 19, 1]])
vC = array([2, 5])
positionsC = [2, 4]


def get_koczkodaj_inconsistency_index(matrix: array):
    max_val = 0
    n = len(matrix)
    for i in range(n):
        for j in range(n):
            for k in range(n):
                min_val = min(abs(1 - matrix[i][k] * matrix[k][j] / matrix[i][j]),
                              abs(1 - matrix[i][j] / (matrix[i][k] * matrix[k][j])))
                max_val = max(max_val, min_val)
    return max_val


def create_A_B_matrices(matrix, known_values, positions):
    n = len(matrix)
    k = len(known_values)

    # check if swap is needed
    swap_wrong_columns_and_rows(matrix, n, positions)

    A = matrix[0:(n - k), 0:(n - k)]
    B = matrix[0:(n - k), (n - k):n]

    for i in range(n - k):
        for j in range(n - k):
            if i != j:
                A[i][j] = A[i][j] * (-1 / (n - 1))

    B = array(B).dot(known_values)
    B = B * (1 / (n - 1))
    return [A, B]


def swap_wrong_columns_and_rows(matrix, n, positions):
    do_swap = False
    index = n
    k = len(positions)
    for i in reversed(positions):
        if index != i:
            do_swap = True
            break
        else:
            index = index - 1
    if do_swap:
        # mix matrice so vectors 'from' position are last in matrix
        index = n - 1
        for position in reversed(positions):
            matrix[:, [position-1, index]] = matrix[:, [index, position-1]]
            index = index - 1

        index = n - 1
        for position in reversed(positions):
            matrix[[position - 1, index], :] = matrix[[index, position - 1], :]
            index = index - 1


def insert_on_correct_positions(known_values, positions, w):
    for i in range(len(positions)):
        w.insert(positions[i] - 1, known_values[i])


def HRE_arithmetic(matrix: array, known_values: array, positions: list):
    matrix = matrix.copy()
    koczkodaj_index = get_koczkodaj_inconsistency_index(matrix)
    n = len(matrix)
    k = len(known_values)
    has_solution = koczkodaj_index < 1 - (1 + (4 * (n - 1) * (n - k - 2)) ** 0.5) / (2 * (n - 1))

    [A, B] = create_A_B_matrices(matrix, known_values, positions)
    solution = solve(A, B).tolist()

    insert_on_correct_positions(known_values, positions, solution)

    for s in solution:
        if s < 0:
            has_solution = False
            break
    return [has_solution, solution]


def HRE_geometric(matrix: array, known_values: array, positions: list):
    def foldl(func, acc, xs):
        return functools.reduce(func, xs, acc)

    matrix = matrix.copy()
    n = len(matrix)
    swap_wrong_columns_and_rows(matrix, n, positions)

    k = len(positions)
    b_ = []
    multiplied_know_values = foldl(operator.mul, 1, known_values)
    for i in range(n - k):
        tmp = foldl(operator.mul, multiplied_know_values, matrix[i, :])
        b_.append(log10(tmp))

    A_ = ones((n - k, n - k)) * (-1) + eye(n - k) * n

    W_ = solve(A_, b_)
    solution = [10 ** x for x in W_]

    insert_on_correct_positions(known_values, positions, solution)
    return solution


def foldl(func, acc, xs):
    return functools.reduce(func, xs, acc)


# # A
# [has_solution, solution] = HRE_arithmetic(A, vA, positionsA)
# s = "It probably doesn't have solution -> "
# if has_solution:
#     s = "It has solution -> "
# print("Arithmetic HRE for matrix A.", s, solution)
# print("Geometric HRE for matrix A", HRE_geometric(A, vA, positionsA))
#
# # B
# [has_solution, solution] = HRE_arithmetic(B, vB, positionsB)
# s = "It probably doesn't have solution -> "
# if has_solution:
#     s = "It has solution -> "
# print("Arithmetic HRE for matrix B", s, solution)
# print("Geometric HRE for matrix B", HRE_geometric(B, vB, positionsB))
#
# # C
# [has_solution, solution] = HRE_arithmetic(C, vC, positionsC)
# s = "It probably doesn't have solution -> "
# if has_solution:
#     s = "It has solution -> "
# print("Arithmetic HRE for matrix C", s, solution)
print("Geometric HRE for matrix C", HRE_geometric(C, vC, positionsC))