
unit uDistCorrMath;

{$MODE Delphi}

interface

uses uVecMat;

function DistCorr(n,p,q:integer; var x,y:datamatrix):real;
{n:number of samples; x:n*p; y:n*q}

procedure DistCoh(n,p,q:integer; // input
                  var xReData,xImData,yReData,yImData:datamatrix; // input
                  var TotCoh,ReCoh,ImCoh:real // output
                  );

implementation

var dx,dy,xReDist,xImDist,yReDist,yImDist:distmatrix;




const eps=1e-5;




procedure Data2Dist(n,p:integer; var x:datamatrix; var d:distmatrix);
var i,j,k:integer; r:extended;
begin
  for i:=1 to n do begin
    d[i,i]:=0;
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(x[i,k]-x[j,k]);
      r:=sqrt(r);
      d[i,j]:=r;
      d[j,i]:=r;
    end;
  end;
end;




procedure SzekelyDoubleCenter(n:integer; var d:distmatrix);
var i,j:integer; t,a:real; r,c:nvector;
begin
  for i:=1 to n do begin
    t:=0;
    for j:=1 to n do
      t:=t+d[i,j];
    r[i]:=t/n;
  end; // row mean

  for j:=1 to n do begin
    t:=0;
    for i:=1 to n do
      t:=t+d[i,j];
    c[j]:=t/n;
  end; // column mean; this is stupid since symmetric matrix has same row and col mean

  a:=0;
  for i:=1 to n do
    a:=a+r[i];
  a:=a/n; // global mean

  for i:=1 to n do
    for j:=1 to n do
      d[i,j]:=d[i,j]-r[i]-c[j]+a; // Szekely biased double centering
end;




procedure VarCovar(n:integer; var xy,xx,yy:real; var dx,dy:distmatrix);
var i,j:integer;
begin
  xy:=0;
  xx:=0;
  yy:=0;
  for i:=1 to n do
    for j:=1 to n do begin
      xy:=xy+dx[i,j]*dy[i,j];
      xx:=xx+dx[i,j]*dx[i,j];
      yy:=yy+dy[i,j]*dy[i,j];
    end; // distance var & covar ; could be more efficient due to symmetric matrix
  xy:=xy/sqr(n);
  xx:=xx/sqr(n);
  yy:=yy/sqr(n);
  if xx<eps then xx:=0;
  if yy<eps then yy:=0;
  if sqrt(xx*yy)<eps then begin
    xx:=0;
    yy:=0;
  end;
end;




function DistCorr(n,p,q:integer; var x,y:datamatrix):real;
{n:number of samples; x:n*p; y:n*q}
var xy,xx,yy:real;
begin
  Data2Dist(n,p,x,dx);
  Data2Dist(n,q,y,dy);

  SzekelyDoubleCenter(n,dx);
  SzekelyDoubleCenter(n,dy);

  VarCovar(n,xy,xx,yy,dx,dy);

  DistCorr:=xy/sqrt(xx*yy);
end;




procedure CmplxData2Dist(n,p:integer;
                         var xre,xim:datamatrix; var d:distmatrix);
var i,j,k:integer; r:extended;
begin
  for i:=1 to n do begin
    d[i,i]:=0;
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(xre[i,k]-xre[j,k])+sqr(xim[i,k]-xim[j,k]);
      r:=sqrt(r);
      d[i,j]:=r;
      d[j,i]:=r;
    end;
  end;
end;




procedure DistCoh(n,p,q:integer; // input
                  var xReData,xImData,yReData,yImData:datamatrix; // input
                  var TotCoh,ReCoh,ImCoh:real // output
                  );
var zxryr,zxryi,zxiyr,zxiyi,zden:real;
    totcovar,totvarx,totvary,fre,fim,ftot,t1,t2:real;
begin
  CmplxData2Dist(n,p,xredata,ximdata,xredist);
  CmplxData2Dist(n,q,yredata,yimdata,yredist);
  SzekelyDoubleCenter(n,xredist);
  SzekelyDoubleCenter(n,yredist);
  VarCovar(n,totcovar,totvarx,totvary,xredist,yredist);

  Data2Dist(n,p,xredata,xredist);
  Data2Dist(n,p,ximdata,ximdist);
  Data2Dist(n,q,yredata,yredist);
  Data2Dist(n,q,yimdata,yimdist);

  SzekelyDoubleCenter(n,xredist);
  SzekelyDoubleCenter(n,ximdist);
  SzekelyDoubleCenter(n,yredist);
  SzekelyDoubleCenter(n,yimdist);

  VarCovar(n,zxryr,t1,t2,xReDist,yReDist);
  VarCovar(n,zxryi,t1,t2,xReDist,yImDist);
  VarCovar(n,zxiyi,t1,t2,xImDist,yImDist);
  VarCovar(n,zxiyr,t1,t2,xImDist,yReDist);

  if (totvarx<eps) or (totvary<eps) then begin
    totcoh:=0;
    recoh:=0;
    imcoh:=0;
    exit;
  end;

  fre:=zxryr+zxiyi;
  fim:=zxryi+zxiyr;
  ftot:=fre+fim;
  if ftot<eps then ftot:=1;
  fre:=fre/ftot;
  fim:=fim/ftot;

  zden:=sqrt(totvarx*totvary);
  totcoh:=totcovar/zden;
  recoh:=fre*totcoh;
  imcoh:=fim*totcoh;

end;




end.

