
unit uSimCorrMath;

{$MODE Delphi}

interface

uses uVecMat;

function SimCorr(n,p,q:integer; alpha,sx,sy:real;
                 var x,y:datamatrix;
                 var xy,xx,yy:real
                 ):real;
{n:number of samples; x:n*p; y:n*q; alpha,sx,sy>0:parameters}

procedure SimCoh(n,p,q:integer; alpha,sx,sy:real; // input
                 var xReData,xImData,yReData,yImData:datamatrix; // input
                 var TotCoh,ReCoh,ImCoh,TotCovar:real // output
                 );

implementation

uses math;

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




const eps=1e-5;




procedure Data2Sim(n,p:integer; alpha,s:real; var x:datamatrix; var d:distmatrix);
var i,j,k:integer; r:extended;
begin
  for i:=1 to n do begin
    d[i,i]:=1;
    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:=power(r,alpha/2);
      r:=r/s;
      if r>20 then
        r:=0
      else
        r:=exp(-r);
      d[i,j]:=r;
      d[j,i]:=r;
    end;
  end;
end;




procedure UnbiasedDoubleCenter(n:integer; var d:distmatrix);
var i,j:integer; t,a,r0,r1: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; // first pass, biased double centering

  a:=0;
  for i:=1 to n do
    a:=a+d[i,i]; // trace
  a:=a/(n-1);
  r0:=a*(1.0-1.0/n);
  r1:=a/n;
  for i:=1 to n do begin
    d[i,i]:=d[i,i]-r0;
    for j:=i+1 to n do begin
      d[i,j]:=d[i,j]+r1;
      d[j,i]:=d[i,j];
    end;
  end; // final unbiased 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 SimCorr(n,p,q:integer; alpha,sx,sy:real;
                 var x,y:datamatrix;
                 var xy,xx,yy:real
                 ):real;
{n:number of samples; x:n*p; y:n*q; alpha,sx,sy>0:parameters}
begin
  Data2Sim(n,p,alpha,sx,x,dx);
  Data2Sim(n,q,alpha,sy,y,dy); // similarity matrix

  UnbiasedDoubleCenter(n,dx);
  UnbiasedDoubleCenter(n,dy);

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

  if (xx<eps) or (yy<eps) then
    SimCorr:=0
  else
    SimCorr:=xy/sqrt(xx*yy);
end;




procedure CmplxData2Sim(n,p:integer; alpha,s:real;
                        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]:=1;
    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:=power(r,alpha/2);
      r:=r/s;
      if r>20 then
        r:=0
      else
        r:=exp(-r);
      d[i,j]:=r;
      d[j,i]:=r;
    end;
  end;
end;




procedure CmplxData2Dist(n,p:integer; alpha,s:real;
                         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]:=1;
    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:=power(r,alpha/2);
      r:=r/s;
      if r>20 then
        r:=0
      else
        r:=exp(-r);
      d[i,j]:=r;
      d[j,i]:=r;
    end;
  end;
end;




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

  Data2Sim(n,p,alpha,sx,xredata,xredist);
  Data2Sim(n,p,alpha,sx,ximdata,ximdist);
  Data2Sim(n,q,alpha,sy,yredata,yredist);
  Data2Sim(n,q,alpha,sy,yimdata,yimdist);

  UnbiasedDoubleCenter(n,xredist);
  UnbiasedDoubleCenter(n,ximdist);
  UnbiasedDoubleCenter(n,yredist);
  UnbiasedDoubleCenter(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.

