testsplineinterpolationunit.pas

maths lib with source

unit testsplineinterpolationunit;
interface
uses Math, Ap, Sysutils, spline3;

function TestSplineInterpolation(Silent : Boolean):Boolean;
function testsplineinterpolationunit_test_silent():Boolean;
function testsplineinterpolationunit_test():Boolean;

implementation

procedure LConst(A : Double;
     B : Double;
     const C : TReal1DArray;
     LStep : Double;
     var L0 : Double;
     var L1 : Double;
     var L2 : Double);forward;
function TestUnpack(const C : TReal1DArray;
     const X : TReal1DArray):Boolean;forward;


function TestSplineInterpolation(Silent : Boolean):Boolean;
var
    WasErrors : Boolean;
    CSErrors : Boolean;
    HSErrors : Boolean;
    ASErrors : Boolean;
    LSErrors : Boolean;
    DSErrors : Boolean;
    UPErrors : Boolean;
    CPErrors : Boolean;
    LTErrors : Boolean;
    IErrors : Boolean;
    N : Integer;
    I : Integer;
    K : Integer;
    Pass : Integer;
    PassCount : Integer;
    BLType : Integer;
    BRType : Integer;
    X : TReal1DArray;
    Y : TReal1DArray;
    Y2 : TReal1DArray;
    D : TReal1DArray;
    C : TReal1DArray;
    C2 : TReal1DArray;
    A : Double;
    B : Double;
    BL : Double;
    BR : Double;
    T : Double;
    SA : Double;
    SB : Double;
    V : Double;
    LStep : Double;
    H : Double;
    L10 : Double;
    L11 : Double;
    L12 : Double;
    L20 : Double;
    L21 : Double;
    L22 : Double;
    S : Double;
    DS : Double;
    D2S : Double;
    S2 : Double;
    DS2 : Double;
    D2S2 : Double;
    VL : Double;
    VM : Double;
    VR : Double;
    Err : Double;
begin
    WasErrors := False;
    PassCount := 20;
    LStep := 0.005;
    H := 0.00001;
    LSErrors := False;
    CSErrors := False;
    HSErrors := False;
    ASErrors := False;
    DSErrors := False;
    CPErrors := False;
    UPErrors := False;
    LTErrors := False;
    IErrors := False;
    
    //
    // General test: linear, cubic, Hermite, Akima
    //
    N:=2;
    while N<=8 do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        SetLength(D, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Prepare task
            //
            A := -1-RandomReal;
            B := +1+RandomReal;
            BL := 2*RandomReal-1;
            BR := 2*RandomReal-1;
            I:=0;
            while I<=N-1 do
            begin
                X[I] := 0.5*(B+A)+0.5*(B-A)*Cos(PI*(2*i+1)/(2*n));
                if I=0 then
                begin
                    X[I] := A;
                end;
                if I=N-1 then
                begin
                    X[I] := B;
                end;
                Y[I] := Cos(1.3*Pi*X[I]+0.4);
                D[I] := -1.3*Pi*Sin(1.3*Pi*X[I]+0.4);
                Inc(I);
            end;
            I:=0;
            while I<=N-1 do
            begin
                K := RandomInteger(N);
                if K<>I then
                begin
                    T := X[I];
                    X[I] := X[K];
                    X[K] := T;
                    T := Y[I];
                    Y[I] := Y[K];
                    Y[K] := T;
                    T := D[I];
                    D[I] := D[K];
                    D[K] := T;
                end;
                Inc(I);
            end;
            
            //
            // Build linear spline
            // Test for general interpolation scheme properties:
            // * values at nodes
            // * continuous function
            // Test for specific properties is implemented below.
            //
            BuildLinearSpline(X, Y, N, C);
            Err := 0;
            I:=0;
            while I<=N-1 do
            begin
                Err := Max(Err, AbsReal(Y[I]-SplineInterpolation(C, X[I])));
                Inc(I);
            end;
            LSErrors := LSErrors or (Err>100*MachineEpsilon);
            LConst(A, B, C, LStep, L10, L11, L12);
            LConst(A, B, C, LStep/3, L20, L21, L22);
            LSErrors := LSErrors or (L20/L10>1.2);
            
            //
            // Build cubic spline.
            // Test for interpolation scheme properties:
            // * values at nodes
            // * boundary conditions
            // * continuous function
            // * continuous first derivative
            // * continuous second derivative
            //
            BLType:=0;
            while BLType<=2 do
            begin
                BRType:=0;
                while BRType<=2 do
                begin
                    BuildCubicSpline(X, Y, N, BLType, BL, BRType, BR, C);
                    Err := 0;
                    I:=0;
                    while I<=N-1 do
                    begin
                        Err := Max(Err, AbsReal(Y[I]-SplineInterpolation(C, X[I])));
                        Inc(I);
                    end;
                    CSErrors := CSErrors or (Err>100*MachineEpsilon);
                    Err := 0;
                    if BLType=0 then
                    begin
                        SplineDifferentiation(C, A-H, S, DS, D2S);
                        SplineDifferentiation(C, A+H, S2, DS2, D2S2);
                        T := (D2S2-D2S)/(2*H);
                        Err := Max(Err, AbsReal(T));
                    end;
                    if BLType=1 then
                    begin
                        T := (SplineInterpolation(C, A+H)-SplineInterpolation(C, A-H))/(2*H);
                        Err := Max(Err, AbsReal(BL-T));
                    end;
                    if BLType=2 then
                    begin
                        T := (SplineInterpolation(C, A+H)-2*SplineInterpolation(C, A)+SplineInterpolation(C, A-H))/Sqr(H);
                        Err := Max(Err, AbsReal(BL-T));
                    end;
                    if BRType=0 then
                    begin
                        SplineDifferentiation(C, B-H, S, DS, D2S);
                        SplineDifferentiation(C, B+H, S2, DS2, D2S2);
                        T := (D2S2-D2S)/(2*H);
                        Err := Max(Err, AbsReal(T));
                    end;
                    if BRType=1 then
                    begin
                        T := (SplineInterpolation(C, B+H)-SplineInterpolation(C, B-H))/(2*H);
                        Err := Max(Err, AbsReal(BR-T));
                    end;
                    if BRType=2 then
                    begin
                        T := (SplineInterpolation(C, B+H)-2*SplineInterpolation(C, B)+SplineInterpolation(C, B-H))/Sqr(H);
                        Err := Max(Err, AbsReal(BR-T));
                    end;
                    CSErrors := CSErrors or (Err>1.0E-3);
                    LConst(A, B, C, LStep, L10, L11, L12);
                    LConst(A, B, C, LStep/3, L20, L21, L22);
                    CSErrors := CSErrors or (L20/L10>1.2) and (L10>1.0E-6);
                    CSErrors := CSErrors or (L21/L11>1.2) and (L11>1.0E-6);
                    CSErrors := CSErrors or (L22/L12>1.2) and (L12>1.0E-6);
                    Inc(BRType);
                end;
                Inc(BLType);
            end;
            
            //
            // Build Hermite spline.
            // Test for interpolation scheme properties:
            // * values and derivatives at nodes
            // * continuous function
            // * continuous first derivative
            //
            BuildHermiteSpline(X, Y, D, N, C);
            Err := 0;
            I:=0;
            while I<=N-1 do
            begin
                Err := Max(Err, AbsReal(Y[I]-SplineInterpolation(C, X[I])));
                Inc(I);
            end;
            HSErrors := HSErrors or (Err>100*MachineEpsilon);
            Err := 0;
            I:=0;
            while I<=N-1 do
            begin
                T := (SplineInterpolation(C, X[I]+H)-SplineInterpolation(C, X[I]-H))/(2*H);
                Err := Max(Err, AbsReal(D[I]-T));
                Inc(I);
            end;
            HSErrors := HSErrors or (Err>1.0E-3);
            LConst(A, B, C, LStep, L10, L11, L12);
            LConst(A, B, C, LStep/3, L20, L21, L22);
            HSErrors := HSErrors or (L20/L10>1.2);
            HSErrors := HSErrors or (L21/L11>1.2);
            
            //
            // Build Akima spline
            // Test for general interpolation scheme properties:
            // * values at nodes
            // * continuous function
            // * continuous first derivative
            // Test for specific properties is implemented below.
            //
            if N>=5 then
            begin
                BuildAkimaSpline(X, Y, N, C);
                Err := 0;
                I:=0;
                while I<=N-1 do
                begin
                    Err := Max(Err, AbsReal(Y[I]-SplineInterpolation(C, X[I])));
                    Inc(I);
                end;
                ASErrors := ASErrors or (Err>100*MachineEpsilon);
                LConst(A, B, C, LStep, L10, L11, L12);
                LConst(A, B, C, LStep/3, L20, L21, L22);
                HSErrors := HSErrors or (L20/L10>1.2);
                HSErrors := HSErrors or (L21/L11>1.2);
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    
    //
    // Special linear spline test:
    // test for linearity between x[i] and x[i+1]
    //
    N:=2;
    while N<=10 do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        
        //
        // Prepare task
        //
        A := -1;
        B := +1;
        I:=0;
        while I<=N-1 do
        begin
            X[I] := A+(B-A)*I/(N-1);
            Y[I] := 2*RandomReal-1;
            Inc(I);
        end;
        BuildLinearSpline(X, Y, N, C);
        
        //
        // Test
        //
        Err := 0;
        K:=0;
        while K<=N-2 do
        begin
            A := X[K];
            B := X[K+1];
            Pass:=1;
            while Pass<=PassCount do
            begin
                T := A+(B-A)*RandomReal;
                V := Y[K]+(T-A)/(B-A)*(Y[K+1]-Y[K]);
                Err := Max(Err, AbsReal(SplineInterpolation(C, T)-V));
                Inc(Pass);
            end;
            Inc(K);
        end;
        LSErrors := LSErrors or (Err>100*MachineEpsilon);
        Inc(N);
    end;
    
    //
    // Special Akima test: test outlier sensitivity
    // Spline value at (x[i], x[i+1]) should depend from
    // f[i-2], f[i-1], f[i], f[i+1], f[i+2], f[i+3] only.
    //
    N:=5;
    while N<=10 do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        SetLength(Y2, N-1+1);
        
        //
        // Prepare unperturbed Akima spline
        //
        A := -1;
        B := +1;
        I:=0;
        while I<=N-1 do
        begin
            X[I] := A+(B-A)*I/(N-1);
            Y[I] := Cos(1.3*Pi*X[I]+0.4);
            Inc(I);
        end;
        BuildAkimaSpline(X, Y, N, C);
        
        //
        // Process perturbed tasks
        //
        Err := 0;
        K:=0;
        while K<=N-1 do
        begin
            APVMove(@Y2[0], 0, N-1, @Y[0], 0, N-1);
            Y2[K] := 5;
            BuildAkimaSpline(X, Y2, N, C2);
            
            //
            // Test left part independence
            //
            if K-3>=1 then
            begin
                A := -1;
                B := X[K-3];
                Pass:=1;
                while Pass<=PassCount do
                begin
                    T := A+(B-A)*RandomReal;
                    Err := Max(Err, AbsReal(SplineInterpolation(C, T)-SplineInterpolation(C2, T)));
                    Inc(Pass);
                end;
            end;
            
            //
            // Test right part independence