"""Waldbrand als zellulärer Automat (Modellierung und Simulation, UE1). Kommandozeile: python ue1.py --experiment alle python ue1.py --experiment wind --laeufe 30 --seed 1 python ue1.py --animation Jupyter / Colab: %run -n ue1.py (oder: from ue1 import *) ergebnis = simuliere(Params(wind_staerke=0.8, seed=1)) zeige_verlauf(ergebnis) experiment_wind() Modellregeln (Kurzfassung): - Zustände: Wald, brennend, abgebrannt, gelöscht, Brandschneise. Nur Wald kann sich entzünden; abgebrannt, gelöscht und Brandschneise verhalten sich gleich. - Nachbarschaft: Moore (8 Zellen), Distanzgewicht 1 (direkt) bzw. 0.7 (diagonal). - Rand: die äußerste Zeile bzw. Spalte auf jeder Seite ist Brandschneise, das Feuer kann das Gebiet also nicht verlassen. Jede innere Zelle hat damit iåmmer 8 Nachbarn. - Brandherd: eine zufällige Waldzelle (durch den Seed festgelegt) oder eine vorgegebene Zelle. - Synchrones Update: alle Regeln lesen den Zustand zum Zeitpunkt t und schreiben t+1. - Pro Zeitschritt, für jede brennende Zelle: 1. Mit Löschwahrscheinlichkeit wird sie gelöscht (und breitet sich nicht mehr aus). 2. Sonst versucht sie jeden Nachbarn, der Wald ist, zu entzünden (unabhängige Versuche): p = 0.5 * Distanzgewicht * Brandfaktor(Nachbar) * Windfaktor Windfaktor = 1 + Windstärke * cos(Winkel(Windvektor, Vektor brennende Zelle -> Nachbar)) 3. Danach Restbrenndauer - 1; bei 0 wird sie abgebrannt. Eine Zelle mit Restbrenndauer d breitet sich also in d Zeitschritten aus. - Windvektor (x, y): x zeigt nach rechts (Osten), y nach oben (Norden), der Vektor gibt die Richtung an, in die der Wind weht. """ import argparse import math from dataclasses import dataclass, field, replace from typing import Dict, List, Optional, Sequence, Tuple import matplotlib.pyplot as plt import numpy as np from matplotlib.animation import FuncAnimation from matplotlib.colors import ListedColormap from matplotlib.patches import Patch # Zustände WALD = 0 BRENNEND = 1 ABGEBRANNT = 2 GELOESCHT = 3 SCHNEISE = 4 ZUSTAND_NAMEN = ["Wald", "brennend", "abgebrannt", "gelöscht", "Brandschneise"] ZUSTAND_FARBEN = ["#008300", "#e34948", "#3a3a38", "#86b6ef", "#d9d4c7"] ZUSTAND_CMAP = ListedColormap(ZUSTAND_FARBEN) # Farben für Varianten in Vergleichsplots (feste Reihenfolge) VARIANTEN_FARBEN = ["#2a78d6", "#eb6834", "#1baf7a"] # Moore-Nachbarschaft: (Zeilenversatz, Spaltenversatz, Distanzgewicht) NACHBARN = [ (-1, -1, 0.7), (-1, 0, 1.0), (-1, 1, 0.7), (0, -1, 1.0), (0, 1, 1.0), (1, -1, 0.7), (1, 0, 1.0), (1, 1, 0.7), ] @dataclass class Params: breite: int = 100 hoehe: int = 100 # Anteil natürlicher Brandschneisen (z.B. Wasser, Fels), zufällig verteilt schneise_anteil: float = 0.05 # Künstliche Brandschneisen als volle Zeilen bzw. Spalten schneisen_zeilen: Sequence[int] = () schneisen_spalten: Sequence[int] = () # Brandfaktor ("Trockenheit") je Waldzelle, gleichverteilt brandfaktor_min: float = 0.1 brandfaktor_max: float = 1.0 # Restbrenndauer je Waldzelle, ganzzahlig gleichverteilt (inklusive Grenzen) brenndauer_min: int = 2 brenndauer_max: int = 5 wind_richtung: Tuple[float, float] = (1.0, 0.0) wind_staerke: float = 0.0 loesch_wahrscheinlichkeit: float = 0.05 # Startzelle (Zeile, Spalte), None = zufällige Waldzelle start: Optional[Tuple[int, int]] = None max_schritte: int = 1000 seed: Optional[int] = None @dataclass class Ergebnis: params: Params verlauf: List[np.ndarray] # Zustandsgitter je Zeitschritt, verlauf[0] = Anfangszustand anzahl: np.ndarray # Form (Schritte + 1, 5): Anzahl Zellen je Zustand und Zeitschritt wald_anfang: int = field(init=False) def __post_init__(self): self.wald_anfang = int(self.anzahl[0, WALD] + self.anzahl[0, BRENNEND]) @property def dauer(self) -> int: return len(self.verlauf) - 1 @property def betroffen(self) -> int: """Zellen, die irgendwann gebrannt haben.""" letzte = self.anzahl[-1] return int(letzte[BRENNEND] + letzte[ABGEBRANNT] + letzte[GELOESCHT]) @property def betroffen_anteil(self) -> float: """Anteil der anfänglichen Waldfläche, die gebrannt hat.""" return self.betroffen / self.wald_anfang def windfaktor(dr: int, dc: int, richtung: Tuple[float, float], staerke: float) -> float: """1 + Windstärke * cos(Winkel zwischen Windvektor und Vektor zur Nachbarzelle).""" wx, wy = richtung norm_wind = math.hypot(wx, wy) if staerke == 0 or norm_wind == 0: return 1.0 # Zeilenindex wächst nach unten, y zeigt nach oben nx, ny = dc, -dr cos_winkel = (wx * nx + wy * ny) / (norm_wind * math.hypot(nx, ny)) return 1 + staerke * cos_winkel def verschiebe(a: np.ndarray, dr: int, dc: int) -> np.ndarray: """Liefert b mit b[r, c] = a[r - dr, c - dc]; außerhalb des Gitters 0 (kein Umlauf wie bei np.roll).""" b = np.zeros_like(a) h, w = a.shape b[max(dr, 0):h + min(dr, 0), max(dc, 0):w + min(dc, 0)] = \ a[max(-dr, 0):h + min(-dr, 0), max(-dc, 0):w + min(-dc, 0)] return b def init_gitter(params: Params, rng: np.random.Generator): """Erzeugt Zustand, Brandfaktor und Restbrenndauer und entzündet die Startzelle.""" form = (params.hoehe, params.breite) zustand = np.where(rng.random(form) < params.schneise_anteil, SCHNEISE, WALD).astype(np.int8) if params.start is not None: # Am vorgegebenen Brandherd wächst immer Wald, auch wenn dort zufällig eine natürliche Schneise läge zustand[params.start] = WALD for zeile in params.schneisen_zeilen: zustand[zeile, :] = SCHNEISE for spalte in params.schneisen_spalten: zustand[:, spalte] = SCHNEISE # Rand: die äußerste Zeile bzw. Spalte auf jeder Seite ist Brandschneise zustand[[0, -1], :] = SCHNEISE zustand[:, [0, -1]] = SCHNEISE if params.start is None: # Zufällige Waldzelle; hängt nur vom Seed und den Zelltypen ab, daher bei gleichem Seed # in allen Varianten dieselbe, solange diese die Zelltypen nicht verändern waldzellen = np.argwhere(zustand == WALD) start = tuple(waldzellen[rng.integers(len(waldzellen))]) else: start = params.start if zustand[start] != WALD: raise ValueError(f"Startzelle {start} liegt auf einer künstlichen Brandschneise oder am Rand") brandfaktor = rng.uniform(params.brandfaktor_min, params.brandfaktor_max, form) restbrenndauer = rng.integers(params.brenndauer_min, params.brenndauer_max + 1, form) brandfaktor[zustand == SCHNEISE] = 0 restbrenndauer[zustand == SCHNEISE] = 0 zustand[start] = BRENNEND return zustand, brandfaktor, restbrenndauer def schritt(zustand, brandfaktor, restbrenndauer, params: Params, rng: np.random.Generator): """Ein synchroner Zeitschritt; die Eingabe-Arrays werden nicht verändert.""" brennend = zustand == BRENNEND wald = zustand == WALD geloescht = brennend & (rng.random(zustand.shape) < params.loesch_wahrscheinlichkeit) ausbreitend = brennend & ~geloescht entzuendet = np.zeros_like(wald) for dr, dc, gewicht in NACHBARN: # Zielzelle (r, c) hat eine ausbreitende Nachbarzelle bei (r - dr, c - dc) quelle = verschiebe(ausbreitend, dr, dc) p = 0.5 * gewicht * brandfaktor * windfaktor(dr, dc, params.wind_richtung, params.wind_staerke) entzuendet |= wald & quelle & (rng.random(zustand.shape) < p) neu_zustand = zustand.copy() neu_brandfaktor = brandfaktor.copy() neu_restbrenndauer = restbrenndauer.copy() neu_restbrenndauer[ausbreitend] -= 1 ausgebrannt = ausbreitend & (neu_restbrenndauer == 0) neu_zustand[ausgebrannt] = ABGEBRANNT neu_zustand[geloescht] = GELOESCHT neu_restbrenndauer[geloescht] = 0 neu_brandfaktor[ausgebrannt | geloescht] = 0 neu_zustand[entzuendet] = BRENNEND return neu_zustand, neu_brandfaktor, neu_restbrenndauer def zaehle(zustand: np.ndarray) -> np.ndarray: return np.bincount(zustand.ravel(), minlength=len(ZUSTAND_NAMEN)) def simuliere(params: Params) -> Ergebnis: """Simuliert, bis keine Zelle mehr brennt oder max_schritte erreicht ist.""" rng = np.random.default_rng(params.seed) zustand, brandfaktor, restbrenndauer = init_gitter(params, rng) verlauf = [zustand] anzahl = [zaehle(zustand)] while anzahl[-1][BRENNEND] > 0 and len(verlauf) <= params.max_schritte: zustand, brandfaktor, restbrenndauer = schritt(zustand, brandfaktor, restbrenndauer, params, rng) verlauf.append(zustand) anzahl.append(zaehle(zustand)) return Ergebnis(params, verlauf, np.array(anzahl)) # --------------------------------------------------------------------------- # Visualisierung # --------------------------------------------------------------------------- def zeige_gitter(zustand: np.ndarray, ax=None, titel: str = ""): if ax is None: _, ax = plt.subplots(figsize=(5, 5)) ax.imshow(zustand, cmap=ZUSTAND_CMAP, vmin=0, vmax=len(ZUSTAND_NAMEN) - 1, interpolation="nearest") ax.set_title(titel, fontsize=10) ax.set_xticks([]) ax.set_yticks([]) return ax def zustand_legende(fig): griffe = [Patch(color=farbe, label=name) for name, farbe in zip(ZUSTAND_NAMEN, ZUSTAND_FARBEN)] fig.legend(handles=griffe, loc="lower center", ncol=len(griffe), frameon=False, fontsize=9) def zeichne_kurven(ax, ergebnis: Ergebnis, farbe: str, name: Optional[str] = None): """Betroffene Fläche über die Zeit; ohne Name zusätzlich die brennenden Zellen (gestrichelt).""" a = ergebnis.anzahl betroffen = (a[:, BRENNEND] + a[:, ABGEBRANNT] + a[:, GELOESCHT]) / ergebnis.wald_anfang * 100 ax.plot(betroffen, color=farbe, linewidth=2, label=name or "betroffen gesamt") if name is None: brennend = a[:, BRENNEND] / ergebnis.wald_anfang * 100 ax.plot(brennend, color=farbe, linewidth=2, linestyle="--", label="brennend") def kurven_achse(ax): ax.set_xlabel("Zeitschritt") ax.set_ylabel("% der Waldfläche") ax.set_ylim(0, 100) ax.grid(alpha=0.3) ax.spines[["top", "right"]].set_visible(False) ax.legend(frameon=False, fontsize=9) def zeige_verlauf(ergebnis: Ergebnis, anzahl_bilder: int = 5): """Momentaufnahmen des Gitters und Zeitverlauf eines einzelnen Laufs.""" zeitpunkte = np.linspace(0, ergebnis.dauer, anzahl_bilder).astype(int) fig = plt.figure(figsize=(3 * anzahl_bilder, 6.5)) raster = fig.add_gridspec(2, anzahl_bilder, height_ratios=[1, 0.8]) for i, t in enumerate(zeitpunkte): zeige_gitter(ergebnis.verlauf[t], fig.add_subplot(raster[0, i]), f"t = {t}") ax = fig.add_subplot(raster[1, :]) zeichne_kurven(ax, ergebnis, VARIANTEN_FARBEN[0]) kurven_achse(ax) fig.suptitle(f"Betroffen: {ergebnis.betroffen_anteil:.0%} der Waldfläche, Dauer: {ergebnis.dauer} Schritte") zustand_legende(fig) fig.tight_layout(rect=(0, 0.05, 1, 1)) return fig def zeige_vergleich(ergebnisse: Dict[str, Ergebnis]): """Endzustand jeder Variante nebeneinander plus gemeinsamer Zeitverlauf.""" n = len(ergebnisse) fig = plt.figure(figsize=(4 * n + 5, 4.8)) raster = fig.add_gridspec(1, n + 1, width_ratios=[1] * n + [1.4]) ax_kurven = fig.add_subplot(raster[0, n]) for i, (name, ergebnis) in enumerate(ergebnisse.items()): titel = f"{name}\nbetroffen {ergebnis.betroffen_anteil:.0%}, {ergebnis.dauer} Schritte" zeige_gitter(ergebnis.verlauf[-1], fig.add_subplot(raster[0, i]), titel) zeichne_kurven(ax_kurven, ergebnis, VARIANTEN_FARBEN[i % len(VARIANTEN_FARBEN)], name) kurven_achse(ax_kurven) zustand_legende(fig) fig.tight_layout(rect=(0, 0.07, 1, 1)) return fig def animiere(ergebnis: Ergebnis, intervall_ms: int = 80): """Animation eines Laufs. In Jupyter: from IPython.display import HTML; HTML(anim.to_jshtml())""" fig, ax = plt.subplots(figsize=(6, 6.4)) bild = ax.imshow(ergebnis.verlauf[0], cmap=ZUSTAND_CMAP, vmin=0, vmax=len(ZUSTAND_NAMEN) - 1, interpolation="nearest") ax.set_xticks([]) ax.set_yticks([]) zustand_legende(fig) def update(t): bild.set_data(ergebnis.verlauf[t]) ax.set_title(f"t = {t}") return (bild,) anim = FuncAnimation(fig, update, frames=len(ergebnis.verlauf), interval=intervall_ms) return anim # --------------------------------------------------------------------------- # Experimente # --------------------------------------------------------------------------- def vergleiche(varianten: Dict[str, Params], laeufe: int = 20, seed: int = 0, zeigen: bool = True): """Simuliert jede Variante mit denselben Seeds und gibt Mittelwerte aus. Gleiche Seeds bedeuten: alle Varianten starten mit demselben Wald und Brandherd (gleiche Brandfaktoren, Brenndauern und natürlichen Schneisen), soweit die Variante die Initialisierung nicht selbst ändert. """ statistik = {} beispiele = {} for name, params in varianten.items(): ergebnisse = [simuliere(replace(params, seed=seed + i)) for i in range(laeufe)] anteile = np.array([e.betroffen_anteil for e in ergebnisse]) dauern = np.array([e.dauer for e in ergebnisse]) statistik[name] = { "betroffen_mittel": anteile.mean(), "betroffen_std": anteile.std(), "dauer_mittel": dauern.mean(), "dauer_std": dauern.std(), "grossbrand_anteil": (anteile > 0.5).mean(), } # Beispiel für die Grafik: der Lauf mit dem Median der betroffenen Fläche beispiele[name] = ergebnisse[int(np.argsort(anteile)[len(anteile) // 2])] print(f"{'Variante':<28}{'betroffen (Mittel ± Std)':>26}{'Dauer (Mittel ± Std)':>24}{'Läufe > 50 %':>14}") for name, s in statistik.items(): print(f"{name:<28}{s['betroffen_mittel']:>17.1%} ± {s['betroffen_std']:>5.1%}" f"{s['dauer_mittel']:>16.1f} ± {s['dauer_std']:>5.1f}{s['grossbrand_anteil']:>14.0%}") print() if zeigen: zeige_vergleich(beispiele) return statistik, beispiele def experiment_wind(laeufe: int = 20, seed: int = 0, zeigen: bool = True): print("Experiment Wind: Windstärke 0 vs. 0.8 (Wind nach Osten)") basis = Params() return vergleiche({ "ohne Wind": basis, "Wind 0.8 nach Osten": replace(basis, wind_richtung=(1, 0), wind_staerke=0.8), }, laeufe, seed, zeigen) def experiment_schneise(laeufe: int = 20, seed: int = 0, zeigen: bool = True): print("Experiment Brandschneise: Wind nach Osten, ohne vs. mit Schneise östlich des Brandherds") # Fester Brandherd westlich der Schneise, damit der Vergleich nicht vom zufälligen Startpunkt abhängt basis = Params(wind_richtung=(1, 0), wind_staerke=0.8, start=(50, 50)) # Eine durchgehende gerade Schneise ist in der Moore-Nachbarschaft schon mit 1 Zelle Breite # unüberwindbar, die Breite spielt in diesem Modell also keine Rolle return vergleiche({ "ohne Schneise": basis, "Schneise in Spalte 65": replace(basis, schneisen_spalten=(65,)), }, laeufe, seed, zeigen) def experiment_trockenheit(laeufe: int = 20, seed: int = 0, zeigen: bool = True): print("Experiment Trockenheit: Brandfaktor feucht (0.1-0.5) vs. normal (0.1-1.0) vs. trocken (0.5-1.0)") basis = Params() return vergleiche({ "feucht (0.1-0.5)": replace(basis, brandfaktor_min=0.1, brandfaktor_max=0.5), "normal (0.1-1.0)": basis, "trocken (0.5-1.0)": replace(basis, brandfaktor_min=0.5, brandfaktor_max=1.0), }, laeufe, seed, zeigen) EXPERIMENTE = { "wind": experiment_wind, "schneise": experiment_schneise, "trockenheit": experiment_trockenheit, } def main(): parser = argparse.ArgumentParser(description="Waldbrand als zellulärer Automat") parser.add_argument("--experiment", choices=[*EXPERIMENTE, "alle"], help="Experiment ausführen (ohne Angabe: einzelner Lauf mit Standardparametern)") parser.add_argument("--laeufe", type=int, default=20, help="Läufe je Variante") parser.add_argument("--seed", type=int, default=0) parser.add_argument("--wind", type=float, default=0.0, help="Windstärke für den Einzellauf") parser.add_argument("--animation", action="store_true", help="Einzellauf animieren") parser.add_argument("--speichern", metavar="PRAEFIX", help="Grafiken als PRAEFIX_.png speichern statt anzeigen") args = parser.parse_args() anim = None if args.experiment: namen = list(EXPERIMENTE) if args.experiment == "alle" else [args.experiment] for name in namen: EXPERIMENTE[name](args.laeufe, args.seed) else: ergebnis = simuliere(Params(wind_staerke=args.wind, seed=args.seed)) if args.animation: anim = animiere(ergebnis) # Referenz halten, sonst stoppt die Animation else: zeige_verlauf(ergebnis) if args.speichern: for i in plt.get_fignums(): plt.figure(i).savefig(f"{args.speichern}_{i}.png", dpi=150) else: plt.show() return anim if __name__ == "__main__": main()