complete recoding based on an implementation by Maik Riedl (much faster calculations)

twinsurr.m

@@ -14,6 +14,11 @@ function [y nTw] = twinsurr(varargin)
 %    [Y,NTWINS]=TWINSURR(...) where NTWINS is the total number
 %    of twins in the RP.
 %
+%    Note: Please check before use the recurrence parameters M, T, and
+%          E by visual inspection of the recurrence plot! Please also
+%          check the number of twins NTWINS whether enough twins have
+%          been found.
+%
 %    Example: x = rand(3,1);
 %             a = [.8 .3 -.25 .9]';
 %             for i = 4:1000,
@@ -28,8 +33,8 @@ function [y nTw] = twinsurr(varargin)
 %    Twin Surrogates to Test for Complex Synchronisation, 
 %    Europhys. Lett., 75, 2006.
 
-% Copyright (c) 2008-2009
-% Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
+% Copyright (c) 2008-2016
+% Norbert Marwan, Maik Riedl, Potsdam Institute for Climate Impact Research, Germany
 % http://www.pik-potsdam.de
 %
 % Copyright (c) 2008
@@ -40,6 +45,9 @@ function [y nTw] = twinsurr(varargin)
 % $Revision$
 %
 % $Log$
+% Revision 5.4  2009/03/24 08:33:47  marwan
+% copyright address changed
+%
 % Revision 5.3  2009/03/17 09:18:17  marwan
 % serious bug fix, zero columns were also considered as twins!
 %
@@ -60,8 +68,8 @@ sil = 0;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input
 
-error(nargchk(1,7,nargin));
-if nargout>2, error('Too many output arguments'), end
+narginchk(1,7);
+nargoutchk(0,2)
 
 
 
@@ -88,10 +96,10 @@ end
 varargin{8}=[];
 i_double=find(cellfun('isclass',varargin,'double'));
 i_char=find(cellfun('isclass',varargin,'char'));
-check_meth={'ma','eu','mi','nr','rr','fa','in','om','op','di'}; 	% maxnorm, euclidean, nrmnorm,  fan, distance
+check_meth={'ma','eu','mi'}; 	% maxnorm, euclidean, nrmnorm,  fan, distance
 check_norm={'non','nor'};                        % nonormalize, normalize
 check_sil={'ve','si'};                         % verbose, silent
-
+nonorm = 0;
 if isnumeric(varargin{1}) 		% read commandline input
    % check the text input parameters for method, gui 
     temp_meth=0;
@@ -110,7 +118,7 @@ if isnumeric(varargin{1}) 		% read commandline input
          for i=1:length(i_char); temp2(i,:)=varargin{i_char(i)}(1:3); end
 
          if isempty(sil), sil=0; end
-         if isempty(nonorm), nonorm=1; end
+         if isempty(nonorm), nonorm=0; end
          if nonorm>1, nonorm=1; end
          if isempty(method_n), method_n=1; end
          if method_n>length(check_meth), method0=length(check_meth); end
@@ -140,52 +148,161 @@ end
 
 
 if size(x,1) < size(x,2), x = x'; end
-m0 = size(x,2);
+if size(x,2) > 1
+    warning('Multicolumn vectors are not supported. Only first column will be used.')
+    x = x(:,1);
+end
+
 N = length(x); 
 
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% make surrogates
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% main part: create surrogates
+
+
+surr = zeros(N,nsur);
 
-X = crp(x,m,t,e,method,norm_str,'sil');
-X = double(X) .* (1-eye(size(X)));
-NX = length(X);
-%% find twins
-if ~sil, h = waitbar(0,'Searching Twins'); end
-for i = 1:NX, if ~sil, waitbar(i/NX); end
-    A = repmat(X(:,i),1,NX);
-%    S{i} = find(all(X == A & X == 1));
-    S{i} = [i find(all(X == A) & any(X(:,i)))];
-    nTwins(i) = numel(S{i});
-end
-if ~sil, delete(h); end
 
-nTwins = sum(nTwins) - NX;
+%% find twins
+h = [];
+if ~sil & N > 3000, h = waitbar(0,'Find recurrences'); end
+twinList = findTwins(x,m,t,e,method,nonorm,sil,h);
+nTwins = length(twinList); % number of twins 
 
 if nTwins < 1
-    warning(['No twins found!',10,'The derived surrogates are identical with the original data.']')
+    error(['No twins found!',10,'The derived surrogates would identical with the original data.']')
+end
+
+%% marker list for number of twins
+markerTwin = zeros(N,1);
+for i = 1:nTwins
+    markerTwin(twinList{i}) = i;
 end
 
+%% remove the last twin from twinList                  
+temp_ind = find(markerTwin);
+temp_ind2 = find(twinList{markerTwin(temp_ind(end))} == temp_ind(end));
+twinList{markerTwin(temp_ind(end))}(temp_ind2) = [];
 
-for k = 1:nsur
-    %% chose randomly first point
-    i = ceil(NX * rand);
-    xs = x(i,:);
-
-    %% jump randomly between twins until length of surrogate is reached
-    while length(xs) < N
-        j = S{i};
-        j_i = ceil(length(j) * rand);
-        i = j(j_i);
-        i = i + 1;
-        if i > NX
-            i = ceil(NX * rand);
-            continue 
+%% create twin surrogate
+cnt = 0;
+while cnt < nsur
+    idxSurr = ceil(rand(1)*N); % random starting point
+    r_ind = rand(N,1);
+    i = 2;
+    while (idxSurr(i-1)+1) <= N % just go through the time series (the next point after the last found point)
+        if markerTwin(idxSurr(i-1)+1) ~= 0 % we have twins: let's jump!
+            nr_sel = ceil(r_ind(i) * length(twinList{markerTwin(idxSurr(i-1)+1)})); % select randomly a new point
+            idxSurr(i)=twinList{markerTwin(idxSurr(i-1)+1)}(nr_sel); % jump to this new point
+        else
+           idxSurr(i) = idxSurr(i-1)+1; % just use the next point from the time series
+        end
+        i = i + 1; % just go through the time series
+        if i > N
+            break
+        end
+    end
+    if length(idxSurr)==N
+        cnt = cnt + 1;
+        surr(:,cnt)=x(idxSurr);
+        if ~sil & exist('h') & mod(cnt,100)==0
+            waitbar((2*N+N*cnt/nsur)/(3*N),h,'Create surrogates')
         end
-        xs = [xs; x(i,:)];
     end
-    y(:,k,1:m0) = reshape(xs(:),N,1,m0);
 end
 
+if ~sil & exist('h') & ishandle(h), close(h), end
+
+
+
+if nargout > 0;
+   y = surr;
+else
+   surr
+end
+
+
 if nargout == 2;
     nTw = nTwins;
 end
+
+
+
+
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+
+function twinList=findTwins(x,m,t,thr,method,nonorm,sil,h)
+% FINDTWINS   Find twins in recurrence plot.
+%
+%    twinList=TWIN_MAIK(X,M,TAU,THR)
+%
+%The function look for the twins in a time series X.
+%Therefore the phase space is reconstructed by time delay embedding
+%approach where M is the dimension of the phase space and TAU is
+%the time delay.
+%Neighbouring points are two points with a distance smaller than the
+%threshold THR.
+%Twins are points with the same neighbourhood.
+
+%% Reconstruct phase space
+if nonorm
+    if ~sil, disp('Normalize data.'), end
+    x = zscore(x);
+end
+
+updateTime = 100;
+N = length(x);
+Nx = N - (m-1)*t;
+xEmb = zeros(Nx,m);                
+for i = 1:m
+    xEmb(:,i) = x(((i-1)*t+1):(end-t*(m-i)));
+end
+
+%% Find recurrence points (= neighbours) for each state
+% this is a list of neighbours for each time point/state
+idxNeighbours = [];
+nNeighbours=[];
+for i = 1:Nx, if ~sil & ~isempty(h) & mod(i,updateTime)==0, waitbar(i/(3*N)), end
+    % distance between current state and all other states
+    switch method
+        case 'eu'
+           d = sqrt(sum((xEmb(1:end,:)-xEmb(i,:)).^2,2));
+        case 'mi'
+           d = sum(abs(xEmb(1:end,:)-xEmb(i,:)),2);
+        otherwise
+           d = max(abs(xEmb(1:end,:)-xEmb(i,:)),[],2);
+    end
+    % indices of the Neighbours
+    idxNeighbours{i} = find(d<thr);
+    % number of Neighbours found (useful for acceleration of the algorithm)
+    nNeighbours(i) = length(idxNeighbours{i});
+end
+
+%% Find twins
+ind_p = 1:Nx; % set of states
+twinList = []; % list of twins 
+cnt = 1;
+while length(ind_p) > 1,if ~sil & ~isempty(h) & mod(i,updateTime)==0, waitbar((2*N - length(ind_p))/(3*N),h,'Find twins'), end
+    currentTwins = 1; % start with the first entry in the (remaining) set 
+    % find only those columns in the RP that have the same number of neighbours
+    potentialTwins = find(nNeighbours == nNeighbours(1)); % this will speed up the code
+    for i = 2:length(potentialTwins)
+        if sum(abs(idxNeighbours{1}-idxNeighbours{potentialTwins(i)})) == 0 & any(idxNeighbours{1})
+            currentTwins = [currentTwins, potentialTwins(i)];
+        end
+    end
+    if length(currentTwins) > 1
+        twinList{cnt} = ind_p(currentTwins);
+        cnt = cnt + 1;
+    end
+    % remove all states that are already found
+    ind_p(currentTwins) = [];
+    idxNeighbours(currentTwins) = [];
+    nNeighbours(currentTwins) = [];
+end
+
+if ~sil & ~isempty(h), waitbar(2/3,h), end