Abstand von Vektoren berechnen

Hallo Leute,

ich habe zur Zeit ein Problem, wo ich bereits seit ein paar Stunden nicht richtig weiterkomme.
Es geht um folgendes: Ich habe einen Zeilenvektor 3. Dimension gegeben mit tr(v)=1, sowie eine fixe stochastische Matrix, nennen wir sie M. Nun möchte ich die Verteilungen berechnen, das geht so: v*M = v'. Darauf folgend soll der Abstand der Vektoren v zu v' berechnet werden, Stichwort ist da die Euklidsche Metrik, kurz: d(v,v')=sqrt( (v-v')^2 ).
Ich poste einfach mal meinen Code an dieser Stelle, und erläutere im Nachhinein:

double mult_matrix_vector(double *matrix, double *vector, int i) {
    double tmp;
    tmp=vector[0]*matrix[i]  +  vector[1]*matrix[i+3]  +  vector[2]*matrix[i+6];
    return tmp;
}

Diese Routine soll den neuen Eintrag an der Stelle i berechnen, konkreter: v=(i=0 i=1 i=2)*M=(i=0 i=1 i=2)=v' Dabei ist die Matrix M zeilenweise in einem Array gespeichert.

double d_euklid (double* vector_old, double* vector_new) {
    int i; double abstand; double array[3]={0};

    for(i=0;i<3;i++) {
        array[i] = pow((vector_old[i]-vector_new[i]),2);    //      (x-y)^2       //
        array[i] = (int) (array[i]*100000);                    // Hier wird gerundet //
        array[i] = array[i]/100000;                            // Hier wird gerundet //
        }
    abstand = sqrt(array[0]+array[1]+array[2]);             //      ||x-y||       //
    abstand = (int) (abstand*1000000);                         // Hier wird gerundet //
    abstand = abstand/1000000;                                 // Hier wird gerundet //
    return abstand;
}

Nicht viel zu erklären. Hier wird der Abstand der Vektoren v zu v' berechnet. Dabei runde ich großzügig, da die Vektoren v,v' Einträge mit einigen Nachkommastellen haben können, und sonst ein underflow auftritt, bzw die Ausgabe nicht mehr schön ist.
Kommen wir zum Automatismus:

void plot (double* matrix, double* vector_new, int n)   {
    double vector_old[3]={0};
    double *array_abstand;
    array_abstand = (double*) malloc(n*sizeof(double));

    FILE *file;
    file = fopen("abstandplot.txt", "w");

    int i,k;
    for(k=0;k<n;k++){
        for(i=0;i<3;i++) {
            vector_old[i]=vector_new[i];
            vector_new[i]=mult_matrix_vector(matrix,vector_new,i);
        }
        array_abstand[k]=d_euklid(&vector_old,&vector_new);
            printf("\nVerteilung nach Iteration %d:(%f  %f  %f)\n",k+1, vector_new[0],vector_new[1],vector_new[2]);
            printf ("d(%d,%d)=%g\n",k,k+1,array_abstand[k]);
            fprintf( file, "(%d,%g)\n", k, array_abstand[k]);
    }
    fclose(file);
}

n wird in der main definiert als die Anzahl der Wiederholungen, oder anders gesagt, wie viele von zwei aufeinanderfolgenden Vektoren verglichen werden sollen.
Hier habe ich ein Bild der Ausgabe hochgeladen.http://s14.directupload.net/images/120414/k9zua8tn.png
Was fällt also auf: Offenbar wird v' für k>1 nicht mehr berechnet. Ich kann mir aber nicht erklären, wieso. Vielleicht kann mir da einer von euch Tipps geben oder gar zeigen, wo der Fehler liegt und bedanke mich jetzt schon mal für eure Hilfe.

Grüße,
Yarox

Yarox

Ich bins nochmal, registriert und wollte eine Demo-Version unten posten:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double mult_matrix_vector(double *matrix, double *vector, int i)
{
    double tmp;
    tmp=vector[0]*matrix[i]  +  vector[1]*matrix[i+3]  +  vector[2]*matrix[i+6];
    return tmp;
}

double d_euklid (double* vector_old, double* vector_new)
{
    int i;
    double abstand;
    double array[3]= {0};

    for(i=0; i<3; i++)
    {
        array[i] = pow((vector_old[i]-vector_new[i]),2);    //      (x-y)^2       //
        array[i] = (int) (array[i]*100000);                    // Hier wird gerundet //
        array[i] = array[i]/100000;                            // Hier wird gerundet //
    }
    abstand = sqrt(array[0]+array[1]+array[2]);             //      ||x-y||       //
    abstand = (int) (abstand*1000000);                         // Hier wird gerundet //
    abstand = abstand/1000000;                                 // Hier wird gerundet //
    return abstand;
}

void plot (double* matrix, double* vector_new, int n)
{
    double vector_old[3]= {0};
    double *array_abstand;
    array_abstand = (double*) malloc(n*sizeof(double));
    int i,k;
    for(k=0; k<n; k++)
    {
        for(i=0; i<3; i++)
        {
            vector_old[i]=vector_new[i];
            vector_new[i]=mult_matrix_vector(matrix,vector_new,i);
        }
        array_abstand[k]=d_euklid(&vector_old,&vector_new);
        printf("\nVerteilung nach Iteration %d:(%f  %f  %f)\n",k+1, vector_new[0],vector_new[1],vector_new[2]);
        printf ("d(%d,%d)=%g\n",k,k+1,array_abstand[k]);
    }
}

int main()
{
    double matrix_P[9]= {0.5,0.3333333333333333333,0.166666667,
                         0.75,0,0.25,
                         0,1.,0
                        };
    double startverteilung[3]= {1,0,0};

    plot (matrix_P,startverteilung,5);
    return 0;
}

Vielleicht hilfts

SeppJ

In Zeile 45 gehören die & dort nicht hin.

Der Cast vom malloc gehört dort nicht hin.

Das malloc hat kein free.

Zum Thema Runden: WTF? Wenn du in der Ausgabe runden willst, dann runde in der Ausgabe (guck mal die Formatoptionen von printf an), aber doch nicht in der Rechnung! Wie kommst du auf so eine Idee?

Deine Matrixmultiplikation ist, zumindest nach gängiger Konvention, falsch. Du multiplizierst den Vector mit einer Spalte der Matrix, nicht mit einer Zeile. Du meinst vermutlich

tmp=vector[0]*matrix[3*i]  +  vector[1]*matrix[3*i+1]  +  vector[2]*matrix[3*i+2];

Das sind die direkten Fehler, die mir aufgefallen sind. Andernfalls ist das Programm natürlich unendlich umständlich, aber das ist dir vermutlich bewusst. Mit den Fehlern behoben sieht die Ausgabe so aus:

Verteilung nach Iteration 1:(0.500000  0.375000  0.375000)
d(0,1)=0.728869

Verteilung nach Iteration 2:(0.437500  0.421875  0.421875)
d(1,2)=0.0911086

Verteilung nach Iteration 3:(0.429688  0.427734  0.427734)
d(2,3)=0.0113886

Verteilung nach Iteration 4:(0.428711  0.428467  0.428467)
d(3,4)=0.00142357

Verteilung nach Iteration 5:(0.428589  0.428558  0.428558)

Ist dies das, was du suchtest?

Yarox

Hallo.
Ja, danke dir. Das ist genau das, was ich haben wollte. Der Abstand geht gegen 0, hervorragend Mir ists ja schon ein wenig peinlich, die Matrixmultiplikation falsch implementiert zu haben.
Ich danke dir auch für deine weiteren Hinweise. Ich bin noch Neuling im Programmieren, da ist sowas sehr hilfreich.
-Der Cast von malloc gehört da nicht rein, weil das array bereits ein double ist, oder?
-free - mal wieder- vergessen. Hätte ich aber in valgrind spätestens gemerkt, danke trotzdem.
-Zum Thema runden: Das Problem ist, dass wenn ich

printf ("d(%d,%d)=%.6g\n",k,k+1,array_abstand[k]);

schreibe, die Ausgabe so aussieht:

0.5             0.333           0.167
        0.75            0               0.25
        0               1               0
Zufaellig generierte Startverteilung:(0.803345  0.096380  0.100276)

Bestimme n: 10

Verteilung nach Iteration 1:(0.450511  0.362953  0.362953)
d(0,1)=0.514346

...

Verteilung nach Iteration 5:(0.400490  0.400469  0.400469)
d(4,5)=0.000124647

Verteilung nach Iteration 6:(0.400479  0.400477  0.400477)
d(5,6)=1.55808e-005

...

Verteilung nach Iteration 10:(0.400478  0.400478  0.400478)
d(9,10)=3.84806e-009

Das Programm ist in der Tat ein wenig umständlich, wollte es anfangs etwas anders handhaben, bin aber daran gescheitert, dass eine double foo() kein array zurückgeben kann. Darum hab ich nach einer alternative gesucht. Schlägt sich das denn auf die Laufzeit des Programmes negativ aus, so, wie ich das Programm realisiert habe?

edit:
ok, mit

printf ("d(%d,%d)=%.6f\n",k,k+1,array_abstand[k]);

funktionierts einwandfrei. Somit hat sich das auch geklärt. Nochmals danke dir für deine Hilfe.

SeppJ

Yarox schrieb:

-Der Cast von malloc gehört da nicht rein, weil das array bereits ein double ist, oder?

Nein, der gehört da nicht hin, weil malloc einen void* zurückliefert, der implizit in alle anderen Datenzeiger überführt werden kann. In C++ wäre solch ein Cast nötig (aber dann auch besser ganz anders und malloc benutzt man dort sowieso nicht). Durch den Cast verdeckst du bloß Compilerfehler, falls du mal malloc falsch benutzt haben solltest, hast aber keine Vorteile.

-Zum Thema runden: Das Problem ist, dass wenn ich
[...]
Schlägt sich das denn auf die Laufzeit des Programmes negativ aus, so, wie ich das Programm realisiert habe?

Ein bisschen. Viel schlimmer ist, dass dein Ergebnis nicht stimmt und du machst ja anscheinend irgendetwas numerisches. Man rundet niemals Zwischenergebnisse.

edit:
ok, mit
printf ("d(%d,%d)=%.6f\n",k,k+1,array_abstand[k]);
funktionierts einwandfrei. Somit hat sich das auch geklärt. Nochmals danke dir für deine Hilfe.

So ist's recht.
printf ist sehr mächtig, da bleibt kein (normaler) Wunsch bei der Formatierung offen.

Wutz

Yarox schrieb:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* hier kannst du jeweils noch ein 'const' spendieren, du kaufst du damit kostenlose Prüfungen des Compilers ein und hilfst ihm beim Optimieren */
double mult_matrix_vector(double *matrix, double *vector, int i)
double mult_matrix_vector(const double *matrix, const double *vector, int i)

/* hier ebenso */
double d_euklid (const double* vector_old, const double* vector_new)

/* hier suggerierst du dem Quellcodeleser eine dir unbekannte Genauigkeit */
    double matrix_P[9]= {0.5,0.3333333333333333333,0.166666667,
/* lasse den Compiler für dich rechnen, also besser */
    double matrix_P[9]= {0.5,1.0/3,1.0/6,

SeppJ schrieb:

printf ist sehr mächtig, da bleibt kein (normaler) Wunsch bei der Formatierung offen

Die Unterdrückung der führenden Nullen im Exponenten bei %e (und %g) hätte ich gerne.

SeppJ

Yarox schrieb:

Das Programm ist in der Tat ein wenig umständlich, wollte es anfangs etwas anders handhaben, bin aber daran gescheitert, dass eine double foo() kein array zurückgeben kann.

Das macht man üblicherweise, indem man das Zielarray als Funktionsparameter mit übergibt und dieses wird dann in der Funktion befüllt.