testldaunit.pas

maths lib with source

unit testldaunit;
interface
uses Math, Ap, Sysutils, blas, rotations, tdevd, sblas, reflections, tridiagonal, sevd, cholesky, spdinverse, lda;

function TestLDA(Silent : Boolean):Boolean;
function testldaunit_test_silent():Boolean;
function testldaunit_test():Boolean;

implementation

procedure GenSimpleSet(NFeatures : Integer;
     NClasses : Integer;
     NSamples : Integer;
     Axis : Integer;
     var XY : TReal2DArray);forward;
procedure GenDeg1Set(NFeatures : Integer;
     NClasses : Integer;
     NSamples : Integer;
     Axis : Integer;
     var XY : TReal2DArray);forward;
function GenerateNormal(Mean : Double; Sigma : Double):Double;forward;
function TestWN(const XY : TReal2DArray;
     const WN : TReal2DArray;
     NS : Integer;
     NF : Integer;
     NC : Integer;
     NDeg : Integer):Boolean;forward;
function CalcJ(NF : Integer;
     const ST : TReal2DArray;
     const SW : TReal2DArray;
     const W : TReal1DArray;
     var P : Double;
     var Q : Double):Double;forward;
procedure FisherS(const XY : TReal2DArray;
     NPoints : Integer;
     NFeatures : Integer;
     NClasses : Integer;
     var ST : TReal2DArray;
     var SW : TReal2DArray);forward;


function TestLDA(Silent : Boolean):Boolean;
var
    MaxNF : Integer;
    MaxNS : Integer;
    MaxNC : Integer;
    PassCount : Integer;
    LDANErrors : Boolean;
    LDA1Errors : Boolean;
    WasErrors : Boolean;
    NF : Integer;
    NC : Integer;
    NS : Integer;
    I : Integer;
    Info : Integer;
    Pass : Integer;
    Axis : Integer;
    XY : TReal2DArray;
    WN : TReal2DArray;
    W1 : TReal1DArray;
begin
    
    //
    // Primary settings
    //
    MaxNF := 10;
    MaxNS := 1000;
    MaxNC := 5;
    PassCount := 1;
    WasErrors := False;
    LDANErrors := False;
    LDA1Errors := False;
    
    //
    // General tests
    //
    NF:=1;
    while NF<=MaxNF do
    begin
        NC:=2;
        while NC<=MaxNC do
        begin
            Pass:=1;
            while Pass<=PassCount do
            begin
                
                //
                // Simple test for LDA-N/LDA-1
                //
                Axis := RandomInteger(NF);
                NS := MaxNS div 2+RandomInteger(MaxNS div 2);
                GenSimpleSet(NF, NC, NS, Axis, XY);
                FisherLDAN(XY, NS, NF, NC, Info, WN);
                if Info<>1 then
                begin
                    LDANErrors := True;
                    Inc(Pass);
                    Continue;
                end;
                LDANErrors := LDANErrors or  not TestWN(XY, WN, NS, NF, NC, 0);
                LDANErrors := LDANErrors or (WN[Axis,0]<=0.90);
                FisherLDA(XY, NS, NF, NC, Info, W1);
                I:=0;
                while I<=NF-1 do
                begin
                    LDA1Errors := LDA1Errors or (W1[I]<>WN[I,0]);
                    Inc(I);
                end;
                
                //
                // Degenerate test for LDA-N
                //
                if NF>=3 then
                begin
                    Axis := RandomInteger(NF);
                    NS := MaxNS div 2+RandomInteger(MaxNS div 2);
                    GenDeg1Set(NF, NC, NS, Axis, XY);
                    FisherLDAN(XY, NS, NF, NC, Info, WN);
                    if Info<>2 then
                    begin
                        LDANErrors := True;
                        Inc(Pass);
                        Continue;
                    end;
                    LDANErrors := LDANErrors or  not TestWN(XY, WN, NS, NF, NC, 1);
                    LDANErrors := LDANErrors or (Sqrt(Sqr(WN[NF-1,NF-1])+Sqr(WN[NF-2,NF-1]))<=0.90);
                    LDANErrors := LDANErrors or (AbsReal(WN[NF-1,NF-1]+WN[NF-2,NF-1])>=0.10);
                    FisherLDA(XY, NS, NF, NC, Info, W1);
                    I:=0;
                    while I<=NF-1 do
                    begin
                        LDA1Errors := LDA1Errors or (W1[I]<>WN[I,0]);
                        Inc(I);
                    end;
                end;
                Inc(Pass);
            end;
            Inc(NC);
        end;
        Inc(NF);
    end;
    
    //
    // Final report
    //
    WasErrors := LDANErrors or LDA1Errors;
    if  not Silent then
    begin
        Write(Format('LDA TEST'#13#10'',[]));
        Write(Format('FISHER LDA-N:                            ',[]));
        if  not LDANErrors then
        begin
            Write(Format('OK'#13#10'',[]));
        end
        else
        begin
            Write(Format('FAILED'#13#10'',[]));
        end;
        Write(Format('FISHER LDA-1:                            ',[]));
        if  not LDA1Errors 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;


(*************************************************************************
Generates 'simple' set - a sequence of unit 'balls' at (0,0), (1,0), (2,0)
and so on.
*************************************************************************)
procedure GenSimpleSet(NFeatures : Integer;
     NClasses : Integer;
     NSamples : Integer;
     Axis : Integer;
     var XY : TReal2DArray);
var
    I : Integer;
    J : Integer;
    C : Integer;
    V : Double;
begin
    Assert((Axis>=0) and (Axis<NFeatures), 'GenSimpleSet: wrong Axis!');
    SetLength(XY, NSamples-1+1, NFeatures+1);
    I:=0;
    while I<=NSamples-1 do
    begin
        J:=0;
        while J<=NFeatures-1 do
        begin
            XY[I,J] := GenerateNormal(0.0, 1.0);
            Inc(J);
        end;
        C := I mod NClasses;
        XY[I,Axis] := XY[I,Axis]+C;
        XY[I,NFeatures] := C;
        Inc(I);
    end;
end;


(*************************************************************************
Generates 'degenerate' set #1.
NFeatures>=3.
*************************************************************************)
procedure GenDeg1Set(NFeatures : Integer;
     NClasses : Integer;
     NSamples : Integer;
     Axis : Integer;
     var XY : TReal2DArray);
var
    I : Integer;
    J : Integer;
    C : Integer;
    V : Double;
begin
    Assert((Axis>=0) and (Axis<NFeatures), 'GenDeg1Set: wrong Axis!');
    Assert(NFeatures>=3, 'GenDeg1Set: wrong NFeatures!');
    SetLength(XY, NSamples-1+1, NFeatures+1);
    if Axis>=NFeatures-2 then
    begin
        Axis := NFeatures-3;
    end;
    I:=0;
    while I<=NSamples-1 do
    begin
        J:=0;
        while J<=NFeatures-2 do
        begin
            XY[I,J] := GenerateNormal(0.0, 1.0);
            Inc(J);
        end;
        XY[I,NFeatures-1] := XY[I,NFeatures-2];
        C := I mod NClasses;
        XY[I,Axis] := XY[I,Axis]+C;
        XY[I,NFeatures] := C;
        Inc(I);
    end;
end;


(*************************************************************************
Normal random number
*************************************************************************)
function GenerateNormal(Mean : Double; Sigma : Double):Double;
var
    U : Double;
    V : Double;
    S : Double;
    Sum : Double;
begin
    Result := Mean;
    while True do
    begin
        u := (2*RandomInteger(2)-1)*RandomReal;
        v := (2*RandomInteger(2)-1)*RandomReal;
        sum := u*u+v*v;
        if (sum<1) and (sum>0) then
        begin
            sum := sqrt(-2*ln(sum)/sum);
            Result := Sigma*u*sum+Mean;
            Exit;
        end;
    end;
end;


(*************************************************************************
Tests WN for correctness
*************************************************************************)
function TestWN(const XY : TReal2DArray;
     const WN : TReal2DArray;
     NS : Integer;
     NF : Integer;
     NC : Integer;
     NDeg : Integer):Boolean;
var
    ST : TReal2DArray;
    SW : TReal2DArray;
    A : TReal2DArray;
    Z : TReal2DArray;
    TX : TReal1DArray;
    JP : TReal1DArray;
    JQ : TReal1DArray;
    WORK : TReal1DArray;
    I : Integer;