maxRevMatrixEl.py

import itertools
import numpy as np
from numpy.linalg import inv
from time import time
from sympy.matrices import Matrix
import sympy

N = 6
# Все ненулевые строки
zero_one_rows = [np.matrix(row, dtype=np.float64) for row in itertools.product([0, 1], repeat=N) if any(row)]
# Текущие лучшие значения
max_el = float('-inf')
max_row = float('-inf')
max_mat = float('-inf')
# Последовательность улучшений
el_seq = []
row_seq = []
mat_seq = []
t = 0
st = time()
for row in itertools.combinations(zero_one_rows, 6):
    t += 1
    if t % 100000 == 0:
        print(f'{100 * t / 67945521:0.3f}%, {time() - st:0.3f}с со старта')
        print(f'макс. элемент {max_el}, макс. сумма в строке {max_row}, макс.сумма всех {max_mat}')
        print(len(el_seq), len(row_seq), len(mat_seq))
    m = np.vstack(row)
    try:
        minv = inv(m)
    except np.linalg.LinAlgError:
        continue
    m_el = minv.max()
    m_row = minv.sum(axis=1).max()
    m_mat = minv.sum()
    # Если матрица «интересная», то считаем точно
    if m_el > max_el or m_row > max_row or m_mat > max_mat:
        M = Matrix(m.astype(np.int))
        try:
            Minv = M.inv()
        except sympy.matrices.common.NonInvertibleMatrixError:
            # Оказалось, что если посчитать точно, то матрица таки необратима
            continue
        minv = np.array(Minv)  # Теперь здесь точные значения
        m_el = minv.max()
        m_row = minv.sum(axis=1).max()
        m_mat = minv.sum()
    # Сохраняем
    if m_el > max_el:
        max_el = m_el
        el_seq.append(m)
    if m_row > max_row:
        max_row = m_row
        row_seq.append(m)
    if m_mat > max_mat:
        max_mat = m_mat
        mat_seq.append(m)

print('total matrixes', t)
print(f'макс. элемент {max_el}, макс. сумма в строке {max_row}, макс.сумма всех {max_mat}')
print(len(el_seq), len(row_seq), len(mat_seq))
for seq in el_seq, row_seq, mat_seq:
    print('Интересные матрицы:')
    for m in el_seq[-10]:
        minv = np.array(Matrix(m.astype(np.int)).inv())
        print(m)
        print(minv)
        print(minv.max(), minv.sum(axis=1).max(), minv.sum())
        print('*' * 50)