/*
 * 'test-korelation.c' created: Mon Oct 28 16:44:08 2002
 *
 * Copyright 2002 by Lubomir Host 'rajo' <8host AT pauli.fmph.uniba.sk>
 *
 * License: GNU GPL (GNU Public License)
 *
 */

/* $Platon$ */

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>

#include "main.h"
#include "devices.h"
#include "tests.h"
#include "test-korelation.h"


RNG_test_t RNG_test_korel = {
	"Korelation test of generators", 2E9, test_korel
};

#define MAX 5

int test_korel(int *counts)
{
	int retval = RNG_OK;
	int i = 0;
	int j = 0;
	int k = 0;
	Real_t *x_i;			/* remembered serie of random numbers */
	Real_t *R_k;			/* korelation koeficients */
	Real_t *Sum_xi;			/* Sum[i=1...n-k] x_i */
	Real_t *Sum_xi_k;		/* Sum[i=1...n-k] x_(i+k) */
	Real_t *Sum_xi_xik;		/* Sum[i=1...n-k] x_(i+k) * x_i */

	while (rng_devices[j] != NULL) {
		int cur = 1;
		int end = counts[0];
		Real_t random;

		printf("# @@ running test for R.N.G. %s\n", (*rng_devices[j]).desc);
		printf("# ** <count>	<seconds>	<microseconds>	<numbers per second>\n");

		/* memory alocations {{{ */
		x_i = (Real_t *) malloc(MAX * sizeof(Real_t));
		if (x_i == NULL) {
			fprintf(stderr, "Allocation of %d Real_t's x_i[] failed!\n", MAX);
			return RNG_FAIL;
		}
		R_k = (Real_t *) malloc(MAX * sizeof(Real_t));
		if (R_k == NULL) {
			fprintf(stderr, "Allocation of %d Real_t's R_k[] failed!\n", MAX);
			return RNG_FAIL;
		}
		Sum_xi = (Real_t *) malloc(MAX * sizeof(Real_t));
		if (Sum_xi == NULL) {
			fprintf(stderr, "Allocation of %d Real_t's Sum_xi[] failed!\n", MAX);
			return RNG_FAIL;
		}
		Sum_xi_k = (Real_t *) malloc(MAX * sizeof(Real_t));
		if (Sum_xi_k == NULL) {
			fprintf(stderr, "Allocation of %d Real_t's Sum_xi_k[] failed!\n", MAX);
			return RNG_FAIL;
		}
		Sum_xi_xik = (Real_t *) malloc(MAX * sizeof(Real_t));
		if (Sum_xi_xik == NULL) {
			fprintf(stderr, "Allocation of %d Real_t's Sum_xi_xik[] failed!\n", MAX);
			return RNG_FAIL;
		}
		/* }}} */

		/* initialize arrays */
		for (k = 0; k < MAX; k++) {
				x_i[k]         = 0.0;
				R_k[k]         = 0.0;
				Sum_xi[k]      = 0.0;
				Sum_xi_k[k]    = 0.0;
				Sum_xi_xik[k]  = 0.0;
		}
		/* generate first MAX random numbers and remember it */
		for (k = 0; k < MAX; k++)
			x_i[k] = (*rng_devices[j]).read();

		/* create Sum's */
		Sum_xi[MAX - 1]     += x_i[0];
		Sum_xi_k[MAX - 1]   += x_i[MAX - 1];
		Sum_xi_xik[MAX - 1] += x_i[MAX - 1] * x_i[0];
		for (k = MAX - 1; k > 0; k--) {
			int l;
			Sum_xi[k - 1]     += Sum_xi[k]     + x_i[MAX - k];
			Sum_xi_k[k - 1]   += Sum_xi_k[k]   + x_i[k - 1];
			for (l = 0; l < (MAX - k + 1); l++)
				Sum_xi_xik[k - 1] +=  x_i[l] * x_i[l + k - 1];
		}
		printf("#$ k	x_i[k]	Sum_xi[k]	Sum_xi_k[k]	Sum_xi_xik[k]	R_k[k]\n");
		for (k = 0; k < MAX; k++) {
			printf("%d	%Lg	%Lg	%Lg	%Lg	%Lg\n",
					k, x_i[k], Sum_xi[k], Sum_xi_k[k], Sum_xi_xik[k], R_k[k]);
		}
		/* look for counts[i] >= MAX */
		while (counts[i] <= MAX && counts[i] != -1)
			i++;
		if (counts[i] != -1) {
			i--;
			cur = MAX;
		}
		while (counts[i] != -1
				&& end != -1
				&& counts[i] <= RNG_test_korel.max_count
				&& counts[i] <= (*rng_devices[j]).max_count
				&& retval == RNG_OK) {
			end = counts[i];
			while ((cur - 1) != end && end != -1) {
				/* printf("cur = %d end = %d\n", cur, end); */

				random = (*rng_devices[j]).read();

				for (k = 0; k < (MAX - 1); k++)
					x_i[k] = x_i[k + 1];
				x_i[MAX - 1] = random;

				for (k = MAX - 1; k > -1; k--) {
					Sum_xi[k]     += x_i[MAX - k - 1];
					Sum_xi_k[k]   += x_i[MAX - 1];
					Sum_xi_xik[k] += x_i[MAX - 1] * x_i[MAX - k - 1];
				}

				cur++;
			}

			/* here we are looking for parameter:
			 *
			 *            1           n-k
			 *         -------  *  Sum[ (X_i - X_avg) * (X_(i+k) - X_avg) ]
			 *         (n - k)        i=1
			 * R_k =  ------------------------------------------------------   =
			 *            1           n
			 *         -------  *  Sum[ (X_i - X_avg) * (X_i - X_avg) ]
			 *            n           i=1
			 *
			 *
			 *     n     Sum_xi_xik[k] - x_avg*(Sum_xi[k] + Sum_xi_k[k]) - (n-k) * x_avg^2
			 * = ----- * -----------------------------------------------------------------
			 *   (n-k)                Sum_xi_xik[0]   -   n * x_avg * x_avg
			 */
			printf("## next count\n%d", cur);
			for (k = 0; k < MAX; k++) {
				Real_t x_avg;

				x_avg = Sum_xi[0] / cur;

				R_k[k] = (cur / (cur - k)) *
					(Sum_xi_xik[k]
					 - x_avg * (Sum_xi[k] + Sum_xi_k[k])
					 - (cur - k) * x_avg * x_avg) /
					(Sum_xi_xik[0] - cur * x_avg * x_avg);

				printf("\t%Lg", R_k[k]);
			}

			printf("\n");
			end = counts[++i + 1];
		}
		j++; /* go to the next device ... */


		free(x_i);
		free(R_k);
		free(Sum_xi);
		free(Sum_xi_k);
		free(Sum_xi_xik);
	}

	return retval;
}

/* Modeline for ViM {{{
 * vim: set ts=4:
 * vim600: fdm=marker fdl=0 fdc=3:
 * }}} */