1 /* splinef.i
2 * $Id: splinef.i,v 1.2 2007-04-07 01:16:12 dhmunro Exp $
3 * piecewise cubic interpolation functions
4 */
5 /* Copyright (c) 2007, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 /* This package extends the original spline.i functions, concentrating
12 * on how arbitrary piecewise cubic functions may be used as a natural
13 * extension of the piecewise linear interp() built-in function.
14 * Given a sequence of abcissa values x(i), the corresponding function
15 * values y(i) uniquely determine a piecewise linear function. Knowing
16 * both the function and its first derivatives, y(i) and dydx(i),
17 * uniquely determines a piecewise cubic function.
18 *
19 * Unlike a true "spline", a general piecewise cubic curve has a
20 * discontinuous second derivative. True spline curves depend only on
21 * the function values y(i) at x(i); the derivatives dydx(i) are
22 * computed to force continuity of the second derivative at x(i).
23 * The less continuous general piecewise cubic permits you to make a
24 * more robust fit, which is still smoother than piecewise linear, at
25 * the cost of having to provide dydx(i) values by some other means.
26 * The splinelsq function can select y(i) and dydx(i) to have least
27 * square residuals of the piecewise cubic from a larger set of data,
28 * like the fitlsq function in fitlsq.i does for piecewise linear fits.
29 *
30 * Note that a pair of piecewise cubic curves is a Bezier curve; that is,
31 * if x(t) and y(t) are piecewise cubic, then (x,y) is a (cubic) Bezier.
32 * These are the curves generated by the Postscript arcto operator.
33 * If (x0,y0) is an initial point, (x1,y1) and (x2,y2) are the Bezier
34 * control points, and (x3,y3) is the final point, then the spline
35 * with 0<t<1 has dxdt0=3*(x1-x0), dydt0=3*(y1-y0), dxdt3=3*(x3-x2),
36 * and dydt3=3*(y3-y2).
37 */
38
splinef(dydx,y,x,xp)39 func splinef(dydx, y, x, xp)
40 /* DOCUMENT yp = splinef(dydx, y, x, xp)
41 * or yp = splinef(x_y_dydx, xp)
42 * returns piecewise cubic function specified by DYDX, Y, X at
43 * the points XP. Extrapolation beyond the extreme endpoints of X
44 * is linear, with slope equal to the final value of DYDX. The
45 * return value dimensions are the same as the dimensions of XP.
46 *
47 * In the second form, X_Y_DYDX is a 3-by-nknots array of [x,y,dydx]
48 * values. The values of X in either case must either increase or
49 * decrease monotonically.
50 *
51 * SEE ALSO: interp, splined, splinei, splinelsq
52 */
53 {
54 local c0, c1, c2, c3, cm1;
55 spline_coef;
56 return poly(x, c0, c1, c2, c3);
57 }
58
splined(dydx,y,x,xp)59 func splined(dydx, y, x, xp)
60 /* DOCUMENT yp = splined(dydx, y, x, xp)
61 * or yp = splined(x_y_dydx, xp)
62 * returns the derivative of the piecewise cubic function specified
63 * by DYDX, Y, X at the points XP. Extrapolation beyond the extreme
64 * endpoints of X is linear, so splined gives the final value of DYDX.
65 * The return value dimensions are the same as the dimensions of XP.
66 *
67 * In the second form, X_Y_DYDX is a 3-by-nknots array of [x,y,dydx]
68 * values. The values of X in either case must either increase or
69 * decrease monotonically.
70 *
71 * SEE ALSO: splined, splinei
72 */
73 {
74 local c0, c1, c2, c3, cm1;
75 spline_coef;
76 return poly(x, c1, c2+c2, 3.*c3);
77 }
78
splinei(dydx,y,x,xp)79 func splinei(dydx, y, x, xp)
80 /* DOCUMENT yp = splinei(dydx, y, x, xp)
81 * or yp = splinei(x_y_dydx, xp)
82 * returns the integral of the piecewise cubic function specified
83 * by DYDX, Y, X at the points XP. The integral is quadratic beyond
84 * the extreme endpoints of X, and zero at X(1). The dimensions of
85 * the return value are the same as the dimensions of XP.
86 * This is the cubic analog of the integ function.
87 *
88 * In the second form, X_Y_DYDX is a 3-by-nknots array of [x,y,dydx]
89 * values. The values of X in either case must either increase or
90 * decrease monotonically.
91 *
92 * SEE ALSO: integ, splinef, splined
93 */
94 {
95 local c0, c1, c2, c3, cm1;
96 cm1 = 1;
97 spline_coef;
98 return poly(x, cm1, c0, 0.5*c1, (1./3.)*c2, 0.25*c3);
99 }
100
101 func spline_coef
102 /* DOCUMENT spline_coef
103 * is the worker for the splinef, splined, and splinei functions.
104 * If you need to compute both function and derivative or integral,
105 * you will improve performance using spline_coef. See the source
106 * code for those functions for usage.
107 *
108 * SEE ALSO: splinef, splined, splinei
109 */
110 {
111 extern c0, c1, c2, c3, cm1, dydx, y, x;
112 if (is_void(xp)) {
113 xp = y;
114 x = dydx(1,);
115 y = dydx(2,);
116 dydx = dydx(3,);
117 }
118 u = digitize(xp, x);
119 l = max(u - 1, 1);
120 u = min(u, numberof(x));
121
122 xl = x(l);
123 c0 = y(l);
124 c1 = dydx(l);
125
126 if (cm1) {
127 c2 = x(dif);
128 cm1 = ((y(zcen) - (1./12.)*dydx(dif)*c2)*c2)(cum)(l);
129 }
130
131 x = x(u) - xl; /* 0 outside endpoints of x */
132 c2 = double(!x);
133 x = (1. - c2)/(x + c2);
134 c2 = (y(u) - c0) * x;
135 c3 = dydx(u) + c1 - c2 - c2;
136 c2 = (c2 - c1 - c3) * x;
137 c3 *= x * x;
138
139 x = xp - xl;
140 }
141
142 /* question: any simple way to get continuous 2nd derivative fit? */
143 func splinelsq(y, x, xfit, weight=, y0=, dydx0=, y1=, dydx1=, constrain=)
144 /* DOCUMENT x_y_dydx = splinelsq(y, x, xfit)
145 ...
146 yp = splinef(x_y_dydx, xp)
147 performs a least squares fit to the data points (X, Y). The input
148 XFIT are the abcissas of the piecewise cubic function with knot
149 points XFIT which is the least squares best fit to the data (X,Y).
150 The XFIT must be strictly increase or decrease.
151
152 Any points in XFIT with no data points in the intervals on
153 either side will be removed.
154
155 A weight= keyword of the same length as X and Y may be supplied in
156 order to weight the various data points differently; a typical
157 WEIGHT function is 1/sigma^2 where sigma are the standard deviations
158 associated with the Y values.
159
160 You can specify y0=, dydx0=, y1=, and dydx1= keywords to fix the
161 value of the function or its derivative at the first (0) or last (1)
162 endpoint. Be sure there is at least one point in the final
163 interval so that the XFIT at the endpoint is not removed.
164
165 More generally, you can specify a constrain= keyword. The value
166 of constrain is a hook function which will be called just before
167 the matrix solve. Your constrain subroutine will be passed no
168 arguments, but it can access and modify the mat and rhs variables.
169
170 SEE ALSO: splinef
171 */
172 {
173 xfit = double(xfit);
174 np1 = numberof(xfit)+1;
175
176 /* bin the input data into the xfit */
177 l = digitize(x, xfit);
178 /* remove any points between adjacent empty bins */
179 list = histogram(l, top=np1);
180 mask = list(2:0) + list(1:-1);
181 list = where(mask);
182 nfit = numberof(xfit);
183 if (numberof(list) < nfit) {
184 if (numberof(list) < 2) error, "all data points outside xfit";
185 xfit = xfit(list);
186 np1 = numberof(xfit)+1;
187 l = digitize(x, xfit);
188 }
189 /* create extended xfit for which l and l+1 may be used directly */
190 dx = xfit(0) - xfit(1);
191 xl = min(min(x), min(xfit)) - abs(dx);
192 xu = max(max(x), max(xfit)) + abs(dx);
193 xx = (dx>0.)? grow(xl, xfit, xu) : grow(xu, xfit, xl);
194
195 xl = xx(l);
196 xu = xx(l+1);
197 dx = xu - xl;
198 rdx = 1./dx;
199
200 ld = xhi = x-xl;
201 lx = ld * rdx;
202 ud = lx*lx;
203 lx = 1.-lx;
204 ld *= lx*lx;
205 uf = (1.+lx+lx)*ud;
206 ud *= xlo = -lx*dx;
207 lf = 1.-uf;
208
209 mist = where(l==1);
210 if (numberof(mist)) {
211 lf(mist) = ld(mist) = 0.;
212 uf(mist) = 1.;
213 ud(mist) = xlo(mist);
214 xlo = [];
215 }
216 mist = where(l==np1);
217 if (numberof(mist)) {
218 uf(mist) = ud(mist) = 0.;
219 lf(mist) = 1.;
220 ld(mist) = xhi(mist);
221 xhi = [];
222 }
223
224 if (is_void(weight)) weight = 1.0;
225 lfy = histogram(l, lf*y*weight, top=np1);
226 ldy = histogram(l, ld*y*weight, top=np1);
227 ufy = histogram(l, uf*y*weight, top=np1);
228 udy = histogram(l, ud*y*weight, top=np1);
229
230 lflf = histogram(l, lf*lf*weight, top=np1);
231 lfld = histogram(l, lf*ld*weight, top=np1);
232 lfuf = histogram(l, lf*uf*weight, top=np1);
233 lfud = histogram(l, lf*ud*weight, top=np1);
234 ldld = histogram(l, ld*ld*weight, top=np1);
235 lduf = histogram(l, ld*uf*weight, top=np1);
236 ldud = histogram(l, ld*ud*weight, top=np1);
237 ufuf = histogram(l, uf*uf*weight, top=np1);
238 ufud = histogram(l, uf*ud*weight, top=np1);
239 udud = histogram(l, ud*ud*weight, top=np1);
240
241 n = 2*numberof(xfit);
242 rhs = array(0., n);
243 mat = array(rhs, n);
244 n2 = n*n;
245 dn = 2*(n+1);
246 mat(1:n2:dn) = lflf(2:0) + ufuf(1:-1);
247 mat(2:n2:dn) = mat(n+1:n2:dn) = lfld(2:0) + ufud(1:-1);
248 mat(n+2:n2:dn) = ldld(2:0) + udud(1:-1);
249 mat(3:n2-n:dn) = mat(dn-1:n2:dn) = lfuf(2:-1);
250 mat(4:n2-n:dn) = mat(n+dn-1:n2:dn) = lfud(2:-1);
251 mat(n+3:n2:dn) = mat(dn:n2:dn) = lduf(2:-1);
252 mat(n+4:n2:dn) = mat(n+dn:n2:dn) = ldud(2:-1);
253 rhs(1:n:2) = lfy(2:0) + ufy(1:-1);
254 rhs(2:n:2) = ldy(2:0) + udy(1:-1);
255
256 if (!is_void(y0)) {
257 if (list(1)!=1) error, "y0= keyword with no data in first interval";
258 rhs(1) = y0;
259 mat(1,) = 0.0;
260 mat(1,1) = 1.0;
261 }
262 if (!is_void(dydx0)) {
263 if (list(1)!=1) error, "dydx0= keyword with no data in first interval";
264 rhs(2) = dydx0;
265 mat(2,) = 0.0;
266 mat(2,2) = 1.0;
267 }
268 if (!is_void(y1)) {
269 if (list(0)!=nfit) error, "y1= keyword with no data in last interval";
270 rhs(-1) = y1;
271 mat(-1,) = 0.0;
272 mat(-1,-1) = 1.0;
273 }
274 if (!is_void(dydx1)) {
275 if (list(0)!=nfit) error, "dydx1= keyword with no data in last interval";
276 rhs(0) = dydx1;
277 mat(0,) = 0.0;
278 mat(0,0) = 1.0;
279 }
280 if (!is_void(constrain)) constrain;
281
282 fd = LUsolve(mat, rhs);
283 xfd = xfit(-:1:3,);
284 xfd(2,) = fd(1:0:2);
285 xfd(3,) = fd(2:0:2);
286 return xfd;
287 }
288