Outdated egg!

This is an egg for CHICKEN 4, the unsupported old release. You're almost certainly looking for the CHICKEN 5 version of this egg, if it exists.

If it does not exist, there may be equivalent functionality provided by another egg; have a look at the egg index. Otherwise, please consider porting this egg to the current version of CHICKEN.

  1. Outdated egg!
  2. bvsp-spline
    1. Description
    2. Installation notes
    3. Library procedures
    4. Example
    5. Version History
    6. License

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

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