testldaunit.pas

maths lib with source

    J : Integer;
    V : Double;
    WPrev : Double;
    TOL : Double;
    P : Double;
    Q : Double;
    i_ : Integer;
begin
    TOL := 10000;
    Result := True;
    FisherS(XY, NS, NF, NC, ST, SW);
    
    //
    // Test for decreasing of J
    //
    SetLength(TX, NF-1+1);
    SetLength(JP, NF-1+1);
    SetLength(JQ, NF-1+1);
    J:=0;
    while J<=NF-1 do
    begin
        for i_ := 0 to NF-1 do
        begin
            TX[i_] := WN[i_,J];
        end;
        V := CalcJ(NF, ST, SW, TX, P, Q);
        JP[J] := P;
        JQ[J] := Q;
        Inc(J);
    end;
    I:=1;
    while I<=NF-1-NDeg do
    begin
        Result := Result and (JP[I-1]/JQ[I-1]>=(1-TOL*MachineEpsilon)*JP[I]/JQ[I]);
        Inc(I);
    end;
    I:=NF-1-NDeg+1;
    while I<=NF-1 do
    begin
        Result := Result and (JP[I]<=TOL*MachineEpsilon*JP[0]);
        Inc(I);
    end;
    
    //
    // Test for J optimality
    //
    for i_ := 0 to NF-1 do
    begin
        TX[i_] := WN[i_,0];
    end;
    V := CalcJ(NF, ST, SW, TX, P, Q);
    I:=0;
    while I<=NF-1 do
    begin
        WPrev := TX[I];
        TX[I] := WPrev+0.01;
        Result := Result and (V>=(1-TOL*MachineEpsilon)*CalcJ(NF, ST, SW, TX, P, Q));
        TX[I] := WPrev-0.01;
        Result := Result and (V>=(1-TOL*MachineEpsilon)*CalcJ(NF, ST, SW, TX, P, Q));
        TX[I] := WPrev;
        Inc(I);
    end;
    
    //
    // Test for linear independence of W
    //
    SetLength(WORK, NF+1);
    SetLength(A, NF-1+1, NF-1+1);
    MatrixMatrixMultiply(WN, 0, NF-1, 0, NF-1, False, WN, 0, NF-1, 0, NF-1, True, 1.0, A, 0, NF-1, 0, NF-1, 0.0, WORK);
    if SMatrixEVD(A, NF, 1, True, TX, Z) then
    begin
        Result := Result and (TX[0]>TX[NF-1]*1000*MachineEpsilon);
    end;
    
    //
    // Test for other properties
    //
    J:=0;
    while J<=NF-1 do
    begin
        V := 0.0;
        for i_ := 0 to NF-1 do
        begin
            V := V + WN[i_,J]*WN[i_,J];
        end;
        V := Sqrt(V);
        Result := Result and (AbsReal(V-1)<=1000*MachineEpsilon);
        V := 0;
        I:=0;
        while I<=NF-1 do
        begin
            V := V+WN[I,J];
            Inc(I);
        end;
        Result := Result and (V>=0);
        Inc(J);
    end;
end;


(*************************************************************************
Calculates J
*************************************************************************)
function CalcJ(NF : Integer;
     const ST : TReal2DArray;
     const SW : TReal2DArray;
     const W : TReal1DArray;
     var P : Double;
     var Q : Double):Double;
var
    TX : TReal1DArray;
    I : Integer;
    J : Integer;
    V : Double;
begin
    SetLength(TX, NF-1+1);
    I:=0;
    while I<=NF-1 do
    begin
        V := APVDotProduct(@ST[I][0], 0, NF-1, @W[0], 0, NF-1);
        TX[I] := V;
        Inc(I);
    end;
    V := APVDotProduct(@W[0], 0, NF-1, @TX[0], 0, NF-1);
    P := V;
    I:=0;
    while I<=NF-1 do
    begin
        V := APVDotProduct(@SW[I][0], 0, NF-1, @W[0], 0, NF-1);
        TX[I] := V;
        Inc(I);
    end;
    V := APVDotProduct(@W[0], 0, NF-1, @TX[0], 0, NF-1);
    Q := V;
    Result := P/Q;
end;


(*************************************************************************
Calculates ST/SW
*************************************************************************)
procedure FisherS(const XY : TReal2DArray;
     NPoints : Integer;
     NFeatures : Integer;
     NClasses : Integer;
     var ST : TReal2DArray;
     var SW : TReal2DArray);
var
    I : Integer;
    J : Integer;
    K : Integer;
    V : Double;
    C : TInteger1DArray;
    Mu : TReal1DArray;
    MuC : TReal2DArray;
    NC : TInteger1DArray;
    TF : TReal1DArray;
    WORK : TReal1DArray;
begin
    
    //
    // Prepare temporaries
    //
    SetLength(TF, NFeatures-1+1);
    SetLength(WORK, NFeatures+1);
    
    //
    // Convert class labels from reals to integers (just for convenience)
    //
    SetLength(C, NPoints-1+1);
    I:=0;
    while I<=NPoints-1 do
    begin
        C[I] := Round(XY[I,NFeatures]);
        Inc(I);
    end;
    
    //
    // Calculate class sizes and means
    //
    SetLength(Mu, NFeatures-1+1);
    SetLength(MuC, NClasses-1+1, NFeatures-1+1);
    SetLength(NC, NClasses-1+1);
    J:=0;
    while J<=NFeatures-1 do
    begin
        Mu[J] := 0;
        Inc(J);
    end;
    I:=0;
    while I<=NClasses-1 do
    begin
        NC[I] := 0;
        J:=0;
        while J<=NFeatures-1 do
        begin
            MuC[I,J] := 0;
            Inc(J);
        end;
        Inc(I);
    end;
    I:=0;
    while I<=NPoints-1 do
    begin
        APVAdd(@Mu[0], 0, NFeatures-1, @XY[I][0], 0, NFeatures-1);
        APVAdd(@MuC[C[I]][0], 0, NFeatures-1, @XY[I][0], 0, NFeatures-1);
        NC[C[I]] := NC[C[I]]+1;
        Inc(I);
    end;
    I:=0;
    while I<=NClasses-1 do
    begin
        V := 1/NC[I];
        APVMul(@MuC[I][0], 0, NFeatures-1, V);
        Inc(I);
    end;
    V := 1/NPoints;
    APVMul(@Mu[0], 0, NFeatures-1, V);
    
    //
    // Create ST matrix
    //
    SetLength(ST, NFeatures-1+1, NFeatures-1+1);
    I:=0;
    while I<=NFeatures-1 do
    begin
        J:=0;
        while J<=NFeatures-1 do
        begin
            ST[I,J] := 0;
            Inc(J);
        end;
        Inc(I);
    end;
    K:=0;
    while K<=NPoints-1 do
    begin
        APVMove(@TF[0], 0, NFeatures-1, @XY[K][0], 0, NFeatures-1);
        APVSub(@TF[0], 0, NFeatures-1, @Mu[0], 0, NFeatures-1);
        I:=0;
        while I<=NFeatures-1 do
        begin
            V := TF[I];
            APVAdd(@ST[I][0], 0, NFeatures-1, @TF[0], 0, NFeatures-1, V);
            Inc(I);
        end;
        Inc(K);
    end;
    
    //
    // Create SW matrix
    //
    SetLength(SW, NFeatures-1+1, NFeatures-1+1);
    I:=0;
    while I<=NFeatures-1 do
    begin
        J:=0;
        while J<=NFeatures-1 do
        begin
            SW[I,J] := 0;
            Inc(J);
        end;
        Inc(I);
    end;
    K:=0;
    while K<=NPoints-1 do
    begin
        APVMove(@TF[0], 0, NFeatures-1, @XY[K][0], 0, NFeatures-1);
        APVSub(@TF[0], 0, NFeatures-1, @MuC[C[K]][0], 0, NFeatures-1);
        I:=0;
        while I<=NFeatures-1 do
        begin
            V := TF[I];
            APVAdd(@SW[I][0], 0, NFeatures-1, @TF[0], 0, NFeatures-1, V);
            Inc(I);
        end;
        Inc(K);
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testldaunit_test_silent():Boolean;
begin
    Result := TestLDA(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testldaunit_test():Boolean;
begin
    Result := TestLDA(False);
end;


end.