Help me translate this easy Phyton code to C, please

[John_Cro]

Hello, people

I only registered here for this matter (although it is a great forum).
I am a ph.d. student and need to investigate eigenvector centrality measure for my homework.
I found Phyton code on net and would like to see does it produce the same results as mine, written in Matlab.
But I don't use Phyton, would someone be so nice and translate it, it is so painfully simple and wright now I don't have time to start learning it from scratch.

Thank you, really

Here is the code:

from numpy import *

def norm2(v):
return sqrt(v.T*v).item()

def PowerMethod(A, y, e):
while True:
v = y/norm2(y)
y = A*v
t = dot(v.T,y).item()
if norm2(y – t*v) <= e*abs(t):
return (t, v)

[Gribouillis]

I guess it could be something like this in C

#include <math.h>

double dot(int sz, double *u, double* v){
    int i;
    double res = 0.0;
    for(i = 0; i < sz; ++i)
        res += u[i] * v[i];
    return res;
}

void cumul(int sz, double *u, double* res){
    int i;
    for(i = 0; i < sz; ++i)
        res[i] += u[i];
}

double norm2(int sz, double *u){
    return sqrt(dot(sz, u, u));
}

void prod(int sz, double **A, double *v, double *res){
    int i;
    for(i = 0; i < sz; ++i)
        res[i] = dot(sz, A[i], v);
}

void zoom(int sz, double *y, double fact, double *res){
    int i;
    for(i = 0; i < sz; ++i){
        res[i] = y[i] * fact;
    }
}

void PowerMethod(int sz, double **A, double* y, double e, int *t, double *v){
    int stop = 0;
    while(1){
        zoom(sz, y, 1.0/norm2(sz, y), v);
        prod(sz, A, v, y);
        *t = dot(sz, v, y);
        if(stop)
            return;
        zoom(sz, v, -*t, v);
        cumul(sz, y, v);
        if(norm2(sz, v) <= e * abs(*t))
            stop = 1;
    }    
}

I don't have many opportunities to code C, so you can certainly improve it :)

Edit: you still need to write the main function and allocate and initialize a few arrays of double ...

[John_Cro]

Thank you really, very good job, I would have lost days, than you once again for your effort

[Gribouillis]

Thank you really, very good job, I would have lost days, than you once again for your effort

Hm, I already found a bug, at line 34 it should be double *t instead of int *t.

[Gribouillis]

Hm, I already found a bug, at line 34 it should be double *t instead of int *t.

Also abs() should be fabs() http://www.gnu.org/software/libc/manual/html_node/Absolute-Value.html .

[John_Cro]

I only don't get line 40. where it says

if(stop) return

.

It will always stop here as variable stop is being declared before

while(1)

and noone is changing it's value, so on line 40. it will still be 0.
I guess some function shoul check stop variable before this?

Thanks

[Gribouillis]

Stop is changed at line 45. You could add printf() statements to see what happens, also could you post a typical matrix A and vector y, so that I can test the algorithm ?

[John_Cro]

Sure, here you go:
A =
[
0 1 1 0
1 0 0 1
1 1 0 1
0 0 1 0]

y= [1 2 3 4]

It should output
y=[2,5 2 3 1,5]

[John_Cro]

But stop at line 40 is happening before the line 45 so it will never reach 45. I am sorry am such a pain in the ***

[Gribouillis]

The algorithm seems to work well in python

from numpy import *

def norm2(v):
    return sqrt(dot(v.T, v)).item()

def PowerMethod(A, y, e):
    while True:
        v = y/norm2(y)
        y = dot(A, v.T)
        t = dot(v.T,y).item()
        if norm2(y - t*v) <= e*abs(t):
            return (t, v)
A =array([
    0., 1., 1., 0.,
    1., 0., 0., 1.,
    1., 1., 0., 1.,
    0., 0., 1., 0.,])
A = A.reshape((4, 4))
y = array([1, 2, 3, 4])

t, v = PowerMethod(A, y, 1.0e-12)
print "eigenvalue", t
print "eigenvector", v
expected=array([2.5, 2, 3, 1.5])
print "expected   ",expected/norm2(expected)

""" my output -->
eigenvalue 2.0
eigenvector [ 0.53916387  0.43133109  0.64699664  0.32349832]
expected    [ 0.53916387  0.43133109  0.64699664  0.32349832]
"""

I'm trying the C version.

[Gribouillis]

But stop at line 40 is happening before the line 45 so it will never reach 45.

Lines 40 and 45 are in a loop, so line 40 is can still be reached at line 45.

Here is a fully working C version with a main function:

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

double dot(int sz, double u[], double v[]){
    int i;
    double res = 0.0;
    for(i = 0; i < sz; ++i)
        res += u[i] * v[i];
    return res;
}

void cumul(int sz, double u[], double res[]){
    int i;
    for(i = 0; i < sz; ++i)
        res[i] += u[i];
}

double norm2(int sz, double u[]){
    return sqrt(dot(sz, u, u));
}

void prod(int sz, double A[4][4], double v[], double res[]){
    int i;
    for(i = 0; i < sz; ++i)
        res[i] = dot(sz, A[i], v);
}

void zoom(int sz, double y[], double fact, double res[]){
    int i;
    for(i = 0; i < sz; ++i){
        res[i] = y[i] * fact;
    }
}

void PowerMethod(int sz, double A[4][4], double y[], double e, double *t, double v[]){
    int stop = 0;
    while(1){
        zoom(sz, y, 1.0/norm2(sz, y), v);
        prod(sz, A, v, y);
        *t = dot(sz, v, y);
        if(stop)
            return;
        zoom(sz, v, -*t, v);
        cumul(sz, y, v);
        if(norm2(sz, v) <= e * fabs(*t))
            stop = 1;
    }    
}

int main(){
    double A[4][4] = {
        {0., 1., 1., 0.},
        {1., 0., 0., 1.},
        {1., 1., 0., 1.},
        {0., 0., 1., 0.}};
    double y[4] = {1., 2., 3., 4.};
    double v[4] = {0., 0., 0., 0.};
    double t;
    PowerMethod(4, A, y, 1.e-12, &t, v);
    printf("eigenvalue %.7f\n", t);
    printf("eigenvector %.7f %.7f %.7f %.7f\n", v[0], v[1], v[2], v[3]); 
}

/**********
 * my output  -->
 * eigenvalue 2.0000000
 * eigenvector 0.5391639 0.4313311 0.6469966 0.3234983
 * 
 * the compilation command was
 *      gcc -o centrality centrality.c -lm
***********/

Note that my C is terrible, so perhaps you could post this code to the C forum and ask if they can improve the code.