diff --git a/twinsurr.m b/twinsurr.m
index 1ba09252fa2a61fd5cf579b14f0e39a9792d49ac..ec94fc0ce524d19e6cf3e8fc5a73ebb92ed65ee7 100644
--- a/twinsurr.m
+++ b/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