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/c1.c | 396 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 396 insertions(+) create mode 100644 rspl/c1.c (limited to 'rspl/c1.c') diff --git a/rspl/c1.c b/rspl/c1.c new file mode 100644 index 0000000..8172cef --- /dev/null +++ b/rspl/c1.c @@ -0,0 +1,396 @@ + +/************************************************/ +/* Investigate various curve approximations */ +/************************************************/ + +/* Discrete regularized spline versions */ +/* Standard test + Random testing */ + +/* Author: Graeme Gill + * Date: 4/10/95 + * Date: 5/4/96 + * + * Copyright 1995, 1996 Graeme W. Gill + * + * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :- + * see the License.txt file for licencing details. + */ + +#undef DIAG +#undef DIAG2 +#undef GLOB_CHECK +#undef RES2 /* Do multiple test at various resolutions */ +#define AVGDEV 0.0 /* Average deviation of function data */ + +#include +#include +#include +#include +#include "copyright.h" +#include "aconfig.h" +#include "numlib.h" +#include "plot.h" +#include "rspl.h" + +double lin(); +void usage(void); + +#define TRIALS 20 /* Number of random trials */ +#define SKIP 0 /* Number of random trials to skip */ + +#define MIN_PNTS 5 +#define MAX_PNTS 40 + +#define MIN_RES 20 +#define MAX_RES 2000 + +double xa[MAX_PNTS]; +double ya[MAX_PNTS]; + +#define XRES 100 + +#define PNTS1 10 +#define GRES1 400 +//#define GRES 800 +double t1xa[PNTS1] = { 0.2, 0.25, 0.30, 0.35, 0.40, 0.44, 0.48, 0.51, 0.64, 0.75 }; +double t1ya[PNTS1] = { 0.3, 0.35, 0.4, 0.41, 0.42, 0.46, 0.5, 0.575, 0.48, 0.75 }; + +#ifndef NEVER + +// Reverse in x */ +#define PNTS2 10 +#define GRES2 400 +double t2xa[PNTS2] = { 0.25, 0.36, 0.49, 0.52, 0.56, 0.60, 0.65, 0.70, 0.75, 0.8 }; +double t2ya[PNTS2] = { 0.75, 0.48, 0.575, 0.5, 0.46, 0.42, 0.41, 0.4, 0.35, 0.3 }; + +#else + +#define PNTS2 10 +#define GRES2 400 +// reverse in y +double t2xa[PNTS2] = { 0.2, 0.25, 0.30, 0.35, 0.40, 0.44, 0.48, 0.51, 0.64, 0.75 }; +double t2ya[PNTS2] = { 0.7, 0.65, 0.6, 0.59, 0.58, 0.54, 0.5, 0.425, 0.52, 0.25 }; + +#endif /* NEVER */ + +//#define PNTS2 2 +//#define GRES2 5 +//double t2xa[PNTS2] = { 0.0, 1.0 }; +//double t2ya[PNTS2] = { 0.33, 0.66 }; + +co test_points[MAX_PNTS]; + +double lin(double x, double xa[], double ya[], int n); + +void usage(void) { + fprintf(stderr,"Test 1D rspl interpolation\n"); + fprintf(stderr,"Author: Graeme W. Gill\n"); + fprintf(stderr,"usage: c1 [options]\n"); + fprintf(stderr," -s smooth Use given smoothness (default 1.0)\n"); + fprintf(stderr," -2 Use two pass smoothing\n"); + fprintf(stderr," -x Use extra fitting\n"); + exit(1); +} + +int main(int argc, char *argv[]) { + int fa,nfa; /* argument we're looking at */ + int i,j, n; + double x; + double xx[XRES]; + double yy[6][XRES]; + rspl *rss; /* incremental solution version */ + datai low,high; + int gres[MXDI]; + double smooth = 1.0; + int twopass = 0; + int extra = 0; + double avgdev[MXDO]; + + low[0] = 0.0; + high[0] = 1.0; + avgdev[0] = AVGDEV; + + error_program = "c1"; + check_if_not_interactive(); + + /* 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(); + + /* smoothness */ + else if (argv[fa][1] == 's' || argv[fa][1] == 'S') { + fa = nfa; + if (na == NULL) usage(); + smooth = atof(na); + } + + else if (argv[fa][1] == '2') { + twopass = 1; + } + else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') { + extra = 1; + } + else + usage(); + } else + break; + } + + for (n = 0; n < TRIALS; n++) { + double lrand = 0.0; /* Amount of level randomness */ + int pnts; + int fres; + + if (n == 0) { /* Standard versions */ +#ifdef NEVER /* Doubled up points */ + pnts = 2 * PNTS; + fres = GRES; + for (i = 0; i < pnts; i++) { + xa[i * 2 + 0] = t1xa[i] - 0.01; + ya[i * 2 + 0] = t1ya[i]; + xa[i * 2 + 1] = t1xa[i] + 0.01; + ya[i * 2 + 1] = t1ya[i]; + } +#else + pnts = PNTS1; + fres = GRES1; + for (i = 0; i < pnts; i++) { + xa[i] = t1xa[i]; + ya[i] = t1ya[i]; + } +#endif + printf("Trial %d, points = %d, res = %d, level randomness = %f\n",n,pnts,fres,lrand); + } else if (n == 1) { /* Second test versions */ + pnts = PNTS2; + fres = GRES2; + for (i = 0; i < pnts; i++) { + xa[i] = t2xa[i]; + ya[i] = t2ya[i]; + } + printf("Trial %d, points = %d, res = %d, level randomness = %f\n",n,pnts,fres,lrand); + } else { /* Random versions */ + lrand = d_rand(0.0,0.1); /* Amount of level randomness */ + pnts = i_rand(MIN_PNTS,MAX_PNTS); + fres = i_rand(MIN_RES,MAX_RES); + + printf("Trial %d, points = %d, res = %d, level randomness = %f\n",n,pnts,fres,lrand); + + /* Create X values */ + xa[0] = d_rand(0.5,1.0); + for (i = 1; i < pnts; i++) + xa[i] = xa[i-1] + d_rand(0.5,1.0); + for (i = 0; i < pnts; i++) /* Divide out */ + xa[i] = (xa[i]/xa[pnts-1]); + + /* Create y values */ + ya[0] = xa[0]; + for (i = 0; i < pnts; i++) + ya[i] = ya[i-1] + d_rand(0.2,1.0) + d_rand(-0.2,0.3) + d_rand(-0.2,0.3); + for (i = 0; i < pnts; i++) /* Divide out */ + ya[i] = (ya[i]/ya[pnts-1]); + } + + if (n < SKIP) + continue; + + /* Create the object */ + rss = new_rspl(RSPL_NOFLAGS, + 1, /* di */ + 1); /* fdi */ + + for (i = 0; i < pnts; i++) { + test_points[i].p[0] = xa[i]; + test_points[i].v[0] = ya[i]; + } + gres[0] = fres; + +#ifdef RES2 + if (n != 0) { +#endif + /* Fit to scattered data */ + rss->fit_rspl(rss, + 0 | (twopass ? RSPL_2PASSSMTH : 0) + | (extra ? RSPL_EXTRAFIT2 : 0) , + test_points, /* Test points */ + pnts, /* Number of test points */ + low, high, gres, /* Low, high, resolution of grid */ + NULL, NULL, /* Default data scale */ + smooth, /* Smoothing */ + avgdev, /* Average deviation */ + NULL); /* iwidth */ + + + /* Display the result */ + for (i = 0; i < XRES; i++) { + co tp; /* Test point */ + x = i/(double)(XRES-1); + xx[i] = x; + yy[0][i] = lin(x,xa,ya,pnts); + tp.p[0] = x; + rss->interp(rss, &tp); + yy[1][i] = tp.v[0]; + if (yy[1][i] < -0.2) + yy[1][i] = -0.2; + else if (yy[1][i] > 1.2) + yy[1][i] = 1.2; + } + + do_plot(xx,yy[0],yy[1],NULL,XRES); + +#ifdef RES2 + } else { /* Multiple resolution version */ + int gresses[5]; + for (j = 0; j < 5; j++) { +#ifndef NEVER + if (j == 0) + gres[0] = fres/8; + else if (j == 1) + gres[0] = fres/4; + else if (j == 2) + gres[0] = fres/2; + else if (j == 3) + gres[0] = fres; + else + gres[0] = fres * 2; +#else /* Check sensitivity to griding of data points */ + if (j == 0) + gres[0] = 192; + else if (j == 1) + gres[0] = 193; + else if (j == 2) + gres[0] = 194; + else if (j == 3) + gres[0] = 195; + else + gres[0] = 196; +#endif + gresses[j] = gres[0]; + + rss->fit_rspl(rss, + 0 | (twopass ? RSPL_2PASSSMTH : 0) + | (extra ? RSPL_EXTRAFIT2 : 0) , + 0, + test_points, /* Test points */ + pnts, /* Number of test points */ + low, high, gres, /* Low, high, resolution of grid */ + NULL, NULL, /* Default data scale */ + smooth, /* Smoothing */ + avgdev, /* Average deviation */ + NULL); /* iwidth */ + + /* Get the result */ + for (i = 0; i < XRES; i++) { + co tp; /* Test point */ + x = i/(double)(XRES-1); + xx[i] = x; + yy[0][i] = lin(x,xa,ya,pnts); + tp.p[0] = x; + rss->interp(rss, &tp); + yy[1+j][i] = tp.v[0]; + if (yy[1+j][i] < -0.2) + yy[1+j][i] = -0.2; + else if (yy[1+j][i] > 1.2) + yy[1+j][i] = 1.2; + } + } + + printf("Black = lin, Red = %d, Green = %d, Blue = %d, Yellow = %d, Purple = %d\n", + gresses[0], gresses[1], gresses[2], gresses[3], gresses[4]); + do_plot6(xx,yy[0],yy[1],yy[2],yy[3],yy[4],yy[5],XRES); + } +#endif /* RES2 */ + } /* next trial */ + return 0; +} + + +/* Simple linear interpolation */ +double +lin( +double x, +double xa[], +double ya[], +int n) { + int i; + double y; + + if (x < xa[0]) + return ya[0]; + else if (x > xa[n-1]) + return ya[n-1]; + + for (i = 0; i < (n-1); i++) + if (x >=xa[i] && x <= xa[i+1]) + break; + + x = (x - xa[i])/(xa[i+1] - xa[i]); + + y = ya[i] + (ya[i+1] - ya[i]) * x; + + return y; +} + + +/******************************************************************/ +/* Error/debug output routines */ +/******************************************************************/ + +/* Next u function done with optimization */ + +/* Structure to hold data for optimization function */ +struct _edatas { + rspl *rss; + int j; + }; typedef struct _edatas edatas; + +#ifdef GLOB_CHECK +/* Overall Global optimization method */ +/* Definition of the optimization function handed to powell() */ +double efunc2(void *edata, double p[]) + { + int j; + double rv; + rspl *rss = (rspl *)edata; + for (j = 0; j < rss->nig; j++) /* Ugg */ + rss->u[j].v = p[j]; + rv = rss->efactor(rss); +#ifdef DIAG2 + /* printf("%c%e",cr_char,rv); */ + printf("%e\n",rv); +#endif + return rv; + } + +solveu(rss) +rspl *rss; + { + int j; + double *cp; + double *s; + + cp = dvector(0,rss->nig); + s = dvector(0,rss->nig); + for (j = 0; j < rss->nig; j++) /* Ugg */ + { + cp[j] = rss->u[j].v; + s[j] = 0.1; + } + powell(rss->nig,cp,s,1e-7,1000,efunc2,(void *)rss); + } +#endif /* GLOB_CHECK */ -- cgit v1.2.3