#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "sparse_matrix.h"
#include "matrix_operation.h"
#include "vector.h"

double* power_algorithm_step(const SparseMatrix *matrix, const double *pi) {
    double* result = malloc(matrix->num_nodes*sizeof(double));
    multiply_vector_matrix_parallel(pi, matrix, result);
    return result;
}

double* power_algorithm(const SparseMatrix *matrix, double epsilon) {
    int N = matrix->num_nodes;
    int vecsize = N*sizeof(double);
    double* pi = malloc(vecsize);
    double* pi2 = malloc(vecsize);
    init_vector(pi, N, 1.0/(double)N);
    pi2 = power_algorithm_step(matrix, pi);
    while (diff_norm_vector(pi, pi2, N)>epsilon) {
        printf("step\n");
        memcpy(pi, pi2, vecsize);
        pi2 = power_algorithm_step(matrix, pi);
    }
    return pi2;
}

double* pagerank(const SparseMatrix *matrix, double epsilon, double alpha) {
    int N = matrix->num_nodes;
    size_t vec_size = N * sizeof(double);

    double* pi = malloc(vec_size);
    double* pi_new = malloc(vec_size);
    double* f = malloc(vec_size);
    double right_const = (1.0 - alpha) / N;

    init_vector(pi, N, 1.0 / N);

    generate_f(matrix, f);

    double diff;
    int iter = 0;

    do {
        // 1. pi * M
        double* temp = power_algorithm_step(matrix, pi);

        // 2. alpha/N * (pi * f)
        double right_var = (alpha/(double)N) * vec_product(pi, f, N);

        // 3. alpha*(pi*M) + (right_const+alpha/N * (pi * f))*e
        for (int i = 0; i < N; i++) {
            pi_new[i] = alpha * temp[i] + right_const + right_var;
        }

        // 4. Normalize
        normalize_vector(pi_new, N);

        // 5. Calculate convergence
        diff = diff_norm_vector(pi, pi_new, N);

        // 6. Update for next iteration
        free(pi);
        pi = pi_new;
        pi_new = malloc(vec_size);

        if ((++iter)%1 == 0) {
            printf("Iteration %d: diff = %.16f\n", iter, diff);
        }
        free(temp);
    } while (diff > epsilon);

    free(pi_new);
    free(f);
    return pi;
}