Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/642.f
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/642.f	(revision 15442)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/642.f	(revision 15442)
@@ -0,0 +1,912 @@
+C     ALGORITHM 642 COLLECTED ALGORITHMS FROM ACM.
+C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.12, NO. 2,
+C     JUN., 1986, P. 150.
+C   SUBROUTINE NAME     - CUBGCV
+C
+C--------------------------------------------------------------------------
+C
+C   COMPUTER            - VAX/DOUBLE
+C
+C   AUTHOR              - M.F.HUTCHINSON
+C                         CSIRO DIVISION OF MATHEMATICS AND STATISTICS
+C                         P.O. BOX 1965
+C                         CANBERRA, ACT 2601
+C                         AUSTRALIA
+C
+C   LATEST REVISION     - 15 AUGUST 1985
+C
+C   PURPOSE             - CUBIC SPLINE DATA SMOOTHER
+C
+C   USAGE               - CALL CUBGCV (X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
+C
+C   ARGUMENTS    X      - VECTOR OF LENGTH N CONTAINING THE
+C                           ABSCISSAE OF THE N DATA POINTS
+C                           (X(I),F(I)) I=1..N. (INPUT) X
+C                           MUST BE ORDERED SO THAT
+C                           X(I) .LT. X(I+1).
+C                F      - VECTOR OF LENGTH N CONTAINING THE
+C                           ORDINATES (OR FUNCTION VALUES)
+C                           OF THE N DATA POINTS (INPUT).
+C                DF     - VECTOR OF LENGTH N. (INPUT/OUTPUT)
+C                           DF(I) IS THE RELATIVE STANDARD DEVIATION
+C                           OF THE ERROR ASSOCIATED WITH DATA POINT I.
+C                           EACH DF(I) MUST BE POSITIVE.  THE VALUES IN
+C                           DF ARE SCALED BY THE SUBROUTINE SO THAT
+C                           THEIR MEAN SQUARE VALUE IS 1, AND UNSCALED
+C                           AGAIN ON NORMAL EXIT.
+C                           THE MEAN SQUARE VALUE OF THE DF(I) IS RETURNED
+C                           IN WK(7) ON NORMAL EXIT.
+C                           IF THE ABSOLUTE STANDARD DEVIATIONS ARE KNOWN,
+C                           THESE SHOULD BE PROVIDED IN DF AND THE ERROR
+C                           VARIANCE PARAMETER VAR (SEE BELOW) SHOULD THEN
+C                           BE SET TO 1.
+C                           IF THE RELATIVE STANDARD DEVIATIONS ARE UNKNOWN,
+C                           SET EACH DF(I)=1.
+C                N      - NUMBER OF DATA POINTS (INPUT).
+C                           N MUST BE .GE. 3.
+C                Y,C    - SPLINE COEFFICIENTS. (OUTPUT) Y
+C                           IS A VECTOR OF LENGTH N. C IS
+C                           AN N-1 BY 3 MATRIX. THE VALUE
+C                           OF THE SPLINE APPROXIMATION AT T IS
+C                           S(T)=((C(I,3)*D+C(I,2))*D+C(I,1))*D+Y(I)
+C                           WHERE X(I).LE.T.LT.X(I+1) AND
+C                           D = T-X(I).
+C                IC     - ROW DIMENSION OF MATRIX C EXACTLY
+C                           AS SPECIFIED IN THE DIMENSION
+C                           STATEMENT IN THE CALLING PROGRAM. (INPUT)
+C                VAR    - ERROR VARIANCE. (INPUT/OUTPUT)
+C                           IF VAR IS NEGATIVE (I.E. UNKNOWN) THEN
+C                           THE SMOOTHING PARAMETER IS DETERMINED
+C                           BY MINIMIZING THE GENERALIZED CROSS VALIDATION
+C                           AND AN ESTIMATE OF THE ERROR VARIANCE IS
+C                           RETURNED IN VAR.
+C                           IF VAR IS NON-NEGATIVE (I.E. KNOWN) THEN THE
+C                           SMOOTHING PARAMETER IS DETERMINED TO MINIMIZE
+C                           AN ESTIMATE, WHICH DEPENDS ON VAR, OF THE TRUE
+C                           MEAN SQUARE ERROR, AND VAR IS UNCHANGED.
+C                           IN PARTICULAR, IF VAR IS ZERO, THEN AN
+C                           INTERPOLATING NATURAL CUBIC SPLINE IS CALCULATED.
+C                           VAR SHOULD BE SET TO 1 IF ABSOLUTE STANDARD
+C                           DEVIATIONS HAVE BEEN PROVIDED IN DF (SEE ABOVE).
+C                JOB    - JOB SELECTION PARAMETER. (INPUT)
+C                         JOB = 0 SHOULD BE SELECTED IF POINT STANDARD ERROR
+C                           ESTIMATES ARE NOT REQUIRED IN SE.
+C                         JOB = 1 SHOULD BE SELECTED IF POINT STANDARD ERROR
+C                           ESTIMATES ARE REQUIRED IN SE.
+C                SE     - VECTOR OF LENGTH N CONTAINING BAYESIAN STANDARD
+C                           ERROR ESTIMATES OF THE FITTED SPLINE VALUES IN Y.
+C                           SE IS NOT REFERENCED IF JOB=0. (OUTPUT)
+C                WK     - WORK VECTOR OF LENGTH 7*(N + 2). ON NORMAL EXIT THE
+C                           FIRST 7 VALUES OF WK ARE ASSIGNED AS FOLLOWS:-
+C
+C                           WK(1) = SMOOTHING PARAMETER (= RHO/(RHO + 1))
+C                           WK(2) = ESTIMATE OF THE NUMBER OF DEGREES OF
+C                                   FREEDOM OF THE RESIDUAL SUM OF SQUARES
+C                           WK(3) = GENERALIZED CROSS VALIDATION
+C                           WK(4) = MEAN SQUARE RESIDUAL
+C                           WK(5) = ESTIMATE OF THE TRUE MEAN SQUARE ERROR
+C                                   AT THE DATA POINTS
+C                           WK(6) = ESTIMATE OF THE ERROR VARIANCE
+C                           WK(7) = MEAN SQUARE VALUE OF THE DF(I)
+C
+C                           IF WK(1)=0 (RHO=0) AN INTERPOLATING NATURAL CUBIC
+C                           SPLINE HAS BEEN CALCULATED.
+C                           IF WK(1)=1 (RHO=INFINITE) A LEAST SQUARES
+C                           REGRESSION LINE HAS BEEN CALCULATED.
+C                           WK(2) IS AN ESTIMATE OF THE NUMBER OF DEGREES OF
+C                           FREEDOM OF THE RESIDUAL WHICH REDUCES TO THE
+C                           USUAL VALUE OF N-2 WHEN A LEAST SQUARES REGRESSION
+C                           LINE IS CALCULATED.
+C                           WK(3),WK(4),WK(5) ARE CALCULATED WITH THE DF(I)
+C                           SCALED TO HAVE MEAN SQUARE VALUE 1.  THE
+C                           UNSCALED VALUES OF WK(3),WK(4),WK(5) MAY BE
+C                           CALCULATED BY DIVIDING BY WK(7).
+C                           WK(6) COINCIDES WITH THE OUTPUT VALUE OF VAR IF
+C                           VAR IS NEGATIVE ON INPUT.  IT IS CALCULATED WITH
+C                           THE UNSCALED VALUES OF THE DF(I) TO FACILITATE
+C                           COMPARISONS WITH A PRIORI VARIANCE ESTIMATES.
+C
+C                IER    - ERROR PARAMETER. (OUTPUT)
+C                         TERMINAL ERROR
+C                           IER = 129, IC IS LESS THAN N-1.
+C                           IER = 130, N IS LESS THAN 3.
+C                           IER = 131, INPUT ABSCISSAE ARE NOT
+C                             ORDERED SO THAT X(I).LT.X(I+1).
+C                           IER = 132, DF(I) IS NOT POSITIVE FOR SOME I.
+C                           IER = 133, JOB IS NOT 0 OR 1.
+C
+C   PRECISION/HARDWARE  - DOUBLE
+C
+C   REQUIRED ROUTINES   - SPINT1,SPFIT1,SPCOF1,SPERR1
+C
+C   REMARKS      THE NUMBER OF ARITHMETIC OPERATIONS REQUIRED BY THE
+C                SUBROUTINE IS PROPORTIONAL TO N.  THE SUBROUTINE
+C                USES AN ALGORITHM DEVELOPED BY M.F. HUTCHINSON AND
+C                F.R. DE HOOG, 'SMOOTHING NOISY DATA WITH SPLINE
+C                FUNCTIONS', NUMER. MATH. (IN PRESS)
+C
+C-----------------------------------------------------------------------
+C
+      SUBROUTINE CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
+C
+C---SPECIFICATIONS FOR ARGUMENTS---
+      INTEGER N,IC,JOB,IER
+      DOUBLE PRECISION X(N),F(N),DF(N),Y(N),C(IC,3),SE(N),VAR,
+     .                 WK(0:N+1,7)
+C
+C---SPECIFICATIONS FOR LOCAL VARIABLES---
+      DOUBLE PRECISION DELTA,ERR,GF1,GF2,GF3,GF4,R1,R2,R3,R4,TAU,RATIO,
+     .                 AVH,AVDF,AVAR,ZERO,ONE,STAT(6),P,Q
+C
+      DATA RATIO/2.0D0/
+      DATA TAU/1.618033989D0/
+      DATA ZERO,ONE/0.0D0,1.0D0/
+C
+C---INITIALIZE---
+      IER = 133
+      IF (JOB.LT.0 .OR. JOB.GT.1) GO TO 140
+      CALL SPINT1(X,AVH,F,DF,AVDF,N,Y,C,IC,WK,WK(0,4),IER)
+      IF (IER.NE.0) GO TO 140
+      AVAR = VAR
+      IF (VAR.GT.ZERO) AVAR = VAR*AVDF*AVDF
+C
+C---CHECK FOR ZERO VARIANCE---
+      IF (VAR.NE.ZERO) GO TO 10
+      R1 = ZERO
+      GO TO 90
+C
+C---FIND LOCAL MINIMUM OF GCV OR THE EXPECTED MEAN SQUARE ERROR---
+   10 R1 = ONE
+      R2 = RATIO*R1
+      CALL SPFIT1(X,AVH,DF,N,R2,P,Q,GF2,AVAR,STAT,Y,C,IC,WK,WK(0,4),
+     .            WK(0,6),WK(0,7))
+   20 CALL SPFIT1(X,AVH,DF,N,R1,P,Q,GF1,AVAR,STAT,Y,C,IC,WK,WK(0,4),
+     .            WK(0,6),WK(0,7))
+      IF (GF1.GT.GF2) GO TO 30
+C
+C---EXIT IF P ZERO---
+      IF (P.LE.ZERO) GO TO 100
+      R2 = R1
+      GF2 = GF1
+      R1 = R1/RATIO
+      GO TO 20
+
+   30 R3 = RATIO*R2
+   40 CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4),
+     .            WK(0,6),WK(0,7))
+      IF (GF3.GT.GF2) GO TO 50
+C
+C---EXIT IF Q ZERO---
+      IF (Q.LE.ZERO) GO TO 100
+      R2 = R3
+      GF2 = GF3
+      R3 = RATIO*R3
+      GO TO 40
+
+   50 R2 = R3
+      GF2 = GF3
+      DELTA = (R2-R1)/TAU
+      R4 = R1 + DELTA
+      R3 = R2 - DELTA
+      CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4),
+     .            WK(0,6),WK(0,7))
+      CALL SPFIT1(X,AVH,DF,N,R4,P,Q,GF4,AVAR,STAT,Y,C,IC,WK,WK(0,4),
+     .            WK(0,6),WK(0,7))
+C
+C---GOLDEN SECTION SEARCH FOR LOCAL MINIMUM---
+   60 IF (GF3.GT.GF4) GO TO 70
+      R2 = R4
+      GF2 = GF4
+      R4 = R3
+      GF4 = GF3
+      DELTA = DELTA/TAU
+      R3 = R2 - DELTA
+      CALL SPFIT1(X,AVH,DF,N,R3,P,Q,GF3,AVAR,STAT,Y,C,IC,WK,WK(0,4),
+     .            WK(0,6),WK(0,7))
+      GO TO 80
+
+   70 R1 = R3
+      GF1 = GF3
+      R3 = R4
+      GF3 = GF4
+      DELTA = DELTA/TAU
+      R4 = R1 + DELTA
+      CALL SPFIT1(X,AVH,DF,N,R4,P,Q,GF4,AVAR,STAT,Y,C,IC,WK,WK(0,4),
+     .            WK(0,6),WK(0,7))
+   80 ERR = (R2-R1)/ (R1+R2)
+      IF (ERR*ERR+ONE.GT.ONE .AND. ERR.GT.1.0D-6) GO TO 60
+      R1 = (R1+R2)*0.5D0
+C
+C---CALCULATE SPLINE COEFFICIENTS---
+   90 CALL SPFIT1(X,AVH,DF,N,R1,P,Q,GF1,AVAR,STAT,Y,C,IC,WK,WK(0,4),
+     .            WK(0,6),WK(0,7))
+  100 CALL SPCOF1(X,AVH,F,DF,N,P,Q,Y,C,IC,WK(0,6),WK(0,7))
+C
+C---OPTIONALLY CALCULATE STANDARD ERROR ESTIMATES---
+      IF (VAR.GE.ZERO) GO TO 110
+      AVAR = STAT(6)
+      VAR = AVAR/ (AVDF*AVDF)
+  110 IF (JOB.EQ.1) CALL SPERR1(X,AVH,DF,N,WK,P,AVAR,SE)
+C
+C---UNSCALE DF---
+      DO 120 I = 1,N
+         DF(I) = DF(I)*AVDF
+  120 CONTINUE
+C
+C--PUT STATISTICS IN WK---
+      DO 130 I = 0,5
+         WK(I,1) = STAT(I+1)
+  130 CONTINUE
+      WK(5,1) = STAT(6)/ (AVDF*AVDF)
+      WK(6,1) = AVDF*AVDF
+      GO TO 150
+C
+C---CHECK FOR ERROR CONDITION---
+  140 CONTINUE
+C     IF (IER.NE.0) CONTINUE
+  150 RETURN
+      END
+      SUBROUTINE SPINT1(X,AVH,Y,DY,AVDY,N,A,C,IC,R,T,IER)
+C
+C INITIALIZES THE ARRAYS C, R AND T FOR ONE DIMENSIONAL CUBIC
+C SMOOTHING SPLINE FITTING BY SUBROUTINE SPFIT1.  THE VALUES
+C DF(I) ARE SCALED SO THAT THE SUM OF THEIR SQUARES IS N
+C AND THE AVERAGE OF THE DIFFERENCES X(I+1) - X(I) IS CALCULATED
+C IN AVH IN ORDER TO AVOID UNDERFLOW AND OVERFLOW PROBLEMS IN
+C SPFIT1.
+C
+C SUBROUTINE SETS IER IF ELEMENTS OF X ARE NON-INCREASING,
+C IF N IS LESS THAN 3, IF IC IS LESS THAN N-1 OR IF DY(I) IS
+C NOT POSITIVE FOR SOME I.
+C
+C---SPECIFICATIONS FOR ARGUMENTS---
+      INTEGER N,IC,IER
+      DOUBLE PRECISION X(N),Y(N),DY(N),A(N),C(IC,3),R(0:N+1,3),
+     .                 T(0:N+1,2),AVH,AVDY
+C
+C---SPECIFICATIONS FOR LOCAL VARIABLES---
+      INTEGER I
+      DOUBLE PRECISION E,F,G,H,ZERO
+      DATA ZERO/0.0D0/
+C
+C---INITIALIZATION AND INPUT CHECKING---
+      IER = 0
+      IF (N.LT.3) GO TO 60
+      IF (IC.LT.N-1) GO TO 70
+C
+C---GET AVERAGE X SPACING IN AVH---
+      G = ZERO
+      DO 10 I = 1,N - 1
+         H = X(I+1) - X(I)
+         IF (H.LE.ZERO) GO TO 80
+         G = G + H
+   10 CONTINUE
+      AVH = G/ (N-1)
+C
+C---SCALE RELATIVE WEIGHTS---
+      G = ZERO
+      DO 20 I = 1,N
+         IF (DY(I).LE.ZERO) GO TO 90
+         G = G + DY(I)*DY(I)
+   20 CONTINUE
+      AVDY = DSQRT(G/N)
+C
+      DO 30 I = 1,N
+         DY(I) = DY(I)/AVDY
+   30 CONTINUE
+C
+C---INITIALIZE H,F---
+      H = (X(2)-X(1))/AVH
+      F = (Y(2)-Y(1))/H
+C
+C---CALCULATE A,T,R---
+      DO 40 I = 2,N - 1
+         G = H
+         H = (X(I+1)-X(I))/AVH
+         E = F
+         F = (Y(I+1)-Y(I))/H
+         A(I) = F - E
+         T(I,1) = 2.0D0* (G+H)/3.0D0
+         T(I,2) = H/3.0D0
+         R(I,3) = DY(I-1)/G
+         R(I,1) = DY(I+1)/H
+         R(I,2) = -DY(I)/G - DY(I)/H
+   40 CONTINUE
+C
+C---CALCULATE C = R'*R---
+      R(N,2) = ZERO
+      R(N,3) = ZERO
+      R(N+1,3) = ZERO
+      DO 50 I = 2,N - 1
+         C(I,1) = R(I,1)*R(I,1) + R(I,2)*R(I,2) + R(I,3)*R(I,3)
+         C(I,2) = R(I,1)*R(I+1,2) + R(I,2)*R(I+1,3)
+         C(I,3) = R(I,1)*R(I+2,3)
+   50 CONTINUE
+      RETURN
+C
+C---ERROR CONDITIONS---
+   60 IER = 130
+      RETURN
+
+   70 IER = 129
+      RETURN
+
+   80 IER = 131
+      RETURN
+
+   90 IER = 132
+      RETURN
+      END
+      SUBROUTINE SPFIT1(X,AVH,DY,N,RHO,P,Q,FUN,VAR,STAT,A,C,IC,R,T,U,V)
+C
+C FITS A CUBIC SMOOTHING SPLINE TO DATA WITH RELATIVE
+C WEIGHTING DY FOR A GIVEN VALUE OF THE SMOOTHING PARAMETER
+C RHO USING AN ALGORITHM BASED ON THAT OF C.H. REINSCH (1967),
+C NUMER. MATH. 10, 177-183.
+C
+C THE TRACE OF THE INFLUENCE MATRIX IS CALCULATED USING AN
+C ALGORITHM DEVELOPED BY M.F.HUTCHINSON AND F.R.DE HOOG (NUMER.
+C MATH., IN PRESS), ENABLING THE GENERALIZED CROSS VALIDATION
+C AND RELATED STATISTICS TO BE CALCULATED IN ORDER N OPERATIONS.
+C
+C THE ARRAYS A, C, R AND T ARE ASSUMED TO HAVE BEEN INITIALIZED
+C BY THE SUBROUTINE SPINT1.  OVERFLOW AND UNDERFLOW PROBLEMS ARE
+C AVOIDED BY USING P=RHO/(1 + RHO) AND Q=1/(1 + RHO) INSTEAD OF
+C RHO AND BY SCALING THE DIFFERENCES X(I+1) - X(I) BY AVH.
+C
+C THE VALUES IN DF ARE ASSUMED TO HAVE BEEN SCALED SO THAT THE
+C SUM OF THEIR SQUARED VALUES IS N.  THE VALUE IN VAR, WHEN IT IS
+C NON-NEGATIVE, IS ASSUMED TO HAVE BEEN SCALED TO COMPENSATE FOR
+C THE SCALING OF THE VALUES IN DF.
+C
+C THE VALUE RETURNED IN FUN IS AN ESTIMATE OF THE TRUE MEAN SQUARE
+C WHEN VAR IS NON-NEGATIVE, AND IS THE GENERALIZED CROSS VALIDATION
+C WHEN VAR IS NEGATIVE.
+C
+C---SPECIFICATIONS FOR ARGUMENTS---
+      INTEGER IC,N
+      DOUBLE PRECISION X(N),DY(N),RHO,STAT(6),A(N),C(IC,3),R(0:N+1,3),
+     .                 T(0:N+1,2),U(0:N+1),V(0:N+1),FUN,VAR,AVH,P,Q
+C
+C---LOCAL VARIABLES---
+      INTEGER I
+      DOUBLE PRECISION E,F,G,H,ZERO,ONE,TWO,RHO1
+      DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/
+C
+C---USE P AND Q INSTEAD OF RHO TO PREVENT OVERFLOW OR UNDERFLOW---
+      RHO1 = ONE + RHO
+      P = RHO/RHO1
+      Q = ONE/RHO1
+      IF (RHO1.EQ.ONE) P = ZERO
+      IF (RHO1.EQ.RHO) Q = ZERO
+C
+C---RATIONAL CHOLESKY DECOMPOSITION OF P*C + Q*T---
+      F = ZERO
+      G = ZERO
+      H = ZERO
+      DO 10 I = 0,1
+         R(I,1) = ZERO
+   10 CONTINUE
+      DO 20 I = 2,N - 1
+         R(I-2,3) = G*R(I-2,1)
+         R(I-1,2) = F*R(I-1,1)
+         R(I,1) = ONE/ (P*C(I,1)+Q*T(I,1)-F*R(I-1,2)-G*R(I-2,3))
+         F = P*C(I,2) + Q*T(I,2) - H*R(I-1,2)
+         G = H
+         H = P*C(I,3)
+   20 CONTINUE
+C
+C---SOLVE FOR U---
+      U(0) = ZERO
+      U(1) = ZERO
+      DO 30 I = 2,N - 1
+         U(I) = A(I) - R(I-1,2)*U(I-1) - R(I-2,3)*U(I-2)
+   30 CONTINUE
+      U(N) = ZERO
+      U(N+1) = ZERO
+      DO 40 I = N - 1,2,-1
+         U(I) = R(I,1)*U(I) - R(I,2)*U(I+1) - R(I,3)*U(I+2)
+   40 CONTINUE
+C
+C---CALCULATE RESIDUAL VECTOR V---
+      E = ZERO
+      H = ZERO
+      DO 50 I = 1,N - 1
+         G = H
+         H = (U(I+1)-U(I))/ ((X(I+1)-X(I))/AVH)
+         V(I) = DY(I)* (H-G)
+         E = E + V(I)*V(I)
+   50 CONTINUE
+      V(N) = DY(N)* (-H)
+      E = E + V(N)*V(N)
+C
+C---CALCULATE UPPER THREE BANDS OF INVERSE MATRIX---
+      R(N,1) = ZERO
+      R(N,2) = ZERO
+      R(N+1,1) = ZERO
+      DO 60 I = N - 1,2,-1
+         G = R(I,2)
+         H = R(I,3)
+         R(I,2) = -G*R(I+1,1) - H*R(I+1,2)
+         R(I,3) = -G*R(I+1,2) - H*R(I+2,1)
+         R(I,1) = R(I,1) - G*R(I,2) - H*R(I,3)
+   60 CONTINUE
+C
+C---CALCULATE TRACE---
+      F = ZERO
+      G = ZERO
+      H = ZERO
+      DO 70 I = 2,N - 1
+         F = F + R(I,1)*C(I,1)
+         G = G + R(I,2)*C(I,2)
+         H = H + R(I,3)*C(I,3)
+   70 CONTINUE
+      F = F + TWO* (G+H)
+C
+C---CALCULATE STATISTICS---
+      STAT(1) = P
+      STAT(2) = F*P
+      STAT(3) = N*E/ (F*F)
+      STAT(4) = E*P*P/N
+      STAT(6) = E*P/F
+      IF (VAR.GE.ZERO) GO TO 80
+      STAT(5) = STAT(6) - STAT(4)
+      FUN = STAT(3)
+      GO TO 90
+
+   80 STAT(5) = DMAX1(STAT(4)-TWO*VAR*STAT(2)/N+VAR,ZERO)
+      FUN = STAT(5)
+   90 RETURN
+      END
+      SUBROUTINE SPERR1(X,AVH,DY,N,R,P,VAR,SE)
+C
+C CALCULATES BAYESIAN ESTIMATES OF THE STANDARD ERRORS OF THE FITTED
+C VALUES OF A CUBIC SMOOTHING SPLINE BY CALCULATING THE DIAGONAL ELEMENTS
+C OF THE INFLUENCE MATRIX.
+C
+C---SPECIFICATIONS FOR ARGUMENTS---
+      INTEGER N
+      DOUBLE PRECISION X(N),DY(N),R(0:N+1,3),SE(N),AVH,P,VAR
+C
+C---SPECIFICATIONS FOR LOCAL VARIABLES---
+      INTEGER I
+      DOUBLE PRECISION F,G,H,F1,G1,H1,ZERO,ONE
+      DATA ZERO,ONE/0.0D0,1.0D0/
+C
+C---INITIALIZE---
+      H = AVH/ (X(2)-X(1))
+      SE(1) = ONE - P*DY(1)*DY(1)*H*H*R(2,1)
+      R(1,1) = ZERO
+      R(1,2) = ZERO
+      R(1,3) = ZERO
+C
+C---CALCULATE DIAGONAL ELEMENTS---
+      DO 10 I = 2,N - 1
+         F = H
+         H = AVH/ (X(I+1)-X(I))
+         G = -F - H
+         F1 = F*R(I-1,1) + G*R(I-1,2) + H*R(I-1,3)
+         G1 = F*R(I-1,2) + G*R(I,1) + H*R(I,2)
+         H1 = F*R(I-1,3) + G*R(I,2) + H*R(I+1,1)
+         SE(I) = ONE - P*DY(I)*DY(I)* (F*F1+G*G1+H*H1)
+   10 CONTINUE
+      SE(N) = ONE - P*DY(N)*DY(N)*H*H*R(N-1,1)
+C
+C---CALCULATE STANDARD ERROR ESTIMATES---
+      DO 20 I = 1,N
+         SE(I) = DSQRT(DMAX1(SE(I)*VAR,ZERO))*DY(I)
+   20 CONTINUE
+      RETURN
+      END
+      SUBROUTINE SPCOF1(X,AVH,Y,DY,N,P,Q,A,C,IC,U,V)
+C
+C CALCULATES COEFFICIENTS OF A CUBIC SMOOTHING SPLINE FROM
+C PARAMETERS CALCULATED BY SUBROUTINE SPFIT1.
+C
+C---SPECIFICATIONS FOR ARGUMENTS---
+      INTEGER IC,N
+      DOUBLE PRECISION X(N),Y(N),DY(N),P,Q,A(N),C(IC,3),U(0:N+1),
+     .                 V(0:N+1),AVH
+C
+C---SPECIFICATIONS FOR LOCAL VARIABLES---
+      INTEGER I
+      DOUBLE PRECISION H,QH
+C
+C---CALCULATE A---
+      QH = Q/ (AVH*AVH)
+      DO 10 I = 1,N
+         A(I) = Y(I) - P*DY(I)*V(I)
+         U(I) = QH*U(I)
+   10 CONTINUE
+C
+C---CALCULATE C---
+      DO 20 I = 1,N - 1
+         H = X(I+1) - X(I)
+         C(I,3) = (U(I+1)-U(I))/ (3.0D0*H)
+         C(I,1) = (A(I+1)-A(I))/H - (H*C(I,3)+U(I))*H
+         C(I,2) = U(I)
+   20 CONTINUE
+      RETURN
+      END
+C CUBGCV TEST DRIVER
+C ------------------
+C
+C AUTHOR          - M.F.HUTCHINSON
+C                   CSIRO DIVISION OF WATER AND LAND RESOURCES
+C                   GPO BOX 1666
+C                   CANBERRA ACT 2601
+C                   AUSTRALIA
+C
+C LATEST REVISION - 7 AUGUST 1986
+C
+C COMPUTER        - VAX/DOUBLE
+C
+C USAGE           - MAIN PROGRAM
+C
+C REQUIRED ROUTINES - CUBGCV,SPINT1,SPFIT1,SPCOF1,SPERR1,GGRAND
+C
+C REMARKS   USES SUBROUTINE CUBGCV TO FIT A CUBIC SMOOTHING SPLINE
+C           TO 50 DATA POINTS WHICH ARE GENERATED BY ADDING A RANDOM
+C           VARIABLE WITH UNIFORM DENSITY IN THE INTERVAL [-0.3,0.3]
+C           TO 50 POINTS SAMPLED FROM THE CURVE  Y=SIN(3*PI*X/2).
+C           RANDOM DEVIATES IN THE INTERVAL [0,1] ARE GENERATED BY THE
+C           DOUBLE PRECISION FUNCTION GGRAND (SIMILAR TO IMSL FUNCTION
+C           GGUBFS).  THE ABSCISSAE ARE UNEQUALLY SPACED IN [0,1].
+C
+C           POINT STANDARD ERROR ESTIMATES ARE RETURNED IN SE BY
+C           SETTING JOB=1.  THE ERROR VARIANCE ESTIMATE IS RETURNED
+C           IN VAR.  IT COMPARES FAVOURABLY WITH THE TRUE VALUE OF 0.03.
+C           SUMMARY STATISTICS FROM THE ARRAY WK ARE WRITTEN TO
+C           UNIT 6.  DATA VALUES AND FITTED VALUES WITH ESTIMATED
+C           STANDARD ERRORS ARE ALSO WRITTEN TO UNIT 6.
+C
+      PARAMETER (N=50, IC=49)
+C
+      INTEGER            JOB,IER
+      DOUBLE PRECISION   X(N),F(N),Y(N),DF(N),C(IC,3),WK(7*(N+2)),
+     *                   VAR,SE(N)
+      DOUBLE PRECISION   GGRAND,DSEED
+C
+C---INITIALIZE---
+      DSEED=1.2345D4
+      JOB=1
+      VAR=-1.0D0
+C
+C---CALCULATE DATA POINTS---
+      DO 10 I=1,N
+      X(I)=(I - 0.5)/N + (2.0*GGRAND(DSEED) - 1.0)/(3.0*N)
+      F(I)=DSIN(4.71238*X(I)) + (2.0*GGRAND(DSEED) - 1.0)*0.3
+      DF(I)=1.0D0
+  10  CONTINUE
+C
+C---FIT CUBIC SPLINE---
+      CALL CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
+C
+C---WRITE OUT RESULTS---
+      WRITE(6,20)
+  20  FORMAT(' CUBGCV TEST DRIVER RESULTS:')
+      WRITE(6,30) IER,VAR,WK(3),WK(4),WK(2)
+  30  FORMAT(/' IER =',I4/' VAR =',F7.4/
+     *        ' GENERALIZED CROSS VALIDATION =',F7.4/
+     *        ' MEAN SQUARE RESIDUAL         =',F7.4/
+     *        ' RESIDUAL DEGREES OF FREEDOM  =',F7.2)
+      WRITE(6,40)
+  40  FORMAT(/' INPUT DATA',17X,'OUTPUT RESULTS'//
+     *         '   I    X(I)    F(I)',6X,'    Y(I)   SE(I)',
+     *          '      C(I,1)      C(I,2)      C(I,3)')
+      DO 60 I=1,N-1
+      WRITE(6,50) I,X(I),F(I),Y(I),SE(I),(C(I,J),J=1,3)
+  50  FORMAT(I4,2F8.4,6X,2F8.4,3E12.4)
+  60  CONTINUE
+      WRITE(6,50) N,X(N),F(N),Y(N),SE(N)
+      STOP
+      END
+      DOUBLE PRECISION FUNCTION GGRAND(DSEED)
+C
+C DOUBLE PRECISION UNIFORM RANDOM NUMBER GENERATOR
+C
+C CONSTANTS: A = 7**5
+C            B = 2**31 - 1
+C            C = 2**31
+C
+C REFERENCE: IMSL MANUAL, CHAPTER G - GENERATION AND TESTING OF
+C                                     RANDOM NUMBERS
+C
+C---SPECIFICATIONS FOR ARGUMENTS---
+      DOUBLE PRECISION DSEED
+C
+C---SPECIFICATIONS FOR LOCAL VARIABLES---
+      DOUBLE PRECISION A,B,C,S
+C
+      DATA A,B,C/16807.0D0, 2147483647.0D0, 2147483648.0D0/
+C
+      S=DSEED
+      S=DMOD(A*S, B)
+      GGRAND=S/C
+      DSEED=S
+      RETURN
+      END
+
+CUBGCV TEST DRIVER RESULTS:
+
+IER =   0
+VAR = 0.0279
+GENERALIZED CROSS VALIDATION = 0.0318
+MEAN SQUARE RESIDUAL         = 0.0246
+RESIDUAL DEGREES OF FREEDOM  =  43.97
+
+INPUT DATA                 OUTPUT RESULTS
+
+  I    X(I)    F(I)          Y(I)   SE(I)      C(I,1)      C(I,2)      C(I,3)
+  1  0.0046  0.2222        0.0342  0.1004  0.3630E+01  0.0000E+00  0.2542E+02
+  2  0.0360 -0.1098        0.1488  0.0750  0.3705E+01  0.2391E+01 -0.9537E+01
+  3  0.0435 -0.0658        0.1767  0.0707  0.3740E+01  0.2175E+01 -0.4233E+02
+  4  0.0735  0.3906        0.2900  0.0594  0.3756E+01 -0.1642E+01 -0.2872E+02
+  5  0.0955  0.6054        0.3714  0.0558  0.3642E+01 -0.3535E+01  0.2911E+01
+  6  0.1078  0.3034        0.4155  0.0549  0.3557E+01 -0.3428E+01 -0.1225E+02
+  7  0.1269  0.7386        0.4822  0.0544  0.3412E+01 -0.4131E+01  0.2242E+02
+  8  0.1565  0.4616        0.5800  0.0543  0.3227E+01 -0.2143E+01  0.6415E+01
+  9  0.1679  0.4315        0.6165  0.0543  0.3180E+01 -0.1923E+01 -0.1860E+02
+ 10  0.1869  0.5716        0.6762  0.0544  0.3087E+01 -0.2985E+01 -0.3274E+02
+ 11  0.2149  0.6736        0.7595  0.0542  0.2843E+01 -0.5733E+01 -0.4435E+02
+ 12  0.2356  0.7388        0.8155  0.0539  0.2549E+01 -0.8486E+01 -0.5472E+02
+ 13  0.2557  1.1953        0.8630  0.0537  0.2139E+01 -0.1180E+02 -0.9784E+01
+ 14  0.2674  1.0299        0.8864  0.0536  0.1860E+01 -0.1214E+02  0.9619E+01
+ 15  0.2902  0.7981        0.9225  0.0534  0.1322E+01 -0.1149E+02 -0.7202E+01
+ 16  0.3155  0.8973        0.9485  0.0532  0.7269E+00 -0.1203E+02 -0.1412E+02
+ 17  0.3364  1.2695        0.9583  0.0530  0.2040E+00 -0.1292E+02  0.2796E+02
+ 18  0.3557  0.7253        0.9577  0.0527 -0.2638E+00 -0.1130E+02 -0.3453E+01
+ 19  0.3756  1.2127        0.9479  0.0526 -0.7176E+00 -0.1151E+02  0.3235E+02
+ 20  0.3881  0.7304        0.9373  0.0525 -0.9889E+00 -0.1030E+02  0.4381E+01
+ 21  0.4126  0.9810        0.9069  0.0525 -0.1486E+01 -0.9977E+01  0.1440E+02
+ 22  0.4266  0.7117        0.8842  0.0525 -0.1756E+01 -0.9373E+01 -0.8925E+01
+ 23  0.4566  0.7203        0.8227  0.0524 -0.2344E+01 -0.1018E+02 -0.2278E+02
+ 24  0.4704  0.9242        0.7884  0.0524 -0.2637E+01 -0.1112E+02 -0.4419E+01
+ 25  0.4914  0.7345        0.7281  0.0523 -0.3110E+01 -0.1140E+02 -0.3562E+01
+ 26  0.5084  0.7378        0.6720  0.0524 -0.3500E+01 -0.1158E+02  0.5336E+01
+ 27  0.5277  0.7441        0.6002  0.0525 -0.3941E+01 -0.1127E+02  0.2479E+02
+ 28  0.5450  0.5612        0.5286  0.0527 -0.4310E+01 -0.9980E+01  0.2920E+02
+ 29  0.5641  0.5049        0.4429  0.0529 -0.4659E+01 -0.8309E+01  0.3758E+02
+ 30  0.5857  0.4725        0.3390  0.0531 -0.4964E+01 -0.5878E+01  0.5563E+02
+ 31  0.6159  0.1380        0.1850  0.0531 -0.5167E+01 -0.8307E+00  0.4928E+02
+ 32  0.6317  0.1412        0.1036  0.0531 -0.5157E+01  0.1499E+01  0.5437E+02
+ 33  0.6446 -0.1110        0.0371  0.0531 -0.5091E+01  0.3614E+01  0.3434E+02
+ 34  0.6707 -0.2605       -0.0927  0.0532 -0.4832E+01  0.6302E+01  0.1164E+02
+ 35  0.6853 -0.1284       -0.1619  0.0533 -0.4640E+01  0.6812E+01  0.1617E+02
+ 36  0.7064 -0.3452       -0.2564  0.0536 -0.4332E+01  0.7834E+01  0.4164E+01
+ 37  0.7310 -0.5527       -0.3582  0.0538 -0.3939E+01  0.8141E+01 -0.2214E+02
+ 38  0.7531 -0.3459       -0.4415  0.0540 -0.3611E+01  0.6674E+01 -0.9205E+01
+ 39  0.7686 -0.5902       -0.4961  0.0541 -0.3410E+01  0.6245E+01 -0.2193E+02
+ 40  0.7952 -0.7644       -0.5828  0.0541 -0.3125E+01  0.4494E+01 -0.4649E+02
+ 41  0.8087 -0.5392       -0.6242  0.0541 -0.3029E+01  0.2614E+01 -0.3499E+02
+ 42  0.8352 -0.4247       -0.7031  0.0539 -0.2964E+01 -0.1603E+00  0.2646E+01
+ 43  0.8501 -0.6327       -0.7476  0.0538 -0.2967E+01 -0.4132E-01  0.1817E+02
+ 44  0.8726 -0.9983       -0.8139  0.0538 -0.2942E+01  0.1180E+01 -0.6774E+01
+ 45  0.8874 -0.9082       -0.8574  0.0542 -0.2911E+01  0.8778E+00 -0.1364E+02
+ 46  0.9139 -0.8930       -0.9340  0.0566 -0.2893E+01 -0.2044E+00 -0.8094E+01
+ 47  0.9271 -1.0233       -0.9723  0.0593 -0.2903E+01 -0.5258E+00 -0.1498E+02
+ 48  0.9473 -0.8839       -1.0313  0.0665 -0.2942E+01 -0.1433E+01  0.4945E+01
+ 49  0.9652 -1.0172       -1.0843  0.0766 -0.2989E+01 -0.1168E+01  0.1401E+02
+ 50  0.9930 -1.2715       -1.1679  0.0998
+Documentation:
+C   COMPUTER            - VAX/DOUBLE
+C
+C   AUTHOR              - M.F.HUTCHINSON
+C                         CSIRO DIVISION OF MATHEMATICS AND STATISTICS
+C                         P.O. BOX 1965
+C                         CANBERRA, ACT 2601
+C                         AUSTRALIA
+C
+C   LATEST REVISION     - 15 AUGUST 1985
+C
+C   PURPOSE             - CUBIC SPLINE DATA SMOOTHER
+C
+C   USAGE               - CALL CUBGCV (X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
+C
+C   ARGUMENTS    X      - VECTOR OF LENGTH N CONTAINING THE
+C                           ABSCISSAE OF THE N DATA POINTS
+C                           (X(I),F(I)) I=1..N. (INPUT) X
+C                           MUST BE ORDERED SO THAT
+C                           X(I) .LT. X(I+1).
+C                F      - VECTOR OF LENGTH N CONTAINING THE
+C                           ORDINATES (OR FUNCTION VALUES)
+C                           OF THE N DATA POINTS (INPUT).
+C                DF     - VECTOR OF LENGTH N. (INPUT/OUTPUT)
+C                           DF(I) IS THE RELATIVE STANDARD DEVIATION
+C                           OF THE ERROR ASSOCIATED WITH DATA POINT I.
+C                           EACH DF(I) MUST BE POSITIVE.  THE VALUES IN
+C                           DF ARE SCALED BY THE SUBROUTINE SO THAT
+C                           THEIR MEAN SQUARE VALUE IS 1, AND UNSCALED
+C                           AGAIN ON NORMAL EXIT.
+C                           THE MEAN SQUARE VALUE OF THE DF(I) IS RETURNED
+C                           IN WK(7) ON NORMAL EXIT.
+C                           IF THE ABSOLUTE STANDARD DEVIATIONS ARE KNOWN,
+C                           THESE SHOULD BE PROVIDED IN DF AND THE ERROR
+C                           VARIANCE PARAMETER VAR (SEE BELOW) SHOULD THEN
+C                           BE SET TO 1.
+C                           IF THE RELATIVE STANDARD DEVIATIONS ARE UNKNOWN,
+C                           SET EACH DF(I)=1.
+C                N      - NUMBER OF DATA POINTS (INPUT).
+C                           N MUST BE .GE. 3.
+C                Y,C    - SPLINE COEFFICIENTS. (OUTPUT) Y
+C                           IS A VECTOR OF LENGTH N. C IS
+C                           AN N-1 BY 3 MATRIX. THE VALUE
+C                           OF THE SPLINE APPROXIMATION AT T IS
+C                           S(T)=((C(I,3)*D+C(I,2))*D+C(I,1))*D+Y(I)
+C                           WHERE X(I).LE.T.LT.X(I+1) AND
+C                           D = T-X(I).
+C                IC     - ROW DIMENSION OF MATRIX C EXACTLY
+C                           AS SPECIFIED IN THE DIMENSION
+C                           STATEMENT IN THE CALLING PROGRAM. (INPUT)
+C                VAR    - ERROR VARIANCE. (INPUT/OUTPUT)
+C                           IF VAR IS NEGATIVE (I.E. UNKNOWN) THEN
+C                           THE SMOOTHING PARAMETER IS DETERMINED
+C                           BY MINIMIZING THE GENERALIZED CROSS VALIDATION
+C                           AND AN ESTIMATE OF THE ERROR VARIANCE IS
+C                           RETURNED IN VAR.
+C                           IF VAR IS NON-NEGATIVE (I.E. KNOWN) THEN THE
+C                           SMOOTHING PARAMETER IS DETERMINED TO MINIMIZE
+C                           AN ESTIMATE, WHICH DEPENDS ON VAR, OF THE TRUE
+C                           MEAN SQUARE ERROR, AND VAR IS UNCHANGED.
+C                           IN PARTICULAR, IF VAR IS ZERO, THEN AN
+C                           INTERPOLATING NATURAL CUBIC SPLINE IS CALCULATED.
+C                           VAR SHOULD BE SET TO 1 IF ABSOLUTE STANDARD
+C                           DEVIATIONS HAVE BEEN PROVIDED IN DF (SEE ABOVE).
+C                JOB    - JOB SELECTION PARAMETER. (INPUT)
+C                         JOB = 0 SHOULD BE SELECTED IF POINT STANDARD ERROR
+C                           ESTIMATES ARE NOT REQUIRED IN SE.
+C                         JOB = 1 SHOULD BE SELECTED IF POINT STANDARD ERROR
+C                           ESTIMATES ARE REQUIRED IN SE.
+C                SE     - VECTOR OF LENGTH N CONTAINING BAYESIAN STANDARD
+C                           ERROR ESTIMATES OF THE FITTED SPLINE VALUES IN Y.
+C                           SE IS NOT REFERENCED IF JOB=0. (OUTPUT)
+C                WK     - WORK VECTOR OF LENGTH 7*(N + 2). ON NORMAL EXIT THE
+C                           FIRST 7 VALUES OF WK ARE ASSIGNED AS FOLLOWS:-
+C
+C                           WK(1) = SMOOTHING PARAMETER (= RHO/(RHO + 1))
+C                           WK(2) = ESTIMATE OF THE NUMBER OF DEGREES OF
+C                                   FREEDOM OF THE RESIDUAL SUM OF SQUARES
+C                           WK(3) = GENERALIZED CROSS VALIDATION
+C                           WK(4) = MEAN SQUARE RESIDUAL
+C                           WK(5) = ESTIMATE OF THE TRUE MEAN SQUARE ERROR
+C                                   AT THE DATA POINTS
+C                           WK(6) = ESTIMATE OF THE ERROR VARIANCE
+C                           WK(7) = MEAN SQUARE VALUE OF THE DF(I)
+C
+C                           IF WK(1)=0 (RHO=0) AN INTERPOLATING NATURAL CUBIC
+C                           SPLINE HAS BEEN CALCULATED.
+C                           IF WK(1)=1 (RHO=INFINITE) A LEAST SQUARES
+C                           REGRESSION LINE HAS BEEN CALCULATED.
+C                           WK(2) IS AN ESTIMATE OF THE NUMBER OF DEGREES OF
+C                           FREEDOM OF THE RESIDUAL WHICH REDUCES TO THE
+C                           USUAL VALUE OF N-2 WHEN A LEAST SQUARES REGRESSION
+C                           LINE IS CALCULATED.
+C                           WK(3),WK(4),WK(5) ARE CALCULATED WITH THE DF(I)
+C                           SCALED TO HAVE MEAN SQUARE VALUE 1.  THE
+C                           UNSCALED VALUES OF WK(3),WK(4),WK(5) MAY BE
+C                           CALCULATED BY DIVIDING BY WK(7).
+C                           WK(6) COINCIDES WITH THE OUTPUT VALUE OF VAR IF
+C                           VAR IS NEGATIVE ON INPUT.  IT IS CALCULATED WITH
+C                           THE UNSCALED VALUES OF THE DF(I) TO FACILITATE
+C                           COMPARISONS WITH A PRIORI VARIANCE ESTIMATES.
+C
+C                IER    - ERROR PARAMETER. (OUTPUT)
+C                         TERMINAL ERROR
+C                           IER = 129, IC IS LESS THAN N-1.
+C                           IER = 130, N IS LESS THAN 3.
+C                           IER = 131, INPUT ABSCISSAE ARE NOT
+C                             ORDERED SO THAT X(I).LT.X(I+1).
+C                           IER = 132, DF(I) IS NOT POSITIVE FOR SOME I.
+C                           IER = 133, JOB IS NOT 0 OR 1.
+C
+C   PRECISION/HARDWARE  - DOUBLE
+C
+C   REQUIRED ROUTINES   - SPINT1,SPFIT1,SPCOF1,SPERR1
+C
+C   REMARKS      THE NUMBER OF ARITHMETIC OPERATIONS REQUIRED BY THE
+C                SUBROUTINE IS PROPORTIONAL TO N.  THE SUBROUTINE
+C                USES AN ALGORITHM DEVELOPED BY M.F. HUTCHINSON AND
+C                F.R. DE HOOG, 'SMOOTHING NOISY DATA WITH SPLINE
+C                FUNCTIONS', NUMER. MATH. (IN PRESS)
+C
+C-----------------------------------------------------------------------
+C
+
+
+ALGORITHM
+
+CUBGCV calculates a natural cubic spline curve which smoothes a given set
+of data points, using statistical considerations to determine the amount
+of smoothing required, as described in reference 2.  If the error variance
+is known, it should be supplied to the routine in VAR.  The degree of
+smoothing is then determined by minimizing an unbiased estimate of the true
+mean square error.  On the other hand, if the error variance is not known, VAR
+should be set to -1.0.  The routine then determines the degree of smoothing
+by minimizing the generalized cross validation.  This is asymptotically the
+same as minimizing the true mean square error (see reference 1).  In this
+case, an estimate of the error variance is returned in VAR which may be
+compared with any a priori approximate estimates.  In either case, an
+estimate of the true mean square error is returned in WK(5).  This estimate,
+however, depends on the error variance estimate, and should only be accepted
+if the error variance estimate is reckoned to be correct.
+
+If job is set to 1, bayesian estimates of the standard error of each
+smoothed data value are returned in the array SE.  These also depend on
+the error variance estimate and should only be accepted if the error
+variance estimate is reckoned to be correct.  See reference 4.
+
+The number of arithmetic operations and the amount of storage required by
+the routine are both proportional to N, so that very large data sets may be
+analysed.  The data points do not have to be equally spaced or uniformly
+weighted.  The residual and the spline coefficients are calculated in the
+manner described in reference 3, while the trace and various statistics,
+including the generalized cross validation, are calculated in the manner
+described in reference 2.
+
+When VAR is known, any value of N greater than 2 is acceptable.  It is
+advisable, however, for N to be greater than about 20 when VAR is unknown.
+If the degree of smoothing done by CUBGCV when VAR is unknown is not
+satisfactory, the user should try specifying the degree of smoothing by
+setting VAR to a reasonable value.
+
+References:
+
+1.  Craven, Peter and Wahba, Grace, "Smoothing noisy data with spline
+    functions", Numer. Math. 31, 377-403 (1979).
+2.  Hutchinson, M.F. and de Hoog, F.R., "Smoothing noisy data with spline
+    functions", Numer. Math. (in press).
+3.  Reinsch, C.H., "Smoothing by spline functions", Numer. Math. 10,
+    177-183 (1967).
+4.  Wahba, Grace, "Bayesian 'confidence intervals' for the cross-validated
+    smoothing spline", J.R.Statist. Soc. B 45, 133-150 (1983).
+
+
+Example
+
+A sequence of 50 data points are generated by adding a random variable with
+uniform density in the interval [-0.3,0.3] to the curve y=sin(3*pi*x/2).
+The abscissae are unequally spaced in [0,1].  Point standard error estimates
+are returned in SE by setting JOB to 1.  The error variance estimate is
+returned in VAR.  It compares favourably with the true value of 0.03.
+The IMSL function GGUBFS is used to generate sample values of a uniform
+variable on [0,1].
+
+
+INPUT:
+
+      INTEGER          N,IC,JOB,IER
+      DOUBLE PRECISION X(50),F(50),Y(50),DF(50),C(49,3),WK(400),
+     *                 VAR,SE(50)
+      DOUBLE PRECISION GGUBFS,DSEED,DN
+      DATA DSEED/1.2345D4/
+C
+      N=50
+      IC=49
+      JOB=1
+      VAR=-1.0D0
+      DN=N
+C
+      DO 10 I=1,N
+      X(I)=(I - 0.5)/DN + (2.0*GGUBFS(DSEED) - 1.0)/(3.0*DN)
+      F(I)=DSIN(4.71238*X(I)) + (2.0*GGUBFS(DSEED) - 1.0)*0.3
+      DF(I)=1.0D0
+  10  CONTINUE
+      CALL CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
+       .
+       .
+       .
+      END
+
+OUTPUT:
+
+IER = 0
+VAR = 0.0279
+GENERALIZED CROSS VALIDATION = 0.0318
+MEAN SQUARE RESIDUAL = 0.0246
+RESIDUAL DEGREES OF FREEDOM = 43.97
+FOR CHECKING PURPOSES THE FOLLOWING OUTPUT IS GIVEN:
+
+X(1)  = 0.0046    F(1)  =  0.2222     Y(1)  =  0.0342     SE(1)  = 0.1004
+X(21) = 0.4126    F(21) =  0.9810     Y(21) =  0.9069     SE(21) = 0.0525
+X(41) = 0.8087    F(41) = -0.5392     Y(41) = -0.6242     SE(41) = 0.0541
+
Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs	(revision 15441)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/GAM.cs	(revision 15442)
@@ -218,9 +218,10 @@
           product = product.Zip(problemData.Dataset.GetDoubleValues(inputVars[i], problemData.TrainingIndices), (pi, vi) => pi * vi).ToArray();
         }
-        // Umständlich!
-        return Splines.CalculatePenalizedRegressionSpline(
+        double optTolerance, looRMSE;
+        return Splines.CalculateSmoothingSplineReinsch(
           product,
-          (double[])target.Clone(), lambda,
-          problemData.TargetVariable, inputVars
+          (double[])target.Clone(), inputVars,
+          problemData.TargetVariable, 
+          out optTolerance, out looRMSE
           );
       } else return new ConstantModel(target.Average(), problemData.TargetVariable);
Index: /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs	(revision 15441)
+++ /branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/Splines.cs	(revision 15442)
@@ -21,19 +21,14 @@
 
 using System;
-using System.Collections.Concurrent;
 using System.Collections.Generic;
 using System.Linq;
 using System.Threading;
-using System.Threading.Tasks;
 using HeuristicLab.Common;
 using HeuristicLab.Core;
 using HeuristicLab.Data;
-using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
 using HeuristicLab.Optimization;
 using HeuristicLab.Parameters;
 using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
 using HeuristicLab.Problems.DataAnalysis;
-using HeuristicLab.Problems.DataAnalysis.Symbolic;
-using HeuristicLab.Problems.DataAnalysis.Symbolic.Regression;
 using HeuristicLab.Random;
 using MathNet.Numerics.Interpolation;
@@ -71,5 +66,5 @@
           "Smoothing Spline (Mine)",
           "Smoothing Spline (Reinsch)",
-          "Smoothing Spline (Reinsch with automatic tolerance determination using CV)",
+          "Smoothing Spline (Reinsch with automatic tolerance determination using LOOCV)",
           "B-Spline Smoothing",
           "Penalized Regression Spline (alglib)"
@@ -150,7 +145,10 @@
               break;
             }
-          case "Smoothing Spline (Reinsch with automatic tolerance determination using CV)": {
-              var model = CalculateSmoothingSplineReinschWithAutomaticTolerance(x, y, inputVars, Problem.ProblemData.TargetVariable);
-              Results.Add(new Result("Solution", new RegressionSolution(model, (IRegressionProblemData)Problem.ProblemData.Clone())));
+          case "Smoothing Spline (Reinsch with automatic tolerance determination using LOOCV)": {
+              double optTol, looRMSE;
+              var model = CalculateSmoothingSplineReinsch(x, y, inputVars, Problem.ProblemData.TargetVariable, out optTol, out looRMSE);
+              Results.Add(new Result("Solution", model.CreateRegressionSolution((IRegressionProblemData)Problem.ProblemData.Clone())));
+              Results.Add(new Result("Optimal tolerance", new DoubleValue(optTol)));
+              Results.Add(new Result("RMSE (LOO-CV)", new DoubleValue(looRMSE)));
               break;
             }
@@ -159,6 +157,5 @@
               break;
             }
-          case "Penalized Regression Spline (alglib)":
-          {
+          case "Penalized Regression Spline (alglib)": {
               var lambda = ((IValueParameter<DoubleValue>)Parameters["Lambda"]).Value.Value; // controls extent of smoothing
               var model = CalculatePenalizedRegressionSpline(x, y, lambda, Problem.ProblemData.TargetVariable, inputVars);
@@ -168,5 +165,5 @@
 
               break;
-          }
+            }
           default: throw new NotSupportedException();
         }
@@ -176,7 +173,5 @@
 
 
-
-    public static IRegressionModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, string targetVar, string[] inputVars)
-    {
+    public static IRegressionModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, string targetVar, string[] inputVars) {
       int info;
       alglib.spline1dinterpolant interpolant;
@@ -189,29 +184,44 @@
     }
 
+    public static IRegressionEnsembleModel CalculatePenalizedRegressionSpline(double[] x, double[] y, double lambda, string targetVar, string[] inputVars,
+      out double avgTrainRMSE,
+      out double loocvRMSE) {
+      int info;
+      alglib.spline1dinterpolant interpolant;
+      alglib.spline1dfitreport rep;
+      int n = x.Length;
+      n = (int)Math.Max(50, 3 * Math.Sqrt(n));
+
+      double[] x_loo = new double[x.Length - 1];
+      double[] y_loo = new double[y.Length - 1];
+
+      var models = new List<IRegressionModel>();
+      var sumTrainRmse = 0.0;
+      var sse = 0.0;
+
+      for (int i = 0; i < x.Length; i++) {
+        Array.Copy(x, 0, x_loo, 0, i);
+        Array.Copy(x, i + 1, x_loo, i, x.Length - (i + 1));
+        Array.Copy(y, 0, y_loo, 0, i);
+        Array.Copy(y, i + 1, y_loo, i, y.Length - (i + 1));
+        alglib.spline1dfitpenalized(x_loo, y_loo, n, lambda, out info, out interpolant, out rep);
+        var y_pred = alglib.spline1dcalc(interpolant, x[i]);
+        double res = y[i] - y_pred;
+        sse += res * res;
+        sumTrainRmse += rep.rmserror;
+
+        var model = new Spline1dModel(interpolant, targetVar, inputVars);
+        models.Add(model);
+      }
+
+      loocvRMSE = Math.Sqrt(sse / x.Length);
+      avgTrainRMSE = sumTrainRmse / x.Length;
+
+      return new RegressionEnsembleModel(models);
+    }
+
     public void BSplineSmoothing(double[] x, double[] y, string[] inputVars) {
       Array.Sort(x, y);
-
       CubicSpline2(x, y, inputVars);
-
-      //  // see Elements of Statistical Learning Appendix: Computations for Splines
-      //  int M = 4; // order of spline (cubic)
-      //  int K = x.Length;
-      //  int n = x.Length;
-      // 
-      //  double[] t = new double[K + 2 * M];
-      //  for (int i = 0; i < M; i++) t[i] = x[0];
-      //  for (int i = K; i < K + 2 * M; i++) t[i] = x[K - 1];
-      // 
-      //  double[,] B = new double[n, n + 4]; // basis function matrix
-      //  double[,] W = new double[n + 4, n + 4]; // penalty matrix (W[j,k] = 
-      //  for (int i = 0; i < M; i++) B[i] = new double[K + 2 * M];
-      //  for (int j = 0; j < K; j++) {
-      //    for (int i = 0; i < K + 2M - 1; i++) {
-      //      if (t[i] <= x[j] && x[j] < t[i + 1])
-      //        B[1][i] = 1.0;
-      //    }
-      //  }
-
-
     }
 
@@ -312,56 +322,78 @@
     }
 
-    public void CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars) {
-      double s = xOrig.Length + 1;
+    public void CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars) {
+      double s = x.Length + 1;
       double stdTol = ((IValueParameter<DoubleValue>)Parameters["StdDev (noise)"]).Value.Value; // controls extent of smoothing
 
-      var model = CalculateSmoothingSplineReinsch(xOrig, yOrig, inputVars, stdTol, s, Problem.ProblemData.TargetVariable);
-
-      Results.Add(new Result("Solution", new RegressionSolution(model, (IRegressionProblemData)Problem.ProblemData.Clone())));
-
-
-    }
-
-
-
-    public static IRegressionModel CalculateSmoothingSplineReinschWithAutomaticTolerance(double[] xOrig, double[] yOrig, string[] inputVars, string targetVar) {
-      int n = xOrig.Length;
-      var rows = Enumerable.Range(0, n);
-      var rand = new FastRandom(1234);
-      var trainRows = rows.Shuffle(rand).Take((int)(n * 0.66)).ToArray();
-      var testRows = rows.Except(trainRows).ToArray();
-
-      var trainX = trainRows.Select(i => xOrig[i]).ToArray();
-      var trainY = trainRows.Select(i => yOrig[i]).ToArray();
-      var testX = testRows.Select(i => xOrig[i]).ToArray();
-      var testY = testRows.Select(i => yOrig[i]).ToArray();
-
-      var shuffledDs = new Dataset(inputVars.Concat(new string[] { targetVar }), new[] { trainX.Concat(testX).ToArray(), trainY.Concat(testY).ToArray() });
-      var shuffeledProblemData = new RegressionProblemData(shuffledDs, inputVars, targetVar);
-      shuffeledProblemData.TrainingPartition.Start = 0;
-      shuffeledProblemData.TrainingPartition.End = trainX.Length;
-      shuffeledProblemData.TestPartition.Start = trainX.Length;
-      shuffeledProblemData.TestPartition.End = shuffledDs.Rows + 1;
-
-      double curTol = trainY.StandardDeviation();
-      double prevTestRMSE = double.PositiveInfinity;
-      double testRMSE = double.PositiveInfinity;
-      IRegressionModel prevModel = null;
-      IRegressionModel model = null;
-      do {                                     
-        prevModel = model;
-        prevTestRMSE = testRMSE;
-        model = CalculateSmoothingSplineReinsch(trainX, trainY, inputVars, curTol, s: trainX.Length + 1, targetVar: targetVar);
-
-        var sol = model.CreateRegressionSolution(shuffeledProblemData);
-        var trainRMSE = sol.TrainingRootMeanSquaredError;
-        testRMSE = sol.TestRootMeanSquaredError;
-        curTol = Math.Max(curTol * 0.5, 1e-12 * trainY.StandardDeviation());
-      } while (testRMSE < prevTestRMSE);
-      return prevModel;
-    }
-
-
-    public static IRegressionModel CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars, double stdTol, double s, string targetVar) {
+      var model = CalculateSmoothingSplineReinsch(x, y, inputVars, stdTol, s, Problem.ProblemData.TargetVariable);
+
+      Results.Add(new Result("Solution", model.CreateRegressionSolution((IRegressionProblemData)Problem.ProblemData.Clone())));
+    }
+
+
+    public static IRegressionEnsembleModel CalculateSmoothingSplineReinsch(double[] x, double[] y, double tol, string targetVar, string[] inputVars,
+      out double avgTrainRMSE,
+      out double loocvRMSE) {
+      int info;
+
+      double[] x_loo = new double[x.Length - 1];
+      double[] y_loo = new double[y.Length - 1];
+
+      var models = new List<IRegressionModel>();
+      var sumTrainRmse = 0.0;
+      var sse = 0.0;
+
+      for (int i = 0; i < x.Length; i++) {
+        Array.Copy(x, i + 1, x_loo, i, x.Length - (i + 1));
+        Array.Copy(y, i + 1, y_loo, i, y.Length - (i + 1));
+
+        var model = CalculateSmoothingSplineReinsch(x_loo, y_loo, inputVars, tol, x_loo.Length + 1, targetVar);
+
+        var y_pred = model.GetEstimatedValue(x[i]);
+        double res = y[i] - y_pred;
+        sse += res * res;
+
+        models.Add(model);
+
+        // add current row for next iteration
+        if (i < x_loo.Length) {
+          x_loo[i] = x[i];
+          y_loo[i] = y[i];
+        }
+      }
+
+      loocvRMSE = Math.Sqrt(sse / x.Length);
+      avgTrainRMSE = sumTrainRmse / x.Length;
+
+      return new RegressionEnsembleModel(models);
+    }
+
+
+    // automatically determines tolerance using loo cv
+    public static IRegressionModel CalculateSmoothingSplineReinsch(double[] x, double[] y, string[] inputVars, string targetVar,
+      out double optTolerance, out double looRMSE) {
+      var maxTolerance = y.StandardDeviation();
+      var minTolerance = maxTolerance * 1E-6;
+      optTolerance = maxTolerance;
+      var tol = maxTolerance;
+      looRMSE = double.PositiveInfinity;
+      IRegressionModel bestModel = new ConstantModel(y.Average(), targetVar);
+      while (tol > minTolerance) {
+        double loo_rmse;
+        double avgTrainRmse;
+        var model = Splines.CalculateSmoothingSplineReinsch(
+          x, y, tol, targetVar, inputVars, out avgTrainRmse, out loo_rmse);
+
+        if (loo_rmse < looRMSE) {
+          looRMSE = loo_rmse;
+          optTolerance = tol;
+          bestModel = model;
+        }
+        tol *= 0.85;
+      }
+      return bestModel;
+    }
+
+    public static ReinschSmoothingSplineModel CalculateSmoothingSplineReinsch(double[] xOrig, double[] yOrig, string[] inputVars, double stdTol, double s, string targetVar) {
       var minX = xOrig.Min();
       var maxX = xOrig.Max();
@@ -605,5 +637,5 @@
 
   // UNFINISHED
-  internal class ReinschSmoothingSplineModel : NamedItem, IRegressionModel {
+  public class ReinschSmoothingSplineModel : NamedItem, IRegressionModel {
     private MyArray<double> a;
     private MyArray<double> b;
@@ -631,5 +663,5 @@
       this.offset = orig.offset;
     }
-    public ReinschSmoothingSplineModel(
+    internal ReinschSmoothingSplineModel(
       MyArray<double> a,
       MyArray<double> b,
@@ -650,5 +682,6 @@
       // extrapolate for xx > x[n2]
       b[b.Length] = b[b.Length - 1];
-      d[b.Length] = d[d.Length - 1];
+      d[1] = 0;
+      // d[b.Length] = -d[d.Length - 1];
     }
 
@@ -661,31 +694,39 @@
     }
 
+    public double GetEstimatedValue(double xx) {
+      int n = x.Length;
+      if (xx <= x[1]) {
+        double h = xx - x[1];
+        return a[1] + h * (b[1] + h * (c[1] + h * d[1]));
+      } else if (xx >= x[n]) {
+        double h = xx - x[n];
+        return a[n] + h * (b[n] + h * (c[n] + h * d[n]));
+      } else {
+        // binary search
+        int lower = 1;
+        int upper = n;
+        while (true) {
+          if (upper < lower) throw new InvalidProgramException();
+          int i = lower + (upper - lower) / 2;
+          if (x[i] <= xx && xx < x[i + 1]) {
+            double h = xx - x[i];
+            return a[i] + h * (b[i] + h * (c[i] + h * d[i]));
+          } else if (xx < x[i]) {
+            upper = i - 1;
+          } else {
+            lower = i + 1;
+          }
+        }
+      }
+    }
+
     public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
       int n = x.Length;
-      foreach (var xx in dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows).Select(xi => (xi - offset) * scale)) {
-        if (xx <= x[1]) {
-          double h = xx - x[1];
-          yield return a[1] + h * (b[1] + h * (c[1] + h * d[1]));
-        } else if (xx >= x[n]) {
-          double h = xx - x[n];
-          yield return a[n] + h * (b[n] + h * (c[n] + h * d[n]));
-        } else {
-          // binary search
-          int lower = 1;
-          int upper = n;
-          while (true) {
-            if (upper < lower) throw new InvalidProgramException();
-            int i = lower + (upper - lower) / 2;
-            if (x[i] <= xx && xx < x[i + 1]) {
-              double h = xx - x[i];
-              yield return a[i] + h * (b[i] + h * (c[i] + h * d[i]));
-              break;
-            } else if (xx < x[i]) {
-              upper = i - 1;
-            } else {
-              lower = i + 1;
-            }
-          }
-        }
+      var product = dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows).ToArray();
+      foreach (var factor in VariablesUsedForPrediction.Skip(1)) {
+        product = product.Zip(dataset.GetDoubleValues(factor, rows), (pi, fi) => pi * fi).ToArray();
+      }
+      foreach (var xx in product.Select(xi => (xi - offset) * scale)) {
+        yield return GetEstimatedValue(xx);
       }
     }
@@ -764,5 +805,5 @@
     public IEnumerable<double> GetEstimatedValues(IDataset dataset, IEnumerable<int> rows) {
       var product = dataset.GetDoubleValues(VariablesUsedForPrediction.First(), rows).ToArray();
-      foreach(var factor in VariablesUsedForPrediction.Skip(1)) {
+      foreach (var factor in VariablesUsedForPrediction.Skip(1)) {
         product = product.Zip(dataset.GetDoubleValues(factor, rows), (pi, fi) => pi * fi).ToArray();
       }
Index: /branches/MathNetNumerics-Exploration-2789/Test/Test.csproj
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/Test/Test.csproj	(revision 15441)
+++ /branches/MathNetNumerics-Exploration-2789/Test/Test.csproj	(revision 15442)
@@ -36,4 +36,9 @@
   </PropertyGroup>
   <ItemGroup>
+    <Reference Include="ALGLIB-3.7.0, Version=3.7.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec, processorArchitecture=MSIL">
+      <SpecificVersion>False</SpecificVersion>
+      <HintPath>..\..\..\trunk\sources\bin\ALGLIB-3.7.0.dll</HintPath>
+    </Reference>
+    <Reference Include="HeuristicLab.Algorithms.DataAnalysis-3.4, Version=3.4.0.0, Culture=neutral, PublicKeyToken=ba48961d6f65dcec" />
     <Reference Include="HeuristicLab.Algorithms.GeneticAlgorithm-3.3">
       <HintPath>..\..\..\trunk\sources\bin\HeuristicLab.Algorithms.GeneticAlgorithm-3.3.dll</HintPath>
Index: /branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs
===================================================================
--- /branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs	(revision 15441)
+++ /branches/MathNetNumerics-Exploration-2789/Test/UnitTest1.cs	(revision 15442)
@@ -1,5 +1,9 @@
 ﻿using System;
+using System.Diagnostics;
+using System.Linq;
 using HeuristicLab.Algorithms.DataAnalysis.Experimental;
 using HeuristicLab.Algorithms.GeneticAlgorithm;
+using HeuristicLab.Common;
+using HeuristicLab.Problems.DataAnalysis;
 using HeuristicLab.SequentialEngine;
 using Microsoft.VisualStudio.TestTools.UnitTesting;
@@ -22,3 +26,90 @@
     }
   }
+
+  [TestClass]
+  public class PenalizedRegressionSplineTests {
+    [TestMethod]
+    public void TestPenalizedRegressionSplinesCrossValidation() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToArray();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToArray();
+
+      double bestRho = -15;
+      double best_loo_rmse = double.PositiveInfinity;
+      var clock = new Stopwatch();
+      clock.Start();
+      int iters = 0;
+      for (int rho = -15; rho <= 15; rho++) {
+        double loo_rmse;
+        double avgTrainRmse;
+        Splines.CalculatePenalizedRegressionSpline(xs, ys_noise, rho, "y", new string[] { "x" }, out avgTrainRmse, out loo_rmse);
+        iters++;
+        Console.WriteLine("{0} {1} {2}", rho, avgTrainRmse, loo_rmse);
+        if (loo_rmse < best_loo_rmse) {
+          best_loo_rmse = loo_rmse;
+          bestRho = rho;
+        }
+      }
+      clock.Stop();
+      Console.WriteLine("Best rho {0}, RMSE (LOO): {1}, ms/run {2}", bestRho, best_loo_rmse, clock.ElapsedMilliseconds / (double)iters);
+    }
+
+    [TestMethod]
+    public void TestReinschSplineCrossValidation() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToList();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToList();
+
+      var tol = ys_noise.StandardDeviation();
+
+      var ds = new Dataset(new string[] { "x", "y" }, new[] { xs, ys_noise });
+      var rows = Enumerable.Range(0, xs.Count);
+
+      var bestTol = double.PositiveInfinity;
+      var bestLooRmse = double.PositiveInfinity;
+      var clock = new Stopwatch();
+      clock.Start();
+      var iters = 0;
+      while (tol > 0.0001) {
+        double loo_rmse;
+        double avgTrainRmse;
+        var model = Splines.CalculateSmoothingSplineReinsch(
+          xs.ToArray(), ys_noise.ToArray(), tol, "y", new string[] { "x" }, out avgTrainRmse, out loo_rmse);
+        var y_pred = model.GetEstimatedValues(ds, rows);
+
+        Console.WriteLine("{0} {1} {2}", tol, avgTrainRmse, loo_rmse);
+        if (loo_rmse < bestLooRmse) {
+          bestLooRmse = loo_rmse;
+          bestTol = tol;
+        }
+        tol *= 0.8;
+        iters++;
+      }
+      clock.Stop();
+      Console.WriteLine("Best tolerance {0}, RMSE (LOO): {1}, ms/run {2}", bestTol, bestLooRmse, clock.ElapsedMilliseconds / (double)iters);
+    }
+
+    [TestMethod]
+    public void TestReinschSplineAutomaticTolerance() {
+      var xs = HeuristicLab.Common.SequenceGenerator.GenerateSteps(-3.5, 3.5, 0.02, includeEnd: true).ToList();
+      var ys = xs.Select(xi => alglib.normaldistr.normaldistribution(xi)); // 1.0 / (Math.Sqrt(2 * Math.PI) * Math.Exp(-0.5 * xi * xi))).ToArray();
+
+      alglib.hqrndstate state;
+      alglib.hqrndseed(1234, 5678, out state);
+      var ys_noise = ys.Select(yi => yi + alglib.hqrndnormal(state) * 0.01).ToList();
+
+      double optTol;
+      double looRMSE;
+      Splines.CalculateSmoothingSplineReinsch(xs.ToArray(), ys_noise.ToArray(), new string[] { "x" }, "y", out optTol, out looRMSE);
+
+      Console.WriteLine("Best tolerance {0}, RMSE (LOO): {1}", optTol, looRMSE);
+
+    }
+  }
 }
