Numerische Bibliothek in C++

headmyshoulder

knivil schrieb:

...

Oha, find ich ja gut, daß es gleich konkret wird. Danke für die Kritik, viele der Sachen die Du angespricht, haben wir auch schon diskutiert.

Geschwindigkeit: Schau Dir mal die Speedtests in https://yanl.svn.sourceforge.net/svnroot/yanl/trunk/examples/odeint an. Das Klassendesign mit virtuellen Funktionen ist Faktor 2 langsamer^^. Sry, falls der Code dort leicht unübersichtlich ist.

Zum Design: Keiner hindert Dich dadran trotzdem Klassen zu verwenden. Hauptsache der (,,) operator ist definiert. Zum Beispiel kannst Du das obige Beispiel auch mit einer Klasse schreiben:

class Lorenz
{
public:

    void operator()( state_type &x , state_type &dxdt , double t )
    {
	dxdt[0] = sigma * ( x[1] - x[0] );
	dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
	dxdt[2] = x[0]*x[1] - b * x[2];
    }
};

und das sieht dann angewandt so aus:

ode_step_euler< vector<double> > euler;
Lorenz lorenzo;
euler.next_step( lorenzo , x , t , dt );

Das ist doppelt so flexibel als ein Klassenansatz, weil man Funktion und auch Klassen benutzen kann.

Iteratoren oder Container: Ich denke das es mindestens eine Variante mit Iteratoren geben wird. Das obige Beispiel in dem anderen Thread ist nicht die Endfassung. Die Diskussion ist noch nicht zu Ende, beide Varianten haben Vor- und Nachteile. Wir können diese Diskussion ja woanders hinverschieben. Der grosste Nachteil ist, daß man die Anzahl der Elemente zwischen zwei Iteratoren nicht notwendigerweise mit O(0) rauskriegt. Bei Listen dauerts beispielsweise O(n). Falls irgendwann der Stepper weiterbenutzt werden soll, muß dort dann halt ein template Parameter stehen, statt einer abstrakten Klasse. Der entsprechende Typ muss die Stepper Eigenschaften erfüllen, genauso wie das dynamische System den (,,) Operator haben muss.

Schrittweitensteuerung Es gibt noch gar keine Schrittweitensteuerung. Fakt ist aber, daß falls die Schrittweite gefunden wurde, der Schritt auch ausgeführt werden soll und dafür ist die Funktionalität da mit beispielsweise ode_step_euler . Die eigentliche Iteration soll dann auch so ähnlich wie Du vorgeschlagen hast aussehen, nur daß noch Zustandsbeobachter hingezugefügt werden, welche in jeder Iteration irgendwas mit dem Zustand machen, Ausgabe, Berechnungen, Mitteln...

knivil

Sieht interessant aus, vom Ansatz her. Ausprobiert habe ich den Code noch nicht, nur gelesen. Aber etwas zum Codestyle: Beginnt man nach einer geschweiften Klammer eine neue Zeile, so rueckt man ein. Das macht das Lesen angenehmer. Ihr koennt euch ja mal einige Syleguides durchlesen.

Auch solltet ihr euch um geeignete Matrix und Vektorklassen bemuehen, die solche Zeilen kapseln:

for( size_t i=0 ; i<dxdt.size() ; ++i ) state[i] += dt * dxdt[i];

und in etwas Vertrauteres verwandeln:

x_next = x +  dt * dx;

Ach verdammt, dann waere man nicht mehr containerunabhaengig. Auch scheint dieses Beispiel nicht mit std::list als Container zu funktionieren, da ihr random access auf die Elemente benoetigt ... getestet habe ich es nicht.

Der grosste Nachteil ist, daß man die Anzahl der Elemente zwischen zwei Iteratoren nicht notwendigerweise mit O(0) rauskriegt.

Wenn man die Methode size einer Liste aufruft, dann garantiert der Standard auch keine O(1) Komplexitaet, sondern auch nur O(n). Moderne Listenimplementationen von std::list achten aber auf ihre Groesse, so dass es zumindest fuer diesen Datentyp keine Probleme gibt. Bei RandomAccessIteratorn sollte es moeglich sein, die Anzahl der Elemente dazwischen in O(1) zu bestimmen.

Das Klassendesign mit virtuellen Funktionen ist Faktor 2 langsamer

Das kann gut sein, wenn es sich um eher "kleine" Differentialgleichungen handelt. Eine virtuelle Funktion in Einfachvererbung kostet halt eine Zeigerdereferenzierung mehr. Ist das im Vergleich zum Rest der Funktion viel, so kann sich das nachteilig auf die Gesamtperformance niederschlagen. Die Funktion wird ja staendig ausgewertet. Die Kosten sind aber pro Funktionsauswertung konstant. Auch hat der Compiler beim Templateeinsatz mehr Moeglichkeiten zu optimieren.

Ich werde mal etwas mit dem neuen Ansatz experimentieren.

headmyshoulder

knivil schrieb:

Sieht interessant aus, vom Ansatz her. Ausprobiert habe ich den Code noch nicht, nur gelesen. Aber etwas zum Codestyle: Beginnt man nach einer geschweiften Klammer eine neue Zeile, so rueckt man ein. Das macht das Lesen angenehmer. Ihr koennt euch ja mal einige Syleguides durchlesen.

Die aktuellen Sourcen befinden sich in https://boost.org/svn/boost/sandbox/odeint nicht daß es da Verwechslungen mit dem alten Code gibt. Dort haben wir auch auf vernünftige Konvertierung geachtet, hoffe ich zumindest.

knivil schrieb:

Auch solltet ihr euch um geeignete Matrix und Vektorklassen bemuehen, die solche Zeilen kapseln:
for( size_t i=0 ; i<dxdt.size() ; ++i ) state[i] += dt * dxdt[i];
und in etwas Vertrauteres verwandeln:
x_next = x +  dt * dx;
Ach verdammt, dann waere man nicht mehr containerunabhaengig. Auch scheint dieses Beispiel nicht mit std::list als Container zu funktionieren, da ihr random access auf die Elemente benoetigt ... getestet habe ich es nicht.

Der code in der boost sandbox sollte mit std::list laufen, man kann halt in der ODE Funktion (oder Klasse) dann keine RandomAccess-Iteratoren benutzen.

knivil schrieb:

Wenn man die Methode size einer Liste aufruft, dann garantiert der Standard auch keine O(1) Komplexitaet, sondern auch nur O(n). Moderne Listenimplementationen von std::list achten aber auf ihre Groesse, so dass es zumindest fuer diesen Datentyp keine Probleme gibt. Bei RandomAccessIteratorn sollte es moeglich sein, die Anzahl der Elemente dazwischen in O(1) zu bestimmen.

Joa, da hast Du Recht. Ich persönlich würde es am besten finden wenn beide Varianten unterstüzt werden. Also einmal, daß die ODE Klasse mit Containern aufgerufen wird, und einmal mit Iteratoren. Allerdings muss man dann wahrscheinlich von jedem Algorithmus zwei Varianten schreiben, einem der die ODE mit Containern aufruft und einem der die ODE mit iteratoren aufruft, oder man abstrahiert das irgendwie weg.

knivil schrieb:

Ich werde mal etwas mit dem neuen Ansatz experimentieren.

Schön, wenn Du Ergebnisse oder Anregungen/Kritik/Vorschläge hast, lass es uns wissen.

knivil

Jap, hatte die anderen Sourcen aus dem uebergeordneten Ordner mit den Speedtests.

PS: Ich habe gerade festgestellt, dass die Welt klein ist. Ich bin auch oefters im Goldenen Kaefig von Golm ... Aber ich dachte, dass es bereits ein Programmpaket fuer Statistik, Odes und den ganzen Kram in der Nichtlinearen Dynamik entwickelt wurde.

headmyshoulder

knivil schrieb:

Jap, hatte die anderen Sourcen aus dem uebergeordneten Ordner mit den Speedtests.

PS: Ich habe gerade festgestellt, dass die Welt klein ist. Ich bin auch oefters im Goldenen Kaefig von Golm ... Aber ich dachte, dass es bereits ein Programmpaket fuer Statistik, Odes und den ganzen Kram in der Nichtlinearen Dynamik entwickelt wurde.

Hehe, ja manchmal ist die Welt wirklich klein.

Ich glaub, hier bei uns wurden schon mehrere solcher und ähnlicher Pakete entwickelt. In der boost gibt es sowas noch nicht, aber ich denke, die Chancen das dort reinzubekommen sind nicht so schlecht. Es wurde mal eine Zeitreihenbibliothek entwickelt, aber die hat es nicht bis zur Veröffentlichung geschafft. Naja, der Sinn ist auch endlich mal was zu schaffen, was Standards benutzt. Da gibt es echt nicht viel. Vielleicht benutzen dann auch wieder mehr Leute C++ statt Python. Das verbreitet sich in der Wissenschaft gerade rasant^^.

FelixKlein

Hallo headmyshoulder!

Ich habe erst gestern das entsprechende Kapitel in den Num. Recipes gelesen und bin heute auf deine Idee von dem C++ Odeint Framework gestoßen. Ich hoffe das Projekt läuft noch!

Ich selber habe Mathematik studiert. Mit der ODE-Solvern und Schrittweitensteurungen beschäftige mich erst seit Kurzem. Auch in C++ bin ich noch relativ unerfahren. Aber ich werde muss mich sowieso auf beiden Gebieten in nächster Zeit einarbeiten. Über jeden Austausch und gegebenenfalls eine Zusammenarbeit würde ich mich freuen!

Viele Grüße
FelixKlein

headmyshoulder

Hi FelixKlein,

das odeint Projekt läuft noch und ist auch schon recht fortgeschritten, du hast quasi eine menge aufzuholen. Zu finden ist das Ganze in der boost sandbox

https://www.boost.org/svn/boost/sandbox/odeint

Akutell könnten wir die Dormand-Prince Methoden gebrauchen. Wenn Du Bock kannst Du dir ja ansehen wie unsere lib funktioneirt und dann die Methoden implementieren.

Viele Grüße,

headmyshoulder

FelixKlein

Also DOPRI interessiert mich persönlich weniger.
Habt ihr euch eigentlich schon Gedanken über Extrapolationsverfahren gemacht?

In jedem Fall muss ich erst einmal noch einiges lesen. Erst wenn ich die theoretischen Grundlagen durchblicke will ich mir eure Implementierungen im Detail ansehen und erst wenn ich die verstehe mache ich mich selber ans Werk.

Soweit ich das bisher abschätzen kann ist euer Projekt sehr strukturiert und sauber programmiert. Ob ich mit dem Code tatsächlich etwas anfangen kann weiß ich noch nicht. Toll jedenfalls, dass ihr das als Open Source ins Netz stellt.

Ich würde gerne weiter in Kontakt bleiben. Gerne auch persönlicher als nur über das Forum (E-Mail o.Ä.)

Gruß
FelixKlein

headmyshoulder

FelixKlein schrieb:

Also DOPRI interessiert mich persönlich weniger.

Ok, kein Problem. Es gibt noch ne Menge anderer Sachen...

FelixKlein schrieb:

Habt ihr euch eigentlich schon Gedanken über Extrapolationsverfahren gemacht?

Ja, haben wir, und das steht auch auf unserer Wunschliste. Allerdings gibt es nicht für jede Methode ein Extrapolationsverfahren, glaube ich zumindest. Deshalb haben wir das erstmal nicht mit aufgenommen. Implizite Methoden wäre prinzipiell auch sehr interessant. Allerdings braucht man dann lineare Algebra (LU-Zerlegung,...) was nicht trivial ist. Und symplektische Verfahren mit Fehlerabschätzung haben wir auch noch nicht implementiert.

FelixKlein schrieb:

In jedem Fall muss ich erst einmal noch einiges lesen. Erst wenn ich die theoretischen Grundlagen durchblicke will ich mir eure Implementierungen im Detail ansehen und erst wenn ich die verstehe mache ich mich selber ans Werk.

Welche Bücher hast Du dazu? Die von Hairer kann ich empfehlen.

FelixKlein schrieb:

Soweit ich das bisher abschätzen kann ist euer Projekt sehr strukturiert und sauber programmiert. Ob ich mit dem Code tatsächlich etwas anfangen kann weiß ich noch nicht. Toll jedenfalls, dass ihr das als Open Source ins Netz stellt.

Wir sind im Moment eigentlich zu zweit. Und das Ziel ist es, das Ganze in die boost bibliotheken zu submitten. Zu den technischen Details haben wir einen Artikel geschrieben:
http://www.codeproject.com/KB/recipes/odeint.aspx, obwohl der schon leicht veraltet ist. Ich werde den aber demnächst aktualisieren.

FelixKlein schrieb:

Ich würde gerne weiter in Kontakt bleiben. Gerne auch persönlicher als nur über das Forum (E-Mail o.Ä.)

Jo, ich schick Dir mal eine pn. Wir haben auch eine Mailingliste zu diesem Thema bei sf: https://lists.sourceforge.net/lists/listinfo/yanl-develop. Wenn Du willst können wir auch darüber kommunizieren.

Viele Grüße,

headmyshoulder

FelixKlein

headmyshoulder schrieb:

Welche Bücher hast Du dazu? Die von Hairer kann ich empfehlen.

Momentan arbeite ich mich durch "Numerische Mathematik (2)" von Bulirsch/Stoer. Ich habe aber auch schon andere Lektüre gesichtet. Unter anderem "Solving ODEs I/II" von Hairer/Norsett/Wanner. Meintest du das oder andere Bücher von ihm?

headmyshoulder schrieb:

Zu den technischen Details haben wir einen Artikel geschrieben:
http://www.codeproject.com/KB/recipes/odeint.aspx, obwohl der schon leicht veraltet ist. Ich werde den aber demnächst aktualisieren.

Den Artikel hatte ich schon entdeckt. Ist interessant. Ich freue mich auf die Aktualisierung.

headmyshoulder schrieb:

Wir haben auch eine Mailingliste zu diesem Thema bei sf: https://lists.sourceforge.net/lists/listinfo/yanl-develop. Wenn Du willst können wir auch darüber kommunizieren.

Was genau ist der Zweck der Mailingliste? Ich möchte weder vollgespammt werden noch, dass meine E-Mails öffentlich gemacht werden.

Ich melde mich heute Abend wieder per E-Mail.

Viele Grüße
FelixKlein

headmyshoulder

FelixKlein schrieb:

Momentan arbeite ich mich durch "Numerische Mathematik (2)" von Bulirsch/Stoer. Ich habe aber auch schon andere Lektüre gesichtet. Unter anderem "Solving ODEs I/II" von Hairer/Norsett/Wanner. Meintest du das oder andere Bücher von ihm?

Jepp, die Solving ODEs I/II" mein ich. Für symplektische Verfahren haben wir auch mit "B. Leimkuhler and S. Reich, Simulating Hamiltonian Dynamics." gearbeitet.

FelixKlein schrieb:

Was genau ist der Zweck der Mailingliste? Ich möchte weder vollgespammt werden noch, dass meine E-Mails öffentlich gemacht werden.

Soviel passiert im Moment nicht auf der Liste, also keine Angst vor Spam. Man kann da halt diskutieren, alles was odeint betrifft. Ist wichtig, falls mehrere Leute programmieren, daß alle halbwegs auf dem Laufenden sind. Einige Sachen sollte man bei solch einem Projekt öffentlich ankündigen.

Wenns um Details geht, kann man natürlich auch direkt via Email kommunizieren.

Bis dahin,

headmyshoulder

headmyshoulder

Nur mal als kleines Update: Was hier als Bibliothek zum Lösen von Differentialgleichungen (odeint) angefangen hat ist jetzt akzeptiert als Boost-Bibliothek. Hat ziemlich lange gedauert, und es war eine Menge Arbeit, aber es hat sich gelohnt. Und natürlich ist der Umfang enorm gewachsen. Ein Hauptvorteil ist die Unabhängigkeit der Algorithmen von arithmetischen Operationen und Containern, so dass man mit relativ einfachem Aufwand auch Differentialgleichungen auf Grafikkarten lösen kann.

rapso

ihr freaks
ich gratuliere, schoen mal ein erfolgreiches projekt zu sehen!

knivil

Ich habe die Entwicklung aus der ferne verfolgt. Habt ja auch einige Vortraege gehalten. Sieht gut aus ... wenn ich nur mehr Zeit haette.

Kellerautomat

Wo finde ich die Bibliothek?

headmyshoulder

Kellerautomat schrieb:

Wo finde ich die Bibliothek?

Die Bibliothek liegt unter odeint.com. Es wird wahrscheinlich bis Release 1.53 dauern bis die direkt in der Boost ist.