/* Test the sobol sub-random vector sequencer */

#include <stdlib.h>
#include "numlib.h"

/* Reference vectors */
/* First 10 of dimension 40 */
static double refv[10][40] = {
	{
	0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000,
	0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000,
	0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000,
	0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000 
	}, {
	0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.2500, 0.7500,
	0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.7500, 0.2500, 0.7500, 0.2500,
	0.7500, 0.2500, 0.7500, 0.2500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500,
	0.2500, 0.7500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500
	}, {
	0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.7500, 0.2500,
	0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.2500, 0.7500, 0.2500, 0.7500,
	0.2500, 0.7500, 0.2500, 0.7500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500,
	0.7500, 0.2500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500, 0.2500, 0.7500
	}, {
	0.3750, 0.3750, 0.6250, 0.1250, 0.8750, 0.8750, 0.1250, 0.6250, 0.1250, 0.8750,
	0.3750, 0.6250, 0.1250, 0.3750, 0.6250, 0.1250, 0.6250, 0.3750, 0.3750, 0.8750,
	0.8750, 0.6250, 0.1250, 0.8750, 0.1250, 0.8750, 0.8750, 0.1250, 0.6250, 0.6250,
	0.3750, 0.3750, 0.3750, 0.3750, 0.6250, 0.1250, 0.8750, 0.8750, 0.1250, 0.6250
	}, {
	0.8750, 0.8750, 0.1250, 0.6250, 0.3750, 0.3750, 0.6250, 0.1250, 0.6250, 0.3750,
	0.8750, 0.1250, 0.6250, 0.8750, 0.1250, 0.6250, 0.1250, 0.8750, 0.8750, 0.3750,
	0.3750, 0.1250, 0.6250, 0.3750, 0.6250, 0.3750, 0.3750, 0.6250, 0.1250, 0.1250,
	0.8750, 0.8750, 0.8750, 0.8750, 0.1250, 0.6250, 0.3750, 0.3750, 0.6250, 0.1250
	}, {
	0.6250, 0.1250, 0.3750, 0.3750, 0.1250, 0.6250, 0.8750, 0.8750, 0.3750, 0.1250,
	0.1250, 0.3750, 0.3750, 0.6250, 0.8750, 0.8750, 0.3750, 0.1250, 0.6250, 0.6250,
	0.1250, 0.8750, 0.8750, 0.6250, 0.3750, 0.1250, 0.6250, 0.8750, 0.8750, 0.3750,
	0.1250, 0.6250, 0.6250, 0.1250, 0.3750, 0.3750, 0.1250, 0.6250, 0.8750, 0.8750
	}, {
	0.1250, 0.6250, 0.8750, 0.8750, 0.6250, 0.1250, 0.3750, 0.3750, 0.8750, 0.6250,
	0.6250, 0.8750, 0.8750, 0.1250, 0.3750, 0.3750, 0.8750, 0.6250, 0.1250, 0.1250,
	0.6250, 0.3750, 0.3750, 0.1250, 0.8750, 0.6250, 0.1250, 0.3750, 0.3750, 0.8750,
	0.6250, 0.1250, 0.1250, 0.6250, 0.8750, 0.8750, 0.6250, 0.1250, 0.3750, 0.3750 
	}, {
	0.1875, 0.3125, 0.3125, 0.6875, 0.5625, 0.1875, 0.0625, 0.9375, 0.1875, 0.0625,
	0.6875, 0.8125, 0.5625, 0.6875, 0.1875, 0.6875, 0.1875, 0.0625, 0.0625, 0.8125,
	0.9375, 0.3125, 0.5625, 0.3125, 0.4375, 0.4375, 0.6875, 0.4375, 0.8125, 0.5625,
	0.9375, 0.8125, 0.1875, 0.3125, 0.3125, 0.6875, 0.5625, 0.1875, 0.0625, 0.9375
	}, {
	0.6875, 0.8125, 0.8125, 0.1875, 0.0625, 0.6875, 0.5625, 0.4375, 0.6875, 0.5625,
	0.1875, 0.3125, 0.0625, 0.1875, 0.6875, 0.1875, 0.6875, 0.5625, 0.5625, 0.3125,
	0.4375, 0.8125, 0.0625, 0.8125, 0.9375, 0.9375, 0.1875, 0.9375, 0.3125, 0.0625,
	0.4375, 0.3125, 0.6875, 0.8125, 0.8125, 0.1875, 0.0625, 0.6875, 0.5625, 0.4375
	}, {
	0.9375, 0.0625, 0.5625, 0.9375, 0.3125, 0.4375, 0.8125, 0.6875, 0.4375, 0.8125,
	0.9375, 0.0625, 0.8125, 0.4375, 0.4375, 0.4375, 0.9375, 0.3125, 0.8125, 0.5625,
	0.1875, 0.0625, 0.3125, 0.0625, 0.1875, 0.6875, 0.9375, 0.6875, 0.5625, 0.3125,
	0.6875, 0.0625, 0.9375, 0.0625, 0.5625, 0.9375, 0.3125, 0.4375, 0.8125, 0.6875
	}
};

int main() {
	sobol *s;
	int i;
	double vec[40];
	int fail = 0;

	printf("Starting sobol test\n");

	if ((s = new_sobol(40)) == NULL) {
		printf("new_sobol failed\n");
		exit(1);
	}

	for (i = 0; i < 10; i++) {
		int j;
		if (s->next(s, vec)) {
			printf("Next failed\n");
			exit(1);
		}
		printf("Vector %d = ",i);
		for (j = 0; j < 40; j++) {
			printf("%f ",vec[j]);
		}
		printf("\n");

		for (j = 0; j < 40; j++) {
			if (refv[i][j] != vec[j]) {
				printf("Warning: failed at vector %d, dim %d (%f != %f)\n",i,j,refv[i][j],vec[j]);
				fail = 1;
			}
		}
	}

	s->del(s);
	if (fail) {
		printf("Sobol test FAILED\n");
		return 1;
	} else {
		printf("Sobol test done OK\n");
		return 0;
	}
}