415 lines
17 KiB
Python
415 lines
17 KiB
Python
"""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_<n>.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()
|