path: root/rspl/c1.c

/************************************************/
/* 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 <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <math.h>
#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 */