simplex.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Симплекс-метод (максимизация) для задач вида:
  max Z = c^T x
  A x <= b,  x >= 0

Печатает:
- приведение к канонической форме
- пошаговые симплекс-таблицы как на тетрадных листах
- выбор ведущего столбца и строки, отношения θ
- итоговый оптимум

Вычисления ведутся в поле рациональных чисел (fractions.Fraction),
чтобы красиво выводились дроби как 1/3, 2/7 и т.п.

Как использовать:
1) Измените блок INPUT_DATA ниже.
2) Запустите:  python3 simplex.py
"""
from __future__ import annotations
from fractions import Fraction as F
from typing import List, Tuple

# ======== INPUT_DATA ========
# Пример из задания:
# F = x1 - x2 -> max
# -2x1 + x2 <= 2
#  x1 - 2x2 <= 2
#  x1 +  x2 <= 5
# x1, x2 >= 0
c = [1, -1]  # коэффициенты целевой функции (x1, x2)
A = [
    [-2, 1],
    [1, -2],
    [1, 1],
]
b = [2, 2, 5]
# ============================


def to_frac_matrix(A: List[List[float]]) -> List[List[F]]:
    return [[F(v) for v in row] for row in A]


def to_frac_vec(v: List[float]) -> List[F]:
    return [F(x) for x in v]


def fmt(x: F) -> str:
    if x == 0:
        return "0"
    if x.denominator == 1:
        return str(x.numerator)
    return f"{x.numerator}/{x.denominator}"


def fmt_term(coeff: F, name: str, first: bool) -> str:
    # формат одного члена: знаки и коэффициенты, 1 не пишем
    if coeff == 0:
        return ""
    sign = "-" if coeff < 0 else "+"
    abs_coeff = -coeff if coeff < 0 else coeff
    if abs_coeff == 1:
        piece = name
    else:
        piece = f"{fmt(abs_coeff)}{name}"
    if first:
        return ("-" if coeff < 0 else "") + piece
    else:
        return f" {sign} {piece}"


def format_lincomb(coeffs: List[F], names: List[str]) -> str:
    # собирает строку вида: -2x1 + x2 - 3/2 x3
    out = []
    first = True
    for c, nm in zip(coeffs, names):
        term = fmt_term(c, nm, first)
        if term:
            out.append(term)
            first = False
    if not out:
        return "0"
    return "".join(out)


def draw_tableau(basis_names: List[str], var_names: List[str], A: List[List[F]], b: List[F],
                 delta: List[F], Z_val: F, step_title: str) -> str:
    # рамки
    LEFT = 6
    W = 7
    lines: List[str] = []
    lines.append(step_title)
    def hline():
        return "+" + "-"*LEFT + "+" + "+".join(["-"*W for _ in range(len(var_names)+1)]) + "+"
    # заголовок с переменными и свободным членом b
    lines.append(hline())
    header = f"{'Базис':<{LEFT}}|" + "|".join(f"{name:^{W}}" for name in var_names) + f"|{'b':^{W}}|"
    lines.append(header)
    lines.append(hline())
    # строки ограничений (базис)
    for r, row in enumerate(A):
        line = f"{basis_names[r]:<{LEFT}}|" + "|".join(f"{fmt(v):^{W}}" for v in row) + f"|{fmt(b[r]):^{W}}|"
        lines.append(line)
        lines.append(hline())
    # строка дельт (используем соглашение Zj - Cj)
    dline = f"{'Δ':<{LEFT}}|" + "|".join(f"{fmt(v):^{W}}" for v in delta) + f"|{fmt(Z_val):^{W}}|"
    lines.append(dline)
    lines.append(hline())
    return "\n".join(lines)


def simplex_max(A_in: List[List[float]], b_in: List[float], c_in: List[float]) -> None:
    # к рациональным
    A0 = to_frac_matrix(A_in)
    b0 = to_frac_vec(b_in)
    c0 = to_frac_vec(c_in)

    m, n = len(A0), len(A0[0])

    # проверка знаков правых частей
    for i in range(m):
        if b0[i] < 0:
            # умножим строку на -1, меняя знак неравенства (здесь все неравенства <=, так что допустимо)
            A0[i] = [ -aij for aij in A0[i] ]
            b0[i] = -b0[i]

    # приведение к канонической форме: добавляем m свободных переменных s_i (обозначим их как x_{n+1}..x_{n+m})
    I = [[F(1) if i==j else F(0) for j in range(m)] for i in range(m)]
    A = [ A0[i] + I[i] for i in range(m) ]  # m × (n+m)
    c = c0 + [F(0)]*m                       # (n+m)
    var_names = [f"x{i+1}" for i in range(n+m)]
    basis = [n+i for i in range(m)]
    basis_names = [var_names[idx] for idx in basis]

    # Печать исходной задачи и затем канонической формы как на фото
    print("Исходная задача:")
    # F = ... -> max
    expr_orig = format_lincomb(c0, [f"x{j+1}" for j in range(n)])
    print(f"F = {expr_orig} -> max")
    for i in range(m):
      left = format_lincomb(A0[i], [f"x{j+1}" for j in range(n)])
      print(f"{left} <= {fmt(b0[i])}")
    print("x_j >= 0")
    print()
    print("=> Каноническая форма (табличный вид):")
    # F - c^T x = 0
    expr_obj = format_lincomb([-c0[j] for j in range(n)], [f"x{j+1}" for j in range(n)])
    print(f"F {(' ' + expr_obj) if expr_obj and not expr_obj.startswith('-') else expr_obj} = 0")
    for i in range(m):
        left = format_lincomb(A0[i], [f"x{j+1}" for j in range(n)]) + f" + x{n+i+1}"
        print(f"{left} = {fmt(b0[i])}")
    print(f"x_i >= 0; i = 1..{n+m}")
    print()

    # вспомогательные функции
    def compute_delta_and_Z(A: List[List[F]], b: List[F], basis: List[int]) -> Tuple[List[F], F, List[F]]:
        # Zj = sum(c_Bi * a_ij)
        cB = [c[idx] for idx in basis]
        Zj = [F(0) for _ in range(len(c))]
        for j in range(len(c)):
            Zj[j] = sum(cB[i]*A[i][j] for i in range(m))
        # Используем соглашение: Δ = Zj - Cj (как на ваших фото)
        delta = [Zj[j] - c[j] for j in range(len(c))]
        Z_val = sum(cB[i]*b[i] for i in range(m))
        return delta, Z_val, Zj

    step = 1
    while True:
        delta, Z_val, Zj = compute_delta_and_Z(A, b0, basis)
        print(draw_tableau(basis_names, var_names, A, b0, delta, Z_val, f"Симплекс-таблица — шаг {step}"))

        # выбор ведущего столбца (entering)
        min_delta = min(delta)
        if min_delta >= 0:
            print("Оптимальность достигнута (все Δ >= 0).")
            break
        j_enter = min(range(len(delta)), key=lambda j: delta[j])
        print(f"Ведущий столбец: {var_names[j_enter]} (Δ = {fmt(delta[j_enter])} = min)")

        # отношения θ = b_i / a_ij (печатаем для всех строк; выбираем минимум среди положительных a_ij)
        ratios_all_lines = []
        ratios_candidates = []
        for i in range(m):
            aij = A[i][j_enter]
            if aij == 0:
                ratios_all_lines.append((basis_names[i], None, aij))
            else:
                val = b0[i]/aij
                ratios_all_lines.append((basis_names[i], val, aij))
                if aij > 0:
                    ratios_candidates.append((val, i))
        if not ratios_candidates:
            print("Задача неограничена: в ведущем столбце нет положительных элементов.")
            return
        theta, i_leave = min(ratios_candidates)
        # печать как у вашей подруги: x_k = b/a = значение, с пометкой -> min у выбранной строки
        for name, val, denom in ratios_all_lines:
            if denom == 0:
                print(f"{name} = —")
            else:
                tail = ""
                if val == theta and denom > 0 and name == basis_names[i_leave]:
                    tail = " -> min"
                print(f"{name} = {fmt(b0[basis_names.index(name)])}/{fmt(denom)} = {fmt(val)}{tail}")
        print(f"Ведущая строка: {basis_names[i_leave]} (минимум θ = {fmt(theta)})")
        print(f"Опорный элемент: A[{i_leave+1},{j_enter+1}] = {fmt(A[i_leave][j_enter])}")

        # поворот (Гаусс-Жордан)
        pivot = A[i_leave][j_enter]
        # нормируем ведущую строку
        A[i_leave] = [v/pivot for v in A[i_leave]]
        b0[i_leave] = b0[i_leave]/pivot
        # зануляем столбец в остальных строках
        for i in range(m):
            if i == i_leave:
                continue
            factor = A[i][j_enter]
            if factor != 0:
                A[i] = [A[i][k] - factor*A[i_leave][k] for k in range(len(c))]
                b0[i] = b0[i] - factor*b0[i_leave]
        # обновляем базис
        basis[i_leave] = j_enter
        basis_names[i_leave] = var_names[j_enter]
        print()
        step += 1

    # вывод оптимума
    x_full = [F(0) for _ in range(len(c))]
    for r, idx in enumerate(basis):
        x_full[idx] = b0[r]
    x = x_full[:n]
    #print("Итоговая таблица/план:")
    #print(draw_tableau(basis_names, var_names, A, b0, delta, Z_val, "Финальная симплекс-таблица"))
    print("Оптимум:")
    for j in range(n):
        print(f"x{j+1} = {fmt(x[j])}")
    print(f"Z* = {fmt(Z_val)}")


if __name__ == "__main__":
    simplex_max(A, b, c)