Mittwoch, 21. Dezember 2011

Noch ein Zugang

Als Zwölfjähriger bekam ich zu Weihnachten ein Fahrrad geschenkt. Es war mir ein Vergnügen, dieses Fahrrad einmal vollständig bis zur letzten Schraube auseinanderzunehmen und danach wieder zusammenzubauen. Es blieben beunruhigenderweise am Schluss einige Schrauben übrig, aber der Fahrtüchtigkeit meines "Highrisers" tat das keinen Abbruch.

Diese Anlage - alles auseinanderzunehmen und wieder zusammenzubauen - könnte ich auch beim swepar-Projekt einbringen. Das ergibt - neben den vielen kleinen Refaktorisierungen, mit denen ich schon begonnen habe, einen weiteren Zugang zu diesem Projekt. Ich könnte damit beginnen, die Berechnung der baryzentrischen Koordinaten der Planeten aus der Swiss Ephemeris zu reproduzieren. Danach könnte das ganze Fleisch hinzugefügt werden: Umrechnung heliozentrisch -> geozentrisch, Berücksichtigung der Präzession, schliesslich Moshier als Ausweich-Ephemeride...

Die folgende kleine Funktion enthält alle nötigen Aufrufe, um die Position eines einzigen Planeten aus der Swiss Ephemeris zu berechnen. Das sieht doch eigentlich überschaubar aus...
 int test_direct() {

double tjd = 2452545.0;
int ipl = 5;
char qfname[255],fname[255],serr[255];
double xp[6];

// File name
swi_gen_filename(tjd, ipl, fname);
strcpy(qfname, "\\sweph\\ephe\\");
strcat(qfname, fname);

// Open the file
FILE* fp = fopen(qfname, BFILE_R_ACCESS);

// Put the minimum data necessary into swed
swed.fidat[0].fptr = fp;
strcpy(swed.fidat[0].fnam, fname);

// Read the necessary stuff
read_const(0, serr, &swed);

struct plan_data *pdp = & swed.pldat[ipl];
get_new_segment(tjd, ipl, 0, serr, &swed);
if (pdp->iflg & SEI_FLG_ROTATE)
rot_back(ipl,&swed);
else
pdp->neval = pdp->ncoe;

/* evaluate chebyshew polynomial for tjd */
double t = (tjd - pdp->tseg0) / pdp->dseg;
t = 2 * t - 1;
for (int i = 0; i <= 2; i++) {
xp[i] = swi_echeb (t, pdp->segp+(i*pdp->ncoe), pdp->neval);
xp[i+3] = swi_edcheb(t, pdp->segp+(i*pdp->ncoe), pdp->neval)
/ pdp->dseg * 2;
}

for (int i=0;i<6;i++) printf("xp[%d] = %15.10f\n",i,xp[i]);
return OK;
}

Montag, 28. November 2011

Long Jumps eingeführt

Wie man in C++ oder Java für verschiedene Ausnahmesituationen eigene Ausnahmeklassen definiert, kann man in C für verschiedene fatale Ausnahmesituationen eigene Ausnahmefunktionen definieren. Dabei bezeichne ich eine Ausnahmesituation als fatal, wenn die Funktion in keiner sinnvollen Weise mehr fortgesetzt werden kann, sondern sofort abgebrochen werden muss, mit Ausnahme einer allfälligen Bereinigung von Ressourcen (free, fclose) und Rückgabeparametern.

Mit der C-Technik der "langen Sprünge" - Long Jumps - ist es möglich, zu einer weiter oben im aktuellen Aufrufstack definierten Stelle zurückzuspringen und dabei alle darunterliegenden Stackebenen aufzulösen. Es handelt sich um eine altbewährte Technik, die sicher bereits in den ersten Versionen von C vorhanden war und auch im Buch von Kernighan und Ritchie dokumentiert ist (Abschnitt B8 - nichtlokale Sprünge).

Das erlaubt es, grundlegende Ausnahmen "oben" zu behandeln, ohne sich in den detaillierteren Stackebenen darum kümmern zu müssen. Es ist also nicht nötig, bei einem Stack der Tiefe N auf jeder Ebene einenen Rückgabewert abfragen zu müssen:
int retc = deep_function(t,i,&x,&y,&z);
if ( (retc == ERR_FILE_NOT_FOUND) ||
(retc == ERR_FILE_DAMAGED) ||
(retc == ERR_OUT_OF_RANGE) ) return retc;

Stattdessen schreibt man nur, was man tun möchte:
deep_function(t,i,&x,&y,&z);

Die Behandlung solcher schwerwiegenden Fehlersymptome findet an einer zentralen Stelle weit oben im Stack statt, zu der man mittels eines Long Jump direkt von der fehlerverursachenden Stelle springt.

Eine solche zentrale Stelle wäre hier die neue Funktion se_calc(). Die Stelle in se_calc(), die dem try-Block entspricht - also die Folge von Anweisungen, bei deren Ausführung potentiell ein Long Jump ausgelöst wird - muss in C in einem if-Block stehen:

int se_calc( 
double tjd,
int ipl,
int iflag,
double *xx,
char *serr,
struct swe_data *swed) {

int iflag_act = iflag,
ephemeris_requested = iflag & SEFLG_EPHMASK,
ipl_act = ipl,
retc;

if ((retc = setjmp((swed->state).env))==0) {

se_calc_prepare( tjd, &ipl_act, &iflag_act, xx, serr, swed);
if (iflag_act==ERR) return ERR;

double xint[24]; // Default result area for ephemeris computation
double *xp = xint;

int iflag_ret = swecalc(tjd, ipl_act, iflag_act, & xp, serr, swed);
return iflag_ret >= 0 ?
se_calc_map_results(
ipl_act,
iflag_ret,
ephemeris_requested,
xp,
xx )
: ERR;

}
else {
return retc;
}
}

Der Kontext, der mit setjmp erzeugt wird, muss im Aufrufstack weitergereicht werden (hier im neuen Member swed->state). Wenn die Ausnahme mit longjmp ausgelöst wird, muss dieser Kontext der Funktion longjmp als Argument mitgegeben werden. Es werden dann alle Stackebenen bis zu der Stelle abgebaut, in der setjmp aufgerufen wurde, und die Funktion kehrt mit dem als zweites Argument angegebenen Rückgabewert (der natürlich ungleich Null sein muss) zurück.

Hier ein Beispiel einer Ausnahme, die über drei Funktionen in steigender Allgemeinheit schliesslich den longjmp in die oberste Stackstufe ausführt:

static void throw_file_damaged(
struct file_data *fdp,
char* serr,
struct swe_data *swed) {

char *serr_file_damage = "Ephemeris file %s is damaged.";
char msg[AS_MAXCH];
int errmsglen = strlen(serr_file_damage) + strlen(fdp->fnam);
sprintf(msg, serr_file_damage,
errmsglen < AS_MAXCH ? fdp->fnam : "" );
throw_file_error(fdp->fptr,msg,serr,swed);
}

static void throw_file_error(
FILE *fp,
char* msg,
char* serr,
struct swe_data *swed) {

fclose(fp);
throw(ERR,msg,serr,swed);

}

static inline void throw(
int retc,
char* msg,
char *serr,
struct swe_data* swed) {
if (serr != msg) msg_append( serr, msg );
longjmp( swed->env, retc);
}


Natürlich betrifft diese ganze Springerei nur zentrale Ausnahmen, die die weitere Bearbeitung direkt verhindern. Auch spricht natürlich nichts grundsätzlich gegen Rückgabewerte von Funktionen - solange diese Rückgabewerte nur direkt nach dem Aufruf verarbeitet werden und nicht durch N Stackebenen abgefragt und weitergereicht werden müssen.

Um den Zauber der Long Jumps zu ermöglichen, benötigt man ein implementierungsabhängiges Datenhäppchen vom Typ jmp_buf - einen Array von einer nicht festgelegten Anzahl von Bytes. Es empfiehlt sich, diesen Array in eine Struktur jmp_state zu packen, um ihn bei Bedarf kopieren zu können (denn Arrays als Members werden ja durch eine Zuweisung komplett kopiert, ohne dass deren Grösse explizit verwendet werden muss):
typedef struct {
jmp_buf env;
} jmp_state;

Wann sollte es nötig sein, diese Rücksprunginformation zu kopieren? In den seltenen Fällen einer "ausnahmsweise abweichenden Ausnahmebehandlung". Die folgende Funktion zum Beispiel benutze ich, solange im aktuellen Mischcode (Altes und Neues) auch noch fatale Ausnahmen per Returncode behandelt werden sollen:
static int sweph_nothrow(
double tjd,
int ipli,
int ifno,
int iflag,
double *xsunb,
AS_BOOL do_save,
double *xpret,
char *serr,
struct swe_data *swed) {
jmp_state state_save = swed->state;
int retc;
if ((retc = setjmp((swed->state).env)) == 0) {
sweph(tjd,ipli,ifno,iflag,xsunb,do_save,xpret,serr,swed);
}
swed->state = state_save;
return retc;
}
Diese Funktion soll noch auf die alte Art arbeiten, also einen Returncode zurückgeben. Innere Funktionen sind aber bereits auf die Long Jump-Technik umgestellt. Was tun?

Ich merke mir die Information über den globalen Rücksprung in der lokalen Variablen state_save. Durch einen neuen setjmp-Aufruf wird dann ein neues Zustands-Datenhäppchen erzeugt, das im if-Zweig der setjmp-Anweisung gültig ist. Danach kehre ich wieder zum alten Zustand zurück. Gibt es keine Ausnahme, ist retc = 0. Andernfalls haben wir den bei Auslösung des Sprungs gesetzten Returncode. In allen Fällen wird dieser retc mit dem erwarteten Wert zurückgegeben.

Freitag, 25. November 2011

github-Repository eingerichtet

Habe nun für das Projekt ein github Repository eingerichtet (mein letztes freies Repository bei github, denn man hat nur fünf kostenlose).

https://github.com/rplantiko/swepar/

Grund dafür war vor allem die Sorge um Datenverlust: Da ich mittlerweile schon recht viele Änderungen gemacht habe, würde mich ein Datenverlust weit zurückwerfen.

Als nächstes werde ich mich mit der Logik für den Dateinamen beschäftigen. Auch da finde ich für meinen Geschmack zu viel Code vor. Es darf nicht schwer sein, aus einer Planetennummer und einem Datum das relevante Ephemeridenfile zu ermitteln.

Was mich noch stört: Das Durchschleusen fataler Fehler durch die Aufrufhierarchie. Ein fataler Fehler ist es z.B., wenn ein Ephemeridenfile nicht geöffnet werden kann. Nur um diesen Fall abzubilden, enthät der Aufruf jeder Funktion lästigen Boilerplate-Code: Nämlich das Entgegennehmen des Return-Codes und das Verlassen der aktuellen Stack-Ebene, falls es sich um einen fatalen Fehler handelt.

Alternative wäre der C-Vorläufer zu throw: Die Funktion longjmp(), mit der ich bei fatalen Fehlern direkt zu der höchstmöglichen Stackebene springen kann. Dies werde ich gewissermassen als eine private Version von throw für sweph.c einführen.

Die globale Datenstruktur struct swe_data swed, die so ziemlich alles enthält, was in sweph.c Rang und Namen hat (und auch so manche No-Name-Daten, die man vielleicht guten Gewissens herausschmeissen kann), wird nun noch auf oberster Stackebene verwendet - beim Einsprung aus den API-Funktionen swe_calc etc., die aus Kompatibilitätsgründen nicht die Möglichkeit bieten können, einen Zeiger auf diese Struktur als Parameter mitzugeben. Ansonsten gibt es in meiner Version von sweph.c (also in test_calc_reduce.c) keine direkten Zugriffe mehr auf diese globale Struktur.

Mittwoch, 16. November 2011

Beginn der Refaktorisierungen

Um einmal zu probieren, wie weit man mit Refaktorisierungen kommen kann, habe ich mir, wie bereits beschrieben, das Programm sweph.c vorgenommen, das allen wesentlichen Code für die Ephemeridenberechnung auf Basis der Swiss Ephemeris enthält.

Um das Programm zunächst mit automatischen Tests abzusichern, habe ich einen Testgenerator testgen entwickelt, der nun aus Dateien wie der folgenden
TESTCASE
description:Swiss Eph, geocentric positions
planets:0-13,15-20,21,22,40,43,47
flags:2,258
dates:1.1.1000-31.12.2000,5,randomly # 5 random dates for each object

TESTCASE
description:Swiss Eph, helio/barycentric positions
planets:1-9,14-20,40,43,47
flags:10,266,16642
dates:1.1.1000-31.12.2000,5,randomly # 5 random dates for each object

eine binäre Datei mit einigen Hundert Referenzdaten berechnet - nämlich die Ergebnisse der swe_calc()-Aufrufe durch die unveränderte swedll32.dll (mit den durch die Fixture definierten einzelnen Aufrufen).

Für meine Kopie von sweph.c mit dem Namen test_calc_reduce.c, mit der ich testhalber die Refaktorisierung ausprobieren möchte, habe ich ein eigenes make-Target eingerichtet, an dessen Schlusspunkt der Aufruf der Tests steht (insgesamt sind es 585 Berechnungen verschiedener Objekte zu verschiedenen Zeiten und Koordinatensystemen). Im Fehlerfall weiss ich, dass meine letzten Änderungen etwas zerstört haben müssen. Die umgekehrte Garantie, dass bei Fehlerfreiheit auch alles wie früher funktioniert, habe ich zwar leider nicht. Aber das Streben nach einer 100%igen Testabdeckung würde mich zuviel Zeit kosten. Lieber füge ich in die obige Fixture noch ein paar weitere Testfälle ein, wenn ich auf einen bislang noch unentdeckten Fehler stossen sollte.

Was fällt auf bei der Betrachtung von sweph.c? Die Datei enthält 59 Funktionen bei 6241 Zeilen, also im Schnitt 105 Zeilen pro Funktion. Die Verteilung ist allerdings sehr ungleichförmig: Die Funktion swecalc() hat beispielsweise fast 600 Zeilen! Clean Code-Profis wie Robert C. Martin empfehlen Funktionsgrössen in der Grössenordnung von 10 Zeilen und darunter. Das sind Erfahrungswerte.

Wenn die Zeilenzahlen deutlich über 10 liegen, erschwert das die Lesbarkeit und das Verständnis des Programms. In der Regel sind alle möglichen Übel die Folge: Lokale Variablen werden für aufeinderfolgenden mehrfachen Gebrauch zweckentfremdet, die Verzweigungslogik wird nicht mehr durchschaubar - und, am unangenehmsten: identische Aufgaben werden mehrfach codiert (Verstoss gegen das elementare Prinzip Don't Repeat Yourself).

Die Funktionsnamen sind zum grossen Teil gut gewählt, man versteht, was die Funktion tun soll. Für die Variablen gilt das nur noch zum Teil, und fast vollends ins Unlesbare taucht man bei der Frage nach der Bedeutung der einzelnen Komponenten der Strukturen swe_data ab. Es wäre zu prüfen, welche Komponenten dieser Struktur überhaupt noch benötigt werden.

In sweph.c finde ich verschiedene Dinge, die mich zu Verbesserungen reizen:

  • Micro-Wiederholungen, z.B. immer wieder das Kopieren von Vektoren mit einer Mini-for-Schleife,

  • Wiederholungen kleinerer Codeblöcke,

  • Immer wieder Entscheidungen anhand von Flagleisten, die die Lesbarkeit erschweren.

  • Sprünge innerhalb von Funktionen mittels goto - es gibt insgesamt 97 goto's in sweph.c - die immer der deutliche Ruf nach Auslagerung von Code in eine eigene Funktion sind.

  • Eine klarere Trennung nach Schichten täte dem Aufbau gut. Wenn man Aufrufe puffern möchte, sollte diese Pufferung möglichst "weit oben" erfolgen – oder "transversal"; jedenfalls so, dass der Teil, der die eigentliche, die ungepufferte Rechnung durchführt, nicht durch die Abfragen nach Pufferung verdunkelt wird.

  • Konstanten täte es gut, wenn sie wirklich als const deklariert und somit dem Compiler bekanntgemacht werden. Das bringt eine Menge Vorteile. Das ist aber für den riesigen Haufen an #define's in diesem Projekt wohl nicht (mehr) praktikabel und auch nicht (mehr) sinnvoll. Für einzelne Konstanten und insbesondere für neu eingeführte Konstanten werde ich die Umwandlung jedoch vornehmen.



Es gibt verschiedene Präfixe wie swe_, se_, swi_. Für die Kenntlichmachung der externen Funktionen sehe ich die Nützlichkeit eines solchen Präfixes ein, das nach aussen, im Code des Aufrufers, wie ein "Namensraum" fungiert. Aber für interne Funktionsnamen hat es sich eingebürgert, keine besonderen Präfixe mehr zu verwenden. Dasselbe gilt für Variablennamen. Es empfehlen sich einfach sprechende, ruhig lange Funktionsnamen, die wie ein Kommentar ihres Inhalts funktionieren - ein Kommentar, der diese Funktion immer begleiten wird!

Die dateiinternen Funktionen static zu deklarieren, ist sicher eine gute Idee. Es ist das private des C-Programmierers!

Insgesamt habe ich beim Lesen den Eindruck: "So schwer kann das doch alles nicht sein!" Sicher, die Ephemeridenrechnung ist kein Deckchensticken. Aber auch komplizierteste Probleme lassen sich in kleine Teilaufgaben zerlegen, die jede für sich in klarer, nachvollziehbarer Weise hingeschrieben werden können. "Was man überhaupt sagen kann, das kann man auch klar sagen", sagt Wittgenstein. Das gilt auch für Programme.

Micro-Funktionen werden vom Compiler meist automatisch in den Code hineingeneriert, es entsteht also bei der Ausführung keine neue Stackebene. In C99, das ich verwende, kann man dem Compiler dies mit dem Qualifier inline nahelegen. Solche Funktionen erzeugen keinen Performance-Overhead. Um der Lesbarkeit des Code willens kann man daher ruhig grosszügig solche Funktionen einführen - das Programm wird dadurch um keine Nanosekunde langsamer und gewinnt in der Regel an Lesbarkeit.

Einige Kandidaten für solche Inline-Funktionen wären [1],[2]:

// Utility functions (might move to swephlib)

static inline void* vec_clear( void* x, int size) {
memset((void*)x,0,size*sizeof(double));
return x;
}

static inline void* vec_copy( void* to, void* from, int n) {
memcpy(to,from,n*sizeof(double));
return to;
}

static inline void* vec_add( void* to, void* from, int n) {
double *fto = (double *) to,
*ffrom = (double *) from;
for (int i=0;i<n;i++) fto[i] += ffrom[i];
return fto;
}

static inline void* smul( void* to, double s, int n) {
double *fto = (double *) to;
for (int i=0;i<n;i++) fto[i] *= s;
return fto;
}

static inline char* msg_append(char* serr, char* text) {
int maxLengthAppendable = AS_MAXCH - strlen(serr);
if (maxLengthAppendable > 0) {
if (strlen(text) >= maxLengthAppendable)
*(text+maxLengthAppendable-1) = '\0';
strcat(serr, text);
}
return serr;
}

static inline bool is_set( int iflag, int flags ) {
return iflag & flags;
}

static inline bool is_not_set( int iflag, int flags ) {
return ! (iflag & flags);
}

static inline void set( int* iflag, int flags ) {
(*iflag) |= flags;
}

static void set_ephemeris( int* iflag, int ephemeris) {
*iflag = ( (*iflag) & ~SEFLG_EPHMASK ) |
( ephemeris & SEFLG_EPHMASK )
;
}


[1] Das void* bei den Vektorfunktionen habe ich mit Bedacht gewählt, weil ich überlege, ob der Code an manchen Stellen nicht mit einer geschachtelten Struktur von 2*3 Doubles lesbarer würde, anstelle eines Double-Arrays, etwa so:
// Using named members instead of a double[6] 
// might lead to more readable code in some cases
typedef struct {
double x;
double y;
double z;
} rect_coord;

typedef struct {
double l;
double b;
double r;
} polar_coord;

typedef struct {
rect_coord pos;
rect_coord speed;
} rect_data;

typedef struct {
polar_coord pos;
polar_coord speed;
} polar_data;


Ein void*-Argument für die Vektorfunktionen hat den Vorteil, dass man ebensogut auch einen Zeiger auf eine solche Struktur übergeben könnte. Das (compiler- und plattformspezifische) Alignment von Strukturen spuckt uns nicht in die Suppe, weil die Struktur aus lauter double's besteht. Das führt, soweit ich das erforscht habe, auf keiner Plattform zum Padding.

[2] Die Präfixe vec für die Vektorfunktionen widersprechen nicht dem oben über Präfixe Gesagten, denn es sind keine Präfixe aufgrund von generellen Namenskonventionen, sondern Teile des konkreten Namens, genau wie bei set_ephemeris()

Montag, 7. November 2011

Ein "paralleler" Zugang

Parallel zum Testframework habe ich begonnen, mir das Include sweph.c vorzuknöpfen, das den Löwenanteil der Funktionalität für die Swiss Ephemeris enthält: Es gibt nur Absprünge in die Moshierfunktionen (swemmoon, swemplan), wenn Moshier explizit angefordert wurde oder das Datum aus dem Range der Swiss Ephemeris herausreicht, einige allgemeine Bibliotheksfunktionen (swephlib) und einen Aufruf von swe_revjul (swedate) bei der Namensbildung der Ephemeridendatei. Alles andere, was zur Berechnung eines Planeten mit Swiss Ephemeris nötig ist, ist in sweph.c enthalten.

Ich habe mir daher sweph.c in ein neues Programm test_calc_reduce.c hineinkopiert, eine main()-Funktion vorangestellt und im Makefile die Compilierung ohne Verwendung von sweph.o definiert:


test_calc_reduce.exe: test_calc_reduce.o
$(CC) $(CFLAGS) -DUSE_DLL -o test_calc_reduce.exe test_calc_reduce.o
swemmoon.o swemplan.o swephlib.o swedate.o
test_calc_reduce.o: swephexp.h swedll.h sweodef.h test_calc_reduce.h

Nun arbeite ich mich von oben durch die Funktionsebenen durch, um die Programmideen klarer hervortreten zu lassen und die Abhängigkeiten mit globalen Daten zu entfernen oder zu untersuchen. Wichtig sind auch hier Tests, um nach jeder (auch kleinen) Refaktorisierung sicherzustellen, dass die grundlegenden Aufrufe alle noch funktionieren. Im Moment mache ich das noch von Hand, plane aber einen Querschnittstest, der für jeden mögliche Berechnungstyp ein oder zwei Beispiele enthält. Eine interessante Aufgabe.

Die switch/cases, Sprungmarken und goto's sind, ebenso wie Kommentarblöcke, die mit einer langen Reihe von Sternchen beginnen, Hinweise darauf, dass man die Lesbarkeit verbessern kann, wenn man den nachfolgenden Code in eine parametrisierte Funktion auslagert. Der Stack ist ein grosser Helfer: Er hilft, die Geltung von Variablen auf kleinere Codeblöcke zu begrenzen. Sie werden einem zu Beginn der Routine automatisch auf- und beim Verlassen automatisch abgebaut. Änderungen an Parametern nimmt die Funktion wegen der Wertübergabe nur vor, wenn ich das explizit wünsche. Ansonsten sind nach dem Rücksprung die Daten automatisch in dem Zustand, den sie vor dem Aufruf hatten. Auch das ist ein hilfreiches Systemverhalten.

Die verschiedenen Arten der Berechung - z.B. Moshier-Mond, Asteroid, Ekliptikschiefe, habe ich in einzelne Funktionen der jeweils gleichen Signatur ausgelagert. Nach dem Strategiemuster vorgehend, kann ich aus iflag und ipl in einer separaten Funktion die Funktionsreferenz ermitteln, die für die Berechnung zuständig ist. Dadurch vereinfacht sich der Code der aufrufenden Funktion (swecalc()) erheblich. Es wird nicht sehr viel weniger Code, aber was übrig bleibt, ist klarer einzelnen Aufgaben zuzuordnen.

Was die globalen Daten angeht, so scheint das meiste in der globalen Struktur swed enthalten zu sein. Intern wird neu allen Funktionen, die auf swed zugreifen, ein Zeiger auf eine Struktur vom Typ swe_data übergeben, und die Zugriffe erfolgen über diesen Zeiger. So bereite ich vor, dass diese Daten im Stack verwaltet werden können statt im globalen Datenbereich.

So schlecht das Schlüsselwort static für Variablen ist, weil es neue versteckte globale Daten einführt, so hilfreich ist es für die Kapselung, weil es die Sichtbarkeit von Variablen und Funktionen auf ein File einschränkt. Eine mit static deklarierte Variable ist immer noch besser als eine überall sichtbare globale Variable. In dieser Hinsicht ist der static-Qualifier vergleichbar dem private in objektorientierten Sprachen.

Donnerstag, 3. November 2011

1000 Testrechnungen in 80 Millisekunden

Die Bresche ist nun geschlagen:

  • Ein Programm testgen – das mit einer Referenzversion der swedll32.dll aus dem Release 1.77 arbeitet – berechnet zu zehn verschiedenen astrologischen Planeten und je 100 verschiedenen Daten die geozentrischen ekliptikalen Positionen und speichert die Ergebnisrecords binär (mittels fwrite) in einer Datei testdata.bin ab.
  • Das Programm tests enthält eine Testklasse TestGeneratedData, die diese Daten einliest, die Berechnungen mit der aktuell in Bearbeitung befindlichen DLL ausführt und die Ergebnisse mit den Erwartungen vergleicht. Abweichungen werden in detaillierter, nachvollziehbarer Form in eine Logdatei testFailures.txt geschrieben. Für 1000 Rechnungen benötigt das Testprogramm 80 Millisekunden. Eine erträgliche Zeit, um sie bei Builds mit auszuführen.

In der Datei folgen ein allgemeiner testHeader, ein Header swissCalcHeader für die Tests, die speziell der Funktion swe_calc gewidmet sind, sowie eine Anzahl von Records swissCalcData - für jede auszuführende Rechnung einer.
// A test
typedef struct {
// Specifies the particular kind of the test (e.g. a test of swe_calc):
int type;
char description[100];
// The number of bytes belonging to this test (for skips):
int length;
} testHeader;

// A set of swe_calc computations
typedef struct {
// flags for swe_calc
int flags;
// precision exponent 10**p as factor to DEFAULT_PRECISION :
int precision[6];
int numberOfRecords;
} swissCalcHeader;

// A single swe_calc computation
typedef struct {
double jd;
int planet;
double result[6];
} swissCalcData;

Hier der Code der Testklasse für die Funktion swe_calc:
class TestGeneratedData : public Test {
public :
virtual void run() {
testHeader header;
swissCalcHeader scHeader;
FILE *testdata = fopen( "testdata.bin","rb" );
if (!testdata) {
fail( "Can't open testdata.bin" );
}
fread( &header, sizeof(header), 1, testdata );
fread( &scHeader, sizeof(scHeader), 1, testdata );

swissCalcData scData[scHeader.numberOfRecords];
fread( &scData,
sizeof(swissCalcData),
scHeader.numberOfRecords,
testdata );
SwissCalc sc;
bool ok = true;
for (int i=0;i<scHeader.numberOfRecords;i++) {
if (!sc.calcsAsExpected(
scData[i].jd, scData[i].planet, scHeader.flags, scData[i].result
)) {
fail( "Differences", sc.diffs );
}
}
fclose( testdata );
}
};

Der File I/O muss noch in andere Funktionen ausgelagert werden (anderer Concern) und gegen Fehler abgesichert werden (Rückgabewerte von fread auswerten). Die eigentliche Rechenarbeit steckt in der Funktion SwissCalc::calcsAsExpected, die ich ja früher schon diskutiert hatte.
Noch offen:

  • Das Programm bricht momentan bei der ersten Rechnung ab, die Differenzen ergab. Der Programmfluss soll so geändert werden, dass trotz Differenzen weitergerechnet wird: Es soll alles protokolliert werden, und am Schluss bekommt dieser Test in der Standardausgabe den Status not ok.
  • Damit grössere Datenmengen eingelesen werden, soll eine Paketbildung gemacht werden, wonach die Daten z.B. in 1000er Päckchen eingelesen und ausgeführt werden.
  • Der Name der Testdatendatei soll im Programm als Parameter angegeben werden können. Er wird dann vom Testrunner an die Testklasse durchgereicht (oder die Testklasse kann diesen Dateinamen beim Testrunner erfragen).



Weiteres Vorgehen: Nun beginnt die eher systematische Arbeit, viele Testfälle zu erfassen. Die Tests sollten eine Zeitreise durch die Jahrhunderte machen und auch weniger verwendete (Klein-)Planeten umfassen. Die verschiedenen Koordinatensysteme, die man anfordern darf, sollten ausgetestet werden - ebenso andere Optionen wie das Speed-Flag.

Solange ein solcher Rundumschlag durch die möglichen Aufrufe noch fehlt, sind Änderungen am Quellcode tabu (sosehr es mich danach gelüstet, schon jetzt Hand anzulegen...).

Montag, 24. Oktober 2011

Differenzenermittlung - Problem ist gelöst

Das Problem mit den Differenzen ist nun gelöst. Das folgende Template kann eine beliebige einzelne, benannte Abweichung festhalten und formatiert ausgeben.

Aus Clientsicht vereinfacht sich der Aufruf erheblich: Hier ein Beispiel, der Aufruf einer Planetenberechnung via swe_calc(). Beim Aufruf gibt man die Argumente sowie die erwarteten Rückgabewerte an. Wenn die Methode calcsAsExpected feststellt, dass einige Werte abweichen, gibt sie false zurück. Man kann dann das Member diffs verwenden, um alle festgestellten Differenzen auszugeben.
class TestSwissEphemeris : public Test {
public :
virtual void run() {
SweCalc se;
if (!se.calcsAsExpected(
2451545.0,
0, // <-- SE_SUN
0, // <-- Flags
(double[6]) { 280.368, 0.000227827, 0.983328, 0.,0.,0. })
) {
fail( "Differences for Sun", se.diffs );
}
}
};


Dieser Teil des Programms wird automatisch generiert werden, wobei dann die Erwartungswerte mit hineinkommen. Auch der Text (hier Differences for Sun kann vom Automaten so produziert werden, dass man nachher bei Betrachtung des Logs den Fall nachvollziehen kann. Es könnte z.B. eine swetest-Kommando sein, mit dem man den Unterschied reproduzieren kann.

Hier die oben verwendete Klasse, die einen Aufruf von swe_calc() ausführt und die Differenzen sammelt:

class SweCalc : public TestCalculation {
public :
bool calcsAsExpected(
double juldate,
int planet,
int flags,
const double expectedResult[],
const string expectedMessage = "" ) {

double actualResult[6];
char actualMessage[255];

int return_flags = swe_calc(
juldate,
planet,
flags,
actualResult,
actualMessage );

diffs.check( "xx", 6, actualResult, expectedResult );
diffs.check( "msg", actualMessage, expectedMessage );
diffs.check( "flags", return_flags, flags );

return ! hasDiffs();

}
};


Was für andere Funktionsaufrufe der Swiss Ephemeris wiederverwendbar ist, kommt in eine Oberklasse TestCalculation. Im Moment stehen hier nur die Differenzen-Sätze:

class TestCalculation {
public :
bool hasDiffs();
Diffs diffs;
};


Die ganz allgemeine Schicht, die diese Prüfungen durchführt und in der ich den Fehler hatte, folgt hier:

template<typename T> 
class Diff : public DiffBase {
public :
Diff(string name, T actual, T expected)
: actual(actual), expected(expected) {
this->var = new Var(name);
}
Diff(string name, int index, T actual, T expected)
: actual(actual), expected(expected) {
this->var = new VarMember(name,index);
}
virtual string toString() const {
stringstream s;
s << var->getName()
<< "\t: actual=" << actual
<< ", "
<< "expected=" << expected ;
return s.str();
}
~Diff() {
delete var;
}
private :
Var* var;
T actual;
T expected;
};

Die Klasse, die das Problem machte, war Diffs - eine Sammlung mehrerer solcher Abweichungen. Da ich sie auch im System herumreichte, hat mir der Destruktor in den Fuss geschossen, den ich zuerst vorgesehen hatte. Nach Löschung des Destruktors und Umstellung auf einen boost::shared_ptr (der das delete ja automatisch ausführt, wenn der letzte Besitzer des Pointers sein Leben aushaucht) lief alles problemlos.

class Diffs {
public :
Diffs() {
diffs = boost::shared_ptr<vector<DiffBase*>>
(new vector<DiffBase*>);
}
void getLog( vector<string>& log ) const {
log.clear();
for (vector<DiffBase*>::const_iterator i=diffs->begin();
i!=diffs->end();
++i) {
log.push_back( (*i)->toString() );
}
}
bool hasDiffs() { return diffs->size() > 0; }
template<typename T> void check( const string& name,
const T actual,
const T expected ) {
if (actual != expected) {
diffs->push_back( new Diff<T>(name,actual,expected) );
}
}
template<typename T> void check( const string& name,
const int size,
const T actual[],
const T expected[] ) {
for (int i=0; i<size;i++) {
if (actual[i] != expected[i]) {
diffs->push_back( new Diff<T>(name,i,actual[i],expected[i]) );
}
}
}
string toString() const {
string s;
for (vector<DiffBase*>::const_iterator i=diffs->begin();
i!=diffs->end();
++i) {
s.append((*i)->toString());
s.push_back('\n');
}
return s;
}
private :
boost::shared_ptr<vector<DiffBase*>> diffs;
};


Ein isoliert lauffähiges Beispiel habe ich codepad.org gestellt, wo es auch ausgeführt werden kann (ist sogar ISO-C++, obwohl ich mir für das Testprogramm immer den "Standard" C++0x erlaube).

Auch heute gibt es ein
Todo: Für double sollte ein Test mit einer vorgegebenen Genaugkeit durchgeführt werden. Wie ist das in das Template einzubauen?

Sonntag, 23. Oktober 2011

Problem bei Differenzen-Ermittlung

Die automatischen Tests werden von einem in C++ geschriebenen Programm namens tests durchgeführt werden. Teile dieses Programms werden anhand von erwarteten Testdaten, die z.B. als Ouptput vom bestehenden Programm swetest produziert werden, generiert.

Bin gerade dabei, ein System von Klassen zu schreiben, das über allfällig festgestellte Abweichungen von Ist- und Sollwerten Buch führt: Die entscheidende Klasse Diffs ist im wesentlichen eine Sammlung (Vektor) von einzelnen Abweichungssätzen. Eine einzelne Abweichung ist vom Typ IntegerDiff, DoubleDiff, DoubleMemberDiff, StringDiff o.ä. - alles Subtypen einer gemeinsamen Oberklasse Diff. Gibt es Abweichungen, so wird eine solche Diffs-Instanz als Member der entsprechenden Failure-Ausnahme an den TestRunner gesendet. Dieser protokolliert in der Standardausgabe die Summenzeile des TestAnythingProtocols und gibt die Details in ein Logfile aus.

ToDo - Code-Duplizierungen verhindern durch Templates: Diff<double>, Diff<int> etc.


Hier gibt es bei der Vererbung noch das Problem, dass mir der Vektor irgendwo auf eine Weise kopiert wird, durch die die Typinformationen der aktuellen Typen verlorengeht und alle Einzelabweichungen nur noch den Typ der gemeinsamen Oberklasse Diff haben.

Das Vorhaben wird dargelegt

Vor kurzem fasste ich den Entschluss, die Bibliothek Swiss Ephemeris so umzuschreiben, dass sie nebenläufig lauffähig ist.

Dies ist ein grösseres Refactoring-Projekt, das folgende Schritte beinhaltet:

  • Bevor auch nur eine Zeile Code geändert wird, muss ein dichtes Netz von Tests bereitstehen, die automatisch, z.B. mit einer bestimmten make-Regel, alle relevanten Funktionen im aktuellen Stand der Bibliothek aufrufen und mit den erwarteten Ergebnissen vergleichen.

  • Zwar ist es theoretisch möglich, globalen Speicher in mehreren Threads wiederzuverwenden, der gelesen und geschrieben wird. Wo es möglich ist, ist es jedoch besser, wenn der Aufrufer den Speicherplatz für einen allfällig benötigten Kontext bereitstellt (also Daten "auf dem Stack statt auf dem Heap" verwendet werden). Ausnahme sind natürlich Ressourcen, die wirklich von allen Threads in identischer Form benötigt werden (z.B. die Daten von Ephemeridenfiles). Es soll ein Set von API-Funktionen angeboten werden, die genau dies ermöglichen.

  • Die bestehenden Funktionen werden an die neuen API-Funktionen angeschlossen (die dann den benötigten Speicher intern allokieren), so dass Kompatibilität mit den bisherigen Verwendungen der Swiss Ephemeris erreicht wird.


Angesichts der vielen Arbeit, die in diese Software gesteckt wurde, und angesichts der geringen Freizeit, die mir dafür zur Verfügung steht, habe ich mir für dieses Vorhaben rund ein Jahr Zeit gegeben. Ich kann damit scheitern - oder auch nicht.

Ein Blog scheint mir eine passende Form, um die kontinuierliche Arbeit an diesem Thema zu dokumentieren. Auch diesen Blog schreibe ich - wie meinen persönlichen Blog eigentlich nur für mich selbst: Durch das systematische Aufschreiben von Themen bekomme ich Klarheit in meinem Gedankenkasten. In dieser Funktion ist der Blog identisch mit einem gewöhnlichen Tagebuch. Abweichend ist nur die elektronische und überall verfügbare Form - das Medium.

Dieser Blog bleibt bis auf weiteres nur einem begrenzten Nutzerkreis zugänglich, den ich persönlich einlade.