[gelöst Danke]Parallelisierung einer For-Schleife
-
Ok, ich habe dann mal die Datei umgeschrieben. Alles ist weitesgehend ausdokumentiert. Korrekt müsste eigentlich alles sein. Danke für Hilfe bei der Optimierung.
#include<iostream> #include<fstream> #include<time.h> #include<math.h> #include<vector> #include<sstream> #include<limits> #include<stdlib.h> #include<windows.h> using namespace std; class HMM { private: vector<string> data; vector<string> states; vector<string> alphabet; vector<double> start; vector<double> transition; vector<vector<double> > emission; vector<double> forward_matrix; vector<double> backward_matrix; double absolut_min ; double epsilon; int my_hash(char a, vector<string> alphabet) { for(int i = 0; i < alphabet.size(); i++) { if(alphabet[i][0] == a) { return i; } } return 0; } /* Die Logarithmen muessen Berechnet werden, damit die Rechnung auch fuer groessere Sequenzen durchgeführt werden kann. Der "find_max" - Teil dient einer Scalierung, um große Werte zu Berechnen. Dies ist notwendig und kann nicht umgangen werden. Hoechstwahrscheinlich macht dies eine Matrix-Multiplikation nicht moeglich. */ void forward_procedure() { for(int n = 0; n < this->data.size(); n++) { for(int i = 0; i < this->states.size(); i++) { this->forward_matrix[n * this->states.size() * this->data[0].size() + i] = log(this->start[i]) + log(this->emission[i][this->my_hash(data[n][0], this->alphabet)]); } for(int i = 1; i < this->data[n].size(); i++) { for(int j = 0; j < this->states.size(); j++) { double find_max = absolut_min; for(int k = 0; k < this->states.size(); k++) { if((this->forward_matrix[n * this->states.size() * this->data[0].size() + (i-1) * this->states.size() + k] + log(this->transition[k * this->states.size() + j])) > find_max ) { find_max = this->forward_matrix[n * this->states.size() * this->data[0].size() + (i-1) * this->states.size() + k] + log(this->transition[k * this->states.size() + j]);} } double tmp_prob = 0; for(int k = 0; k < this->states.size(); k++) { tmp_prob += exp(forward_matrix[n * this->states.size() * this->data[0].size() + (i-1) * this->states.size() + k] + log(this->transition[k * this->states.size() + j]) - find_max); } this->forward_matrix[n * this->states.size() * this->data[0].size() + i * this->states.size() + j] = log(tmp_prob) + find_max + log(this->emission[j][this->my_hash(data[n][i], this->alphabet)]); } } } } /* Die Logarithmen muessen Berechnet werden, damit die Rechnung auch fuer groessere Sequenzen durchgeführt werden kann. Der "find_max" - Teil dient einer Scalierung, um große Werte zu Berechnen. Dies ist notwendig und kann nicht umgangen werden. Hoechstwahrscheinlich macht dies eine Matrix-Multiplikation nicht moeglich. */ void backward_procedure() { for(int n = 0; n < this->data.size(); n++) { for(int i = 0; i < this->states.size(); i++) { this->backward_matrix[n * this->states.size() * this->data[0].size() + (this->data[n].size()-1) * this->states.size() + i] = log(1); } for(int i = this->data[n].size()-2; i >= 0; i--) { for(int j = 0; j < this->states.size(); j++) { double find_max = this->absolut_min; for(int k = 0; k < this->states.size(); k++) { if(this->backward_matrix[n * this->states.size() * this->data[0].size() + (i+1) * this->states.size() + k] + log(this->transition[j * this->states.size() + k]) + log(this->emission[k] [this->my_hash(data[n][i+1], this->alphabet)]) > find_max ) { find_max = this->backward_matrix[n * this->states.size() * this->data[0].size() + (i+1) * this->states.size() + k] + log(this->transition[j * this->states.size() + k]) + log(this->emission[k] [this->my_hash(data[n][i+1], this->alphabet)]);} } double tmp_prob = 0; for(int k = 0; k < this->states.size(); k++) { tmp_prob += exp(this->backward_matrix[n * this->states.size() * this->data[0].size() + (i+1) * this->states.size() + k] + log(this->transition[j * this->states.size() + k]) + log(this->emission[k][this->my_hash(data[n][i+1],this->alphabet)]) - find_max); } this->backward_matrix[n * this->states.size() * this->data[0].size() + i * this->states.size() + j] = log(tmp_prob) + find_max; } } } } public : HMM(int argc, char *argv[]) { // fuer Auswertung der Fliesskomma Zahlen interessant this->absolut_min = -1 * numeric_limits<double>::max(); this->epsilon = numeric_limits<double>::epsilon(); // Alphabet this->alphabet.push_back("A"); this->alphabet.push_back("C"); this->alphabet.push_back("G"); this->alphabet.push_back("T"); // Data-Strings string a = "ATCGTATCGATCGATCGCTAGCTA"; string b = "ATGCTGACTGACTGACTGACTGAC"; string c = "AAATGTGTGTCATGCATCGTACGT"; this->data.push_back(a); this->data.push_back(b); this->data.push_back(c); // Zustande this->states.push_back("X"); this->states.push_back("A"); // Start- und Transitions-Wkt. this->start.push_back(0.75); this->start.push_back(0.25); this->transition.push_back(0.9); this->transition.push_back(0.1); this->transition.push_back(0.2); this->transition.push_back(0.8); // Emissions-Wkt. vector<double> tmp (this->alphabet.size(), 0); vector<vector<double> > tmp2 (this->states.size(), tmp); this->emission = tmp2; this->emission[0][0] = 0.25; this->emission[0][1] = 0.25; this->emission[0][2] = 0.25; this->emission[0][3] = 0.25; this->emission[1][0] = 0.1; this->emission[1][1] = 0.1; this->emission[1][2] = 0.1; this->emission[1][3] = 0.1; }; void Baum_Welsh (int iteration) { vector<double> forward_matrix (this->data.size() * this->data[0].size() * this->states.size(),0); this->forward_matrix = forward_matrix; vector<double> backward_matrix(this->data.size() * this->data[0].size() * this->states.size(),0); this->backward_matrix = backward_matrix; vector<double> seq_prob (data.size(), 0); for(int it = 0; it < iteration; it++) { // Berechnung der Forward-Backward-Sequenz-Wkt. this->forward_procedure(); this->backward_procedure(); for(int n = 0; n < this->data.size(); n++) { seq_prob[n] = 0; for(int i = 0; i < this->states.size(); i++) {seq_prob[n] += exp(this->backward_matrix[n * this->states.size() * this->data[0].size() + (this->data[n].size()/2) * this->states.size() + i] + this->forward_matrix[n * this->states.size() * this->data[0].size() + (this->data[n].size()/2) * this->states.size() + i]);} seq_prob[n] = log(seq_prob[n]); } // E-M-Schritt for(int s = 0; s < this->states.size(); s++) { double Sum_Msb = 0; vector<double> Msa(this->alphabet.size(), 0); for(int n = 0; n < this->data.size(); n++) { for(int k = 0; k < this->data[n].size(); k++) { for(int b = 0; b < this->alphabet.size(); b++) { if( this->data[n][k] == this->alphabet[b][0] ) { double tmp = exp(this->forward_matrix[n * this->states.size() * this->data[0].size() + k * this->states.size() + s] + this->backward_matrix[n * this->states.size() * this->data[0].size() + k * this->states.size() + s] - seq_prob[n]); Sum_Msb += tmp; Msa[b] += tmp; } } } } double Sum_Nsv = 0; vector<double> Nsr(this->states.size(), 0); for(int n = 0; n < this->data.size(); n++) { for(int j = 1; j < this->data[n].size(); j++) { for(int v = 0; v < this->states.size(); v++) { // eigentlich emission[s] statt emission[v], aber R machst es mit v double tmp = exp(this->forward_matrix[n * this->states.size() * this->data[0].size() + (j-1)* this->states.size() + s] + this->backward_matrix[n * this->states.size() * this->data[0].size() + j * this->states.size() + v] + log(this->transition[s * this->states.size() + v])+log(this->emission[v][this->my_hash(data[n][j], this->alphabet)])-seq_prob[n]); Sum_Nsv += tmp; Nsr[v] += tmp; } } } for(int a = 0; a < this->alphabet.size(); a++) { this->emission[s][a] = Msa[a] / Sum_Msb; } for(int r = 0; r < this->states.size(); r++) { this->transition[s * this->states.size() + r] = Nsr[r] / Sum_Nsv; } } } } }; int main(int argc, char *argv[]) { HMM *tmp_HMM = new HMM(argc, argv); tmp_HMM->Baum_Welsh(100); return 0; }
-
Achso bevor die Frage kommt. Meine Eingabe mit realistischen Werten ist natürlich wesentlich größer. Und hier wird die Laufzeit erst unerträglich.
Meine eigentlichen Daten:
ca. 1000 DNA Sequenzen der Länge 500 mit etwa 30 Zuständen (und Alphabetgröße 4).
-
Dieses Programm hat eine Ausführungszeit von ungefähr 0. Wieviel besser soll das noch werden?
Außerdem hat das Programm jede Menge Speicherlöcher und du selber scheinst an Includitis zu leiden.
P.S.: Und du hast nicht einmal grundlegende Tipps umgesetzt, die man die schon vor Wochen gegeben hat.
-
Geht doch.
(Auch wenns schrecklich formatiert ist..)
Aber wie lange soll die Ausführung denn dauern für 100, 1000, 10000? Also wie lange dauert sie bei dir? Weil 100 ist bei mir in 0.05 Sekunden fertig..
-
Ein Testdatensatz, bei dem das Programm nicht nach dem Bruchteil einer Sekunde fertig ist, wäre noch schön.
-
...hatte mich vertippt...
-
ComputerCarl schrieb:
Achso bevor die Frage kommt. Meine Eingabe mit realistischen Werten ist natürlich wesentlich größer. Und hier wird die Laufzeit erst unerträglich.
Meine eigentlichen Daten:
ca. 1000 DNA Sequenzen der Länge 500 mit etwa 30 Zuständen (und Alphabetgröße 4).Bitte mit den Daten. Die Ausführung sollte etwa 4 - 8 Sekunden dauern, um ordentlich testen zu können.
(Und alle Verhältnisse von Daten <-> Iterationen etc. sollten möglichst dem Originaldatensatz entsprechen.)
-
Mir war so das ich die Tipps umgesetzt habe. Nur in Emission ist noch ein 2-dim Vector, aber den habe ich einfach noch nicht geschaft zu ändern. Was sind diese Speicherlöcher? Das sagt mir gerade nichts.
Wie lang dieses Programm bei mir dauert weiß ich nichtt genau. Auch nicht mehr als eine Sekunde oder so. Nur wie gesagt bei großen Daten wird es schlimm. Einen kleinen Tippfehler hatte ich noch hier:this->emission[1][0] = 0.1;
this->emission[1][1] = 0.1;
this->emission[1][2] = 0.1;
this->emission[1][3] = 0.1; // <- 0.7 !!! sind ja WahrscheinlichkeitenSonst mit diesen Strings müsste es (jedenfalls bei mir) mehr als 4 Sekunden dauern.
-
Das ist dann gut 0.9 Sekunden bei mir. Gibt es eine Möglichkeit, das Ergebnis zu überprüfen? Muss man dir eigentlich alles so aus der Nase ziehen? Deine gesamte Einstellung deutet darauf hin, dass du keine Hilfe möchtest. Sieben Seiten lang. Warum sollte man dir noch helfen?
-
Neue Strings + neue Wahrscheinlichkeit -> Laufzeit 1.33 Sekunden. Na ja zumindest kann man da irgendwie was messen. Bitte die Strings hier wegeditieren und auf pastebin.com oder einer ähnlichen Seite Posten und dann verlinken. (Die machen das arme Forum kaputt. :))
-
Sorry für Forum kaputt machen. Sind schon wieder weg.
"Deine gesamte Einstellung deutet darauf hin, dass du keine Hilfe möchtest" ? Ähm. Das sehe ich anders. Ich weiss auch nicht wie du darauf kommst. Aber wenn es dich so stört, dann überlies dieses Thema einfach. Ich hatte nicht vor jemand zur Last zu fallen. Trotzdem Danke für dein Bemühen.
Meine Ergebnisse aus diesem Programm habe ich mit einer vorimplementierten R-Klasse überprüft. Diese kann aber nur für eine Sequenz Messungen durchführen. Da ich aber identische Likelihoods wie mein Übungsleiter habe, denke ich dass das Programm korrekt arbeiten. Wie man es anders überprüfen kann, weiß aus dem stehgreif nicht.
-
Ich habe gerade die Laufzeit etwas mehr als halbiert, dadurch dass ich "int my_hash(char a, vector<string> alphabet)" zu "int my_hash(char a, const vector<string>& alphabet)" geändert habe..
Mal gucken was da noch so kommt.
Edit:
Ist dir eigentlich bewusst, dass du i zwei mal verwendest und das eine damit verdeckst? Zudem berechnest du einige Sachen gleich zwei mal direkt hintereinander, da hattest du gleich doppelt Glück, dass der Compiler das gerafft hat.
Edit2:
Na ja, der Profiler mag gerade vor allem log und exp nicht mehr, so wie es sein soll. Wahrscheinlich kann man noch einiges rausholen, aber bei dem Durcheinander habe ich da nicht wirklich mehr Lust drauf. (Bitte rücke Code vernünftig ein! Und lass Zeilen nicht länger als ~120 Zeichen lang werden! Nutze mal ein paar Leerzeilen!)(Nicht gerade optimal, aber bei 100 Threads geht das noch) parallelisieren kannste so:
#include <thread> void Baum_Welsh (int iteration) { // ... vector<thread> threads; for (int i = 0; i != iteration; ++i) { threads.push_back(std::thread([&] (int it) -> void { // ... }, i)); } for (auto i = threads.begin(); i != threads.end(); ++i) i->join();
-
cooky451 schrieb:
Ich habe gerade die Laufzeit etwas mehr als halbiert, dadurch dass ich "int my_hash(char a, vector<string> alphabet)" zu "int my_hash(char a, const vector<string>& alphabet)" geändert habe..
Der Parameter wird sogar gar nicht gebraucht, weil my_hash eh Zugriff aufs Alphabet hat. Warum alphabet ein vector von std::string ist, ist mir auch noch nicht klar. Ein einfacher std::string tuts auch.
-
ComputerCarl schrieb:
Meine Ergebnisse aus diesem Programm habe ich mit einer vorimplementierten R-Klasse überprüft. Diese kann aber nur für eine Sequenz Messungen durchführen. Da ich aber identische Likelihoods wie mein Übungsleiter habe, denke ich dass das Programm korrekt arbeiten. Wie man es anders überprüfen kann, weiß aus dem stehgreif nicht.
Du gibst mit deinem Beispielprogramm etwas aus, was vom Ergebnis abhangt oder das Ergebnis ist. Dann können wir nämlich am Programm herumfummeln und sehen, ob immer noch das gleiche herauskommt. dann hat man viel mehr Freiheit im aufräumen deines Codes.
Aber vergiss es einfach. Ich habe schon lange keinen Bock mehr und nur noch cooky schont interessiert. Und der scheint das Ergebnis nicht zu brauchen, sonst hätte er gefragt.
-
i verwende ich nacheinander in 2 verschiedenen Schleifen. Die kommen direkt nacheinander und sehen wohl nur wegen der Formatierung so aus als ob diese ineinander greifen. Oder hab ich mich jetzt geirrt?
Alphabet braucht ich in der Übergabe wirklich nicht. Das ist noch ein Relikt aus einer früheren Version. Sorry dafür. Aber ein einfacher String würde es wohl nicht tun, da ich eigentlich noch das Alphabet variabel setzen wollte. Also unter Umständen vielleicht ein Protein-Alphabet bei denen eine 3-Buchstaben-Codierung vorkommt. Nötig ist das nicht, könnte aber vielleicht von Vorteil werden. Dafür müsste natürlich noch mehr geändert werden, aber alles muss ja erstmal im Groben laufen...
-
SeppJ schrieb:
und nur noch cooky schont interessiert. Und der scheint das Ergebnis nicht zu brauchen, sonst hätte er gefragt.
Ich hab' halt den Code kopiert und das weggemacht, was der Profiler mir gesagt hat. Viel mehr kann man in dem Code eh nicht kopfschmerzfrei machen.
