summaryrefslogtreecommitdiff
path: root/rspl/cw1.c
blob: 2a026c8558e2400b2a47c3f8f5ec80b20e2e9cad (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309

/************************************************/
/* Investigate various curve approximations     */
/************************************************/

/* Discrete regularized spline versions */
/* Test variable grid spacing in 1D */

/* Author: Graeme Gill
 * Date:   4/10/95
 * Date:   5/4/96
 *
 * Copyright 1995, 2013 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
#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 "rspl.h"
#include "plot.h"
#include "ui.h"

double lin();
void usage(void);

#define TRIALS 1	/* Number of random trials */
#define XRES 100		/* Plotting res */


#define PNTS1 10
#define GRES1 100

/* Function that we're approximating */
static double func(double val) {
	double out;
	int sgn = 0;

	if (val < 0.0) {
		sgn = 1;
		val = -val;
	}

	out = pow(val, 2.5);

	if (sgn)
		out = -out;

	return out;
}

/* Grid spacing power */
#define GRIDPOW 1.5

/* Scattered data sample spacing power */
#define SAMPPOW 1.3


co test_points[PNTS1];

double ipos[GRES1];
double *iposes[1] = { ipos };

double powlike(double vv, double pp);

void usage(void) {
	fprintf(stderr,"Test 1D rspl interpolation variable grid spacing\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");
	exit(1);
}

int main(int argc, char *argv[]) {
	int fa,nfa;				/* argument we're looking at */
	int i,j, n;
	double x;
	double xx1[XRES];
	double yy1[6][XRES];
	double xx2[XRES];
	double yy2[6][XRES];
	rspl *rss;		/* incremental solution version */
	datai low,high;
	int gres[MXDI];
	double smooth = 1.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 
				usage();
		} else
			break;
	}

	/* Only one trial */
	for (n = 0; n < TRIALS; n++) {
		double lrand = 0.0;		/* Amount of level randomness */
		int pnts;
		int fres;

		pnts = PNTS1;
		fres = GRES1; 
		gres[0] = fres;

		/* Create the object */
		rss =  new_rspl(RSPL_NOFLAGS,
		                1,				/* di */
		                1);				/* fdi */

		/* For each grid spacing test */
		for (j = 0; j < 2; j++) {

			for (i = 0; i < pnts; i++) {
				double val, out;

				val = i/(pnts - 1.0);

				if (j == 0) {			/* Linear */
//					out = func(val + d_rand(-0.001, 0.001));
					out = func(val);
					test_points[i].p[0] = val;	
					test_points[i].v[0] = out;

				} else {					/* power 2.0 grid spacing */
					val = pow(val, SAMPPOW);	/* sample point spacing */ 

//					out = func(val + d_rand(-0.001, 0.001));
					out = func(val);
					val = powlike(val, 1.0/GRIDPOW);	/* To grid spacing */ 
					test_points[i].p[0] = val;
					test_points[i].v[0] = out;
				}
			}

			/* Set grid position info */
			for (i = 0; i < gres[0]; i++) {
				double val, out;

				val = i/(gres[0] - 1.0);
				if (j == 1)
					val = powlike(val, 1.0/GRIDPOW);	/* To grid spacing */ 
				ipos[i] = val; 
			}
			
			rss->fit_rspl(rss,
		           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 */
		           j == 0 ? NULL : iposes);				/* iwidth */

			/* Show the result full scale */
			for (i = 0; i < XRES; i++) {
				co tp;	/* Test point */
				x = i/(double)(XRES-1);
				xx1[i] = x;
				yy1[0][i] = func(x);
				tp.p[0] = x;
				if (j == 1)
					tp.p[0] = powlike(tp.p[0], 1.0/GRIDPOW);	/* Grid spacing input */ 
				rss->interp(rss, &tp);
				yy1[1+j][i] = tp.v[0];
				if (yy1[1+j][i] < -0.2)
					yy1[1+j][i] = -0.2;
				else if (yy1[1+j][i] > 1.2)
					yy1[1+j][i] = 1.2;
			}

			/* Show the result magnified */
			for (i = 0; i < XRES; i++) {
				co tp;	/* Test point */
				x = i/(double)(XRES-1);
				x *= 0.1;
				xx2[i] = x;
				yy2[0][i] = func(x);
				tp.p[0] = x;
				if (j == 1)
					tp.p[0] = powlike(tp.p[0], 1.0/GRIDPOW);	/* Grid spacing input */ 
				rss->interp(rss, &tp);
				yy2[1+j][i] = tp.v[0];
				if (yy2[1+j][i] < -0.2)
					yy2[1+j][i] = -0.2;
				else if (yy2[1+j][i] > 1.2)
					yy2[1+j][i] = 1.2;
			}
		}
	
		printf("Full scale: Black = ref, Red = even, Green = pow2\n");
		do_plot6(xx1,yy1[0],yy1[1],yy1[2],NULL,NULL,NULL,XRES);

		printf("Magnified: Black = ref, Red = even, Green = pow2\n");
		do_plot6(xx2,yy2[0],yy2[1],yy2[2],NULL,NULL,NULL,XRES);
	}
	return 0;
}


/* A power-like function, based on Graphics Gems adjustment curve. */
/* Avoids "toe" problem of pure power. */
/* Adjusted so that "power" 2 and 0.5 agree with real power at 0.5 */

double powlike(double vv, double pp) {
	double tt, g;

	if (pp >= 1.0) {
		g = 2.0 * (pp - 1.0);
		vv = vv/(g - g * vv + 1.0);
	} else {
		g = 2.0 - 2.0/pp;
		vv = (vv - g * vv)/(1.0 - g * vv);
	}

	return vv;
}


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