testsplineinterpolationunit.pas

maths lib with source

            //
            if K+3<=N-2 then
            begin
                A := X[K+3];
                B := +1;
                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;
            Inc(K);
        end;
        ASErrors := ASErrors or (Err>100*MachineEpsilon);
        Inc(N);
    end;
    
    //
    // Differentiation, unpack test
    //
    N:=2;
    while N<=10 do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        
        //
        // Prepare cubic spline
        //
        A := -1-RandomReal;
        B := +1+RandomReal;
        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;
        BuildCubicSpline(X, Y, N, 2, 0.0, 2, 0.0, C);
        
        //
        // Test diff
        //
        Err := 0;
        Pass:=1;
        while Pass<=PassCount do
        begin
            T := A+(B-A)*RandomReal;
            SplineDifferentiation(C, T, S, DS, D2S);
            VL := SplineInterpolation(C, T-H);
            VM := SplineInterpolation(C, T);
            VR := SplineInterpolation(C, T+H);
            Err := Max(Err, AbsReal(S-VM));
            Err := Max(Err, AbsReal(DS-(VR-VL)/(2*H)));
            Err := Max(Err, AbsReal(D2S-(VR-2*VM+VL)/Sqr(H)));
            Inc(Pass);
        end;
        DSErrors := DSErrors or (Err>0.001);
        
        //
        // Test copy
        //
        SplineCopy(C, C2);
        Err := 0;
        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;
        CPErrors := CPErrors or (Err>100*MachineEpsilon);
        
        //
        // Test unpack
        //
        UPErrors := UPErrors or  not TestUnpack(C, X);
        
        //
        // Test lin.trans.
        //
        Err := 0;
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // LinTransX, general A
            //
            SA := 4*RandomReal-2;
            SB := 2*RandomReal-1;
            T := A+(B-A)*RandomReal;
            SplineCopy(C, C2);
            SplineLinTransX(C2, SA, SB);
            Err := Max(Err, AbsReal(SplineInterpolation(C, T)-SplineInterpolation(C2, (T-SB)/SA)));
            
            //
            // LinTransX, special case: A=0
            //
            SB := 2*RandomReal-1;
            T := A+(B-A)*RandomReal;
            SplineCopy(C, C2);
            SplineLinTransX(C2, 0, SB);
            Err := Max(Err, AbsReal(SplineInterpolation(C, SB)-SplineInterpolation(C2, T)));
            
            //
            // LinTransY
            //
            SA := 2*RandomReal-1;
            SB := 2*RandomReal-1;
            T := A+(B-A)*RandomReal;
            SplineCopy(C, C2);
            SplineLinTransY(C2, SA, SB);
            Err := Max(Err, AbsReal(SA*SplineInterpolation(C, T)+SB-SplineInterpolation(C2, T)));
            Inc(Pass);
        end;
        LTErrors := LTErrors or (Err>100*MachineEpsilon);
        Inc(N);
    end;
    
    //
    // Testing integration
    //
    Err := 0;
    N:=20;
    while N<=35 do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Prepare cubic spline
            //
            A := -1-0.2*RandomReal;
            B := +1+0.2*RandomReal;
            I:=0;
            while I<=N-1 do
            begin
                X[I] := A+(B-A)*I/(N-1);
                Y[I] := Sin(Pi*X[I]+0.4)+Exp(X[I]);
                Inc(I);
            end;
            BL := Pi*Cos(Pi*A+0.4)+Exp(A);
            BR := Pi*Cos(Pi*B+0.4)+Exp(B);
            BuildCubicSpline(X, Y, N, 1, BL, 1, BR, C);
            
            //
            // Test
            //
            T := A+(B-A)*RandomReal;
            V := -Cos(Pi*A+0.4)/Pi+Exp(A);
            V := -Cos(Pi*T+0.4)/Pi+Exp(T)-V;
            V := V-SplineIntegration(C, T);
            Err := Max(Err, AbsReal(V));
            Inc(Pass);
        end;
        Inc(N);
    end;
    IErrors := IErrors or (Err>0.001);
    
    //
    // report
    //
    WasErrors := LSErrors or CSErrors or HSErrors or ASErrors or DSErrors or CPErrors or UPErrors or LTErrors or IErrors;
    if  not Silent then
    begin
        Write(Format('TESTING SPLINE INTERPOLATION'#13#10'',[]));
        
        //
        // Normal tests
        //
        Write(Format('LINEAR SPLINE TEST:                      ',[]));
        if LSErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('CUBIC SPLINE TEST:                       ',[]));
        if CSErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('HERMITE SPLINE TEST:                     ',[]));
        if HSErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('AKIMA SPLINE TEST:                       ',[]));
        if ASErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('DIFFERENTIATION TEST:                    ',[]));
        if DSErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('COPY TEST:                               ',[]));
        if CPErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('UNPACK TEST:                             ',[]));
        if UPErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('LIN.TRANS. TEST:                         ',[]));
        if LTErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('INTEGRATION TEST:                        ',[]));
        if IErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        if WasErrors then
        begin
            Write(Format('TEST FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST PASSED'#13#10'',[]));
        end;
        Write(Format(''#13#10''#13#10'',[]));
    end;
    
    //
    // end
    //
    Result :=  not WasErrors;
end;


(*************************************************************************
Lipschitz constants for spline inself, first and second derivatives.
*************************************************************************)
procedure LConst(A : Double;
     B : Double;
     const C : TReal1DArray;
     LStep : Double;
     var L0 : Double;
     var L1 : Double;
     var L2 : Double);
var
    T : Double;
    VL : Double;
    VM : Double;
    VR : Double;
    PrevF : Double;
    PrevD : Double;
    PrevD2 : Double;
    F : Double;
    D : Double;
    D2 : Double;
begin
    L0 := 0;
    L1 := 0;
    L2 := 0;
    T := A-0.1;
    VL := SplineInterpolation(C, T-2*LStep);
    VM := SplineInterpolation(C, T-LStep);
    VR := SplineInterpolation(C, T);
    F := VM;
    D := (VR-VL)/(2*LStep);
    D2 := (VR-2*VM+VL)/Sqr(LStep);
    while T<=B+0.1 do
    begin
        PrevF := F;
        PrevD := D;
        PrevD2 := D2;
        VL := VM;
        VM := VR;
        VR := SplineInterpolation(C, T+LStep);
        F := VM;
        D := (VR-VL)/(2*LStep);
        D2 := (VR-2*VM+VL)/Sqr(LStep);
        L0 := Max(L0, AbsReal((F-PrevF)/LStep));
        L1 := Max(L1, AbsReal((D-PrevD)/LStep));
        L2 := Max(L2, AbsReal((D2-PrevD2)/LStep));
        T := T+LStep;
    end;
end;


(*************************************************************************
Lipschitz constants for spline inself, first and second derivatives.
*************************************************************************)
function TestUnpack(const C : TReal1DArray; const X : TReal1DArray):Boolean;
var
    I : Integer;
    N : Integer;
    Err : Double;
    T : Double;
    V1 : Double;
    V2 : Double;
    Pass : Integer;
    PassCount : Integer;
    Tbl : TReal2DArray;
begin
    PassCount := 20;
    Err := 0;
    SplineUnpack(C, N, Tbl);
    I:=0;
    while I<=N-2 do
    begin
        Pass:=1;
        while Pass<=PassCount do
        begin
            T := RandomReal*(Tbl[I,1]-Tbl[I,0]);
            V1 := Tbl[I,2]+T*Tbl[I,3]+Sqr(T)*Tbl[I,4]+T*Sqr(T)*Tbl[I,5];
            V2 := SplineInterpolation(C, Tbl[I,0]+T);
            Err := Max(Err, AbsReal(V1-V2));
            Inc(Pass);
        end;
        Inc(I);
    end;
    I:=0;
    while I<=N-2 do
    begin
        Err := Max(Err, AbsReal(X[I]-Tbl[I,0]));
        Inc(I);
    end;
    I:=0;
    while I<=N-2 do
    begin
        Err := Max(Err, AbsReal(X[I+1]-Tbl[I,1]));
        Inc(I);
    end;
    Result := Err<100*MachineEpsilon;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testsplineinterpolationunit_test_silent():Boolean;
begin
    Result := TestSplineInterpolation(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testsplineinterpolationunit_test():Boolean;
begin
    Result := TestSplineInterpolation(False);
end;


end.