Freitag, 26. Februar 2016

Ende

Mit diesem Post ist dieser Blog offiziell beendet. Die Swiss Ephemeris unterstützt seit der Version 2.03 das Multithreading.

Die Implementierungsidee war sehr einfach und kommt ganz ohne das in diesem Blog skizzierte grundsätzliche Refactoring aus: Globale Variablen, die pro Thread separat behandelt werden müssen, bekommen einfach bei der Deklaration den Storage class specifier _thread und gehören somit zur Speicherklasse "TLS" (für thread-local storage), der von C-Compilern schon seit gut anderthalb Jahrzehnten für die gängigen Zielarchitekturen unterstützt wird. Daten von dieser Speicherklasse werden bei der Erzeugung eines neuen Threads dupliziert, so dass der Thread über eine eigene Version von ihnen verfügt und sich beim Ändern der Inhalte nicht mit anderen Threads in die Quere kommt. In der Swiss Ephemeris war das im wesentlichen die globale Variable swed - eine Sammelstruktur, die als Zwischenspeicher für alle Art von Daten verwendet wird, die den einzelnen Aufruf einer Funktion der Swiss Ephemeris überleben sollen.

Wäre ich selbst mehr als nur freizeitmässig dem Thema zugeneigt gewesen, hätte mich schon 2013 die folgende Antwort des Users skylendar in der Swiss Ephemeris-Mailingliste nachdenklich stimmen müssen. Ich hatte gefragt, ob sich schon einmal andere Entwickler dieses Themas angenommen hätten und bekam von jenem skylendar die lapidare Antwort, er hätte schon vor langer Zeit eine mehrfädige Version der Swiss Ephemeris geschrieben, diese jedoch mangels Interesse in der Mailingliste nicht weiter verbreitet. Da ich damals nur den Lösungsweg gesehen hatte, das Programm total umzubauen, um die globalen Variablen aus dem Code zu eliminieren, fragte ich ihn, wie er dieses Kunststück vollbracht habe. Seine - wiederum lapidare - Antwort hätte mich damals nachdenklich stimmen müssen. Stattdessen überlas ich sie:

Delivered-To: ruediger.plantiko@gmail.com
Received: by 10.194.135.203 with SMTP id pu11csp79435wjb;
        Fri, 1 Mar 2013 08:04:14 -0800 (PST)
Message-ID: 
In-Reply-To: 
User-Agent: eGroups-EW/0.82
From: "skylendar" 
Sender: swisseph@yahoogroups.com
Mailing-List: list swisseph@yahoogroups.com; contact swisseph-owner@yahoogroups.com
Delivered-To: mailing list swisseph@yahoogroups.com
Date: Fri, 01 Mar 2013 16:04:07 -0000
Subject: [SWISSEPH] Re: Is libswe thread-safe?

--- In swisseph@yahoogroups.com, Rüdiger Plantiko  wrote:
>
> Skylendar,
>
> I missed this. Did you rewrite
> sweph.c?
> If not, how did you avoid the usage of global and static variables?
>
> Regards,
> Rüdiger

I used the thread keyword wherever it is necessary.

Da das Multithreading in der Mailingliste aber auch weiterhin immer wieder zum Thema wurde, postete skylendar am 4. September 2015 schliesslich doch seinen Patch, den er damals noch auf der Version 1.71 aufgesetzt hatte.

Es zeigte sich: wenn man diesen Weg verfolgt, bringt man sich zwar um eine an sich gute Gelegenheit, seinen Code zu refaktorisieren, aber man kann mit vergleichsweise wenigen Änderungen die Mehrfädigkeit erreichen.

Ich folgte seinem Weg für das damals aktuelle Release 2.02.01:

  • Ich änderte alle globalen Variablen in const, auf die im Programm nur lesend zugegriffen wurde. Davon gab es eine ganze Reihe (etwa Koeffizienten oder sonstige Konstanten aus Reihenentwicklungen).
  • Alle globalen Variablen, die nicht nur gelesen, sondern auch verändert wurden, bekamen das neue Attribut TLS (ein Makro, hinter dem sich das Schlüsselwort __thread verbarg.
  • Dasselbe für die zahlreichen statischen Variablen, die sich ja nur durch die eingeschränkte Sichtbarkeit (File oder Function) von den globalen unterscheiden, aber derselben Verletzlichkeit bei mehrfädiger Programmausführung unterliegen wie globale Variablen.
Mit folgendem kleinen C99-Testprogramm, das die POSIX Schnittstelle <pthread.h> verwendet, konnte ich auf meinem Rechner mit 2*4 Prozessorkernen eine Laufzeitverbesserung um den Faktor 5 für die zufällige Berechnung von Planetenpositionen erreichen. Wie ich aber schon früher hier begründete, geht es nicht nur um die Laufzeitverbesserung durch effizientere Ressourcennutzung. Ein anderer Gewinn ist, dass mehrfädige Clientprogramme (und Plattformen wie z.B. mod_perl oder Java Server Pages) die Swiss Ephemeris nur dann sicher und ohne Extraprogrammierung zur Synchronisierung (etwa mit Mutexen) benutzen können, wenn diese selbst mehrfägig ist. Daher erweitert die Mehrfädigkeit auch die Einsatzmöglichkeiten der Swiss Ephemeris.
/*
 * Test program for multithreading
 * Performs many swe_calc() executions in several threads simultaneously
 *
 * Syntax: 'test <number of threads> [ <calcs per thread> [ <repetitions> [ <trace mode> ] ] ]'
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <sys/time.h>
#include <stdbool.h>
#include "swephexp.h"
#include "sweph.h"

typedef struct {
  double jd,l;
  int ipl;
  bool ok;
  } calc_data;

typedef struct {
  calc_data *m;
  int length;
} calc_data_array;

double get_time();
void *get_longitude( void* arg );
void *get_longitudes( void* arg );
void prepare_random_arguments( int nthreads, int calcs_per_thread, calc_data_array* data );
void parse_cmdline_args( int argc, char* argv[], int* nthreads, int *calcs_per_thread, int* repetitions, bool* trc);
void print_trace( int index,int nthreads, int calcs_per_thread, calc_data_array *data);


int main( int argc, char* argv[] ) {
  
  int nthreads         = 1;
  int calcs_per_thread = 1;
  int repetitions      = 100;
  bool trc             = false;
  double total_time    = 0;

// Inititialize pseudo random generator with system time as seed
  srand( time(NULL) );

  parse_cmdline_args( argc, argv, &nthreads, &calcs_per_thread, &repetitions, &trc);
 
// Allocate calc data array on stack
  calc_data _data[nthreads][calcs_per_thread];
  calc_data_array data[nthreads];
  for (int i=0;i<nthreads;i++) {
    data[i].m      = _data[i];
    data[i].length = calcs_per_thread;
  }

  pthread_t tid[nthreads];
 
  for (int i=0; i<repetitions; i++ ) {

// Prepare random argument chunks for swe_calc
    prepare_random_arguments( nthreads, calcs_per_thread, data );

// Compute a chunk of nthread*calcs_per_thread data
    double begin = get_time( );
    if (nthreads == 1) {
      // Direct call for comparison
      get_longitudes( &data[0] );
      }
    else {
      // ... or multithreaded
      for (int j=0;j<nthreads;j++) {
        pthread_create( &tid[j], NULL, get_longitudes, (void *) &data[j] );
      }
      for (int j=0;j<nthreads;j++) {
        pthread_join( tid[j], NULL);
      }
    }
    double end = get_time();
    total_time += (end - begin);

// Details, if wanted
    if (trc) print_trace(i,nthreads,calcs_per_thread,data);

  }  

// Execution times summary
  printf( 
      "%d times (%d threads, in chunks of %d), "
      "done in %lf msec total (=%lf /each) \n",
      repetitions,
      nthreads,
      calcs_per_thread,
      total_time*1000,
      total_time*1000/(repetitions*nthreads*calcs_per_thread)
      );

} 

// The "worker" function for a chunk of calculations, executed simultaneously in several threads
void *get_longitudes( void* arg) {
  calc_data_array *data = (calc_data_array*) arg;
  for (int i=0;i<data->length;i++) {
    get_longitude(&data->m[i]);
  }
  return NULL;
}

// The "worker" function for a single calculation
void *get_longitude( void* arg ) {
  calc_data *data = (calc_data*) arg;
  char serr[255];
  int  iflag = 0, iret = 0;
  double xx[6] = {};
  iret = swe_calc( data->jd, data->ipl, iflag, xx, serr );
  data->ok = (iret == iflag);
  data->l = xx[0];
  return NULL;
  }


// Use pseudo random generator to generate some Julian dates and planet numbers
void prepare_random_arguments( int nthreads, int calcs_per_thread, calc_data_array* data ) {
  for (int i=0; i<nthreads;i++)
    for (int j=0; j<calcs_per_thread;j++){
      data[i].m[j].jd  = J1900 + (rand() % 10000); // range: 10000 days from 1.1.1900 onwards
      data[i].m[j].ipl = rand( ) % 10;
    }
}


void parse_cmdline_args( int argc, char* argv[], int* nthreads, int *calcs_per_thread, int* repetitions, bool* trc) {

  if (argc == 1) {
    printf( "Syntax: 'test <number of threads> [ <calcs per thread> [ <repetitions> [ <trace mode> ] ] ]'\n");
    exit(4);
  }
  if (argc > 1) sscanf( argv[1], "%d", nthreads    );
  if (argc > 2) sscanf( argv[2], "%d", calcs_per_thread );
  if (argc > 3) sscanf( argv[3], "%d", repetitions );
  if (argc > 4) *trc = true; // *Any* fourth argument activates trace mode

}

double get_time() {
  struct timeval tv;
  gettimeofday(&tv,NULL);
  return tv.tv_sec + tv.tv_usec / 1000000.0 ;
  }

void print_trace( int index, int nthreads, int calcs_per_thread, calc_data_array *data) {
for (int k=0,kmax = nthreads*calcs_per_thread;k<kmax;k++) {
  int i = k / calcs_per_thread,
      j = k % calcs_per_thread;
  printf( "%d / %d / %d : %lf %d -> %lf %s\n",
          index,i,j,
          data[i].m[j].jd,
          data[i].m[j].ipl,
          data[i].m[j].l,
          data[i].m[j].ok ? "" : "(problem)");
  }
}

Freitag, 24. Juli 2015

Ein Exkurs - Fast CGI

Das Projekt swepar, dem dieser Blog gewidmet ist, ist ein hartes Brot, an dem ich noch lange zu kauen haben werde.

"Die Kunst ist lang, das Leben schwierig." Beruflich Software entwickelnd, ist meine Lust oft begrenzt, nach einem anstrengenden Entwicklungstag auch abends noch weiter an Code herumzudrechseln - vor allem wenn es um ein Megaprojekt geht, in das man besser eine Serie voller Arbeitstage investieren sollte.

So bin ich auch immer auf Ausschau nach Alternativen. Der Reiz einer parallelisierten Version der Swiss Ephemeris liegt ja darin, dass sie, auf einem Server angeboten, eine hohe Zahl von Anfragen (scheinbar) gleichzeitig abarbeiten könnte. Dieses Ziel mag aber auch auf andere Weise erreichbar sein.

Der Ephemeridenservice

Auf den Webseiten astrotexte.ch verwende ich derzeit Anfragen in einem einfachen Query-Response-Stil:

Eine Anfrage der Form

ephemeris?planets&houses&jd=2451544.5&lon=7.2833&lat=52.3333
liefert dabei die Antwort
planets:279.8592,217.2933,271.1118,240.9614,327.5755,25.2331,40.4059,314.7841,303.1753,251.4372
houses:192.0288,217.1332,248.6104,285.8994,321.3361,349.8952,12.0288,37.1332,68.6104,105.8994,141.3361,169.8952

Die Syntax des Querystrings wäre natürlich nach Bedarf erweiterbar. In dieser Form benutze ich aber solche Anfragen bereits in einigen Anwendungen.

Das einfache Schema mit einem Querystring, der die gewünschten Daten angibt, und einem Klartextformat (text/plain) als Antwort, strukturiert in Form von Zeilen, denen jeweils ein Bezeichner vorangestellt ist, lässt sich mit wenig Aufwand in jeder beliebigen Programmiersprache weiterverarbeiten.

Natürlich ist ein solcher Aufruf nur als Low-Level-Funktion vorgesehen, um in einer Applikation den Ephemeridenservice aufzurufen.

Eine Webanwendung, in der durch Operationen mit der Maus die Parameter geändert werden können, indem beispielsweise:

  • ein Zeitstrahl als Slider Control verwendet wird, um parallel die entsprechenden Planetenstände im Horoskop anzuzeigen,
  • die Planetensymbole in einem Horoskop mittels Drag und Drop auf eine andere Position gebracht werden können, wobei das Horoskop auf den nächsten Zeitpunkt aktualisiert, zu dem der Planet diese Position erreicht,
  • eine frei definierbare Reihe von abgeleiteten Horoskopen errechnet und angezeigt wird (Auslösungsketten),
erfordert pro Benutzer eine Vielzahl von Ephemerideaufrufen in kurzer Zeit. Hier wäre es wichtig, dass der Server möglichst gute Antwortzeiten für viele Requests liefert.

Alternativen zu Multithreading

Für eine Bibliothek wie SwissEph, die sich so sehr gegen das Parallelisiertwerden sträubt, stehen auch andere Lösungen zur Auswahl.

Wenn der Server ein paar freie Ressourcen hat, ist z.B. FastCGI eine gute Idee und eine interessante Alternative. Es ist eine dieser gut abgehangenen, soliden Lösungen, die uns aus dem letzten Jahrtausend überliefert wurden.

Folgendes Statement am Ende ihrer Dokumentation hat mich gleich für die Entwickler eingenommen:

There is not much development on FastCGI because it is a very stable protocol / application.

Bei klar definierten Infrastrukturentwicklungen geht es eben anders zu als bei Business-Software! Während letztere praktisch nie fertig wird und bis zum Tag der Abschaltung Defects und Änderungswünsche hereinschneien, gibt es für technische Komponenten, sind sie einmal implementiert, praktisch nichts mehr hinzuzufügen - auch die Zahl der sinnvoll einbaubaren Features ist begrenzt.

Die Idee

Die Idee von fastcgi ist einfach erklärt. Statt - wie beim klassischen CGI - für jeden Request einen Prozess zu starten, ein Programm darin auszuführen und dessen Standardausgabe als Antwort an den Client zurückzusenden, steht bei fastcgi eine festgelegte, frei wählbare Anzahl von Prozessen bereit, in einer Warteschleife hängend, um die eingehenden Requests zu beantworten, die reihum an diese Prozesse verteilt werden. Der Start der Prozesse wird dabei auf den Start des gesamten Servers vorverlegt - dann bleiben sie "up and running" bis zum Herunterfahren. So müssen Bibliotheken und Daten nicht pro Prozess neu gelesen werden, sondern stehen in jedem Prozess zur Verfügung.

Dasselbe gilt für Pufferung: Gepufferte Daten bleiben über den Request hinaus erhalten - sogar Daten im Stack können beim nächsten Request weiterverwendet werden (in der Funktion, die die zentrale Warteschleife enthält, oder mit Zeigen auf Daten in darüberliegenden Stackebenen).

Das Programm, das die Anfragen entgegennimmt, ist meist mit einer interpretierten Sprache wie Perl oder PHP geschrieben. Aber auch dies ist ein vermeidbarer Overhead. Es ist selbstverständlich auch ein Programm in Maschinencode möglich, womit man neben Ausführungszeit auch Speicher spart, da kein kompletter Interpreter in den Prozesspeicher geladen werden muss.

So erwuchs in mir der Plan, den obigen Ephemeridenservice als C-Programm zu erstellen.

Installation von fastcgi

Man kann sich die letzte Version von der Seite fastcgi.com herunterladen, ein Link befindet sich unter Application Libraries And Development Kits / Development Kit / Current und trägt den schönen Namen fcgi-2.4.1-SNAP-0311112127. Ich habe unter Ubuntu 14.04 nach Entpacken des heruntergeladenen Verzeichnisses den üblichen Weg probiert, ein Softwarepaket flott zu bekommen: Ich navigierte in die oberste Ebene des Verzeichnisses und führt die Befehlsfolge
./configure
make
make install
aus. Leider scheiterte schon der make-Schritt mit einigen Fehlermeldungen. Nach kurzer Google-Recherche wurde ich auf dem Blog von Muriel Salvan fündig, der vor drei Jahren exakt die gleichen Fehlermeldungen bekamn. Er hatte das analysiert und zwei Korrekturen angegeben, die man nach dem ./configure-Schritt ausführen kann.

Zur Bequemlichkeit für zukünftige Leser habe ich diesen Patch als Datei muriels.patch aufgezeichnet. Wenn man nach dem Download in die oberste Ebene des heruntergeladenen Verzeichnisses wechselt und den Befehl

patch -p1 <../muriels.patch
ausführt (ggf. ../muriels.patch durch den vollen Pfadnamen der Patchdatei ersetzen), so werden die beiden Änderungen automatisch ausgeführt. Danach kann man mit make und make install fortfahren.

Das Gerüst

Nun kann man mit der eigentlich interessanten Entwicklungsarbeit beginnen. Der grobe Programmaufbau ist sehr einfach: zunächst müssen gewisse einmalige Initialisierungsarbeiten vorgenommen werden; danach folgt die endlose Schleife, in der auf eingehende Requests gewartet wird.

Sehr rudimentär sieht das Programmgerüst also folgendermassen aus:

#include "fcgi_stdio.h"
#include ...

int main (int argc, char** argv) {

    initializations( );  

    while (FCGI_Accept() >= 0) {

// Parsen des Requests
      const char *qs = getenv("QUERY_STRING");

// Ausführung der angeforderten Aktion

// Antwort mit printf() ausgeben

     
    } /* while */

  return 0;

}
Statt im QUERY_STRING (wie für meinen Ephemeridenservice) könnte die Daten der Anfrage natürlich auch im HTTP-Body stehen. Dieser wird - wie bei CGI - über stdin zugänglich gemacht.

Das Gerüst des Ephemeridenservices

Der main()-Funktion des Ephemeridenservice habe ich folgende konkrete Gestalt gegeben:
/* ephemeris via fastcgi */

#include "ephemeris.h"

int main (int argc, char** argv) {

    query q = {};
    result r = {};
    
    char ephepath[MAXLEN_EPHE_PATH];
    initializations( ephepath );  

    swe_set_ephe_path( ephepath );

    while (FCGI_Accept() >= 0) {
      init_request(&q,&r);
      TRY {  
        const char *qs = getenv("QUERY_STRING");
        if (qs) {
          parse(qs, &q);
          compute(&q,&r);
          }
        output(&q,&r);
        }
      CATCH {
        do_header( );
        printf("error: %s\n",r.error);
        }  
     
    } /* while */

    return 0;

}
  • In den Initialisierungen (vor Beginn der while-Schleife) lese ich aus einer Konfigurationsdatei .ephemeris den Pfad ephepath des Verzeichnisses ein, in dem die Swiss Ephemeris die Ephemeridendateien vorfinden wird, und teile diesen Pfad der Bibliothek mit.
  • Geht ein Request ein, so initialisiere ich zunächst die Query-Struktur q und die Ergebnisstruktur r, damit nicht etwa noch irgendwelche Relikte aus dem letzten Aufruf übrigbleiben; auch gebe ich Speicher wieder frei, der durch den letzten Request auf dem Heap reserviert worden war.
  • parse(qs,&q): Nun wird der eingehende Querystring analysiert und in das strukturierte Datenobjekt q geschrieben.
  • compute(&q,&r): In dieser Funktion wird die eigentliche Rechnung ausgeführt. Die Ergebnisse landen in der result-Struktur r.
  • output(&q,&r): Hier - und nur hier - werden mit printf() die Ergebnisse in die HTTP-Antwort zurückgeschrieben.

TRY - CATCH in C

Das Sprachfeature try - catch ermöglicht es, Fehler auf einer höheren Stackebene abzufangen, als dort, wo sie aufgetreten sind. Die Sprache C bietet kein try-catch-Idiom an, aber wir können es uns leicht mit den sogenannten "Long Jumps" emulieren. Hinter den Befehlen TRY und CATCH im obigen Listing verbergen sich zwei bescheidene Macros:
jmp_buf _state;
#define TRY if (setjmp(_state)==0) 
#define CATCH else 
#define THROW(errmsg,...) { if (errmsg) sprintf(errmsg,__VA_ARGS__); longjmp(_state, 1); }
Für die Details dieser Konstruktion verweise ich auf einen Blogpost von Francesco Nidito, dessen Idee ich mir hier passend zurechtgekürzt habe: Es funktioniert wunderbar: Mit dem dritten Macro THROW kann ich auf beliebiger Ebene unterhalb von main() alle Stackebenen bis main() löschen und die Programmausführung im CATCH-Block fortsetzen.

Testen

Man kann das Programm auch in der Konsole starten. FCGI erkennt dann, dass es nicht in der Web-Umgebung aufgerufen wird und durchläuft in diesem Fall die Schleife genau einmal. Zum Testen kann man vor dem Programmaufruf die Umgebungsvariable QUERY_STRING setzen.

So kann ich verschiedene erwartete Reaktionen des Programms in einer Datei ./t/test.pl notieren:

#!/usr/bin/perl

use strict;
use warnings;

use Test::More;

ephemeris_call( 
  "Planets for given julian date 2451545",
  "planets&jd=2451545",
  <<EXP );
Content-type: text/plain
planets:280.3689,223.3238,271.8893,241.5658,327.9633,25.2531,40.3957,314.8092,303.1930,251.4548
EXP

ephemeris_call( 
  "Sun And Moon for given julian date 2451545",
  "planets=0,1&jd=2451545",
  <<EXP );
Content-type: text/plain
planets:280.3689,223.3238
EXP

... usw ...

done_testing( );

sub ephemeris_call {
  my ($test,$query_string,$exp) = @_;
  $ENV{'SE_EPHE_PATH'} = "";
  $ENV{'QUERY_STRING'} = $query_string;
  my $act = `ephemeris`;    
  $exp =~ s/\s*$//gm;
  $act =~ s/\s*$//gm;
  is( $act, $exp, $test );
}

Konfiguration im Apache

Für fcgi-Programme habe ich mir in den Webressourcen ein eigenes Verzeichnis /var/www/fcgi eingerichtet: dorthin kam nun das fertige Binary ephemeris.

Um nun das Programm mit dem Apache Server zu "verheiraten", muss es ihm mit der Direktive FastCgiServer bekanntgemacht werden:

<IfModule mod_fastcgi.c>
    
    Alias /ephemeris /var/www/fcgi/ephemeris
    FastCgiServer /var/www/fcgi/ephemeris -processes 5
    
</IfModule>
Nach
  • Einrichtung im Verzeichnis /etc/apache2/mods-available,
  • Stoppen des Servers mit service apache2 stop,
  • Aktivierung der Einstellung mit dem Befehl a2enmod fastcgi
  • Neustart des Servers mit service apache2 start
sieht man mit dem Befehl ps -ef | grep ephemeris die fünf gestarteten Prozesse, wie sie in einer öden Schleife vor sich hindümpeln und, wie morgens auf dem Marktplatz die Tagelöhner, auf Arbeit warten.

Der %lf Bug

Einen Bug musste ich dennoch feststellen. Die Ausgabe von double-Werten mit dem in C dafür vorgesehenen Formatierer %lf im Formatstring von printf funktionierte leider nicht. Die so ausgegebenen Werte wurden rundweg verschluckt.

Eine diesbezügliche Frage in die Entwicklergemeinde bei StackOverflow blieb bislang unbeantwortet.

Vollständiger Quellcode

Hier sind das Programm ephemeris.c und seine Headerdatei ephemeris.h.

Zusammenfassung

Die ganze Technik funktioniert hervorragend und sehr performant. Wenn ich im Zugriffs-Log von Apache die Requestbearbeitungszeiten mit dem Formatierer %D protokolliere, liegen diese praktisch "unter der Nachweisgrenze", auch bei gleichem Input stark schwankend zwischen 100 und 1000 Mikrosekunden pro Request. Da es für 10 Planeten ähnlich lange dauert wie für einen einzelnen (den ich in meiner leicht erweiterten Syntax mit z.B. planets=5 aufrufen kann), vermute ich, dass der Löwenanteil dieser 100 bis 1000 Mikrosekunden in der Apache-Laufzeitumgebung liegt. Aber eine Bearbeitungszeit von einer Millisekunde ist bei Netzwerk-Antwortzeiten von 10-50 Millisekunden sowieso vernachlässigbar.

Mittwoch, 6. August 2014

Refaktorisierung von sweph()

Aufs Geratewohl herausgegriffen aus den verbliebenen knapp 5'000 Zeilen sei die Funktion sweph(). Ich stiess auf sie, weil sie von oben die erste Funktion ist, die die goto-Anweisung verwendet. Sicher ist goto considered harmful (Dijkstra) und sollte daher vermieden werden. Wichtiger ist aber die Verletzung des Single Responsibility Principle (SRP), dass eine Funktion nur für eine einzelne Aufgabe zuständig sein sollte. Das gilt ebenso für alle anderen Funktionen dieses Codes, aber irgendwo muss man einmal anfangen.

Die Funktion sweph() berechnet aus den Files der Swiss Ephemeris die baryzentrische Position des gewünschten Objekts. Gemäss ihrem einleitenden Kommentar hat sie folgende Aufgaben:

/*
 * this function looks for an ephemeris file,
 * opens it, if not yet open,
 * reads constants, if not yet read,
 * computes a planet, if not yet computed
 * attention: asteroids are heliocentric
 *            other planets barycentric
*/
Bei genauerer Inspektion des Codes führt sie die folgenden Aufgaben aus:
  • Sie puffert den letzten Aufruf pro Planet. Die Struktur swed enthält für jeden Planeten einen "Slot" mit den Ergebnissen des letzten Aufrufs.
  • Sie invalidiert den Zeiger auf das Ephemeris File, falls dieses von Datum oder Planetennummer nicht mehr gültig ist. In diesem Fall wird das File auch mit fclose() geschlossen und allokiertes Memory freigegeben.
  • Ist das File ungültig, werden drei mögliche Filenamen für den korrekten Kandidaten ermittelt:
    1. Standardnamen wie sepl_18.se1 oder ast0/se000302.se1,
    2. sogenannte Kurzversionen mit einem "s" vor der File Extension wie in ast0/se00005s.se1,
    3. und schliesslich, falls der Pfad ein Verzeichnis enthielt wie in ast0/se000302.se1, wird ein letzter Versuch für den Filenamen ohne Verzeichnis gemacht, im Beispiel also se000302.se1.
    Die Kandidaten 2 und 3 werden dabei nur für Asteroiden probiert. Es wird versucht, Files mit diesem Namen zu öffnen. Das erste File, für das dies gelingt, wird das neue "aktuelle Ephemeris File".
  • Es werden die für das ganze File gültigen Konstanten in die swed-Struktur eingelesen.
  • Es wird am Anfang und Ende des Gültigkeitsbereichs der Swiss Ephemeris eine Randbedingung fürs Datum geprüft.
  • Nun werden die Koeffizienten des für das Julianische Datum gültigen Interpolationspolynoms eingelesen.
  • Das Polynom wird evaluiert.
  • Für den Sonderfall der baryzentrischen Sonne wird die gesamte Funktion sweph() noch einmal für das Erde-Mond-Baryzentrum aufgerufen und das Ergebnis vektoriell zur Position der heliozentrischen Sonne addiert.
  • Die Asteroidenpositionen sind heliozentrisch und werden in die baryzentrische Positionen konvertiert.
  • Zum Kennzeichen, dass die berechneten Positionen gemerkt werden, wird die Zahl pdp->teval auf das Julianische Datum der Anfrage gesetzt.
Für eine Funktion mit 185 Zeilen sind das ziemlich viele Aufgaben. Auch gehören die Aufgaben zu verschiedenen Ebenen:
  • Ebene File: Datei öffnen und schliessen, Filenamen
  • Abbildung von Filedaten in Datenstruktur: Dateiweite Konstanten und Koeffizienten des zum aktuellen Objekt / Datum gültigen Interpolationspolynoms.
  • Berechnung der baryzentrischen Position: Auswertung des Interpolationspolynoms, Umrechnung von Asteroidenpositionen von helio- in baryzentrisch.
  • Pufferung der Abfrage
Nach Refaktorisierung dokumentiert der Code in der obersten Ebene die einzelnen Aufgaben, aus denen die Funktion besteht:
static int sweph(double tjd, int ipli, int ifno, int iflag, 
                 double *xsunb, AS_BOOL do_save, double *xpret, 
                 char *serr, struct swe_data *swed) {

  int ipl = is_main_object(ipli) ? ipli : SEI_ANYBODY;
  struct plan_data *pedp = &swed->pldat[SEI_EARTH];
  struct plan_data *psdp = &swed->pldat[SEI_SUNBARY];
  struct plan_data *pdp = &swed->pldat[ipl];

  /* In do_save mode, use swe_data memory slot for results */
  double* xp = do_save ? pdp->x : xpret;

  /* if planet has already been computed for this date, return */
  if (sweph_isPositionInBuffer(tjd, ipl, iflag, pdp, xpret))
    return (OK);

  /* Get ephemeris file handle */
  struct file_data* fdp = 
    sweph_getEphemerisFile(tjd, ipli, ifno, pdp, swed, serr);
  if (fdp == NULL )
    return NOT_AVAILABLE;

  if (sweph_isDateOutOfBound(tjd, fdp, serr, swed))
    return NOT_AVAILABLE;

  /* Read coefficient data */
  sweph_getCoefficientSegment(tjd, ipl, ifno, pdp, serr, swed);

  /* Also compute speed, if
   * 1. true position is being computed before applying light-time etc.
   *    this is the position saved in pdp->x.
   *    in this case, speed is needed for light-time correction.
   * 2. the speed flag has been specified.
   */
  bool need_speed = do_save || is_set(iflag, SEFLG_SPEED);

  /* Compute position */
  sweph_evaluateChebyshevPolynomial(tjd, pdp, need_speed, xp);

  /* Special case barycentric Sun */
  if (ipl == SEI_SUNBARY && (pdp->iflg & SEI_FLG_EMBHEL))
    sweph_barycentricSunFromEarth(tjd, ifno, iflag, need_speed, 
                                  pedp, xp, serr, swed);

  /* asteroids are heliocentric.
   * if JPL or SWISSEPH, convert to barycentric */
  if ((ipl >= SEI_ANYBODY) && 
      is_set(iflag, SEFLG_JPLEPH | SEFLG_SWIEPH))
    vec_add(xp, xsunb, need_speed ? 6 : 3);

  /* Mark result in save area as valid buffer */
  if (do_save) {
    sweph_validateBuffer(tjd, ifno, pdp, psdp);
    vec_copy(xpret, xp, 6);
  }

  return (OK);

}
Zwei Ebenen tiefer finden wir die zuvor mit goto again gelöste Aufgabe, das Ephemeridenfile anhand der drei Kandidaten für den Namen zu ermitteln. Die Code-Ersparnis, die das goto motiviert hatte, erwies sich als nicht besonders hoch: Man sparte sich nur die eine Codezeile mit dem Aufruf von swid_fopen(). Dadurch, dass es nun eine eigene kleine Funktkion, zwei Ebenen unterhalb von sweph() gibt, können wir return wie ein break benutzen, um auszusteigen, sobald swid_fopen() einen gültigen Filepointer lieferte:
static FILE* sweph_openEphemerisFile(int ipli, int ifno, char *fname, 
                                     char *serr, struct swe_data *swed) {
 
  FILE *fptr = swid_fopen(ifno, fname, swed->ephepath, serr, swed);
  if (fptr)
    return fptr;
 
  if (!is_main_object(ipli)) {
 
    char *sp;

    /*
     * if it is a numbered asteroid file:
     * try also for short files (..s.se1)
     */
    char s[AS_MAXCH];
    strcpy(s, fname);
    if ((sp = strrchr(s, '.'))) { 
      strcpy(sp, "s." SE_FILE_SUFFIX); /* insert an 's' */
      if ((fptr = swid_fopen(ifno, s, swed->ephepath, serr, swed)))
        return fptr;
      }
 
    /*
     * if we still have 'ast0' etc. in front of the filename,
     * we try in the main ephemeris directory instead of the
     * asteroid subdirectory.
     */
    if ((sp = strrchr(fname, *DIR_GLUE))) {
      if ((fptr = swid_fopen(ifno, sp+1, swed->ephepath, serr, swed)))
        return fptr;
    }
  
  }
  return NULL ;
}
Um festzustellen, wie Code für die Pufferung durchlaufen wird, schrieb ich einen neuen Test test_identical_call() für die Testsuite, der die Funktion se_calc zweimal nacheinander mit identischen Parametern aufruft. Überraschenderweise lieferte die Funktion sweph_isPositionInBuffer immer false!

Grund dafür war, dass exakt der gleiche Code, um festzustellen, ob das Ergebnis sich bereits im Puffer befindet, auch auf höherer Ebene bereits programmiert war - in der darüberliegenden Funktion sweplan(). Ich konnte diesen duplikativen Code nun durch Aufruf der Bedingungsfunktion ersetzen. Insgesamt gab es 6 identische Codestrecken in sweph.c, die jeweils diese Abfrage auf den Puffer machten, die ich nun alle durch den Aufruf der Bedingungsfunktion ersetzen konnte. Solche kleinen Bedingungsfunktionen werden übrigens vom Compiler meist automatisch ge-inlined. Die Bedenken, dass man durch den Funktionsaufruf eine zusätzliche Stackebene aufbauen und somit kostbare Rechenzeit verlieren würde, sind daher unbegründet.

Samstag, 22. Februar 2014

Die 90% Marke ist erreicht

Mit einer Abdeckung von 90.17% habe ich nun den Löwenanteil meiner Version test_calc_reduce.c (aus einer Kopie von sweph.c hervorgegangen) abgedeckt. Dafür werden insgesamt 1194 Testaufrufe der Swiss Ephemeris mit den verschiedensten Parametern ausgeführt, wie das folgende Kommando ergibt:
test_calc_reduce | perl -ne 'BEGIN { $count = 0;} $count += $1 if /(\d+) OK/; END { print qq(Total number of tests: $count\n); }'
Selbstverständlich sind alle API-Funktionsaufrufe enthalten - für die Planetenberechnungen auch mit einer grossen Bandbreite verschiedener Zeiten, verschiedener Flags, verschiedener Planeten. Um auch den Fall abzufangen, dass bei beschädigten Ephemeridenfiles korrekt ein Fehler gemeldet wird, habe ich sogar die Ephemeridendatei sepl_18.se1 in zwei automatischen Tests mit präparierten Versionen ausgetauscht.

Was nun übrig bleibt an nicht abgedecktem Code, betrifft

  • Ausnutzung von Pufferungen: Wenn ein in zeitlicher Nähe liegender Planetenstand zuvor schon einmal berechnet wurde, können Teile der Daten von der letzten Rechnung wiederverwendet werden. Die Zweige, die auf die Pufferung zugreifen, werden in meinen Tests nicht durchlaufen, da ich zwischen den Aufrufen stets die globale Struktur swed initialisiere. Ich könnte zwar in den Testfällen ein Flag clear_buffer vorsehen (das standardmässig gesetzt ist, aber abgeschaltet werden kann), um auch diese Codezweige zu testen. Aber ich will das gar nicht. Pufferungen werde ich nämlich ignorieren, entfernen und bei Bedarf hinterher einfügen, wenn es eine Performanceanalyse als nötig zeigt.
  • Extrem seltene Fehlersituationen: Beispielsweise ein ungültiger Koeffizient mitten im File sepl_18.se1.
  • Grundsätzlich nie durchlaufene Codestellen: Z.B. eine Codestrecke für den Fall baryzentrischer Moshierephemeriden, die aber gar nicht vorgesehen sind und daher weiter oben schon abgefangen werden. Oder eine Fehlermeldung, falls in einer Routine zur Sonnenberechnung nicht wirklich die Sonne
Bei diesem Stand der Testüberdeckung vom 22.2.2014 werde ich es nun im wesentlichen belassen (die Auswertungsfiles wurden mit dem Programm lcov aus den gcov-Resultaten erzeugt). Das Risiko, einige Randwerte beim Refaktorisieren falsch zu behandeln, ist zwar vorhanden, aber es lähmt mich nicht!

Mittwoch, 22. Januar 2014

Initialisierung ausschalten

Das Ausschalten der Initialisierung
    if (0) {
    swe_close( );
    const struct swe_data swed0 = {
      .ephe_path_is_set = false,
      .jpl_file_is_open = false,
      .fixfp = NULL,
      .ephepath =  SE_EPHE_PATH,
      .jplfnam = SE_FNAME_DFT,    /* JPL file name, default */
      .geopos_is_set = false,   /* geopos is set, for topocentric */
      .ayana_is_set = false     /* ayanamsa is set */
      };
    memcpy(&swed,&swed0,sizeof(struct swe_data));
    }
treibt die Abdeckung um ein halbes Prozent in die Höhe. Vermutlich ist es anderer Code, der in diesen Fällen durchlaufen wird. Die Tests bleiben aber alle auf OK, was beruhigend ist: Es gibt anscheinend keine weiteren verschleppten Abhängigkeiten mehr, die die Ergebnisse der Tests beeinflussen könnten.

Sonntag, 12. Januar 2014

82 Prozent

Ich habe nun den Testgenerator erweitert, so dass auch Testdaten für siderische Positionen generiert werden können. Ebenso kann eine topozentrische Position in der Fixture angegeben werden, so dass topozentrische Planetenberechnungen in der Fixture zugänglich sind. Für die Fixsterne und ihre Grössen, für die Ayanamsa-Namen, Planetennamen und die Zeitgleichung erlaube ich mir, die Tests und Erwartungen in test_test_calc_reduce.c hart zu codieren. Sie kommen im Testprotokoll immer am Schluss, nach den Tests aus der Fixture

Und so sieht das Protokoll nun aus.

ruediger@herschel:~/workspace/swepar/src$ test_calc_reduce
geocentric positions: 250 OK   0 NOT OK
helio/barycentric positions: 285 OK   0 NOT OK
barycentric Sun: 5 OK   0 NOT OK
astrometric positions J2000: 90 OK   0 NOT OK
topocentric positions: 15 OK   0 NOT OK
sidereal positions with Lahiri Ayanamsha: 25 OK   0 NOT OK
sidereal positions / Suryasiddhanta, strictly computed: 25 OK   0 NOT OK
individual ayanamsa: 25 OK   0 NOT OK
sidereal rectangular geocentric: 40 OK   0 NOT OK
Asteroid Urania: 5 OK   0 NOT OK
Hamburger Poseidon: 5 OK   0 NOT OK
different coordinate systems: 4 OK   0 NOT OK
Obliquity of ecliptic and nutation: 1 OK   0 NOT OK
Moshier, geocentric positions: 140 OK   0 NOT OK
Moshier, heliocentric positions: 100 OK   0 NOT OK
Provoking a "file not found" error: 1 OK   0 NOT OK
swe_calc_ut: 1 OK   0 NOT OK
New mode: Barycentric, directly from file: 1 OK   0 NOT OK
Fixed Stars: 7 OK   0 NOT OK
Get Ayanamsa: 6 OK   0 NOT OK
Get Planet Names: 3 OK   0 NOT OK
Equation of Time for J2000: 1 OK   0 NOT OK
Die damit erreichte Testabdeckung liegt bei 82 Prozent. In dem arg zerklüfteten gcov-File gibt es jetzt keine Hauptfunktion mehr, die nicht abgedeckt wäre. Die vielen kleinen roten Einsprengsel entstehen:
  • durch zahllose nicht durchlaufene Ausnahmesituationen (z.B. "Ephemeride nach Moshier darf nicht benutzt werden für baryzentrische Koordinaten"). Einiges davon hat mit nicht mehr unterstützten Optionen zu tun.
  • durch nicht ausgenutzte Puffer: Da ich zwischen den Tests normalerweise nun swe_close() aufrufe, wird der Code zum Wiederverwenden bereits berechneter Positionen nicht durchlaufen. Dies könnte ich noch mit einer weiteren kleinen Testserie erreichen.

Freitag, 10. Januar 2014

Der Nageltest für die Local History

Gleich nach meinem gestrigen Lobpreis auf die Local History konnte ich nun die Probe aufs Exempel machen. In der Konsole hatte ich versehentlich ein rm test_calc* ausgeführt, was mir ohne Sicherungen einige Backup-Arbeit beschert hätte. Nun konnte ich auf die Local History zurückgreifen: Sie erschien wieder, nachdem ich die Files test_calc_reduce.c und test_calc_reduce.h wieder leer angelegt hatte.

Nach einigen kleineren Macken konnte ich nun fünf Fixsterntests und vier Ayanamsa-Tests hinzufügen, was die Abdeckung zwar nicht weltbewegend hochgetrieben hat (sie liegt stieg von 68% auf 71% an), aber mit swe_fixstar() und swe_deltat() sind nun zwei prominente Funktionen mit im Rennen.

Die folgende Übersicht zeigt,

  • dass im Bereich der Berechnung von Fixsternpositionen und -grössen noch besonders viel Luft ist.
  • lunar_osc_elem könnte mit ein paar Lilith-Berechnungen besser ausgebeutet werden.
  • sweph, die zentrale Funktion für die Planetenberechnungen, ist erst zu drei Vierteln abgedeckt, obwohl sie mit den bestehenden Tests bereits über 6000 Mal aufgerufen wird.
  • Der Code, der die Anfrage aus dem Puffer bedient, wird nicht durchlaufen, da ich bislang keine zwei identische Positionsanfragen nacheinander stelle.
  • Ausserdem werden keine Asteroiden abgefragt, für die eines der ast-Subdirectories abgefragt werden müsste.
  • Die Abdeckung lässt sich deutlich steigern, wenn noch ein paar topozentrische Berechnungen dazukommen (vor allem durch swid_get_observer).
  • Die Zweige, die Planetenberechnungen mit Ayanamsa enthalten, werden noch nicht durchlaufen.