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.

Donnerstag, 9. Januar 2014

Local History

Hätte ich die Local History von Eclipse schon 2012 gekannt, dann hätte ich nicht folgende Frage im Forum der Clean Code Developer abgesetzt:

Delta-Backup für alle Änderungen an Textdateien

Genau die Umkehrung des üblichen Master/Slave-Verhältnisses leistet die Local History, die für mich ohne weiteres Zutun von Eclipse verwaltet wird. Alles, was ich brauche, kann ich damit machen: Versionen vergleichen (und zwar jeden beliebigen Stand, ohne dass ich Commits einchecken und mit einem Text versehen muss), und bei Bedarf auch alte Versionen wieder zurückholen.

Wenn ich will, kann ich sogar einzelne Deltas aus früheren Versionen wieder holen. Alles sehr komfortabel und nützlich.

Mittwoch, 8. Januar 2014

Über die Benutzung vorgesehener Funktionen. Coverage jetzt 68%

Stand gestern: In der bash ausgeführt, lieferte das Programm test_calc_reduce unsinnige Testergebnisse, während es in Eclipse mit dem Pfad Run as->Local C/C++ Application korrekt ausgeführt wurde. Wie ist dieser Unterschied zu erklären?

Zwischen den Tests hatte ich ja die Struktur swed explizit initialisiert.

    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));
Indem ich nun dieser Initialisierung noch einen Aufruf von swe_close() voranstellte, funktionierte das Programm mit beiden Aufrufmethoden wieder gleich. Es lohnt sich vermutlich nicht, der genauen Ursache nachzugehen, warum es zu diesem Fehlverhalten kam. Denn selbst wenn ich es vergessen haben sollte, ein globales Datum zu löschen: Warum verhielt sich das Programm dann beim Aufruf aus Eclipse anders als beim Aufruf aus der Konsole?

Anderes Thema: Ich habe einmal den Coverage Analyzer gcov angeworfen und festgestellt, dass ich mit der Testabdeckung für sweph.c (also mein test_calc_reduce.c) schon ganz gut liege:

Die vertikale Seitenleiste in der Abdeckungs-Quelltextanzeige zeigt, dass es vor allem im unteren Teil des Files noch viel Rot gibt: da sind einige Funktionen, die noch gar nicht aufgerufen werden. Das topozentrische Koordinatensystem habe ich noch mit keinem Testfall angefordert. Die vielen kleinen Splitter zwischendurch sind Ausnahmefälle wie: Die Fehlermeldung bei der Chironberechnung für ein Datum kleiner als CHIRON_START. Oder ein File, das behauptet, Big Endian zu sein, aber es in Wahrheit nicht ist, weil schon seine eigene Grösse in falscher Codierung abgespeichert ist. Hier werde ich irgendwo einen Stopp machen. Ich könnte natürlich in einem Test ein präpariertes File in den Ephemeridenpfad schieben, das dort das entsprechende korrekte File ersetzt. So könnte ich auch viele dieser Splitter noch "erschlagen".

Eine 100%-Abdeckung bekomme ich voraussichtlich nicht hin. Aber ab 95% könnte ich mit gutem Gefühl damit beginnen, auch aggressivere Umbaumassnahmen in meinem "Branch" test_calc_reduce.c vorzunehmen. Und mir dabei ziemlich sicher zu sein, dass die Kernfunktionen der Swiss Ephemeris erhalten bleiben.

Dienstag, 7. Januar 2014

Statische Variablen

Gestern abend räsonierte ich noch über rätselhafte Fixsterne. Das Rätsel der zwei verschiedenen Positionen desselben Fixsterns (Aldebaran) war schnell gelöst.

Es ist so, wie ich gestern schon schrieb:

Das kann nur bedeuten, dass es noch irgendeinen anderen, nicht in swed enthaltenen, Zustand im Programm gibt. Das werde ich morgen weiter analysieren.
Und dieser verborgene Zustand bestand aus zwei in swe_fixstar deklarierten statischen Variablen:
  static char slast_stardata[AS_MAXCH];
  static char slast_starname[AS_MAXCH];
die zur Pufferung des zuletzt gelesenen Satzes im Fixsternkatalog verwendet wurden.

Das erklärt, dass der innere Zustand des Programms beim zweiten, identischen Aufruf ein anderer war als beim ersten Aufruf. Nun waren die beiden statischen Variablen - noch vom ersten Aufruf - gefüllt.

Wenn es aber nur eine Pufferung war - warum ergab der zweite Aufruf dann andere Werte als der erste?

Der Grund ist, dass ich durch das Löschen der swed (z.B. mittels eines swe_close()-Aufrufs, oder - wie in meinem Testprogramm: explizit) auch die Information swed.is_old_starfile = TRUE zerstört hatte, die beim ersten Aufruf ermittelt worden war (in meinem Ephemeridenverzeichnis befand sich nämlich nur das "alte" File fixstars.cat - was aber ja nicht verboten ist).

  if (swed.fixfp == NULL) {
    if ((swed.fixfp = swi_fopen(SEI_FILE_FIXSTAR, SE_STARFILE, swed.ephepath, serr)) == NULL) {
      swed.is_old_starfile = TRUE;
      if ((swed.fixfp = swi_fopen(SEI_FILE_FIXSTAR, SE_STARFILE_OLD, swed.ephepath, NULL)) == NULL) {
        ...
Wenn der Eintrag im Puffer gefunden wurde, springt das Programm direkt zum Label found. Nach Erreichen der Sprungmarke gibt es aber noch eine Stelle, in der mit den Koordinaten abhängig vom Flag swed.is_old_starfile unterschiedlich weitergerechnet wird. Das erklärt's:
  /* speed in ra and de, degrees per century */
  if (swed.is_old_starfile == TRUE) {
    ra_pm = ra_pm * 15 / 3600.0;
    de_pm = de_pm / 3600.0;
  } else {
    ra_pm = ra_pm / 10.0 / 3600.0;
    de_pm = de_pm / 10.0 / 3600.0;
    parall /= 1000.0;
  }
Nach Ausbau dieser Pufferung und Entfernung der beiden statischen Variablen slast_stardata und slast_starname liefen die beiden Fixsterntests nun beide mit OK durch. Insbesondere müssen die beiden sukzessiven Aufrufe nun auch dieselben Ergebnisse gebracht haben.

Nachdem dies erledigt war, habe ich auch das neuere File sefstars.txt ins Ephemeridenverzeichnis hinzugefügt. Im Grunde ist fixstars.cat ja ein alter Zopf, den ich irgendwann kappen kann.

Wie ich ja früher auch schon den Compileparameter TRACE eliminiert habe.

Das Tool nm, das das Symbolverzeichnis eines Objektfiles auflistet, zeigt, dass es noch zwei weitere statische Variablen in sweph.c gibt: epheflag_sv und nutflag:

ruediger@herschel:~/workspace/swepar/src$ nm -aC sweph.o | grep ' [bB] ' 
00000000 b .bss
00000204 b epheflag_sv.3938
00000200 b nutflag.5350
00000000 b slast_stardata.5374
00000100 b slast_starname.5375
Um diese werde ich mich morgen kümmern.

Verwirrend kommt nun hinzu, dass der Start des Programms test_calc_reduce in der Bash nicht mehr funktioniert, wohl aber der Start von Eclipse aus mittels Run As -> Local C/C++ Application

Montag, 6. Januar 2014

Rätselhafte Fixsterne

Das Programm testgen.c funktioniert wunderbar für die Erzeugung von Test Fixtures, die den Aufruf der Funktion swe_calc betreffen. Nun gibt es aber noch eine ganze Reihe weiterer öffentlicher Funktionen, die für eine volle Testabdeckung des Codes aufzurufen wären. Hierzu hätte man testgen.c durch andere Testtypen erweitern müssen.

Leider war testgen.c nur halbherzig für solche Erweiterungen vorbereitet. Es gibt zwar im Header einer Testreihe ein Feld type. Aber damit hat es sich auch schon. Um es wirklich erweiterbar zu halten, müsste für jeden Testtyp ein Set von drei Funktionen plus eine Datenstruktur bereitgestellt werden:

  • Eine Funktion zum Einlesen der textförmigen Fixture in die Datenstruktur;
  • Eine Funktion, die mit gegebenen Werten in dieser Datenstruktur mittels der Originalbibliothek libswe.so Ausgabewerte erzeugt;
  • Eine Funktion, die mit gegebenen Werten in dieser Datenstruktur die entsprechende Funktion im Prototyp test_calc_reduce.c aufruft und die Ergebnisse mit den durch die vorherige Funktion berechneten (und in der binären Fixture abgespeicherten) Ergebnissen des Standards vergleicht.
So etwas wäre machbar. Aber aus Zeitgründen verschiebe ich das nach hinten und schreibe lieber erstmal einige ad hoc-Tests für die noch verbliebenen Funktionen. Allen voran swe_fixstar().

Hier mache ich nun die merkwürdige Erfahrung, dass die swe_fixstar(), zweimal nacheinander mit denselben Werten aufgerufen, zu unterschiedlichen Ergebnissen führt, obwohl ich die globale Datenstruktur swed zwischen den beiden Aufrufen lösche.

static bool test_fixstar() {

  result r = {};
  const double prec = 1e-7;
  const struct swe_data swed0 = {};

  testdata_fixstar td[] = {
      {
        .star ="Aldebaran",
        .jd = 2451544.5,
        .result_exp = {
           .pos = {
            l:69.7903072,
            b:-5.4676328
          }
        }
      },
      {
        .star ="Aldebaran",
        .jd = 2451544.5,
        .result_exp = {
           .pos = {
            l:69.7903072,
            b:-5.4676328
          }
        }
      }
    };
  const int TESTCOUNT = 2;

  for (int i=0;i<TESTCOUNT;i++) {
    memcpy(&swed,&swed0,sizeof(struct swe_data));
    bool ok = test_single_fixstar( td[i], prec );
    if (ok) r.ok++; else r.not_ok++;
    }

  summary("Fixed Stars",r);

  return r.not_ok == 0;
  }

static bool test_single_fixstar( testdata_fixstar td, double prec) {

  double xx_act[6];
  char serr_act[255], sname[41];

  strcpy(sname,td.star);

  int iflag_act = swe_fixstar( sname, td.jd, td.iflag_in, xx_act, serr_act );

  td.result_exp.pos.r = 1;
  td.result_exp.speed.r = 0;
  if (!td.serr_exp) td.serr_exp = "";

  bool ok;
  ok = diff_check_int( "iflag", iflag_act, td.iflag_exp  );
  ok = diff_check_degree( "pos.l", xx_act[0], td.result_exp.pos.l, prec ) && ok;
  ok = diff_check_degree( "pos.b", xx_act[1], td.result_exp.pos.b, prec ) && ok;
  ok = diff_check( "pos.r", xx_act[2], td.result_exp.pos.r, prec ) && ok;
  ok = diff_check_degree( "speed.l", xx_act[3], td.result_exp.speed.l, prec ) && ok;
  ok = diff_check_degree( "speed.b", xx_act[4], td.result_exp.speed.b, prec ) && ok;
  ok = diff_check( "speed.r", xx_act[5], td.result_exp.speed.r, prec ) && ok;
  ok = diff_check_char( "serr",  serr_act,  td.serr_exp ) && ok;

  return ok;

  }
Das kann nur bedeuten, dass es noch irgendeinen anderen, nicht in swed enthaltenen, Zustand im Programm gibt. Das werde ich morgen weiter analysieren.

Sonntag, 5. Januar 2014

Flags mit GNU bc aufschlüsseln

Mehrmals schon hatte ich beim Testen oder Debuggen das Bedürfnis, die Bestandteile eines iflag-Wertes in seine einzelnen Optionen zu zerlegen, ohne ihn erst (z.B. im Kopf) binär umrechnen zu müssen.

Für solche kleinen Rechnereien empfiehlt sich bc(1), ursprünglich ein kleiner Rechner für die Kommandozeile, und von Entwicklern weiter aufgepeppt. bc(1) kann auch Funktionen, hat Kontrollstrukturen und Schleifen, kennt Arrays, symbolische Variable usw.

Hier mein Tool iflag.

#!/usr/bin/bc -q 

scale = 0

print "iflag = "
iflag = read()

define bit( i, p ) {
  return ((i/p) % 2 == 1)
  }

if (bit( iflag, 1)) print "seflg_jpleph\n";
if (bit( iflag, 2)) print "seflg_swieph\n";
if (bit( iflag, 4)) print "seflg_moseph\n" ;
if (bit( iflag, 8)) print "seflg_helctr\n";
if (bit( iflag, 16)) print "seflg_truepos\n";
if (bit( iflag, 32)) print "seflg_j2000\n";
if (bit( iflag, 64)) print "seflg_nonut\n";
if (bit( iflag, 128)) print "seflg_speed3\n";
if (bit( iflag, 256)) print "seflg_speed\n";
if (bit( iflag, 512)) print "seflg_nogdefl\n";
if (bit( iflag, 1024)) print "seflg_noaberr\n";
if (bit( iflag, (2*1024))) print "seflg_equatorial\n";
if (bit( iflag, (4*1024))) print "seflg_xyz\n";
if (bit( iflag, (8*1024))) print "seflg_radians\n";
if (bit( iflag, (16*1024))) print "seflg_baryctr\n";
if (bit( iflag, (32*1024))) print "seflg_topoctr\n";
if (bit( iflag, (64*1024))) print "seflg_sidereal\n";
if (bit( iflag, (128*1024))) print "seflg_icrs\n";
if (bit( iflag, (256*1024))) print "seflg_new_mode\n";

quit
Und hier eine Beispielsitzung:
ruediger@herschel:~/workspace/swepar/src$ iflag 
iflag = 16912
seflg_truepos
seflg_nogdefl
seflg_baryctr

Nachtrag

Gegen eine Scriptsprache, die sich gewaschen hat, wie z.B. Perl, kann bc natürlich nicht ankommen. Das ist auch gar nicht die Absicht. bc will nicht viel mehr als eine Art fortgeschrittener Taschenrechner sein, der in die Unix-Philosophie integriert ist, d.h. sich komfortabel mit anderen Befehlen in der Kommandozeile kombinieren lässt.

Hier das selbe Beispiel in Perl:

#!/usr/bin/perl

my @SEFLG=qw(SEFLG_JPLEPH  SEFLG_SWIEPH    SEFLG_MOSEPH 
             SEFLG_HELCTR  SEFLG_TRUEPOS   SEFLG_J2000   
             SEFLG_NONUT   SEFLG_SPEED3    SEFLG_SPEED   
             SEFLG_NOGDEFL SEFLG_NOABERR   SEFLG_EQUATORIAL 
             SEFLG_XYZ     SEFLG_RADIANS   SEFLG_BARYCTR 
             SEFLG_TOPOCTR SEFLG_SIDEREAL  SEFLG_ICRS   
             SEFLG_NEW_MODE);
my $iflag = $ARGV[0];
foreach my $flag (@SEFLG) {
  print "$flag\n" if $iflag % 2;
  last unless $iflag >>= 1;
  }

Samstag, 4. Januar 2014

Mythos Top-Down-Schleife

In einem früheren Blog mutmasste ich, dass ich in test_calc_reduce.c die Macros dot_prod und square_sum bereits eliminiert hätte. Diese Mutmassung war falsch. Es passiert mir oft, dass ich in der Erinnerung das intensiv ins Auge gefasste Vorhaben, etwas zu tun, damit verwechsle, ich hätte es bereits getan - wohl ein allgemein-menschliches Problem, das etwas mit der Art zu tun hat, wie unser Gedächtnis funktioniert.[1]

Also führte ich es nun aus: Ich führte zwei neue Inline-Funktionen in test_calc_reduce.h ein:

static inline double vec_dot( void* v, void* w, int n) {
  double *fv = (double *)v,
         *fw = (double *)w,
         s = 0;
  if (v && w) for (int i=n-1;i>=0;--i) s += fv[i]*fw[i];
  return s;
  }

static inline double vec_length( void* v, int n) {
  return v ? sqrt( vec_dot(v,v,n) ) : 0;
  }
Dann ersetzte ich die Aufrufe der Macros dot_prod und square_sum durch die entsprechenden Inline-Funktionen. Während nach der dot_prod-Ersetzung die Tests noch sauber durchliefen, kam es nach der square_sum-Ersetzung (alles schön in einzelnen Versionen der Dropbox nachvollziehbar!) zu einer Verletzung eines Tests. Ein einziger Test schlug fehl - die Berechnung einer Geschwindigkeitskomponente des vermeintlich "wahren" Mondknotens (Swiss Ephemeris Objekt 13) war in einem einzigen Fall zu ungenau:
ruediger@herschel:~/workspace/swepar/src$ test_calc_reduce tcr.bin

          xx[3] :    2.5273368038 <>    2.5273404552
Differences from test case 138
Juldat 2263323.5683292975, planet 13, iflags 258

Swiss Eph, geocentric positions: 249 OK   1 NOT OK
Natürlich prüfte ich zuerst sehr genau, ob ich irgendwelche Fehler bei den Ersetzungen gemacht hatte. Sie waren aber alle korrekt. Als nächstes debuggte ich den Testfall und verglich die Ergebnisse von square_sum und vec_dot bei jedem Aufruf. Dabei kam es zu kleinen Abweichungen. Für die Werte
xpos[2][3]  double  -2.4692434310721423e-05
xpos[2][4]  double  -0.00055779853107697566
xpos[2][5]  double  -4.9608092113099317e-05
ergab die Berechnung
v2 = vec_dot(xpos[i]+3,xpos[i]+3,3);  // 3.1420988038692279e-07
v2 = square_sum((xpos[i]+3));         // 3.1420988038692284e-07 
// Exakter Wert mit Wolfram Alpha:    // 3.14209880386922820362148922449447018e-07
Ärgerlicherweise liefert das Macro sogar den genaueren Wert, obwohl man doch annehmen sollte, dass nach Auflösung der Inline-Funktion praktisch dieselben Maschineninstruktionen entstehen! Die exakte Berechnung führte ich mit den obigen Werten für xpos[2][3..5] mit der Suchmaschine Wolfram Alpha durch.

Wie kommt es zu diesem Unterschied?

Die Erklärung ist banal: Das Rechnen mit Gleitkommazahlen verletzt das Assoziativgesetz aus Rundungsgründen. Es kommt also darauf an, in welcher Reihenfolge man eine Summe mit drei Summanden berechnet. Dagegen brachte Wolfram Alpha, wie zu erwarten, immer dasselbe Ergebnis, unabhängig von der Reihenfolge der Summanden, da sie solche Berechnungen mit arbitrary precision ausführt.

Nun hatte ich aber das Skalarprodukt mit einer Top-Down-Schleife berechnet, wobei ich einem Performance-Mythos aufgesessen war, der möglicherweise in den 80er Jahren noch seine Berechtigung gehabt hatte.[2]

Nach Umstellung auf eine normale Bottom-Up-Schleife,

static inline double vec_dot( void* v, void* w, int n) {
  double *fv = (double *)v,
         *fw = (double *)w,
         s = 0;
  if (v && w) for (int i=0;i<n;i++) s += fv[i]*fw[i];
  return s;
  }
die auch für das Auge angenehmer = lesbarer ausfällt, war wieder alles OK, die Inline-Funktion ergab bis auf die letzte Dezimale denselben Wert wie das Macro.


[1] Das ist wohl ein allgemein-menschliches Problem, Ich habe den Eindruck, bei der Erinnerung spielen die Wertungen "wichtig"/"unwichtig" bzw. "interessant"/"uninteressant" eine Rolle, mit der Tendenz, dass Unwichtiges oder Uninteressantes aus dem Gedächtnis herausfällt. Wenn man sich etwas vornimmt, heisst das, dass man einem Thema Aufmerksamkeit schenkt und sich die vorzunehmende Handlung imaginiert. Dadurch kann dieser Vorsatz in der Erinnerung mit der Erinnerung an die tatsächlich vorgenommene Handlung verwechselt werden.
[2] Erst recht ist natürlich das Dekrementieren des Index mit dem Präfix-Operator statt mit der üblichen Postfix-Notation nicht durch Performancegründe zu rechtfertigen. Auch dies beleidigt nur das Auge - weiter nichts!

Freitag, 3. Januar 2014

Heliozentrisch!

Der Grund dafür, dass meine experimentellen Funktion calcMainPlanetFromFile() für die inneren Planeten fehlschlägt, scheint der zu sein, dass diese im File helio- und nicht baryzentrisch abgespeichert sind. Der Debugger zeigt, dass bei einem Vergleichsaufruf von swetest mit denselben Daten in Zeile 1601ff. von sweph.c in diesem Fall noch die baryzentrische Sonnenposition zur errechneten (heliozentrischen) Position addiert wird:
      /* if planet is heliocentric, it must be transformed to barycentric */
      if (pdp->iflg & SEI_FLG_HELIO) {
 /* now barycentric planet */
 for (i = 0; i <= 2; i++) 
   xp[i] += xps[i];
 if (do_save || (iflag & SEFLG_SPEED))
   for (i = 3; i <= 5; i++) 
     xp[i] += xps[i];
      }
Im Debugger sieht man, dass vor dieser Korrektur die errechneten Werte noch exakt mit denen aus meiner experimentellen Funktion übereinstimmen (im Konsolenfenster im Hintergrund des Screenshots sieht man die Soll/Ist-Abweichungen des Tests).

Die Flagleiste pdp->iflg gehört dabei zur Struktur plan_data, die beim Lesen des Files gefüllt wird. Die Information, ob die Daten heliozentrisch oder baryzentrisch sind, ist also nicht statisch im Programm definiert, sondern mit den Daten zusammen im File abgespeichert.

Habe mich auch wieder daran erinnert, dass gdb-Breakpoints in dynamischen Bibliotheken erst dann gesetzt werden können, wenn die Bibliothek innerhalb der Debuggingsitzung das erste Mal geladen wurde.

pdfposter

Da der Aufrufgraph von sweph.c (siehe letzter Post) mir im Druck immer auf A4 reduziert wurde, so dass selbst ein weniger alterssichtiger Mensch als ich die Texte der einzelnen Knoten nicht mehr lesen kann, suchte ich ein Programm für den Posterdruck und fand Hartmut Goebels pdfposter.

Mit diesem Programm gelang mir die Aufteilung des PDF-Dokuments auf vier A4-Seiten, die zusammen ein A2-Poster im Hochformat ergeben. Beim Ausdruck ergaben sich allerdings Ränder, auf denen die Farbe "schmierte". Ob es am Treiber oder am pdfposter-Programm selbst liegt, war mir zu mühsam herauszufinden. Ich rückte den vier Seiten einfach mit der Schere zu Leibe und habe nun mein sweph.c-Poster.

Aufrufgraph auf vier A4-Seiten

Donnerstag, 2. Januar 2014

Ein Aufrufgraph für sweph.c

Das Kernstück der Swiss Ephemeris ist das Includefile sweph.c. Um die Kohäsion und die Kopplung der einzelnen Funktionen einzuschätzen und das Programm insgesamt besser zu verstehen, hatte ich das Bedürfnis nach einem Aufrufgraphen. Es soll für ein bestimmtes C-Sourcefile (natürlich sweph.c, das fetteste) ein Diagramm erstellt werden, das für jede Funktion zeigt, welche anderen Funktionen sie jeweils aufruft.

Der in Eclipse CDT vorgesehene View Call Graph klingt vom Namen her verheissungsvoll, bricht aber bei Aufruf mit folgendem Stacktrace ab:

org.eclipse.core.runtime.AssertionFailedException: null argument:Action must not be null
 at org.eclipse.core.runtime.Assert.isNotNull(Assert.java:85)
 at org.eclipse.jface.action.ContributionManager.add(ContributionManager.java:75)
 at org.eclipse.linuxtools.internal.callgraph.CallgraphView.createPartControl(CallgraphView.java:532)
 at org.eclipse.ui.internal.e4.compatibility.CompatibilityPart.createPartControl(CompatibilityPart.java:138)
 at org.eclipse.ui.internal.e4.compatibility.CompatibilityView.createPartControl(CompatibilityView.java:155)
 at org.eclipse.ui.internal.e4.compatibility.CompatibilityPart.create(CompatibilityPart.java:313)
        ...
Der Hauptgrund für den Abbruch sind nicht etwa die in meinem System fehlenden Komponenten wie System Tap (ein interessantes Programm, wie mir scheint), sondern dass der View Call Graph nicht statisch, sondern dynamisch gebildet wird. Dieser Call Graph ist ein Werkzeug fürs Profiling, zum Studium des Laufzeitverhaltens – und nicht das statische Analysetool, nach dem ich suche. Also suchte ich nach einer Alternative.

Das Projekt besteht aus einer Reihe von C-Dateien, die wie üblich jede einzeln in eine gleichnamige Objektdatei compiliert wird. Um einen Aufrufgraphen zu erstellen, müssen wir wissen,

  • welche Funktion welche anderen Funktionen aufruft, und
  • welche Funktionen in welchem Include enthalten sind, um die Funktionen als Gruppe markieren zu können.

Listen der in einem Include enthaltenen Funktionen kann man aus dem jeweils generierten Objectfile (*.o) herauslesen, das ja eine Symboltabelle der compilierten Funktionen enthält.

Funktionen, die mittels Compilerschaltern ausgeblendet sind, werden bei diesem Vorgehen allerdings ignoriert, aber das stört mich hier nicht besonders. Falls das in anderen Situationen stören sollte, kann man auch CTags verwenden, um den Quellcode selbst auf Funktionsdefinitionen zu durchforsten.

Mit dem GCC-Tool objdump kann diese Symboltabelle mit der Option -t ausgelesen werden. Für sweph.o beginnt der Symboltabellen-Dump z.B. wie folgt:

ruediger@herschel:~/workspace/swepar/src$ objdump -t sweph.o

sweph.o:     file format elf32-i386

SYMBOL TABLE:
00000000 l    df *ABS* 00000000 sweph.c
00000000 l    d  .text 00000000 .text
00000000 l    d  .data 00000000 .data
00000000 l    d  .bss 00000000 .bss
00000000 l    d  .rodata 00000000 .rodata
00000000 l     O .rodata 000000a8 pla_diam
000000c0 l     O .rodata 000001e0 ayanamsa
00000000 l    d  .data.rel.local 00000000 .data.rel.local
00000000 l     O .data.rel.local 0000006c ayanamsa_name
00000420 l     O .rodata 0000002c pnoint2jpl
00000460 l     O .rodata 00000054 pnoext2int
00000204 l     O .bss 00000004 epheflag_sv.3938
000005b1 l     F .text 00001a56 swecalc
0000f554 l     F .text 00000165 denormalize_positions
0000f6b9 l     F .text 000000d6 calc_speed
0000faa0 l     F .text 000000c4 plaus_iflag
0000395c l     F .text 000004e1 jplplan
00003169 l     F .text 000007f3 sweplan
In der Symboltabelle wird so allerhand Mystisches aufgeführt; aber eben auch alle enthaltenen Funktionen; diese sind dadurch zu erkennen, dass sie in der Liste das Kürzel F haben und ihr Name mit .text eingeleitet wird.

Mit der folgenden Perl-Routine kann ich alle Object Files eines Verzeichnisses auf einmal durchgehen und eine Hashtabelle erstellen, die für jede Funktion die Nummer des Object Files liefert, in der sie vorhanden ist:

sub get_all_symbols {

# Generates a hash containing all symbols 
# of all object files in this directory

  my %symbols;
  my @obj = split /\s+/, `ls *.o`;
  my $count = -1;
  my $exceptions = join "|", @ignore_object_files;
  for my $fname (@obj) {
    next if $fname =~ /$exceptions/; # Skip the exceptions
    $count++;
    foreach (`objdump -t $fname`) {
      chomp;
# This is very GCC specific, expecting the OBJDUMP output format
      next unless /^[0-9a-f]{8}.*?F\s+\.text\s+[0-9a-f]{8}/;
      my $name = $_;
      $name =~ s/^.*?(\w+)$/$1/;
      $symbols{$name} = $count;
      }
    }
  return \%symbols;
  }
Hierbei kann man mit dem Array @ignore_object_files noch eine Reihe von Object Files aufführen, die nicht analysiert werden sollen.

Für die eigentliche Aufrufanalyse gibt es das Programm cflow. Mit der Option -l aufgerufen liefert es die gewünschte Vorwärtsreferenztabelle: Für jede Funktion werden alle von dieser Funktion aufgerufenen Funktionen gelistet. Für sweph.c beginnt die Liste beispielsweise so:

ruediger@herschel:~/workspace/swepar/src$ cflow -l sweph.c | head
{   0} swe_calc() :
{   1}     memset()
{   1}     fopen()
{   1}     fgets()
{   1}     strchr()
{   1}     atol()
{   1}     fclose()
{   1}     swi_open_trace()
{   1}     trace_swe_calc() <void trace_swe_calc (int swtch, double tjd, int ipl, int32 iflag, double *xx, char *serr) at sweph.c:5976>:
{   2}         fputs()
Wieder ist es so, dass bei einem automatischen Auswerten dieser Information bestimmte Aufrufe ignoriert werden sollten, hier z.B. die Aufrufe von Standardfunktionen wie fclose() usw.

Im Grunde sind damit alle Informationen beisammen, um den Graphen zu erstellen, was am besten auf dem Weg über die graphische Sprache dot von graphviz geschieht. Um die GraphViz-Befehle zu erforschen, gibt es übrigens eine schöne Webanwendung, den GraphViz Workspace. Er erzeugt bereits während der Eingabe das resultierende Diagramm.

Ich orientierte mich, was die dot-Syntax und Gestalt der zu erzeugenden Knoten angeht, an dem Script cflow2dot.pl.

Meine Version unterscheidet sich von der Vorlage im wesentlichen nur dadurch, dass es Filtermöglichkeiten für die Objektfiles und die Listen der aufgerufenen Funktionen gibt, und dass die verschiedenen Farbcodes dafür genutzt werden, die Funktionen nach ihrer Zugehörigkeit zu Source-Code-Includes zu gruppieren. Die Formen dagegen entsprechen der jeweiligen Aufruftiefe (so ist es auch im Originalskript).

Hier ist das Ergebnis für sweph.c (noch mit inkscape von SVG ins PDF-Format convertiert):

Die Trace-Funktionen sind gelb markiert, da sie in keinem Objektfile gefunden werden konnten (denn ich kompiliere das Projekt nicht mit dem Traceschalter). Auch konnten die Macro-Funktionen dot_prod und square_sum nicht gefunden werden. Aus Gründen der Compilersicherheit wäre es sowieso besser, diese Funktionen als Inline-Funktionen zu deklarieren (was ich in meinem Prototyp test_calc_reduce.c auch schon umgestellt habe - glaube ich).

Donnerstag, 26. Dezember 2013

Backup und kurzfristige Versionierung mit der Dropbox

Da ich dieses Projekt nicht im Team bearbeite, sondern solo, ist ein "richtiges" Repository ein bisschen mit Kanonen auf Spatzen schiessen. Eigentlich reicht vorerst die auf dreissig Tage befristete Versionierung, wie sie die Dropbox-Gratis-Accounts anbieten.

In Linux gibt es die schönen symbolischen Links. Ich habe einen Ordner angelegt, der symbolische Links auf die Dateien enthält, die von der SwissEphemeris-Distribution abweichen:

#!/usr/bin/perl -w

use strict;

my @syncfiles = qw(
  test_calc_reduce.c
  test_calc_reduce.h
  test_test_calc_reduce.c
  test_test_calc_reduce.h
  testgen.c
  testgen.h
  tcr.fix
  tests.cpp
  tests.hpp
  makefile
  );

foreach my $file (@syncfiles) {
  `ln -s ~/workspace/swepar/src/$file $file`;
  }
... und hier ist der

Link auf den Dropbox-Ordner

Der Mondknoten

Ich habe nun die Änderungen der letzten beiden Releases, soweit sie das Include sweph.c betrafen, in meine experimentelle Version test_calc_reduce hineingearbeitet.

Ein diff -uw mit der neuen Version förderte auch die folgende Modifikation an sweph.c zutage, die ich damals eigenmächtig vorgenommen hatte - und die nun die mysteriösen Differenzen produziert hatte!

@@ -783,8 +789,6 @@
     if (!(iflag & SEFLG_SIDEREAL) && !(iflag & SEFLG_J2000)) {
       ndp->xreturn[1] = 0.0; /* ecl. latitude       */
       ndp->xreturn[4] = 0.0; /*               speed */
-/* For true node, radial speed can be ommitted too. RPL, 17.11.2011 */
-      ndp->xreturn[5] = 0.0;  /* radial        speed */
       ndp->xreturn[8] = 0.0; /* z coordinate        */
       ndp->xreturn[11] = 0.0; /*               speed */
     }
Damit hatte ich eine Änderung am Referenzalgorithmus (sweph.c) vorgenommen. Eine genauere Prüfung zeigt, dass dies meine einzige Änderung am Referenzcode ist.

Die Fixture neu generiert mit

testgen tcr.fix tcr.bin
- und schon kommt das neue Output wie erwartet:
ruediger@herschel:~/workspace/swepar/src$ test_calc_reduce tcr.bin | grep OK
Swiss Eph, geocentric positions: 250 OK
Swiss Eph, helio/barycentric positions: 285 OK
Swiss Eph, different coordinate systems: 4 OK   
Obliquity of ecliptic and nutation: 1 OK   
Moshier, geocentric positions: 140 OK   
Moshier, heliocentric positions: 100 OK   
Provoking a "file not found" error: 1 OK   
New mode: Barycentric, directly from file: 1 OK   

Montag, 23. Dezember 2013

Präzessionsalgorithmen geändert

Die Erklärung der kleinen Abweichungen meiner Tests dürfte darin liegen, dass ich die neueste Version der Swiss Ephemeris für die Berechnung der Testerwartungen verwende. Die Release Notes sagen aber:
1.78  2-aug-2012  Precession calculation updated acc. to Vondrák et alii 2012
Die Vondrák-Codeänderung in sweph.c ist in meinem vorläufigen Ersatzinclude test_calc_reduce.c für sweph.c noch nicht vorhanden. Es muss ge-merged werden.

LD_LIBRARY_PATH ist anders - in Ubuntu

Es gefiel den Entwicklern von Ubuntu, Änderungen der Variable LD_LIBRARY_PATH in der .profile-Datei nicht zuzulassen. Was für beliebige andere Variablen geht – z.B. auch für PATH – wird in Ubuntu für LD_LIBRARY_PATH stillschweigend ignoriert.

Hier die Erklärung - nebst zwei Workarounds, wie man es trotzdem hinbekommt (ich habe mich für den zweiten entschieden).