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

double *calcHistograms (double ***ai, int n1, int n2, int n3)
{
	int i, j, k, l;

	double const d0 = 1;
	double const dx = 2.2;
	double w = d0 * 0.7;

	double *histogram = (double *) malloc(sizeof(double) * n1 * n2 * BB);
	double *x = (double *) malloc(sizeof(double) * n3);
	double v[BB];

	/* Initialize v */
	for (i = 0; i < BB; i++) {
		v[i] = dx - d0 * i;
		v[i] /= w;
	}

    double sumOfSquares = 0.0;

	for (i = 0; i < n1; i++) {
		double r = pow(2, (i+1)/2.0);
		for (j = 0; j < n2; j++) {
			int lengthX = 0;
			for (k = 0; k < n3; k++) {
				if (ai[i][j][k] > 0) {
					x[lengthX++] = log(ai[i][j][k]);
				}
			}

			for (k = 0; k < BB; k++) {
				double sum = 0;
				for (l = 0; l < lengthX; l++) {
					double xk = 1 - (x[l]/w - v[k]) * (x[l]/w - v[k]);
					if (xk > 0) {
						sum += xk * r;
					}
				}

				// To match Gwen's reshaping, we have to unpack a 3-dimensional histogram into
				// a vector. Note that, in Matlab, the array indexes start with *least* significant index,
				// and become more significant to the right.
                sumOfSquares += sum * sum;
				histogram[k * (n1 * n2) + j * n1 + i] = sum;
			}
		}
	}

    sumOfSquares = sqrt( sumOfSquares );
//    printf( "calcHistograms sumOfSquares: %f\n", sumOfSquares );
    int cnt = n1 * n2 * BB;
    int idx;
    for( idx = 0; idx < cnt; ++idx )
        histogram[idx] /= sumOfSquares;

	free(x);

	return histogram;
}

/* Test case -- outputs Matlab code */
void test_main()
{
	int const n1 = 8, n2 = 5, n3 = 500;
	int i, j, k;
	double ***ai;

	/* Initialize 3-D array */
	ai = (double ***) malloc(sizeof(double **) * n1);
	for (i = 0; i < n1; i++) {
		ai[i] = (double **) malloc(sizeof(double *) * n2);
		for (j = 0; j < n2; j++) {
			ai[i][j] = (double *) malloc(sizeof(double) * n3);
			for (k = 0; k < n3; k++) {
				double value = (double) rand() / RAND_MAX;
				ai[i][j][k] = value;
				printf("ai(%d,%d,%d) = %.9f;\n", i+1, j+1, k+1, value);
			}
		}
	}

	double *histogram = calcHistograms(ai, n1, n2, n3);
	printf("cHistograms = [ \n");
	for (i = 0; i < BB * n1 * n2; i++) {
		printf("%.9f ", histogram[i]);
	}
	printf("\n];\n");

	/* Free 3-D array */
	for (i = 0; i < n1; i++) {
		ai[i] = (double **) malloc(sizeof(double *) * n2);
		for (j = 0; j < n2; j++) {
			free(ai[i][j]);
		}
		free(ai[i]);
	}
	free(ai);

	free(histogram);
}

/*main ()
{
	int const n1 = 12, n2 = 23, n3 = 1304;
	int i, j, k;
	double ***ai;

	ai = (double ***) malloc(sizeof(double **) * n1);
	FILE *fp = fopen("ai2.dat", "rt");
	for (i = 0; i < n1; i++) {
		ai[i] = (double **) malloc(sizeof(double *) * n2);
		for (j = 0; j < n2; j++) {
			ai[i][j] = (double *) malloc(sizeof(double) * n3);
			for (k = 0; k < n3; k++) {
				double value;
				fscanf(fp, "%lf", &value);
				//printf("val=%.9f\n", value);
				ai[i][j][k] = value;
			}
		}
	}
	fclose(fp);

	double *histogram = calcHistograms(ai, n1, n2, n3);
	printf("cHistograms = [ \n");
	for (i = 0; i < BB * n1 * n2; i++) {
		printf("%.9f ", histogram[i]);
	}
	printf("\n];\n");

	for (i = 0; i < n1; i++) {
		ai[i] = (double **) malloc(sizeof(double *) * n2);
		for (j = 0; j < n2; j++) {
			free(ai[i][j]);
		}
		free(ai[i]);
	}
	free(ai);

	free(histogram);
}*/