[gelöst Danke]Parallelisierung einer For-Schleife
-
Kann mich jemand bestätigen?
Das hier:
for(int i = 1; i < data[n].size(); i++) { for(int j = 0; j < states.size(); j++) { double tmp_prob = 0; for(int k = 0; k < states.size(); k++) { tmp_prob += forward_matrix[... + (i-1) * states.size() + k] * transition[k][j]; } } }sieht mir sehr wie eine Matrix-Matrix multiplikation aus. Der Code ist leider etwas undurchsichtig, weswegen ich das nicht 100% validieren kann. Aber in dem Fall müsste zum beispiel ATLAS oder Eigen an der Stelle bis zu Faktor 40 raus kriegen.
-
otze schrieb:
Kann mich jemand bestätigen?
Das hier: [...]
sieht mir sehr wie eine Matrix-Matrix multiplikation aus.Jo. Passt zum Algorithmus.
-
Hallo,
viele schreiben, man sollte an den Optionen für den Compiler etwas machen und hoffen,dass er dann gut optimiert. Einfache Optimierungen kann man aber auch selbst machen.
- Methoden-/Funktionsaufrufe vermeiden,
- Teil-Offsets berechnen und verwenden
- 1-dim. Felder verwenden. (Mehrdim. Felder nur bei wirklich großen Datenmenge, mehrere GB)Zum Beispiel:
int numStates = (int)states.size(); int sizeData = (int)data.size(); int statesMulLength = numStates * max_length; int matrixSize = statesMulLength * sizeData; double* forward_matrix = new double[matrixSize]; double* backward_matrix = new double[matrixSize]; double* seq_prob = new double[sizeData]; for(int i=0; i < matrixSize; i++) { forward_matrix[i] = 0; } for(int i=0; i < matrixSize; i++) { backward_matrix[i] = 0; } for(int i=0; i < sizeData; i++) { seq_prob[i] = 0; } //int matIdx; int matRowOffset; int matRowOffsetLoopConst; int data_nsize; for(int n = 0; n < sizeData; n++) { // Berechnung aller Forward-Matrizen matRowOffsetLoopConst = statesMulLength*n; data_nsize = data[n].size(); int hash_value = my_hash(data[n][0], alphabet); for(int i = 0; i < numStates; i++) { forward_matrix[matRowOffsetLoopConst + i] = start[i] * emission[i][hash_value]; } matRowOffset2 = matRowOffsetLoopConst + numStates; for(int i = 1; i < data_nsize; i++) { hash_value = my_hash( data[n][i], alphabet); for(int j = 0; j < numStates; j++) { double tmp_prob = 0; for(int k = 0; k < numStates; k++) { tmp_prob += forward_matrix[matRowOffset + k] * transition[k][j]; } forward_matrix[matRowOffset2+ j] = tmp_prob * emission[j][hash_value]; } matRowOffset += numStates; matRowOffset2 += numStates; } // Berechnung aller Backward-Matrizen matRowOffset = matRowOffsetLoopConst + (data_nsize-1)*numStates; for(int i = 0; i < numStates; i++) { backward_matrix[matRowOffset + i] = 1; } matRowOffset = matRowOffsetLoopConst; int matRowOffset2 = matRowOffsetLoopConst + numStates*(data_nsize-1); for(int i = data_nsize-2; i >= 0; i--) { hash_value = my_hash(data[n][i+1],alphabet); for(int j = 0; j < numStates; j++) { double tmp_prob = 0; for(int k = 0; k < numStates; k++) { tmp_prob += backward_matrix[matRowOffset2 + k] * transition[j][k] * emission[k][hash_value]; } backward_matrix[matRowOffset + j] = tmp_prob; } matRowOffset += numStates; matRowOffset2 -= numStates; } // Berechnung der Sequenz-Wahrscheinlichkeiten matRowOffset = matRowOffsetLoopConst + 1*numStates; double sp= 0; for(int i = 0; i < numStates; i++) { sp += backward_matrix[matRowOffset + i] * forward_matrix[matRowOffset + i]; } seq_prob[n] = sp; }Für parallele Implementierung, falls es hier möglich. Empfehle ich OpenCL oder CUDA (für NVIDIA Grafikkarten). OpenCL funktioniert nicht nur mit GPUs, sondern auch mit CPUs, wenn die entsprechenden Libs vorhanden/installiert sind.
PS: Der Code oben kann Fehler enthalten. Hoffe die Idee ist klar.
-
dimiKL schrieb:
- Methoden-/Funktionsaufrufe vermeiden,
Fail.
dimiKL schrieb:
- Teil-Offsets berechnen und verwenden
Fail.
dimiKL schrieb:
- 1-dim. Felder verwenden.
Hey, der Tipp ist brauchbar! (Nur halt schon gefühlte 20 mal geschrieben worden in diesem Thread.)
-
Oh ich hatte schon nicht damit gerechnet, dass noch jemand etwas dazu schreibt. Ich habe die letzten Tage versucht ein bisschen mit diesen WINAPI threads zu arbeiten. Hat noch ein bisschen Zeit rausgehold, da ich Forward- und Backward- Berechnung trennen kann.
Ich werde mal schauen, ob ich es die nächsten Tage schaffe den Code hier zu posten. Ich müsste viel, was aus einer Datei ausgelesen wird, in den Code schreiben und darauf hatte ich eigentlich nicht soviel Lust. Wenn das Angebot aber noch steht, würde ich die Tage gerne darauf zurückgreifen.
-
naja, für mich ist die Sache klar: schnapp dir eine ordentliche Bibliothek für matrix-Matrixmultiplikation. Damit sollten sich alle Performancesorgen in Luft auflösen. Und die meisten Bibliotheken bieten außerdem noch Unterstützung für mehrere Kerne an, also musst du dich nichtmal mit threads rumschlagen.
-
Du kannst die Dateien ruhig mit hochladen. Es sollten halt nur nicht 5 .hpp/cpp Paare sein.
-
Hallo,
@cooky451: Bin ziemlich neu in diesem Forum. Was bedeutet "Fail" ?
- Methoden-/Funktionsaufrufe vermeiden: Der Compiler kann normalerweise davon nicht ausgehen, dass bspw. "states.size()" immer den gleichen Wert zurückliefert, also wird immer die Methode immer aufgerufen
- Teil-Offsets berechnen und verwenden: Auch hier kann der Compiler nicht davon ausgehen dass bspw. "states.size()" immer den gleichen Wert zurückliefert. Also wird nochmal Methode aufgerufen, multipliziert und (viel) adddiert.
-
genau das kann der Compiler in vielen Fällen sehr wohl. Das magische Wort in dem Zusammenhang heißt inlining. Die Funktion ist für den Compiler keine Blackbox, sondern er kann rein schauen und sehen was passiert. Und wenn er sieht, dass size() nur eine Variable zurück gibt und die Variable nirgendwo in der Schleife verändert wird - dann wird er sie aus der Schleife entfernen.
-
dimiKL schrieb:
Hallo,
@cooky451: Bin ziemlich neu in diesem Forum. Was bedeutet "Fail" ?
FAIL!!! ROFL!!
dimiKL schrieb:
- Methoden-/Funktionsaufrufe vermeiden: Der Compiler kann normalerweise davon nicht ausgehen, dass bspw. "states.size()" immer den gleichen Wert zurückliefert, also wird immer die Methode immer aufgerufen
Ja natürlich! Bist n ganz schlauer, was? Aber zum glück können die geinlined werden...
(Das wurde früher durch das Schlüsselwort inlineimpliziert, falls du dich daran erinnern kannst- heute machen es die meisten Compiler automatisch, das Schlüsselwort wird nur noch in einem anderen Kontext verwendet).
-
ha! du bist so fail das du schon wieder epic bist!
-
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?
