#!/usr/bin/env catnip
# Planification de production par programmation linéaire en nombres entiers.
#
# Trois produits consomment des heures machine, du travail et de la matière. Le
# solveur choisit des quantités entières et peut activer un lot d'heures
# supplémentaires, représenté par une variable binaire. Catnip garde le modèle
# métier lisible (produits, scénarios, statuts et rapports) ; SciPy construit et
# résout le MILP. Une planche compare les décisions sous trois contraintes.
#
# Contrôles : intégralité, capacités et objectif sont recalculés indépendamment
# depuis la solution entière renvoyée par HiGHS.
#
# DEPS: numpy scipy matplotlib

numpy = import('numpy')
import('scipy.optimize', 'milp', 'Bounds', 'LinearConstraint')
mpl = import('matplotlib')
mpl.use('Agg')
plt = import('matplotlib.pyplot')
import('pathlib', 'Path')

script_dir = Path(META.file).parent
output_dir = script_dir / 'output'
output_dir.mkdir(exist_ok=True)

struct Product {
    name: str; profit: float; machine: float; labor: float; material: float;
    minimum: int; maximum: int;
}

struct Scenario {
    name: str; machine: float; labor: float; material: float; overtime_cost: float;

    machine_capacity(self, overtime: bool): float => {
        extra = if overtime { 18.0 } else { 0.0 }
        self.machine + extra
    }

    labor_capacity(self, overtime: bool): float => {
        extra = if overtime { 12.0 } else { 0.0 }
        self.labor + extra
    }
}

union SolveStatus {
    optimal; infeasible; failed

    label(self): str => {
        match self {
            SolveStatus.optimal    => { "optimal" }
            SolveStatus.infeasible => { "infaisable" }
            SolveStatus.failed     => { "échec solveur" }
        }
    }
}

struct Plan {
    scenario: Scenario; units: list[int]; overtime: bool; net_profit: float;
    machine_used: float; labor_used: float; material_used: float;
    status: SolveStatus; integrality_error: float;
    capacity_violation: float; objective_error: float;

    display(self, products: list[Product]) => {
        parts = list()
        i = 0
        while i < len(products) {
            parts.append(f"{products[i].name}={self.units[i]}")
            i = i + 1
        }
        print(f"  {self.scenario.name:<14} → {self.status.label()}, {', '.join(parts)}")
        print(f"    heures sup.={self.overtime}  marge nette={round(self.net_profit, 2)} €")
        print(
            f"    machine={self.machine_used}/{self.scenario.machine_capacity(self.overtime)} h  " +
                f"travail={self.labor_used}/{self.scenario.labor_capacity(self.overtime)} h  " +
                f"matière={self.material_used}/{self.scenario.material} kg"
        )
        print(
            f"    contrôles : intégralité={self.integrality_error}, " +
                f"violation={self.capacity_violation}, objectif={self.objective_error}"
        )
    }
}

products = list(
    Product('capteurs', 38.0, 2.0, 1.0, 3.0, 6, 45),
    Product('contrôleurs', 52.0, 1.0, 3.0, 4.0, 4, 40),
    Product('passerelles', 71.0, 3.0, 2.0, 5.0, 2, 30),
)

scenarios = list(
    Scenario('semaine normale', 72.0, 76.0, 120.0, 180.0),
    Scenario('commande urgente', 52.0, 100.0, 180.0, 230.0),
    Scenario('matière rare', 76.0, 82.0, 88.0, 150.0),
)

status_for = (result): SolveStatus => {
    if result.success { SolveStatus.optimal }
    elif int(result.status) == 2 { SolveStatus.infeasible }
    else { SolveStatus.failed }
}

solve = (scenario: Scenario): Plan => {
    n_products = len(products)

    # SciPy minimise : les marges sont donc niées. La dernière variable est
    # binaire et facture le lot d'heures supplémentaires dans l'objectif.
    objective = products.[(p) => { -p.profit }] + list(scenario.overtime_cost)
    integrality = numpy.ones(n_products + 1, dtype='int32')
    lower = products.[(p) => { p.minimum }] + list(0)
    upper = products.[(p) => { p.maximum }] + list(1)

    # Les coefficients négatifs de la variable binaire ajoutent 18 h machine
    # et 12 h de travail quand le lot supplémentaire est activé.
    matrix = numpy.array(list(
            products.[(p) => { p.machine }] + list(-18.0),
            products.[(p) => { p.labor }] + list(-12.0),
            products.[(p) => { p.material }] + list(0.0),
        ),
        dtype='float64')
    capacities = numpy.array(list(scenario.machine, scenario.labor, scenario.material))

    result = milp(
        c=numpy.array(objective, dtype='float64'),
        integrality=integrality,
        bounds=Bounds(numpy.array(lower), numpy.array(upper)),
        constraints=LinearConstraint(matrix, numpy.full(3, -numpy.inf), capacities),
        options=dict(disp=False),
    )

    status = status_for(result)
    if not result.success {
        Plan(scenario, list(0, 0, 0), False, 0.0, 0.0, 0.0, 0.0, status, 0.0, 0.0, 0.0)
    } else {
        rounded = numpy.rint(result.x).astype('int64')
        units = rounded[:n_products].tolist()
        overtime = int(rounded[n_products]) == 1
        unit_array = numpy.array(units, dtype='float64')

        margins = numpy.array(products.[(p) => { p.profit }])
        machine_coeff = numpy.array(products.[(p) => { p.machine }])
        labor_coeff = numpy.array(products.[(p) => { p.labor }])
        material_coeff = numpy.array(products.[(p) => { p.material }])

        machine_used = float(numpy.dot(machine_coeff, unit_array))
        labor_used = float(numpy.dot(labor_coeff, unit_array))
        material_used = float(numpy.dot(material_coeff, unit_array))
        overtime_charge = if overtime { scenario.overtime_cost } else { 0.0 }
        net_profit = float(numpy.dot(margins, unit_array)) - overtime_charge

        effective = numpy.array(list(
                scenario.machine_capacity(overtime),
                scenario.labor_capacity(overtime),
                scenario.material,
            ))
        used = numpy.array(list(machine_used, labor_used, material_used))
        violation = float(numpy.max(numpy.maximum(used - effective, 0.0)))
        integrality_error = float(numpy.max(numpy.abs(result.x - rounded)))
        objective_error = abs(net_profit + float(result.fun))

        Plan(
            scenario,
            units,
            overtime,
            net_profit,
            machine_used,
            labor_used,
            material_used,
            status,
            integrality_error,
            violation,
            objective_error,
        )
    }
}

print("⇒ Planification MILP de trois scénarios")
plans = scenarios.[(scenario) => { solve(scenario) }]
for plan in plans {
    plan.display(products)
    print()
}

# Oracle global : tous les scénarios doivent être optimaux et respecter les
# trois invariants numériques recalculés ci-dessus.
valid = True
for plan in plans {
    if plan.status != SolveStatus.optimal { valid = False }
    if plan.integrality_error > 1.0e-8 { valid = False }
    if plan.capacity_violation > 1.0e-8 { valid = False }
    if plan.objective_error > 1.0e-8 { valid = False }
}
print(f"⇒ Oracle global (optimal + entier + capacités + objectif) : {valid}")

# Planche : mix de production à gauche, taux d'utilisation des ressources à
# droite. L'astérisque indique l'activation de la variable binaire d'heures sup.
plan_label = (plan: Plan): str => {
    suffix = if plan.overtime { '*' } else { '' }
    plan.scenario.name + suffix
}
labels = plans.[(p) => { plan_label(p) }]
x = numpy.arange(len(plans))
fig, axes = plt.subplots(1, 2, figsize=tuple(13, 5.5))

bottom = numpy.zeros(len(plans))
colors = list('#4c78a8', '#f58518', '#54a24b')
i = 0
while i < len(products) {
    values_list = list()
    for plan in plans {
        values_list.append(plan.units[i])
    }
    values = numpy.array(values_list)
    axes[0].bar(x, values, bottom=bottom, label=products[i].name, color=colors[i])
    bottom = bottom + values
    i = i + 1
}
axes[0].set_xticks(x.tolist())
axes[0].set_xticklabels(labels)
axes[0].set_ylabel("unités produites")
axes[0].set_title("Mix optimal (* = heures supplémentaires)")
axes[0].legend()

machine_pct = plans.[(p) => { 100.0 * p.machine_used / p.scenario.machine_capacity(p.overtime) }]
labor_pct = plans.[(p) => { 100.0 * p.labor_used / p.scenario.labor_capacity(p.overtime) }]
material_pct = plans.[(p) => { 100.0 * p.material_used / p.scenario.material }]
width = 0.24
axes[1].bar(x - width, machine_pct, width, label='machine')
axes[1].bar(x, labor_pct, width, label='travail')
axes[1].bar(x + width, material_pct, width, label='matière')
axes[1].axhline(100.0, color='#444444', linewidth=0.8, linestyle='--')
axes[1].set_xticks(x.tolist())
axes[1].set_xticklabels(labels)
axes[1].set_ylim(0.0, 108.0)
axes[1].set_ylabel("utilisation (%)")
axes[1].set_title("Contraintes actives")
axes[1].legend()

fig.tight_layout()
output_path = output_dir / 'production_planning_milp.png'
fig.savefig(str(output_path), dpi=120, bbox_inches='tight')
print(f"⇒ Figure → {output_path}")

# Aperçu navigateur : sert la figure une fois puis rend la main.
# --no-browser garde un chemin headless (le PNG reste écrit ci-dessus).
if '--no-browser' not in import('sys').argv {
    http = import('http')
    b64 = import('base64').b64encode(output_path.read_bytes()).decode('ascii')
    http.serve(f'<!doctype html><meta charset="utf-8"><title>Production planning MILP</title><body style="margin:0;background:#0d1117"><img style="max-width:100%;display:block;margin:0 auto" src="data:image/png;base64,{b64}">',
        0, 'text/html; charset=utf-8', True)
}