
unit uSimDistCorr;

{$MODE Delphi}

interface

uses
  LCLIntf, LCLType, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, uFunctGeneration, uVecMat, StdCtrls, ExtCtrls, uDistCorrMath,
  uSimCorrMath, uHookeJeeves;

type
  TfSimCorr = class(TForm)
    LineButt: TButton;
    distexpTxt: TLabeledEdit;
    XweightTxt: TLabeledEdit;
    YweightTxt: TLabeledEdit;
    noisetxt: TLabeledEdit;
    dcorrTxt: TLabeledEdit;
    scorrTxt: TLabeledEdit;
    MaxSxyChk: TCheckBox;
    Xbutt: TButton;
    SineLineButt: TButton;
    RandBoxButt: TButton;
    CircleButt: TButton;
    ParaboloidButt: TButton;
    CmplxLinButt: TButton;
    GroupBox1: TGroupBox;
    ReTxt: TLabeledEdit;
    ImTxt: TLabeledEdit;
    GroupBox2: TGroupBox;
    totsimtxt: TLabeledEdit;
    resimtxt: TLabeledEdit;
    imsimtxt: TLabeledEdit;
    GroupBox3: TGroupBox;
    totdisttxt: TLabeledEdit;
    redisttxt: TLabeledEdit;
    imdisttxt: TLabeledEdit;
    CmplxQuadButt: TButton;
    CmplxSqrtButt: TButton;
    CmplxLogButt: TButton;
    Panel1: TPanel;
    procedure LineButtClick(Sender: TObject);
    procedure XbuttClick(Sender: TObject);
    procedure SineLineButtClick(Sender: TObject);
    procedure RandBoxButtClick(Sender: TObject);
    procedure CircleButtClick(Sender: TObject);
    procedure ParaboloidButtClick(Sender: TObject);
    procedure CmplxLinButtClick(Sender: TObject);
    procedure CmplxQuadButtClick(Sender: TObject);
    procedure CmplxSqrtButtClick(Sender: TObject);
    procedure CmplxLogButtClick(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
    n:integer;
    alpha,xx,yy,xy,simTotCoh,simReCoh,simImCoh:real;
    xre,xim,yre,yim,x,y:datamatrix;
  end;

var
  fSimCorr: TfSimCorr;

implementation

{$R *.lfm}

procedure AddNoise1(n:integer; a:real; var y:datamatrix);
var i:integer;
begin
  randomize;
  for i:=1 to n do
    y[i,1]:=y[i,1]+2*(random-0.5)*a;
end;

procedure WriteTxt1(fn:string; n:integer; var x,y:datamatrix);
var txt:textfile; i:integer;
begin
  assignfile(txt,fn);
  filemode:=1;
  rewrite(txt);
  for i:=1 to n do
    writeln(txt,' ',x[i,1]:15,' ',y[i,1]:15);
  closeFile(txt);
end;

function fHJSimCov(var a:vvector):real;
begin
  a[1]:=exp(a[1]); a[2]:=exp(a[2]);
  SimCorr(fSimCorr.n,1,1,fSimCorr.alpha,a[1],a[2],
          fSimCorr.x,fSimCorr.y,fSimCorr.xy,fSimCorr.xx,fSimCorr.yy);
  fHJSimCov:=-fSimCorr.xy;
  a[1]:=ln(a[1]); a[2]:=ln(a[2]);
end;

function fHJSimCorr(var a:vvector):real;
begin
  a[1]:=exp(a[1]); a[2]:=exp(a[2]);
  fHJSimCorr:=-SimCorr(fSimCorr.n,1,1,fSimCorr.alpha,a[1],a[2],
          fSimCorr.x,fSimCorr.y,fSimCorr.xy,fSimCorr.xx,fSimCorr.yy);
  a[1]:=ln(a[1]); a[2]:=ln(a[2]);
end;

procedure TfSimCorr.LineButtClick(Sender: TObject);
var eps,scorr,a,dcorr,sx,sy,r:real; fn:string;
    p1,d1:vvector;
begin
  fn:='Line.txt';
  n:=1000;
  linear(n,x,y);
  a:=strtofloat(noisetxt.Text);
  AddNoise1(n,a,y);
  WriteTxt1(fn,n,x,y);
  dcorr:=DistCorr(n,1,1,x,y);
  dcorrtxt.Text:=floattostr(dcorr);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    scorr:=SimCorr(n,1,1,alpha,sx,sy,x,y,xy,xx,yy);
    scorrtxt.Text:=floattostr(scorr);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJSimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJSimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    scorr:=-r;
    scorrtxt.Text:=floattostr(scorr);
  end;

  showmessage('finished');
end;

procedure TfSimCorr.XbuttClick(Sender: TObject);
var eps,scorr,a,dcorr,sx,sy,r:real; fn:string;
    p1,d1:vvector;
begin
  fn:='Xshape.txt';
  n:=2000;
  xshape(n,x,y);
  a:=strtofloat(noisetxt.Text);
  AddNoise1(n,a,y);
  WriteTxt1(fn,n,x,y);
  dcorr:=DistCorr(n,1,1,x,y);
  dcorrtxt.Text:=floattostr(dcorr);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    scorr:=SimCorr(n,1,1,alpha,sx,sy,x,y,xy,xx,yy);
    scorrtxt.Text:=floattostr(scorr);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJSimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJSimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    scorr:=-r;
    scorrtxt.Text:=floattostr(scorr);
  end;

  showmessage('finished');
end;

procedure TfSimCorr.SineLineButtClick(Sender: TObject);
var eps,scorr,a,dcorr,sx,sy,r:real; fn:string;
    p1,d1:vvector;
begin
  fn:='SinePlusLine.txt';
  n:=1000;
  LinPeriod1(n,x,y);
  a:=strtofloat(noisetxt.Text);
  AddNoise1(n,a,y);
  WriteTxt1(fn,n,x,y);
  dcorr:=DistCorr(n,1,1,x,y);
  dcorrtxt.Text:=floattostr(dcorr);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    scorr:=SimCorr(n,1,1,alpha,sx,sy,x,y,xy,xx,yy);
    scorrtxt.Text:=floattostr(scorr);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJSimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJSimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    scorr:=-r;
    scorrtxt.Text:=floattostr(scorr);
  end;

  showmessage('finished');
end;

procedure TfSimCorr.RandBoxButtClick(Sender: TObject);
var eps,scorr,a,dcorr,sx,sy,r:real; fn:string;
    p1,d1:vvector;
begin
  fn:='RandomDotsBox.txt';
  n:=1000;
  RandBox(n,x,y);
  a:=strtofloat(noisetxt.Text);
  AddNoise1(n,a,y);
  WriteTxt1(fn,n,x,y);
  dcorr:=DistCorr(n,1,1,x,y);
  dcorrtxt.Text:=floattostr(dcorr);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    scorr:=SimCorr(n,1,1,alpha,sx,sy,x,y,xy,xx,yy);
    scorrtxt.Text:=floattostr(scorr);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJSimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJSimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    scorr:=-r;
    scorrtxt.Text:=floattostr(scorr);
  end;

  showmessage('finished');
end;

procedure TfSimCorr.CircleButtClick(Sender: TObject);
var eps,scorr,a,dcorr,sx,sy,r:real; fn:string;
    p1,d1:vvector;
begin
  fn:='Circle.txt';
  n:=1000;
  CircleFull(n,x,y);
  a:=strtofloat(noisetxt.Text);
  AddNoise1(n,a,y);
  WriteTxt1(fn,n,x,y);
  dcorr:=DistCorr(n,1,1,x,y);
  dcorrtxt.Text:=floattostr(dcorr);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    scorr:=SimCorr(n,1,1,alpha,sx,sy,x,y,xy,xx,yy);
    scorrtxt.Text:=floattostr(scorr);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJSimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJSimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    scorr:=-r;
    scorrtxt.Text:=floattostr(scorr);
  end;

  showmessage('finished');
end;

procedure WriteTxt2(fn:string; n:integer; var x,y:datamatrix);
var txt:textfile; i:integer;
begin
  assignfile(txt,fn);
  filemode:=1;
  rewrite(txt);
  for i:=1 to n do
    writeln(txt,' ',x[i,1]:15,' ',y[i,1]:15,' ',y[i,2]:15);
  closeFile(txt);
end;

function fHJ2SimCov(var a:vvector):real;
begin
  a[1]:=exp(a[1]); a[2]:=exp(a[2]);
  SimCorr(fSimCorr.n,1,2,fSimCorr.alpha,a[1],a[2],
          fSimCorr.x,fSimCorr.y,fSimCorr.xy,fSimCorr.xx,fSimCorr.yy);
  fHJ2SimCov:=-fSimCorr.xy;
  a[1]:=ln(a[1]); a[2]:=ln(a[2]);
end;

function fHJ2SimCorr(var a:vvector):real;
begin
  a[1]:=exp(a[1]); a[2]:=exp(a[2]);
  fHJ2SimCorr:=-SimCorr(fSimCorr.n,1,2,fSimCorr.alpha,a[1],a[2],
          fSimCorr.x,fSimCorr.y,fSimCorr.xy,fSimCorr.xx,fSimCorr.yy);
  a[1]:=ln(a[1]); a[2]:=ln(a[2]);
end;

procedure TfSimCorr.ParaboloidButtClick(Sender: TObject);
var eps,scorr,a,dcorr,sx,sy,r:real; fn:string;
    p1,d1:vvector;
begin
  fn:='Paraboloid.txt';
  n:=900;
  paraboloid(n,x,y);
  a:=strtofloat(noisetxt.Text);
  AddNoise1(n,a,y);
  WriteTxt2(fn,n,x,y);
  dcorr:=DistCorr(n,1,2,x,y);
  dcorrtxt.Text:=floattostr(dcorr);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    scorr:=SimCorr(n,1,2,alpha,sx,sy,x,y,xy,xx,yy);
    scorrtxt.Text:=floattostr(scorr);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJ2SimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJ2SimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    scorr:=-r;
    scorrtxt.Text:=floattostr(scorr);
  end;

  showmessage('finished');
end;

procedure WriteTxt3(fn:string; n:integer; var xre,xim,yre,yim:datamatrix);
var txt:textfile; i:integer;
begin
  assignfile(txt,fn);
  filemode:=1;
  rewrite(txt);
  for i:=1 to n do
    writeln(txt,' ',xre[i,1]:15,' ',xim[i,1]:15,' ',yre[i,1]:15,' ',yim[i,1]:15);
  closeFile(txt);
end;

procedure AddNoise2(n:integer; a:real; var y:datamatrix);
var i:integer;
begin
  randomize;
  for i:=1 to n do begin
    y[i,1]:=y[i,1]+2*(random-0.5)*a;
    y[i,2]:=y[i,2]+2*(random-0.5)*a;
  end;
end;

function fHJ3SimCov(var a:vvector):real;
begin
  a[1]:=exp(a[1]); a[2]:=exp(a[2]);
  SimCoh(fSimCorr.n,1,1,fSimCorr.alpha,a[1],a[2],
    fSimCorr.xre,fSimCorr.xim,fSimCorr.yre,fSimCorr.yim,
    fSimCorr.simTotCoh,fSimCorr.simReCoh,fSimCorr.simImCoh,fSimCorr.xy);
  fHJ3SimCov:=-fSimCorr.xy;
  a[1]:=ln(a[1]); a[2]:=ln(a[2]);
end;

function fHJ3SimCorr(var a:vvector):real;
begin
  a[1]:=exp(a[1]); a[2]:=exp(a[2]);
  SimCoh(fSimCorr.n,1,1,fSimCorr.alpha,a[1],a[2],
    fSimCorr.xre,fSimCorr.xim,fSimCorr.yre,fSimCorr.yim,
    fSimCorr.simTotCoh,fSimCorr.simReCoh,fSimCorr.simImCoh,fSimCorr.xy);
  fHJ3SimCorr:=-fSimCorr.simTotCoh;
  a[1]:=ln(a[1]); a[2]:=ln(a[2]);
end;

procedure TfSimCorr.CmplxLinButtClick(Sender: TObject);
var eps,a,sx,sy,r,are,aim,dTotCoh,dReCoh,dImCoh:real; fn:string;
    p1,d1:vvector;
begin
  fn:='CmplxLinear.txt';
  n:=900;
  are:=strtofloat(retxt.Text);
  aim:=strtofloat(imtxt.Text);
  CmplxLinear(n,are,aim,xre,xim,yre,yim);
  a:=strtofloat(noisetxt.Text);
  AddNoise2(n,a,y);
  WriteTxt3(fn,n,xre,xim,yre,yim);
  DistCoh(n,1,1,xre,xim,yre,yim,dTotCoh,dReCoh,dImCoh);
  totdisttxt.Text:=floattostr(dTotCoh);
  redisttxt.Text:=floattostr(dreCoh);
  imdisttxt.Text:=floattostr(dimCoh);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    SimCoh(n,1,1,alpha,sx,sy,xre,xim,yre,yim,simTotCoh,simReCoh,simImCoh,xy);
    totsimtxt.Text:=floattostr(simTotCoh);
    resimtxt.Text:=floattostr(simreCoh);
    imsimtxt.Text:=floattostr(simimCoh);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJ3SimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJ3SimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    totsimtxt.Text:=floattostr(-r);
    resimtxt.Text:=floattostr(simreCoh);
    imsimtxt.Text:=floattostr(simimCoh);
  end;

  showmessage('finished');
end;

procedure TfSimCorr.CmplxQuadButtClick(Sender: TObject);
var eps,a,sx,sy,r,are,aim,dTotCoh,dReCoh,dImCoh:real; fn:string;
    p1,d1:vvector;
begin
  fn:='CmplxQuad.txt';
  n:=900;
  are:=strtofloat(retxt.Text);
  aim:=strtofloat(imtxt.Text);
  CmplxQuad(n,are,aim,xre,xim,yre,yim);
  a:=strtofloat(noisetxt.Text);
  AddNoise2(n,a,y);
  WriteTxt3(fn,n,xre,xim,yre,yim);
  DistCoh(n,1,1,xre,xim,yre,yim,dTotCoh,dReCoh,dImCoh);
  totdisttxt.Text:=floattostr(dTotCoh);
  redisttxt.Text:=floattostr(dreCoh);
  imdisttxt.Text:=floattostr(dimCoh);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    SimCoh(n,1,1,alpha,sx,sy,xre,xim,yre,yim,simTotCoh,simReCoh,simImCoh,xy);
    totsimtxt.Text:=floattostr(simTotCoh);
    resimtxt.Text:=floattostr(simreCoh);
    imsimtxt.Text:=floattostr(simimCoh);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJ3SimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJ3SimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    totsimtxt.Text:=floattostr(-r);
    resimtxt.Text:=floattostr(simreCoh);
    imsimtxt.Text:=floattostr(simimCoh);
  end;

  showmessage('finished');
end;

procedure TfSimCorr.CmplxSqrtButtClick(Sender: TObject);
var eps,a,sx,sy,r,are,aim,dTotCoh,dReCoh,dImCoh:real; fn:string;
    p1,d1:vvector;
begin
  fn:='CmplxSqrt.txt';
  n:=900;
  are:=strtofloat(retxt.Text);
  aim:=strtofloat(imtxt.Text);
  CmplxPrincipleSqrt(n,are,aim,xre,xim,yre,yim);
  a:=strtofloat(noisetxt.Text);
  AddNoise2(n,a,y);
  WriteTxt3(fn,n,xre,xim,yre,yim);
  DistCoh(n,1,1,xre,xim,yre,yim,dTotCoh,dReCoh,dImCoh);
  totdisttxt.Text:=floattostr(dTotCoh);
  redisttxt.Text:=floattostr(dreCoh);
  imdisttxt.Text:=floattostr(dimCoh);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    SimCoh(n,1,1,alpha,sx,sy,xre,xim,yre,yim,simTotCoh,simReCoh,simImCoh,xy);
    totsimtxt.Text:=floattostr(simTotCoh);
    resimtxt.Text:=floattostr(simreCoh);
    imsimtxt.Text:=floattostr(simimCoh);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJ3SimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJ3SimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    totsimtxt.Text:=floattostr(-r);
    resimtxt.Text:=floattostr(simreCoh);
    imsimtxt.Text:=floattostr(simimCoh);
  end;

  showmessage('finished');
end;

procedure TfSimCorr.CmplxLogButtClick(Sender: TObject);
var eps,a,sx,sy,r,are,aim,dTotCoh,dReCoh,dImCoh:real; fn:string;
    p1,d1:vvector;
begin
  fn:='CmplxLog.txt';
  n:=900;
  are:=strtofloat(retxt.Text);
  aim:=strtofloat(imtxt.Text);
  CmplxLog(n,are,aim,xre,xim,yre,yim);
  a:=strtofloat(noisetxt.Text);
  AddNoise2(n,a,y);
  WriteTxt3(fn,n,xre,xim,yre,yim);
  DistCoh(n,1,1,xre,xim,yre,yim,dTotCoh,dReCoh,dImCoh);
  totdisttxt.Text:=floattostr(dTotCoh);
  redisttxt.Text:=floattostr(dreCoh);
  imdisttxt.Text:=floattostr(dimCoh);
  alpha:=strtofloat(distexptxt.Text);
  sx:=strtofloat(xweighttxt.Text);
  sy:=strtofloat(yweighttxt.Text);

  if not MaxSxyChk.Checked then begin
    SimCoh(n,1,1,alpha,sx,sy,xre,xim,yre,yim,simTotCoh,simReCoh,simImCoh,xy);
    totsimtxt.Text:=floattostr(simTotCoh);
    resimtxt.Text:=floattostr(simreCoh);
    imsimtxt.Text:=floattostr(simimCoh);
  end else begin
    eps:=0.1;
    p1[1]:=0; p1[2]:=0;
    d1[1]:=1; d1[2]:=1;
    HookeJeeves(2,p1,d1,fHJ3SimCov,r,eps);
    HookeJeeves(2,p1,d1,fHJ3SimCorr,r,eps);
    p1[1]:=exp(p1[1]); p1[2]:=exp(p1[2]);
    xweighttxt.Text:=floattostr(p1[1]);
    yweighttxt.Text:=floattostr(p1[2]);
    totsimtxt.Text:=floattostr(-r);
    resimtxt.Text:=floattostr(simreCoh);
    imsimtxt.Text:=floattostr(simimCoh);
  end;

  showmessage('finished');
end;

end.

