testlogitunit.pas

maths lib with source

unit testlogitunit;
interface
uses Math, Ap, Sysutils, descriptivestatistics, mlpbase, cholesky, spdsolve, tsort, bdss, logit;

function TestLogit(Silent : Boolean):Boolean;
function testlogitunit_test_silent():Boolean;
function testlogitunit_test():Boolean;

implementation

const
    StrongThreshold = 0.7;
    StrongFreq = 4/5;

procedure DegTests(NF : Integer;
     NC : Integer;
     var DegErrors : Boolean);forward;
procedure STest1(NF : Integer;
     var ConvErrors : Boolean;
     var WeakCorrErrors : Boolean;
     var StrongCorrErrors : Boolean;
     var OtherErrors : Boolean);forward;
procedure STest2(NF : Integer;
     NC : Integer;
     var ConvErrors : Boolean;
     var WeakCorrErrors : Boolean;
     var StrongCorrErrors : Boolean;
     var OtherErrors : Boolean);forward;
procedure STest3(NF : Integer;
     var ConvErrors : Boolean;
     var WeakCorrErrors : Boolean;
     var StrongCorrErrors : Boolean;
     var OtherErrors : Boolean);forward;
procedure AddTest1(NP : Integer;
     NF : Integer;
     NC : Integer;
     const XY : TReal2DArray;
     var W : LogitModel;
     var OtherErrors : Boolean);forward;
procedure AddTest2(NF : Integer;
     NC : Integer;
     var OtherErrors : Boolean);forward;
procedure CVTest(var ConvErrors : Boolean; var OtherErrors : Boolean);forward;
procedure ErrTest(var ConvErrors : Boolean; var OtherErrors : Boolean);forward;
procedure TestProcessing(var ConvErrors : Boolean;
     var ProcErrors : Boolean);forward;
function RNormal():Double;forward;
function RSphere(var XY : TReal2DArray;
     N : Integer;
     I : Integer):Double;forward;
procedure UnsetLM(var LM : LogitModel);forward;


function TestLogit(Silent : Boolean):Boolean;
var
    NF : Integer;
    MaxNF : Integer;
    NC : Integer;
    MaxNC : Integer;
    PassCount : Integer;
    Pass : Integer;
    WasErrors : Boolean;
    DegErrors : Boolean;
    ConvErrors : Boolean;
    OtherErrors : Boolean;
    WeakCorrErrors : Boolean;
    StrongCorrErrors : Boolean;
begin
    
    //
    // Primary settings
    //
    MaxNF := 4;
    MaxNC := 4;
    PassCount := 3;
    WasErrors := False;
    ConvErrors := False;
    OtherErrors := False;
    DegErrors := False;
    WeakCorrErrors := False;
    StrongCorrErrors := False;
    
    //
    // Tests
    //
    NF:=1;
    while NF<=MaxNF do
    begin
        NC:=2;
        while NC<=MaxNC do
        begin
            
            //
            // Degenerate tests
            //
            DegTests(NF, NC, DegErrors);
            
            //
            // General tests
            //
            Pass:=1;
            while Pass<=PassCount do
            begin
                
                //
                // Simple test #1
                //
                if NC=2 then
                begin
                    STest1(NF, ConvErrors, WeakCorrErrors, StrongCorrErrors, OtherErrors);
                    STest3(NF, ConvErrors, WeakCorrErrors, StrongCorrErrors, OtherErrors);
                end;
                
                //
                // Simple test #2
                //
                STest2(NF, NC, ConvErrors, WeakCorrErrors, StrongCorrErrors, OtherErrors);
                Inc(Pass);
            end;
            
            //
            // Additional tests
            //
            AddTest2(NF, NC, OtherErrors);
            Inc(NC);
        end;
        Inc(NF);
    end;
    CVTest(ConvErrors, OtherErrors);
    ErrTest(ConvErrors, OtherErrors);
    TestProcessing(ConvErrors, OtherErrors);
    
    //
    // Final report
    //
    WasErrors := WeakCorrErrors or ConvErrors or StrongCorrErrors or DegErrors or OtherErrors;
    if  not Silent then
    begin
        Write(Format('LOGIT TEST'#13#10'',[]));
        Write(Format('TOTAL RESULTS:                           ',[]));
        if  not WasErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* CONVERGENCE:                           ',[]));
        if  not ConvErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* DEGENERATE CASES:                      ',[]));
        if  not DegErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* WEAK CORRELATION:                      ',[]));
        if  not WeakCorrErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* STRONG CORRELATION:                    ',[]));
        if  not StrongCorrErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('* OTHER PROPERTIES:                      ',[]));
        if  not OtherErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        if WasErrors then
        begin
            Write(Format('TEST SUMMARY: FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST SUMMARY: PASSED'#13#10'',[]));
        end;
        Write(Format(''#13#10''#13#10'',[]));
    end;
    Result :=  not WasErrors;
end;


(*************************************************************************
Degenerate tests
*************************************************************************)
procedure DegTests(NF : Integer; NC : Integer; var DegErrors : Boolean);
var
    I : Integer;
    J : Integer;
    Info : Integer;
    XY : TReal2DArray;
    E : TReal2DArray;
    X : TReal1DArray;
    Y : TReal1DArray;
    NP : Integer;
    W : LogitModel;
    Rep : MNLReport;
begin
    SetLength(X, NF-1+1);
    SetLength(Y, NC-1+1);
    
    //
    // Test #2
    //
    NP:=NF+2;
    while NP<=NF+NF do
    begin
        SetLength(XY, NP-1+1, NF+1);
        I:=0;
        while I<=NP-1 do
        begin
            J:=0;
            while J<=NF-1 do
            begin
                XY[I,J] := 2*RandomReal-1;
                Inc(J);
            end;
            XY[I,NF] := NP mod NC;
            Inc(I);
        end;
        MNLTrainH(XY, NP, NF, NC, Info, W, Rep);
        if Info<>1 then
        begin
            DegErrors := True;
            Exit;
        end;
        I:=0;
        while I<=NF-1 do
        begin
            X[I] := 2*RandomReal-1;
            Inc(I);
        end;
        MNLProcess(W, X, Y);
        I:=0;
        while I<=NC-1 do
        begin
            if I=NP mod NC then
            begin
                DegErrors := DegErrors or (Y[I]<>1);
            end
            else
            begin
                DegErrors := DegErrors or (Y[I]<>0);
            end;
            Inc(I);
        end;
        Inc(NP);
    end;
end;


(*************************************************************************
Simple test 1.
Two well-separated balls.
*************************************************************************)
procedure STest1(NF : Integer;
     var ConvErrors : Boolean;
     var WeakCorrErrors : Boolean;
     var StrongCorrErrors : Boolean;
     var OtherErrors : Boolean);
var
    I : Integer;
    J : Integer;
    Info : Integer;
    XY : TReal2DArray;
    X : TReal1DArray;
    Y : TReal1DArray;
    C : TReal1DArray;
    W : LogitModel;
    NC : Integer;
    NP : Integer;
    S : Double;
    V : Double;
    StrongCnt : Integer;
    WeakCnt : Integer;
    Rep : MNLReport;
begin
    NP := 2*(10+2*NF+RandomInteger(NF));
    NC := 2;
    SetLength(X, NF-1+1);
    SetLength(C, NF-1+1);
    SetLength(Y, NC-1+1);
    SetLength(XY, NP-1+1, NF+1);
    
    //
    // Fill
    //
    J:=0;
    while J<=NF-1 do
    begin
        C[J] := 2*RandomReal-1;
        Inc(J);
    end;
    S := APVDotProduct(@C[0], 0, NF-1, @C[0], 0, NF-1);
    S := Sqrt(S);
    I:=0;
    while I<=NP-1 do
    begin
        RSphere(XY, NF, I);
        V := 0.2*S;
        APVMul(@XY[I][0], 0, NF-1, V);
        if I mod 2=0 then
        begin
            
            //
            // label 0
            //
            APVAdd(@XY[I][0], 0, NF-1, @C[0], 0, NF-1);
            XY[I,NF] := 0;
        end
        else
        begin
            
            //
            // label 1
            //
            APVSub(@XY[I][0], 0, NF-1, @C[0], 0, NF-1);
            XY[I,NF] := 1;
        end;
        Inc(I);
    end;
    
    //
    // Train
    //
    MNLTrainH(XY, NP, NF, 2, Info, W, Rep);
    if Info<>1 then
    begin
        ConvErrors := True;
        Exit;
    end;
    
    //
    // Test
    //
    WeakCnt := 0;
    StrongCnt := 0;
    I:=0;
    while I<=NP-1 do
    begin
        APVMove(@X[0], 0, NF-1, @XY[I][0], 0, NF-1);
        MNLProcess(W, X, Y);
        if Y[I mod 2]>0.5 then
        begin
            WeakCnt := WeakCnt+1;
        end;
        if Y[I mod 2]>StrongThreshold then
        begin
            StrongCnt := StrongCnt+1;
        end;
        J:=0;
        while J<=NC-1 do
        begin
            OtherErrors := OtherErrors or (Y[J]<0) or (Y[J]>1);
            Inc(J);
        end;
        OtherErrors := OtherErrors or (AbsReal(Y[0]+Y[1]-1)>100*MachineEpsilon);
        Inc(I);
    end;
    StrongCorrErrors := StrongCorrErrors or (StrongCnt/NP<StrongFreq);
    WeakCorrErrors := WeakCorrErrors or (WeakCnt/NP<2/3);
    
    //
    // Additional tests
    //
    AddTest1(NP, NF, NC, XY, W, OtherErrors);
    OtherErrors := OtherErrors or (MNLClsError(W, XY, NP)<>NP-WeakCnt);
end;


(*************************************************************************
Simple test 2:
NC unit balls separated by 3 units wide gaps
*************************************************************************)
procedure STest2(NF : Integer;
     NC : Integer;
     var ConvErrors : Boolean;
     var WeakCorrErrors : Boolean;
     var StrongCorrErrors : Boolean;
     var OtherErrors : Boolean);
var
    I : Integer;
    J : Integer;
    Info : Integer;
    XY : TReal2DArray;
    X : TReal1DArray;
    Y : TReal1DArray;
    C : TReal1DArray;
    W : LogitModel;
    NP : Integer;
    S : Double;
    V : Double;
    WeakCnt : Integer;
    StrongCnt : Integer;
    Rep : MNLReport;
begin
    NP := NC*(10+2*NF+RandomInteger(NF));
    SetLength(X, NF-1+1);
    SetLength(C, NF-1+1);
    SetLength(Y, NC-1+1);
    SetLength(XY, NP-1+1, NF+1);
    
    //
    // Fill
    //
    J:=0;
    while J<=NF-1 do
    begin
        C[J] := RNormal;
        Inc(J);
    end;
    S := APVDotProduct(@C[0], 0, NF-1, @C[0], 0, NF-1);
    S := 1/Sqrt(S);
    APVMul(@C[0], 0, NF-1, S);
    I:=0;
    while I<=NP-1 do
    begin
        RSphere(XY, NF, I);
        V := 2*(I mod NC);
        APVAdd(@XY[I][0], 0, NF-1, @C[0], 0, NF-1, V);
        XY[I,NF] := I mod NC;
        Inc(I);
    end;
    
    //
    // Train
    //
    MNLTrainH(XY, NP, NF, NC, Info, W, Rep);
    if Info<>1 then
    begin
        ConvErrors := True;
        Exit;
    end;
    
    //
    // Test
    //
    WeakCnt := 0;
    StrongCnt := 0;
    I:=0;
    while I<=NP-1 do
    begin
        APVMove(@X[0], 0, NF-1, @XY[I][0], 0, NF-1);
        MNLProcess(W, X, Y);
        if Y[I mod NC]>StrongThreshold then
        begin
            StrongCnt := StrongCnt+1;
        end;
        if Y[I mod NC]>0.5 then
        begin
            WeakCnt := WeakCnt+1;
        end;
        S := 0;
        J:=0;
        while J<=NC-1 do
        begin
            S := S+Y[J];
            OtherErrors := OtherErrors or (Y[J]<0) or (Y[J]>1);
            Inc(J);
        end;
        OtherErrors := OtherErrors or (AbsReal(S-1)>100*MachineEpsilon);
        Inc(I);
    end;
    StrongCorrErrors := StrongCorrErrors or (StrongCnt/NP<StrongFreq);
    WeakCorrErrors := WeakCorrErrors or (WeakCnt/NP<2/3);
    
    //
    // Additional tests
    //
    AddTest1(NP, NF, NC, XY, W, OtherErrors);
    OtherErrors := OtherErrors or (MNLClsError(W, XY, NP)<>NP-WeakCnt);
end;


(*************************************************************************
Simple test 3.
Two badly-separated balls.
*************************************************************************)
procedure STest3(NF : Integer;
     var ConvErrors : Boolean;
     var WeakCorrErrors : Boolean;
     var StrongCorrErrors : Boolean;
     var OtherErrors : Boolean);
var
    I : Integer;
    J : Integer;
    Info : Integer;
    XY : TReal2DArray;
    X : TReal1DArray;
    Y : TReal1DArray;
    C : TReal1DArray;
    W : LogitModel;
    NC : Integer;
    NP : Integer;
    S : Double;
    V : Double;
    WeakCnt : Integer;
    Rep : MNLReport;
begin
    NP := 100;
    NC := 2;
    SetLength(X, NF-1+1);
    SetLength(C, NF-1+1);
    SetLength(Y, NC-1+1);
    SetLength(XY, NP-1+1, NF+1);
    
    //
    // Fill
    //
    J:=0;
    while J<=NF-1 do
    begin
        C[J] := 2*RandomReal-1;
        Inc(J);
    end;
    S := APVDotProduct(@C[0], 0, NF-1, @C[0], 0, NF-1);
    S := Sqrt(S);
    I:=0;
    while I<=NP-1 do
    begin
        RSphere(XY, NF, I);
        V := 1.5*S;
        APVMul(@XY[I][0], 0, NF-1, V);
        if I mod 2=0 then
        begin
            
            //
            // label 0
            //
            APVAdd(@XY[I][0], 0, NF-1, @C[0], 0, NF-1);