[gelöst Danke]Parallelisierung einer For-Schleife

[nurf]

Computer Carl schrieb:

Also warum soll den noch Code gepostet werden, wenn keine Verbindung besteht?

Damit wir damit ein Schweine-Geld verdienen können. Wird dann direkt an SAP, Oracle und Microsoft verkauft und fließt direkt in das aktuelle Produkt ein.
Einmal kurz das Forum durchstöbern, Code abgrasen und schon wandert das Geld aufs Konto. Einfacher kann man einfach kein Geld verdienen.

[otze]

Computer Carl schrieb:

Also warum soll den noch Code gepostet werden, wenn keine Verbindung besteht?

Vielleicht, aber auch nur vielleicht, weil die Leute, die sich professionell mit hochoptimiertem numerischen Code auseinander setzen, sich bereit erklärt haben, das für dich zu profilen und genau zu sehen, wo das Problem ist? Um dir dann mit ihrem Wissen helfen zu können, das zu verbessern?

ich geb dir mal ein bekanntes Beispiel: die naive Matrix-Matrix-Multiplikation kann ohne Parallelisierung und irgendwelche algorithmischen Tricks um Faktor 40 beschleunigt werden. Mit parallelisierung holst du im Vergleich < Faktor 4.

Wenn ich auf deinen Code schaue, sehe ich zum Beispiel ganz spontan sowas:

for(int k = 0; k < states.size(); k++){ 
    tmp_prob +=  forward_matrix[...+ k] * transition[k][j];
}

ist states.size() groß, kann das richtig richtig langsam sein, weil du für jeden lookup einen Cachemiss kriegst (Faktor 5-40 langsamer). Da ichs net profilen kann, kann ich jetzt auch nicht überprüfen, ob man durch Schleifentausch noch was reissen kann.

[Athar]

Wenn du einfach nur einen Weg suchst, etwas parallel auszuführen, dann ist #pragma omp parallel for vor einer for-Schleife einer der simpelsten (dann mit -fopenmp kompilieren und linken). Es gibt eine Reihe von Regeln, die du dabei einhalten musst, insbesondere müssen die Iterationen unabhängig voneinander sein (gleichzeitig ausführbar und in beliebiger Reihenfolge). So etwas wie tmp_prob += ... ginge so schon einmal nicht, allerdings kann dir OpenMP auch da helfen (reduction). Jedoch muss die Anzahl der Iterationen groß sein, damit sich das lohnt ("groß" ist allerdings auch wieder relativ).

[otze]

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.

[Michael E.]

otze schrieb:

Kann mich jemand bestätigen?

Das hier: [...]
sieht mir sehr wie eine Matrix-Matrix multiplikation aus.

Jo. Passt zum Algorithmus.

[dimiKL]

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.

[cooky451]

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.)

[ComputerCarl]

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.

[otze]

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.

[cooky451]

Du kannst die Dateien ruhig mit hochladen. Es sollten halt nur nicht 5 .hpp/cpp Paare sein.

[dimiKL]

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.

[otze]

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.

[Sone]

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 inline impliziert, falls du dich daran erinnern kannst- heute machen es die meisten Compiler automatisch, das Schlüsselwort wird nur noch in einem anderen Kontext verwendet).

[lecker]

ha! du bist so fail das du schon wieder epic bist!

[Computer Carl]

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;
}

[ComputerCarl]

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).

[SeppJ]

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.

[cooky451]

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..

[Michael E.]

Ein Testdatensatz, bei dem das Programm nicht nach dem Bruchteil einer Sekunde fertig ist, wäre noch schön.

[ComputerCarl]

...hatte mich vertippt...