IncrementalPCA.pas

// ###################################################################
// #### This file is part of the mathematics library project, and is
// #### offered under the licence agreement described on
// #### http://www.mrsoft.org/
// ####
// #### Copyright:(c) 2011, Michael R. . All rights reserved.
// ####
// #### Unless required by applicable law or agreed to in writing, software
// #### distributed under the License is distributed on an "AS IS" BASIS,
// #### WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// #### See the License for the specific language governing permissions and
// #### limitations under the License.
// ###################################################################

unit IncrementalPCA;

// ###############################################
// #### Incremental version of the PCA algorithm
// ###############################################

interface

uses SysUtils, PCA, Matrix, BaseMathPersistence, Types;

type
  TIncrPCAData = record
    Props : TPCAData;
    Weights : TDoubleDynArray;
    NumEigenVecs : integer;
    A : TDoubleMatrix;
  end;

// ###############################################
// #### incremental version of the base pca algorithm
// based on: Danijel Skocaj
// Robust Subspace Approaches to Visual Learning and Recognition
type
  TIncrementalPCA = class(TMatrixPCA)
  private
    fNumEigenvectorsToKeep: integer;
    fA : IMatrix;
    fWeights : TDoubleDynArray;
    fNumWeights : integer;

    procedure AddWeight(const weight : double);
    function GetA: TDoubleMatrix;
  protected
    procedure DefineProps; override;
    function PropTypeOfName(const Name: string): TPropType; override;
    class function ClassIdentifier : String; override;
    function OnLoadObject(const Name : String; obj : TBaseMathPersistence) : boolean; override;
    procedure OnLoadIntProperty(const Name : String; Value : integer); override;
    procedure OnLoadDoubleArr(const Name : String; const Value : TDoubleDynArray); override;
  public
    property NumEigenvectorsToKeep : integer read fNumEigenvectorsToKeep write fNumEigenvectorsToKeep;
    property A : TDoubleMatrix read GetA;

    function UpdateEigenspace(mtx : TDoubleMatrix) : boolean; overload;
    // this function is an extension to the base incremental pca algorithm -> it uses the same blending algorithm
    // as the weighted incremental version (this is works since it is only a special case of the weighted version!)
    function UpdateEigenspace(mtx : TDoubleMatrix; const spacialWeights : Array of double) : boolean; overload;
    // note the spacial weights MUST have values between 0 and 1. 1 means no weighting and 0 means that pixel is not taken into account.
    // implementation of Daniejl S. weighted incremental PCA function (p 69ff):
    // this function supports spacial weights as well as temporal weights. If the function without spacial weights is called
    // all weights are set to 1 (no blending with the original image)
    function UpdateEigenspaceWeighted(mtx : TDoubleMatrix; const weight : double) : boolean; overload;
    function UpdateEigenspaceWeighted(mtx : TDoubleMatrix; const weight : double; const spacialWeights : Array of double) : boolean; overload;

    function PCA(Examples : TDoubleMatrix; CutEps : double; IsRelativeValue : boolean) : boolean; override;
    function TemporalWeightPCA(Examples : TDoubleMatrix; CutEps : double; IsRelativeValue : boolean; const Weights : Array of Double) : boolean; override;

    constructor Create(const Data : TIncrPCAData); overload;
    constructor Create(KeepFlags : TPCASaveDataSet); overload;
    destructor Destroy; override;
  end;

type
  TFastRobustPCADataEx = record
    RData : TFastRobustPCAData;
    A : TDoubleMatrix;
    Weights : TDoubleDynArray;
    NumEigenVecs : integer;
    SubA : IMatrixDynArr;
    SubWeights : Array of TDoubleDynArray;
  end;

// ###############################################
// #### incremental version of the fast robust pca algorithm
// a combination of the base fast robust pca method and the incremental pca algorithm
type
  TFastRobustIncrementalPCA = class(TFastRobustPCA)
  private
    fNumEigenvectorsToKeep : integer;
    fA : IMatrix;
    fSubA : IMatrixDynArr;
    fWeights : Array of TDoubleDynArray;
    fNumWeights : TIntegerDynArray;
    fSubEigVals : IMatrixDynArr;

    procedure AddWeight(idx : integer; const weight : double);
    procedure BuildSubspaces(Examples : TDoubleMatrix; origKeepFlags : TPCASaveDataSet);
  protected
    type
      TListLoadType = (llNone, llSubA, llWeights, llSubEigVals);
  private
    fIdx : integer;
    fListType : TListLoadType;
    function getA: TDoubleMatrix;
  protected
    procedure DefineProps; override;
    function PropTypeOfName(const Name : string) : TPropType; override;

    class function ClassIdentifier : String; override;
    procedure OnLoadBeginList(const Name : String; count : integer); override;
    procedure OnLoadEndList; override;
    function OnLoadObject(Obj : TBaseMathPersistence) : boolean; overload; override;
    function OnLoadObject(const Name : String; obj : TBaseMathPersistence) : boolean; overload; override;
    procedure OnLoadListDoubleArr(const Value : TDoubleDynArray); override;
    procedure OnLoadIntArr(const Name : String; const Value : TIntegerDynArray); override;
    procedure OnLoadIntProperty(const Name : String; Value : integer); overload; override;
  protected
    procedure Clear; override;
  public
    property NumEigenvectorsToKeep : integer read fNumEigenvectorsToKeep write fNumEigenvectorsToKeep;
    property A : TDoubleMatrix read getA;

    function UpdateEigenspace(mtx : TDoubleMatrix) : boolean; overload;
    // this function is an extension to the base incremental pca algorithm -> it uses the same blending algorithm
    // as the weighted incremental version (this is works since it is only a special case of the weighted version!)
    function UpdateEigenspace(mtx : TDoubleMatrix; const spacialWeights : Array of double) : boolean; overload;
    // note the spacial weights MUST have values between 0 and 1. 1 means no weighting and 0 means that pixel is not taken into account.
    // implementation of Daniejl S. weighted incremental PCA function (p 69ff):
    // this function supports spacial weights as well as temporal weights. If the function without spacial weights is called
    // all weights are set to 1 (no blending with the original image)
    function UpdateEigenspaceWeighted(mtx : TDoubleMatrix; const weight : double) : boolean; overload;
    function UpdateEigenspaceWeighted(mtx : TDoubleMatrix; const weight : double; const spacialWeights : Array of double) : boolean; overload;

    function PCA(Examples : TDoubleMatrix; CutEps : double; IsRelativeValue : boolean) : boolean; override;
    function TemporalWeightPCA(Examples : TDoubleMatrix; CutEps : double; IsRelativeValue : boolean; const Weights : Array of Double) : boolean; override;

    constructor Create(const Props : TFastRobustPCADataEx); overload;
    constructor Create(KeepFlags : TPCASaveDataSet); overload;
    destructor Destroy; override;
  end;


implementation

uses Math, MatrixConst, MathUtilFunc;

// ####################################################################
// #### Local functions - used to update the eigenspaces
// ####################################################################

// implements the update Eigenspace functions and holds references to all necessary data structures
type
  TSubspaceUpdater = class(TMatrixClass)
  private
    fNumEigenvectorsToKeep : integer;
    fEigVecs : IMatrix;
    fEigVals : IMatrix;
    fMeanElem : IMatrix;
    fA : IMatrix;
    fWeights : TDoubleDynArray;
    fNumWeights : integer;

    function InternalUpdateEigenspace(mtx : TDoubleMatrix; weighted : boolean) : boolean;
    function PrepareBlendExample(mtx : TDoubleMatrix; const spacialWeights : Array of double) : IMatrix;
  public
    function UpdateEigenspace(mtx : TDoubleMatrix; var EigVecs, EigVals, MeanElem, A : IMatrix) : boolean; overload;
    // this function is an extension to the base incremental pca algorithm -> it uses the same blending algorithm
    // as the weighted incremental version (this is works since it is only a special case of the weighted version!)
    function UpdateEigenspace(mtx : TDoubleMatrix; const spacialWeights : Array of double; var EigVecs, EigVals, MeanElem, A : IMatrix) : boolean; overload;
    // note the spacial weights MUST have values between 0 and 1. 1 means no weighting and 0 means that pixel is not taken into account.
    // implementation of Daniejl S. weighted incremental PCA function (p 69ff):
    // this function supports spacial weights as well as temporal weights. If the function without spacial weights is called
    // all weights are set to 1 (no blending with the original image)
    function UpdateEigenspaceWeighted(mtx : TDoubleMatrix; var weights : TDoubleDynArray; var numWeights : integer; var EigVecs, EigVals, MeanElem, A : IMatrix) : boolean; overload;
    function UpdateEigenspaceWeighted(mtx : TDoubleMatrix; var weights : TDoubleDynArray; var numWeights : integer; const spacialWeights : Array of double; var EigVecs, EigVals, MeanElem, A : IMatrix) : boolean; overload;

    constructor Create(NumVectorsToKeep : integer);
  end;

{ TSubspaceUpdater }

constructor TSubspaceUpdater.Create(NumVectorsToKeep: integer);
begin
     inherited Create;

     fNumEigenvectorsToKeep := NumVectorsToKeep;
end;

function TSubspaceUpdater.InternalUpdateEigenspace(mtx: TDoubleMatrix; weighted: boolean): boolean;
var normExample : IMatrix;
    a : IMatrix;
    r : IMatrix;
    a1 : IMatrix;
    u1 : IMatrix;
    normR : double;
    y : Integer;
    x : Integer;
    pca : TMatrixPCA;
    temp : IMatrix;
    eigVecsT : IMatrix;
begin
     Result := False;

     // project current example to the current feature space
     normExample := mtx.Sub(fMeanElem);
     eigVecsT := fEigVecs.Transpose;
     a := EigVecsT.Mult(normExample);
     eigVecsT := nil;

     // reconstruct example
     normExample := fEigVecs.Mult(a);
     normExample.AddInplace(fMeanElem);

     // compute residual vector and normalize
     r := mtx.Sub(normExample);
     normR := 0;
     for y := 0 to r.Height - 1 do
         normR := normR + sqr(r[0, y]);
     normR := sqrt(normR);
     // check error -> todo: better Boundary check
     normExample := nil;
     if normR < 100*MinDouble then
        exit;

     // normalize vector
     r.ScaleInPlace(1/normR);

     // append new basis vector to u
     u1 := MatrixClass.Create(fEigVecs.Width + 1, fEigVecs.Height);
     u1.AssignSubMatrix(fEigVecs);
     u1.SetColumn(u1.Width - 1, r);

     r := nil;

     // append new value in feature space
     a1 := MatrixClass.Create(fA.Width + 1, fA.Height + 1);
     a1.SetSubMatrix(0, 0, a1.Width, a1.Height - 1);
     for x := 0 to fA.Width - 1 do
         a1.SetColumn(x, fA, x);
     a1.SetColumn(a1.Width - 1, a);
     a1.UseFullMatrix;
     a1[a1.width - 1, a1.Height - 1] := normR;

     a := nil;

     // perform full PCA on the new matrix a1
     pca := TMatrixPCA.Create([pcaEigVals]);
     try
        if not weighted
        then
            Result := pca.PCA(a1.GetObjRef, 1, True)
        else
            Result := pca.TemporalWeightPCA(a1.GetObjRef, 1, True, Copy(fWeights, 0, fNumWeights));

        if not Result then
           exit;

        // project the coefficient vectors to the new basis
        a1.SubVecInPlace(pca.Mean, False);
        fA := pca.EigVecsT.Mult(a1);
        a1 := nil;

        // update the mean
        temp := u1.Mult(pca.Mean);
        fMeanElem.AddInplace(temp);
        temp := nil;

        // rotate the subspace
        u1.MultInPlace(pca.EigVecs);
        fEigVecs := u1;

        // discard the last eigenvector -> keep the basis constant
        if (fNumEigenvectorsToKeep > 0) and (fEigVecs.Width > fNumEigenvectorsToKeep) then
        begin
             fEigVecs.Resize(fEigVecs.Width - 1, fEigVecs.Height);

             // remove the last row -> it contains the values of the least significant
             // elements
             fA.SetSubMatrix(0, 0, fA.Width, fA.Height - 1);
        end;

        if Assigned(fEigVals) then
           fEigVals.Assign(pca.EigVals);
     finally
            pca.Free;
     end;
end;

function TSubspaceUpdater.PrepareBlendExample(mtx: TDoubleMatrix;
  const spacialWeights: array of double): IMatrix;
var uPinv : IMatrix;
    x, y : integer;
    rootWeights : IMatrix;
    residual : IMatrix;
    a : IMatrix;
begin
     // #########################################################
     // #### Calculate the coefficients considering the spacial weights
     rootWeights := MatrixClass.Create(1, Length(spacialWeights));
     for y := 0 to rootWeights.Height - 1 do
         rootWeights[0, y] := sqrt(spacialWeights[y]);

     uPinv := fEigVecs.Clone;

     // elementwise multiply the elements of U with the spacial weights
     for x := 0 to uPinv.Width - 1 do
     begin
          uPinv.SetSubMatrix(x, 0, 1, uPinv.Height);
          uPinv.ElementWiseMult(rootWeights);
     end;
     uPinv.UseFullMatrix;
     if uPinv.PseudoInversionInPlace <> srOk then
        exit;

     // calculated weighted residual
     residual := mtx.Sub(fMeanElem);
     residual.ElementWiseMult(rootWeights);
     rootWeights := nil;

     // now calculate the projected elements in feature space
     a := uPinv.Mult(residual);
     uPinv := nil;

     // reconstruct example
     residual := fEigVecs.Mult(a);
     residual.AddInplace(fMeanElem);

     // ############################################################
     // #### Blend the input and the reconstructed image considering the spatial weights
     Result := MatrixClass.Create(1, mtx.Height);

     for y := 0 to Result.Height - 1 do
         Result[0, y] := mtx[0, y]*spacialWeights[y] + residual[0, y]*(1 - spacialWeights[y]);
end;

function TSubspaceUpdater.UpdateEigenspace(mtx: TDoubleMatrix;
  const spacialWeights: array of double; var EigVecs, EigVals, MeanElem,
  A: IMatrix): boolean;
var blendMtx : IMatrix;
begin
     // #######################################################
     // #### Execute the generic eigenspace update
     fEigVecs := EigVecs;
     fEigVals := EigVals;
     fMeanElem := MeanElem;
     fA := A;

     blendMtx := PrepareBlendExample(mtx, spacialWeights);
     Result := InternalUpdateEigenspace(blendMtx.GetObjRef, False);

     EigVecs := fEigVecs;
     EigVals := fEigVals;
     MeanElem := fMeanElem;
     A := fA;
end;

function TSubspaceUpdater.UpdateEigenspace(mtx: TDoubleMatrix; var EigVecs,
  EigVals, MeanElem, A: IMatrix): boolean;
begin
     // #######################################################
     // #### Execute the generic eigenspace update
     assert(Assigned(mtx), 'Error missing data');
     fEigVecs := EigVecs;
     fEigVals := EigVals;
     fMeanElem := MeanElem;
     fA := A;

     Result := InternalUpdateEigenspace(mtx, False);

     EigVecs := fEigVecs;
     EigVals := fEigVals;
     MeanElem := fMeanElem;
     A := fA;
end;

function TSubspaceUpdater.UpdateEigenspaceWeighted(mtx: TDoubleMatrix;
  var weights: TDoubleDynArray; var NumWeights : integer; var EigVecs, EigVals, MeanElem,
  A: IMatrix): boolean;
begin
     // #######################################################
     // #### Execute the generic eigenspace update
     fEigVecs := EigVecs;
     fEigVals := EigVals;
     fMeanElem := MeanElem;
     fWeights := Weights;
     fNumWeights := numWeights;
     fA := A;

     Result := InternalUpdateEigenspace(mtx, True);

     EigVecs := fEigVecs;
     EigVals := fEigVals;
     MeanElem := fMeanElem;
     A := fA;
     weights := fWeights;
     numWeights := fNumWeights;
end;

function TSubspaceUpdater.UpdateEigenspaceWeighted(mtx: TDoubleMatrix;
  var weights : TDoubleDynArray; var numWeights : integer; const spacialWeights: array of double; var EigVecs, EigVals, MeanElem,
  A: IMatrix): boolean;
var blendMtx : IMatrix;
begin
     // #######################################################
     // #### Execute the generic eigenspace update
     fEigVecs := EigVecs;
     fEigVals := EigVals;
     fMeanElem := MeanElem;
     fWeights := weights;
     fNumWeights := numWeights;
     fA := A;

     blendMtx := PrepareBlendExample(mtx, spacialWeights);
     assert(Assigned(blendMtx), 'No blend matrix created');
     Result := InternalUpdateEigenspace(blendMtx.GetObjRef, True);

     EigVecs := fEigVecs;
     EigVals := fEigVals;
     MeanElem := fMeanElem;
     A := fA;
     numWeights := fNumWeights;
     weights := fWeights;
end;


// ####################################################################
// #### Object Implementations
// ####################################################################

{ TIncrementalPCA }

procedure TIncrementalPCA.AddWeight(const weight: double);
begin
     if Length(fWeights) <= fNumWeights + 1 then
        SetLength(fWeights, Max(10, Min(Length(fWeights) + 1000, 2*fNumWeights)));

     fWeights[fNumWeights] := weight;
     inc(fNumWeights);
end;

constructor TIncrementalPCA.Create(const Data: TIncrPCAData);
begin
     inherited Create(Data.Props);

     fA := Data.A;
     fWeights := Data.Weights;
     fNumWeights := Length(fWeights);
     fNumEigenvectorsToKeep := Data.NumEigenVecs;
end;

constructor TIncrementalPCA.Create(KeepFlags: TPCASaveDataSet);
begin
     inherited Create(KeepFlags);

     fNumEigenvectorsToKeep := -1;
     fNumWeights := 0;
end;

destructor TIncrementalPCA.Destroy;
begin
     fA := nil;

     inherited;
end;

function TIncrementalPCA.GetA: TDoubleMatrix;
begin
     if not Assigned(fA) then
        fA := MatrixClass.Create;

     Result := fA.GetObjRef;
end;

// ###############################################
// #### Persistence functions

class function TIncrementalPCA.ClassIdentifier: String;
begin
     Result := 'IPCA';
end;

procedure TIncrementalPCA.DefineProps;
begin
     inherited;

     if Assigned(fA) then
        AddObject('A', fA.GetObjRef);
     AddIntProperty('NumEigenvectorsToKeep', fNumEigenvectorsToKeep);
     AddDoubleArr('IncrWeights', fWeights);
end;

function TIncrementalPCA.PropTypeOfName(const Name: string): TPropType;
begin
     if CompareText(Name, 'A') = 0
     then
         Result := ptObject
     else if CompareText(Name, 'NumEigenvectorsToKeep') = 0
     then
         Result := ptInteger
     else if CompareText(Name, 'IncrWeights') = 0
     then
         Result := ptDouble
     else
         Result := inherited PropTypeOfName(Name);
end;

procedure TIncrementalPCA.OnLoadDoubleArr(const Name: String;
  const Value: TDoubleDynArray);
begin
     if CompareText(Name, 'IncrWeights') = 0 then
     begin
          fNumWeights := Length(value);
          fWeights := Value;
     end
     else
         inherited;
end;

function TIncrementalPCA.OnLoadObject(const Name: String;
  obj: TBaseMathPersistence): boolean;
begin
     Result := True;
     if CompareText(Name, 'A') = 0
     then
         fA := obj as TDoubleMatrix
     else
         Result := inherited OnLoadObject(Name, Obj);
end;

procedure TIncrementalPCA.OnLoadIntProperty(const Name: String;
  Value: integer);
begin
     if CompareText(Name, 'NumEigenvectorsToKeep') = 0
     then
         fNumEigenvectorsToKeep := Value
     else
         inherited;
end;

function TIncrementalPCA.PCA(Examples: TDoubleMatrix; CutEps: double;
  IsRelativeValue: boolean): boolean;
var origKeepFlags : TPCASaveDataSet;
begin
     fNumWeights := 0;
     fWeights := nil;

     origKeepFlags := fKeepFlags;
     include(fKeepFlags, pcaMeanNormalizedData);

     Result := inherited PCA(Examples, CutEps, IsRelativeValue);

     if Result then
     begin
          fA := EigVecT.Mult(fMeanNormExamples);

          if not (pcaMeanNormalizedData in origKeepFlags) then
             fMeanNormExamples := nil;

          if not (pcaTransposedEigVec in origKeepFlags) then
             fEigVecsT := nil;

          fKeepFlags := origKeepFlags;
     end;
end;

function TIncrementalPCA.TemporalWeightPCA(Examples: TDoubleMatrix;
  CutEps: double; IsRelativeValue: boolean;
  const Weights: array of Double): boolean;
var origKeepFlags : TPCASaveDataSet;
begin
     SetLength(fWeights, Length(Weights));
     fNumWeights := Length(fWeights);
     SetLength(fWeights, Length(Weights));
     if Length(Weights) > 0 then
        Move(Weights[0], fWeights[0], sizeof(double)*Length(Weights));

     origKeepFlags := fKeepFlags;
     include(fKeepFlags, pcaMeanNormalizedData);

     Result := inherited TemporalWeightPCA(Examples, CutEps, IsRelativeValue, Weights);

     if Result then
     begin
          fA := EigVecT.Mult(fMeanNormExamples);

          if not (pcaMeanNormalizedData in origKeepFlags) then
             fMeanNormExamples := nil;

          if not (pcaTransposedEigVec in origKeepFlags) then
             fEigVecsT := nil;

          fKeepFlags := origKeepFlags;
     end;
end;

function TIncrementalPCA.UpdateEigenspace(mtx: TDoubleMatrix) : boolean;
var updater : TSubspaceUpdater;
begin
     Result := True;
     assert(mtx.Width = 1, 'Error only one example allowed');

     if not Assigned(fMeanElem) then
     begin
          fMeanElem := mtx.Clone;
          fEigVecs := MatrixClass.Create(1, mtx.Height);
          fA := MatrixClass.Create(1, 1);

          if pcaEigVals in fKeepFlags then
             fEigVals := MatrixClass.Create(1, 1);

          exit;
     end;

     // #######################################################
     // #### Execute the generic eigenspace update
     updater := TSubspaceUpdater.Create(fNumEigenvectorsToKeep);
     try
        updater.MatrixClass := MatrixClass;
        Result := updater.UpdateEigenspace(mtx, fEigVecs, fEigVals, fMeanElem, fA);
     finally
            updater.Free;
     end;
end;

function TIncrementalPCA.UpdateEigenspace(mtx: TDoubleMatrix;
  const spacialWeights: array of double): boolean;
var updater : TSubspaceUpdater;
begin
     Result := True;
     assert(mtx.Width = 1, 'Error only one example allowed');
     assert(mtx.Height = length(spacialWeights), 'Dimension error, spatial weights lenght differs from matrix height');

     if not Assigned(fMeanElem) then
     begin
          fMeanElem := mtx.Clone;
          fEigVecs := MatrixClass.Create(1, mtx.Height);
          fA := MatrixClass.Create(1, 1);

          if pcaEigVals in fKeepFlags then
             fEigVals := MatrixClass.Create(1, 1);

          exit;
     end;

     // ########################################################
     // #### Use the subspace updater:
     updater := TSubspaceUpdater.Create(fNumEigenvectorsToKeep);
     try
        updater.MatrixClass := MatrixClass;
        Result := updater.UpdateEigenspace(mtx, spacialWeights, fEigVecs, fEigVals, fMeanElem, fA);
     finally
            updater.Free;
     end;
end;

function TIncrementalPCA.UpdateEigenspaceWeighted(mtx: TDoubleMatrix;
  const weight: double; const spacialWeights: array of double): boolean;
var updater : TSubspaceUpdater;
begin
     Result := True;
     assert(mtx.Width = 1, 'Error only one example allowed');
     assert(mtx.Height = length(spacialWeights), 'Dimension error, spatial weights lenght differs from matrix height');

     AddWeight(weight);

     if not Assigned(fMeanElem) then
     begin
          fMeanElem := mtx.Clone;
          fEigVecs := MatrixClass.Create(1, mtx.Height);
          fA := MatrixClass.Create(1, 1);

          if pcaEigVals in fKeepFlags then
             fEigVals := MatrixClass.Create(1, 1);

          exit;
     end;

     // ########################################################
     // #### Use the subspace updater:
     updater := TSubspaceUpdater.Create(fNumEigenvectorsToKeep);
     try
        updater.MatrixClass := MatrixClass;
        Result := updater.UpdateEigenspaceWeighted(mtx, fWeights, fNumWeights, spacialWeights, fEigVecs, fEigVals, fMeanElem, fA);
     finally
            updater.Free;
     end;
end;

function TIncrementalPCA.UpdateEigenspaceWeighted(mtx: TDoubleMatrix;
  const weight: double): boolean;
var updater : TSubspaceUpdater;
begin
     Result := True;
     assert(mtx.Width = 1, 'Error only one example allowed');

     AddWeight(weight);

     if not Assigned(fMeanElem) then
     begin
          fMeanElem := MatrixClass.Create;
          fMeanElem.Assign(mtx, True);
          fEigVecs := MatrixClass.Create(1, mtx.Height);
          fA := MatrixClass.Create(1, 1);

          if pcaEigVals in fKeepFlags then
             fEigVals := MatrixClass.Create(1, 1);

          exit;
     end;

     // ########################################################
     // #### Use the subspace updater:
     updater := TSubspaceUpdater.Create(fNumEigenvectorsToKeep);
     try
        updater.MatrixClass := MatrixClass;
        Result := updater.UpdateEigenspaceWeighted(mtx, fWeights, fNumWeights, fEigVecs, fEigVals, fMeanElem, fA);
     finally
            updater.Free;
     end;
end;

{ TFastRobustIncrementalPCA }

procedure TFastRobustIncrementalPCA.AddWeight(idx : integer; const weight: double);
begin
     if Length(fWeights[idx]) <= fNumWeights[idx] + 1 then
        SetLength(fWeights[idx], Max(10, Min(Length(fWeights[idx]) + 1000, 2*fNumWeights[idx])));

     fWeights[idx][fNumWeights[idx]] := weight;
     inc(fNumWeights[idx]);
end;

procedure TFastRobustIncrementalPCA.BuildSubspaces(Examples: TDoubleMatrix; origKeepFlags : TPCASaveDataSet);
var i, idx : integer;
    j : integer;
    x : IMatrix;
    subT : IMatrix;
begin
     // main subspace
     fA := EigVecT.Mult(fMeanNormExamples);

     if not (pcaMeanNormalizedData in origKeepFlags) then
        fMeanNormExamples := nil;

     if not (pcaTransposedEigVec in origKeepFlags) then
        fEigVecsT := nil;

     // subsubspaces
     SetLength(fSubA, fProps.NumSubSubSpaces);
     for idx := 0 to Length(fSubEigVecs) - 1 do
     begin
          x := MatrixClass.Create(Examples.Width, fSubMeanElem[idx].Height);
          // project to feature space
          for j := 0 to Examples.Width - 1 do
          begin
               for i := 0 to x.Height - 1 do
                   x[j, i] := Examples[j, fSubItemIdx[idx][i]];

               x.SetSubMatrix(j, 0, 1, x.Height);
               x.SubInPlace(fSubMeanElem[idx]);
               x.UseFullMatrix;
          end;

          subT := fSubEigVecs[idx].Transpose;
          fSubA[idx] := subT.Mult(x);
     end;
end;

procedure TFastRobustIncrementalPCA.Clear;
begin
     inherited;

     if not fIsLearning then
     begin
          fA := nil;
          fSubA := nil;
          fWeights := nil;
          fNumWeights := nil;
          fSubEigVals := nil;
          fSubEigVecsT := nil;
     end;
end;

constructor TFastRobustIncrementalPCA.Create(KeepFlags: TPCASaveDataSet);
begin
     inherited Create(KeepFlags);
end;

constructor TFastRobustIncrementalPCA.Create(const Props: TFastRobustPCADataEx);
var i : integer;
begin
     inherited Create(props.RData);

     fA := props.A;
     fNumEigenvectorsToKeep := Props.NumEigenVecs;
     fSubA := Props.SubA;
     SetLength(fWeights, 1 + Length(Props.SubWeights));
     fWeights[0] := Props.Weights;
     SetLength(fNumWeights, Length(fWeights));
     for i := 0 to Length(fNumWeights) - 1 do
         fNumWeights[i] := Length(fWeights[i]);
     fSubA := Props.SubA;
     for i := 0 to Length(Props.SubWeights) - 1 do
         fWeights[1 + i] := Props.SubWeights[i];
end;

destructor TFastRobustIncrementalPCA.Destroy;
begin
     // inherited calls Clear!
     inherited;
end;

procedure TFastRobustIncrementalPCA.DefineProps;
var i : Integer;
begin
     inherited;

     if Assigned(fA) then
        AddObject('A', fA.GetObjRef);
     AddIntProperty('NumEigenvectorsToKeep', fNumEigenvectorsToKeep);
     if Length(fWeights) > 0 then
     begin
          BeginList('IPCAWeights', Length(fWeights));
          for i := 0 to Length(fWeights) - 1 do
              AddListDoubleArr(fWeights[i]);
          EndList;
     end;
     if Length(fNumWeights) > 0 then
        AddIntArr('IPCANumWeights', fNumWeights);
     if Length(fSubA) > 0 then
     begin
          BeginList('IPCASUBA', Length(fSubA));
          for i := 0 to Length(fSubA) - 1 do
              AddObject(fSubA[i].GetObjRef);
          EndList;
     end;
     if Length(fSubEigVals) > 0 then
     begin
          BeginList('IPCASUBEIGVALS', Length(fSubEigVals));
          for i := 0 to Length(fSubEigVals) - 1 do
              AddObject(fSubEigVals[i].GetObjRef);
          EndList;
     end;
end;

function TFastRobustIncrementalPCA.PropTypeOfName(
  const Name: string): TPropType;
begin
     if (CompareText(Name, 'IPCASUBEIGVALS') = 0) or (CompareText(Name, 'IPCASUBA') = 0) or
        (CompareText(Name, 'A') = 0)
     then
         Result := ptObject
     else if CompareText(Name, 'IPCAWeights') = 0
     then
         Result := ptObject
     else if CompareText(Name, 'NumEigenvectorsToKeep') = 0
     then
         Result := ptInteger
     else
         Result := inherited PropTypeOfName(Name);
end;


class function TFastRobustIncrementalPCA.ClassIdentifier: String;
begin
     Result := 'FRIPCA';
end;

procedure TFastRobustIncrementalPCA.OnLoadBeginList(const Name: String;
  count: integer);
begin
     fIdx := 0;
     if CompareText(Name, 'IPCAWeights') = 0 then
     begin
          fListType := llWeights;
          SetLength(fWeights, Count);
     end
     else if CompareText(Name, 'IPCASUBA') = 0 then
     begin
          fListType := llSubA;
          SetLength(fSubA, Count);
     end
     else if CompareText(Name, 'IPCASUBEIGVALS') = 0 then
     begin
          fListType := llSubEigVals;
          SetLength(fSubEigVals, Count);
     end
     else
         inherited;
end;

procedure TFastRobustIncrementalPCA.OnLoadEndList;
begin
     inherited;

     fListType := llNone;     
end;

procedure TFastRobustIncrementalPCA.OnLoadIntArr(const Name: String;
  const Value: TIntegerDynArray);
begin
     if CompareText(Name, 'IPCANumWeights') = 0
     then
         fNumWeights := Value
     else
         inherited;
end;

procedure TFastRobustIncrementalPCA.OnLoadListDoubleArr(
  const Value: TDoubleDynArray);
begin
     if fListType = llWeights then
     begin
          fWeights[fIdx] := Value;
          inc(fIdx);
     end
     else
         inherited;
end;

function TFastRobustIncrementalPCA.OnLoadObject(const Name: String;
  obj: TBaseMathPersistence): boolean;
begin
     Result := True;
     if CompareText(Name, 'A') = 0
     then
         fA := obj as TDoubleMatrix
     else
         Result := inherited OnLoadObject(Name, Obj);
end;

procedure TFastRobustIncrementalPCA.OnLoadIntProperty(const Name: String;
 Value: integer);
begin
     if CompareText(Name, 'NumEigenvectorsToKeep') = 0
     then
         fNumEigenvectorsToKeep := Value
     else
         inherited;
end;

function TFastRobustIncrementalPCA.OnLoadObject(
  Obj: TBaseMathPersistence): boolean;
begin
     Result := True;
     if fListType = llSubA then
     begin
          fSubA[fIdx] := obj as TDoubleMatrix;
          inc(fIdx);
     end
     else if fListType = llSubEigVals then
     begin
          fSubEigVals[fIdx] := obj as TDoubleMatrix;
          inc(fIdx);
     end
     else
         Result := inherited OnLoadObject(obj);
end;

function TFastRobustIncrementalPCA.PCA(Examples: TDoubleMatrix; CutEps: double;
  IsRelativeValue: boolean): boolean;
var origKeepFlags : TPCASaveDataSet;
begin
     Clear;

     fNumWeights := nil;
     fWeights := nil;

     origKeepFlags := fKeepFlags;
     include(fKeepFlags, pcaMeanNormalizedData);

     Result := inherited PCA(Examples, CutEps, IsRelativeValue);

     if Result then
     begin
          BuildSubspaces(Examples, origKeepFlags);
           
          fKeepFlags := origKeepFlags;
     end;
end;

function TFastRobustIncrementalPCA.TemporalWeightPCA(Examples: TDoubleMatrix;
  CutEps: double; IsRelativeValue: boolean;
  const Weights: array of Double): boolean;
var origKeepFlags : TPCASaveDataSet;
    i : Integer;
    sortWeights : TDoubleDynArray;
begin
     Clear;

     SetLength(fNumWeights, fProps.NumSubSubSpaces + 1);
     SetLength(fWeights, fProps.NumSubSubSpaces + 1);

     // build the list as it would be done incrementaly (sorted!)
     SetLength(sortWeights, Length(Weights));
     if Length(Weights) > 0 then
        Move(Weights[0], sortWeights[0], Length(Weights)*sizeof(double));
     QuickSort(sortWeights[0], sizeof(double), Length(Weights), @DoubleSortFunc);
     
     for i := 0 to Length(fNumWeights) - 1 do
     begin
          fWeights[i] := Copy(sortWeights, 0, Min(fNumEigenvectorsToKeep, Length(Weights)));
          fNumWeights[i] := Length(fWeights[i]);
     end;

     // ###############################################################
     // #### Build Sub + SubSubspaces
     origKeepFlags := fKeepFlags;
     include(fKeepFlags, pcaMeanNormalizedData);

     Result := inherited TemporalWeightPCA(Examples, CutEps, IsRelativeValue, Weights);

     // build the A matrices of all subspaces
     if Result then
     begin
          BuildSubspaces(Examples, origKeepFlags);

          fKeepFlags := origKeepFlags;
     end;
end;

function TFastRobustIncrementalPCA.UpdateEigenspace(
  mtx: TDoubleMatrix): boolean;
var spacialWeights : TDoubleDynArray;
begin
     spacialWeights := nil;
     Result := UpdateEigenspace(mtx, spacialWeights);
end;

function TFastRobustIncrementalPCA.UpdateEigenspace(mtx: TDoubleMatrix;
  const spacialWeights: array of double): boolean;
begin
     Result := UpdateEigenspaceWeighted(mtx, -1);
end;

function TFastRobustIncrementalPCA.UpdateEigenspaceWeighted(mtx: TDoubleMatrix;
  const weight: double; const spacialWeights: array of double): boolean;
var i, j : integer;
    subExample : IMatrix;
    pcaUpdater : TSubspaceUpdater;
    subSpaceWeights : TDoubleDynArray;
    example : IMatrix;
    eigValMtx : IMatrix;
begin
     Result := True;
     assert(mtx.Width = 1, 'Error only one example allowed');

     subSpaceWeights := nil;
     // #######################################################
     // #### Update sub sub spaces
     fIsLearning := True;
     try
        if not Assigned(fMeanElem) then
        begin
             SetLength(fSubItemIdx, fProps.NumSubSubSpaces);
             SetLength(fSubEigVecs, fProps.NumSubSubSpaces);
             SetLength(fSubMeanElem, fProps.NumSubSubSpaces);
             SetLength(fSubEigVecsT, fProps.NumSubSubSpaces);
             SetLength(fSubEigVals, fProps.NumSubSubSpaces);
             SetLength(fSubA, fProps.NumSubSubSpaces);
        end;
        for i := 0 to fProps.NumSubSubSpaces - 1 do
        begin
             if fSubItemIdx[i] = nil then
                fSubItemIdx[i] := GenerateSampleList(mtx.Height, i);

             SubExample := MatrixClass.Create(mtx.Width, Length(fSubItemIdx[i]));

             if Length(spacialWeights) > 0 then
             begin
                  SetLength(subSpaceWeights, SubExample.Height);
                  for j := 0 to Length(subSpaceWeights) - 1 do
                      subSpaceWeights[j] := spacialWeights[fSubItemIdx[i][j]];
             end;

             for j := 0 to Length(fSubItemIdx[i]) - 1 do
                 SubExample.SetRow(j, mtx, fSubItemIdx[i][j]);

             if weight > 0 then
                AddWeight(i + 1, weight);

             if not Assigned(fMeanElem) then
             begin
                  fSubMeanElem[i] := SubExample.Clone;
                  fSubEigVecs[i] := MatrixClass.Create(1, SubExample.Height);
                  fSubA[i] := MatrixClass.Create(1, 1);

                  if pcaEigVals in fKeepFlags then
                     fSubEigVals[i] := MatrixClass.Create(1, 1);
                  if pcaTransposedEigVec in fKeepFlags then
                     fSubEigVecsT[i] := MatrixClass.Create(1, mtx.Height);
             end
             else
             begin
                  pcaUpdater := TSubspaceUpdater.Create(fNumEigenvectorsToKeep);
                  try
                     pcaUpdater.MatrixClass := MatrixClass;
                     example := SubExample;
                     if Length(spacialWeights) <> 0 then
                        example := pcaUpdater.PrepareBlendExample(SubExample.GetObjRef, subSpaceWeights);

                     eigValMtx := nil;
                     if Length(fSubEigVals) > 0 then
                        eigValMtx := fSubEigVals[i];

                     if weight > 0
                     then
                         pcaUpdater.UpdateEigenspaceWeighted(Example.GetObjRef, fWeights[i + 1], fNumWeights[i + 1], fSubEigVecs[i], eigValMtx, fsubMeanElem[i], fSubA[i])
                     else
                         pcaUpdater.UpdateEigenspace(Example.GetObjRef, fSubEigVecs[i], eigValMtx, fsubMeanElem[i], fSubA[i]);

                     fSubEigVecsT[i] := nil;

                     if pcaTransposedEigVec in fKeepFlags then
                        fSubEigVecsT[i] := fSubEigVecs[i].Transpose;
                  finally
                         pcaUpdater.Free;
                  end;
             end;
        end;

        // #######################################################
        // #### Update full sub space
        if weight > 0 then
           AddWeight(0, weight);

        if not Assigned(fMeanElem) then
        begin
             fMeanElem := mtx.Clone;
             fEigVecs := MatrixClass.Create(1, mtx.Height);
             fA := MatrixClass.Create(1, 1);

             fEigVals := nil;
             if pcaEigVals in fKeepFlags then
                fEigVals := MatrixClass.Create(1, 1);

             exit;
        end;

        // #######################################################
        // #### Execute the generic eigenspace update
        pcaUpdater := TSubspaceUpdater.Create(fNumEigenvectorsToKeep);
        try
           pcaUpdater.MatrixClass := MatrixClass;
           example := mtx.Clone;
           if Length(spacialWeights) <> 0 then
              example := pcaUpdater.PrepareBlendExample(mtx, spacialWeights);

           if weight > 0
           then
               Result := pcaUpdater.UpdateEigenspaceWeighted(example.GetObjRef, fWeights[0], fNumWeights[0], fEigVecs, fEigVals, fMeanElem, fA)
           else
               Result := pcaUpdater.UpdateEigenspace(example.GetObjRef, fEigVecs, fEigVals, fMeanElem, fA);
        finally
               pcaUpdater.Free;
        end;
     finally
            fIsLearning := False;
     end;
end;

function TFastRobustIncrementalPCA.UpdateEigenspaceWeighted(mtx: TDoubleMatrix;
  const weight: double): boolean;
var spacialWeights : TDoubleDynArray;
begin
     spacialWeights := nil;
     Result := UpdateEigenspaceWeighted(mtx, weight, spacialWeights);
end;

function TFastRobustIncrementalPCA.getA: TDoubleMatrix;
begin
     Result := nil;
     if Assigned(fA) then
        Result := fA.GetObjRef;
end;

initialization
  RegisterMathIO(TIncrementalPCA);
  RegisterMathIO(TFastRobustIncrementalPCA);

end.