foldername = '~/Documents/MATLAB/ecgdata/dataBSPMforOlivier';
load(fullfile(foldername, 'invalidElectrodes.mat'));
listing = dir(fullfile(foldername, '*_QRST.mat'));
% listing = dir(fullfile(foldername, '*_PP.mat'));
filenameList = fullfile(foldername, {listing.name}');
shortFilenameList = {listing.name}';

load BSPMoutcome


Index_to_CRF=[];
for k=1:75,
st=shortFilenameList{k};

po=0;
while (length(str2num(st(10:10+po))) ~= 0)
    po=po+1;
end
Index_to_CRF(k)=str2num(st(10:10+po-1));
end


indrecur=[];
indnonrecur=[];

for k=1:75,
if (BSPMOutcome(Index_to_CRF(k),3) == 1) && (BSPMOutcome(Index_to_CRF(k),5) == 1)
%if (BSPMOutcome(Index_to_CRF(k),3) == 1) 
%       indsucCV=indsucCV+1;
        if (BSPMOutcome(Index_to_CRF(k),4) == 1)
            indrecur=[indrecur;k];
        end
        if (BSPMOutcome(Index_to_CRF(k),4) == 0)
            indnonrecur=[indnonrecur;k];
        end 
    end
end

load RT3

RE=[];
REi=[];
PER=[];
POS=[];
PIS=[];

jj=2;
RF=[];
for k=1:75,
    for l=1:48/jj,
        RF(k,:,l)=mean(RT(k,:,jj*l-(jj-1):jj*l),3);
    end;
end

for subj=1:75,
for k=1:size(RF,3),
[I,J]=max(abs(fft(RF(subj,1:500,k)-mean(RF(subj,1:500,k)))));
Per=1/(((J-1)/250)*0.5);
PER(subj,k)=Per;
posmax=1;
posmin=1;
    for ll=1:7,
    [I,J]=max(RF(subj,posmax+round(Per/2):posmax+3*round(Per/2),k));
    [I,Ji]=min(RF(subj,posmax:posmax+round(Per/2)+J-1,k));
    PIS(subj,ll,k)=posmax+Ji-1;
    
    y=RF(subj,PIS(subj,ll,k)-1:PIS(subj,ll,k)+1,k);
    c=y(2);a=(y(1)+y(3)-2*c)/2;b=(y(3)-y(1))/2;
    deltatemps=-b/(2*a);
    PIS(subj,ll,k)=PIS(subj,ll,k)+deltatemps;
    
    %REi(subj,ll,k)=RF(subj,posmax+Ji-1,k);
    REi(subj,ll,k)=-b^2/a+c;
    
    posmax=posmax+round(Per/2)+J-1;
    
    y=RF(subj,posmax-1:posmax+1,k);
    c=y(2);a=(y(1)+y(3)-2*c)/2;b=(y(3)-y(1))/2;
    deltatemps=-b/(2*a);
    
    POS(subj,ll,k)=posmax+deltatemps;
    
    %RE(subj,ll,k)=RF(subj,posmax,k);
    RE(subj,ll,k)=-b^2/a+c;
    
    
  
    end
end
end

RFF=[];
for l=1:size(RF,3),
    for k=1:75,
%        RF(k,:,l)=filtfilt(B,A,RT(k,:,l));
        RFF(k,:,l)=(abs(hilbert(RF(k,:,l))));
    end;
end

C=median(RFF(:,1:500,:),3);
C=diag(1./median(C(:,150:450),2))*C;
CC=cumsum(C,2);

c=[];
for k=1:75,c(k)=C(k,round(mean(POS(k,1,:),3)));end
[p,h]=ranksum(c(indrecur),c(indnonrecur))
%plot(mean(PER(:,:),2),c,'o')

in=[indrecur;indnonrecur];
%in=[1:75];
labels=zeros(63,1);
labels(1:length(indnonrecur))=1;

%[T,TT]=corrcoef(mean(PER(in,:),2),c(in))

disp('with first max');

[T,TT]=corrcoef(mean(POS(in,1,:),3),mean(RE(in,1,:),3)./mean(mean(abs(RF(in,250:450,:)),2),3))
[p,h]=ranksum(mean(RE(indrecur,1,:),3)./mean(mean(abs(RF(indrecur,150:450,:)),2),3),mean(RE(indnonrecur,1,:),3)./mean(mean(abs(RF(indnonrecur,150:450,:)),2),3))
%[T,TT]=corrcoef(mean(POS(in,1,:),3),mean(mean(abs(RF(in,150:450,:)),2),3))
[X,Y,T,AUC,OPTROCPT]=perfcurve(labels,mean(RE([indnonrecur;indrecur],1,:),3)./mean(mean(abs(RF([indnonrecur;indrecur],150:450,:)),2),3),1);
AUC
acc=((1-OPTROCPT(1))*32+OPTROCPT(2)*31)/(31+32)
%[X,Y,T,AUC,OPTROCPT]=perfcurve(labels,c([indnonrecur;indrecur]),1);
%AUC
%[X,Y,T,AUC,OPTROCPT]=perfcurve(labels,mean(PER([indnonrecur;indrecur],:),2),1);
%AUC

disp('with first min');

[T,TT]=corrcoef(mean(PIS(in,1,:),3),abs(mean(REi(in,1,:),3))./mean(mean(abs(RF(in,250:450,:)),2),3))
[p,h]=ranksum(abs(mean(REi(indrecur,1,:),3)./mean(mean(abs(RF(indrecur,150:450,:)),2),3)),abs(mean(REi(indnonrecur,1,:),3)./mean(mean(abs(RF(indnonrecur,150:450,:)),2),3)))
[X,Y,T,AUC,OPTROCPT]=perfcurve(labels,abs(mean(REi([indnonrecur;indrecur],1,:),3))./mean(mean(abs(RF([indnonrecur;indrecur],150:450,:)),2),3),1);
AUC
acc=((1-OPTROCPT(1))*32+OPTROCPT(2)*31)/(31+32)

for k=1:63, 
    [T,TT]=corrcoef(squeeze(REi(in(k),1,:)),squeeze(PIS(in(k),1,:)));
    t(k)=T(1,2);
    pt(k)=TT(1,2);
end


%[T,TT]=corrcoef(mean(PIS(in,1,:),3),mean(mean(abs(RF(in,250:450,:)),2),3))
%for k=1:63, [T,TT]=corrcoef(squeeze(POS(in(k),1,:)),squeeze(RE(in(k),1,:)));t(k)=T(1,2);pt(k)=TT(1,2);end
%for k=1:63, [T,TT]=corrcoef(squeeze(PIS(in(k),1,:)),squeeze(abs(REi(in(k),1,:))));t(k)=T(1,2);pt(k)=TT(1,2);end
