llstestunit.pas

maths lib with source

unit llstestunit;
interface
uses Math, Ap, Sysutils, spline3, reflections, lq, bidiagonal, rotations, bdsvd, qr, blas, svd, leastsquares;

function TestLLS(Silent : Boolean):Boolean;
function IsGLSSolution(N : Integer;
     M : Integer;
     const Y : TReal1DArray;
     const W : TReal1DArray;
     const FMatrix : TReal2DArray;
     const C : TReal1DArray):Boolean;
function llstestunit_test_silent():Boolean;
function llstestunit_test():Boolean;

implementation

function TestLLS(Silent : Boolean):Boolean;
var
    WasErrors : Boolean;
    GLSErrors : Boolean;
    PLSErrors : Boolean;
    LinLSErrors : Boolean;
    Threshold : Double;
    N : Integer;
    M : Integer;
    MaxN : Integer;
    MaxM : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    Pass : Integer;
    PassCount : Integer;
    CTask : Integer;
    XScale : Double;
    X : TReal1DArray;
    Y : TReal1DArray;
    W : TReal1DArray;
    C : TReal1DArray;
    AC : TReal1DArray;
    C1 : Double;
    C2 : Double;
    V : Double;
    S1 : Double;
    S2 : Double;
    S3 : Double;
    Delta : Double;
    T : Double;
    VG : Double;
    VK : Double;
    VK1 : Double;
    VK2 : Double;
    VP : Double;
    A : TReal2DArray;
begin
    WasErrors := False;
    GLSErrors := False;
    PLSErrors := False;
    LinLSErrors := False;
    Threshold := 10000*MachineEpsilon;
    MaxN := 10;
    MaxM := 5;
    PassCount := 10;
    Delta := 0.001;
    
    //
    // Testing general least squares
    //
    N:=1;
    while N<=MaxN do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        SetLength(W, N-1+1);
        SetLength(AC, N-1+1);
        M:=1;
        while M<=MaxM do
        begin
            SetLength(A, N-1+1, M-1+1);
            Pass:=1;
            while Pass<=PassCount do
            begin
                
                //
                // Prepare task.
                // Use Chebyshev basis. Its condition number is very good.
                //
                XScale := 0.9+0.1*RandomReal;
                I:=0;
                while I<=N-1 do
                begin
                    if N=1 then
                    begin
                        X[I] := 2*RandomReal-1;
                    end
                    else
                    begin
                        X[I] := XScale*(2*I/(N-1)-1);
                    end;
                    Y[I] := 3*X[I]+Exp(X[I]);
                    W[I] := 1+RandomReal;
                    A[I,0] := 1;
                    if M>1 then
                    begin
                        A[I,1] := X[I];
                    end;
                    J:=2;
                    while J<=M-1 do
                    begin
                        A[I,J] := 2*X[I]*A[I,J-1]-A[I,J-2];
                        Inc(J);
                    end;
                    Inc(I);
                end;
                
                //
                // Solve General Least Squares task
                //
                BuildGeneralLeastSquares(Y, W, A, N, M, C);
                GLSErrors := GLSErrors or  not IsGLSSolution(N, M, Y, W, A, C);
                Inc(Pass);
            end;
            Inc(M);
        end;
        Inc(N);
    end;
    
    //
    // Test polynomial least squares
    //
    N:=1;
    while N<=MaxN do
    begin
        SetLength(X, N-1+1);
        SetLength(Y, N-1+1);
        SetLength(W, N-1+1);
        SetLength(AC, N-1+1);
        M:=1;
        while M<=MaxM do
        begin
            SetLength(A, N-1+1, M-1+1);
            Pass:=1;
            while Pass<=PassCount do
            begin
                
                //
                // Prepare task.
                // Use power basis.
                //
                XScale := 0.9+0.1*RandomReal;
                I:=0;
                while I<=N-1 do
                begin
                    if N=1 then
                    begin
                        X[I] := 2*RandomReal-1;
                    end
                    else
                    begin
                        X[I] := XScale*(2*I/(N-1)-1);
                    end;
                    Y[I] := 3*X[I]+Exp(X[I]);
                    W[I] := 1;
                    A[I,0] := 1;
                    J:=1;
                    while J<=M-1 do
                    begin
                        A[I,J] := X[I]*A[I,J-1];
                        Inc(J);
                    end;
                    Inc(I);
                end;
                
                //
                // Solve polynomial least squares task
                //
                BuildPolynomialLeastSquares(X, Y, N, M-1, C);
                PLSErrors := PLSErrors or  not IsGLSSolution(N, M, Y, W, A, C);
                Inc(Pass);
            end;
            Inc(M);
        end;
        Inc(N);
    end;
    
    //
    // Linear approximation.
    // These tests are easy to do, but I think it will be enough
    //
    N:=1;
    while N<=MaxN do
    begin
        SetLength(X, 2*N-1+1);
        SetLength(Y, 2*N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // Generate y = C1 + C2*x
            // Generate N pairs of points: (x, y(x)+s1) and (x, y(x)-s1)
            // C1 and C2 must be calculated exactly
            //
            C1 := 2*RandomReal-1;
            C2 := 2*RandomReal-1;
            S1 := 1;
            I:=0;
            while I<=N-1 do
            begin
                Inc(I);
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    
    //
    //
    // report
    //
    WasErrors := GLSErrors or PLSErrors;
    if  not Silent then
    begin
        Write(Format('TESTING LINEAR LEAST SQUARES'#13#10'',[]));
        Write(Format('GENERAL LLS                              ',[]));
        if GLSErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('POLYNOMIAL LLS                           ',[]));
        if PLSErrors 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;


function IsGLSSolution(N : Integer;
     M : Integer;
     const Y : TReal1DArray;
     const W : TReal1DArray;
     const FMatrix : TReal2DArray;
     const C : TReal1DArray):Boolean;
var
    I : Integer;
    J : Integer;
    AC : TReal1DArray;
    V : Double;
    S1 : Double;
    S2 : Double;
    S3 : Double;
    Delta : Double;
begin
    Delta := 0.001;
    SetLength(AC, N-1+1);
    
    //
    // Test result
    //
    Result := True;
    I:=0;
    while I<=N-1 do
    begin
        V := APVDotProduct(@FMatrix[I][0], 0, M-1, @C[0], 0, M-1);
        AC[I] := V;
        Inc(I);
    end;
    S1 := 0;
    I:=0;
    while I<=N-1 do
    begin
        S1 := S1+Sqr(W[I]*(AC[I]-Y[I]));
        Inc(I);
    end;
    J:=0;
    while J<=M-1 do
    begin
        S2 := 0;
        S3 := 0;
        I:=0;
        while I<=N-1 do
        begin
            S2 := S2+Sqr(W[I]*(AC[I]+FMatrix[I,J]*Delta-Y[I]));
            S3 := S3+Sqr(W[I]*(AC[I]-FMatrix[I,J]*Delta-Y[I]));
            Inc(I);
        end;
        Result := Result and (S2>=S1) and (S3>=S1);
        Inc(J);
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function llstestunit_test_silent():Boolean;
begin
    Result := TestLLS(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function llstestunit_test():Boolean;
begin
    Result := TestLLS(False);
end;


end.