1 subroutine percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp, 2 * wrk,lwrk,iwrk,ier) 3c given the set of data points (x(i),y(i)) and the set of positive 4c numbers w(i),i=1,2,...,m-1, subroutine percur determines a smooth 5c periodic spline approximation of degree k with period per=x(m)-x(1). 6c if iopt=-1 percur calculates the weighted least-squares periodic 7c spline according to a given set of knots. 8c if iopt>=0 the number of knots of the spline s(x) and the position 9c t(j),j=1,2,...,n is chosen automatically by the routine. the smooth- 10c ness of s(x) is then achieved by minimalizing the discontinuity 11c jumps of the k-th derivative of s(x) at the knots t(j),j=k+2,k+3,..., 12c n-k-1. the amount of smoothness is determined by the condition that 13c f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, with s a given non- 14c negative constant, called the smoothing factor. 15c the fit s(x) is given in the b-spline representation (b-spline coef- 16c ficients c(j),j=1,2,...,n-k-1) and can be evaluated by means of 17c subroutine splev. 18c 19c calling sequence: 20c call percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk, 21c * lwrk,iwrk,ier) 22c 23c parameters: 24c iopt : integer flag. on entry iopt must specify whether a weighted 25c least-squares spline (iopt=-1) or a smoothing spline (iopt= 26c 0 or 1) must be determined. if iopt=0 the routine will start 27c with an initial set of knots t(i)=x(1)+(x(m)-x(1))*(i-k-1), 28c i=1,2,...,2*k+2. if iopt=1 the routine will continue with 29c the knots found at the last call of the routine. 30c attention: a call with iopt=1 must always be immediately 31c preceded by another call with iopt=1 or iopt=0. 32c unchanged on exit. 33c m : integer. on entry m must specify the number of data points. 34c m > 1. unchanged on exit. 35c x : real array of dimension at least (m). before entry, x(i) 36c must be set to the i-th value of the independent variable x, 37c for i=1,2,...,m. these values must be supplied in strictly 38c ascending order. x(m) only indicates the length of the 39c period of the spline, i.e per=x(m)-x(1). 40c unchanged on exit. 41c y : real array of dimension at least (m). before entry, y(i) 42c must be set to the i-th value of the dependent variable y, 43c for i=1,2,...,m-1. the element y(m) is not used. 44c unchanged on exit. 45c w : real array of dimension at least (m). before entry, w(i) 46c must be set to the i-th value in the set of weights. the 47c w(i) must be strictly positive. w(m) is not used. 48c see also further comments. unchanged on exit. 49c k : integer. on entry k must specify the degree of the spline. 50c 1<=k<=5. it is recommended to use cubic splines (k=3). 51c the user is strongly dissuaded from choosing k even,together 52c with a small s-value. unchanged on exit. 53c s : real.on entry (in case iopt>=0) s must specify the smoothing 54c factor. s >=0. unchanged on exit. 55c for advice on the choice of s see further comments. 56c nest : integer. on entry nest must contain an over-estimate of the 57c total number of knots of the spline returned, to indicate 58c the storage space available to the routine. nest >=2*k+2. 59c in most practical situation nest=m/2 will be sufficient. 60c always large enough is nest=m+2*k,the number of knots needed 61c for interpolation (s=0). unchanged on exit. 62c n : integer. 63c unless ier = 10 (in case iopt >=0), n will contain the 64c total number of knots of the spline approximation returned. 65c if the computation mode iopt=1 is used this value of n 66c should be left unchanged between subsequent calls. 67c in case iopt=-1, the value of n must be specified on entry. 68c t : real array of dimension at least (nest). 69c on successful exit, this array will contain the knots of the 70c spline,i.e. the position of the interior knots t(k+2),t(k+3) 71c ...,t(n-k-1) as well as the position of the additional knots 72c t(1),t(2),...,t(k+1)=x(1) and t(n-k)=x(m),..,t(n) needed for 73c the b-spline representation. 74c if the computation mode iopt=1 is used, the values of t(1), 75c t(2),...,t(n) should be left unchanged between subsequent 76c calls. if the computation mode iopt=-1 is used, the values 77c t(k+2),...,t(n-k-1) must be supplied by the user, before 78c entry. see also the restrictions (ier=10). 79c c : real array of dimension at least (nest). 80c on successful exit, this array will contain the coefficients 81c c(1),c(2),..,c(n-k-1) in the b-spline representation of s(x) 82c fp : real. unless ier = 10, fp contains the weighted sum of 83c squared residuals of the spline approximation returned. 84c wrk : real array of dimension at least (m*(k+1)+nest*(8+5*k)). 85c used as working space. if the computation mode iopt=1 is 86c used, the values wrk(1),...,wrk(n) should be left unchanged 87c between subsequent calls. 88c lwrk : integer. on entry,lwrk must specify the actual dimension of 89c the array wrk as declared in the calling (sub)program. lwrk 90c must not be too small (see wrk). unchanged on exit. 91c iwrk : integer array of dimension at least (nest). 92c used as working space. if the computation mode iopt=1 is 93c used,the values iwrk(1),...,iwrk(n) should be left unchanged 94c between subsequent calls. 95c ier : integer. unless the routine detects an error, ier contains a 96c non-positive value on exit, i.e. 97c ier=0 : normal return. the spline returned has a residual sum of 98c squares fp such that abs(fp-s)/s <= tol with tol a relat- 99c ive tolerance set to 0.001 by the program. 100c ier=-1 : normal return. the spline returned is an interpolating 101c periodic spline (fp=0). 102c ier=-2 : normal return. the spline returned is the weighted least- 103c squares constant. in this extreme case fp gives the upper 104c bound fp0 for the smoothing factor s. 105c ier=1 : error. the required storage space exceeds the available 106c storage space, as specified by the parameter nest. 107c probably causes : nest too small. if nest is already 108c large (say nest > m/2), it may also indicate that s is 109c too small 110c the approximation returned is the least-squares periodic 111c spline according to the knots t(1),t(2),...,t(n). (n=nest) 112c the parameter fp gives the corresponding weighted sum of 113c squared residuals (fp>s). 114c ier=2 : error. a theoretically impossible result was found during 115c the iteration process for finding a smoothing spline with 116c fp = s. probably causes : s too small. 117c there is an approximation returned but the corresponding 118c weighted sum of squared residuals does not satisfy the 119c condition abs(fp-s)/s < tol. 120c ier=3 : error. the maximal number of iterations maxit (set to 20 121c by the program) allowed for finding a smoothing spline 122c with fp=s has been reached. probably causes : s too small 123c there is an approximation returned but the corresponding 124c weighted sum of squared residuals does not satisfy the 125c condition abs(fp-s)/s < tol. 126c ier=10 : error. on entry, the input data are controlled on validity 127c the following restrictions must be satisfied. 128c -1<=iopt<=1, 1<=k<=5, m>1, nest>2*k+2, w(i)>0,i=1,...,m-1 129c x(1)<x(2)<...<x(m), lwrk>=(k+1)*m+nest*(8+5*k) 130c if iopt=-1: 2*k+2<=n<=min(nest,m+2*k) 131c x(1)<t(k+2)<t(k+3)<...<t(n-k-1)<x(m) 132c the schoenberg-whitney conditions, i.e. there 133c must be a subset of data points xx(j) with 134c xx(j) = x(i) or x(i)+(x(m)-x(1)) such that 135c t(j) < xx(j) < t(j+k+1), j=k+1,...,n-k-1 136c if iopt>=0: s>=0 137c if s=0 : nest >= m+2*k 138c if one of these conditions is found to be violated,control 139c is immediately repassed to the calling program. in that 140c case there is no approximation returned. 141c 142c further comments: 143c by means of the parameter s, the user can control the tradeoff 144c between closeness of fit and smoothness of fit of the approximation. 145c if s is too large, the spline will be too smooth and signal will be 146c lost ; if s is too small the spline will pick up too much noise. in 147c the extreme cases the program will return an interpolating periodic 148c spline if s=0 and the weighted least-squares constant if s is very 149c large. between these extremes, a properly chosen s will result in 150c a good compromise between closeness of fit and smoothness of fit. 151c to decide whether an approximation, corresponding to a certain s is 152c satisfactory the user is highly recommended to inspect the fits 153c graphically. 154c recommended values for s depend on the weights w(i). if these are 155c taken as 1/d(i) with d(i) an estimate of the standard deviation of 156c y(i), a good s-value should be found in the range (m-sqrt(2*m),m+ 157c sqrt(2*m)). if nothing is known about the statistical error in y(i) 158c each w(i) can be set equal to one and s determined by trial and 159c error, taking account of the comments above. the best is then to 160c start with a very large value of s ( to determine the least-squares 161c constant and the corresponding upper bound fp0 for s) and then to 162c progressively decrease the value of s ( say by a factor 10 in the 163c beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the 164c approximation shows more detail) to obtain closer fits. 165c to economize the search for a good s-value the program provides with 166c different modes of computation. at the first call of the routine, or 167c whenever he wants to restart with the initial set of knots the user 168c must set iopt=0. 169c if iopt=1 the program will continue with the set of knots found at 170c the last call of the routine. this will save a lot of computation 171c time if percur is called repeatedly for different values of s. 172c the number of knots of the spline returned and their location will 173c depend on the value of s and on the complexity of the shape of the 174c function underlying the data. but, if the computation mode iopt=1 175c is used, the knots returned may also depend on the s-values at 176c previous calls (if these were smaller). therefore, if after a number 177c of trials with different s-values and iopt=1, the user can finally 178c accept a fit as satisfactory, it may be worthwhile for him to call 179c percur once more with the selected value for s but now with iopt=0. 180c indeed, percur may then return an approximation of the same quality 181c of fit but with fewer knots and therefore better if data reduction 182c is also an important objective for the user. 183c 184c other subroutines required: 185c fpbacp,fpbspl,fpchep,fpperi,fpdisc,fpgivs,fpknot,fprati,fprota 186c 187c references: 188c dierckx p. : algorithms for smoothing data with periodic and 189c parametric splines, computer graphics and image 190c processing 20 (1982) 171-184. 191c dierckx p. : algorithms for smoothing data with periodic and param- 192c etric splines, report tw55, dept. computer science, 193c k.u.leuven, 1981. 194c dierckx p. : curve and surface fitting with splines, monographs on 195c numerical analysis, oxford university press, 1993. 196c 197c author: 198c p.dierckx 199c dept. computer science, k.u. leuven 200c celestijnenlaan 200a, b-3001 heverlee, belgium. 201c e-mail : Paul.Dierckx@cs.kuleuven.ac.be 202c 203c creation date : may 1979 204c latest update : march 1987 205c 206c .. 207c ..scalar arguments.. 208 real*8 s,fp 209 integer iopt,m,k,nest,n,lwrk,ier 210c ..array arguments.. 211 real*8 x(m),y(m),w(m),t(nest),c(nest),wrk(lwrk) 212 integer iwrk(nest) 213c ..local scalars.. 214 real*8 per,tol 215 integer i,ia1,ia2,ib,ifp,ig1,ig2,iq,iz,i1,i2,j1,j2,k1,k2,lwest, 216 * maxit,m1,nmin 217c ..subroutine references.. 218c perper,pcheck 219c .. 220c we set up the parameters tol and maxit 221 maxit = 20 222 tol = 0.1e-02 223c before starting computations a data check is made. if the input data 224c are invalid, control is immediately repassed to the calling program. 225 ier = 10 226 if(k.le.0 .or. k.gt.5) go to 50 227 k1 = k+1 228 k2 = k1+1 229 if(iopt.lt.(-1) .or. iopt.gt.1) go to 50 230 nmin = 2*k1 231 if(m.lt.2 .or. nest.lt.nmin) go to 50 232 lwest = m*k1+nest*(8+5*k) 233 if(lwrk.lt.lwest) go to 50 234 m1 = m-1 235 do 10 i=1,m1 236 if(x(i).ge.x(i+1) .or. w(i).le.0.) go to 50 237 10 continue 238 if(iopt.ge.0) go to 30 239 if(n.le.nmin .or. n.gt.nest) go to 50 240 per = x(m)-x(1) 241 j1 = k1 242 t(j1) = x(1) 243 i1 = n-k 244 t(i1) = x(m) 245 j2 = j1 246 i2 = i1 247 do 20 i=1,k 248 i1 = i1+1 249 i2 = i2-1 250 j1 = j1+1 251 j2 = j2-1 252 t(j2) = t(i2)-per 253 t(i1) = t(j1)+per 254 20 continue 255 call fpchep(x,m,t,n,k,ier) 256 if (ier.eq.0) go to 40 257 go to 50 258 30 if(s.lt.0.) go to 50 259 if(s.eq.0. .and. nest.lt.(m+2*k)) go to 50 260 ier = 0 261c we partition the working space and determine the spline approximation. 262 40 ifp = 1 263 iz = ifp+nest 264 ia1 = iz+nest 265 ia2 = ia1+nest*k1 266 ib = ia2+nest*k 267 ig1 = ib+nest*k2 268 ig2 = ig1+nest*k2 269 iq = ig2+nest*k1 270 call fpperi(iopt,x,y,w,m,k,s,nest,tol,maxit,k1,k2,n,t,c,fp, 271 * wrk(ifp),wrk(iz),wrk(ia1),wrk(ia2),wrk(ib),wrk(ig1),wrk(ig2), 272 * wrk(iq),iwrk,ier) 273 50 return 274 end 275