
{$R-}

unit uLaplacianCorrMath;

interface

uses uMatToolsDynSingle33;

var dCorrel,REdCorrel,IMdCorrel,dCovarXY,dVarX,dVarY:real;
// read only! contains dCorr and dCov after call to dCorr0 for reals (no coherence stuff)

var DistExp:real=2.0; // write only. this is the distance exponent, default=2
    eps3:extended=1e-5;

procedure Data2Dist(var x,d:matrix; s2:real); overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; creates d if d.nr=0

procedure Data2Dist(var x:vector; var d:matrix; s2:real); overload;
// in data x:n*1; n objects; p vars; out d:n*n distance; creates d if d.nr=0

procedure CenterDist(var d:matrix); overload;
// makes rows and cols have zero mean

function dCov(var CenterDist1,CenterDist2:matrix):real; overload;
// computes dist covar from centered distance matrices

function dCorr(var CenterDist1,CenterDist2:matrix):real; overload;
// computes dist corr from centered distance matrices

// computes dist corr from data matrices or vectors:
function dCorr0(var datax,datay:matrix; xs2,ys2:real):real; overload;
function dCorr0(var datax,datay:vector; xs2,ys2:real):real; overload;

// computes DISTANCE COHERENCE from data matrices or vectors:
procedure dCorr0a(var xReData,xImData,yReData,yImData:matrix;
                  var TotCoh,ReCoh,ImCoh:real;
                  xs2,ys2:real); overload;
procedure dCorr0a(var xReData,xImData,yReData,yImData:matrix;
                  var TotCoh,ReCoh,ImCoh:real
                  ); overload;

function dCorr0TimeSeries(var datax,datay:matrix; xs2,ys2:real; TimeWidth:integer):real; overload;

function dCovar0(var datax,datay:vector; xs2,ys2:real):real; overload;
function dCovar0TimeSeries(var datax,datay:matrix; xs2,ys2:real; TimeWidth:integer):real; overload;
function dCorr0TimeSeries(var datax,datay:matrix; TimeWidth:integer):real; overload;
function dCorr0(var datax,datay:vector):real; overload;

implementation

uses math;




procedure CenterDist(var d:matrix); overload;
// makes rows and cols have zero mean
var i,j,n:integer; a,r0,r1:real;
begin
  n:=d.nr;
  subtractrowmeans(d);
  subtractcolumnmeans(d);
  a:=0;
  for i:=1 to n do
    a:=a+d.m[i,i];
  a:=a/(n-1);
  r0:=a*(1.0-1.0/n);
  r1:=a/n;
  for i:=1 to n do begin
    d.m[i,i]:=d.m[i,i]-r0;
    for j:=i+1 to n do begin
      d.m[i,j]:=d.m[i,j]+r1;
      d.m[j,i]:=d.m[i,j];
    end;
  end;
end;




function dCov(var CenterDist1,CenterDist2:matrix):real; overload;
// computes dist covar from centered distance matrices
var r:extended; n,i,j:integer;
begin
  r:=0;
  n:=CenterDist1.nr;
  for i:=1 to n do
    for j:=i+1 to n do
      r:=r+CenterDist1.m[i,j]*CenterDist2.m[i,j];
  r:=r*2;
  for i:=1 to n do
    r:=r+CenterDist1.m[i,i]*CenterDist2.m[i,i];
  dCov:=r/sqr(n);
end;




function dCorr(var CenterDist1,CenterDist2:matrix):real; overload;
// computes dist corr from centered distance matrices
var xy,xx,yy:extended;
begin
  xy:=dCov(CenterDist1,CenterDist2);
  dCovarXY:=xy;
  xx:=dCov(CenterDist1,CenterDist1);
  yy:=dCov(CenterDist2,CenterDist2);
  if xx<eps3 then xx:=0;
  if yy<eps3 then yy:=0;
  if sqrt(xx*yy)<eps3 then begin
    xx:=0;
    yy:=0;
  end;
  if (xx<eps3) or (yy<eps3) then
    dCorr:=0
  else
    dCorr:=xy/sqrt(xx*yy);
  dvarx:=xx;
  dvary:=yy;
end;




function dCorr0(var datax,datay:matrix; xs2,ys2:real):real; overload;
var distx,disty:matrix;
begin
  distx.nr:=0; Data2Dist(datax,distx,xs2);
  disty.nr:=0; Data2Dist(datay,disty,ys2);
  CenterDist(distx);
  CenterDist(disty);
  dCorrel:=dCorr(distx,disty);
  dCorr0:=dCorrel;
  destroy1(distx);
  destroy1(disty);
end;




function dCorr0(var datax,datay:vector; xs2,ys2:real):real; overload;
var distx,disty:matrix;
begin
  distx.nr:=0; Data2Dist(datax,distx,xs2);
  disty.nr:=0; Data2Dist(datay,disty,ys2);
  CenterDist(distx);
  CenterDist(disty);
  dCorrel:=dCorr(distx,disty);
  dCorr0:=dCorrel;
  destroy1(distx);
  destroy1(disty);
end;




function dCovar0(var datax,datay:vector; xs2,ys2:real):real; overload;
var distx,disty:matrix;
begin
  distx.nr:=0; Data2Dist(datax,distx,xs2);
  disty.nr:=0; Data2Dist(datay,disty,ys2);
  CenterDist(distx);
  CenterDist(disty);
  dCovar0:=dCov(distx,disty);
  destroy1(distx);
  destroy1(disty);
end;




procedure Data2Dist(var x,d:matrix; s2:real); overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,p,i,j,k:integer; r:extended;
begin
  n:=x.nr;
  p:=x.nc;
  if d.nr=0 then
    create1(n,n,d);
  for i:=1 to n do begin
    d.m[i,i]:=1;
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(x.m[i,k]-x.m[j,k]);
      r:=power(r,distexp/2);
      r:=r/s2;
      if r>20 then
        r:=0
      else
        r:=exp(-r);
      d.m[i,j]:=r;
      d.m[j,i]:=r;
    end;
  end;
end;




function Data2AveDist(var x:matrix):real; overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,p,i,j,k:integer; r,a:extended;
begin
  n:=x.nr;
  p:=x.nc;
  a:=0;
  for i:=1 to n do
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(x.m[i,k]-x.m[j,k]);
      r:=power(r,distexp/2);
      a:=a+r;
    end;
  result:=2*a/(n*(n-1));
end;




procedure Data2Dist(var x:vector; var d:matrix; s2:real); overload;
// in data x:n*1; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,i,j:integer; r:extended;
begin
  n:=x.n;
  if d.nr=0 then
    create1(n,n,d);
  for i:=1 to n do begin
    d.m[i,i]:=1;
    for j:=i+1 to n do begin
      r:=sqr(x.v[i]-x.v[j]);
      r:=power(r,distexp/2);
      r:=r/s2;
      if r>20 then
        r:=0
      else
        r:=exp(-r);
      d.m[i,j]:=r;
      d.m[j,i]:=r;
    end;
  end;
end;




function Data2AveDist(var x:vector):real; overload;
// in data x:n*1; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,i,j:integer; r,a:extended;
begin
  n:=x.n;
  a:=0;
  for i:=1 to n do
    for j:=i+1 to n do begin
      r:=sqr(x.v[i]-x.v[j]);
      r:=power(r,distexp/2);
      a:=a+r
    end;
  result:=2*a/(n*(n-1));
end;




procedure CmplxData2Dist(var xre,xim,d:matrix; s2:real); overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,p,i,j,k:integer; r:extended;
begin
  n:=xre.nr;
  p:=xre.nc;
  if d.nr=0 then
    create1(n,n,d);
  for i:=1 to n do begin
    d.m[i,i]:=1;
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(xre.m[i,k]-xre.m[j,k])+sqr(xim.m[i,k]-xim.m[j,k]);
      r:=power(r,distexp/2);
      r:=r/s2;
      if r>20 then
        r:=0
      else
        r:=exp(-r);
      d.m[i,j]:=r;
      d.m[j,i]:=r;
    end;
  end;
end;




function CmplxData2AveDist(var xre,xim:matrix):real; overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,p,i,j,k:integer; r,a:extended;
begin
  n:=xre.nr;
  p:=xre.nc;
  a:=0;
  for i:=1 to n do
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(xre.m[i,k]-xre.m[j,k])+sqr(xim.m[i,k]-xim.m[j,k]);
      r:=power(r,distexp/2);
      a:=a+r;
    end;
  result:=2*a/(n*(n-1));
end;




procedure dCorr0a(var xReData,xImData,yReData,yImData:matrix;
                  var TotCoh,ReCoh,ImCoh:real;
                  xs2,ys2:real); overload;
var xReDist,xImDist,yReDist,yImDist:matrix;
    zxryr,zxryi,zxiyr,zxiyi,zden:extended;
    totcovar,totvarx,totvary,fre,fim,ftot:extended;
begin
  xReDist.nr:=0; CmplxData2Dist(xReData,xImData,xReDist,xs2);
  yReDist.nr:=0; CmplxData2Dist(yReData,yImData,yReDist,ys2);
  CenterDist(xReDist);
  CenterDist(yReDist);
  totcovar:=dCov(xReDist,yReDist);
  totvarx:=dCov(xReDist,xReDist);
  totvary:=dCov(yReDist,yReDist);
  destroy1(xReDist);
  destroy1(yReDist);

  xReDist.nr:=0; Data2Dist(xReData,xReDist,xs2);
  xImDist.nr:=0; Data2Dist(xImData,xImDist,xs2);
  yReDist.nr:=0; Data2Dist(yReData,yReDist,ys2);
  yImDist.nr:=0; Data2Dist(yImData,yImDist,ys2);

  CenterDist(xReDist);
  CenterDist(xImDist);
  CenterDist(yReDist);
  CenterDist(yImDist);

  zxryr:=dCov(xReDist,yReDist);
  zxryi:=dCov(xReDist,yImDist);
  zxiyr:=dCov(xImDist,yReDist);
  zxiyi:=dCov(xImDist,yImDist);

  if totvarx<eps3 then totvarx:=0;
  if totvary<eps3 then totvary:=0;
  if sqrt(totvary*totvarx)<eps3 then begin
    totvarx:=0;
    totvary:=0;
  end;

  if (totvarx<eps3) or (totvary<eps3) then begin
    totcoh:=0;
    recoh:=0;
    imcoh:=0;
    dCorrel:=0;
    REdCorrel:=0;
    IMdCorrel:=0;
    dCovarXY:=totcovar;
    dVarX:=totvarx;
    dVarY:=totvary;
    destroy1(xRedist);
    destroy1(xImdist);
    destroy1(yRedist);
    destroy1(yImdist);
    exit;
  end;

  if totvarx<eps3 then totvarx:=0;
  if totvary<eps3 then totvary:=0;
  if sqrt(totvary*totvarx)<eps3 then begin
    totvarx:=0;
    totvary:=0;
  end;

  if (totvarx<eps3) or (totvary<eps3) then begin
    totcoh:=0;
    recoh:=0;
    imcoh:=0;
    dCorrel:=0;
    REdCorrel:=0;
    IMdCorrel:=0;
    dCovarXY:=totcovar;
    dVarX:=totvarx;
    dVarY:=totvary;
    destroy1(xRedist);
    destroy1(xImdist);
    destroy1(yRedist);
    destroy1(yImdist);
    exit;
  end;

  zden:=sqrt(totvarx*totvary);
  fre:=zxryr+zxiyi;
  fim:=zxryi+zxiyr;
  ftot:=fre+fim;
  if ftot<1e-7 then ftot:=1;
  fre:=fre/ftot;
  fim:=fim/ftot;

  totcoh:=totcovar/zden;
  recoh:=fre*totcoh;
  imcoh:=fim*totcoh;

  dCorrel:=totcoh;
  REdCorrel:=recoh;
  IMdCorrel:=imcoh;
  dCovarXY:=totcovar;
  dVarX:=totvarx;
  dVarY:=totvary;

  destroy1(xRedist);
  destroy1(xImdist);
  destroy1(yRedist);
  destroy1(yImdist);

end;




procedure CmplxData2Dist(var xre,xim,d:matrix); overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,p,i,j,k:integer; r:extended;
begin
  n:=xre.nr;
  p:=xre.nc;
  if d.nr=0 then
    create1(n,n,d);
  for i:=1 to n do begin
    d.m[i,i]:=0;
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(xre.m[i,k]-xre.m[j,k])+sqr(xim.m[i,k]-xim.m[j,k]);
      r:=sqrt(r);
      d.m[i,j]:=r;
      d.m[j,i]:=r;
    end;
  end;
end;




procedure Data2Dist(var x,d:matrix); overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,p,i,j,k:integer; r:extended;
begin
  n:=x.nr;
  p:=x.nc;
  if d.nr=0 then
    create1(n,n,d);
  for i:=1 to n do begin
    d.m[i,i]:=0;
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(x.m[i,k]-x.m[j,k]);
      r:=sqrt(r);
      d.m[i,j]:=r;
      d.m[j,i]:=r;
    end;
  end;
end;




procedure dCorr0a(var xReData,xImData,yReData,yImData:matrix;
                  var TotCoh,ReCoh,ImCoh:real
                  ); overload;
var xReDist,xImDist,yReDist,yImDist:matrix;
    zxryr,zxryi,zxiyr,zxiyi,zden:extended;
    totcovar,totvarx,totvary,fre,fim,ftot:extended;
begin
  xReDist.nr:=0; CmplxData2Dist(xReData,xImData,xReDist);
  yReDist.nr:=0; CmplxData2Dist(yReData,yImData,yReDist);
  CenterDist(xReDist);
  CenterDist(yReDist);
  totcovar:=dCov(xReDist,yReDist);
  totvarx:=dCov(xReDist,xReDist);
  totvary:=dCov(yReDist,yReDist);
  destroy1(xReDist);
  destroy1(yReDist);

  xReDist.nr:=0; Data2Dist(xReData,xReDist);
  xImDist.nr:=0; Data2Dist(xImData,xImDist);
  yReDist.nr:=0; Data2Dist(yReData,yReDist);
  yImDist.nr:=0; Data2Dist(yImData,yImDist);

  CenterDist(xReDist);
  CenterDist(xImDist);
  CenterDist(yReDist);
  CenterDist(yImDist);

  zxryr:=dCov(xReDist,yReDist);
  zxryi:=dCov(xReDist,yImDist);
  zxiyr:=dCov(xImDist,yReDist);
  zxiyi:=dCov(xImDist,yImDist);

  if totvarx<eps3 then totvarx:=0;
  if totvary<eps3 then totvary:=0;
  if sqrt(totvary*totvarx)<eps3 then begin
    totvarx:=0;
    totvary:=0;
  end;

  if (totvarx<eps3) or (totvary<eps3) then begin
    totcoh:=0;
    recoh:=0;
    imcoh:=0;
    dCorrel:=0;
    REdCorrel:=0;
    IMdCorrel:=0;
    dCovarXY:=totcovar;
    dVarX:=totvarx;
    dVarY:=totvary;
    destroy1(xRedist);
    destroy1(xImdist);
    destroy1(yRedist);
    destroy1(yImdist);
    exit;
  end;

  zden:=sqrt(totvarx*totvary);
  fre:=zxryr+zxiyi;
  fim:=zxryi+zxiyr;
  ftot:=fre+fim;
  if ftot<1e-7 then ftot:=1;
  fre:=fre/ftot;
  fim:=fim/ftot;

  totcoh:=totcovar/zden;
  recoh:=fre*totcoh;
  imcoh:=fim*totcoh;

  dCorrel:=totcoh;
  REdCorrel:=recoh;
  IMdCorrel:=imcoh;
  dCovarXY:=totcovar;
  dVarX:=totvarx;
  dVarY:=totvary;

  destroy1(xRedist);
  destroy1(xImdist);
  destroy1(yRedist);
  destroy1(yImdist);
end;




procedure Data2DistTimeSeries(var x,d:matrix; s2:real; ii:integer); overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; must pre-create d
var n,p,i,j,k,kk:integer; r:extended;
begin
  n:=d.nr;
  p:=x.nc;
  kk:=(ii-1)*n;
  for i:=1 to n do begin
    d.m[i,i]:=1;
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(x.m[kk+i,k]-x.m[kk+j,k]);
      r:=power(r,distexp/2);
      r:=r/s2;
      if r>20 then
        r:=0
      else
        r:=exp(-r);
      d.m[i,j]:=r;
      d.m[j,i]:=r;
    end;
  end;
end;




function Data2AveDistTimeSeries(var x:matrix; WindowSize:integer):real; overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var ii,nw,p,i,j,k,kk:integer; r,a:extended;
begin
  p:=x.nc;
  nw:=x.nr div WindowSize;

  a:=0;

  for ii:=1 to nw do begin
    kk:=(ii-1)*WindowSize;
    for i:=1 to WindowSize do
      for j:=i+1 to WindowSize do begin
        r:=0;
        for k:=1 to p do
          r:=r+sqr(x.m[kk+i,k]-x.m[kk+j,k]);
        r:=power(r,distexp/2);
        a:=a+r;
      end;
  end;

  a:=2*a/(WindowSize*(WindowSize-1));
  result:=a/nw;
end;




function dCorr0TimeSeries(var datax,datay:matrix; xs2,ys2:real; TimeWidth:integer):real; overload;
var distx,disty:matrix; n,nr,i:integer; xy,xx,yy:extended;
begin
  create1(timewidth,timewidth,distx);
  create1(timewidth,timewidth,disty);
  n:=datax.nr;
  nr:=n div timewidth;
  xy:=0;
  xx:=0;
  yy:=0;
  for i:=1 to nr do begin
    Data2DistTimeSeries(datax,distx,xs2,i);
    Data2DistTimeSeries(datay,disty,ys2,i);
    CenterDist(distx);
    CenterDist(disty);
    xy:=xy+dCov(distx,disty);
    xx:=xx+dCov(distx,distx);
    yy:=yy+dCov(disty,disty);
  end;
  xy:=xy/nr;
  xx:=xx/nr;
  yy:=yy/nr;

  dCovarXY:=xy;
  if xx<eps3 then xx:=0;
  if yy<eps3 then yy:=0;
  if sqrt(xx*yy)<eps3 then begin
    xx:=0;
    yy:=0;
  end;
  if (xx<eps3) or (yy<eps3) then
    dCorrel:=0
  else
    dCorrel:=xy/sqrt(xx*yy);
  dCorr0TimeSeries:=dcorrel;
  dvarx:=xx;
  dvary:=yy;
  destroy1(distx);
  destroy1(disty);
end;




function dCovar0TimeSeries(var datax,datay:matrix; xs2,ys2:real; TimeWidth:integer):real; overload;
var distx,disty:matrix; n,nr,i:integer; xy:extended;
begin
  create1(timewidth,timewidth,distx);
  create1(timewidth,timewidth,disty);
  n:=datax.nr;
  nr:=n div timewidth;
  xy:=0;
  for i:=1 to nr do begin
    Data2DistTimeSeries(datax,distx,xs2,i);
    Data2DistTimeSeries(datay,disty,ys2,i);
    CenterDist(distx);
    CenterDist(disty);
    xy:=xy+dCov(distx,disty);
  end;
  xy:=xy/nr;
  dCovar0TimeSeries:=xy;
  destroy1(distx);
  destroy1(disty);
end;




procedure Data2DistTimeSeries(var x,d:matrix; ii:integer); overload;
// in data x:n*p; n objects; p vars; out d:n*n distance; must pre-create d
var n,p,i,j,k,kk:integer; r:extended;
begin
  n:=d.nr;
  p:=x.nc;
  kk:=(ii-1)*n;
  for i:=1 to n do begin
    d.m[i,i]:=0;
    for j:=i+1 to n do begin
      r:=0;
      for k:=1 to p do
        r:=r+sqr(x.m[kk+i,k]-x.m[kk+j,k]);
      r:=sqrt(r);
      d.m[i,j]:=r;
      d.m[j,i]:=r;
    end;
  end;
end;




function dCorr0TimeSeries(var datax,datay:matrix; TimeWidth:integer):real; overload;
var distx,disty:matrix; n,nr,i:integer; xy,xx,yy:extended;
begin
  create1(timewidth,timewidth,distx);
  create1(timewidth,timewidth,disty);
  n:=datax.nr;
  nr:=n div timewidth;
  xy:=0;
  xx:=0;
  yy:=0;
  for i:=1 to nr do begin
    Data2DistTimeSeries(datax,distx,i);
    Data2DistTimeSeries(datay,disty,i);
    CenterDist(distx);
    CenterDist(disty);
    xy:=xy+dCov(distx,disty);
    xx:=xx+dCov(distx,distx);
    yy:=yy+dCov(disty,disty);
  end;
  xy:=xy/nr;
  xx:=xx/nr;
  yy:=yy/nr;

  dCovarXY:=xy;
  if xx<eps3 then xx:=0;
  if yy<eps3 then yy:=0;
  if sqrt(xx*yy)<eps3 then begin
    xx:=0;
    yy:=0;
  end;
  if (xx<eps3) or (yy<eps3) then
    dCorrel:=0
  else
    dCorrel:=xy/sqrt(xx*yy);
  dCorr0TimeSeries:=dcorrel;
  dvarx:=xx;
  dvary:=yy;
  destroy1(distx);
  destroy1(disty);
end;




procedure Data2Dist(var x:vector; var d:matrix); overload;
// in data x:n*1; n objects; p vars; out d:n*n distance; creates d if d.nr=0
var n,i,j:integer; r:extended;
begin
  n:=x.n;
  if d.nr=0 then
    create1(n,n,d);
  for i:=1 to n do begin
    d.m[i,i]:=0;
    for j:=i+1 to n do begin
      r:=abs(x.v[i]-x.v[j]);
      d.m[i,j]:=r;
      d.m[j,i]:=r;
    end;
  end;
end;




function dCorr0(var datax,datay:vector):real; overload;
var distx,disty:matrix;
begin
  distx.nr:=0; Data2Dist(datax,distx);
  disty.nr:=0; Data2Dist(datay,disty);
  CenterDist(distx);
  CenterDist(disty);
  dCorrel:=dCorr(distx,disty);
  dCorr0:=dCorrel;
  destroy1(distx);
  destroy1(disty);
end;




end.

