{*************************************}
{ Clustering algorithms               }
{                                     }
{ Copyright Steema Software 2012-2020 }
{ www.steema.com                      }
{*************************************}
unit TeeClustering;
{$I TeeDefs.inc}

interface

{$IFDEF NEXTGEN}
// Not possible to use our TFastList with Android / iOS
{$DEFINE TEEGENERICS}
{$ELSE}
{$DEFINE TEEFASTLIST}
{$ENDIF}

uses
  {$IFDEF MSWINDOWS}
  Windows,
  {$ENDIF}
  {$IFDEF TEEGENERICS}
  System.Generics.Collections,
  {$ENDIF}
  Classes;

type
  {$IFNDEF TEEGENERICS}

  {$IFDEF D16}
  PPointerList=^TPointerList;
  TPointerList = Array[0..(Maxint div 16) - 1] of Pointer;
  {$ENDIF}

  {$IFNDEF TEEFASTLIST}
  TFastList=TList;
  {$ELSE}
  TFastList=class
  private
    FCapacity : Integer;
    FCount : Integer;
    List   : PPointerList;

    function Get(Index:Integer):Pointer; {$IFDEF D9}inline;{$ENDIF}
    procedure Grow;
    procedure SetCapacity(const NewCapacity: Integer);
  public
    Destructor Destroy; override;

    procedure Add(const Item:Pointer);
    procedure Clear; virtual;
    procedure Delete(const Index:Integer);

    property Count:Integer read FCount;
    property Items[Index:Integer]:Pointer read Get; default;
  end;
  {$ENDIF}
  {$ENDIF}

  TCluster=class;

  TCluster=class({$IFDEF TEEGENERICS}TList<TCluster>{$ELSE}TFastList{$ENDIF})
  private
    Flag : Boolean;

    function Get(Index:Integer):TCluster; {$IFDEF D9}inline;{$ENDIF}
  protected
    OwnsData : Boolean;
  public
    Data  : Pointer;

    Constructor Create(AOwnsData:Boolean=False);
    Destructor Destroy; override;

    function AddItem(const AData:Pointer; const AOwnsData:Boolean=False):TCluster;
    procedure Clear; {$IFNDEF TEEGENERICS}override;{$ENDIF}

    property Item[Index:Integer]:TCluster read Get; default;
  end;

  TSingleArray=Array of Single;

  TBaseClustering=class;

  // Abstract Cluster Provider

  TClusterProvider=class {$IFDEF D11}abstract{$ENDIF}
  protected
    procedure BeginClustering; virtual;
    function CalcDistanceData(const A,B:Integer):Single; overload; virtual; abstract;
    function CalcDistanceData(const A,B:TCluster):Single; overload; virtual; abstract;
    function CalcDistanceData(const A:Integer; const B:TCluster):Single; overload; virtual; abstract;
    function CloneData(const AData:Pointer):Pointer; overload; virtual; abstract;
    function EmptyData:Pointer; {$IFDEF D9}inline;{$ENDIF}
    procedure EndClustering; virtual;
    function GetData(const Index:Integer):Pointer; virtual; abstract;
    procedure RecalcValues(const ACluster:TCluster); virtual; abstract;
    function SourceCount:Integer; virtual; abstract;
  end;

  TClusterEvaluation=record
    DaviesBouldin,
    Dunn : Double;
  end;

  TLinkageMethod=(lmNone, lmSingle, lmComplete, lmAverage, lmWard);

  // Abstract Cluster

  TBaseClustering=class(TPersistent)
  private
    FLinkage    : TLinkageMethod;
    FThreads    : Integer;
    FProvider   : TClusterProvider;

    SourceCount : Integer;
    Matrix      : Array of TSingleArray;

    procedure BuildMatrix;
    procedure CalcMatrix(Index:Integer);
  protected
    function CalcDistance(const A,B:TCluster):Single;
  public
    Clusters : TCluster;

    Constructor Create(const AProvider:TClusterProvider); virtual;
    Destructor Destroy; override;

    function Execute:TCluster; virtual;
    class function IsUndefined(const ACluster:TCluster):Boolean;
    class function Evaluation(AProvider:TClusterProvider; ACluster:TCluster):TClusterEvaluation;

    property Threads:Integer read FThreads write FThreads default 1;
  published
    property Linkage:TLinkageMethod read FLinkage write FLinkage default lmNone;
    property Provider:TClusterProvider read FProvider;
  end;

{$DEFINE MATRICES}

{$IFDEF MATRICES}
  PSingleList=^TSingleList;
  TSingleList = array[0..(MaxInt div 16) - 1] of Single;

  TFloatList=class
  private
    FCapacity : Integer;
    FCount : Integer;
    List   : PSingleList;

    function Get(Index:Integer):Single; {$IFDEF D9}inline;{$ENDIF}
    procedure Grow;
    procedure SetCapacity(const NewCapacity: Integer);
  public
    Destructor Destroy; override;

    procedure Add(const Item:Single);
    procedure Clear; virtual;
    procedure Delete(const Index:Integer);

    property Count:Integer read FCount;
    property Items[Index:Integer]:Single read Get; default;
  end;

  PMatrixItem=^TMatrixItem;
  TMatrixItem=record
    Value : Single;
    Next  : PMatrixItem;
  end;

  TMatrix=class(TFloatList)
  private
    Min : Single;
    MinIndex : Integer;
    NeedsRecalcMin : Boolean;

    procedure RecalcMin;
  public
    procedure Clear; override;
  end;

  TMatrixList=class({$IFDEF TEEGENERICS}TList<TMatrix>{$ELSE}TFastList{$ENDIF})
  private
  public
    procedure Clear; {$IFNDEF TEEGENERICS}override;{$ENDIF}
  end;
{$ENDIF}

  THierarchicalClustering=class(TBaseClustering)
  private
    FCount : Integer;

    {$IFDEF MATRICES}
    Matrices : TMatrixList;
    {$ENDIF}

    procedure BuildMatrix(Index:Integer);
  public
    Constructor Create(const AProvider:TClusterProvider); override;

    function Execute:TCluster; override;
  published
    property NumClusters:Integer read FCount write FCount default 5;
  end;

  TKMeansClustering=class(TBaseClustering)
  private
    FKCount  : Integer;
    FMaxIter : Integer;
  public
    Constructor Create(const AProvider:TClusterProvider); override;

    function Execute:TCluster; override;
  published
    property NumClusters:Integer read FKCount write FKCount default 5;
    property MaxIterations:Integer read FMaxIter write FMaxIter default 1000;
  end;

  TDistance=record
    Index : Integer;
    Value : Single;
  end;

  TDistances=Array of TDistance;

  TQTClustering=class(TBaseClustering)
  private
    FMinCount    : Integer;
    FMaxDiameter : Double;

    Distances    : Array of TDistances;
    procedure SortedDists(n:Integer);
  public
    Constructor Create(const AProvider:TClusterProvider); override;

    function Execute:TCluster; override;
  published
    property MinCount:Integer read FMinCount write FMinCount default 1;
    property MaxDiameter:Double read FMaxDiameter write FMaxDiameter;
  end;

implementation

uses
  Math, SysUtils,
  {$IFDEF FMX}
  FMXTee.Procs
  {$ELSE}
  TeeProcs
  {$ENDIF};

{$IFNDEF FMX}
{$DEFINE TEEPARALLEL}
{$ENDIF}

{ TCluster }

Constructor TCluster.Create(AOwnsData:Boolean=False);
begin
  inherited Create;
  OwnsData:=AOwnsData;
end;

Destructor TCluster.Destroy;
begin
  if OwnsData then
     FreeMem(Data);

  inherited;
end;

procedure TCluster.Clear;
var t : Integer;
begin
  for t:=0 to Count-1 do
      {$IFDEF AUTOREFCOUNT}
      Item[t].DisposeOf;
      {$ELSE}
      Item[t].Free;
      {$ENDIF}

  inherited;
end;

function TCluster.AddItem(const AData: Pointer; const AOwnsData:Boolean=False): TCluster;
begin
  result:=TCluster.Create(AOwnsData);
  result.Data:=AData;
  Add(result);
end;

function TCluster.Get(Index: Integer): TCluster;
begin
  result:=List[Index];
end;

{ TBaseClustering }

Constructor TBaseClustering.Create(const AProvider:TClusterProvider);
begin
  inherited Create;

  FThreads:=1;
  FProvider:=AProvider;

  Clusters:=TCluster.Create(True);
  Clusters.Data:=Provider.EmptyData;
end;

Destructor TBaseClustering.Destroy;
begin
  Clusters.Free;
  inherited;
end;

class function TBaseClustering.IsUndefined(const ACluster:TCluster):Boolean;
begin
  result:=ACluster.Flag;
end;

procedure TBaseClustering.CalcMatrix(Index:Integer);
var t   : Integer;
    tmp : TSingleArray;
begin
  SetLength(Matrix[Index], SourceCount);
  tmp:=Matrix[Index];

  for t:=Index+1 to SourceCount-1 do
      tmp[t]:=Provider.CalcDistanceData(Index, t);
end;

procedure TBaseClustering.BuildMatrix;
var t : Integer;
begin
  SetLength(Matrix, SourceCount);

  {$IFDEF TEEPARALLEL}
  if (FThreads>1) and (TTeeCPU.NumberOfProcessors>1) then
     TTeeCPU.ParallelFor(0, SourceCount-1, FThreads, CalcMatrix)
  else
  {$ENDIF}
  for t:=0 to SourceCount-1 do
      CalcMatrix(t);
end;

function TBaseClustering.Execute: TCluster;
begin
  Provider.BeginClustering;

  Clusters.Clear;

  if Assigned(Provider) then
     SourceCount:=Provider.SourceCount
  else
     SourceCount:=0;

  result:=Clusters;
end;

function TBaseClustering.CalcDistance(const A, B: TCluster): Single;
var
  tmpA,
  tmpB : Integer;

  procedure SingleLinkage;
  var t, tt : Integer;
      tmp   : Single;
  begin
    result:=MaxSingle;

    if (tmpA>0) and (tmpB>0) then
    begin
      for t:=0 to tmpA-1 do
          for tt:=0 to tmpB-1 do
          begin
            tmp:=Provider.CalcDistanceData(A[t], B[tt]);

            if tmp<result then
               result:=tmp;
          end
    end
    else
    if tmpA>0 then
    begin
      for t:=0 to tmpA-1 do
      begin
        tmp:=Provider.CalcDistanceData(A[t], B);

        if tmp<result then
           result:=tmp;
      end;
    end
    else
      for t:=0 to tmpB-1 do
      begin
        tmp:=Provider.CalcDistanceData(A, B[t]);

        if tmp<result then
           result:=tmp;
      end;
  end;

  procedure CompleteLinkage;
  var t, tt : Integer;
      tmp : Single;
  begin
    result:=MinSingle;

    if (tmpA>0) and (tmpB>0) then
    begin
      for t:=0 to tmpA-1 do
          for tt:=0 to tmpB-1 do
          begin
            tmp:=Provider.CalcDistanceData(A[t], B[tt]);

            if tmp>result then
               result:=tmp;
          end;
    end
    else
    if tmpA>0 then
    begin
      for t:=0 to tmpA-1 do
      begin
        tmp:=Provider.CalcDistanceData(A[t], B);

        if tmp>result then
           result:=tmp;
      end;
    end
    else
      for t:=0 to tmpB-1 do
      begin
        tmp:=Provider.CalcDistanceData(A, B[t]);

        if tmp>result then
           result:=tmp;
      end;
  end;

  procedure AverageLinkage;
  var t, tt : Integer;
  begin
    result:=0;

    if (tmpA>0) and (tmpB>0) then
    begin
      for t:=0 to tmpA-1 do
          for tt:=0 to tmpB-1 do
              result:=result+Provider.CalcDistanceData(A[t], B[tt]);

      result:=result/(tmpA*tmpB);
    end
    else
    if tmpA>0 then
    begin
      for t:=0 to tmpA-1 do
          result:=result+Provider.CalcDistanceData(A[t], B);

      result:=result/tmpA;
    end
    else
    begin
      for t:=0 to tmpB-1 do
          result:=result+Provider.CalcDistanceData(A, B[t]);

      result:=result/tmpB;
    end;
  end;

  procedure WardLinkage;

    function ErrorSumOfSquares(const ACluster:TCluster):Double;
    var t : Integer;
    begin
      result:=0;

      for t:=0 to ACluster.Count-1 do
          result:=result+Sqr(Provider.CalcDistanceData(ACluster[t], ACluster));
    end;

  var tmp : TCluster;
      t   : Integer;
  begin

    tmp:=TCluster.Create(True);
    tmp.Data:=Provider.EmptyData;

    if tmpA=0 then
       tmp.Add(A)
    else
    for t:=0 to tmpA-1 do
        tmp.Add(A[t]);

    if tmpB=0 then
       tmp.Add(B)
    else
    for t:=0 to tmpB-1 do
        tmp.Add(B[t]);

    Provider.RecalcValues(tmp);

    result:=ErrorSumOfSquares(tmp);

    //TFastList(tmp).Clear;
    tmp.Free;

    if (tmpA>0) and (tmpB>0) then
       result:=result-(ErrorSumOfSquares(A)+ErrorSumOfSquares(B))
    else
    if tmpA>0 then
       result:=result-ErrorSumOfSquares(A)
    else
       result:=result-ErrorSumOfSquares(B);
  end;

begin
  if FLinkage=lmNone then
     result:=Provider.CalcDistanceData(A,B)
  else
  begin
    tmpA:=A.Count;
    tmpB:=B.Count;

    if (tmpA=0) and (tmpB=0) then
       result:=Provider.CalcDistanceData(A,B)
    else
    begin
      case FLinkage of
        lmSingle: SingleLinkage;
      lmComplete: CompleteLinkage;
       lmAverage: AverageLinkage;
      else
        WardLinkage;
      end;
    end;
  end;
end;

class function TBaseClustering.Evaluation(AProvider:TClusterProvider; ACluster:TCluster): TClusterEvaluation;

  function ClusterScatter(Index:Integer):Double;
  var tmp : TCluster;
      tmpSum : Double;
      t      : Integer;
  begin
    tmpSum:=0;

    tmp:=ACluster[Index];

    for t:=0 to tmp.Count-1 do
        tmpSum:=tmpSum+Sqr(AProvider.CalcDistanceData(tmp[t], tmp));

    result:=Sqrt(tmpSum/tmp.Count);
  end;

  function ClustersScatter(const A,B:TCluster):Double;
  var t,tt : Integer;
      tmpSum : Double;
  begin
    tmpSum:=0;

    for t:=0 to A.Count-1 do
        for tt:=0 to B.Count-1 do
            tmpSum:=tmpSum+Sqr(AProvider.CalcDistanceData(A[t], B[tt]));

    result:=Sqrt(tmpSum/(A.Count*B.Count));
  end;

var t, tt : Integer;
    tmpScatter,
    tmpMax,
    tmp,
    tmpSum : Double;
    //tmpData : Pointer;
begin
  if ACluster.Count=0 then
  begin
    result.DaviesBouldin:=0;
    result.Dunn:=0;
  end
  else
  begin

    tmpSum:=0;

    for t:=0 to ACluster.Count-1 do
    begin
      //tmpData:=ACluster[t].Data;
      tmpScatter:=ClusterScatter(t);

      tmpMax:=0;

      for tt:=0 to ACluster.Count-1 do
      if tt<>t then
      begin
        tmp:=(tmpScatter+ClusterScatter(tt))/ClustersScatter(ACluster[t],ACluster[tt]);

        if tmp>tmpMax then
           tmpMax:=tmp;
      end;

      tmpSum:=tmpSum+tmpMax;
    end;

    result.DaviesBouldin:=tmpSum/ACluster.Count;

    result.Dunn:=0;
  end;
end;

{ THierarchicalClustering }

Constructor THierarchicalClustering.Create(const AProvider: TClusterProvider);
begin
  inherited;
  FCount:=5;
end;

procedure THierarchicalClustering.BuildMatrix(Index:Integer);
var tmpMatrix : TMatrix;
    t         : Integer;
    tmpValue  : Single;
begin
  tmpMatrix:=Matrices[Index];

  tmpMatrix.SetCapacity(SourceCount-Index);

  tmpMatrix.Min:=MaxSingle;
  tmpMatrix.MinIndex:=-1;

  for t:=Index+1 to SourceCount-1 do
  begin
    tmpValue:=Provider.CalcDistanceData(Index,t);
    tmpMatrix.Add(tmpValue);

    if tmpValue<tmpMatrix.Min then
    begin
      tmpMatrix.Min:=tmpValue;
      tmpMatrix.MinIndex:=tmpMatrix.Count-1;
    end;
  end;
end;

function THierarchicalClustering.Execute:TCluster;

  {$IFNDEF MATRICES}
  // Brute force O(n^2)
  // TODO: Speed opt, multithread candidate
  procedure FindNearestPair(out A,B:Integer); register;
  var tmpDistance : Single;
      t,tt : Integer;
  begin
    tmpDistance:=MaxSingle;

    for t:=0 to SourceCount-1 do
        for tt:=t+1 to SourceCount-1 do
          if tmpDistance>Matrix[t,tt] then
          begin
            tmpDistance:=Matrix[t,tt];
            A:=t;
            B:=tt;
          end;
  end;
  {$ELSE}
  procedure FindNearestPairMatrices(out A,B:Integer); register;
  var tmpDistance : Single;
      t : Integer;
  begin
    tmpDistance:=MaxSingle;

    for t:=0 to SourceCount-1 do
    with TMatrix(Matrices.List[t]) do
    if tmpDistance>Min then
    begin
      tmpDistance:=Min;
      A:=t;
      B:=MinIndex+A+1;
    end;
  end;
  {$ENDIF}

  procedure InitAllClusters;
  var t : Integer;
  begin
    for t:=0 to SourceCount-1 do
        Clusters.AddItem(Provider.GetData(t));
  end;

  {$IFDEF MATRICES}
  procedure BuildMatrices;
  var t : Integer;
  begin
    for t:=0 to SourceCount-1 do
        Matrices.Add(TMatrix.Create);

    {$IFDEF TEEPARALLEL}
    if (FThreads>1) and (TTeeCPU.NumberOfProcessors>1) then
       TTeeCPU.ParallelFor(0, SourceCount-1, FThreads, BuildMatrix)
    else
    {$ENDIF}
    for t:=0 to SourceCount-1 do
        BuildMatrix(t);

  end;
  {$ENDIF}

var a,b,
    t  : Integer;
    tmpCluster,
    CA : TCluster;
    {$IFDEF MATRICES}
    tmpMatrixA  : PSingleList;
    tmpVal      : Single;
    {$ELSE}
    tt : Integer;
    {$ENDIF}
begin
  result:=inherited Execute;

  if SourceCount>1 then
  begin
    InitAllClusters;

    {$IFDEF MATRICES}
    Matrices:=TMatrixList.Create;
    BuildMatrices;
    {$ELSE}
    BuildMatrix;
    {$ENDIF}

    repeat
      {$IFDEF MATRICES}
      FindNearestPairMatrices(a,b);
      {$ELSE}
      FindNearestPair(a,b);
      {$ENDIF}

      CA:=Clusters[a];

      // Create new Cluster to group "a" and "b"
      Clusters.List[a]:=TCluster.Create(True);

      tmpCluster:=Clusters[a]; // share "a"

      with tmpCluster do
      begin
        Data:=Provider.EmptyData;

        Add(CA);
        Add(Clusters[b]);
      end;

      Provider.RecalcValues(tmpCluster);

      Dec(SourceCount);

      Clusters.Delete(b);

      {$IFNDEF MATRICES}
      // Remove Matrix "b" rows and columns
      for t:=b to SourceCount-1 do
          Matrix[t]:=Matrix[t+1];

      for t:=0 to SourceCount-1 do
          for tt:=b to SourceCount-1 do
              Matrix[t,tt]:=Matrix[t,tt+1];

      {$ELSE}
      // Matrices

      for t:=0 to b-1 do
      with TMatrix(Matrices[t]) do
      begin
        NeedsRecalcMin:=((b-t-1)=MinIndex);

        Delete(b-t-1);

        if (b-t-1)<MinIndex then
           Dec(MinIndex);
      end;

      {$IFDEF AUTOREFCOUNT}
      TMatrix(Matrices[b]).DisposeOf;
      {$ELSE}
      TMatrix(Matrices[b]).Free;
      {$ENDIF}

      Matrices.Delete(b);
      {$ENDIF}

      // Recalc Matrix new "a" rows and columns
      {$IFDEF MATRICES}
      tmpMatrixA:=TMatrix(Matrices[a]).List;

      for t:=a+1 to SourceCount-1 do
          tmpMatrixA[t-a-1]:=CalcDistance(tmpCluster, Clusters[t]);

      TMatrix(Matrices[a]).RecalcMin;

      for t:=0 to a-1 do
      begin
        tmpVal:=CalcDistance(Clusters[t], tmpCluster);

        with TMatrix(Matrices.List[t]) do
        begin
          List[a-t-1]:=tmpVal;

          if not NeedsRecalcMin then
             if tmpVal<Min then
             begin
               Min:=tmpVal;
               MinIndex:=a-t-1;
             end;
        end;
      end;

      for t:=0 to Max(a,b)-1 do
      with TMatrix(Matrices.List[t]) do
      if NeedsRecalcMin then
         RecalcMin;

      {$ELSE}
      for t:=a+1 to SourceCount-1 do
          Matrix[a,t]:=CalcDistance(tmpCluster, Clusters[t]);

      for t:=0 to a-1 do
          Matrix[t,a]:=CalcDistance(Clusters[t], tmpCluster);
      {$ENDIF}

    until SourceCount<=FCount;

    {$IFDEF MATRICES}
    Matrices.Free;
    {$ELSE}
    Matrix:=nil;
    {$ENDIF}

    Provider.RecalcValues(Clusters);
  end;

  Provider.EndClustering;
end;

{ TKMeansClustering }

Constructor TKMeansClustering.Create;
begin
  inherited;
  FKCount:=5;
  FMaxIter:=1000;
end;

function TKMeansClustering.Execute: TCluster;

  procedure RemoveEmptyClusters;
  var t : Integer;
  begin
    t:=0;

    repeat
      if Clusters[t].Count=0 then
      begin
        {$IFDEF AUTOREFCOUNT}
        Clusters[t].DisposeOf;
        {$ELSE}
        Clusters[t].Free;
        {$ENDIF}

        Clusters.Delete(t);
      end
      else
        Inc(t);

    until t=Clusters.Count;
  end;

  procedure InitRandomClusters;
  var Randoms : Array of Integer;

    function UniqueRandom(AMax:Integer):Integer;
    var t : Integer;
        Found : Boolean;
    begin
      repeat
        result:=Random(SourceCount);

        Found:=False;

        for t:=0 to AMax-1 do
            if Randoms[t]=result then
            begin
              Found:=True;
              break;
            end;

      until not Found;
    end;

  var t : Integer;
  begin
    SetLength(Randoms,FKCount);

    // Init random K-Clusters:
    for t:=0 to FKCount-1 do
        Randoms[t]:=UniqueRandom(t-1);

    for t:=0 to FKCount-1 do
        Clusters.AddItem(Provider.CloneData(Provider.GetData(Randoms[t])), True);

    Randoms:=nil;
  end;

  {$DEFINE TEEOPT}

type
  TIntegerArray=Array of Integer;

var t,tt,
    tmpI,
    Changes,
    tmpOld,
    tmpOldK,
    Iterations : Integer;
    tmp,tmp2 : Single;

    {$IFDEF TEEOPT}
    ClusterIndex : TIntegerArray;
    ClusterItemIndex : Array of TIntegerArray;
    tmpList : TIntegerArray;
    {$ENDIF}
begin
  result:=inherited Execute;

  {$IFNDEF D8}
  tmpList:=nil;
  {$ENDIF}
  
  if SourceCount>0 then
  begin
    InitRandomClusters;

    {$IFDEF TEEOPT}
    SetLength(ClusterIndex,SourceCount);
    SetLength(ClusterItemIndex,FKCount,SourceCount);

    for t:=0 to FKCount-1 do
        for tt:=0 to SourceCount-1 do
            ClusterItemIndex[t,tt]:=-1;

    for t:=0 to SourceCount-1 do
        ClusterIndex[t]:=-1;
    {$ENDIF}

    // Loop:
    Iterations:=0;

    repeat
      Changes:=0;

      for t:=0 to SourceCount-1 do
      begin
        // Find nearest K-Cluster:

        tmp:=MaxSingle;
        tmpI:=-1;

        for tt:=0 to FKCount-1 do
        begin
          tmp2:=Provider.CalcDistanceData(t, Clusters[tt]);

          if tmp2=0 then
          begin
            tmpI:=tt;
            break;
          end
          else
          if tmp2<tmp then
          begin
            tmp:=tmp2;
            tmpI:=tt;
          end;
        end;

        // Find existing K-Cluster:

        tmpOld:=-1;

        if Iterations>0 then
        begin
          {$IFDEF TEEOPT}
          tmpOld:=ClusterIndex[t];
          tmpOldK:=ClusterItemIndex[tmpOld,t];
          {$ELSE}
          tmpOldK:=Clusters[tmpI].IndexOfData(tmpData);

          if tmpOldK<>-1 then
             tmpOld:=tmpI
          else
          for tt:=0 to FKCount-1 do
          if tt<>tmpI then
          begin
            tmpOldK:=Clusters[tt].IndexOfData(tmpData);

            if tmpOldK<>-1 then
            begin
              tmpOld:=tt;
              break;
            end;
          end;
          {$ENDIF}
        end
        else
          tmpOldK:=-1;

        // If changing Cluster:

        if tmpOld<>tmpI then
        begin
          Inc(Changes);

          // Add to new cluster:
          Clusters[tmpI].AddItem(Provider.GetData(t));

          // Remove from old cluster:
          if tmpOld<>-1 then
          with Clusters[tmpOld] do
          begin
            {$IFDEF AUTOREFCOUNT}
            Item[tmpOldK].DisposeOf;
            {$ELSE}
            Item[tmpOldK].Free;
            {$ENDIF}

            Delete(tmpOldK);

            {$IFDEF TEEOPT}
            // Update indexes:
            tmpList:=ClusterItemIndex[tmpOld];

            for tt:=0 to SourceCount-1 do
                if tmpList[tt]>tmpOldK then
                   Dec(tmpList[tt]);
            {$ENDIF}
          end;

          {$IFDEF TEEOPT}
          ClusterIndex[t]:=tmpI;
          ClusterItemIndex[tmpI,t]:=Clusters[tmpI].Count-1;
          {$ENDIF}
        end;
      end;

      // Recalculate K clusters "centers":

      for t:=0 to FKCount-1 do
          Provider.RecalcValues(Clusters[t]);

      Inc(Iterations);

    until (Changes=0) or (Iterations>=MaxIterations);

    {$IFDEF TEEOPT}
    ClusterIndex:=nil;
    ClusterItemIndex:=nil;
    {$ENDIF}

    RemoveEmptyClusters;
  end;

  Provider.EndClustering;
end;

{ TQTClustering }

Constructor TQTClustering.Create;
begin
  inherited;
  FMinCount:=1;
  FMaxDiameter:=1;
end;

procedure TQTClustering.SortedDists(n:Integer);
var
  CurrentDist : TDistances;

  procedure PrivateSort(const l,r:Integer);
  var i : Integer;
      j : Integer;
      x : Integer;
      tmp : TDistance;
  begin
    i:=l;
    j:=r;
    x:=(i+j) shr 1;

    while i<j do
    begin
      while CurrentDist[i].Value<CurrentDist[x].Value do inc(i);
      while CurrentDist[x].Value<CurrentDist[j].Value do dec(j);

      if i<j then
      begin
        tmp:=CurrentDist[i];
        CurrentDist[i]:=CurrentDist[j];
        CurrentDist[j]:=tmp;

        if i=x then x:=j else if j=x then x:=i;
      end;

      if i<=j then
      begin
        inc(i);
        dec(j)
      end;
    end;

    if l<j then PrivateSort(l,j);
    if i<r then PrivateSort(i,r);
  end;

var t : Integer;
begin
  CurrentDist:=Distances[n];

  for t:=0 to n-1 do
  with CurrentDist[t] do
  begin
    Value:=Matrix[t,n];
    Index:=t;
  end;

  for t:=n+1 to SourceCount-1 do
  with CurrentDist[t] do
  begin
    Value:=Matrix[n,t];
    Index:=t;
  end;

  CurrentDist[n].Value:=MaxSingle;
  CurrentDist[n].Index:=n;

  PrivateSort(0,SourceCount-1);
end;

function TQTClustering.Execute:TCluster;
var n,t,
    tmpDistIndex : Integer;
    tmp,tmpRes,tmp2 : Integer;
    NotConsider : Array of Boolean;
    HasCandidate : Boolean;
    tmpNew : TCluster;
begin
  result:=inherited Execute;

  if SourceCount>0 then
  begin
    BuildMatrix;

    SetLength(Distances, SourceCount, SourceCount);

    {$IFDEF TEEPARALLEL}
    if (FThreads>1) and (TTeeCPU.NumberOfProcessors>1) then
       TTeeCPU.ParallelFor(0, SourceCount-1, FThreads, SortedDists)
    else
    {$ENDIF}
    for t:=0 to SourceCount-1 do
        SortedDists(t);

    Matrix:=nil;

    SetLength(NotConsider, SourceCount);

    repeat

      // Find bigger cluster (the one with max diameter)

      tmp:=0;
      tmpRes:=-1;

      for n:=0 to SourceCount-1 do
      if not NotConsider[n] then
      begin
        // Count sub items within diameter:

        tmp2:=0;

        for t:=0 to SourceCount-1 do
        begin
          tmpDistIndex:=Distances[n,t].Index;

          if (not NotConsider[tmpDistIndex]) and (tmpDistIndex<>n) then
          begin
            if Distances[n,t].Value<MaxDiameter then
               Inc(tmp2)
            else
               break;
          end;
        end;

        // Remember maximum

        if tmp2>tmp then
        begin
          tmp:=tmp2;
          tmpRes:=n;
        end;

      end;

      if tmpRes<>-1 then
      begin
        HasCandidate:= (tmp >= MinCount);

        if HasCandidate then
        begin
          NotConsider[tmpRes]:=True;

          // Create a new cluster:

          tmpNew:=Clusters.AddItem(Provider.EmptyData,True);

          // Add main point:

          tmpNew.AddItem(Provider.GetData(tmpRes));

          // Add items:

          for t:=0 to SourceCount-1 do
          begin
            tmpDistIndex:=Distances[tmpRes,t].Index;

            if (not NotConsider[tmpDistIndex]) and (tmpDistIndex<>tmpRes) then
            begin
              if Distances[tmpRes,t].Value<MaxDiameter then
              begin
                NotConsider[tmpDistIndex]:=True;

                tmpNew.AddItem(Provider.GetData(tmpDistIndex));
              end
              else
                break;
            end;
          end;

          Provider.RecalcValues(tmpNew);
        end;

      end
      else
        HasCandidate:=False;

    until not HasCandidate;

    Distances:=nil;

    // Count undefined items:

    tmp:=0;
    for t:=0 to SourceCount-1 do
        if not NotConsider[t] then
           Inc(tmp);

    // Add undefined items to new cluster, if any:

    if tmp>0 then
    with Clusters.AddItem(Provider.EmptyData, True) do
    begin
      Flag:=True;

      for t:=0 to SourceCount-1 do
          if not NotConsider[t] then
             AddItem(Provider.GetData(t));
    end;

    // Finish

    NotConsider:=nil;
  end;

  Provider.EndClustering;
end;

{ TClusterProvider }

procedure TClusterProvider.BeginClustering;
begin
end;

procedure TClusterProvider.EndClustering;
begin
end;

function TClusterProvider.EmptyData:Pointer;
begin
  result:=CloneData(nil);
end;

{$IFDEF MATRICES}
{ TMatrixList }

procedure TMatrixList.Clear;
var t : Integer;
begin
  for t:=0 to Count-1 do
      TMatrix(List[t]).Free;

  inherited;
end;

{ TMatrix }

procedure TMatrix.Clear;
begin
  Min:=0;
  inherited;
end;

procedure TMatrix.RecalcMin;
var t : Integer;
begin
  Min:=MaxSingle;
  MinIndex:=-1;

  for t:=0 to Count-1 do
      if List[t]<Min then
      begin
        Min:=List[t];
        MinIndex:=t;
      end;

  NeedsRecalcMin:=False;
end;

{$ENDIF}

{$IFNDEF TEEGENERICS}

{$IFDEF TEEFASTLIST}
{ TFastList }

procedure TFastList.Add(const Item: Pointer);
begin
  if FCount=FCapacity then
     Grow;

  List[FCount]:=Item;
  Inc(FCount);
end;

procedure TFastList.Clear;
begin
  {$IFDEF POSIX}
  if List<>nil then
  {$ENDIF}
     ReallocMem(List, 0);

  FCapacity:=0;
  FCount:=0;
end;

procedure TFastList.SetCapacity(const NewCapacity: Integer);
begin
  if NewCapacity <> FCapacity then
  begin
    ReallocMem(List, NewCapacity * SizeOf(Pointer));
    FCapacity := NewCapacity;
  end;
end;

procedure TFastList.Delete(const Index: Integer);
begin
  Dec(FCount);
  if Index < FCount then
    System.Move(List^[Index + 1], List^[Index], (FCount - Index) * SizeOf(Pointer));
end;

Destructor TFastList.Destroy;
begin
  Clear;
  inherited;
end;

procedure TFastList.Grow;
var Delta: Integer;
begin
  if FCapacity > 64 then
     Delta := FCapacity div 4
  else
  if FCapacity > 8 then
     Delta := 16
  else
     Delta := 4;

  SetCapacity(FCapacity + Delta);
end;

function TFastList.Get(Index: Integer): Pointer;
begin
  result:=List[Index];
end;
{$ENDIF}
{$ENDIF}

{$IFDEF MATRICES}
{ TFloatList }

procedure TFloatList.Add(const Item: Single);
begin
  if FCount=FCapacity then
     Grow;
  List^[FCount]:=Item;
  Inc(FCount);
end;

procedure TFloatList.Clear;
begin
  {$IFDEF POSIX}
  if List<>nil then
  {$ENDIF}
     ReallocMem(List, 0);

  FCount:=0;
end;

procedure TFloatList.Delete(const Index: Integer);
begin
  Dec(FCount);
  if Index < FCount then
     System.Move(List^[Index + 1], List^[Index], (FCount - Index) * SizeOf(Single));
end;

Destructor TFloatList.Destroy;
begin
  Clear;
  inherited;
end;

function TFloatList.Get(Index: Integer): Single;
begin
  result:=List[Index];
end;

procedure TFloatList.Grow;
var Delta: Integer;
begin
  if FCapacity > 64 then
     Delta := FCapacity div 4
  else
  if FCapacity > 8 then
     Delta := 16
  else
     Delta := 4;

  SetCapacity(FCapacity + Delta);
end;

procedure TFloatList.SetCapacity(const NewCapacity: Integer);
begin
  if NewCapacity <> FCapacity then
  begin
    ReallocMem(List, NewCapacity * SizeOf(Single));
    FCapacity := NewCapacity;
  end;
end;
{$ENDIF}

end.
