From 22f703cab05b7cd368f4de9e03991b7664dc5022 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Mon, 1 Sep 2014 13:56:46 +0200 Subject: Initial import of argyll version 1.5.1-8 --- rspl/smtnd.c | 1171 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1171 insertions(+) create mode 100644 rspl/smtnd.c (limited to 'rspl/smtnd.c') diff --git a/rspl/smtnd.c b/rspl/smtnd.c new file mode 100644 index 0000000..cdef0a2 --- /dev/null +++ b/rspl/smtnd.c @@ -0,0 +1,1171 @@ + +/*****************************************************/ +/* Smoothness factor tuning of RSPL in N Dimensions. */ +/*****************************************************/ + +/* Author: Graeme Gill + * Date: 28/11/2005 + * Derived from cmatch.c + * Copyright 1995 - 2005 Graeme W. Gill + * + * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :- + * see the License.txt file for licencing details. + * + * Test set for tuning smoothness factor for optimal interpolation + * with respect to dimension, number of sample points, and uncertainty + * of the sample points. + */ + + +#undef DEBUG +#undef DETAILED + +#include +#include +#include +#include "rspl.h" +#include "numlib.h" +#include "xicc.h" /* For mpp support */ +#include "plot.h" +#include "rspl_imp.h" +#include "counters.h" /* Counter macros */ + +/* rspl flags */ +#define MXCHPARAMS 8 + +#define PLOTRES 256 + +/* Function being modeled by rspl */ +/* Similar to MPP model */ +typedef struct { + int di; /* Number of dimensions */ + double ip[MXDI][MXCHPARAMS]; /* Input channel parameters */ + double shape[MXDI][1 << MXDI]; /* Channel interaction shape parameters */ + double op[1 << MXDI]; /* Output channel combination parameters */ +} funcp; + + +/* Setup a random function in the given dimensions */ +static void setup_func(funcp *p, int di) { + double mn,mx; + int i, j; + + p->di = di; + + /* Setup random input parameters */ + /* (This is the one that effects smoothness of function the most) */ + for (j = 0; j < di; j++) { + for (mx = 3.0, i = 0; i < MXCHPARAMS; i++, mx *= 0.6) { + p->ip[j][i] = d_rand(-mx, mx); + } + } + + /* Setup random shape parameters */ + for (j = 0; j < di; j++) { + for (i = 0; i < (1 << di); i++) { /* Initially random */ + p->shape[j][i] = d_rand(-1.0, 1.0); + } + } + + /* Setup the random output parameters */ + mn = 2.0; + mx = -1.0; + for (i = 0; i < (1 << di); i++) { /* Initially random */ + p->op[i] = d_rand(0.0, 1.0); + if (p->op[i] < mn) + mn = p->op[i]; + if (p->op[i] > mx) + mx = p->op[i]; + } + for (i = 0; i < (1 << di); i++) { /* Then scale to between 0.0 and 1.0 */ + p->op[i] = (p->op[i] - mn)/(mx - mn); + } +} + +/* Lookup the function value */ +static double lookup_func(funcp *p, double *v) { + int m, k; + int di = p->di; + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + double ov; /* Output value */ + + /* Input curve lookup */ + for (m = 0; m < di; m++) { + tcnv[m] = icxTransFunc(p->ip[m],MXCHPARAMS,v[m]); + tcnv1[m] = 1.0 - tcnv[m]; + } + + for (m = 0; m < di; m++) + ww[m] = 0.0; + + /* Lookup the shape values */ + for (k = 0; k < (1 << di); k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < di; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + for (m = 0; m < di; m++) { + ww[m] += p->shape[m][k & ~(1<= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + tcnv[m] = vv; + tcnv1[m] = 1.0 - vv; + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < (1 << di); k++) { + double vv = p->op[k]; + for (m = 0; m < di; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + return ov; +} + +/* Do one set of tests and return the results */ +static void do_test( + double *trmse, /* RETURN total RMS error */ + double *tmaxe, /* RETURN total maximum error */ + double *tavge, /* RETURN total average error */ + int verb, /* Verbosity */ + int plot, /* Plot graphs */ + int di, /* Dimensions */ + int its, /* Number of function tests */ + int res, /* RSPL grid resolution */ + int ntps, /* Number of sample points */ + double noise, /* Sample point noise volume */ + int unif, /* NZ if uniform rather than standard deistribution noise */ + double smooth, /* Smoothness to test */ + int twopass, /* Two pass flag */ + int extra /* Extra fit flag */ +); + +/* Compute smoothness of function */ +static double do_stest( + int verb, /* Verbosity */ + int di, /* Dimensions */ + int its, /* Number of function tests */ + int res /* RSPL grid resolution */ +); + +/* ---------------------------------------------------------------------- */ +/* Locate minimum of smoothness series result */ + +#define MXMSS 50 /* Maximum smoothness series */ + +/* Return the optimal smoothness value, based on the */ +/* minimum RMS value. */ +static double best(int n, double *rmse, double *smv) { + int i; + rspl *curve; + co *tps = NULL; + int ns = 500; /* Number of samples */ + datai low,high; + int gres[1]; + datai dlow,dhigh; + double avgdev[1]; + double brmse; /* best solution value */ + double blsmv = 0.0; /* best solution location */ + double rv; /* Return value */ + + /* Create interpolated curve */ + if ((curve = new_rspl(RSPL_NOFLAGS,1, 1)) == NULL) + error ("New rspl failed"); + + /* Create the list of sampling points */ + if ((tps = (co *)malloc(n * sizeof(co))) == NULL) + error ("malloc failed"); + + for (i = 0; i < n; i++) { + tps[i].p[0] = log10(smv[i]); + tps[i].v[0] = rmse[i]; + } + + gres[0] = 100; + low[0] = log10(smv[0]); + high[0] = log10(smv[n-1]); + dlow[0] = 0.0; + dhigh[0] = 1.0; + avgdev[0] = 0.0; + + curve->fit_rspl(curve, + 0, /* Non-mon and clip flags */ + tps, /* Test points */ + n, /* Number of test points */ + NULL, NULL, gres, /* Low, high, resolution of grid */ + NULL, NULL, /* Default data scale */ + -0.00001, /* Underlying smoothing */ + avgdev, /* Average deviation */ + NULL); /* iwidth */ + +#ifdef NEVER + /* Check the fit */ + for (i = 0; i < n; i++) { + co tp; + + tp.p[0] = log10(smv[i]); + curve->interp(curve, &tp); + + printf("Point %d at %f, should be %f is %f\n",i,log10(smv[i]),rmse[i],tp.v[0]); + } + +#define TPRES 100 + /* Plot the result */ + { + double xx[TPRES], yy[TPRES]; + + for (i = 0; i < TPRES; i++) { + co tp; + double vi = i/(TPRES-1.0); + + tp.p[0] = log10(smv[0]) + (log10(smv[n-1]) - log10(smv[0])) * vi; + curve->interp(curve, &tp); + xx[i] = tp.p[0]; + yy[i] = tp.v[0]; + } + do_plot(xx,yy,NULL,NULL,TPRES); + } +#endif + + /* Choose a solution */ + brmse = 1e38; + for (i = 0; i < ns ; i++) { + co tp; + double vi; + + vi = i/(ns-1.0); + tp.p[0] = log10(smv[0]) + (log10(smv[n-1]) - log10(smv[0])) * vi; + curve->interp(curve, &tp); + + if (tp.v[0] < brmse) { + blsmv = tp.p[0]; + brmse = tp.v[0]; + } + } + + rv = pow(10.0, blsmv); + return rv; +} + +/* ---------------------------------------------------------------------- */ +/* Test series */ + +/* Explore ideal smoothness change with test point number and noise volume */ +/* If tdi != 0, just do the given dimension */ +/* If tntps != 0, just do the given number of points */ +/* If tnlev != 0, just do the given noise level */ +static void do_series_1(int unif, int tdi, int tntps, int tnlev, int twopass, int extra) { + int verb = 0; + int plot = 0; + int sdi = 1, edi = 4, di; + int its; + int res = 0; + int ntps = 0; + double noise = 0.0; + double smooth = 0.0; + double trmse, tavge, tmaxe; + int m, i, j, k; + + /* Number of trials to do for each dimension */ + int trials[4] = { + 8, + 6, + 4, + 3 + }; + + /* Resolution of grid for each dimension */ + int reses[4][4] = { + { 257, 129, 65, 33 }, + { 128, 65, 33, 17 }, + { 65, 33, 17, 9 }, + { 33, 17, 9, 5 } + }; + + /* Set of sample points to explore */ + int nset[4][20] = { + { + 5, 10, 20, 50, 100, 200, /* di = 1 */ + }, + { + 25, 100, 400, 2500, 10000, 40000, /* di = 2 */ + }, + { + 25, 50, 75, 125, 250, 500, 1000, 2000, 8000, 125000, /* di = 3 */ + }, + { + 50, 100, 200, 450, 625, 900, 1800, 3600, 10000, 160000, 1000000, /* di = 4 */ + } + }; + + /* Set of smoothnesses to explore */ + double smset[4][20] = { + { + -0.0000001, + -0.0000010, + -0.0000100, + -0.0001000, + -0.0010000, + -0.0100000, + -0.1000000, + -1.0000000, + 0.0 + }, + { + -0.0000001, + -0.0000010, + -0.0000100, + -0.0001000, + -0.0010000, + -0.0100000, + -0.1000000, + -1.0000000, + 0.0 + }, + { + -0.0000010, + -0.0000100, + -0.0001000, + -0.0010000, + -0.0100000, + -0.1000000, + -1.0000000, + 0.0 + }, + { + -0.0000100, + -0.0001000, + -0.0010000, + -0.0100000, + -0.1000000, + -1.0000000, + -10.000000, + 0.0 + } + }; + + /* Set of smoothnesses for twopass smoothing */ + double smset2[4][20] = { + { + -0.0000001, + -0.0000010, + -0.0000100, + -0.0001000, + -0.0010000, + -0.0100000, + -0.1000000, + -1.0000000, + 0.0 + }, + { + -0.0000001, + -0.0000010, + -0.0000100, + -0.0001000, + -0.0010000, + -0.0100000, + -0.1000000, + -1.0000000, + 0.0 + }, + { + -0.0000010, + -0.0000100, + -0.0001000, + -0.0010000, + -0.0100000, + -0.1000000, + -1.0000000, + 0.0 + }, + { + -0.0000100, + -0.0001000, + -0.0010000, + -0.0100000, + -0.1000000, + -1.0000000, + -10.000000, + 0.0 + } + }; + + /* Set of noise levels to explore (average deviation * 4) */ + double noiseset[4][20] = { + { + 0.0, /* Perfect data */ + 0.01, /* 1.0 % */ + 0.02, /* 2.0 % */ + 0.05, /* 5.0 % */ + 0.10, /* 10.0 % */ + 0.20, /* 20.0 % */ + -1.0, + }, + { + 0.0, /* Perfect data */ + 0.01, /* 1.0 % */ + 0.02, /* 2.0 % */ + 0.05, /* 5.0 % */ + 0.10, /* 10.0 % */ + 0.20, /* 20.0 % */ + -1.0, + }, + { + 0.0, /* Perfect data */ + 0.01, /* 1.0 % */ + 0.02, /* 2.0 % */ + 0.05, /* 5.0 % */ + 0.10, /* 10.0 % */ + 0.20, /* 20.0 % */ + -1.0, + }, + { + 0.0, /* Perfect data */ + 0.01, /* 1.0 % */ + 0.02, /* 2.0 % */ + 0.03, /* 3.0 % */ + 0.05, /* 5.0 % */ + 0.10, /* 10.0 % */ + 0.20, /* 20.0 % */ + -1.0, + }, + }; + + + printf("Testing effect of underlying smoothness factors\n"); + if (twopass) + printf("Two Pass smoothing\n"); + if (extra) + printf("Extra fitting\n"); + + /* For dimensions */ + if (tdi != 0) + sdi = edi = tdi; + for (di = sdi; di <= edi; di++) { // dimensions + + its = trials[di-1]; + + for (m = 1; m < 2; m++) { // Just 2nd-highest resolution + res = reses[di-1][m]; + + printf("Tests %d\n",its); + printf("Dimensions %d\n",di); + printf("RSPL resolution %d\n",res); + + /* For noise levels */ + for (j = tnlev; j < 20; j++) { // All noise levels + double smv[20]; + double rmse[20]; + double bfit; + + if (tnlev != 0 && j != tnlev) + break; + + noise = noiseset[di-1][j]; + if (noise < 0.0) + break; + printf("\nNoise volume %f%%, average deviation %f%%\n",noise * 100.0, noise * 25.0); + + /* For number of sample points */ + for (i = 0; i < 20; i++) { // All test points + int rpts; + ntps = nset[di-1][i]; + + if (ntps == 0) + break; + + if (tntps != 0 && ntps != tntps) /* Skip any not requested */ + continue; + + /* Make sure at least 100 points are tested */ + rpts = 1 + 100/ntps; + if (rpts > 5) + rpts = 5; + + printf("\nNo. Sample points %d, norm %8.2f, total its %d\n",ntps, pow((double)ntps, 1.0/di),its * rpts); + + /* For smooth factors */ + for (k = 0; k < 20; k++) { // All smoothing levels + if (twopass) + smooth = smset2[di-1][k]; + else + smooth = smset[di-1][k]; + if (smooth == 0.0) + break; + + printf("Underlying smooth %9.7f, ",-smooth); fflush(stdout); + + do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its * rpts, res, ntps, noise, unif,smooth, twopass, extra); + smv[k] = -smooth; + rmse[k] = trmse; + printf("maxerr %f%%, avgerr %f%%, rmserr %f%%\n", + tmaxe * 100.0, tavge * 100.0, trmse * 100.0); + } + + bfit = best(k, rmse, smv); + printf("Best smoothness = %9.7f, log10 = %4.1f\n",bfit,log10(bfit)); + } + } + } + printf("\n"); + } +} + +/* Explore performance of "optimised" smoothness over test point number and noise volume */ +static void do_series_2(int unif, int twopass, int extra) { + int verb = 0; + int plot = 0; + int di = 0; + int its; + int res = 0; + int ntps = 0; + double noise = 0.0; + double smooth = 0.0; + double trmse, tavge, tmaxe; + int i, j, k; + + /* Number of trials to do for each dimension */ + int trials[4] = { + 16, + 12, + 8, + 5 + }; + + + /* Resolution of grid for each dimension */ + int reses[4] = { + 129, + 65, + 33, + 17 + }; + +#ifdef NEVER + /* Set of sample points to explore */ + int nset[4][20] = { + { + 5, 10, 20, 50, 0 + }, + { + 25, 100, 400, 2500, 0 + }, + { + 125, 1000, 8000, 125000, 0 + }, + { + 625, 10000, 160000, 1000000, 0 + } + }; +#else + /* Set of sample points to explore */ + int nset[4][20] = { + { + 5, 10, 20, 50, 0 + }, + { + 25, 100, 400, 2500, 0 + }, + { + 250, 500, 1000, 2000, 0 + }, + { + 450, 900, 1800, 3600, 0 + } + }; +#endif /* NEVER */ + + + /* Set of smoothnesses to explore */ + double smset[5] = { + 00.01, + 00.10, + 01.00, + 10.00, + 100.0 + }; + + /* Set of noise levels to explore (average deviation * 4) */ + double noiseset[6] = { + 0.0, /* Perfect data */ + 0.01, /* 1.0 % */ + 0.02, /* 2.0 % */ + 0.05, /* 5.0 % */ + 0.10, /* 10.0 % */ + 0.20, /* 20.0 % */ + }; + + + printf("Verifying optimised smoothness factors\n"); + + /* For dimensions */ + for (di = 1; di <= 4; di++) { + + its = trials[di-1]; + res = reses[di-1]; + + printf("Tests %d\n",its); + printf("Dimensions %d\n",di); + printf("RSPL resolution %d\n",res); + + /* For number of sample points */ + for (i = 0; i < 20; i++) { + ntps = nset[di-1][i]; + + if (ntps == 0) + break; + + printf("\nNo. Sample points %d, norm %8.2f\n",ntps, pow((double)ntps, 1.0/di)); + + /* For noise levels */ + for (j = 0; j < 6; j++) { + double rmse[20]; + double bfit; + + noise = noiseset[j]; + printf("Noise volume %f%%, average deviation %f%%\n",noise * 100.0, noise * 25.0); + + /* For smooth factors */ + for (k = 0; k < 5; k++) { + smooth = smset[k]; + + printf("Extra smooth %f, ",smooth); fflush(stdout); + + do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, smooth, twopass, extra); + + rmse[k] = trmse; + printf("maxerr %f%%, avgerr %f%%, rmserr %f%%\n", + tmaxe * 100.0, tavge * 100.0, trmse * 100.0); + } + bfit = best(5, rmse, smset); + printf("Best smoothness = %9.7f, log10 = %4.1f\n",bfit,log10(bfit)); + } + } + printf("\n"); + } +} + +/* ---------------------------------------------------------------------- */ +void usage(void) { + fprintf(stderr,"Test smoothness factor tuning of RSPL in N Dimensions\n"); + fprintf(stderr,"Author: Graeme W. Gill\n"); + fprintf(stderr,"usage: smtnd [options]\n"); + fprintf(stderr," -v Verbose\n"); + fprintf(stderr," -p Plot graphs\n"); + fprintf(stderr," -z n Do test series ""n"" else single test\n"); + fprintf(stderr," 1 = Underlying smoothness\n"); + fprintf(stderr," 2 = Verify optimised smoothness\n"); + fprintf(stderr," -S Compute smoothness factor instead\n"); + fprintf(stderr," -u Use uniformly distributed noise\n"); + fprintf(stderr," -d n Test ""d"" dimension, 1-4 (default 1)\n"); + fprintf(stderr," -t n Test ""n"" random functions (default 1)\n"); + fprintf(stderr," -r res Rspl resolution (defaults 129, 65, 33, 17)\n"); + fprintf(stderr," -n no Test ""no"" sample points (default 20, 40, 80, 100)\n"); + fprintf(stderr," -a amnt Add total randomness to function value (default 0.0)\n"); + fprintf(stderr," -A n Just do the n'th noise level of series\n"); + fprintf(stderr," -2 Use two pass smoothing\n"); + fprintf(stderr," -x Use extra fitting\n"); + fprintf(stderr," -s smooth RSPL extra smoothness factor to test (default 1.0)\n"); + fprintf(stderr," -g smooth RSPL underlying smoothness factor to test\n"); + exit(1); +} + +int main(int argc, char *argv[]) { + int fa,nfa; /* argument we're looking at */ + int verb = 0; + int plot = 0; + int series = 0; + int unif = 0; + int di = 0; + int its = 1; + int res = -1; + int ntps = 0; + double noise = 0.0; + int nlev = 0; + double smooth = 1.0; + double gsmooth = 0.0; + int twopass = 0; + int extra = 0; + int smfunc = 0; + double trmse, tavge, tmaxe; + + error_program = "smtnd"; + +#ifdef NEVER + { + double rmse[10], smv[10], rv; + + smv[0] = 0.0000100, rmse[0] = 2.566116; + smv[1] = 0.0001000, rmse[1] = 2.528666; + smv[2] = 0.0010000, rmse[2] = 2.489116; + smv[3] = 0.0100000, rmse[3] = 3.409045; + smv[4] = 0.1000000, rmse[4] = 5.727079; + smv[5] = 1.0000000, rmse[5] = 6.653747; + + rv = best(6,rmse, smv); + printf("~1 best = %f\n",rv); + exit(0); + } +#endif + /* Process the arguments */ + for(fa = 1;fa < argc;fa++) { + nfa = fa; /* skip to nfa if next argument is used */ + if (argv[fa][0] == '-') { /* Look for any flags */ + char *na = NULL; /* next argument after flag, null if none */ + + if (argv[fa][2] != '\000') + na = &argv[fa][2]; /* next is directly after flag */ + else { + if ((fa+1) < argc) { + if (argv[fa+1][0] != '-') { + nfa = fa + 1; + na = argv[nfa]; /* next is seperate non-flag argument */ + } + } + } + + if (argv[fa][1] == '?') { + usage(); + + } else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') { + verb = 1; + + } else if (argv[fa][1] == 'p' || argv[fa][1] == 'P') { + plot = 1; + + } else if (argv[fa][1] == 'u' || argv[fa][1] == 'U') { + unif = 1; + + /* Test series */ + } else if (argv[fa][1] == 'z' || argv[fa][1] == 'Z') { + fa = nfa; + if (na == NULL) usage(); + series = atoi(na); + if (series <= 0) usage(); + + /* Compute smoothness factor */ + } else if (argv[fa][1] == 'S') { + smfunc = 1; + + /* Dimension */ + } else if (argv[fa][1] == 'd' || argv[fa][1] == 'D') { + fa = nfa; + if (na == NULL) usage(); + di = atoi(na); + if (di <= 0 || di > 4) usage(); +printf("~1 Got -d %s = %d\n",na,di); + + /* Number of tests */ + } else if (argv[fa][1] == 't' || argv[fa][1] == 'T') { + fa = nfa; + if (na == NULL) usage(); + its = atoi(na); + if (its <= 0) usage(); + + /* Resolution */ + } else if (argv[fa][1] == 'r' || argv[fa][1] == 'R') { + fa = nfa; + if (na == NULL) usage(); + res = atoi(na); + if (res <= 0) usage(); + + /* Number of sample points */ + } else if (argv[fa][1] == 'n' || argv[fa][1] == 'N') { + fa = nfa; + if (na == NULL) usage(); + ntps = atoi(na); + if (ntps <= 0) usage(); + + /* Randomness */ + } else if (argv[fa][1] == 'a') { + fa = nfa; + if (na == NULL) usage(); + noise = atof(na); + if (noise < 0.0) usage(); + + /* Series Noise Level */ + } else if (argv[fa][1] == 'A') { + fa = nfa; + if (na == NULL) usage(); + nlev = atoi(na); + if (noise < 0) usage(); + + } else if (argv[fa][1] == '2') { + twopass = 1; + + } else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') { + extra = 1; + + /* Extra smooth factor */ + } else if (argv[fa][1] == 's') { + fa = nfa; + if (na == NULL) usage(); + smooth = atof(na); + if (smooth < 0.0) usage(); + + /* Underlying smoothnes factor */ + } else if (argv[fa][1] == 'g') { + fa = nfa; + if (na == NULL) usage(); + smooth = atof(na); + if (gsmooth < 0.0) usage(); + + } else + usage(); + } else + break; + } + + if (series > 0) { + if (series == 1) + do_series_1(unif, di, ntps, nlev, twopass, extra); + else if (series == 2) + do_series_2(unif, twopass, extra); + else + error("Unknown series %d\n",series); + return 0; + } + + if (di == 0) + di = 1; + + if (res < 0) { + if (di == 1) + res = 129; + else if (di == 2) + res = 65; + else if (di == 3) + res = 33; + else + res = 17; + } + + if (ntps <= 0) { + if (di == 1) + ntps = 20; + else if (di == 2) + ntps = 40; + else if (di == 3) + ntps = 60; + else + ntps = 80; + } + + if (smfunc) { + double sm; + + if (verb) { + printf("Dimensions %d\n",di); + printf("Tests %d\n",its); + printf("Grid resolution %d\n",res); + } + + sm = do_stest(verb, di, its, res); + + printf("Results: smoothness factor = %f\n",sm); + + } else { + + if (verb) { + printf("Dimensions %d\n",di); + printf("Tests %d\n",its); + printf("RSPL resolution %d\n",res); + printf("No. Sample points %d (norm %f)\n",ntps, pow((double)ntps, 1.0/di)); + printf("Noise volume %f\n",noise); + if (gsmooth > 0.0) + printf("Underlying smooth %f\n",gsmooth); + else + printf("Extra smooth %f\n",smooth); + } + + if (gsmooth > 0.0) + do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, -gsmooth, twopass, extra); + else + do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, smooth, twopass, extra); + + printf("Results: maxerr %f%%, avgerr %f%%, rmserr %f%%\n", + tmaxe * 100.0, tavge * 100.0, trmse * 100.0); + } + + return 0; +} + +/* Do one set of tests and return the results */ +static void do_test( + double *trmse, /* RETURN total RMS error */ + double *tmaxe, /* RETURN total maximum error */ + double *tavge, /* RETURN total average error */ + int verb, /* Verbosity */ + int plot, /* Plot graphs */ + int di, /* Dimensions */ + int its, /* Number of function tests */ + int res, /* RSPL grid resolution */ + int ntps, /* Number of sample points */ + double noise, /* Sample point noise volume (total = 4 x average deviation) */ + int unif, /* NZ if uniform rather than standard deistribution noise */ + double smooth, /* Smoothness to test, +ve for extra, -ve for underlying */ + int twopass, /* Two pass flag */ + int extra /* Extra fit flag */ +) { + funcp fp; /* Function parameters */ + sobol *so; /* Sobol sequence generator */ + co *tps = NULL; + rspl *rss; /* Multi-resolution regularized spline structure */ + datai low,high; + double avgdev[MXDO]; + int gres[MXDI]; + int i, j, it; + int flags = RSPL_NOFLAGS; + + *trmse = 0.0; + *tmaxe = 0.0; + *tavge = 0.0; + + for (j = 0; j < di; j++) { + low[j] = 0.0; + high[j] = 1.0; + gres[j] = res; + } + + if (twopass) + flags |= RSPL_2PASSSMTH; + + if (extra) + flags |= RSPL_EXTRAFIT2; + + if ((so = new_sobol(di)) == NULL) + error("Creating sobol sequence generator failed"); + + for (it = 0; it < its; it++) { + double rmse, avge, maxe; + + /* Make repeatable by setting random seed before a test set. */ + rand32(0x12345678 + 0x1000 * it); + + /* New function */ + setup_func(&fp, di); + + /* Create the object */ + rss = new_rspl(RSPL_NOFLAGS,di, 1); + + /* Create the list of sampling points */ + if ((tps = (co *)malloc(ntps * sizeof(co))) == NULL) + error ("malloc failed"); + + so->reset(so); + + if (verb) printf("Generating the sample points\n"); + + for (i = 0; i < ntps; i++) { + so->next(so, tps[i].p); + tps[i].v[0] = lookup_func(&fp, tps[i].p); + if (unif) + tps[i].v[0] += d_rand(-0.5 * noise, 0.5 * noise); + else + tps[i].v[0] += noise * 0.25 * 1.2533 * norm_rand(); + } + + /* Fit to scattered data */ + if (verb) printf("Fitting the scattered data\n"); + avgdev[0] = 0.25 * noise; + rss->fit_rspl(rss, + flags, /* Non-mon and clip flags */ + tps, /* Test points */ + ntps, /* Number of test points */ + low, high, gres, /* Low, high, resolution of grid */ + low, high, /* Default data scale */ + smooth, /* Smoothing to test */ + avgdev, /* Average deviation */ + NULL); /* iwidth */ + + /* Plot out function values */ + if (plot) { + int slice; + printf("Black is target, Red is rspl\n"); + for (slice = 0; slice < (di+1); slice++) { + co tp; /* Test point */ + double x[PLOTRES]; + double ya[PLOTRES]; + double yb[PLOTRES]; + double yc[PLOTRES]; + double pp[MXDI], p1[MXDI], p2[MXDI], ss[MXDI]; + int n = PLOTRES; + + /* setup slices on each axis at 0.5 and diagonal */ + if (slice < di) { + for (j = 0; j < di; j++) + p1[j] = p2[j] = 0.5; + p1[slice] = 0.0; + p2[slice] = 1.0; + printf("Slice along axis %d\n",slice); + } else { + for (j = 0; j < di; j++) { + p1[j] = 0.0; + p2[j] = 1.0; + } + printf("Slice along diagonal\n"); + } + + for (j = 0; j < di; j++) { + ss[j] = (p2[j] - p1[j])/n; + pp[j] = p1[j]; + } + + for (i = 0; i < n; i++) { + double vv = i/(n-1.0); + x[i] = vv; + + /* Reference */ + ya[i] = lookup_func(&fp, pp); + + /* RSPL aproximation */ + for (j = 0; j < di; j++) + tp.p[j] = pp[j]; + + if (rss->interp(rss, &tp)) + tp.v[0] = -0.1; + yb[i] = tp.v[0]; + + /* Crude way of setting the scale: */ + yc[i] = 0.0; + if (i == (n-1)) + yc[0] = 1.0; + + for (j = 0; j < di; j++) + pp[j] += ss[j]; + } + + /* Plot the result */ + do_plot(x,ya,yb,yc,n); + } + } + + /* Compute statistics */ + rmse = 0.0; + avge = 0.0; + maxe = 0.0; +// so->reset(so); + + /* Fit to scattered data */ + if (verb) printf("Fitting the scattered data\n"); + for (i = 0; i <100000; i++) { + co tp; /* Test point */ + double aa, bb, err; + + so->next(so, tp.p); + + /* Reference */ + aa = lookup_func(&fp, tp.p); + + /* RSPL aproximation */ + rss->interp(rss, &tp); + bb = tp.v[0]; + + err = fabs(aa - bb); + avge += err; + rmse += err * err; + if (err > maxe) + maxe = err; + } + avge /= (double)i; + rmse /= (double)i; + + if (verb) + printf("Dim %d, res %d, noise %f, points %d, maxerr %f%%, rmserr %f%%, avgerr %f%%\n", + di, res, noise, ntps, maxe * 100.0, sqrt(rmse) * 100.0, avge * 100.0); + + *trmse += rmse; + *tmaxe += maxe; + *tavge += avge; + + rss->del(rss); + free(tps); + } + so->del(so); + + *trmse = sqrt(*trmse/(double)its); + *tmaxe /= (double)its; + *tavge /= (double)its; +} + +/* Do smoothness scaling check & return results */ +static double do_stest( + int verb, /* Verbosity */ + int di, /* Dimensions */ + int its, /* Number of function tests */ + int res /* RSPL grid resolution */ +) { + funcp fp; /* Function parameters */ + DCOUNT(gc, MXDIDO, di, 1, 1, res-1); + int it; + double atse = 0.0; + + /* Make repeatable by setting random seed before a test set. */ + rand32(0x12345678); + + for (it = 0; it < its; it++) { + double tse; + setup_func(&fp, di); /* New function */ + + DC_INIT(gc) + tse = 0.0; + for (; !DC_DONE(gc);) { + double g[MXDI]; + int e, k; + double y1, y2, y3; + double del; + + for (e = 0; e < di; e++) + g[e] = gc[e]/(res-1.0); + y2 = lookup_func(&fp, g); + + del = 1.0/(res-1.0); + + for (k = 0 ; k < di; k++) { + double err; + + g[k] -= del; + y1 = lookup_func(&fp, g); + g[k] += 2.0 * del; + y3 = lookup_func(&fp, g); + g[k] -= del; + + err = 0.5 * (y3 + y1) - y2; + tse += err * err; + } + + DC_INC(gc); + } + /* Apply adjustments and corrections */ + tse *= pow((res-1.0), 4.0); /* Aprox. geometric resolution factor */ + tse /= pow((res-2.0),(double)di); /* Average squared non-smoothness */ + + if (verb) + printf("smf for it %d = %f\n",it,tse); + atse += tse; + } + + return atse/(double)its; +} + + -- cgit v1.2.3