taucrp.m

function X = taucrp(varargin)
%TAUCRP   Creates a close returns plot.
%    R=TAUCRP(X [,Y] [,param1,param2,...]) creates a cross 
%    recurrence plot/ recurrence plot R for a limited range of
%    past and future states, also known as close returns plot.
%
%    R=TAUCRP(X [,Y],M,T,E,W) uses the dimension M, delay T, the size 
%    of neighbourhood E and the range W of past and future time
%    steps.
%
%    If X and Y are multi-column vectors then they will be 
%    considered as phase space vectors (TAUCRP can be used
%    for real phase space vectors without embedding).
%    
%    Parameters: dimension M, delay T, the size of neighbourhood 
%    E and the range W are the first four numbers after the data 
%    series; further parameters can be used to switch between 
%    various methods of finding the neighbours of the phasespace 
%    trajectory and to suppress the normalization of the data.
%
%    Methods of finding the neighbours/ of plot.
%      maxnorm     - Maximum norm.
%      euclidean   - Euclidean norm.
%      minnorm     - Minimum norm.
%      rr          - Maximum norm, fixed recurrence rate.
%      fan         - Fixed amount of nearest neighbours.
%      distance    - Distance coded matrix (global CRP, Euclidean norm).
%
%    Normalization of the data series.
%      normalize   - Normalization of the data.
%      nonormalize - No normalization of the data.
%
%    Parameters not needed to be specified.
%
%
%    Examples: a = sin((1:1000) * 2 * pi/67);
%              w = 160;
%              X = taucrp(a,2,17,.2,w,'nonorm','euclidean');
%              imagesc(1:size(X,2),-w:w,X), colormap([1 1 1; 0 0 0])
%
%    See also CRP, CRP2, CRP_BIG, JRP, CRQA.
%
%    References: 
%    Marwan, N., Romano, M. C., Thiel, M., Kurths, J.: 
%    Recurrence Plots for the Analysis of Complex Systems, Physics 
%    Reports, 438(5-6), 2007.

% Copyright (c) 2008-2009
% Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
% http://www.pik-potsdam.de
%
% Copyright (c) 2008
% Norbert Marwan, Potsdam University, Germany
% http://www.agnld.uni-potsdam.de
%
% $Date$
% $Revision$
%
% $Log$
% Revision 5.2  2009/03/24 08:33:47  marwan
% copyright address changed
%
% Revision 5.1  2008/01/25 12:47:26  marwan
% initial import
%


warning off
global errcode props nonorm

errcode=0;

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

init_properties
m_init = 1;
tau_init = 1;
eps_init = 0.1;
w_init = 100;

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

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

check_meth={'ma','eu','mi','rr','fa','di'}; 	% maxnorm, euclidean, nrmnorm,  fan, distance
check_norm={'non','nor'};				% nonormalize, normalize
check_gui={'gui','nog','sil'};				% gui, nogui, silent

if isnumeric(varargin{1}) 		% read commandline input
   varargin{9}=[];
   % transform any int to double
   intclasses = {'uint8';'uint16';'uint32';'uint64';'int8';'int16';'int32';'int64'};
   flagClass = [];
   for i = 1:length(intclasses)
       i_int=find(cellfun('isclass',varargin,intclasses{i}));
       if ~isempty(i_int)
           for j = 1:length(i_int)
               varargin{i_int(j)} = double(varargin{i_int(j)});
           end
           flagClass = [flagClass; i_int(:)];
       end
   end
   if ~isempty(flagClass)
       disp(['Warning: Input arguments at position [',num2str(flagClass'),'] contain integer values']);
       disp(['(now converted to double).'])
   end
   i_double=find(cellfun('isclass',varargin,'double'));
   i_char=find(cellfun('isclass',varargin,'char'));

   % check the text input parameters for method, gui and normalization
   temp_meth=0;
   temp_norm=0;
   temp_gui=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_norm=temp_norm+strcmpi(varargin{i_char(i)}(1:3),check_norm'); 
         temp_gui=temp_gui+strcmpi(varargin{i_char(i)}(1:3),check_gui'); 
      end
      method=min(find(temp_meth));
      nonorm=min(find(temp_norm))-1;
      nogui=min(find(temp_gui))-1;
      if isempty(method), method=1; end
      if isempty(nonorm), nonorm=1; end
      if isempty(nogui), nogui=0; end
      if method>length(check_meth), method=length(check_meth); end
      if nonorm>1, nonorm=1; end
      if nogui>2, nogui=2; end
   else
      method=1; nonorm=1; nogui=0;
   end
   if nogui==0 & nargout>0, nogui=1; 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 (isnumeric(varargin{2}) & max(size(varargin{2}))==1) | ~isnumeric(varargin{2})
       y=x;
       if ~isempty(varargin{i_double(2)}), m0=varargin{i_double(2)}(1); else m0=m_init; end
       if ~isempty(varargin{i_double(3)}), t=varargin{i_double(3)}(1); else t=tau_init; end
       if ~isempty(varargin{i_double(4)}), e=varargin{i_double(4)}(1); else e=eps_init; end
       if ~isempty(varargin{i_double(5)}), w=varargin{i_double(5)}(1); else w=w_init; end
     else
       if ~isempty(varargin{i_double(3)}), m0=varargin{i_double(3)}(1); else m0=m_init; end
       if ~isempty(varargin{i_double(4)}), t=varargin{i_double(4)}(1); else t=tau_init; end
       if ~isempty(varargin{i_double(5)}), e=varargin{i_double(5)}(1); else e=eps_init; end
       if ~isempty(varargin{i_double(6)}), w=varargin{i_double(6)}(1); else w=w_init; end
     end
     t=round(t); m0=round(m0); mflag=method;
     if e<0, e=1; disp('Warning: The threshold size E cannot be negative and is now set to 1.'), end
     if t<1, t=1; disp('Warning: The delay T cannot be smaller than one and is now set to 1.'), end
     if m0 < 1, m0 = 1; end
     if t < 1, t = 1; end
     if size(x,1)==1, x=x'; end, if size(y,1)==1, y=y'; end 
     m=max([size(x,2) size(y,2)]);
     if w > size(x,1); w = size(x,1); end

     if method==8 & (m*m0) > 1, 
       m0=1; 
       error(['The neighbourhood criterion ''Oder matrix''',10,'is not implemented - use crp or crp_big instead.'])
     end
     if method==9 & (m*m0) == 1, 
       m0=2; 
         disp(['Warning: For order patterns recurrence plots the dimension must',10,...
              'be larger than one. ',...
              'Embedding dimension is set to ',num2str(m0),'.'])
     end
     action='init';

  if ~isempty(find(isnan(x)))
     disp('NaN detected (in first variable) - will be cleared.')
     for k=1:size(x,2),  x(find(isnan(x(:,k))),:)=[]; end
  end
  if ~isempty(find(isnan(y)))
     disp('NaN detected (in second variable) - will be cleared.')
     for k=1:size(y,2),  y(find(isnan(y(:,k))),:)=[]; end
  end
  if size(x,1) < t*(m-1)+1 | size(y,1) < t*(m-1)+1
     error(['Too less data',10,...
            'Either too much NaN or the number of columns in the vectors do not match.'])
  end

    Nx=size(x,1); Ny=size(y,1);
    NX=Nx-t*(m0-1);NY=Ny-t*(m0-1);
    x0=zeros(Nx,m);y0=zeros(Ny,m);
    x0(1:size(x,1),1:size(x,2))=x; 
    y0(1:size(y,1),1:size(y,2))=y; 

    if nonorm==1, 
	 x=(x0-repmat(mean(x0),Nx,1))./repmat(std(x0),Nx,1);
	 y=(y0-repmat(mean(y0),Ny,1))./repmat(std(y0),Ny,1);
    end

  if ~isempty(find(isnan(x))), for k=1:size(x,2),  x(find(isnan(x(:,k))),:)=[]; end, end
  if ~isempty(find(isnan(y))), for k=1:size(y,2),  y(find(isnan(y(:,k))),:)=[]; end, end
  if size(x,1) < t*(m0-1)+1 | size(y,1) < t*(m0-1)+1
     error(['Too less data',10,...
            'Either too much NaN or the number of columns in the vectors do not match.'])
  end

else
  error('No valid input given!')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if m0>1
  x2=x(1:end-t*(m0-1),:);
  y2=y(1:end-t*(m0-1),:);
  for i=1:m0-1,
    x2(:,m*i+1:m*(i+1))=x(1+t*i:end-t*(m0-i-1),:);
    y2(:,m*i+1:m*(i+1))=y(1+t*i:end-t*(m0-i-1),:);
  end
  x=x2; y=y2; Nx=size(x,1); Ny=size(y,1);
  m=m0*m; clear x2 y2
end

x1=repmat(x,1,Ny);
for mi=1:m, x2(:,mi)=reshape(rot90(x1(:,0+mi:m:Ny*m+mi-m)),Nx*Ny,1); end
y1=repmat(y,Nx,1); x1=x2; clear x2


% create embedding vectors
NX = Nx - t*(m0-1); NY = Ny - t*(m0-1);

[idx add] = meshgrid(linspace(0,(m0-1)*t,m0), 1:Nx-(m0-1)*t);
idx = idx + add;
x1 = []; y1 = [];
for i = 1:m
    x1 = [x1 reshape(x(idx,i),size(idx,1),size(idx,2))];
    y1 = [y1 reshape(y(idx,i),size(idx,1),size(idx,2))];
end


% indices for which the recurrences should be computed
[i1 i2]  = meshgrid(1:NX,-w:w);
i2 = i2 + i1;

i_remove = i2 > length(y1) | i2 < 1;
i2(i_remove) = 1;

% compute distance matrix
D = (x1(i1,:) - y1(i2,:));

switch method
    case 1
        D2 = max(abs(D),[],2);
        X = double(D2 < e);
    case 2
        D2 = sqrt(sum(D.^2, 2));
        X = double(D2 < e);
    case 3
        D2 = sum(abs(D), 2);
        X = double(D2 < e);
    case 4
        D2 = max(abs(D),[],2);
        SS = sort(D2(:));
        idx = ceil(e * length(SS));
        e_ = SS(idx);
        X = double(D2 < e_);
    case 5
        D2 = sqrt(sum(D.^2, 2));
        D3 = reshape(D2,size(i1,1),size(i1,2));
        [SS, JJ] = sort(D3,1); JJ = JJ';
        mine = round((2*w+1)*e);
        X1(NX*(2*w+1)) = 0; 
        X1(JJ(:,1:mine)+repmat([0:(2*w+1):NX*(2*w+1)-1]',1,mine)) = 1;
        X = reshape(X1,(2*w+1),NX); 
    case 6
        D2 = sqrt(sum(D.^2, 2));
        X = double(D2);
      
    end

clear X1 SS JJ s px D D_ D2 D3
            
X = reshape(X,size(i1,1),size(i1,2)); X(i_remove) = 0;
%imagesc(X)