#include "geostat.h"
#include "log.h"
#include "memory.h"
#include <time.h>

/*builds the sampled covariance function*/
/*dimensions are even*/
void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) {
    double* used_ram_t0 = malloc(sizeof(double));
    getVirtualMemUsed(used_ram_t0);
    
    clock_t t = clock();
    
    int i, j, k, maille, n2[3], symmetric;
    double di, dj, dk;

    log_info("RESULT = in progress");

    struct cpustat initial[cores];
    struct cpustat final[cores];

    for (int i = 0; i < cores; i++) {
        get_stats(&initial[i], i - 1);
    }

    for (i = 0; i < 3; i++)
        n2[i] = n[i] / 2;

    for (i = 0; i <= n2[0]; i++) {
        for (j = 0; j <= n2[1]; j++) {
            for (k = 0; k <= n2[2]; k++) {

                /*area 1*/
                maille = 1 + i + n[0] * (j + n[1] * k);
                di = (double)i * mesh.DX;
                dj = (double)j * mesh.DY;
                dk = (double)k * mesh.DZ;
                covar[maille] = (double)cov_value(variogram, di, dj, dk, cores);

                if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) {
                    /*area 2*/
                    symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k));
                    covar[symmetric] = covar[maille];
                }

                if (i > 0 && i < n2[0]) {
                    /*area 4*/
                    di = -(double)i * mesh.DX;
                    dj = (double)j * mesh.DY;
                    dk = (double)k * mesh.DZ;
                    maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k);
                    covar[maille] = (double)cov_value(variogram, di, dj, dk, cores);
                }

                if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) {
                    /*area 8*/
                    symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k));
                    covar[symmetric] = covar[maille];
                }

                if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) {
                    /*area 5*/
                    di = -(double)i * mesh.DX;
                    dj = -(double)j * mesh.DY;
                    dk = (double)k * mesh.DZ;
                    maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k);
                    covar[maille] = (double)cov_value(variogram, di, dj, dk, cores);
                }

                if (k > 0 && k < n2[2]) {
                    /*area 6*/
                    symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k));
                    covar[symmetric] = covar[maille];
                }

                if (j > 0 && j < n2[1]) {
                    /*area 3*/
                    di = (double)i * mesh.DX;
                    dj = -(double)j * mesh.DY;
                    dk = (double)k * mesh.DZ;
                    maille = 1 + i + n[0] * (n[1] - j + n[1] * k);
                    covar[maille] = (double)cov_value(variogram, di, dj, dk, cores);
                }

                if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) {
                    /*area 7*/
                    symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k));
                    covar[symmetric] = covar[maille];
                }
            }
        }
    }

    t = clock() - t;
    double time_taken = ((double)t)/CLOCKS_PER_SEC; // calculate the elapsed time

    for (int i = 0; i < cores; i++) {
        get_stats(&final[i], i - 1);
    }

    for (int i = 0; i < cores; i++) {
        log_info("CPU %d: %lf%%", i, calculate_load(&initial[i], &final[i]));
    }

    double* used_ram_tf = malloc(sizeof(double));
    getVirtualMemUsed(used_ram_tf);

    log_info("RESULT = success, di = %f, dj = %f, dk = %f, ELAPSED = %f seconds, DIF USED VIRTUAL MEM = %5.1f MB", di, dj, dk, time_taken, *used_ram_tf - *used_ram_t0);

    free(used_ram_t0);
    free(used_ram_tf);
}