crqad_big.m

function out=crqad_big(varargin)
% CRQAD_BIG   Computes and plots the diagonalwise CRQA measures.
%    Y=CRQAD_BIG(X [,Y] [,param1,param2,...]) 
%    Recurrence quantification analysis of diagonals in the
%    cross recurrence plot of the vectors X and Y as well as
%    X and -Y. The output is a structure (see below).
%
%    The input vectors can be multi-column vectors, where
%    each column will be used as a component of the 
%    phase-space vector. However, if the first column is
%    monotonically increasing, it will be used as an
%    time scale for plotting.
%
%    Y=CRQAD(X,M,T,E,W) computes the recurrence
%    quantification analysis of the recurrence plot
%    of X by using the dimension M, delay T, the
%    size of neighbourhood E, for the diagonals within
%    the range [-W,W] around the main diagonal.
%
%    Parameters: dimension M, delay T, the size of
%    neighbourhood E and the range size W are the first 
%    five numbers after the data series; if W is empty,
%    the whole plot will be calculated. Further parameters 
%    can be used to switch between various methods of finding
%    the neighbours of the phasespace trajectory, to suppress 
%    the normalization of the data and to suppress the GUI
%    (useful in order to use this programme by other programmes).
%
%    Methods of finding the neighbours.
%      maxnorm     - Maximum norm.
%      euclidean   - Euclidean norm.
%      minnorm     - Minimum norm.
%      nrmnorm     - Euclidean norm between normalized vectors
%                    (all vectors have the length one).
%      fan         - Fixed amount of nearest neighbours.
%      inter       - Interdependent neighbours.
%      omatrix     - Order matrix.
%      opattern    - Order patterns recurrence plot.
%
%    Normalization of the data series.
%      normalize   - Normalization of the data.
%      nonormalize - No normalization of the data.
%
%    Parameters not needed to be specified.
%
%    Output:
%      Y.RRp  = RRp
%      Y.RRm  = RRm
%      Y.DETp = DETp
%      Y.DETm = DETm
%      Y.Lp   = Lp
%      Y.Lm   = Lm
%      
%    Examples: a = sin(0:.1:80) + randn(1,801);
%              b = sin(0:.1:80) + randn(1,801);
%              crqad_big(a,b,3,15,.1,100,'fan')
%
%    See also CRQA, CRQAD, CRP, CRP2, CRP_BIG, DL, TT.
%
%    References: 
%    Marwan, N., Kurths, J.:
%    Nonlinear analysis of bivariate data with cross recurrence plots,
%    Phys. Lett. A, 302, 2002.

% Copyright (c) 2002-2006 by AMRON
% Norbert Marwan, Potsdam University, Germany
% http://www.agnld.uni-potsdam.de
%
%
% $Date$
% $Revision$
%
% $Log$
%
%
% This program is part of the new generation XXII series.
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or any later version.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties

global props

init_properties

lmin=3;
w=[]; method='max'; method=1; t=1; m=1; e=.1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input

error(nargchk(1,10,nargin));
if nargout>1, error('Too many output arguments'), end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% splash the GPL

splash_gpl('crp');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check and read the input

   varargin{11}=[];
   i_double=find(cellfun('isclass',varargin,'double'));
   i_char=find(cellfun('isclass',varargin,'char'));
   nogui=0;

   if nargin & isnumeric(varargin{1})

     % check the text input parameters for method, gui 
      check_meth={'ma','eu','mi','nr','fa','in','om','op','di'}; 	% maxnorm, euclidean, nrmnorm,  fan, distance
      check_gui={'gui','nog','sil'};		% gui, nogui, silent
      check_norm={'non','nor'};				% nonormalize, normalize
      temp_meth=0;
      temp_gui=0;
      temp_norm=0;
      if ~isempty(i_char)
         for i=1:length(i_char), 
            varargin{i_char(i)}(4)='0';
            temp_meth=temp_meth+strcmpi(varargin{i_char(i)}(1:2),check_meth'); 
            temp_gui=temp_gui+strcmpi(varargin{i_char(i)}(1:3),check_gui'); 
            temp_norm=temp_norm+strcmpi(varargin{i_char(i)}(1:3),check_norm'); 
         end
         method=min(find(temp_meth));
         nonorm=min(find(temp_norm))-1;
         nogui=min(find(temp_gui))-1;
         for i=1:length(i_char); temp2(i,:)=varargin{i_char(i)}(1:3); end
         i_char(strmatch(check_gui(find(temp_gui)),temp2))=[];
         if isempty(nogui), nogui=0; end
         if isempty(method), method=1; end
         if nonorm>1, nonorm=1; end
         if nogui>2, nogui=1; end
         if method>length(check_meth), method=length(check_meth); end
      else
         nogui=0; nonorm=1; 
         if nargout
            nogui=1;
            action='compute'; 
         end
      end
      if nogui==0
        action='init';
      else
        action='compute'; 
      end
      
      % get the parameters for creating RP
      if max(size(varargin{1}))<=3
         error('To less values in data X.')
      end
      x=double(varargin{1});
      if isempty(varargin{2}) | ~isnumeric(varargin{2}), y=x; else
      y=double(varargin{2}); end
      if sum(double(diff(x(:,1))<=0)), embed_flag=0; end
   
      if (isnumeric(varargin{2}) & max(size(varargin{2}))==1) | ~isnumeric(varargin{2})
        y=x;
        if ~isempty(varargin{i_double(2)}), m=varargin{i_double(2)}(1); else m=1; end
        if ~isempty(varargin{i_double(3)}), t=varargin{i_double(3)}(1); else t=1; end
        if ~isempty(varargin{i_double(4)}), e=varargin{i_double(4)}(1); else e=.1; end
        if ~isempty(varargin{i_double(5)}), w=varargin{i_double(5)}(1); else w=varargin{i_double(5)}; end
%        if ~isempty(varargin{i_double(6)}), wstep=varargin{i_double(6)}(1); else wstep=1; end
      else
        if ~isempty(varargin{i_double(3)}), m=varargin{i_double(3)}(1); else m=1; end
        if ~isempty(varargin{i_double(4)}), t=varargin{i_double(4)}(1); else t=1; end
        if ~isempty(varargin{i_double(5)}), e=varargin{i_double(5)}(1); else e=.1; end
        if ~isempty(varargin{i_double(6)}), w=varargin{i_double(6)}(1); else w=varargin{i_double(6)}; end
%        if ~isempty(varargin{i_double(7)}), wstep=varargin{i_double(7)}(1); else wstep=1; end
      end
    else
      error('No valid arguments.')
    end

   
    Nx=length(x); Ny=length(y);
    if size(x,1)<size(x,2), x=x'; end
    if size(y,1)<size(y,2), y=y'; end
   
   if size(x,2)>=2
      xscale=x(:,1); 
      if ~isempty(find(diff(xscale)<0)), embed_flag=0;end
      if nonorm==1, x=(x(:,2)-mean(x(:,2)))/std(x(:,2)); else x=x(:,2); end
   else
      if nonorm==1, x=(x-mean(x))/std(x); end
      xscale=(1:length(x))'; 
   end
   if size(y,2)>=2
      yscale=y(:,1); 
      if ~isempty(find(diff(yscale)<0)), embed_flag=0;end
      if nonorm==1, y=(y(:,2)-mean(y(:,2)))/std(y(:,2)); else y=y(:,2); end
   else
       if nonorm==1, y=(y-mean(y))/std(y); end
       yscale=(1:length(y))';
   end
      
      if max(size(x))~=max(size(y)),
        if ~nogui, errordlg('Data must have the same length.','Check Data'), else error('Data must have the same length.'), end
      end
      if e<0, 
        e=1; 
	if ~nogui
	   warndlg('The threshold size E can not be negative and is now set to 1.','Check Data')
	   h=findobj('Tag','crqa_eps');
	   set(h(1),'String',str2num(e)) 
        else 
	   disp('The threshold size E can not be negative and is now set to 1.'), 
        end
      end
      if t<1, 
        t=1; 
	if ~nogui
	   warndlg('The delay T can not be smaller than one and is now set to 1.','Check Data')
	   h=findobj('Tag','crqa_maxLag');
	   set(h(1),'String',str2num(t)) 
	else
	   disp('The delay T can not be smaller than one and is now set to 1.')
	end
      end
      if isempty(w) & Nx > 5000, w = 100; wstep=1; end
      if isempty(w), w=.5*Nx; wstep=1; end
%      if w<2, 
%        w=2; 
%	if ~nogui, warndlg('The window size W exceeds the valid range.','Check Data')
%	else, disp('The window size W exceeds the valid range.'), end
%      end
      if w>Nx, 
        w=Nx; wstep=1;; 
	if ~nogui, warndlg('The window size W exceeds the valid range.','Check Data')
	else, disp('The window size W exceeds the valid range.'), end
      end
      t=round(t); m=round(m); w=round(w);% wstep=round(wstep); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute

flag=1;

x1=x; x2=y;
  if method > 3
      disp('Warning: RQA from choosen neighbourhood is not possible!')
      return
  end

warning off
if size(x1,1)<size(x1,2), x1=x1'; end
if size(x2,1)<size(x2,2), x2=x2'; end


  
  % embedding vectors
  NX=Nx-t*(m-1);NY=Ny-t*(m-1);  
  i=(1:NX)';j=0:t:0+(m-1)*t;
  i=reshape(i(:,ones(1,m))+j(ones(NX,1),:),m*NX,1);
  x1=x(i);
  x2=reshape(x1,NX,m);

  i=(1:NY)';j=0:t:0+(m-1)*t;
  i=reshape(i(:,ones(1,m))+j(ones(NY,1),:),m*NY,1);
  y1=y(i);
  y2=reshape(y1,NY,m);


% compute the diagonalwise RQA
clear DET L RR; j = 0;
hw = waitbar(0,'Compute RQA');
for i = 0:w,waitbar(i/(2*w))
clear z z0 z1

    if m > 1

        switch(method)

        %%%%%%%%%%%%%%%%% local CRP, fixed distance

          case 1
        %%%%%%%%%%%%%%%%% maximum norm
            s = max(abs(x2(1:NX-i,:) - y2(i+1:NY,:))');
          case 2
        %%%%%%%%%%%%%%%%% euclidean norm
            errcode=112;
            s = sqrt(sum((x2(1:NX-i,:) - y2(i+1:NY,:)).^2, 2));
          case 3
        %%%%%%%%%%%%%%%%% minimum norm
            errcode=113;
            s = sum(abs(x2(1:NX-i,:) - y2(i+1:NY,:))');
        end
    else
        s = abs(x2(1:NX-i,:) - y2(i+1:NY,:));
    end
    
    X = s(:) < e;
    

    if sum(X) == length(X), l1=length(X); else
        z=diff(X); z0 = [];
        if ~isempty(find(~(z-1))),z0(:,1)=find(~(z-1));else,z0(1:length(X))=0;end, 
        if ~isempty(find(~(z+1))),z1=find(~(z+1));else,z1(1:length(X))=0;end

        if z0(1)>z1(1)
            z0(2:end+1,1)=z0(1:end);z0(1)=0; 
        end
        if length(z0)>length(z1), z0(end)=[]; end
        l=sort(z1-z0); l1=l(find(l>lmin));
    end
    DET(j+w+1)=sum(l1)/sum(X);
    L(j+w+1)=mean(l1);
    RR(j+w+1)=sum(X)/length(X);

    if m > 1

        switch(method)

        %%%%%%%%%%%%%%%%% local CRP, fixed distance

          case 1
        %%%%%%%%%%%%%%%%% maximum norm
            s = max(abs(x2(i+1:NX,:) - y2(1:NY-i,:))');
          case 2
        %%%%%%%%%%%%%%%%% euclidean norm
            errcode=112;
            s = sqrt(sum((x2(i+1:NX,:) - y2(1:NY-i,:)).^2, 2));
          case 3
        %%%%%%%%%%%%%%%%% minimum norm
            errcode=113;
            s = sum(abs(x2(i+1:NX,:) - y2(1:NY-i,:))');
        end
    else
        s = abs(x2(i+1:NX,:) - y2(1:NY-i,:));
    end
    X = s(:) < e;

    if sum(X) == length(X), l1=length(X);else
        z=diff(X); z0 = [];
        if ~isempty(find(~(z-1))),z0(:,1)=find(~(z-1));else,z0(1:length(X))=0;end, 
        if ~isempty(find(~(z+1))),z1=find(~(z+1));else,z1(1:length(X))=0;end

        if z0(1)>z1(1)
            z0(2:end+1,1)=z0(1:end);z0(1)=0; 
        end
        if length(z0)>length(z1), z0(end)=[]; end
        l=sort(z1-z0); l1=l(find(l>lmin));
    end
    DET(w-j+1)=sum(l1)/sum(X);
    L(w-j+1)=mean(l1);
    RR(w-j+1)=sum(X)/length(X);

    j=j+1;

end


L(find(isnan(L)))=0;
RR(find(isnan(RR)))=0;
DET(find(isnan(DET)))=0;

if nargout, XCF=xcf(x1,x2,w,1); end

if ~nargout
subplot(2,2,1)
clim=1;
xcf(x,y,w)
set(gca,'fonta','i')
xlabel('Lag'), axis([-w w -clim clim])
ylabel('Cross Correlation')
h=text(0,0,'A','fontw','b');set(h,'un','pi'),set(h,'pos',[9,16,0])

switch flag
case 1
subplot(2,2,2), plot([-w:w],RR,'k','linew',.7), 
set(gca,'fonta','i'),axis([-w w 0 1])
xlabel('Lag'),ylabel('Recurrence Rate'),grid on
h=text(0,0,'B','fontw','b');set(h,'un','pi'),set(h,'pos',[9,16,0])

subplot(2,2,3), plot([-w:w],DET,'k','linew',.7), 
set(gca,'fonta','i'),xlabel('Lag')
axis([-w w 0 1]),ylabel('Determinism'),grid on
h=text(0,0,'C','fontw','b');set(h,'un','pi'),set(h,'pos',[9,16,0])

subplot(2,2,4), plot([-w:w],L,'k','linew',.7), 
set(gca,'fonta','i'),xlabel('Lag')
Lmax = max(L); if ~Lmax, Lmax = 1; end
axis([-w w 0 Lmax]),ylabel('Averaged Line Length'),grid on
h=text(0,0,'D','fontw','b');set(h,'un','pi'),set(h,'pos',[9,16,0])

case 2
subplot(2,2,2), plot([-w:w],smooth(RR,5,5),'k','linew',.7), 
set(gca,'fonta','i'),axis([-w w 0 1])
xlabel('Lag'),ylabel('Recurrence Rate'),grid on
h=text(0,0,'B','fontw','b');set(h,'un','pi'),set(h,'pos',[9,16,0])

subplot(2,2,3), plot([-w:w],smooth(DET,5,5),'k','linew',.7), 
set(gca,'fonta','i'),xlabel('Lag')
axis([-w w 0 1]),ylabel('Determinism'),grid on
h=text(0,0,'C','fontw','b');set(h,'un','pi'),set(h,'pos',[9,16,0])

subplot(2,2,4), plot([-w:w],smooth(L,5,5),'k','linew',.7), 
set(gca,'fonta','i'),xlabel('Lag')
axis([-w w 0 max(L)]),ylabel('Averaged Line Length'),grid on
h=text(0,0,'D','fontw','b');set(h,'un','pi'),set(h,'pos',[9,16,0])

end
else
  out.XCF=XCF';
  out.RRp=RR;
  out.DETp=DET;
  out.Lp=L;
end







% compute the diagonalwise RQA
clear DET L RR, j = 0;
for i = 0:w,waitbar((i+w)/(2*w))
clear z z0 z1

    if m > 1

        switch(method)

        %%%%%%%%%%%%%%%%% local CRP, fixed distance

          case 1
        %%%%%%%%%%%%%%%%% maximum norm
            s = max(abs(-x2(1:NX-i,:) - y2(i+1:NY,:))');
          case 2
        %%%%%%%%%%%%%%%%% euclidean norm
            errcode=112;
            s = sqrt(sum((-x2(1:NX-i,:) - y2(i+1:NY,:)).^2, 2));
          case 3
        %%%%%%%%%%%%%%%%% minimum norm
            errcode=113;
            s = sum(abs(-x2(1:NX-i,:) - y2(i+1:NY,:))');
        end
    else
        s = abs(-x2(1:NX-i) - y2(i+1:NY));
    end
      
    X = s(:) < e;

    if sum(X) == length(X), l1=length(X); else
        z=diff(X); z0 = [];
        if ~isempty(find(~(z-1))),z0(:,1)=find(~(z-1));else,z0(1:length(X))=0;end, 
        if ~isempty(find(~(z+1))),z1=find(~(z+1));else,z1(1:length(X))=0;end

        if z0(1)>z1(1)
            z0(2:end+1,1)=z0(1:end);z0(1)=0; 
        end
        if length(z0)>length(z1), z0(end)=[]; end
        l=sort(z1-z0); l1=l(find(l>lmin));
    end
    DET(j+w+1)=sum(l1)/sum(X);
    L(j+w+1)=mean(l1);
    RR(j+w+1)=sum(X)/length(X);

    if m > 1

        switch(method)

        %%%%%%%%%%%%%%%%% local CRP, fixed distance

          case 1
        %%%%%%%%%%%%%%%%% maximum norm
            s = max(abs(-x2(i+1:NX,:) - y2(1:NY-i,:))');
          case 2
        %%%%%%%%%%%%%%%%% euclidean norm
            errcode=112;
            s = sqrt(sum((-x2(i+1:NX,:) - y2(1:NY-i,:)).^2, 2));
          case 3
        %%%%%%%%%%%%%%%%% minimum norm
            errcode=113;
            s = sum(abs(-x2(i+1:NX,:) - y2(1:NY-i,:))');
        end
    else
        s = abs(-x2(i+1:NX,:) - y2(1:NY-i,:));
    end
    
    X = s(:) < e;

    if sum(X) == length(X), l1=length(X);else
        z=diff(X); z0 = [];
        if ~isempty(find(~(z-1))),z0(:,1)=find(~(z-1));else,z0(1:length(X))=0;end, 
        if ~isempty(find(~(z+1))),z1=find(~(z+1));else,z1(1:length(X))=0;end

        if z0(1)>z1(1)
            z0(2:end+1,1)=z0(1:end);z0(1)=0; 
        end
        if length(z0)>length(z1), z0(end)=[]; end
        l=sort(z1-z0); l1=l(find(l>lmin));
    end
    
    DET(w-j+1)=sum(l1)/sum(X);
    L(w-j+1)=mean(l1);
    RR(w-j+1)=sum(X)/length(X);

    j=j+1;

end
if ishandle(hw), delete(hw), end

L(find(isnan(L)))=0;
RR(find(isnan(RR)))=0;
DET(find(isnan(DET)))=0;

if ~nargout
switch flag
case 1
subplot(2,2,2), hold on, plot([-w:w],RR,'r','linew',.7), hold off, set(gca,'YLimMode','Auto')
subplot(2,2,3), hold on, plot([-w:w],DET,'r','linew',.7), hold off, hold off, set(gca,'YLimMode','Auto')
subplot(2,2,4), hold on, plot([-w:w],L,'r','linew',.7), hold off, hold off, set(gca,'YLimMode','Auto')
case 2
subplot(2,2,2), hold on, plot([-w:w],smooth(RR,5,5),'r','linew',.7), hold off, hold off, set(gca,'YLimMode','Auto')
subplot(2,2,3), hold on, plot([-w:w],smooth(DET,5,5),'r','linew',.7), hold off, hold off, set(gca,'YLimMode','Auto')
subplot(2,2,4), hold on, plot([-w:w],smooth(L,5,5),'r','linew',.7), hold off, hold off, set(gca,'YLimMode','Auto')
end
else
  out.RRm=RR;
  out.DETm=DET;
  out.Lm=L;
end
warning on