The Chicken bvsp-spline library provides bindings to the routines in the BVSPIS Fortran library. BVSPIS (Boundary-Valued Shape-Preserving Interpolating Splines) is a package for computing and evaluating shape-preserving interpolating splines, and is described in
P. Costantini. Algorithm 770: BVSPIS -- a package for computing boundary-valued shape-preserving interpolating splines. ACM Trans. Math. Softw. 23, No.2, 252-254 (1997)
The Chicken bvsp-spline library must be compiled along with the BVSPIS Fortran source code, however the ACM licensing terms prevent the distribution of BVSPIS via the Chicken package system.
Therefore, the user must download and unpack BVSPIS, and set the BVSPIS_PATH environment variable to indicate the location where the sources can be found. chicken-install can then be invoked as usual.
The following commands perform the necessary operations when invoked from bash in a typical Linux distribution:
wget http://www.netlib.org/toms/770 mkdir bvspis && awk 'NR>4' 770 > bvspis/bvspis.sh && chmod u+x bvspis/bvspis.sh pushd bvspis && ./bvspis.sh && popd BVSPIS_PATH=$PWD/bvspis chicken-install bvsp-spline
Library procedures[procedure] (compute n k x y #!key (shape-constraint 'none) (boundary-condition 'none) (derivative-computation 'order2) (d0 #f) (dnp #f) (d20 #f) (d2np #f) (eps 1e-4) (constr #f) (beta #f) (betainv #f) (rho #f) (rhoinv #f) (kmax #f) (maxstp #f) (d #f) (d2 #f) )
Computes the coefficients of a shape-preserving spline, of continuity class C(k), k=1,2 , which interpolates a set of data points and, if required, satisfies additional boundary conditions.
The result of the routine is a list of the form (D D2 CONSTR ERRC DIAGN). D D2 ERRC provide the input parameters for evaluate, which evaluates the spline and its derivatives along a set of tabulation points. CONSTR is an SRFI-4 s32vector that contains computed constraint information.
The required arguments are:
- the degree of the spline (must be integer >= 3)
- the class of continuity of the spline (first or second derivative). K=1 or K=2 and N >= 3*K
- SRFI-4 f64vector value containing the x coordinates of the data points (must be the same length as Y)
- SRFI-4 f64vector value containing the y coordinates of the data points (must be the same length as X)
The optional arguments are:
- one of 'none, 'monotonicity, 'convexity, 'monotonicity+convexity, 'local. Default is 'none
- one of 'none, 'non-separable, 'separable. Default is 'none
- one of 'order1, 'order2, 'order3, 'classic. Default is 'order2
- left separable boundary condition for the first derivative (only used when boundary-condition is 'separable)
- right separable boundary condition for the first derivative (only used when boundary-condition is 'separable)
- left separable boundary condition for the second derivative (only used when boundary-condition is 'separable and K=2)
- right separable boundary condition for the second derivative (only used when boundary-condition is 'separable and K=2)
- relative tolerance of the method. Default is 1e-4
- if shape-constraint is 'local, this argument containts a s32vector value with the desired constraints on the shape for each subinterval. Each element can be one of 0,1,2,3 (none, monotonicity, convexity, monotonicity and convexity constraint)
- user-supplied procedure of the form (LAMBDA X), which represents non-separable boundary conditions for the first derivatives (only used when boundary-condition is 'non-separable)
- user-supplied procedure of the form (LAMBDA X), which is the inverse of BETA (only used when boundary-condition is 'non-separable)
- user-supplied procedure of the form (LAMBDA X), which represents non-separable boundary conditions for the second derivatives (only used when boundary-condition is 'non-separable and K=2)
- user-supplied procedure of the form (LAMBDA X), which is the inverse of RHO (only used when boundary-condition is 'non-separable and K=2)
- the number of iterations allowed for selecting the minimal set ASTAR (described in the paper)
- the number of iterations allowed for finding the set DSTAR (described in the paper)
- SRFI-4 f64vector value containing the first derivatives at the points in X (only used when derivative-computation is 'classic)
- SRFI-4 f64vector value containing the second derivatives at the points in X (only used when derivative-computation is 'classic and K=2)
Evaluates the given spline at points given by argument XTAB, which must be an SRFI-4 f64vector value. Arguments N K X Y have the same meaning as for the compute routine. Arguments D D2 ERRC are produced by compute.
(use srfi-1 srfi-4 bvsp-spline) ;; Input data (let ((x (f64vector .000000000E+00 .628318500E+00 .125663700E+01 .188495600E+01 .251327400E+01 .314159300E+01 .376991100E+01 .439823000E+01 .502654800E+01 .565486700E+01 .628318500E+01)) (y (f64vector .200000000E+00 .387785300E+00 .115105700E+01 .751056500E+00 .787785200E+00 -.200000100E+00 -.387785300E+00 -.115105700E+01 -.751056500E+00 -.787784900E+00 .200000200E+00))) ;; Set the degree and the class of continuity of the spline. (let ((n 3) (k 1)) (let-values (((d d2 constr errc diagn) (compute n k x y))) (assert (zero? errc)) (let* ((xpn 100) (dxp (/ (- (f64vector-ref x (- (f64vector-length x) 1)) (f64vector-ref x 0)) xpn)) (xp (list->f64vector (list-tabulate 100 (lambda (i) (* i dxp)))))) (let-values (((y0tab y1tab y2tab erre) (evaluate n k x y d d2 xp errc))) (print "erre = " erre) (for-each (lambda (i v) (printf "y0tab(~A) = ~A~%" (f64vector-ref xp i) (f64vector-ref y0tab i)) ) '(32 65 98) '(0.746843749779158 -0.806157453777544 -0.0795180362289622)) ))) ))
- 1.0 Initial release
Copyright 2011 Ivan Raikov. All rights reserved.
This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A full copy of the GPL license can be found at <http://www.gnu.org/licenses/>.