## bvsp-spline

### Description

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)

### Installation notes

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:

- N
- the degree of the spline (must be integer >= 3)
- K
- the class of continuity of the spline (first or second derivative). K=1 or K=2 and N >= 3*K
- X
- SRFI-4
`f64vector`value containing the x coordinates of the data points (must be the same length as Y) - Y
- SRFI-4
`f64vector`value containing the y coordinates of the data points (must be the same length as X)

The optional arguments are:

- shape-constraint
- one of
`'none`,`'monotonicity`,`'convexity`,`'monotonicity+convexity`,`'local`. Default is`'none` - boundary-condition
- one of
`'none`,`'non-separable`,`'separable`. Default is`'none` - derivative-computation
- one of
`'order1`,`'order2`,`'order3`,`'classic`. Default is`'order2` - d0
- left separable boundary condition for the first derivative (only used when
`boundary-condition`is`'separable`) - dnp
- right separable boundary condition for the first derivative (only used when
`boundary-condition`is`'separable`) - d20
- left separable boundary condition for the second derivative (only used when
`boundary-condition`is`'separable`and K=2) - d2np
- right separable boundary condition for the second derivative (only used when
`boundary-condition`is`'separable`and K=2) - eps
- relative tolerance of the method. Default is 1e-4
- constr
- 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) - beta
- 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`) - betainv
- user-supplied procedure of the form
`(LAMBDA X)`, which is the inverse of`BETA`(only used when`boundary-condition`is`'non-separable`) - rho
- 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) - rhoinv
- 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) - kmax
- the number of iterations allowed for selecting the minimal set
`ASTAR`(described in the paper) - maxstp
- the number of iterations allowed for finding the set
`DSTAR`(described in the paper) - d
- SRFI-4
`f64vector`value containing the first derivatives at the points in X (only used when`derivative-computation`is`'classic`) - d2
- SRFI-4
`f64vector`value containing the second derivatives at the points in X (only used when`derivative-computation`is`'classic`and K=2)

*[procedure]*

`(evaluate n k x y d d2 xtab errc #!key (search-method 'binary) (derivatives 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`.

### Example

(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)) ))) ))

### Version History

- 1.0 Initial release

### License

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