#include <stdio.h>
#include <numeric.h>

int main() {

    /* data consists of ten points in two dimensional space */
    array double data[2][10] = {{2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1},
                                {2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9}};
    array double mymean[2];
    array double mean_comp[2][10];
    array double trans_manual_covar[10][10], manual_covar[2][2];
    array double evalues[2], evectors[2][2];
    array double trans_evalues[10], trans_evectors[10][10];
    int i;

    /* computes mean of the data and put the result in mymean */
    mean(data, mymean);
    printf("mean  =\n%f\n", mymean);

    /* this subtracts the mean from the original data */
    mean_comp = data;
    for (i = 0; i < 10; i++) {
        mean_comp[0][i] = mean_comp[0][i] - mymean[0];
        mean_comp[1][i] = mean_comp[1][i] - mymean[1];
    }

    printf("mean comp  =\n%f\n", mean_comp);

    /* just to show how the covariance routine actually works */
    manual_covar = mean_comp * transpose(mean_comp);
    printf("manual covar =\n%f\n", manual_covar);
    /* now we compute eigenvalues and eigenvectors */
    eigensystem(evalues, evectors, manual_covar);

    //printf("evalues  =\n%f\n", evalues/norm(evalues,"2"));
    printf("evalues  =\n%f\n", evalues);
    //printf("evectors  =\n%f\n", evectors);

  
    /* just to show how the covariance routine actually works */
    trans_manual_covar = transpose(mean_comp) * mean_comp;
    //printf("manual covar =\n%f\n", manual_covar);
    /* now we compute eigenvalues and eigenvectors */
    eigensystem(trans_evalues, trans_evectors, trans_manual_covar);

    //printf("trans evalues  =\n%f\n", trans_evalues/norm(trans_evalues,"2"));
    printf("trans evalues  =\n%f\n", trans_evalues);
    //printf("evectors  =\n%f\n", evectors);
  
        
    return 0;
}