Files
argyllcms/numlib/soboltest.c
2026-08-20 20:28:40 +01:00

128 lines
4.1 KiB
C

/* 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;
}
}