#!/usr/bin/env python3

# Big_M_Falle.py
"""
Kapitel MILP: Was ein zu grosses Big-M wirklich anrichtet.

Lehrbuecher warnen vor grossem M mit dem Hinweis "die Relaxation wird
schwach, der Solver braucht mehr Knoten". Das stimmt - ist aber die
harmlosere Haelfte der Wahrheit. Die gefaehrlichere: Bei sehr grossem M
kann der Solver eine Loesung als ganzzahlig ANNEHMEN, in der die
Binaervariablen bei 1e-8 stehen. Dann liefern zugeschaltete Anlagen Ware
aus, waehrend das Modell ihre Fixkosten mit 0 verbucht.

Beispiel: Standortplanung, 12 moegliche Lager, 40 Kunden.
    min  sum_i fix_i * y_i  +  sum_ij kosten_ij * x_ij
    u.d.N.  sum_i x_ij = bedarf_j                    (jeder Kunde beliefert)
            sum_j x_ij <= M_i * y_i                  (Lager offen, wenn es liefert)
            y_i binaer,  x_ij >= 0

Zwei Laeufe:
  1. M knapp gewaehlt (= tatsaechliche Lagerkapazitaet)  -> richtig
  2. M = 1e7-fach zu gross, Presolve abgeschaltet        -> stilles Desaster
Beide werden mit derselben Pruefung kontrolliert, die den Fall auffliegen laesst.

Benoetigt: numpy, highspy
"""

from __future__ import annotations

import time

import numpy as np
import highspy

RNG = np.random.default_rng(11)

N_LAGER, N_KUNDE = 12, 40
FIXKOSTEN = RNG.uniform(2000, 5000, N_LAGER)
TRANSPORT = RNG.uniform(5, 40, (N_LAGER, N_KUNDE))
BEDARF = RNG.uniform(10, 60, N_KUNDE)
KAPAZITAET = BEDARF.sum() * 0.45          # jedes Lager schafft 45 % des Gesamtbedarfs

# Toleranz, ab der ein Solver eine Variable als ganzzahlig durchgehen laesst.
# HiGHS und die meisten anderen verwenden 1e-6 als Standard.
GANZZAHL_TOLERANZ = 1e-6


def loese(big_m: float, presolve: str = "on") -> dict:
    """Baut und loest das Standortmodell. Liefert Loesung und Solverkennzahlen."""
    hochschule = highspy.Highs()
    hochschule.setOptionValue("output_flag", False)
    hochschule.setOptionValue("presolve", presolve)
    hochschule.setOptionValue("time_limit", 300.0)

    anzahl_x = N_LAGER * N_KUNDE
    unendlich = highspy.kHighsInf

    # Spalten: erst alle x_ij, dann die y_i
    hochschule.addVars(anzahl_x, np.zeros(anzahl_x), np.full(anzahl_x, unendlich))
    hochschule.addVars(N_LAGER, np.zeros(N_LAGER), np.ones(N_LAGER))
    for i in range(N_LAGER):
        hochschule.changeColIntegrality(anzahl_x + i,
                                        highspy.HighsVarType.kInteger)
        hochschule.changeColCost(anzahl_x + i, FIXKOSTEN[i])
        for j in range(N_KUNDE):
            hochschule.changeColCost(i * N_KUNDE + j, TRANSPORT[i, j])

    # Jeder Kunde wird genau beliefert
    for j in range(N_KUNDE):
        index = np.array([i * N_KUNDE + j for i in range(N_LAGER)], dtype=np.int32)
        hochschule.addRow(BEDARF[j], BEDARF[j], len(index), index,
                          np.ones(len(index)))

    # Die Kopplung: sum_j x_ij - M * y_i <= 0
    for i in range(N_LAGER):
        index = np.array([i * N_KUNDE + j for j in range(N_KUNDE)] + [anzahl_x + i],
                         dtype=np.int32)
        werte = np.concatenate([np.ones(N_KUNDE), [-big_m]])
        hochschule.addRow(-unendlich, 0.0, len(index), index, werte)

    t0 = time.perf_counter()
    hochschule.run()
    dauer = time.perf_counter() - t0

    loesung = np.array(hochschule.getSolution().col_value)
    info = hochschule.getInfo()
    return {
        "x": loesung[:anzahl_x].reshape(N_LAGER, N_KUNDE),
        "y": loesung[anzahl_x:],
        "zielwert": info.objective_function_value,
        "knoten": info.mip_node_count,
        "dauer": dauer,
    }


def pruefe(ergebnis: dict) -> tuple[bool, list[str]]:
    """Die Pruefung, die in jedes MILP-Auswertungsskript gehoert.

    Sie rechnet die Kosten AUS DER LOESUNG neu aus, statt dem Zielwert des
    Solvers zu glauben - und vergleicht beide. Genau diese Gegenrechnung
    entlarvt eine Loesung, in der Binaervariablen bei 1e-8 haengengeblieben
    sind.
    """
    beanstandungen = []
    x, y = ergebnis["x"], ergebnis["y"]

    # 1. Sind die Binaervariablen wirklich binaer?
    abstand = np.abs(y - np.round(y))
    if abstand.max() > GANZZAHL_TOLERANZ:
        beanstandungen.append(
            f"y ist nicht ganzzahlig: groesster Abstand {abstand.max():.2e}")

    # 2. Liefert ein Lager, dessen Schalter aus ist?
    liefert = x.sum(axis=1) > 1e-6
    geschlossen_aber_aktiv = np.where(liefert & (y < 0.5))[0]
    if geschlossen_aber_aktiv.size:
        beanstandungen.append(
            f"Lager {geschlossen_aber_aktiv.tolist()} liefern Ware, "
            f"gelten im Modell aber als geschlossen")

    # 3. Stimmt der Zielwert mit den echten Kosten ueberein?
    echte_fixkosten = FIXKOSTEN[liefert].sum()
    echte_transportkosten = float((TRANSPORT * x).sum())
    echte_kosten = echte_fixkosten + echte_transportkosten
    if abs(echte_kosten - ergebnis["zielwert"]) > 1e-4 * max(1.0, echte_kosten):
        beanstandungen.append(
            f"Zielwert {ergebnis['zielwert']:,.2f} weicht von den echten Kosten "
            f"{echte_kosten:,.2f} ab (Differenz {echte_kosten - ergebnis['zielwert']:,.2f})")

    return not beanstandungen, beanstandungen


def zeige(titel: str, ergebnis: dict) -> None:
    x, y = ergebnis["x"], ergebnis["y"]
    liefert = x.sum(axis=1) > 1e-6
    print(f"\n{titel}")
    print(f"  Zielwert laut Solver        {ergebnis['zielwert']:>14,.2f} EUR")
    print(f"  Knoten / Zeit               {ergebnis['knoten']:>14,}  "
          f"/ {ergebnis['dauer']:.3f} s")
    print(f"  Lager mit y = 1             {int((y > 0.5).sum()):>14}")
    print(f"  Lager, die tatsaechlich liefern {int(liefert.sum()):>10}")
    unter = y[y < 0.5]
    print(f"  groesster y-Wert unter 0.5  "
          f"{(unter.max() if unter.size else 0.0):>14.3e}")
    print(f"  Fixkosten real / verbucht   {FIXKOSTEN[liefert].sum():>14,.2f} "
          f"/ {float((FIXKOSTEN * y).sum()):,.2f} EUR")

    in_ordnung, beanstandungen = pruefe(ergebnis)
    if in_ordnung:
        print("  PRUEFUNG: bestanden")
    else:
        print("  PRUEFUNG: DURCHGEFALLEN")
        for text in beanstandungen:
            print(f"    - {text}")


if __name__ == "__main__":
    print("=" * 78)
    print("  DIE BIG-M-FALLE: STANDORTPLANUNG, 12 LAGER, 40 KUNDEN")
    print("=" * 78)
    print(f"Tatsaechliche Lagerkapazitaet: {KAPAZITAET:,.1f} Einheiten.")
    print("Genau das ist das kleinstmoegliche gueltige M - mehr kann ein Lager")
    print("ohnehin nicht ausliefern.")

    knapp = loese(KAPAZITAET)
    zeige("[1] M = Kapazitaet (richtig gewaehlt)", knapp)

    gross = loese(1e7 * KAPAZITAET, presolve="off")
    zeige("[2] M = 10 Millionen mal Kapazitaet, Presolve abgeschaltet", gross)

    print("\n" + "=" * 78)
    print("  WAS DA PASSIERT IST")
    print("=" * 78)
    fehlbetrag = knapp["zielwert"] - gross["zielwert"]
    print(f"Lauf [2] meldet {gross['zielwert']:,.2f} EUR und sieht damit um")
    print(f"{fehlbetrag:,.2f} EUR BESSER aus als die richtige Loesung - ein Ergebnis,")
    print("ueber das sich jeder Auftraggeber freuen wuerde.")
    print()
    print("Der Grund steht in der Zeile 'groesster y-Wert unter 0.5': Die")
    print(f"Schaltervariablen stehen bei rund "
          f"{gross['y'][gross['y'] < 0.5].max():.0e}.")
    print(f"Das ist kleiner als die Ganzzahltoleranz {GANZZAHL_TOLERANZ:.0e}, also gilt")
    print("y = 0 - 'Lager geschlossen'. Zugleich ist M so gross, dass")
    print("  sum_j x_ij <= M * 1e-8")
    print("immer noch reichlich Liefermenge erlaubt. Die Lager liefern also,")
    print("ohne dass ihre Fixkosten je bezahlt werden. Der Fachbegriff dafuer")
    print("ist 'trickle flow'.")
    print()
    print("WICHTIG: Mit eingeschaltetem Presolve (Standard) faellt HiGHS hier")
    print("nicht darauf herein - es zieht M selbst zurecht. Verlassen Sie sich")
    print("nicht darauf: Presolve kann das nur, wenn eine implizite Schranke")
    print("herleitbar ist. Die Pruefung aus pruefe() kostet Millisekunden und")
    print("funktioniert immer.")
    print("=" * 78)