Tag Python

Eine neue Metrik für Webseiten: Crapicity

2023-02-08
Die Krapizität der Webseiten, auf die ich hier so verlinkt habe (und noch ein paar mehr): works with netsurf!

Kurz nachdem das Web seine akademische Unschuld verloren hat, Ende der 1990er, habe ich den UNiMUT Schwobifying Proxy geschrieben, ein kleines Skript, das das ganze Web in Schwäbisch, hust, erlebbar machte. Tatsächlich hat mir das gegen 2002 meine 15 Minuten des Ruhms verschafft, inklusive 200'000 Zugriffen pro Tag (fürs damalige Netz rasend viel), Besprechungen auf heise.de, spiegel.de und, für mich offen gestanden am schmeichelndsten, in Forschung aktuell (also gut: Computer und Kommunikation) im Deutschlandfunk.

Ein heimlich verwandtes Web-Experiment war dass Dummschwätzranking von 1998, das – nicht völlig albern – die Dichte von (damals gerade modernen) Heißdampfworten in Webseiten beurteilt hat – dafür interessieren[1] sich bis heute Leute.

Beide Spielereien sind praktisch abgestellt, teils, weil das kommerzielle Internet und SEO solche Sachen rechtlich oder praktisch sprengen. Teils aber auch, weil sie darauf bauen, dass das, was Leute im Browser sehen, auch in etwa das ist, was im HTML drinsteht, das so ein Programm von den Webservern bekommt. Das aber ist leider im Zuge der Javascriptisierung des Web immer weniger der Fall.

Nun ist sicherlich die Habituierung der Menschen an „lass einfach mal alle Leute, von denen du was lesen willst, Code auf deiner Maschine ausführen“ die deutlich giftigere Folge des Post-Web-1-Megatrends. Aber schade ist es trotzdem, dass zumindest im kommerziellen Web kein Mensch mehr in die Seiten schreibt, was nachher in den dicken Javascript-Browsern angezeigt werden wird.

Und weil das schade ist, habe ich einen postmodernen Nachfolger des Dummschwätzrankings geschrieben: Die Crapicity-Maschine (sorry, nicht lokalisiert, aber sprachlich sowieso eher arm). Ihre Metrik, namentlich die Crapicity oder kurz c7y: das Verhältnis der Länge lesbaren Textes (das, was ohne Javascript an Zeichen auf den Bildschirm zu lesen ist) zur Gesamtlänge der Seite mit allem eingebetteten Markup, Javascript und CSS (also: externes Javascript, CSS usf nicht gerechnet). In Python mit dem wunderbaren BeautifulSoup-Modul ist das schnell berechnet:
```
def compute_crapicity(doc):
  """returns the crapicity of html in doc.

  doc really should be a str -- but if len() and BeautifulSoup() return
  something sensible with it, you can get away with something else, too.
  """
  parsed = BeautifulSoup(doc, "html.parser")
  content_length = max(len(parsed.text), 1)
  return len(doc)/content_length
```
Um diese knappe Funktion herum habe fast 700 Zeilen herumgeklöppelt, die Ergebnisse in einer SQLite-Datenbank festhalten und ein Webinterface bereitstellen. Als Debian-und-eine-Datei-Programm sollte das recht einfach an vielen Stellen laufen können – wer mag, bekommt die Software bei codeberg.

Die Web-Schnittstelle bei https://blog.tfiu.de/c7y hat aber den Vorteil, dass sich die Scores sammeln. Spielt gerne damit rum und empfiehlt es weiter – ich fände es putzig, da vielleicht 10'000 Seiten vertreten zu haben. Ich habe selbst schon knapp 200 Web-Ressourcen durchgepfiffen, meistenteils Links aus dem Blog hier.

Im Groben kommt raus, was wohl jedeR erwartet hat: Das kommerzielle Netz stinkt, alter Kram und Techno-Seiten sind meist ganz ok. Allerdings habe ich den aktuellen Spitzenreiter, eine Reddit-Seite mit einem c7y von über 17'000, noch nicht ganz debuggt: Erstaunlicherweise ist die Seite auch im netsurf und ohne Javascript lesbar. Wie sie das macht: nun, das habe ich in 800 kB Wirrnis noch nicht rausgefunden, und der Quellcode der Seite sieht so schrecklich aus, dass der Score sicherlich verdient ist.

Ich nehme mal an, dass derzeit alle youtube-Seiten bei c7y=8222 liegen; dort ist ja durchweg ohne Javascript nichts zu sehen, und so haut auch dieser Score prima hin. Bei taz.de (gerade 892) geht es vielleicht nicht so gerecht zu, denn die Seite funktioniert in der Tat auch ohne Javascript ganz gut. Eventuell ist hier BeautifulSoup schuld. Hochverdient hingegen sind die 682 von nina.no – das ist ohne Javascript leer. Eine Twitter-Seite liegt bei 413, Bandcamp bei 247.

Vernüftige Seiten liegen dagegen zwischen etwas über eins (minimales Markup) und zehn (z.B. wenig Text mit viel CSS). Dishonorable Mention: uni-heidelberg.de liegt trotz Akademia bei 177. Tatsächlich ist die Seite auch in normalen Browsern[2] halbwegs lesbar. Der schlechte Score liegt vor allem an eingebetten SVGs, ist also schon ein ganz klein wenig unfair. Aber ehrlich: wer für ein bisschen Glitz ein paar hundert Zeichen Text auf satte 680k aufbläht, hat eine große crapicity verdient, auch wenn die Seite selbst nicht richtig kaputt ist. Wer unbedingt so viel Glitz haben will, soll externe Bilder verwenden – die muss ich nicht runterladen, wenn ich nicht will.

Wer interessante Krapizitäten findet: Die Kommentarbox wartet auf euch.

[1] Na gut, viel von dem Interesse kam aus erkennbar aus SEO-Kreisen; das war dann auch einer der Gründe, warum ich das eintragen von Links beim Dummschwätzranking abgestellt habe.

[2] Definiert als: Alles außer den Monstren Firefox, Chrome, Webkit und ihren Derivaten.
⎆
Speech Recognition with Whisper.cpp

2023-02-02
Today I stumbled across Whispers of A.I.'s Modular Future by James Somers, a piece that, at least by the standards of publications aimed at the general public, makes an excellent point of why whisper.cpp might finally be some useful and non-patronising output of the current AI hype.

What can I say? I think I'm sold. And perhaps I'm now a little bit scared, too. If you want to understand way and speak a bit of German, you can skip to The Crazy right away.
The Good

You know, so far I've ignored most of the current statistical modelling (“AI”, “Machine Learning“) – if you need a graphics chip with drivers even worse than Intel's, and that then needs 8 GB of video RAM before anything works, I'm out. And I'm also out when the only way I can use some software is on some web page because there's proprietary data behind it.

Not so for whisper.cpp. This is software as it was meant to be: trivial dependencies, compact, works on basically any hardware there is. To build it, you just run:
```
git clone https://github.com/ggerganov/whisper.cpp/
cd whisper.cpp
make
```
– and that's it. No dependency juggling down to incompatible micro versions, no fancy build system, just a few C(++) sources and a Makefile. The thing works in place without a hitch, and it has a sensible command line interface.

Well, you need the language models, of course. There are some reasonably free ones for English. The whisper.cpp distribution's models/README.md explains how to obtain some. I got myself ggml-small.en.bin, recorded a few words of English into a file zw.wav and ran:
```
./main -m models/ggml-small.en.bin ~/zw.wav
```
The machine demanded I use a samplerate of 16 kHz, I made audacity oblige, ran the thing again and was blown away when – admittedly after a surprisingly long time – my words appeared on the screen.

I immediately tried to figure out how to stream in data but then quickly decided that's probably not worth the effort; the software needs to see words in context, and for what I plan to do – transcribing radio shows – having an intermediate WAV file really does not hurt.

I quickly cobbled together a piece of Python wrapping the conversion (using the perennial classic of audio processing, sox) somewhat cleverly, like this:
```
#!/usr/bin/python
# A quick hack to transcribe audio files
#
# Dependencies:
# * sox (would be mpv, but that's somehow broken)
# * a build of whispercpp (https://github.com/ggerganov/whisper.cpp/)
# * a language model (see models/README.md in the whisper source)

import contextlib
import os
import subprocess
import sys
import tempfile

WHISPER_DIR = "/usr/src/whisper.cpp"


@contextlib.contextmanager
def workdir(wd):
        prev_dir = os.getcwd()
        try:
                os.chdir(wd)
                yield
        finally:
                os.chdir(prev_dir)


def transcribe(audio_source, model, lang):
        """transcibes an audio file, creating an in-place .txt.

        model must be the name of a model file in WHISPER_DIR/models;
        lang is the ISO language code in which the output should turn up.
        """
        audio_source = os.path.join(os.getcwd(), audio_source)
        with tempfile.TemporaryDirectory(suffix="transcribe", dir="/var/tmp") as wd:
                with workdir(wd):
                        subprocess.check_call(["sox",
                                audio_source,
                                "-b", "16", "-r", "16000", "-c", "1",
                                "audiodump.wav"])

                        out_name = os.path.splitext(audio_source)[0]
                        subprocess.check_call([WHISPER_DIR+"/main",
                                "-l", lang,
                                "-m", WHISPER_DIR+"/models/"+model,
                                "-otxt", "-of", out_name,
                                "audiodump.wav"])


def parse_command_line():
        import argparse
        parser = argparse.ArgumentParser(description="Wrap whisper.cpp to"
                " bulk-transcribe audio files.")
        parser.add_argument("model", type=str, help="name of ggml language"
                f" model to use, relative to {WHISPER_DIR}/models")
        parser.add_argument("audios", type=str, nargs="+",
                help="Sox-translatable audio file to transliterate.")
        parser.add_argument("--lang", type=str, default="en",
                help="Spoken language to try and recogonise")

        return parser.parse_args()


if __name__=="__main__":
        args = parse_command_line()
        for audio in args.audios:
                transcribe(audio, args.model, args.lang)
```
Nachtrag (2023-06-26)

(Added a --lang option as per ron's feedback below)

I have that as transcribe.py in my path, and I can now enter the rip of an audiobook and say:
```
transcribe.py ggml-small.en.bin *.ogg
```
(provided I have downloaded the model as per whisper.cpp's instructions). After a little while (with high CPU usage), there is a transcript on my disk that's better what I had typed myself even after two rounds of proff-reading, except that whisper.cpp doesn't get the paragraphs right.

For the first time in the current AI hype, I start getting carried away, in particular when I consider how much speech recognition sucked when I last played with it around 2003, using a heap of sorry failure called viavoice.
The Bad

Skip the rant to get to the exciting part.

Trouble is: What I'd mainly like to transcribe is German radio, and whisper.cpp does not come with a German language model. Not to worry, one would think, as whisper.cpp comes with conversion scripts for the pyTorch-based whisper models like those one can get from Hugging Face. I downloaded what I think is the model file and cheerfully ran:
```
$ python convert-h5-to-ggml.py /media/downloads/model.bin
Traceback (most recent call last):
  File "/home/src/whisper.cpp/models/convert-h5-to-ggml.py", line 24, in <module>
    import torch
ModuleNotFoundError: No module named 'torch'
```
Oh bummer. Well, how hard can it be? Turns out: Surprisingly hard. There is no pytorch package Debian stable. Ah… I very much later realised there is, it's just that my main system still has an i386 userland, and pytorch is only available for amd64. But I hadn't figured that out then. So, I enabled a virtual python (never mix your system python and pip) and ran:
```
$ pip install torch
ERROR: Could not find a version that satisfies the requirement torch
ERROR: No matching distribution found for torch
```
Huh? What's that? I ran pip with a couple of -v sprinkled in, which at least yielded:
```
[...]
Skipping link: none of the wheel's tags match: cp38-cp38-win_amd64: https://download.pytorch.org/whl/cpu/torch-1.9.0%2Bcpu-cp38-cp38-win_amd64.whl (from https://download.pytorch.org/whl/cpu/torch/)
[...]
Given no hashes to check 0 links for project 'torch': discarding no candidates
ERROR: Could not find a version that satisfies the requirement torch
ERROR: No matching distribution found for torch
[...]
```
The message with “Given no“ has a certain lyric quality, but other than that from the “Skipping“ messages I concluded they don't have 32 bit builds any more.

Well, how hard can it be? Pypi says the sources are on github, and so I cloned that repo. Oh boy, AI at its finest. The thing pulls in a whopping 3.5 Gigabytes of who-knows-what. Oh, come on.

python setup.py build fails after a short while, complaining about missing typing_extensions. Manually running pip install typing_extensions fixes that. But I killed setup.py build after a few minutes when there were only 50/5719 files built. Has AI written that software?

In the meantime, I had gone to a machine with a 64 bit userland, and to be fair the experience wasn't too bad there, except for the hellish amount of dependencies that pytorch pulls in.

So, my expectations regarding “AI code” were by and large met in that second part of the adventure, including the little detail that the internal links on https://pypi.org/project/torch/ are broken because right now their document processor does not produce id attributes on the headlines. Yeah, I know, they're giving it all away for free and all that. But still, after the brief glimpse into the paradise of yesteryear's software that whisper.cpp afforded, this was a striking contrast.
The Crazy

So, I converted the German language model doing, in effect:
```
git clone https://github.com/openai/whisper.git
git lfs install
git clone https://huggingface.co/bofenghuang/whisper-small-cv11-german
python convert-h5-to-ggml.py whisper-small-cv11-german/ whisper tmp
```
(where I took convert-h5-to-ggml.py from whisper.cpp's repo). Then I moved the resulting tmp/ggml-model.bin to german-small.ggml and ran:
```
transcribe.py german-small.ggml peer_review_wie_objektiv_ist_das_wissenschaftliche_dlf_20221214_1646_8a93e930.mp3
```
with my script above and this German-language mp3 from Deutschlandfunk. From the English experience, I had expected to get an almost flawless transliteration of the German text. What I got instead was (paragraphs inserted by me); listen to the audio in parallel if you can:

Germany. Research is on [that was: Deutschlandfunk Forschung aktuell]

A Nobel Prize for Science is not easy without further ado. They really need to find something out. For example, Vernon Smith, who is now 95 years old, is now the father of the Experimental Economy. In 2002 he won the Nobel Prize for Science.

This made such a prize and renommee also make impression on other Fachleuteen and that actually influenced the unabhängig well-office method for scientific publications. This has recently shown a study of Business Science in the Fachmagazin PNS. Anike Meyer spoke with one of the authors.

When Jürgen Huber and his colleagues thought about the experiment, it was clear to them that this is not fair. The same manuscript was given by two different authors, Vernon …
⎆
Neun Monate Umwelt-CO₂, Teil II: Hochpass, Tiefpass, Spektrum

2023-01-25
Am Ende des Posts verrate ich, welche Bewandnis es mit diesem hübschen, nachgerade frühlingshaften Muster hat. Und auch, wie mensch sowas selbst macht.
Ich habe letztes Jahr neun Monate lang CO₂-Konzentrationen auf meinem Balkon gemessen, in der Hoffnung, ein wenig hinter die Gründe für die doch überraschend großen Konzentrationsschwankungen zu kommen, für die ich bei einer ersten Messung im November 2021 nur sehr spekulative Erklärungen gefunden habe (Stand Februar vor allem: die lokale Heizung und eventuell das Kraftwerk in Rheinau oder die Chemiefabriken in Ladenburg oder Ludwigshafen). Ich habe dann neulich an der Kalibration der Daten gespielt und bin zum Schluss gekommen, dass sie nach Maßstäben von KonsumentInnenelektronik recht ordentlich sein dürften.

Jetzt möchte ich ein wenig in den Daten rumfummeln. „Explorative Datenanalyse“ nennen das Leute, die das mit dem Fummeln nicht zugeben wollen, und in der ernsthaften Wissenschaft wird zurecht etwas die Nase gerümpft, wenn jemand sowas macht: Es stecken in jedem hinreichend großen Datensatz fast beliebig viele Korrelationen, und wer sie sucht, wird sie auch finden. Leider sind (fast) alle davon Eigenschaften des Datensatzes bzw. der Messung, nicht aber des untersuchten Gegenstands.

Andererseits: Ohne Induktion (was in normale Sprache übersetzt „Rumspielen, bis mensch einen Einfall hat“ heißt) gibt es, da muss ich Karl Popper ganz heftig widersprechen, keine Erkenntnis, und Deduktion kann mensch nachher immer noch machen (also, ich jetzt nicht, weil das hier mein Freizeitvergnügen ist und keine Wissenschaft).
Erstmal glätten

Meine ganz große Enttäuschung mit dem Datensatz war ja, dass sich fast kein Jahreszeiteneffekt gezeigt hat; natürlich werde ich hier, mitten im dichtbesiedelten Oberrheingraben, kein so hübsches Signal bekommen wie die Leute mit der klassischen Messung am Mauna Loa – davon, dass ich eine Schwingung von ein paar ppm nachvollziehen könnte, ganz zu schweigen. Aber ich hatte im September 2021 unter 300 ppm gemessen, während es im Herbst auf die global eher aktuellen 400 bis 500 ppm hochging. Ich hatte gehofft, den Weg zurück zu den 300 ppm im Laufe des Sommers beobachten zu können. Doch leider ist in den (Roh-) Daten erstmal nur Gekrakel:

Wenn in einem Datensatz nur Gekrakel ist, es aber einen langfristigen Effekt geben soll, hilft oft Glätten. Dabei ersetzt mensch jeden Punkt durch einen geeignet gebildeten Mittelwert über eine größere Umgebung dieses Punktes, wodurch kurzfristiges Gewackel wie in meinen Rohdaten stark unterdrückt wird und die langfristigen Trends besser sichtbar werden sollte. Es gibt noch einen zweiten Gewinn: Mensch kann das so gewonnene Mittel von den Daten abziehen und hat sozusagen destilliertes Gewackel (den „hochfrequenten Anteil“), der auf diese Weise auch ein klareres Signal zeigen mag.

Im einfachsten Fall lässt sich so eine Glättung in ein paar Zeilen Python schreiben, und das auch noch relativ effizient als gleitendes Mittel unter Verwendung einer zweiendigen Queue aus Pythons großartigem Collections-Modul: Um das Mittel zu berechnen, addiere ich auf einen Akkumulator, und das, was ich da addiere, muss ich wieder abziehen, wenn der entsprechende Wert aus dem Fenster rausläuft. Dazu gibts noch die Subtilität, dass ich Zeit und Konzentration in der Mitte des Fensters zurückgeben will. Das führt auf so eine Funktion:
```
def iter_naively_smoothed(in_file, smooth_over):
    queue = collections.deque()
    accum = 0

    for time, co2 in iter_co2(in_file):
        queue.append((time, co2))
        accum += co2
        while time-queue[0][0]>smooth_over:
            _, old = queue.popleft()
            accum -= old
        time, base_co2 = queue[len(queue)//2]
        yield time, accum/len(queue), base_co2-accum/len(queue)
```
Als Plot sehen die über drei Tage[1] geglättete Konzentrationskurve und die Residuen so aus:

Die Ränder, vor allem der linke, sind dabei mit Vorsicht zu genießen, da dort fast nichts geglättet ist; die Funktion ist so geschrieben, dass dort die Fenster klein sind. So oder so: auch die geglättete Kurve hat keine nennenswerte Tendenz. Allenfalls die beiden großen Ausschläge bei den Tagen 80 und 110 könnten ein Signal sein.
Zeitrechnungen

„Tag 80” ist jetzt keine wirklich intuitive Angabe, aber ich wollte in den Plots nicht mit läsitgen Kalenderdaten operieren. Für die Interpretation wären sie jedoch schon ganz gut. Wie komme ich von Messtag auf ein bürgerliches Datum?

Die Daten kommen original mit Unix-Timestamps, also der Zahl der Sekunden seit dem 1.1.1970, 0:00 Uhr UTC, und sie fangen mit 1640612194 an. Um ausgehend davon bürgerliche Daten auszurechnen, muss mensch wissen, dass ein Tag 86400 Sekunden hat. Pythons datetime-Modul liefert dann:
```
>>> import datetime
>>> datetime.datetime.fromtimestamp(1640612194+80*86400)
datetime.datetime(2022, 3, 17, 14, 36, 34)
>>> datetime.datetime.fromtimestamp(1640612194+110*86400)
datetime.datetime(2022, 4, 16, 15, 36, 34)
```
Mitte März und Mitte April ist also möglicherweise was Interessantes passiert. Da gucke ich nochmal drauf, wenn ich demnächst Wetterdaten einbeziehe.

Die andere Seite, die Residuen, kann ich für eine saubere Fassung des Tagesplots aus Teil 1 verwenden. Weil längerfristige Schwankungen bei ihnen rausgerechnet sind, sollte ein regelmäßiger Tagesgang in den Residuen deutlicher herauskommen als im ganzen Signal. Dazu nehme ich wieder den Rest bei der Division durch 86400 (siehe oben), und die Warnungen wegen Zeitzonen aus dem Kalibrationspost gelten immer noch.

Dass da „am Morgen” (15000 Sekunden sind etwa 4:15 Uhr UTC, also 5:15 MEZ und 6:15 MESZ) irgendwas passiert, dürfte ein robustes Signal sein, und es kommt ohne die niederfrequenten Signale deutlich besser raus. Mit etwas Wohlwollen könnte mensch vielleicht sogar zwei Berge im Abstand von einer Stunde sehen, was dann Normal- und Sommerzeit entspräche. Da der Sensor nicht weit von der Bundesstraße 3 stand, liegt als Erklärung zunächst Berufs- und Pendelverkehr nahe.

Wenn der Tagesgang tatsächlich mit dem Berufsverkehr am Morgen zu tun hat, müsste sich eigentlich ein Signal im Wochenrhythmus zeigen, denn zumindest Sonntags ist es hier am Morgen deutlich ruhiger als an Werktagen. Die Projektion auf die Woche ist einfach: Statt den Rest bei der Division durch 86'400 nehme ich jetzt den Rest bei der Division durch 604'800, die Zahl der Sekunden in einer Woche. Das Ergebnis zeigt, vielleicht etwas überraschend, das Tagessignal eher deutlicher; ein Wochensignal hingegen ist weniger überzeugend erkennbar:

Vielleicht ist an den Tagen drei und vier etwas weniger los – welche Wochentage sind das? Nun, ich dividiere hier Unix-Timestamps; ihr erinnert euch: Die Sekunden seit dem 1.1.1970. Welcher Wochentag war dieser 1.1.?
```
$ cal 1 1970
    January 1970
Su Mo Tu We Th Fr Sa
             1  2  3
 4  5  6  7  8  9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30 31
```
Weil die Unix-Epoche so eine große Rolle im heutigen Leben spielt[2], nehme das mal als erweiterte Kopfzahl: Er war ein Donnerstag. Und so sind die Tage 3 und 4 tatsächlich Samstag und Sonntag. Meine Autoabgas-These ist mir durch diese Prüfung etwas sympathischer geworden. Vermutlich fallen die Abgase nur bei wenig Wind und wenig Thermik (das ist praktisch die berüchtigte Dunkelflaute…) auf, so dass das dunkle Band mit hoher Punktdichte mehr oder minder bei der Null bleibt, jedoch immer wieder Spitzen mit bis zu 100 ppm extra in der Rush Hour auftreten. Wenn das hier Wissenschaft wäre, müsste ich mich jetzt auf die Suche nach Daten aus der Verkehrszählung machen.

Und ich müsste überlegen, ob die Zacken im dunklen Band nicht doch ein echtes Signal sind. Speziell der Durchhänger so um die 180'000 Sekunden (also in der Nacht von Freitag auf Samstag) ist eigentlich kaum wegzudiskutieren. Aber was könnte gerade da plausiblerweise mehr frische Luft heranführen? Oder ist das, weil die Freitagnacht wegen Spätcorona noch besonders ruhig war?

In Sachen Spitzen am Morgen hätte ich eine Alternativhypothese zu den Autoabgasen: Das könnte nämlich wieder der lokale Brenner sein. Dann wäre das, was da zu sehen ist, die Warmwasserbereitung zur Morgendusche. Attraktiv ist diese These schon allein, weil mir 6:15 eigentlich ein wenig früh vorkommt für das Einsetzen der Rush Hour nach Heidelberg rein; es ist ja nicht so, als seien da noch viele nennenswerte Industriebetriebe, in denen die Leute um sieben stechen müssten. Allerdings: in einem Wohnblock mit so vielen Studis wie meinem hier ist 6:15 auch keine plausible Zeit für massenhaftes Duschen…
Spektralanalyse

Was ich gerade gemacht habe, ist ein Spezialfall einer gerade für uns AstrophysikerInnen extrem wichtigen Technik – ich habe nämlich ein Signal nach Frequenzen aufgeteilt. Bis hierher hatte ich nur in zwei Bänder, ein hohes und ein tiefes. Mit der Fourieranalyse geht das viel allgemeiner: Mensch steckt ein Signal rein und bekommt ein Spektrum heraus, eine Kurve, die für jede Frequenz sagt, wie stark Schwingungen auf dieser Frequenz zum Signal beitragen.

Mit dem Universal-Werkzeugkasten scipy kriegt mensch so ein Spektrum in anderthalb Zeilen Python:
```
from scipy import signal
f, power = signal.periodogram(samples[:,1], fs=86400/24/120)
```
Dabei übergebe ich nur meine Konzentrationsmessungen, nicht die zugehörigen Zeiten, weil die periodogram-Funktion nur mit gleichmäßig gesampleten Daten zurechtkommt. Das sind meine Daten zum Glück: zwischen zwei Messungen liegen immer 30 Sekunden, und es sind auch keine Lücken in den Daten.

Die halbe Minute Abstand zwischen zwei Punkten übergebe ich im Argument fs, allerdings als Frequenz (der Kehrwert der Periode) und umgerechnet in die Einheit „pro …
⎆
Neun Monate Umwelt-CO2, Teil I: Taugen die Daten?

2022-12-25

Im vergangenen Jahr habe ich meine CO₂-Messung am Balkon bis Mitte September laufen lassen, vor allem, weil ich sehen wollte, wie die Konzentrationen im Laufe der Zeit auf unter 300 ppm sinkt. So weit unten (relativ zum gegenwärtigen globalen Mittelwert von um die 400 ppm) lag zu meiner damaligen Überraschung die Konzentration mal ganz am Anfang meiner CO₂-Messungen, im September 2021. Ich hatte eigentlich erwartet, dass all das Grün im Sommer nach und nach auch in diesem Jahr dafür sorgen würde.

Daraus ist nichts geworden. Im Groben ist die CO₂-Konzentration über Frühling und Sommer 2022 konstant geblieben:

(Alle Plots in diesem Post von TOPCAT). Mag sein, dass es einfach zu trocken war in diesem Jahr – mag aber auch sein, dass meine ersten Experimente einfach in besonders frischer Luft stattfanden. Schließlich kämen natürlich noch Kalibrationsprobleme in Betracht; ich habe nicht versucht, meine Messungen mit denen anderer zu abzugleichen.

Der Negativbefund hat mich aber dazu gebracht, die Daten im Hinblick auf statistische und vor allem systematische Fehler genauer unter die Lupe zu nehmen. Darum geht es in diesem Post.

Zunächst stellt sich die Frage, ob die generelle Zittrigkeit der Kurve eigentlich Rauschen des Sensors ist oder etwas anderes – von Anfang an haben mich ja die teils erheblichen Konzentrationsschwankungen verblüfft. Nun: angesichts der hohen Korrelation benachbarter Messwerte kommen diese jedenfalls nicht aus statistischen Fehlern im Sensor. Ich greife mal den 26. Mai (einen Donnerstag) heraus:

Wenn das Wackeln ein statistischer Fehler wäre, dann wäre die Linie, die ich durch die Punkte gemalt habe, völliges Gekrakel und nicht erkennbarer Verlauf. Ich will gerne glauben, dass da ein Rauschen irgendwo unterhalb von 10 ppm drin ist. Der Rest ist irgendein Signal.

Welcher Natur das Signal ist, ist eine Frage, die sich allenfalls mit dem ganzen Datensatz beantworten lässt. Angesichts der überragenden Bedeutung der Sonne für Physik und Leben auf der Erde geht der erste Blick auf der Suche nach systematischen Fehlern oder auch echter Physik bei solchen Zeitreihen erstmal auf die Tagesverläufe. Um da eine Idee des Verhaltens des ganzen Datensatzes zu bekommen, habe ich mir die Punktdichte in einem Plot von Tageszeit gegen CO₂-Konzentration angesehen (ich „falte auf den Tag“):

Die Abszisse hier verdient einen kurzen Kommentar: Sie zeigt den Rest bei der Division meiner Timestamps durch 84600, ausgedrückt in Stunden. Meine Timestamps sind in Sekunden, und 84600 ist einfach 24 mal 3600, also die Zahl der Sekunden an einem Tag. Mithin steht hier etwas wie eine Tageszeit.

Ganz so einfach ist das aber nicht, weil meine Timestamps immer in UTC laufen, während die Umgebung der Willkür der Zeitzonen unterworfen ist; die 15 auf der Abszisse entspricht also manchmal 16 Uhr bürgerlicher Zeit (nämlich, wenn die Umgebung MEZ hatte) und manchmal 17 Uhr (wenn ihr Sommerzeit verordnet war). Aber schon rein optisch liegt nicht nahe, dass viel mehr zu sehen wäre, wenn ich die politischen Kapriolen nachvollziehen würde, um so weniger, als die Sonne die ja auch nicht nachvollzieht.

Die Bäuche in der dunkleren Fläche, also einer besonders hohen Dichte von Punkten, entsprechen nun Tageszeiten, zu denen es häufiger mal hohe CO₂-Konzentrationen auf meinem Balkon gab. Das könnte ein Signal der Lüftung unserer Wohnung (oder vielleicht sogar der unserer Nachbarn) sein. Es ist aber auch plausibel, dass es der Reflex der Verkehrsdichte ist. Der Balkon befindet sich etwa 10 Meter über und 20 Meter neben einer recht viel befahrenen Straße, weshalb im September 2021 die Unsichtbarkeit der CO₂-Emissionen der Fahrzeuge mit meine größte Überraschung war. Mit hinreichend viel Statistik und Mitteln über hinreichend viele Wetterlagen zeigen sich die Autos (vielleicht) eben doch.

Die Messanordnung; das wird für die Kalibration noch wichtig…

Wenn Effekte so sehr auf zusammengekniffenen Augen beruhen wie hier beim Tagesverlauf, hilft es nichts: Da braucht es einen zweiten Blick auf innere Korrelationen der Daten, um zu sehen, ob sich da schlicht systematische Fehler zeigen oder ob es wirklich die Autos sind. Klar: besser wäre es natürlich, mit bekannten Konzentrationen oder einfacher einem bekannt guten Messgerät zu kalibrieren, also zu sehen, welche Anzeige das Gerät für welchen wahren Wert hat.

Aber das ist aufwändig, und zumeist zeigen sich Systematiken mit ein paar plausiblen Annahmen („Modelle“) auch schon in den Daten selbst. Zur plausiblen Modellierung lohnt es sich, das Messprinzip des Geräts zu betrachten. Der Sensor ist im Groben eine Infrarot-Leuchtdiode, die in einem der Spektralbereiche sendet, in denen Kohlendioxid stark absorbiert (weswegen es ja den Treibhauseffekt macht). Das Signal wird dann von einer Fotodiode (oder etwas ähnlichem) aufgefangen, und die Schwächung des Signals ist ein Maß für die Konzentration von CO₂ zwischen LED und Fotodiode.

Allerdings sind alle Halbleiter temperaturempfindlich, und irgendwas, das im Infrarotbereich empfängt, wird schon zwei Mal viel Kalibration brauchen, um Temperatursystematik wegzukriegen. Mit Sicherheit tut die eingebaute Software schon viel in der Richtung. Aber ein Dichteplot zwischen Temperatur und Konzentration zeigt durchaus einen ganzen Haufen Struktur:

Manches davon ist ziemlich sicher Physik, so insbesondere, dass die ganz hohen Konzentrationen bei niedrigen Temperaturen auftreten – das ist ein Jahreszeiteneffekt. Anderes ist ganz klar Instrumentensignatur, am klarsten das abgetrennte Schwanzende jenseits von 50 Grad Sensortemperatur. Offenbar[1] ist die Kalibrationskurve (also: Welches Signal von der Fotodiode soll bei welcher Temperatur welche CO₂-Konzentration ausgeben?) abschnittsweise definiert, und beim Abschnitt über 50 Grad wird sich wohl wer bei den Hundertern vertippt haben. Im Tagesplot entspricht dieses Schwänzchen übrigens dem abgesetzten Häubchen am Nachmittag.

Im Neunmonatsplot zeigen sich die Punkte dort in ein paar der Spitzen zwischen dem 18. Juli und dem 4. August, nur dass sie dort mit den „normalen“ Daten mit einer Linie verbunden sind und nicht als abgesetzt auffallen; Grundregel Nummer 312: Vorsicht beim Verbinden mit Linien. In einem Scatterplot, bei dem Punkte in dem abgetrennten Schwanzende rot gefärbt sind, sind die Unstetigkeiten (und mithin die Fehlkalibration des Geräts) offensichtlich:

In meinem Datensatz betrifft das 1027 von 724'424 Datenpunkten – eigentlich sollte ich die wegwerfen, aber wenn mensch einfach 100 von ihnen abzieht, kommen weitgehend glatte Kurven raus, und so wird das schon nicht völlig unvernünftig sein. Ich bin auch ganz glücklich mit meiner Erklärung des Vertippers bei der Hunderterstelle der abschnittsweise definierten Kalibrationskurve.

Mir gefällt aber auch die offensichtliche Korrelation von Temperatur und CO₂ zwischen 30 und 50 Grad nicht[2], die sich in der „Fahne nach Nordosten“ im T-Kalibrationsplot zeigt. Mein, na ja, Modell würde da keine Korrelation[3], also einen ebenen Verlauf geradeaus „nach Osten“ erwarten lassen.

Soweit das Modell zutrifft, ist die ganze Steigung nur eine weitere Fehlkalibration der Temperaturabhängigkeit der Photodiode. Mithin sollte mensch wahrscheinlich oberhalb von 30 Grad etwas wie (T − 30) ⁄ 20⋅50 ppm (weil: über die 20 Grad oberhalb von 30 Grad Celsius geht die Fahne um rund 50 ppm nach oben) abziehen. Gegen Ende des Posts erwähne ich, warum ich die Korrektur am Ende auf (T − 25) ⁄ 25⋅80 ppm erweitert habe.

Zur näheren Untersuchung habe ich die Punkte aus der Fahne in einer normalen CO₂-Zeitreihe rosa eingefärbt, die Temperatur dazugeplottet und bin dabei zunächst auf ein Ereignis gestoßen, das mir sehr merkwürdig vorkam:

Ich habe immer noch keine Ahnung, was hier passiert ist: Wenn nicht einfach nur das Instrument durchgedreht ist, dann muss von irgendwoher ein Schwung kohlendioxidreiche Luft gekommen sein, die aber an der Temperatur unter der Eisdose nichts geändert hat. Wenn der Spike von einer Wohnungslüftung käme, wäre das sehr seltsam, denn wenn die Leute die Fenster zu hatten – und nur dann hätte sich CO₂ anreichern können – wäre die Luft innen je nach Bauphysik fast sicher kühler oder wärmer gewesen als draußen. Hm. Mein bestes Angebot: Luft ist ein schlechter Wärmeleiter, und welches Lüftchen auch immer hier wehte, hat den sonnenbeschienenen Sensor einfach nicht kühlen können.

Ach so: die Kurve ist beim Ergeignis grün, obwohl sich an der Temperatur nichts geändert hat, weil bei den hohen CO₂-Konzentrationen die ensprechenden Punkte aus der kleinen Nase über der Nordost-Fahne im T-Kalibrationsplot zu liegen kommen. Pink ist aber nur gemalt, was in der Fahne selbst ist.

Fruchtbarer ist die Betrachtung des parallelen Verlaufs von CO₂ und Temperatur zwischen 13 und 16 Uhr. Diese Parallelität besteht tatsächlich nur für die pinken Punkte. Es gibt keine plausible Physik, die CO₂ und Temperatur (auf diesen Skalen) gleich schwingen lassen würde. Wenn ich mit meiner Augenmaß-Kalibrationsfunktion von oben korrigiere, verschwindet dieses Signal auch tatsächlich fast vollständig:

Beachtet die obere Achsenbeschriftung; das ist die Nachkalibration, die ich zunächst angebracht habe; mit meiner verbesserten Nachkalibration ab 25°C bleibt etwas mehr Signal übrig:

Dass die Flanken des „Störsignals“ jetzt steiler sind, finde ich eher beruhigend, denn ich glaube, dass das etwas mit direkter Sonneneinstrahlung zu tun hat (die die Infrarotdiode ganz bestimmt stört), und Licht und Schatten gehen natürlich viel schneller als die Erwärmung von Luft und Gerät. In der Tat würde Streulicht von der Sonne so tun, als käme etwas mehr Licht beim Sensor an, als wäre also etwas weniger CO₂ im Strahlengang. Wenn ihr scharf schaut: im Plot sieht es aus, als sei die CO₂-Schätzung niedriger, wenn die Temperatur ansteigt (also vermutlich die Sonne schien) und höher, wenn sie das nicht tat. Wäre das hier Wissenschaft, müsste ich dieser Spur genauer nachgehen. So, wie es ist, kann …

⎆
My First Libreoffice Macro in Python

2022-12-21
This is what I was after: Immediate Wikipedia search from within Libreoffice. In the document, you can also see the result of a Python dir() as produced by the inspect macro discussed below.
While I still believe the creation „office“ software was one of the more fateful turns in the history of computing and this would be a much better world if there hadn't been VisiCalc and WordStar, not to mention all the software they spun off, I do get asked about Libreoffice quite a bit by people I have helped to get off of Windows.

The other day one of them said: „You know, wouldn't it be nifty if I could mark a term, hit F3 and then that'd do a Wikipedia search for that term?“ I could sympathise with that, since the little one-line CLI I have on my desktop has a function pretty much for this, too. That program, however, probably is too mean for people using Libreoffice. But tomorrow is that guy's birthday, and so I thought: how hard can it be to teach this to Libreoffice directly?

Turns out it's harder (for me) than I thought. Which is why I'm writing this post: perhaps a few people will find it among all the partially outdated or (to me) not terribly helpful material. I think I'd have appreciated a post like this before I started investigating Libreoffice's world of macros.

In particular, I'd have liked the reassuring words: „There's no reason to fiddle with the odd dialect of BASIC that's built in, and there's no reason either to use the odd IDE they have.” The way things were, I did fiddle around with both until I couldn't seem to find a way to open a URL from within StarBasic or whatever that thing is called today. At that point I noticed that all it takes for Python support is the installation of a single small package. In addition, for all I can see the BASIC variant has just about as much relevant documentation as the Python API. So… let's use the latter.
Preparations

(a) To enable Python macros in libreoffice version 7 (as in Debian bullseye), you have to install the libreoffice-script-provider-python package.

(b) The extensions go into a directory deep within your XDG .config. So, create and enter this directory:
```
mkdir ~/.config/libreoffice/4/user/Scripts/python/
cd ~/.config/libreoffice/4/user/Scripts/python/
```
I'm calling this directory the script path below.
Figuring out Libreoffice's API

My main problem with this little project has been that I could not figure out Libreoffice's macro-related documentation. The least confusing material still seems to be maintained by openoffice (rather than libreoffice), and what I ended up doing was using Python introspection to discover attribute names and then entering the more promising ones into the search box of the openoffice wiki. I strongly suspect that's not how it's meant to work. If you know about better ways: please drop me a note and I will do an update here.

But: How do you introspect given these macros do not (easily) have a stdout, and there seems to be no support for the Python debugger either?

Based on an example from openoffice, I figured out that to Libreoffice, macros written in Python are just functions in Python modules in the script path, and that the basic entry point to libreoffice is through a global variable that the libreoffice runtime tricks into the interpreter's namespace, namely XSCRIPTCONTEXT. With this realisation, the example code, and some guessing I came up with this (save into the script path as introspect.py):
```
def introspect():
    desktop = XSCRIPTCONTEXT.getDesktop()
    model = desktop.getCurrentComponent()
    text = getattr(model, "Text", None)
    if not text:
        # We're not in writer
        return

    text.End.String = str(dir(model))
```
If all goes well, this will append a string representation of dir(model) to the end of the document, and in this way you can look at just about any part of the API – perhaps a bit clumsily, but well enough.

But first, run Python itself on your new module to make sure there are no syntax errors:
```
python introspect.py
```
That is important because if Python cannot parse your module, you will not find the function in the next step, which is linking your function to a key stroke.

To do that, in libreoffice, create a new text document and do Tools → Customize from the menu. Then, in Range, you should find LibreOffice Macros, then My Macros, and in there the introspect module you just created. Click it, and you should be able to select introspect (or whatever function name you picked) under Function. Then, select the key F4 in the upper part of this dialog and click Modify[1].

After you did that, you can hit F4, and you will see all attributes that the result of getCurrentComponent has. As I said, pasting some of these attribute names into the search box on openoffice's wiki has helped, and of course you can further introspect the values of all these attributes. Thankfully, libreoffice auto-reloads modules, and so traversing all these various objects in this way is relatively interactive.

I freely admit that I have also used this text.End.String = trick to printf-debug when I did the next steps.
The Wikipedia-Opening Script

These next steps included figuring out the CurrentSelection object and, in particular, resisting the temptation to get its Text attribute (which points to its parent, the whole document). Instead, use the String attribute to retrieve what the user has selected. The rest is standard python fare with a dash of what I suppose is cargo-culting on my end (the supportsService thing seeing whether I can subscript the selection; I lifted that from another example I ran into on the openoffice wiki):
```
import webbrowser
from urllib.parse import quote as urlquote

def search_in_wikipedia():
    """Do a wikipedia search for the current selection"""
    desktop = XSCRIPTCONTEXT.getDesktop()
    model = desktop.getCurrentComponent()
    sel = model.CurrentSelection
    if sel.supportsService("com.sun.star.text.TextRanges"):
        selected = sel[0].String.strip()
        if selected:
            webbrowser.open_new_tab("https://de.wikipedia.org/w/index.php"
                "?fulltext=Suchen&search="+urlquote(selected))
```
For all I can see, the Wikipedia search URI is the same across the instances modulo the authority part – so, replace the de in the URL with the code for whatever language you (or persons you need a birthday present for) prefer. Then, save it to the script path and bind it to a function key as before.

I'm sure this can (and should) be made a whole lot more robust with a bit more actual Libreoffice-fu. But it seems to work nicely enough. And me, I have a very personal (in a somewhat twisted sense) birthday present.

[1] I always find that describing operations in GUIs tends to sound like incomprehensible origami instructions. Is there a term for that kind of language already? Gooeynese, perhaps?
⎆
A QR Code Scanner for the Desktop

2022-12-01

qropen.py in action: Here, it has scanned a QR code on a chocolate wrapper and asks for confirmation that you actually want to open the URI contained (of course, it doesn't come with the glitzy background).

When I was investigating the SARS-2 vaccination certificates last year (sorry: that post is in German), I played a bit with QR codes. A cheap by-product of this was a little Python program scanning QR codes (and other symbologies). I cannot say I am a huge fan of these things – I'd take short-ish URIs without cat-on-the-keyboard strings like “?id=508“ instead any day –, but sometimes I get curious, and then this thing comes in handy given that my telephone cannot deal with QR codes.

Yesterday, I have put that little one-file script, qropen.py, on codeberg, complemented by a readme that points to the various ways you could hack the code. I'll gratefully accept merge requests, perhaps regarding a useful CLI for selecting cameras – you see, with an external camera, something like this thing starts being actually useful, as when I used it to cobble together a vaccination certificate checker in such a setup. Or perhaps doing something smart with EAN codes parsed (so far, they just end up on stdout)?

On that last point, I will admit that with the camera on my Thinkpad X240, most product EAN codes do not scan well. The underlying reason has been the biggest challenge with this program even for QR codes: Laptop cameras generally have a wide field of view with a fixed focus.

The wide field of view means that you have to bring the barcodes pretty close to the camera in order to have the features be something like three pixels wide (which is what zbar is most fond of). At that small distance, the fixed focus means that the object is severely out of focus and hence the edges are so blurry that zbar again does not like them.

Qropen.py tries to mitigate that by unsharp masking and potentially steep gammas. But when the lines are very thin to begin with – as with EAN stripes –, that does not really help. Which means that the QR codes, perhaps somewhat surprisingly given their higher information content, in general work a lot better for qropen.py than to the simple and ancient EAN codes.

There clearly is a moral to this part of the story. I'm just not sure which (beyond the triviality that EANs were invented for lasers rather than cameras).

⎆
34 Monate Corona im Film

2022-10-27

Im Sommer 2021 habe ich in zwei Posts Filme zur Visualisierung der Corona-Inzidenzen und einer Art Alters-Scores aus den RKI-Daten der vorherigen anderthalb Jahren vorgestellt und ein wenig das Python-Programm diskutiert, das die generiert.

Ich wollte das über den Sommer immer mal aktualisieren. Aber die Pandemie hat keine erkennbare Pause eingelegt, die eine gute Gelegenheit für eine Art Rückblick gewesen wäre. Jetzt aber ist wohl so in etwa der letzte Moment, in dem so ein Film noch nicht vollständige Clownerei ist, denn schon jetzt testet zumindest anekdotisch kaum mehr jemand. Wenn aber fast niemand PCR-Tests machen lässt, werden die Zahlen, die ich da visualisiere, weitgehend bedeutunglos (weil nicht small data).

Bevor das Ganze zu einer Art aufwändigen Lavalampen-Simulation wird, habe ich die Programme genommen, dem Inzidenzfilm eine dynamische Colorbar gegönnt – angesichts von Inzidenzen von örtlich über 4000/100'000 im letzten Winter war das bisher feste Maximum von 250 nicht mehr sinnvoll – und habe die Filmchen nochmal mit aktuellen Daten gerendert.

Hier also die Inzidenzen per Kreis (Hinweise dazu von vor einem Jahr):

und hier die Alters-Scores (auch dazu Hinweise vom letzten Mal):

Beide Filme sind mit -i 7 gebaut, jeder Tag ist also sieben Frames lang. Inzwischen wären die Filme zwar auch ohne Interpolation lang genug, aber sie hilft wenigstens dem Lavalampen-Effekt – und natürlich verletze ich gerne das Durchhörbarkeits-Zeitlimit von drei Minuten. Auch wenn die Filme gar keinen Ton machen.

Ich habe die Gelegenheit genutzt, um den Code, der das macht – statt ihn nur hier zu verlinken, wo ihn niemand finden wird –, brav auf dem Codeberg abzuladen. Vielleicht hilft er dort ja Leuten, die irgendwann irgendwelche Daten in Kreispolygonen von matplotlib aus plotten wollen. Wenn es dafür bessere Standardverfahren geben sollte, habe ich die zumindest vor einem Jahr nicht gefunden.

⎆
Blog Extensions on Codeberg

2022-10-07
This post takes an odd bend to become an apology for CGI (as in common gateway interface) scripts. This is the netsurf browser communicating with the CGI shell script at the foot of this post.

I have written a few plugins and extensions for this blog, and I have discussed them in a few posts (e.g., feedback form, tag explanations, cited-by links, or the search engine). The code implementing these things has been strewn across the various posts. I have to admit that having that code attached to just a few blog posts has always felt somewhat too early-90iesy to me.

Now that I have created my Codeberg account, I have finally copied together all the various bits and pieces to create a repository on Codeberg that you are welcome to clone if you're running pelican or perhaps some other static blog engine. And of course I appreciate merge requests with improvements.

There is one major file in there I have not previously discussed here: cgiserver.py. You see, I'm a big fan of CGI scripts. They're reasonably simple to write, trivial to deploy, and I have CGIs that have been working with minimal maintenance for more than 20 years. Sure, pulling up an interpreter for every request is not terribly efficient, but for your average CGI that is perhaps called a dozen times per day (depending on how many web crawlers find it interesting) this really doesn't matter. And that's why both the feedback script and the search engine are written as CGIs.

However, in contrast to apache, nginx (which serves this blog) does not support CGI scripts. I even by and large agree with their rationale for that design decision. Still, I would like to run CGIs, and that's why I've written the cgiserver. It is based on Python's built-in HTTP server and certainly will not scale – but for your average blog (or similar site) it should be just fine. And I don't think it has any glaring security problems (that you don't introduce with your CGIs, that is).

Installation is almost trivial: put the file somewhere (the in-source sysvinit script assumes /var/www/blog-media/cgiserver.py, but there's absolutely no magic about this), and then run it with a port number (it will only bind to localhost; the default in the sysvinit script is 6070) and a directory into which you put your CGI scripts (the sysvinit script assumes /var/www/blog-media/cgi).

When you have a cgi script foo, you can dump it in this directory, make it executable and then run it by retrieving http://localhost:6070/foo. In case you have nothing else, you can try a shell script like:
```
#!/bin/sh
echo "content-type: text/plain"
echo
/usr/games/fortune
```
(which of course only works in this form if you have something like fortunes-en installed on a Debian box). That should be enough to give you something like the screenshot opening this post. Even more than 25 years after I have written my first CGI, I am still amazed how simple this is.

Disclaimer: Writing CGI scripts that take input such that they are not trivially exploitable is higher art. So… don't do it, except as a game. Oh, and to debug your scripts, simply let cgiserver run in a terminal – that way, you will see what your scripts emit on stderr. Note, however, that the way the sysvinit script starts cgiserver, it will run as nobody; if things work when you start cgiserver yourself but not when it's running as a daemon, that's the most likely reason.
⎆
Bahnauskuft auf antiken Geräten – und auf Codeberg

2022-09-04

Bahnauskunft von 2022 auf einem Nokia N900 von 2009: Es braucht inzwischen etwas Mühe, um das gebastelt zu kriegen.

Als die Bahn-Webseite nicht mehr ordentlich auf kompakten Browsern wie dillo funktionierte und auch nicht per WAP– also Mitte der 2010er Jahre –, habe ich mir ein ein kleines Skript geschrieben, das die wesentlichen Infos zur Zugauskunft aus dem HTML herausklaubte und dann in einem einfachen Kommandozeilen-Interface darstellte. Das war, worum es im letzten Sommer bei meinem Rant gegen Zwangs-Redirects umittelbar ging. Der weitere Hintergrund: Ich will Zugauskünfte von meinem alten Nokia N900 aus bekommen (und im Übrigen seit der Abschaltung von UMTS auch über eine 2G-Funkverbindung, also etwas wie 10 kB/s)[1].

Nachdem das – jedenfalls nach Maßstäben von Programmen, die HTML auf Webseiten zerpflücken – überraschend lang gut ging, ist das im Rahmen der derzeitigen Verschlimmbesserung der Bahn-Seite neulich kaputt gegangen. Obendrauf ist die Javascript-Soße auf bahn.de damit so undurchsichtig geworden, dass mich die Lust, das Skript zu pflegen, sehr nachhaltig verlassen hat. In dieser Lage kam ein Vortrag über die Bahn-APIs, den jemand bei der Gulasch-Programmiernacht 2019 gehalten hat, gerade recht. Also: Das Video davon.

In diesem Video habe ich gelernt, dass mein „unpromising“ im Rant vor einem Jahr,

I know bahn.de has a proper API, too, and I'm sure it would be a lot faster if I used it, but alas, my experiments with it were unpromising [...],

einen tiefen Hintergrund hat. Die Bahn hat nämlich keine API für die Fahrplanauskunft.

Was es aber stattdessen gibt: die HaFAS-API, auf die die Reiseplanung der Bahn-App selbst aufsetzt. Und es stellt sich heraus, dass Leute schon mit viel Fleiß ausbaldowert haben, wie die so funktioniert, etwa in pyhafas.

Mit pyhafas kann ich all das schreckliche HTML-parsing aus dem alten bahnconn.py durch ein paar Aufrufe in pyhafas rein ersetzen. Aber leider: pyhafas ist echt modernes Python, und weil es viel mehr kann als es für bahnconn.py bräuchte, wäre das Rückportieren davon nach Python 2.5 ein ernsthaftes Projekt; mehr habe ich aber auf meinem N900 nicht. Außerdem bin ich bekennender Fan von ein-Modul-und-stdlib-Programmen: die brauchen keine Installation und laufen zudem mit allem, das irgendwie Python verdauen kann, also etwa auch jython oder sowas, was spätestens dann in Frage steht, wenn Abhängigkeiten C-Code enthalten.

Deshalb habe ich aus pyhafas die Dinge, die bahnconn dringend braucht, abgeschaut und eine minimale, Python-2.5-verträgliche Implementation gebastelt. Das Ergebnis: ein neues bahnconn. Holt es euch, wenn ihr Bahnauskunft auf älteren Geräten haben wollt. Ich habe es jetzt nicht auf Atari TTs probiert, aber ich kann mir gut vorstellen, dass es selbst da noch benutzbar ist.

Codeberg

Gerade, als ich den Code einfach wieder hier auf dem Blog abwerfen wollte, habe ich beschlossen, das könne ein guter Anlass sein, endlich mal einen zweiten Blick auf Codeberg zu werfen.

Bisher habe ich nämlich für allen etwas langlebigeren oder größeren Code (also: nicht einfach nur am Blog abgeworfenen Kram), ganz DIY, ein eigenes Subversion-Repository betrieben. Was in den letzten Jahren neu dazukam, habe ich in git+ssh+cgit gesteckt.

Natürlich hat das niemand mehr gesehen; nicht mal Suchmaschinen gucken mehr auf sowas, seit aller Code der Welt bei github landet. Deshalb, und auch, weil ich Monstren wie gitea und gitlab wirklich nicht auf meiner Maschine haben will (allerdings: cgit ist ok und würde für Publikation auf subversion-Niveau reichen), habe ich mich mit dem Gedanken, dass mein Kram auf einer öffentlichen Plattform besser aufgehoben sein mag, mehr oder minder abgefunden.

Auf Github bin ich beruflich schon so viel zu viel unterwegs, und der Laden ist deutlich zu nah am Surveillance Capitalism. Zwar kenne ich hinreichend Projekte und Firmen, die ihnen Geld geben, so dass sie gewiss ein konventionell-kapitalistisches Geschäftsmodell fahren könnten; aber schon da fehlt mir der Glaube. Obendrauf hat mir Microsoft in meinem Leben schon so viel Kummer bereitet, dass ich ihnen (bzw. ihrem Tochterunternehmen) nicht noch mehr KundInnen zutreiben will.

Codeberg, auf der anderen Seite, wird von einem Verein betrieben und macht generell vieles richtig, bis hin zu Einblendungen von Javascript-Exceptions (warum machen das eigentlich nicht alle?), so dass die Seite nicht einfach heimlich kaputt ist, wenn ich Local Storage verbiete (gitea, die Software, auf der Codeberg soweit ich sehe aufsetzt, kann leider immer noch nicht ohne).

Von dem gitea-Krampf abgesehen hat gestern alles schön funktioniert, nichts an der Anmeldeprozedur war fies oder unzumutbar. Codeberg hat hiermit erstmal das Anselm Flügel Seal of Approval. Ich denke, da werde ich noch mehr Code hinschaffen. Und mal ernsthaft über Spenden nachdenken.

[1] Janaja, und natürlich nervte mich die fette Bahn-Webseite mit all dem Unsinn darauf auch auf dem Desktop und auch schon vor der gegenwärtigen Verschlimmbesserung.

⎆
Quick RST Previews for Posts in Pelican

2022-05-26
In January, I described how I use this blog's engine, pelican, and how I have a “development” and a “production” site (where I will concede any time that it's exceedingly silly to talk about “production” in this context). Part of that was a trivial script, remake.sh, that I would run while writing and revising a post to format it without doing too much unnecessary work. This script was running between a couple and a couple of dozen times until I was happy with an article.

What the script did was call pelican asking to only write the document being processed. When pelican was instructed to cache work on the other articles, that was enough to keep build times around a second on my box; but as the number of posts on this blog approaches 200, build times ended up on the totally wrong side of that second, and I thought: “Well, why don't I run, perhaps, rst2html for formatting while revising?” That would be, essentially, instantaneous.

But pelican does a lot more than rst2html. Especially, having the plugins and the templating available is a good thing when inspecting a post. So, I got to work and figured out how pelican builds a document. The result is a script build-one that only looks at a single (ReStructuredText) article – which it gets from its command line – and ignores everything else.

This is fast enough to be run whenever I save the current file. Therefore, in my pelican directory I now have, together with the script, the following .vimrc enabling just that (% expands to the file currently edited in vim):
```
augroup local
  au!
  autocmd BufWritePost *.rst !python build-one %
augroup END
```
I've briefly considered whether I should also add some trick to automatically reload a browser window when saving but then figured that's probably overdoing things: In all likelihood I want to scroll around in the rendered document, and hence I will have to focus it anyway. If I do that, then effort spent on saving pressing r after focusing feels misplaced.

The script does have an actual drawback, though: Since pelican does not get to scan the file system with build-one, it cannot do file name substitution (as in {filename}2022-05-26.rst) and will instead warn whenever seeing one of these. Since, as described in January, my static files are not managed by pelican, that is not a serious problem in my setup, except I have to watch out for broken substitutions when doing a final make html (or the make install).
Insights into Pelican

It took me a bit to figure out how the various parts of pelican fit together at least to the extent of letting me format a ReStructuredText document with the jinja templates. Let me therefore briefly discuss what the script does.

First, to make pelican do anything remotely resembling what it will do on make html, you have to load its settings; since I assume I am running in pelican's directory and this is building a “draft” version, I can simply do:
```
settings = pelican.read_settings("pelicanconf.py")
```
With that, I already now where to write to, which lets me construct a writer object; that will later arrange for actually placing the files. I can also construct a reader for my ReStructuredText files (and you would have to change that if you are writing in Markdown); these readers decouple the Article class from input formats:
```
writer = writers.Writer(settings["OUTPUT_PATH"], settings)
reader = readers.RstReader(settings)
```
With that, I have to delve deep into pelican's transformation machinery, which consists of various generators – for articles, static files, pages, whatever. The constructors of these generator classes (which are totally unrelated to Python generators) take a lot of arguments, and I cannot say I investigated why they insist on having them passed in when I fill them with data from settings anyway (as does pelican itself); but then I suspect these extra arguments are important for non-Article generators. I only need to generate a single article, and so stereotypically writing:
```
artgen = generators.ArticlesGenerator(
    settings.copy(), settings,
    settings["PATH"], settings["THEME"], settings["OUTPUT_PATH"])
```
does the trick for me.

Article generators will usually collect the articles to generate by looking at the file system. I don't want that; instead, I want to construct an Article instance myself and then restrict the generator's action to that.

The Article class needs to be constructed with content and metadata, which happen to be what readers return. So, to construct an Article from the RST file passed in in source_path, I need to say:
```
content, metadata = reader.read(source_path)
art = contents.Article(content, metadata,
    source_path=source_path, settings=settings)
```
After all that preparation, all that is left to do is overwrite any misguided ideas the article generator might have on what I would like to have processed and then let it run:
```
artgen.translations = []
artgen.articles = [art]
artgen.generate_articles(
    functools.partial(writer.write_file, relative_urls=True))
```
(the currying of the writer's write_file method to make sure it creates relative URLs you can probably do without, but I'm a fan of relative URLs and of almost anything in functools).
⎆
Globaldokumente in Libreoffice zusammenführen

2022-04-16
Ich hatte gerade einer armen Seele zu helfen, die drei Bücher mit Microsoft Word für Windows geschrieben hat und dazu Word-Zentraldokumente verwendet hat, vor allem wohl, weil vor 20 Jahren – als diese Projekte starteten – Word nicht mit mehreren hundert Seiten auf einmal umgehen konnte. Unbenommen, dass mensch mit Office-Software keine richtige Arbeit machen sollte: da die arme Seele auf Linux migrierte, musste der ganze Kram auf Libreoffice. Für eine Migration auf vernünftige Technologien (TeX, ReStructuredText, Docbook oder was auch immer) reichten weder meine Geduld noch der Erlösungswille der armen Seele.

Erste Ernüchterung: Word-Zentraldokumente (wie auch Libreoffice-Globaldokumente) speichern ihre Pfade absolut: sie können also ohne Tricks nicht bewegt werden, insbesondere nicht von einem Windows-Dateisystem in ein Linux-Dateisystem. Wie so oft stellt sich die Frage, warum das mal ein wie ein gutes Design ausgesehen haben mag.

Nach einem kurzen Blick auf das Arbeiten mit Globaldokumenten habe ich beschlossen, dass Rechner jetzt wirklich groß genug sind, um 500-seitige Dokumente im Speicher zu halten und dass der ganze Zentral- und Globaldokumentenzauber mit den begleitenden Komplikationen nicht mehr sein muss.

Nur: Wie kopiert mensch 20, 30 solche Dateien zusammen? Was ich dazu in der verbleibenden Libreoffice-Doku (das Wiki, das das jetzt sein soll, überzeugt mich übrigens nicht) und dem weiteren Netz gefunden habe, fand ich eher… unbefriedigend – Erinnerungen an Windows-Foren („um den Sound-Treiber zu reparieren, musst du den CD-Treiber deinstallieren. Hat bei mir funktioniert“) werden da dann und wann schon wach. Deshalb dachte ich mir, es könnte nützlich sein, wenn ich auch ein paar Rezepte beitrage, auch wenn ich – Disclaimer – selbst keine Office-Software verwende und in dem Sinn selbst höchstens einäugig bin.
Wie machte ich aus einem Satz von ODTs oder DOCs ein einziges ODT?
1. Ein Verzeichnis anlegen und alle Filialdateien reinkopieren.
2. Die Dateien so benennen, dass eine einfache Sortierung sie in die richtige Reihenfolge bringt (ich habe einfach 000_, 001_ usf vor die Namen gesezt).
3. Libreoffice starten, Neu → Globaldokument.
4. F5 drücken, um in den Navigator zu kommen, dort aufs Einfügen-Icon klicken; das poppt eine Auwahlbox auf.
5. In dieser Auswahlbox in das Verzeichnis mit den Filialdateien gehen und diese in die Reihenfolge sortieren lassen, in der sie nachher im Dokument erscheinen sollen.
6. Alle Dateien im Verzeichnis auswählen, z.B. durch Control-A oder per Shift-Klick, dann den Import bestätigen.
7. Datei → Exportieren, dabei ODT als Zielformat auswählen; das ist der erste Schritt, um von der Einbettung im Globaldokument wegzukommen. Ich nenne diese Datei jetzt mal joined.odt.
8. Das so erzeugte ODT ist leider überall schreibgeschützt, und ich habe keinen Weg gefunden, diesen Schreibschutz per Klicken wegzuzaubern, bevor ich die Geduld mit Doku, Menüs und vor allem Foren verloren habe und mit epubedit beigegangen bin (vgl. unten). Mit dem kleinen Skript dort unten könnt ihr Folgendes in einer Shell laufen lassen:
  
  epubedit joined.odt sed -ie 's/text:protected="[^"]*"//g' content.xml
  
  (ihr könnt natürlich auch mit einem Editor oder gar mit dem hervorragenden xmlstarlet die ganzen text:protected-Attribute löschen)[1]. Geht dann aus der Shell vom epubedit wieder raus; das schreibt joined.odt neu.
9. Das neue joined.odt in libreoffice öffnen.
10. Bearbeiten → Verknüpfungen, wieder alle Auswählen (^A), und dann den Lösen-Knopf drücken.
Das Ergebnis dieser Prozedur ist ein zusammenhängendes „Dokument“ (wenn mensch keine großen Ansprüche an Dokumente hat).

Zumindest in meinem Fall fing damit die Arbeit allerdings erst an, weil jedes Filialdokument eigene und verrückte Absatzvorlagen hatte. Ich schreibe gleich, wie ich das aufgeräumt habe, aber zunächst müsen wir geschwind über das erwähnte epubedit reden.
epubedit

Was Open Document-Dateien tatsächlich etwas angenehmer im Umgang macht als einige andere Office-Dateien, die ich hier erwähnen könnte, ist, dass sie eigentlich nur zip-Archive sind, in denen von für sich nicht unvernünftigen Standards (z.B. XML und CSS) beschriebene Textdateien leben. Das hat mir beispielsweise den obigen Trick mit der Ersetzung der text:protected-Attribute erlaubt.

Diese Architektur haben sie mit Ebooks im epub-Format gemein, und um an denen geschwind mal kleine Korrekturen vorzunehmen, habe ich mir vor Jahren ein kleines Shell-Skript geschrieben:
```
#!/bin/bash

if [ $# -ne 1 ]; then
        echo "Usage: $0 <epub> -- unpack an epub, open a shell, pack it again."
        exit 0
fi

workdir=$(mktemp -d /tmp/workXXXXXX)

cleanup() {
        rm -rf $workdir
}
trap cleanup EXIT
if [ ! -f  "$1".bak ]; then
        cp -a "$1" "$1".bak
fi

unzip "$1" -d $workdir
(cd $workdir; bash)
fullpath=$(pwd)/"$1"

cd $workdir
zip -r "$fullpath" *
```
Nehmt das und legt es als – sagen wir – epubedit irgendwo in euren Pfad und macht es ausführbar. Ihr könnt dann für irgendein epub oder odt epubedit datei.odt sagen und landet in einer Shell, die im Wurzelverzeichnis des jeweiligen ZIP-Archivs läuft. Dort könnt ihr nach Herzenslust editieren – bei ODTs ist der Inhalt in content.xml –, und wenn ihr fertig seid, beendet ihr die Shell und habt ein entsprechend verändertes ODT oder epub.

Weil dabei gerne mal was schief geht, legt das Skript ein Backup der Originaldatei an (es sei denn, es gäbe schon so ein Backup; die Erfahrung zeigt, dass mensch in der Regel lieber das ursprüngliche Backup behalten will…).
Stilfragen

Nun ist das vereinte Dokument zwar immerhin nur noch eine einzige Datei, die zudem – wow! – auch bewegt werden kann. Zumindest mit der Genese in meinem Fall, also den vielen Einzel-Word-Dateien, ist sie trotzdem kaum brauchbar, weil Word einige hundert Formatvorlagen erzeugt hat, meist mit so nützlichen Namen wie Formatvorlage_20_16_20_pt_20_Block_20_Erste_20_Zeile_3a__20__20_05_20_cm_20_Zchn oder Fußnotentext1_20_Zchn oder auch apple-converted-space. Dieses Problem ist schlimm, und ich habe schließlich eingesehen, dass es ohne ein kleines Programm und einige Handarbeit nicht lösbar ist.

Das Programm hat am Anfang nur Stilnamen aus dem Dokument rausgeprökelt und auf die Standardausgabe gelegt. Inzwischen ist das zu einer Basis für eine Abbildungsdatei geworden, und auch für die Abbildung als solche haben reguläre Ausdrücke noch gereicht. Wegen der Abhängigkeiten der Stile untereinander blieb jedoch immer noch jede Menge Mist in der Liste der verwendeten Stile zurück. Deshalb musste ich schließlich doch noch ordentliche XML-Verarbeitung anwerfen, um die styles.xml umzufummeln. Das Ergebnis ist das Programm defuse-libreoffice-style.py. Wenn ihr dieses Programm für die Dauer der Verarbeitung in euer Homeverzeichnis legt, würdet ihr die Stile wie folgt vereinheitlichen:
1. epubedit joined.odt; alles Weitere passiert in der Shell, die das öffnet.
2. python3 ~/defuse_libreoffice-style.py > ~/style-map.txt – wenn ihr das Skript nicht in eurem Home lagert, müsst ihr diesen Pfad anpassen. Und ich lege die Stil-Abbildung ins Home und nicht ins aktuelle Verzeichnis, damit die Abbildung (die recht viel Arbeit ist) nicht gleich verloren ist, wenn ihr die Shell verlasst. Ich jedenfalls habe besonders beim ersten Mal ein paar Anläufe gebraucht, bis das Mapping gut gepasst hat.
3. Editiert die Datei ~/style-map.txt mit einem Texteditor (also auf keinen Fall mit libreoffice selbst). Da drin stehen Zeilen wie:
  
  Footnote_20_Symbol -> Footnote_20_Symbol
  
  – in meinem Fall ungefähr 200 davon. Die Aufgabe ist jetzt, die rechten Seiten dieser Zeilen auf eine Handvoll Stile runterzubringen (Textkörper, Überschrift_1, Überschrift_2, Zitat, Fußnotenzeichen und Fußnote waren mein Minimum); die Zeile oben habe ich zum Beispiel zu:
  
  Footnote_20_Symbol -> Fußnotenzeichen
  
  gemacht. Es ist nicht immer einfach, herauszukriegen, was wohl eine Vorlage mal tun sollte; meist hat Word aber doch einen gewissen Hinweis darauf im Namen hinterlassen.
4. Wenn die Abbildung fertig ist, lasst das Python-Skript nochmal laufen; wenn es nämlich ein Argument bekommt, interpretiert es das als Abbildung und passt sowohl content.xml als auch style.xml entsprechend an:
  
  python3 ~/defuse_libreoffice-style.py ~/style-map.txt
5. Um zu sehen, welche Stile noch übrig sind, könnt ihr das Skript ein weiteres Mal ohne Argumente laufen lassen; das gibt dann die noch vorhandenen Stile ins Terminal aus:
  
  python3 ~/defuse_libreoffice-style.py
  
  Wenn noch was dabei ist, das nicht übrig bleiben soll, könnt ihr style-map.txt anpassen und Schritt (4) nochmal laufen lassen (oder nochmal vom Backup des ODT anfangen).
6. Verlasst zum Abschluss die Shell vom epubedit und guckt im libreoffice nach, ob alles geklappt hat. Libreoffice erzählt wahrscheinlich, dass das Dokument beschädigt sei (aber nicht genauer, was eigentlich; hier rächt sich, dass ich die Open Document-Standards nicht gelesen und stattdessen einfach munter drauflosgehackt habe). Das, was es zur Reparatur unternimmt, hat aber bei mir immer gut funktioniert – insofern: Nur Mut.
Und für den Fall, dass jemand in den Python-Code reinguckt: Nein, auch wenn der StyleSanitiser immerhin ordentlich XML bearbeitet (im Gegensatz zu dem RE-Hacks von oben), ist das immer noch nicht Open Document-allgemein, denn ich habe die spezifische Wahl des text:-Präfix von Libreoffice darin hart kodiert, was sich für „richtige“ Software nicht gehören würde. Aber SAX mit richtigen Namespaces macht keinen Spaß, und ich rechne erstmal nicht damit, dass dieser Code je mit ODTs laufen wird, die nicht von Libreoffice kommen.
Und die Stichworte?

Die Bücher hatten auch je ein Stichwortverzeichnis. Bei einem Dokument hat das gut funktioniert, bei den anderen standen im Verzeichnis ein paar ordentliche Begriffe, ein paar Begriffe mit sinnlosen typografischen Anführungszeichen und ganz viele Einträge für das leere Wort. Ich habe keine Ahnung, wie es dazu kam.

Bei der Reparatur hilft erneut der Umstand, dass ODT im Kern ein nicht ganz unvernünftiges XML ist. Dabei sieht das Markup für ein Stichwort beispielsweise so aus:
```
Karl Valentin …
```
⎆
Von Geburtstagen und /etc/papersize

2022-04-11
Ich bin ein stiller Fan des Debian-Pakets installation-birthday. Das hat mir vorhin eine Mail geschrieben:
```
Date: Mon, 11 Apr 2022 11:04:11 +0200
From: Anacron <root@hostname-withheld>
To: root@hostname-withheld
Subject: Anacron job 'cron.daily' on hostname-withheld

/etc/cron.daily/installation-birthday:

                  0   0
                  |   |
              ____|___|____
           0  |~ ~ ~ ~ ~ ~|   0
           |  |           |   |
        ___|__|___________|___|__
        |/\/\/\/\/\/\/\/\/\/\/\/|
    0   |       H a p p y       |   0
    |   |/\/\/\/\/\/\/\/\/\/\/\/|   |
   _|___|_______________________|___|__
  |/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/|
  |                                   |
  |         B i r t h d a y! ! !      |
  | ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ |
  |___________________________________|


Congratulations, your Debian system "hostname-withheld" was installed
15 year(s) ago today!


Best wishes,

Your local system administrator
```
Oh wow. So ein glatter Geburtstag ist doch eigentlich ein Grund zu feiern.

Er ist aber auch ein Grund zum Nachdenken. Irgendwas kann da nicht stimmen. Bevor meine Haupt-Arbeitsmaschine (deren Filesystem so etwa aus dem Jahr 1996 stammt) ein Debian-System wurde, habe ich das selbst gestrickt und gebaut. Das habe ich mir 2007 ganz sicher nicht mehr angetan, denn die Zahl der auf einem durchschnittlichen Desktop-Linux verbauten Zeilen und irgendwie zusammenzufummelnden Komponenten war schon ein paar Jahre vorher deutlich jenseits meiner persönlichen Schmerzgrenze. Was aber ist dann 2007 passiert?

Nun: wozu gibt es die Quellen? Im Fall von installation-birthday führt ein schneller Blick in die cron-Datei, die in der Mail erwähnt ist, auf ein Python-Skript, das das Installations-Datum berechnet als:
```
for candidate in self.gen_installation_datetimes():
  # Use the oldest mtime we can find
  if dt is None or candidate < dt:
    self.log.debug("Preferring %s over %s", candidate, dt)
    dt = candidate
```
gen_installation_datetime nun geht durch die Timestamps des Root-Filesystems, von /var/log/installer (gibts bei mir nicht), /var/log/bootstrap.log (auch nicht), /var/lib/vim (da kommt das Datum her und reflektiert so wohl eher irgendein Datum in der Entwicklung des Pakets, wahrscheinlich die Einführung formaler Addons), /root und /etc/machine-id (diese Datei ist von 2014, und eine kurze Recherche zeigt: Da habe ich zum ersten Mal versucht, mich mit systemd anzufreunden).

So romatisch es ist, mein Systemalter anhand von vi-Hilfsdateien zu schätzen: Das will ich schon anders haben. Wie aber bekomme ich das Migrationsdatum raus?

In meinem Dateisystem sind noch viele Dateien, die auf das selbstgestrickte Linux zurückgehen (und weitere, die von noch älteren Maschinen stammen, aber die beschränken sich auf /home), ich kann also nicht einfach nur die älteste Datei suchen. Es gibt aber eine interessante Häufung von Timestamps in /etc/ am 7. Juli 2005, darunter adduser.conf und papersize, die sehr nach etwas klingen, das Debian beim Überbügeln des Selbstbausystems angelegt haben könnte.

Umgekehrt gibts davor eigentlich nichts, das irgendwie nach Debian aussieht; der Timestamp von /etc/papersize wird es dann wohl sein.

Jetzt wollte ich es genauer wissen, und weil ich zu faul bin, aus der ls-Manpage oder sonstwoher rauszupopeln, wie ich ohne Code an die Uhrzeit der letzten Änderung herankomme, muss es ein Python-Einzeiler tun (ja, sowas gibts):
```
$ python -c "import os, datetime; print(datetime.datetime.utcfromtimestamp(os.path.getmtime('papersize')))
2005-07-07 12:48:37
```
Also: am Donnerstag, 7.7.2005, so gegen drei am Nachmittag, habe ich offenbar meine Arbeitsmaschine auf Debian umgestelt. Uiuiui. Dass ich das Selbstbauen doch so lange durchgehalten habe, hatte ich gar nicht mehr in Erinnerung.

Damit in Zukunft installation-birthday das richtige Datum nimmt, habe ich gerade noch:
```
$ sudo touch -r /etc/papersize /var/log/installer
```
gesagt. Damit gibt es ein /var/log/installer mit dem Timestamp von /etc/papersize, das jetzt installation-birthday das richtige Signal gibt:
```
$ installation-birthday --verbosity 1
I: Installation date: 2005-07-07
```
Ich freue mich schon auf den 7.7. Das ist übrigens in diesem Jahr, wie schon 2005, wieder ein Donnerstag.
⎆
A Feedback Form in Pelican

2022-04-09
I realise that the great days of discussions on blogs are over, as Sam Hartman blogged the other day – at least for now. Still, I'd like to make it somewhat more straightforward to send me feedback on the posts here than having to get the contact address and dropping me a mail. Hence, I've written a little Python script, feedback, that lets people comment from within their web browsers.

Nachtrag (2022-10-07)

Don't take it from here; rather, see https://codeberg.org/AnselmF/pelican-ext

While the script itself is perfectly general and ought to work with any static blog engine, the form template I give in the module docstring is geared towards pelican and jinja, although only in very few places.

To make it work, this needs to become a CGI (the template assumes it will show up in /bin/feedback according to the server configuration). The notes on deployment from my post on the search engine apply here, too, except that in addition the host has to be able to deliver mail. Most Unix boxes do locally, but whether anyone reads such mail is a different question.

Is it ethical to check “ok to publish” by default?

To configure where it sends mail to (by default, that's root, which may make sense if you have your own VM), you can set the CONTACT_ADDRESS environment variable (see the search engine post in case you're unsure how to do that for a web context). If your machine is set up to deliver mail to remote addresses – be it with a full mail server or using a package like nullmailer –, you could use your “normal” mail address here. In that case, you probably should inform people in your privacy policy that their comments will be sent by unencrypted mail, in particular if that “normal“ e-mail is handled by one of the usual rogues (Google still gets about a half of the mail I send – sigh).

If you look below (or perhaps right if you run your browser full-screen), you will see that there is a checkbox “feel free to publish“ that, right now, I have checked by default. I had some doubts about that in terms of creepy antipatterns. Of course I am as annoyed by most contemporary cookie banners as anyone else, which in violation of the GDPR usually have practical defaults – sure: not what you get when you say “configure” – set at the maximum creep level the operators believe they can get away with. On the other hand, defaults should also be expectable, and I'd guess the expectation when someone fills out a reply form on a blog is that the response will be published with the article. If you disagree: well, the comment form is there for you.

In terms of spam protection, I do an incredibly dumb textcha. Even if this script got deployed to a few dozen sites (which of course would be charming), I cannot see some spam engine bothering to figure it out; since it just sends a mail to the operator, there is basically nothing to be gained from spamming using the CGI. I fully expect this will be enough to keep out the dumb spambots that blindly send whatever forms they can find – it has worked on many similar services.
Security Considerations

The feedback script does at least two things that could be exploited:
1. It enters remotely controlled values into rendered HTML and
2. It calls a local binary with content controlled by the remote user.
In case (1) – that's when I put the URI of the originating article into the reply message to send people back to where they came from –, this can potentially be exploited in cross-site attacks. Suppose you trust my site on only execute benign javascript (I give you that's close to an oxymoron these days), someone could trick you into clicking on a link that goes to my site but executes their, presumably adversarial, javascript.

Bugs aside, the script is resilient against that, as it properly escapes any user input that gets copied into the output. That is thanks to my “micro templating“ that I keep around to paste into such one-script wonders. Have a look into the code if you're interested in how that works. And totally feel free to paste that into any Python code producing HTML or XML templated in any way – sure, it's not jinja or stan, but it has covered 80% of my templating needs at much less than 20% of the effort (counted in code lines of whatever dependency you'd pull in otherwise), which is a good deal in my book.

Case (2) is probably a lot more interesting. In the evaluate_form function, I am doing:
```
mail_text = MAIL_TEMPLATE.format(**locals())
```
Code like this usually is reason for alarm, as far too many text formats can be used to execute code or cause other havoc – the cross-site thing I've discussed for HTML above being one example, the totally bizarre Excel CSV import exploit another (where I really cannot see how this doesn't immediately bring every Windows machine on this planet to a grinding halt). In this case, people could for example insert \ncc: victim@address into anything that gets into headers naively and turn the form into a spam engine.

There are exactly 10000 lines if Python's email module in version 3.9.

In addition, there is a concrete risk creating some way of locally executing code, as the template being filled out is then used as an input for a local program – in this case, whatever you use as sendmail. In theory, I'm pretty sure this is not a problem here, as no user-controlled input goes into the headers. If you change this, either sanitise the input, probably by clamping everything down to printable ASCII and normalising whitespace, or by parsing them yourself. The message content, on the other hand, gets properly MIME-encapsulated. In practice, I can't say I trust Python's email package too much, as by Python stdlib standards, it feels not terribly mature and is probably less widely used than one may think.

But that's a risk I'm willing to take; even if someone spots a problem in the email module, shodan or a similar service still has no way to automatically figure out that it is in use in this form, and my page's insignificance makes it extremely unlikely that someone will do a targeted attack on day 0. Or even day 10.

But then, perhaps this is a good occasion to read through email's source code? Fun fact: in python 3.9, a find . -name "*.py" | xargs wc -l gives exactly 10000 lines. And my instinct that headers are the trickiest part is probably right, too: 3003 of those are in _header_value_parser.py.
⎆

Jost Bürgi, der Sinus und Umberto Eco

2022-03-21

Kryptische, leicht mathematisch aussehende Zeicnung

„Reimers' Diagramm“: Für 400 Jahre der einzige Hinweis darauf, wie Jost Bürgi wohl seine Sinustabelle berechnet hat. Nicht mal Kepler hat das Rätsel lösen können.

Bürgis geniale Methode
Ein numerischer Traum

Ein Geheimnis, das im antiken Griechenland ein wenig angekratzt wurde, über das dann Gelehrte in Indien und Arabien nachgedacht haben, für das in der zweiten Hälfte des 16. Jahrhunderts ein wandernder Schweizer Uhrmacher eine genial einfache Lösung fand, von der aber die Nachwelt nur ein paar kryptische Referenzen hat, bis ein unerschrockener Wissenschaftler in den Tiefen längst vergessener Bibliotheken ein Manuskript entdeckt, das des Rätsels Lösung enthält: Es gibt Geschichten, die klingen nach einem Roman von Umberto Eco (oder, je nach Temperament und/oder Geschmack, Dan Brown) und sind doch wahr.

Auf die Geschichte von Jost Bürgis Sinusberechnung bin ich über die DLF-Sternzeit vom 27.2.2022 gekommen, und sie geht ungefähr so: Nachdem Hipparchos von Nicäa so um die 150 vdcE[1] nicht nur den ersten brauchbaren Sternkatalog vorgelegt, sondern auch die ersten ordentlichen Rezepte angegeben hatte, wie mensch für jede Menge Winkel den zugehörigen Sinus [2] ausrechnet, gab es zur Berechnung trigonometrischer Funktionen sehr lange nicht viel Neues.

Klar, der große Ptolomaios, genau, der mit dem Weltbild, hat Hipparchos' Methode zur Berechnung des Sinus über regelmäßige Vielecke kanonisiert. In Indien gab es einige Fortschritte – etwa den Übergang von der Sehne zum Sinus –, in Arabien wurden verschiedene Ergebnisse zusammengetragen und systematisiert, aber immer war es eine mühsame, geometrisch insprierte, endlose Rechnerei.

Und dann kommen wir in die frühe Neuzeit in Europa, genauer die zweite Hälfte des 16. Jahrhunderts. Kopernikus hat noch einmal ganz klassisch den Sinus mit Vielecken berechnet, während er die Konflikte zwischen Ptolomaios und der Realität untersuchte. In Italien macht sich allmählich Galileo bereit, die Physik als experimentelle Naturwissenschaft zu etablieren. Und in Kassel, beim wissenschaftsbegeisterten hessischen Landgraf Wilhelm IV, sammeln sich ein paar Mathe- und Astro-Nerds, die beim ebenso berühmten wie fiesen Tycho gelernt haben, unter ihnen Nicolaus Reimers, der das kryptische Bild über diesem Post veröffentlicht hat, vermutlich, weil er versprochen hatte, nicht mehr zu verraten.

Bürgis geniale Methode

Es weist auf ein Verfahren zur Berechnung von Werten der Sinusfunktion hin, das nichts mehr mit den umschriebenen Polygonen des Hipparchos zu tun hat. Sein Erfinder, der Toggenburger Uhrmacher-Astronom-Erfinder Jost Bürgi, hatte zu dieser Zeit ein großes Tabellenwerk vorgelegt, mit dem mensch auch ohne Taschenrechner rausbekommen konnte, wie viel wohl sin(27^∘ 32’ 16’’) sei[3]. Offensichtlich funktionierte die Methode. Doch hat Bürgi – Autodidakt und vielleicht etwas verschroben – die Methode nie richtig veröffentlicht, und so brüteten MathematikerInnen, unter ihnen wie gesagt Johannes Kepler, der immerhin die Sache mit den Ellipsenbahnen im Planetensystem rausbekommen hat, lang über der eigenartigen Grafik. Und kamen und kamen nicht weiter.

Das war der Stand der Dinge bis ungefähr 2014, als der (emeritierte) Münchner Wissenschaftshistoriker Menso Folkerts im Regal IV Qu. 38ª der Universitätsbibliothek in Wrocław auf eine lange übersehene gebundene Handschrift stieß. Ein wenig konnte er ihre Geschichte nachvollziehen: Jost Bürgi persönlich hatte das Werk Kaiser Rudolf II – dem Mäzen von Tycho und Kepler – am 22. Juli 1592 (gregorianisch) in Prag übergeben, was ihm eine Zuwendung von 3000 Talern eingebracht hat. Ich habe leider nicht die Spur eines Gefühls, wie sich der Betrag mit modernen Drittmittelanträgen vergleicht. Die Form des Antrags jedenfalls ist aus heutiger Sicht als unkonventionell einzustufen.

Das Werk fand seinen Weg in das Augustinerkloster im unterschlesischen Sagan (heute Żagań). Wie es dort hinkam, ist nicht überliefert, aber mensch mag durchaus eine Verbindung sehen zu Keplers Aufenthalt in Sagan in den Jahren von 1628 bis 1630. Vielleicht hat er das Buch ja nach seinen Diensten in Prag mitgenommen, auf seinen verschiedenen Wanderungen mitgenommen und schließlich selbst im Kloster gelassen? Aber warum hätte er dann über Bürgis Methode gerätselt?

Wie auch immer: Im Gefolge des Reichsdeputationshauptschlusses wurde das Kloster 1810 aufgelöst – ich stelle mir das ein wenig vor wie in Poes „Die Grube und das Pendel“ –, und der Bestand der Bibliothek fiel an die Universität Breslau, die wiederum nach dem zweiten Weltkrieg zur polnischen Uni Wrocław wurde.

In diesem geschichtsträchtigen Manuskript erklärt Bürgi seinen Algorithmus. Dargestellt ist das in der Abhandlung von Folkerts et al (arXiv:1510.03180), in der sich auf den Seiten 11 und 12 auch die Auflösung für Reimers' Rätsel findet. Etwas schöner beschreibt das Verfahren Denis Roegel in seinem Aufsatz Jost Bürgi's skillful computation of sines. Dort wird auch Bürgis mutmaßliche Grundeinsicht besser erläutert, nach der der Sinus einfach das Ding ist, das, modern gesprochen, nach zweifacher Ableitung sich selbst (mal minus eins) ergibt. Das ist der mathematische Hintergrund dafür, dass folgendes Stück Python tatsächlich relativ schnell zu einer Tabelle der Sinuswerte von n im ersten Quadranten gleichverteilten Winkeln konvergiert:

tot_sines = list(range(n+1))
for iter_index in range(n_iter):

  intermediates = [tot_sines[-1]//2]
  for tot in reversed(tot_sines[1:-1]):
    intermediates.append(intermediates[-1]+tot)

  tot_sines = [0]
  for diff in reversed(intermediates):
    tot_sines.append(tot_sines[-1]+diff)

return dict((k*math.pi/2/n,  v/tot_sines[-1])
  for k, v in enumerate(tot_sines))

– es ist, glaube ich, unmöglich, zu verstehen, was hier passiert (und warum), ohne den Roegel oder zumindest den Folkerts gelesen zu haben. Aber ich könnte andächtig werden angesichts so simpler Manipulationen, die so schnell zu richtig vielen Stellen einer transzendenten Funktion wie des Sinus führen.

Ein numerischer Traum

Wie schnell das mit den vielen Stellen bei Bürgis Algorithmus geht, zeigt folgende Grafik:

Hier läuft horizontal der Winkel – und der Algorithmus funktioniert wirklich nur, wenn das einen rechten Winkel einschließt –, vertikal die Iterationen von Bürgis Algorithmus. Die Farben entsprechen dem dekadischen Logarithmus der Abweichung der Bürgi-Ergebnisse von dem, was die Python-Standardbibliothek gibt, im Groben also die Zahl der Stellen, die der Algorithmus richtig hinbekommt. Mehr als 18 geht da schon mal nicht, weil die Zahlen von Python in 64-bittigen IEEE-Fließkommazahlen („double precision“) kommen, und mehr als 18 Dezimalstellen sind da nicht drin (in der Grafik steckt die Zusatzannahme, dass wir von Zahlen in der Größenordnung von eins sprechen).

Mithin gewinnt der Algorithmus pro Iteration ungefähr eine Dezimalstelle, und das gleichmäßig über den ganzen Quadranten. DemoprogrammiererInnen: billiger kommt ihr, glaube ich, nicht an eine beliebig präzise Sinustabelle ran.

Spannend fand ich noch die kleinen dunkelblauen Klötzchen ganz unten in der Grafik: Hier werden sich Bürgi und Python offenbar auf Dauer nicht einig. So, wie ich das geschrieben habe, würde ich fast eher Bürgi vertrauen, denn bei den Ganzzahlen, die da verwendet werden, kann eigentlich nichts schief gehen. Andererseits sind Fließkommazahlen eine heikle Angelegenheit, insbesondere, wenn es ums letzte Bit geht. Um mich zu überzeugen, dass es nur um genau jenes unheimliche letzte Bit geht, habe ich mir geschwind eine Funktion geschrieben, die die Fließkommazahlen vinär ausgibt, und der Code gefällt mir so gut, dass ich sie hier kurz zeigen will:

import struct

_BYTE_LUT = dict((v, "{:08b}".format(v)) for v in range(256))
def float_to_bits(val):
  return "".join(_BYTE_LUT[v] for v in struct.pack(">d", val))

Mit anderen Worten lasse ich mir geschwind ausrechnen, wie jedes Byte in binär aussehen soll (_BYTE_LUT), wobei die Python-Bibliothek mit dem 08b-Format die eigentliche Arbeit macht, und dann lasse ich mir die Bytes der Fließkommazahl vom struct-Modul ausrechnen. Der einzige Trick ist, dass ich das Big-end-first bestellen muss, also mit dem signfikantesten Byte am „linken“ Ende. Tue ich das nicht, ist z.B. auf Intel-Geräten alles durcheinander, weil die Bits in der konventionellen Schreibweise daherkommen, die Bytes aber (wie bei Intel üblich) umgedreht, was ein furchtbares Durcheinander gibt.

Jedenfalls: Tatsächlich unterscheiden sich die Werte schon nach 20 Iterationen nur noch im letzten bit, was für 45 Grad alias π/4 z.B. so aussieht:

45
  0011111111100110101000001001111001100110011111110011101111001101
  0011111111100110101000001001111001100110011111110011101111001100

Ich lege mich jetzt mal nicht fest, was das „bessere“ Ergebnis ist; ich hatte kurz überlegt, ob ich z.B. mit gmpy2 einfach noch ein paar Stellen mehr ausrechnen sollte und sehen, welches Ergebnis näher dran ist, aber dann hat mich meine abergläubische Scheu vor dem letzten Bit von Fließkommazahlen davon abgehalten.

Wer selbst spielen will: Meine Implementation des Bürgi-Algorithmus, der Code zur Erzeugung der Grafik und die Bitvergleicherei sind alle enthalten in buergi.py.

[1]

Das vdcE bringe ich hiermit als Übertragung von BCE, before the Christian era, in Gebrauch. Und zwar, weil v.Chr völlig albern ist (es ist ja nicht mal klar, ob es irgendeine konkrete Figur „Christus“ eigentlich gab; wenn, ist sie jedenfalls ganz sicher nicht zur aktuellen Epoche – also dem 1. Januar 1 – geboren) und „vor unserer Zeitrechnung“ ist auch Quatsch, denn natürlich geht Zeitrechnung (wenn auch mangels Jahr 0 etwas mühsam) auch vor der Epoche. „Vor der christlichen Epoche“ hingegen bringt ganz schön auf den Punkt, was das ist, denn die „Epoche“ ist in der Zeitrechnung einfach deren Nullpunkt (oder halt, vergurkt wie das alles ist, hier der Einspunkt).

[2]	Na ja, in Wirklichkeit hat er mit der Länge der Sehne gerechnet, die ein Winkel in einem Kreis aufspannt, aber das ist im Wesentlichen das Gleiche wie der Sinus, der ja gerade der Hälfte dieser Sehne entspricht.

[3]	Ich bleibe natürlich bei meiner Verurteilung …

Werkstattbericht: Kohlendioxid auf dem Balkon

2022-02-06

Im November hatte ich mich gefragt, was wohl die recht deutlichen Spitzen der CO₂-Konzentration auf meinem Balkon verursachen mag, die sich da immer mal wieder zeigen. Um Antworten zu finden, habe ich seit Ende Dezember eine längere Messreihe laufen lassen und derweil vierstündlich Windrichtungen von der Open Weathermap aufgenommen. Das, so hoffte ich, sollte zeigen, woher der Wind weht, wenn die Konzentration auffällige Spitzen hat.

Leider gibt ein schlichter optischer Vergleich von Konzentration (oben) und Windrichtung (unten; hier als Cosinus, damit das Umschlagen von 0 auf 360 Grad nicht so hässlich aussieht) nicht viel her:

Zwei unterbrochene Kurven, die jeweils recht munter vor sich hinwackeln

CO₂-Konzentration auf meinem Balkon und Windrichtung für Heidelberg aus der Open Weathermap zwischen Ende Dezember 2021 und Anfang Februar 2022. Die Lücken ergeben sich aus fehlenden Daten zur Windrichtung.

Tatsächlich hilft es ein wenig, wenn mensch das anders plottet. Unten bespreche ich kurz das Programm, das Wind- und CO₂-Daten zusammenbringt. Dieses Programm produziert auch folgenden Plot in Polarkoordinaten:

Scatterplot in Polarkoordinaten: Im Wesentlichen ein oranger Ring

CO₂-Konzentration auf meinem Balkon gegen meteorologische Windrichtung (also: Herkunft des Windes, hier gezählt ab Nord über Ost, so dass das orientiert ist wie eine Landkarte) und farbkodierte Windgeschwindigkeit (in Meter pro Sekunde). Das ist ein PNG und kein SVG, weil da doch viele Punkte drauf sind und Browser mit so großen SVGs immer noch ins Schlingern kommen.

Ich hatte mich seit einem Monat auf diesen Plot gefreut, weil ich erwartet habe, darin eine ordentliche „Beule“ zu sehen dort, wo die CO₂-Emission herkommt. Gemessen daran ist wirkliche Ergebnis eher ernüchternd. Dort, wo ich die Abgasfahne des Großkraftwerk Mannheim sehen würde, etwas unterhalb der 270°-Linie, ist allenfalls ein kleines Signälchen und jedenfalls nichts, was ich wirklich ernst nehmen würde.

Etwas deutlicher zeichnet sich etwas zwischen 280 und 305 Grad ab, also Westnordwest. Das könnte die Ladenburger Chemieindustrie oder die BASF in Ludwigshafen sein; zu letzterer haben die kritischen Aktionäre im letzten Jahr angesagt, sie emittiere als Konzern 20 Megatonnen Kohlendioxid im Jahr. Wenn, was nicht unplausibel ist, die Hälfte davon am Standort Ludwigshafen anfällt, würden sich diese 10 Mt ganz gut vergleichen mit den 8 Mt, die ich neulich fürs Großkraftwerk gefunden hatte – die Abschätzung von dort, so eine Abgasfahne könne durchaus die Konzentrationsspitzen erklären, kommt also auch für die BASF hin. Allerdings wird deren Emission angesichts des riesigen Werksgeländes natürlich auch verteilter sein…

Also: Überzeugend ist das alles nicht. Ein anderes Feature ist jedoch schlagend, wegen weniger Übermalung – die bei beiden Plots ein echtes Problem ist; nächstes Mal muss ich mit halbtransparenten Punkten arbeiten – noch mehr, wenn ich den Polarplot „ausrolle“, also den Winkel von unten nach oben laufen lasse:

Scatterplot kartesisch: ein starker dunkler Klops bei 230 Grad

In dieser Darstellung fällt ins Auge, dass die CO₂-Konzentration bei starken (dunkle Töne) Südwest- (um die 225°) -strömungen recht drastisch fällt. Das passt sehr gut zu meinen Erwartungen: Südwestwind schafft hier in der Rheinebene Luft durch die Burgundische Pforte, hinter der im Mittelmeerraum auch jetzt im Winter eifrig Photosynthese stattfindet. Wer drauf aufpasst, sieht die Entsprechungen auch im Polarplot von oben, in dem dann sogar auffällt, dass reiner Südwind gelegentlich noch besser photosynthetisierte Luft heranführt, auch wenn der Wind nicht ganz so stark bläst.

Demgegenüber ist mir eigentlich alles, was sich im nordöstlichen Quadranten des Polarplots (und hier zwischen 0 und 90 Grad) abspielt, eher rätselhaft. Der doppelseitige Sporn bei genau 90 Grad ist vermutlich auf Datenmüll der Wetterstation zurückzuführen: Wahrscheinlich hat die einen Bias, der bei wenig Wind diese 90 Grad ausspuckt. Selbst nach meiner Interpolation (vgl. unten) ist das noch zu ahnen, wenn mensch die Verteilung der Geschwindigkeiten insgesamt (in rot) und die der Geschwindigkeiten rund um einen auffälligen Hügel rund um 90° Windrichtung herum (in blau) ansieht:

Zwei Histogramme über Geschwindigkeiten, bei dem das blaue nur im linken Bereich ist

Die elegante Schleife, die von (0, 500) über (70, 540) nach (90, 510) führt und die im Polarplot ganz alleine außen vor sich hinläuft, dürfte ziemlich sicher teils physikalisch sein. Dass das da so einen Ring macht, dürfte zwar ein Artefakt meiner gewagten Interpolation sein (vgl. Technics). Der Anstieg als solcher und wohl auch die grobe Verortung dürften aber ganz gut hinkommen. Sieht mensch sich das im zeitlichen Verlauf an, entspricht die Schleife der höchsten Spitze in der ganzen Zeitreihe.

Nur leider ist im Nordosten von meinem Balkon nicht mehr viel: Ein paar Dutzend Häuser und dann der Odenwald, also für fast 10 km nur Bäume. Na gut, und ein Ausflugsrestaurant.

Die aus meiner Sicht plausibelste Interpretation für diese Stelle basiert auf der Beobachtung, dass in der fraglichen Zeit (am 10.1.) wenig Wind wehte, die Temperaturen aber ziemlich niedrig lagen. Vielleicht schauen wir hier wirklich auf die Heizungen der Umgebung? Der Schlot unserer lokalen Gemeinschafts-Gasheizung ist in der Tat so in etwa im Nordosten des Balkons – und vielleicht wurde ja sonst nicht so viel geheizt?

Technics

Die wesentliche Schwierigkeit in diesem Fall war, dass ich viel engmaschiger CO₂-Konzentrationen (alle paar Minuten) habe als Windrichtungen (bestenfalls alle vier Stunden), und zudem viele Windrichtungen aus welchen Gründen auch immer (offensichtlich wäre etwa: zu wenig Wind) fehlen. Auf der positiven Seite erzeugt mein Open Weathermap-Harvester weathercheck.py eine SQLite-Datenbank, so dass ich, wenn es nicht furchtbar schnell gehen muss, recht bequem interessante Anfragen laufen lassen kann.

Mein Grundgedanke war, die beiden einem CO₂-Wert nächsten Wind-Werte zu bekommen und dann linear zu interpolieren[1]. Das ist schon deshalb attraktiv, weil die Zeit (als Sekunden seit 1.1.1970) als Primärschlüssel der Tablle deklariert ist und deshalb ohnehin ein Index darauf liegt.

Dabei sind aber je nach Datenverfügbarkeit ein Haufen Fälle zu unterscheiden, was zu hässlichen if-else-Ketten führt:

def get_for_time(self, time, col_name, default=None):
  res = list(self.conn.execute(f"SELECT timestamp, {col_name} FROM climate"
    " WHERE TIMESTAMP BETWEEN ? AND ?"
    " ORDER BY ABS(timestamp-?) LIMIT 2",
    (time-40000, time+40000, time)))

  if len(res)!=2:
    if default is not None:
      return default
    raise ValueError(f"No data points close to {time}")

  elif abs(res[0][0]-time)<200 and res[0][1] is not None:
    return res[0][1]

  elif res[0][1] is None or res[1][1] is None:
    if default is not None:
      return default
    raise ValueError("One or more limits missing.  Cannot interpolate.")

  else:
    t1, v1 = res[0]
    t2, v2 = res[1]
    return (v1*(t2-time)+v2*(time-t1))/(t2-t1)

Die Fallunterscheidung ist:

Es gibt überhaupt keine Daten innerhalb von einem halben Tag. Dann kann ich nur einen Fehler werfen; zumindest in unseren Breiten sind Windrichtungen eigentlich schon über kürzere Zeiträume hinweg nur lose korreliert.
Innerhalb von 200 Sekunden der gesuchten Zeit gibt es einen tatsächlichen Messwert, und dieser ist nicht NULL. Dann gebe ich den direkt zurück.
Einer der beiden Werte, die um die gesuchte Zeit herum liegen, fehlt (also ist NULL). Dann kann ich nicht interpolieren und muss wieder einen Fehler werfen. Hier wäre es nicht viel unplausibler als die Interpolation, wenn ich einfach einen nicht-NULL-Wert nehmen würde; aber es wäre doch nochmal ein Stückchen spekulativer.
Ansonsten mache ich einfach eine lineare Interpolation.

NULL-Werte machen die Dinge immer komplex. Aber wenn ihr euch überlegt, wie viel Stress sowas ohne SQL wäre, ist das, finde ich, immer noch ganz elegant. Im echten Code kommt noch etwas Zusatzkomplexität dazu, weil ich Winkel interpolieren will und dabei immer die Frage ist, wie mensch die Identität von 360 und 0 Grad einrührt.

Eine vorsorgliche Warnung: aus der Art, wie ich den Spaltennamen hier reinfummele, folgt, dass, wer den Parameter kontrolliert, beliebiges SQL ausführen kann. Sprich: wer diesen Code irgendwie Web-zugänglich macht, darf keine unvalidierte Eingabe in col_name reinlassen.

Eingestandenermaßen ist diese Sorte von datenbankbasierter Interpolation nicht furchtbar effizient, aber für die 100000 Punkte, die ich im Augenblick plotten will, reicht es. Siehe: Den Code.

[1]	Klar: Windrichtungen über Stunden linear zu interpolieren ist in den meisten Wetterlagen eher zweifelhaft. So, wie ich meine Plots mache, ist es aber nicht wesentlich verschieden davon, die Punkte über den Bereich zu verschmieren. Das wiederum wäre konzeptionell gar nicht so arg falsch.

Inlining xs:include in XML Schema

2022-01-31

Please don't do it like this: for users of a schema, having to pull it in a dozen fragments is just pain and no gain. See below for a program that lets you heal this particular disease.

While I'm a big fan of XML – which is governed by a very well-written standard and is (DTDs aside) about as easy to process as something context-free can be –, I have always been a lot more skeptical about XML Schema, which is horrendously complex, has a few nasty misfeatures[1] and generally has had a major role in giving XML a bad name.

But well, it's there, and it won't go away. That ought to be reason enough to not encumber it with further and totally avoidable pain. As, for instance, splitting up a single schema into fragments of a couple of lines and then using xs:include liberally to re-assemble the fragments at the client side. Datacite, I'm looking at you. Regrettably, they're not the only ones doing that. And as opposed to splitting a domain mapping into different schemas – which might improve re-usability, this lexical splitting really helps nobody except perhaps the authors.

The use of xs:include is a pain in particular when one tries to implement redistributable validators, as these then need to keep a lot of files in a defined hierarchy. Just pointing to the vendor's site is not an option, because the software would hit that every time it validates something, which is, if nothing else, a privacy and stability problem.

Well, today I had another case of XSD splititis, and this time it was bad enough that I decided to merge the fragments. I had expected people had written “inliners” expanding xs:include in XSDs into standalone XSDs. After all, it's basically a lexical thing (well, excepting namespace mappings and perhaps re-indentation). Five minutes of operating a search engine didn't bring up anything, though, and so I wrote a quick cure: expand-xsd-include.py.

I'll be the first to admit that it is a hack at this point, mainly because I blindly discard the root tags for the included documents. That's wrong because these might add or, worse, change the mapping from prefixes to XML namespaces. In its current state, this will fail badly if included documents use different or extra mappings and declare them in the root element; declarations further down are ok.

Another problem resulting from keeping namespace processing off on the parser is that I hardcode the prefix for the XSD schema to xs. If your schema uses something else, change XSD_PREFIX in the script.

Mending these deficiencies wouldn't be an undue effort, and if you have XSDs that need it, let me know and I'll do proper namespace processing. Or perhaps, in addition, teach the thing to pull the input files via http. Meanwhile, I suspect that the large majority of atomised XSDs can be merged with this code, and so I thought I might as well put it online in its slightly embarrassing shape.

Let me know if you use it. And if you distribute fragmented XSDs: Why not use the script to assemble your XSD before publishing it?

[1]

The worst XSD misfeature IMHO are the namespaced attribute values; where XML has been designed to be parsable without external DTDs (ok, not generally, but under well-defined conditions, and it's been a long time since I saw a document that didn't meet those), parsing results with namespaced attribute values depend on whether or not the parser knows the XSD. And that would even be bad without the ugly schemaLocation hacks in both schema and schema instance.

Explaining Tags in Pelican

2022-01-19
Right after I had celebrated the first anniversary of this blog with the post on my Pelican setup, I decided to write another plugin I've been planning to write for a while: taginfo.py.

Nachtrag (2022-10-07)

Don't take it from here; rather, see https://codeberg.org/AnselmF/pelican-ext

This is for:

that is, including explanations in on pages for tags, telling people what the tag is supposed to mean.

To use taginfo, put the file into your plugins folder, add taginfo to the PLUGINS list in your pelicanconf.py, and then create a folder taginfo next to your content folder. In there, for each tag you want to comment, create a file <tagname>.rstx (or just rst). Such a file has to contain reStructuredText, where pelican's extensions (e.g., {filename} links) do not work (yet). I suppose it wouldn't be hard to support them; if you're interested in this plugin, feel free to poke me in case you'd like to see the extra pelican markup.

To make the descriptions visible, you need to change your tag.html template (typically in theme/templates/tag.html) in order to arrange for tag.make_description() to be callsed when rendering the document. Me, I'm doing it like this:
```
{% block content_title %}
<h1>Tag <em>{{ tag }}</em></h1>
<div id="taginfo">
        {{ tag.make_description() }}
</div>
{% endblock %}
```
(And I still find jinja templates exceptionally ugly).
⎆
Fernseh-Livestreams mit Python und mpv

2021-12-05
- Technics
- Sicherheitsfragen
Nachtrag (2022-10-11)

Ich habe das Programm jetzt auf codeberg untergebracht

Die Unsitte, alles in den Browser zu verlegen, ist ja schon aus einer Freiheitsperspektive zu verurteilen – „die Plattform“ gibt die Benutzerschnittstelle vor, kann die Software jederzeit abschalten und sieht (potenziell) noch den kleinsten Klick der NutzerIn (vgl. WWWorst App Store). Bei Mediatheken und Livestreams kommt noch dazu, dass videoabspielende Browser jedenfalls in der Vergangenheit gerne mal einen Faktor zwei oder drei mehr Strom verbraucht haben als ordentliche Videosoftware, von vermeidbaren Hakeleien aufgrund von schlechter Hardwarenutzung und daraus resultierendem Elektroschrott ganz zu schweigen.

Es gibt also viele Gründe, speziell im Videobereich den Web-Gefängnissen entkommen zu wollen. Für übliche Videoplattformen gibt es dafür Meisterwerke der EntwicklerInnengeduld wie youtube-dl oder streamlink.

Für die Live-Ströme der öffentlich-rechtlichen Anstalten hingegen habe ich zumindest nichts Paketiertes gefunden. Vor vielen Jahren hatte ich von einem Freund ein paar screenscrapende Zeilen Python erledigt dazu. In diesen Zeiten müssen allerdings Webseiten alle paar Monate komplett umgeschrieben („jquery ist doch total alt, angular.js hat auch schon bessere Tage gesehen“) und regelauncht werden, und das Skript ging mit einem Relaunch ca. 2015 kaputt. Als ich am Freitag die Tagesschau ansehen wollte, ohne DVB-Hardware zu haben, habe ich mich deshalb nach einer Neufassung des Skripts umgesehen.

Das Ergebnis war eine uralte Seite mit mpv-Kommandozeilen und ein Verweis auf ein von den MediathekView-Leuten gepflegtes Verzeichnis von Live-Strömen. Da stehen zwar oben alt aussehende Timestamps drin, das Log aber zeigt, dass der Kram durchaus gepflegt wird.

Aus letzterem habe ich livetv.py (ja, das ist ein Download-Link) gestrickt, ein weiteres meiner Ein-Datei-Programme. Installiert mpv (und Python, klar), macht chmod +x livetv.py und sagt dann ./livetv.py ARD Livestream – fertig. Bequemer ist es natürlich, das Skript einfach irgendwo in den Pfad zu legen.

Das Argument, das das Programm haben will, kann irgendeine Zeichenfolge sein, die eindeutig einen Sender identifiziert, also nur in einem Sendernamen vorkommt. Welche Sender es gibt, gibt das Programm aus, wenn es ohne Argumente aufgerufen wird:
```
$ livetv.py
3Sat Livestream
...
PHOENIX Livestream
```
Mit der aktuellen Liste könnt ihr z.B. livetv.py Hamburg sagen, weil „Hamburg“ (auch nach Normalisierung auf Kleinbuchstaben) nur in einer Stationsbezeichnung vorkommt, während „SWR“ auf eine Rückfrage führt:
```
$ livetv.py SWR
SWR BW Livestream? SWR RP Livestream?
```
„SWR BW“ (mit oder ohne Quotes auf der Kommandozeile) ist dann eindeutig, woraufhin livetv an den mpv übergibt.

Ich gehe davon aus, dass die Anstalten die URLs ihrer Streams auch weiterhin munter verändern werden. Deshalb kann sich das Programm neue URLs von den MediathekView-Leuten holen, und zwar durch den Aufruf:
```
$ livety.py update
```
Das schreibt, wenn alles gut geht, die Programmdatei neu und funktioniert mithin nur, wenn ihr Schreibrechte auf das Verzeichnis habt, in dem livetv.py liegt – was besser nicht der Fall sein sollte, wenn ihr es z.B. nach /usr/local/bin geschoben habt.

A propos Sicherheitsüberlegungen: Der update-Teil vertraut gegenwärtig ein wenig den MediathekView-Repo – ich entschärfe zwar die offensichtlichsten Probleme, die durch Kopieren heruntergeladenen Materials in ausführbaren Code entstehen, aber ich verspreche nicht, raffinierteren Angriffen zu widerstehen. Abgesehen vom update-Teil halte ich das Programm für sicherheits-unkritisch. Es redet selbst auch nicht mit dem Netz, sondern überlässt das dem mpv.

Livetv.py sagt per Voreinstellung dem mpv, es solle einen „vernünftigen“ Stream aussuchen, was sich im Augenblick zu „2 Mbit/s oder weniger“ übersetzt. Wer eine andere Auffassung von „vernünftig“ hat, kann die --max-bitrate-Option verwenden, die einfach an mpvs --hls-bitrate weitergereicht wird. Damit könnt ihr
```
$ livetv.py --max-bitrate min arte.de
```
für etwas sagen, das für die Sender, die ich geprüft habe, auch auf sehr alten Geräten noch geht,
```
$ livetv.py --max-bitrate max arte.fr
```
für HD-Wahnsinn oder
```
$ livetv.py --max-bitrate 4000000 dw live
```
für einen Stream, der nicht mehr als 4 MB/s verbraucht.
Technics

Die größte Fummelei war, die Kanalliste geparst zu bekommen, denn aus Gründen, für die meine Fantasie nicht ausreicht (MediathekView-Leute: Ich wäre echt neugierig, warum ihr das so gemacht habt), kommen die Sender in einem JSON-Objekt (statt einer Liste), und jeder Sender hat den gleichen Schlüssel:
```
"X" : [ "3Sat", "Livestream", ...
"X" : [ "ARD", "Livestream", ...
```
– ein einfaches json.loads liefert also ein Dictionary, in dem nur ein Kanal enthalten ist.

Auch wenn ich sowas noch nie gesehen habe, ist es offenbar nicht ganz unüblich, denn der json-Parser aus der Python-Standardbibliothek ist darauf vorbereitet. Wer ein JSONDecoder-Objekt konstruiert, kann in object_pairs_hook eine Funktion übergeben, die entscheiden kann, was mit solchen mehrfach besetzen Schlüsseln pasieren soll. Sie bekommt vom Parser eine Sequenz von Schlüssel-Wert-Paaren übergeben.

Für meine spezielle Anwendung will ich lediglich ein Mapping von Stationstiteln (in Element 3 der Kanaldefinition) zu Stream-URLs (in Element 8) rausziehen und den Rest der Information wegwerfen. Deshalb reicht mir Code wie dieser:
```
def load_stations():
  channels = {}
  def collect(args):
    for name, val in args:
      if name=="X":
        channels[val[2]] = val[8]

  dec = json.JSONDecoder(object_pairs_hook=collect)
  dec.decode(LIST_CACHE)

  return channels
```
– das channels-Dictionary, das collect nach und nach füllt, ist wegen Pythons Scoping-Regeln das, das load_stations definiert. Die collect-Funktion ist also eine Closure, eine Funktion, die Teile ihres Definitionsumfelds einpackt und mitnehmt. So etwas macht das Leben von AutorInnen von Code sehr oft leichter – aber vielleicht nicht das Leben der späteren LeserInnen. Dass die collect-Funktion als ein Seiteneffekt von dec.decode(...) aufgerufen wird und dadurch channels gefüllt wird, braucht jedenfalls erstmal etwas Überlegung.

Der andere interessante Aspekt am Code ist, dass ich die Liste der Live-Streams nicht separat irgendwo ablegen wollte. Das Ganze soll ja ein Ein-Datei-Programm sein, das einfach und ohne Installation überall läuft, wo es Python und mpv gibt. Ein Blick ins Commit-Log der Kanalliste verrät, dass sich diese allein im letzten Jahr über ein dutzend Mal geändert hat (herzlichen Dank an dieser Stelle an die Maintainer!). Es braucht also eine Möglichkeit, sie aktuell zu halten, wenn ich die Liste nicht bei jedem Aufruf erneut aus dem Netz holen will. Das aber will ich auf keinen Fall, weniger, um github zu schonen, mehr, weil sonst github sehen kann, wer so alles wann livetv.py verwendet.

Ich könnte die Liste beim ersten Programmstart holen und irgendwo im Home (oder gar unter /var/tmp) speichern. Aber dann setzt zumindest der erste Aufruf einen Datenpunkt bei github, und zwar für neue NutzerInnen eher überraschend. Das kann ich verhindern, wenn ich die Liste einfach im Programm selbst speichere, also selbstverändernden Code schreibe.

Das ist in interpretierten Sprachen eigentlich nicht schwierig, da bei ihnen Quellcode und ausgeführtes Programm identisch sind. Zu den großartigen Ideen in Unix gehört weiter, dass (das Äquivalent von) sys.argv[0] den Pfad zur gerade ausgeführten Datei enthält. Und so dachte ich mir, ich ziehe mir einfach den eigenen Programmcode und ersetze die Zuweisung des LIST_CACHE (das json-Literal von github) per Holzhammer, also regulärem Ausdruck. In Code:
```
self_path = sys.argv[0]
with open(self_path, "rb") as f:
  src = f.read()

src = re.sub(b'(?s)LIST_CACHE = """.*?"""',
  b'LIST_CACHE = """%s"""'%(in_bytes.replace(b'"', b'\\"')),
  src)

with open(self_path, "wb") as f:
  f.write(src)
```
Dass das Schreiben ein eigener, fast atomarer Schritt ist, ist Vorsicht: Wenn beim Ersetzen etwas schief geht und das Programm eine Exception wirft, ist das open(... "wb") noch nicht gelaufen. Es leert ja die Programmdatei, und solange es das nicht getan hat, hat mensch eine zweite Chance. Ähnlich übrigens meine Überlegung, das alles in Binärstrings zu bearbeiten: beim Enkodieren kann es immer mal Probleme geben, die am Ende zu teilgeschriebenen Dateien führen können. Vermeide ich Umkodierungen, kann zumindest die Sorte von Fehler nicht auftreten.

Wie dem auch sei: Dieser Code funktioniert nicht. Und zwar in recht typischer Weise innerhalb der Familie von Quines und anderen Selbstanwendungsproblemen: das re.sub erwischt auch seine beiden ersten Argumente, denn beide passen auf das Muster LIST_CACHE = """.*?""". Deshalb würden von livetv.py update auch diese beiden durch das json-Literal mit den Senderdefinitionen ersetzt. Das so geänderte Programm hat zwei Syntaxfehler, weil das json natürlich nicht in die String-Literale passt, und selbst wenn es das täte, gingen keine weiteren Updates mehr, da die Such- und Ersatzpatterns wegersetzt wären.

Eine Lösung in diesem Fall ist geradezu billig: in Python kann mensch ein Leerzeichen auch als '\x20' schreiben (das ASCII-Zeichen Nummer 0x20 oder 32), und schon matcht der reguläre Ausdruck nicht mehr sich selbst:
```
re.sub(b'(?s)LIST_CACHE\x20= """.*?"""',
  b'LIST_CACHE\x20= """%s"""'...
```
Sicherheitsfragen

Ein Programm, das Daten aus dem Netz in sich selbst einbaut, muss eigentlich eine Ecke vorsichtiger vorgehen als dieses hier. Stellt euch vor, irgendwer bekommt etwas wie:
```
{ "Filmliste": [....,
  "X": ["...
  "Igore": ['"""; os.system("rm -r ~"); """']
}
```
in das MediathekView-Repo committet; das würde für die MediathekView immer noch prima funktionieren, das Objekt mit dem Schlüssel Ignore würde fast sicher tatsächlich einfach ignoriert.

Wer dann allerdings livetv.py update laufen lässt, bekommt den ganzen Kram in Python-Quelltext gepackt, und der Inhalt des Ignore-Schlüssels wird vom Python-Parser gelesen. Der sieht, wie der lange String mit den drei Anführungszeichen geschlossen wird. Danach kommt eine normale Python-Anweisung. Die hier das Home-Verzeichnis der NutzerIn löscht. Python wird die treu ausführen. Bumm.

So funktioniert das in Wirklichkeit zum Glück nicht, denn ich escape im realen Code Anführungszeichen (das .replace(b'"', b'\\"')). Damit …
⎆
Kohlendioxid auf dem Balkon

2021-11-20
CO₂-Konzentrationen auf meinem Straßenbalkon, zusammen mit Windgeschwindigkeiten und Temperaturen. Das ist ein SVG, es lohnt sich also durchaus, in einem separaten Browserfenster in den Plot zu zoomen.
Ich habe neulich eine längere Zeitreihe mit CO₂-Konzentrationen auf meinem „vorderen” Balkon genommen. Zur Einordnung: Das Messgerät steht so etwa 10 Meter über und 15 Meter neben einer halbwegs viel befahrenen Straße. Ob das wohl etwas mit den wilden Schwankungen zu tun hat, die in der Kurve oben vor allem um den 9.11. herum zu sehen sind? Muss ich meine Einschätzung von neulich, einzelne Autos seien selbst im mittleren Nahbereich im CO₂ kaum nachzuweisen (nun: an der frischen Luft, natürlich), revidieren?

Verheizt jemand 100000 Tonnen Kohlenstoff am Tag?

Wer die Kurven von Windgeschwindigkeit[1] und CO₂-Konzentration vergleicht, könnte schon glauben wollen, ohne externe Frischluftzufuhr (also bei niedrigen Windgeschwindigkeiten) gehe das CO₂ lokal merklich nach oben. Wirklich überzeugen kann mich aber keine Korrelation zwischen den verschiedenen geplotteten Größen.

Darum gehe ich die Frage zunächst deduktiv an: woher könnten die enormen Schwankungen der CO₂-Konzentration wohl kommen? Wir reden hier von einer Spanne zwischen 260 ppm und über 400 ppm, wobei es vorkommen kann, dass ich innerhalb von wenigen Stunden 100 ppm mehr CO₂ sehe. Der langfristig ansteigende Trend macht mir übrigens weniger Sorgen: Wenn die Photosyntheserate Richtung Winter dramatisch sinkt, die Emission aber z.B. wegen Heizung eher zunimmt, ist das angesichts der beschränkten globalen Durchmischung der Atmosphäre auf der Erde zu erwarten[2], auch wenn das vielleicht nicht gerade innerhalb von zwei Wochen vonstatten gehen sollte.

Mit den Werkzeugen aus dem Artikel zu meiner Heizleistung von neulich kann mensch abschätzen, was so eine Konzentrationsschwankung in einer lokal gut durchmischten Atmosphäre in, sagen wir, verbranntem Kohlenstoff bedeuten würde.

Dafür muss ich erst überlegen, wie viele CO₂-Teilchen ΔN_CO2, oder der Bequemlichkeit halber eher welche CO₂-Stoffmenge Δn_CO2 = N_CO2 ⁄ A („in mol”) es braucht, um die Konzentration (in ppm, also CO₂-Molekülen pro Million Teilchen insgesamt) innerhalb eines angenommenen Volumens V um das Δc_ppm zu erhöhen, das ich aus dem Plot ablese. Gemäß meinen Rezepten von neulich ist das:

Δn_CO2 = (V)/(V_m)⋅Δc_ppm⋅10⁶,

wobei V_m wieder das Normvolumen ist (22.4 Liter pro mol); das A von oben war die Avogadro-Konstante. Um herauszukriegen, wie viel Kohlenstoff (sagen wir, in Kilogramm) ich verbrennen muss, um diese Änderung quasi durch „frisches“ CO₂ hinzukriegen, muss ich das nur noch mit dem Atomgewicht von Kohlenstoff u_C multiplizieren.

Das Atomgewicht ist, weil Kohlenstoffkerne meist 6 Protonoen und 6 Neutronen enthalten, mit 12 g/mol gut abgeschätzt (ganz genau ist das nicht, vor allem weil in der Atmosphäre auch etwas C-13 und sogar ein wenig C-14 herumschwebt). In dieser Kopfzahl steht das Gramm aus historischen Gründen. Das Mol wurde so definiert, dass die Zahl der Nukleonen im Kern so in etwa das Atomgewicht liefert, als in der Wissenschaft das cgs-System (aus Zentimeter, Gramm und Sekunde) seine große Zeit hatte. Würde mensch das Mol in den heutigen SI-Zeiten (na gut: die meisten AstronomInnen bleiben dem cgs verhaftet und reden zum Beispiel über Energien in erg) definieren, wäre die Avogadro-Konstante um einen Faktor 1000 (nämlich den Faktor zur SI-Einheit Kilogramm) größer.

Wie auch immer: Wenn ich mir mal vorstelle, dass das, was ich da auf meinem Balkon messe, repräsentativ für den Umkreis von 10 km und bis in eine Höhe von 2 km wäre (mensch ahnt schon: Ich eröffne hier eine Reductio ad absurdum), komme ich auf ein Volumen von

V = 2⋅π⋅(10 km)²⋅(2 km) ≈ 1.3⋅10¹² m³

was mit V_m ≈ 0.02 m³ ⁄ mol, einer Änderung von 100 ppm, die mensch als Sprung am 9. und 10.11. sehen kann, sowie der Formel oben auf

Δm_C = u_C⋅(V)/(V_m)⋅Δc_ppm⋅10⁶ ≈ 0.012 kg ⁄ mol⋅(1.3⋅10¹² m³)/(0.02 m³ ⁄ mol)⋅100⋅10^− 6 ≈ 8⋅10⁷ kg

oder achzigtausend Tonnen verbrannten Kohlenstoff führt. Das klingt nach richtig viel und ist es auch. Aber das Volumen, das ich hier betrachte, sind eben auch 1200 Kubikkilometer, und wer sich erinnert, dass ein Kubikmeter eines normalen Gase bei Normalbedingungen um die 1 kg wiegt, kann leicht ausrechnen, dass die Luft in diesem Volumen 1.2⋅10¹² kg (oder 1.2 Milliarden Tonnen – Luft in großen Mengen ist überhaupt nicht leicht) wiegen wird. Dieser ganze Kohlenstoff macht also ungefähr 0.07 Promille (oder 70 Milionstel) der Masse der Atmosphäre aus, was ganz gut mit den 100 ppm in Teilchen zusammengeht, die wir in die ganze Rechnung reingesteckt haben.

Andersrum gerechnet

Tatsächlich kann mensch die Kohlenstoffmasse, die eine Erhöhung der Teilchenkonzentration in einem Gasvolumen bewirkt, auch so herum abschätzen. Der Umrechnungsfaktor von Teilchen- zu Massenkonzentration ist der Faktor zwischen den Dichten von CO₂ und Luft. Das Verhältnis dieser Dichten ist wiederum das der jeweiligen Atommassen, solange jedes Teilchen das gleiche Volumen einnimmt; das schließlich folgt aus der Annahme, dass die Gase ideal sind, was wiederum für unsere Abschätzungen überallhin gut genug ist.

Für CO₂ ist das mit den überwiegend vorkommenden Isotopen von Sauerstoff und Kohlenstoff 16 + 16 + 12 = 44, für Luft, wenn wir nur auf den Stickstoff N₂ schauen, 14 + 14 = 28. Demnach macht 1 ppm in der Teilchenzahl von CO₂ 44 ⁄ 28 ≈ 1.6 ppm in der Masse aus, solange die CO₂-Konzentration so gering ist, dass tatsächlich N₂ die Dichte dominiert.

Andererseits macht Kohlenstoff nur 12 ⁄ 44 ≈ 0.3 an der Masse im CO₂ aus, die Zunahme an Kohlenstoff ist demnach nur ein Drittel von dem gerade berechneten 1.6, also etwas wie 0.5. Folglich werden aus 100 ppm Änderung in der Teilchenzahl etwas wie 100⋅0.5 = 50 ppm Änderung in der Masse; wer das genauer rechnet, bekommt auf diese Weise natürlich das gleiche Resultat wie oben raus.

Wie herum mensch das auch rechnet, es ist klar, dass niemand in der kurzen Zeit so viel Kohlenstoff verbrennt. Ein schneller Reality Check: Meine Kohlendioxid-Kopfzahl war, dass die BRD 2/3 Gigatonnen im Jahr emittiert, was mit dem C/CO₂-Verhältnis von 0.3 von oben ungefähr 200 Megatonnen Kohlenstoff entspricht, oder irgendwas wie gut 500000 Tonnen am Tag. Damit wäre die Zunahme, die ich hier sehe, rund ein Sechstel des gesamten Kohlenstoffbudgets der BRD, und mehr, wenn der Anstieg schneller als in einem Tag vonstatten geht: Das ist (fast) natürlich Quatsch.

Aber was ist es dann? Noch immer gefällt mir die These ganz lokaler Schwankungen nicht. Wenn hier wirklich nur das CO₂ von Autos und Heizungen nicht mehr weggepustet würde, müsste die Korrelation zwischen CO₂ und Wind viel deutlicher sein.

Ist es eine die Abgasfahne des GKM?

Nächster Versuch: Rund 12 km westlich von meiner Wohnung läuft das Großkraftwerk Mannheim („GKM“). Wenn das Volllast fährt und meine Wohnung in seine Abgasfahne kommt, könnte das so ein Signal geben?

Nun, so ein Kraftwerk liefert ungefähr 1 Gigawatt elektrische Leistung (wie mir der Wikipedia-Artikel dazu gerade verrät: darunter 15% des deutschen Bahnstroms), was bei einem Wirkungsgrad von 1/3 (ok, bei modernen Kohlekraftwerken ist das noch ein wenig mehr, aber als Kopfzahl taugt es) auf 3 Gigawatt thermische Leistung führt (tatsächlich nennt die Wikpedia eine Bruttoleistung von 2146 MW für das GKM).

Aus den 394 kJ/mol, die bei der Verbrennung von Kohlenstoff frei werden (vgl. den Artikel zu meiner thermischen Leistung) könnte mensch jetzt die CO₂-Emission aus der Bruttoleistung ableiten, aber ich bin mal faul und sehe beim WWF nach, der für Kraftwerke dieser Größenordnung ansagt, für eine Kilowattstunde Strom (wir sind dann also wieder bei der Nutzleistung) werde rund ein Kilogramm CO₂ emittiert.

Wenn das Kraftwerk also Volldampf (rund ein GW netto) macht, wird es etwa

10⁹ W⋅0.001 kg ⁄ Wh = 10⁶ kg ⁄ h

CO₂ emittieren, also etwa 1000 Tonnen, was wiederum mit unserem 0.3-Faktor zwischen Kohlenstoff und CO₂ zu einem Kohleverbrauch von 300 Tonnen pro Stunde führt.

Damit leert das Kraftwerk unter Vollast ein Großes Rheinschiff in zehn Stunden – das scheint mir zwar schon sehr schnell zu gehen, ist aber auch nicht gänzlich unplausibel. Gegenrechnung: Das WWF-Dokument von oben nennt 7.7⋅10⁹ kg ⁄ a als CO₂-Emission des GKM im Jahr 2006. Mit der Ur-Kopfzahl π ⋅ 1e7 Sekunden pro Jahr übersetzt sich das in eine mittlere Emission von etwa 200 kg pro Sekunde oder gut 1000 Tonnen pro Stunde. Das passt fast zu gut, denn als jemand, der das Kraftwerk von seiner Leseecke aus sehen kann, kann ich zuverlässig sagen, dass das Ding keineswegs durchläuft. Andererseits hatte das Kraftwerk 2006 auch noch einen Block weniger, und überhaupt ist in der Rechnung genug Luft für Stillstandszeiten.

Nehmen wir …
⎆
Stemming for the Search Engine

2021-11-13
First off, here is a quick reference for the search syntax on this site (the search form links here):
- Phrase searches ("this is a phrase")
- Exclusions (-dontmatch)
- Matches only when two words appear within 10 tokens of each other (matches NEAR appear)
- Trailing wildcard as in file patterns (trail*)
- Searches don't use stemming by default, but stem for German when introduced with l:de and for English when introduced with l:en
- See also the Xapian syntax.
If you only came here for the search syntax, that's it, and you can stop reading here.

Otherwise, if you have read the previous post on my little search engine, you will remember I was a bit unhappy that I completely ignored the language of the posts and had wanted to support stemming so that you can find, ideally, documents containing any of "search", "searches", "searching", and "searched" when searching for any of these. Being able to do that (without completely ruining precision) is obviously language-dependent, which means the first step to make it happen is to properly declare the languague of your posts.

As discussed in the previous post, my blogsearch script only looks at elements with the CSS class indexable, and so I decided to have the language declaration there, too. In my templates, I hence now use:
```
<div class="indexable" lang="{{ article.lang }}">
```
or:
```
<div class="indexable" lang="{{ page.lang }}">
```
as appropriate.

This is interpreted by the indexer rather straightforwardly by pulling the value out of the attribute and asking xapian for a stemmer for the named language. That works for at least most European two-letter country codes, because those happen to coincide with what's legal in HTML's lang universal attribute. It does not work for the more complex BCP 47 language tags like de-AT (where no actually existing stemmer would give results different from plain de anyway) or even sr-Latn-RS (for which, I think, no stemmer exists).

On searching, I was worried that enabling stemming would blow unstemmed searches, but xapian's indexes are clever enough that that's not a problem. But I still cannot stem queries by default, because it is hard to guess their language from just a word or two. Hence, I have defined a query syntax extension: If you prefix your query with l:whatever, blogsearch will try to construct a xapian stemmer from whatever. If that fails, you'll get an error, if it succeeds, it will stem the query in that language.

As an aside, I considered for a moment whether it is a terribly good idea to hand through essentially unfiltered user input to a C++ API like xapian's. I eventually settled for just making it a bit harder to craft buffer overflows by saying:
```
lang = parts[0][2:30]
```
– that is, I'm only allowing through up to 28 characters of language code. Not that I expect that anything in between my code and xapian's core has an overflow problem, but this is a cheap defensive measure that would also limit the amount of code someone could smuggle in in case some vulnerability did sneak in. Since it's essentially free, I'd say that's reasonable defensive programming.

In closing, I do not think stemmed searches will be used a lot, and as usual with these very simple stemmers, they leave a lot to be desired from a linguistic point of view. Compare, for instance, a simple search for going with the result l:en going to see where this is supposed to go (and compare with the result when stemming as German). And then compare with l:en went, which should return the same as l:en going in an ideal world but of course doesn't: Not with the simple snowball stemmer that xapian employs.

I'm still happy the feature's there, and I'm sure I'll need it one of these days.

And again, if you need a CGI that can index and query your static HTML collection with low deployment effort: you're welcome.
⎆
Der hundertste Post

2021-11-10
Vor 10 Monaten habe ich den ersten Artikel für dieses Blog geschrieben, und siehe da: Mit diesem sind es jetzt 100 Posts geworden.

Das wäre ein guter Vorwand für ein paar Statistiken, aber da ich ja generell ein Feind von Metriken bin, die mensch ohne konkrete Fragestellung sammelt (das ist ein wenig wie beim statistischen Testen: Wenn du nicht von vorneherein weißt, worauf du testest, machst du es falsch), bestätige ich mir nur, dass meine Posts viel länger sind als ich das eigentlich will. Insgesamt nämlich habe ich nach Zählung von wc -l auf den Quelldateien fast 93000 Wörter in diesen Artikeln. Zur Fehlerabschätzung: xapian (vgl. unten) zählt nur 89000.

Die Länge der Artikel ist nach wc-Wörtern so verteilt:

Ich weiß auch nicht recht, warum ich mich nicht kürzer fassen kann. Oder will. Der überlange Post mit 3244 Wörtern ist übrigens der über die Konfiguration eines Mailservers – und das ist wieder ein gutes Beispiel für die Fragwürdigkeit von Metriken, denn erstens hat Englisch fast keine Komposita und ist von daher im Nachteil beim Wörterzählen und zweitens ist in dem Artikel ziemlich viel Material, das in Wirklichkeit Rechner lesen, und das sollte wirklich anders zählen als natürlichsprachiger Text.

Na gut, und einem Weiteren kann ich nicht widerstehen: Wie viele verschiedene Wörter („Paradigmata“) kommen da eigentlich vor? Das ist natürlich auch Mumpitz, denn die Definition, wann zwei Wörter verschieden sind („die Token verschiedenen Paradigmata angehören“), ist alles andere als tivial. So würde ich beispielsweise behaupten, dass die Wörter Worte und Wörter praktisch nichts miteinander zu tun haben, während im Deuschen z.B. auf, schaute und aufschauen besser alle zusammen ein einziges Paradigma bilden sollten (zusammen mit allerlei anderem).

Aber ist ja egal, sind ja nur Metriken, ist also eh Quatsch. Und es gibt die Daten auch schon, was für die Nutzung von und die Liebe zu Kennzahlen immer ein Vorteil ist. Ich habe nämlich den xapian-Index über dem Blog, und mit dem kann ich einfach ein paar Zeilen Python schreiben:
```
import xapian
db = xapian.Database("output/.xapian_db")
print(sum(1 for w in db.allterms()))
```
(Beachtet die elegante Längenbestimmung mit konstantem Speicherbedarf – db.allterms() ist nämlich ein Iterator).

Damit bekomme ich – ich stemme nach wie vor nicht – 16540 raus. Klar, diese 16540 für die Zahl der verschiedenen Wörter ist selbst nach den lockeren Maßstäben von Metriken ganz besonders sinnlos, weil es ja eine wilde Mischung von Deutsch und Englisch ist.

Um so mehr Spaß macht es, das mit den 100'000 Wörtern zu vergleichen, die schließlich mal im Goethe-Wörterbuch sein sollen, wenn es fertig ist. Eine schnelle Webrecherche hat leider nichts zur Frage ergeben, wie entsprechende Schätzungen für Thomas Mann aussehen. Einmal, wenn ich gerne Kennzahlen vergleichen würde…
⎆
Halbwegs gute Nachrichten

2021-11-06
Aus dieser Kurve lässt sich ablesen, dass wir in acht Tagen knapp 3500 SARS-2-Fälle auf Intensivstationen haben werden. Wie, verrate ich in diesem Artikel.
In meinen Corona-Überlegungen gestern habe ich mich mal wieder gefragt, wie lange wohl die „Intensiv-Antwort“ dem Inzidenzsignal nachläuft, wie viele Tage es also dauert, bis sich ein Anstieg in den Inzidenzen in der Intensivbelegung reflektiert. Diese Frage ist, wie ich unten ausführe, derzeit ziemlich relevant im Hinblick auf Überlegungen, wie lange wir eigentlich noch Zeit haben, um massenhafte Triage in (oder vor) unseren Intensivstationen abzuwenden.

Meine Null-Annahme für diese Verzögerung war seit Mai 2020 – nach einer entsprechenden Ansage eines befreundeten Anästhesisten – „eher so drei Wochen“. Bei nährem Nachdenken ist mir gestern aber aufgefallen, dass er wohl eher „von Infektion bis Intensiv“ gemeint haben wird, und dann ist die Antwort sicher schneller, denn von Infektion bis Meldung vergeht normalerweise wohl mindestens eine Woche. Aber: Ich muss nicht raten. Wir spielen hier mit Kennzahlen, die vielleicht in der Realität nicht immer viel bedeuten, aber zumindest klar definiert sind. Daher lässt sich der Verzug auf der Basis von RKI- und DIVI-Zahlen nachrechnen.

Ich verrate gleich mal das Ergebnis: ich komme für die zweite Welle auf gut fünf Tage Verzug der Intensiv-Antwort, für die dritte auf etwa sieben Tage, für die vierte Welle auf acht bis neun Tage. Wie ich unten ausführe, sind das relativ gute Nachrichten.
Zeitreihen

Wie habe ich das gerechnet? Nun, die Korrelation von Zeitreihen ist ein ganzer Satz von Wissenschaften, bei denen es meist darum geht, ungleiche Abdeckungen, unregelmäßig gesetzte Messpunkte sowie allerlei Rauschen und Schmutz weggefummelt zu kriegen, ohne allzu viele Informationen zu verlieren oder, vielleicht schlimmer, damit Artefakte einzubauen.

Die so gereinigten Daten lassen sich dann zum Beispiel geeignet skaliert und verschoben übereinanderlegen. Tatsächlich sind die Parameter dieser Transformationen im Regelfall (und gewissermaßen auch hier) viel interessanter als die Zeitreihen selbst[1]. Es gibt daher zahlreiche mathematische Verfahren, die das von einer Augenmaß-Übung in etwas verwandeln, das reproduzierbar und auch quantifizierbar ist. Vorsicht: ich bin da kein Experte und gehe hier nur mit nicht allzu schwer erkranktem Menschenverstand ran, verwende also (zumindest in Summe) gerade kein wirklich wohldurchdachtes Verfahren. Wer sowas wie das hier für Hausaufgaben oder Hausarbeiten verwendet, tut das auf eigene Gefahr.

Andererseits bin ich recht zuversichtlich, dass der Kram insgesamt schon stimmt.

Ausgangsdaten sind meine aus RKI-Berichten gescrapten Intensivbelegungen und ein RKI-Sheet (in, ach weh, „Office Open XML“ a.k.a. XSLX – muss das sein?) mit den Inzidenzen. Dabei ist das erste Problem, dass ich aus verschiedenen Gründen nicht für jeden Tag Belegungszahlen habe, und so ist mein erster Schritt, fehlende Punkte zu interpolieren. Das Scipy-Paket macht es leicht, aus einem Satz von Zeit/Wert-Paaren (die Zeit wird in den Lesefunktionen auf Sekunden seit einer Epoche gewandelt) ein handliches Array zu rechnen:
```
grid_points = numpy.arange(
  raw_crits[0][0],
  raw_crits[-1][0],
  86400)  # ein Tag in Sekunden
critnums = interpolate.griddata(
  raw_crits[:,0], raw_crits[:,1], grid_points)
```
Ich wollte diese Interpolation eigentlich visualisieren, habe aber keine gute Stelle gefunden, an der sie einen sichtbaren Unterschied gemacht hätte, und entscheidend ist eigentlich nur, dass ich ab diesem Schritt blind mit Arrays arbeiten kann; deren Index ist zunächst die Zahl der Tage seit dem ersten Tag mit Daten, hier speziell dem 11.8.2020, denn damals habe ich mit dem Screenscrapen der DIVI-Daten angefangen.
Übergeplottete Kurven

Wenn ich diese interpolierten Kurven übereinanderlege, kommt das hier heraus:

Weil das, wie gesagt, erst im August 2020 anfängt, fehlt die erste Welle.

Das bloße Auge reicht für die Bestätigung der Erwartung, dass die Intensivbelegung der Inzidenzkurve meist etwas hinterherläuft, wenn auch nicht so, dass mensch hoffen könnte, die beiden durch etwas Schieben global übereinanderzubekommen. Die großen Zacken in den Inzidenzen durch Weihnachten und Ostern finden sich in der Intensivbelegung gar nicht wieder, was ein klares Zeichen ist, dass sie weitgehend Erfassungsartefakte sind. In der Hinsicht wären die „großen“ RKI-Zahlen mit Referenzdaten (vgl. die Film-Geschichte) bestimmt besser, aber ich wollte für dieses Ding nicht die 200 Megabyte durchkämmen, zumal die „kleinen“ RKI-Daten besser auf die Meldedaten aus den Tagesberichten passen, um die es mir hier ja geht.

Vor allem fällt auf, dass meine Jammerei darüber, dass der Impffortschritt die Intensivantwort nicht wesentlich abgeflacht hat, unzutreffend ist: Mit der Skalierung aus dem Plot liegen Inzidenz und Belegung in der zweiten und dritten Welle ziemlich übereinander, während in der vierten Welle doch ein knapper Faktor zwei dazwischenliegt; da scheint die Impfung doch ein wenig gegenüber Delta zu gewinnen (aber, klar, nirgendwo hinreichend).

Ich hatte erwartet, dass die abfallenden Flanken der Intensivkurven deutlich flacher sind als die der Inzidenzkurven, weil Leute unter Umständen lang auf Intensiv liegen und es entsprechend lang dauern sollte, bis sich die Stationen wieder leeren. Das sieht im Abfall der zweiten Welle auch ein wenig so aus, nicht jedoch bei dem der dritten Welle. Bei ihr fällt die Intensivbelegung sehr treu mit der Inzidenz. So viele LangzeitpatientInnen gibt es glücklicherweise wohl doch nicht.
Nach Glätten differenzierbar

Wie kann ich jetzt den Verzug zu quantifizieren? Mein (wie gesagt eher intuitiv gefasster) Plan ist, über die Ableitung der jeweiligen Kurven zu gehen, und zwar aus der Überlegung heraus, dass mich ja Veränderungen viel mehr als Pegel interessieren. Nun habe ich aber keine differenzierbaren Funktionen, sondern lediglich Arrays, bei denen ich Ableitungen allenfalls durch Subtrahieren benachbarter Elemente simulieren kann. Solche numerischen „Ableitungen“ reagieren ziemlich empfindlich auf das Gewackel („Rauschen“), das es in realen Daten immer gibt („Subtraktion ist in der Regel numerisch schlecht konditioniert“). Deshalb will ich meine rohen Daten glätten, bevor ich die „Ableitung“ ausrechne, sprich: das Rauschen rausnehmen, ohne das Signal wesentlich zu verzerren.

Glättung heißt eigentlich immer, Kurvenpunkte in einem Zelle für Zelle über die Daten laufenden Fenster zu mitteln, also z.B., indem mensch je fünf Nachbarpunkte rechts und links auf den aktuellen Punkt draufaddiert und ihn dann durch die durch elf geteilte Summe ersetzt. Mit diesem ganz naiven Rezept bekommen relativ weit entfernte Werte aber genauso viel Einfluss auf den aktuellen Punkt wie die unmittelbaren Nachbarn. Das modelliert meist die sachlichen Grundlagen des Rauschens nur schlecht. Es ist auch aus theoretischeren Gründen normalerweise eher ungünstig.

Der eher theoretische Hintergrund ist grob, dass bei der Glättung letztlich zwei Funktionen gefaltet werden. Das ist äquivalent dazu, die Spektren der beiden Signale (ihre Fouriertransformierten) zu multiplizieren und dann wieder zurückzutransformieren. Bei einem einfachen Fenster des oben skizierten Typs („Rechteckfunktion“) gibts nun sehr scharfe Kanten, die wiederum zu einem sehr breiten Spektrum führen, so dass auch das Spektrum der geglätteten Funktion allerlei unwillkommene Features bekommen kann (manchmal ist das aber auch genau das, was mensch haben will – hier nicht).

Wie auch immer: scipy macht es einfach, mit einer lärmarmen Funktion – die sieht ein wenig gaußig aus und ist im nächsten Plot im kleinen Inset zu sehen – zu glätten, nämlich etwa so:
```
import numpy
from scipy import signal

smoothing_kernel = numpy.kaiser(20, smoothing_width)
smoothing_kernel = smoothing_kernel/sum(smoothing_kernel)
convolved = signal.convolve(arr, smoothing_kernel, mode="same"
  )[smoothing_width:-smoothing_width]
```
Die Division durch die Summe von smoothing_kernel in der zweiten Zeile sorgt dafür, dass sich an der „Höhe“ der geglätteten Funktion insgesamt nichts ändert: Sozusagen ausgeklammert ist die Faltung eine Multiplikation mit eins. Das Wegschneiden der Ränder in der letzten Zeile wiederum entfernt Punkte, bei denen fehlende Werte am Anfang und Ende der Zeitreihe im Fenster waren. Scipy ersetzt die durch Nullen, so dass die geglätteten Werte am Rand steil abfallen. Was ich hier mache, entspricht technisch dem valid-Mode der convolve-Funktion. Nur weiß ich hier zuverlässig, wie lang das Array am Schluss ist.

Der Effekt, hier auf Indzidenzdaten der vierten Welle:

Damit kann ich jetzt meine numerische „Ableitung“ bilden, ohne dass mir das Ergebnis furchtbar rauscht:
```
diff = convolved[1:]-convolved[:-1]
```
Die nächsten beiden Grafiken zeigen diese Pseudo-Ableitungen für Inzidenz und Intensivbelegung. Ich habe jeweils in blau reingemalt, wie es ohne Glättung aussehen würde, um deutlich zu machen, warum diese eine gute Idee ist und was ich mit „furchtbar rauschen“ meine. Weiterverwendet werden natürlich die orangen Verläufe:

In der Inzidenzkurve fallen wieder Weihnachten und Ostern besonders auf, weil sie selbst in den geglätteten Graphen noch wilde Ausschläge verursachen.
In der Zeit verschieben

Schon der optische Eindruck aus dem Rohdaten-Plot legt nahe, die einzelnen Wellen getrennt zu untersuchen, und das ist angesichts von über die Zeit stark veränderlichen Viren, Testverhältnissen, demographischen Gegebenheiten und nichtpharamzeutischen Maßnahmen sicher auch sachlich geboten.

Deshalb hier zunächst der Verlauf der Ableitungen von Inzidenzen und Intensivbelegung in der aktuellen vierten Welle (das ist eine Kombination der rechten Enden der orangen Graphen der letzten beiden Plots):

In dem Graphen steckt noch eine weitere Normalisierung. Und zwar habe ich für beide Arrays etwas wie:
```
incs /= sum(abs(incs))
```
laufen lassen. Damit ist die Fläche zwischen den Kurven …
⎆
A Local Search Engine for Pelican-based Blogs

2021-10-31
As the number of posts on this blog approaches 100, I figured some sort of search functionality would be in order. And since I'm wary of “free” commercial services and Free network search does not seem to go anywhere[1], the only way to offer that that is both practical and respectful of the digital rights of my readers is to have a local search engine. True, having a search engine running somewhat defeats the purpose of a static blog, except that there's a lot less code necessary for doing a simple search than for running a CMS, and of course you still get to version-control your posts.

I have to admit that the “less code” argument is a bit relative given that I'm using xapian as a full-text indexer here. But I've long wanted to play with it, and it seems reasonably well-written and well-maintained. I have hence written a little CGI script enabling search over static collections of HTML files, which means in particular pelican blogs. In this post, I'll tell you first a few things about how this is written and then how you'd run it yourself.
Using Xapian: Indexing

At its core, xapian is not much more than an inverted index: Essentially, you feed it words (“tokens”), and it will generate a database pointing from each word to the documents that contain it.

The first thing to understand when using xapian is that it doesn't really have a model of what exactly a document is; the example indexer code, for instance, indexes a text file such that each paragraph is treated as a separate document. All xapian itself cares about is a string („data“, but usually rather metadata) that you associate with a bunch of tokens. This pair receives a numeric id, and that's it.

There is a higher-level thing called omega built on top of xapian that does identify files with xapian documents and can crawl and index a whole directory tree. It also knows (to some extent) how to pull tokens from a large variety of file types. I've tried it, and I wasn't happy; since pelican creates all those ancillary HTML files for tags, monthly archives, and whatnot, when indexing with omega, you get lots of really spurious matches as soon as people enter a term that's in an article title, and entering a tag or a category will yield almost all the files.

So, I decided to write my own indexer, also with a view to later extending it to language detection (this blog has articles in German and English, and they eventually should be treated differently). The core is rather plain in Python:
```
for dir, children, names in os.walk(document_dir):
  for name in fnmatch.filter(names, "*.html"):
    path = os.path.join(dir, name)
    doc = index_html(indexer, path, document_dir)
```
That's enough for iterating over all HTML files in a pelican output directory (which document_dir should point to).

In the code, there's a bit of additional logic in the do_index function. This code enables incremental indexing, i.e., only re-indexing a file if it has changed since the last indexing run (pelican fortunately manages the file timestamps properly).

Nachtrag (2021-11-13)

It didn't, actually; see the search engine update post for how to fix that.

What I had to learn the hard way is that since xapian has no built-in relationship between what it considers a document and an operating system file, I need to explicitly remove the previous document matching a particular file. The function get_indexed_paths produces a suitable data structure for that from an existing database.

The indexing also defines my document model; as said above, as far as xapian is concerned, a document is just some (typically metadata) string under user control (plus the id and the tokens, obviously). Since I want structured metadata, I need to structure that string, and these days, json is the least involved thing to have structured data in a flat string. That explains the first half of the function that actually indexes one single document, the path of which comes in in f_name:
```
def index_html(indexer, f_name, document_dir):
  with open(f_name, encoding="utf-8") as f:
    soup = bs4.BeautifulSoup(f, "lxml")
  doc = xapian.Document()
  meta = {
    "title": soup_to_text(soup.find("title")),
    "path": remove_prefix(f_name, document_dir),
    "mtime": os.path.getmtime(f_name),}
  doc.set_data(json.dumps(meta))

  content = soup.find(class_="indexable")
  if not content:
    # only add terms if this isn't some index file or similar
    return doc
  print(f"Adding/updating {meta['path']}")

  indexer.set_document(doc)
  indexer.index_text(soup_to_text(content))

  return doc
```
– my metadata thus consists of a title, a path relative to pelican's output directory, and the last modification time of the file.

The other tricky part in here is that I only index children of the first element with an indexable class in the document. That's the key to keeping out all the tags, archive, and category files that pelican generates. But it means you will have to touch your templates if you want to adopt this to your pelican installation (see below). All other files are entered into the database, too, in order to avoid needlessly re-scanning them, but no tokens are associated with them, and hence they will never match a useful query.

Nachtrag (2021-11-13)

When you add the indexable class to your, also declare the language in order to support stemming; this would look like lang="{{ page.lang }} (substituting article for page as appropriate).

There is a big lacuna here: the recall, i.e., the ratio between the number of documents actually returned for a query and the number of documents that should (in some sense) match, really suffers in both German and English if you don't do stemming, i.e., fail to strip off grammatical suffixes from words.

Stemming is of course highly language-dependent. Fortunately, pelican's default metadata includes the language. Less fortunately, my templates don't communicate that metadata yet – but that would be quick to fix. The actual problem is that when I stem my documents, I'll also have to stem the incoming queries. Will I stem them for German or for English?

I'll think about that problem later and for now don't stem at all; if you remember that I don't stem, you can simply append an asterisk to your search term; that's not exactly the same thing, but ought to be good enough in many cases.
Using xapian: Searching

Running searches using xapian is relatively straightforward: You open the database, parse the query, get the set of matches and then format the metadata you put in during indexing into links to the matches. In the code, that's in cgi_main; one could do paging here, but I figure spitting out 100 matches will be plenty, and distributing 100 matches on multiple HTML pages is silly (unless you're trying to optimise your access statistics; since I don't take those, that doesn't apply to me).

The part with the query parser deserves a second look, because xapian supports a fairly rich query language, where I consider the most useful features:
- Phrase searches ("this is a phrase")
- Exclusions (-dontmatch)
- Matches only when two words appear within 10 tokens of each other (matches NEAR appear)
- Trailing wildcard as in file patterns (trail*)
That last feature needs to be explicitly enabled, and since I find it somewhat unexpected that keyword arguments are not supported here, and perhaps even that the flag constant sits on the QueryParser object, here's how enabling wildcards in xapian looks in code:
```
qp = xapian.QueryParser()
parsed = qp.parse_query(query, qp.FLAG_WILDCARD)
```
Deploying this on your Pelican Installation

You can re-use my search script on your site relatively easily. It's one file, and if you're running an apache or something else that can run CGIs[2], making it run first is close to trivial: Install your equivalents of the Debian python3-xapian, python3-bs4, and python3-lxml packages. Perhaps you also need to explicitly allow CGI execution on your web server. In Debian's apache, that would be a2enmod cgi, elsewhere, you may need to otherwise arrange for mod_cgi or its equivalent to be loaded.

Then you need to dump blogsearch somewhere in the file system.

Nachtrag (2022-10-07)

Don't take it from here; rather, see https://codeberg.org/AnselmF/pelican-ext

While Debian has a default CGI directory defined, I'd suggest to put blogsearch somewhere next to your blog; I keep everything together in /var/blog (say), have the generated output in /var/blog/generated and would then keep the script in a directory /var/blog/cgi. Assuming this and apache, You'd then have something like:
```
DocumentRoot /var/blog/generated
ScriptAlias /bin /var/blog/cgi
```
in your configuration, presumably in a VirtualHost definition. In addition, you will have to tell the script where your pelican directory is. It expects that information in the environment variable BLOG_DIR; so, for apache, add:
```
SetEnv BLOG_DIR /var/blog/generated
```
to the VirtualHost.

After restarting your web server, the script would be ready (with the configuration above …
⎆
Math with ReStructuredText and Pelican

2021-10-24
I recently wrote a piece on estimating my power output from CO₂ measurements (in German) and for the first time in this blog needed to write at least some not entirely trivial math. Well: I was seriously unhappy with the way formulae came out.

Ugly math of course is very common as soon as you leave the lofty realms of LaTeX. This blog is made with ReStructuredText (RST) in pelican. Now, RST at least supports the math interpreted text role (“inline”) and directive (“block“ or in this case rather “displayed“) out of the box. To my great delight, the input syntax is a subset of LaTeX's, which remains the least cumbersome way to input typeset math into a computer.

But as I said, once I saw how the formulae came out in the browser, my satifsfaction went away: there was really bad spacing, fractions weren't there, and things were really hard to read.

In consequence, when writing the post I'm citing above, rather than reading the docutils documentation to research whether the ugly rendering was a bug or a non-feature, I wrote a footnote:

Sorry für die hässlichen Formeln. Vielleicht schreibe ich mal eine Erweiterung für ReStructuredText, die die ordentlich mit TeX formatiert. Oder zumindest mit MathML. Bis dahin: Danke für euer Verständnis.

(Sorry for the ugly formulae. Perhaps one of these days I'll write an RST extension that properly formats using TeX. Or at least MathML. Until then: thanks for your understanding.)

This is while the documentation clearly said, just two lines below the example that was all I had initially bothered to look at:

For HTML, the math_output configuration setting (or the corresponding --math-output command line option) selects between alternative output formats with different subsets of supported elements.

Following the link at least would have told me that MathML was already there, saving me some public embarrassment.

Anyway, when yesterday I thought I might as well have a look at whether someone had already written any of the code I was talking about in the footnote, rather than properly reading the documentation I started operating search engines (shame on me).

Only when those lead me to various sphinx and pelican extensions and I peeked into their source code I finally ended up at the docutils documentation again. And I noticed that the default math rendering was so ugly just because I didn't bother to include the math.css stylesheet. Oh, the miracles of reading documentation!

With this, the default math rendering suddenly turns from ”ouch” to “might just do”.

But since I now had seen that docutils supports MathML, and since I have wanted to have a look at it at various times in the past 20 years, I thought I might as well try it, too. It is fairly straightforward to turn it on; just say:
```
[html writers]
math_output: MathML
```
in your ~/.docutils (or perhaps via a pelican plugin).

I have to say I am rather underwhelmed by how my webkit renders it. Here's what the plain docutils stylesheet works out to in my current luakit:

$Screenshot with ok formulae.$

And here's how it looks like via MathML:

$Screenshot with less ok formulae.$

For my tastes, the spacing is quite a bit worse in the MathML case; additionally, the Wikipedia article on MathML mentions that the Internet Explorer never supported it (which perhaps wouldn't bother me too much) and that Chromium withdrew support at some point (what?). Anyway: plain docutils with the proper css is the clear winner here in my book.

I've not evaluated mathjax, which is another option in docutils math_output and is what pelican's render_math plugin uses. Call me a luddite, but I'll file requiring people to let me execute almost arbitrary code on their box just so they see math into the big folder labelled “insanities of the modern Web”.

So, I can't really tell whether mathjax would approach TeX's quality, but the other two options clearly lose out against real TeX, which using dvipng would render the example to:

$Screenshot with perfect formulae$

– the spacing is perfect, though of course the inline equation has a terrible break (which is not TeX's fault). It hence might still be worth hacking a pelican extension that collects all formulae, returns placeholder image links for them and then finally does a big dvipng run to create these images. But then this will mean dealing with a lot of files, which I'm not wild about.

What I'd like to ideally use for the small PNGs we are talking about here would be inline images using the data scheme, as in:
```
<img src="data:image/png;base64,AAA..."/>
```
But since I would need to create the data string when docutils calls my extension function, I in that scheme cannot collect all the math rendering for a single run of LaTeX and dvipng. That in turn would mean either creating a new process for TeX and dvipng each for each piece of math, which really sounds bad, or hacking some wild pipeline involving both, which doesn't sound like a terribly viable proposition either.

While considering this, I remembered that matplotlib renders quite a bit of TeX math strings, too, and it lets me render them without any fiddling with external executables. So, I whipped up this piece of Python:
```
import base64
import io
import matplotlib
from matplotlib import mathtext

matplotlib.rcParams["mathtext.fontset"] = "cm"

def render_math(tex_fragment):
    """returns self-contained HTML for a fragment of TeX (inline) math.
    """
    res = io.BytesIO()
    mathtext.math_to_image(f"${tex_fragment}$",
      res, dpi=100, format="png")
    encoded = base64.b64encode(res.getvalue()).decode("ascii")
    return (f'<img src="data:image/png;base64,{encoded}"'
        f' alt="{tex_fragment}" class="math-png"/>')

if __name__=="__main__":
    print(render_math("\int_0^\infty \sin(x)^2\,dx"))
```
This prints the HTML with the inline formula, which with the example provided looks like this: $\int_0^\infty \sin(x)^2\,dx$ – ok, there's a bit too much cropping, I'd have to trick in transparency, there's no displayed styles as far as I can tell, and clearly one would have to think hard about CSS rules to make plausible choices for scale and baseline – but in case my current half-satisfaction with docutils' text choices wears off: This is what I will try to use in a docutils extension.
⎆

Seite 1 / 2 »

Tag Python

The Good

The Bad

The Crazy

Erstmal glätten

Zeitrechnungen

Spektralanalyse

Preparations

Figuring out Libreoffice's API

The Wikipedia-Opening Script

Codeberg

Insights into Pelican

Wie machte ich aus einem Satz von ODTs oder DOCs ein einziges ODT?

epubedit

Stilfragen

Und die Stichworte?

Security Considerations

Bürgis geniale Methode

Ein numerischer Traum

Technics

Technics

Sicherheitsfragen

Verheizt jemand 100000 Tonnen Kohlenstoff am Tag?

Andersrum gerechnet

Ist es eine die Abgasfahne des GKM?

Zeitreihen

Übergeplottete Kurven

Nach Glätten differenzierbar

In der Zeit verschieben

Using Xapian: Indexing

Using xapian: Searching

Deploying this on your Pelican Installation

Letzte Ergänzungen