Untitled

import sympy as sp
import pandas as pd

def M_method_symbolic(A, b, c, maximize=True, M_large_value=10**6):
    """
    M-метод (символічно з урахуванням M як "дуже великого" позитивного числа).
    Параметри:
      A, b, c - стандартні матриці/вектори (list або sympy-compatible)
      maximize - True для максимізації, False для мінімізації
      M_large_value - числова підстановка для порівнянь (fallback)
    Повертає (x_vector, objective_value) або (None, None) при необмеженості.
    """
    A = sp.Matrix(A)
    b = sp.Matrix(b)
    c = list(c)  # зробимо список, будемо дописувати штучні коефіцієнти
    m, n = A.shape
    M = sp.Symbol('M', real=True, positive=True)

    # === 1) Формуємо розширену матрицю A_aug і початковий базис B ===
    A_aug = A.copy()
    artificial_indices = []   # індекси доданих штучних змінних
    B = []                    # індекси базисних змінних, повинно вийти len(B) == m

    for i in range(m):
        unit_vector = [1 if k == i else 0 for k in range(m)]
        found = None
        for j in range(n):
            if list(A[:, j]) == unit_vector:
                found = j
                break
        if found is not None:
            B.append(found)
        else:
            # додаємо штучну змінну (новий стовпець)
            new_col = [0]*m
            new_col[i] = 1
            A_aug = A_aug.col_insert(A_aug.cols, sp.Matrix(new_col))
            new_idx = A_aug.cols - 1
            artificial_indices.append(new_idx)
            B.append(new_idx)
            c.append(-M)  # штучна змінна має коеф -M

    # якщо вхідний вектор c мав довжину n, ми вже дописали -M для штучних
    n_new = A_aug.shape[1]
    # захист: якщо c має менше елементів (наприклад, не дорівнює n), привести до n_new
    if len(c) < n_new:
        c = list(c) + [0]*(n_new - len(c))
    c_new = sp.Matrix(c)

    # Побудуємо початковий розв'язок x0: значення базисних змінних = b (за рядками)
    x = sp.Matrix([0]*n_new)
    for row_i in range(m):
        basis_var_idx = B[row_i]  # базисна змінна, що відповідає цьому рядку
        x[basis_var_idx] = b[row_i]

    N = [i for i in range(n_new) if i not in B]

    print(f"Початкова база: {[f'x{i+1}' for i in B]}")
    print(f"Початкова точка x0 = {x.T}")

    # === Допоміжні функції для визначення знаку з урахуванням M ===
    def sign_with_M(expr):
        """Повертає 1 якщо позитивне, -1 якщо негативне, 0 якщо нуль або невизначено (намагатись чисельно)."""
        expr = sp.simplify(expr)
        # якщо присутній M, дивимось коефіцієнт при M
        if expr.has(M):
            # коефіцієнт при M (припускаємо лінійність по M)
            coeff = sp.expand(expr).coeff(M)
            if coeff.is_positive is True:
                return 1
            if coeff.is_negative is True:
                return -1
            # невизначено за символікою — використовуємо числову підстановку
            try:
                val = float(sp.N(expr.subs(M, M_large_value)))
                return 1 if val > 0 else (-1 if val < 0 else 0)
            except Exception:
                return 0
        # якщо M не присутній
        if expr.is_number:
            try:
                val = float(sp.N(expr))
                return 1 if val > 0 else (-1 if val < 0 else 0)
            except Exception:
                return 0
        if expr.is_positive is True:
            return 1
        if expr.is_negative is True:
            return -1
        # числова підстановка як останній fallback
        try:
            val = float(sp.N(expr))
            return 1 if val > 0 else (-1 if val < 0 else 0)
        except Exception:
            return 0

    def is_positive(expr):
        return sign_with_M(expr) == 1

    def is_negative(expr):
        return sign_with_M(expr) == -1

    step = 0
    # === Основний цикл ===
    while True:
        step += 1

        A_B = A_aug[:, B]   # m x m (за умовою B має довжину m)
        A_N = A_aug[:, N]   # m x (n_new - m)

        # c_B та c_N як вектори
        c_B = sp.Matrix([c_new[i] for i in B])
        c_N = sp.Matrix([c_new[i] for i in N])

        # Обчислення потенціалів π: A_B^T * π = c_B
        try:
            pi = A_B.T.LUsolve(c_B)
        except Exception:
            # fallback на псевдообернену
            pi = A_B.T.pinv() * c_B

        # Δ = c_N^T - π^T * A_N
        z_N = pi.T * A_N
        delta = (c_N.T - z_N)  # 1 x len(N)

        # β = A_B^{-1} * b
        try:
            beta = A_B.LUsolve(b)
        except Exception:
            beta = A_B.pinv() * b

        z_value = (c_B.T * beta)[0]

        print("\n" + "="*70)
        print(f"🔹 КРОК {step}")
        print("="*70)
        print(f"База: {[f'x{i+1}' for i in B]}")
        print(f"π = {sp.simplify(pi.T)}")
        print(f"Δ = {sp.simplify(delta)}")
        print(f"β = {beta.T}")
        print(f"f(x) = {sp.simplify(z_value)}")

        # === Формуємо таблицю (A_B^{-1} * A_aug) у вигляді list-of-lists ===
        try:
            table_mat = A_B.inv() * A_aug
        except Exception:
            table_mat = A_B.pinv() * A_aug
        table = table_mat.tolist()  # m x n_new як звичайний nested list

        # формуємо рядки DataFrame: [β] + ряд таблиці
        df_rows = []
        for i_row in range(m):
            # безпечний доступ: якщо table має менше рядків — це помилка (але має бути m)
            if i_row >= len(table):
                raise IndexError("internal error: table has fewer rows than m (це не повинно відбуватися)")
            df_rows.append([beta[i_row]] + table[i_row])

        df = pd.DataFrame(
            df_rows,
            columns=["β"] + [f"x{i+1}" for i in range(n_new)],
            index=[f"x{i+1}" for i in B]
        )
        print("\nПоточна таблиця:")
        print(df)

        # === Умова оптимальності ===
        delta_list = [delta[0, j] for j in range(delta.shape[1])] if delta.shape[1] > 0 else []
        if maximize:
            enter_candidates = [j for j, val in enumerate(delta_list) if is_positive(val)]
        else:
            enter_candidates = [j for j, val in enumerate(delta_list) if is_negative(val)]

        if not enter_candidates:
            print("\n✅ Оптимум досягнуто.")
            print(f"x* = {x.T}")
            print(f"f(x*) = {sp.simplify(z_value)}")
            return x, sp.simplify(z_value)

        # вибираємо перший кандидат на вхід (можна змінити політику)
        j_in = enter_candidates[0]
        entering = N[j_in]  # індекс змінної в повному наборі

        # обчислення d: A_B * d = A_aug[:, entering]  => d = A_B^{-1} * column
        try:
            d = A_B.LUsolve(A_aug[:, entering])
        except Exception:
            d = A_B.pinv() * A_aug[:, entering]

        # перевірка необмеженості: якщо немає додатних елементів у d -> необмежений
        if not any(sign_with_M(d[i]) == 1 for i in range(m)):
            print("⚠️ Розв’язок не обмежений.")
            return None, None

        # обчислення θ: ratios тільки для d[i] > 0 (враховуємо символіку через sign_with_M)
        ratios = [sp.oo] * m
        ratios_eval = [float('inf')] * m
        for i_row in range(m):
            if sign_with_M(d[i_row]) == 1:
                r = sp.simplify(beta[i_row] / d[i_row])
                ratios[i_row] = r
                try:
                    # числова апроксимація для порівняння — підставляємо велике M
                    r_eval = float(sp.N(r.subs(M, M_large_value)))
                except Exception:
                    r_eval = float('inf')
                ratios_eval[i_row] = r_eval

        i_out = int(min(range(m), key=lambda i: ratios_eval[i]))
        leaving = B[i_out]

        print(f"\n➡️ Входить змінна x{entering+1}, виходить x{leaving+1}")
        print(f"d = {d.T}")
        print(f"θ (символічні) = {ratios}")
        print(f"θ (числ. апрокс) = {ratios_eval}")
        print(f"Мінімальне θ = {ratios[i_out]} (рядок {i_out})")

        # Оновлення вектора x
        theta = ratios[i_out]
        for i_row in range(m):
            x[B[i_row]] = sp.simplify(beta[i_row] - theta * d[i_row])
        x[entering] = sp.simplify(theta)
        x[leaving] = 0

        # Оновлюємо базис
        B[i_out] = entering
        N = [i for i in range(n_new) if i not in B]

        print(f"🔁 Новий базис: {[f'x{i+1}' for i in B]}")
        print(f"x = {x.T}")
        print(f"f(x) = {sp.simplify(c_new.dot(x))}")

# === ПРИКЛАД ===
if __name__ == "__main__":
    A = [
        [3, 1, 1, 1, 1],
        [2, -1, 3, 0, 0],
        [0, 5, 6, 1, 0]
    ]
    b = [5, 4, 11]
    c = [5, -1, -1, 0, 0]

    print("=== АНАЛІТИЧНИЙ M-МЕТОД ===")
    x_star, f_star = M_method_symbolic(A, b, c, maximize=False)
    print("\nРЕЗУЛЬТАТ:")
    print("x* =", x_star)
    print("f(x*) =", f_star)