adding Hirata's reconstruction method to recons.m

mi.m

@@ -7,13 +7,13 @@ function varargout=mi(varargin)
 %    between all pairs of X1,X2,...,Xn using 10 bins. The result I
 %    is a N x N matrix.
 %
-%    I=MI(X1,X2,...,Xn,N), where N is a scalar, determines the number
-%    of bins N.
-%
 %    I=MI(X1,X2,...,Xn,L) or I=MI(X1,X2,...,Xn,N,L), where L is a 
 %    scalar, computes mutual informations for shifting of the 
 %    components of each pair of X1,X2,...,Xn until a maximal lag L. 
 %
+%    I=MI(X1,X2,...,Xn,N,L), where N is a scalar determining the number
+%    of bins N.
+%
 %    [I S]=MI(...) computes the mutual information and the 
 %    standard error (only for one- and two-dimensional data).
 %
@@ -60,6 +60,10 @@ function varargout=mi(varargin)
 % $Revision$
 %
 % $Log$
+% Revision 3.18  2016/03/03 14:57:40  marwan
+% updated input/output check (because Mathworks is removing downwards compatibility)
+% bug in crqad_big fixed (calculation of xcf).
+%
 % Revision 3.17  2016/02/05 17:16:37  marwan
 % fix compatibility issues (R2014b)
 %

recons.m

@@ -3,12 +3,20 @@ function varargout=recons(varargin)
 %    Y = RECONS(X) reconstructs a time series Y from the 
 %    recurrence plot in the matrix X.
 %
-%    Y = RECONS(X,NAME) reconstructs the time series using
+%    Y = RECONS(X,...,METHOD) specifies the reconstruction method
+%    where string METHOD can be either 'thiel' or 'hirata', using
+%    the method by Marco Thiel (defeault) or Yoshito Hirata.
+%
+%    Y = RECONS(X,...,NAME) reconstructs the time series using
 %    the named cumulative distribution function, which can
-%    be 'norm' or 'Normal' (defeault), 'unif' or 'Uniform'.
+%    be 'norm' or 'Normal' (defeault), 'unif' or 'Uniform'. This
+%    is only necessary for the method by Thiel. The reconstruction
+%    from the Hirata will not be scaled.
 %
-%    Y = RECONS(X,P) reconstructs the time series using
-%    the cumulative distribution function given by vector P.
+%    Y = RECONS(X,P,...) reconstructs the time series using the
+%    cumulative distribution function given by vector P (should be
+%    2nd argument). This is only necessary for the method by Thiel. 
+%    The reconstruction from the Hirata will not be scaled.
 %
 %    See also CRP, CRP2, JRP, TWINSURR.
 %
@@ -16,6 +24,10 @@ function varargout=recons(varargin)
 %    Thiel, M., Romano, M. C., Kurths, J.: 
 %    How much information is contained in a recurrence plot?, 
 %    Phys. Lett. A, 330, 2004.
+%    Hirata, Y., Horai, S., Aihara, K.: 
+%    Reproduction of distance matrices from recurrence plots 
+%    and its applications, 
+%    Eur. Phys. J. ST, 164, 2008.
 
 % Copyright (c) 2008-
 % Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
@@ -29,6 +41,10 @@ function varargout=recons(varargin)
 % $Revision$
 %
 % $Log$
+% Revision 5.4  2016/03/03 14:57:40  marwan
+% updated input/output check (because Mathworks is removing downwards compatibility)
+% bug in crqad_big fixed (calculation of xcf).
+%
 % Revision 5.3  2009/03/24 08:34:10  marwan
 % copyright address changed
 %
@@ -50,7 +66,7 @@ function varargout=recons(varargin)
 
 
 errcode = 0; %xout = [];
-narginchk(1,2)
+narginchk(1,3)
 nargoutchk(0,1)
 
 if isnumeric(varargin{1})==1 		% read RP
@@ -65,119 +81,172 @@ N = size(X); maxN = prod(N)+1;
 if N(1) ~= N(2)
     error('Input must be square matrix.')
 end
+N = N(1);
 
 try
+         
 % final distribution        
 errcode = 1;
-i=(-5:10/(N(1)-1):5)';
-fs=cumsum(10*(1/(sqrt(2*pi))).*exp(-(i).^2./2));
-P=interp1(fs,i,[1:length(fs)],'nearest')';
+i = (-5:10/(N-1):5)';
+fs = cumsum(10*(1/(sqrt(2*pi))).*exp(-(i).^2./2));
+P = interp1(fs,i,[1:length(fs)],'nearest')';
+
+% check for input arguments of type char
+method = 1; % default reconstruction method by Thiel
+i_char = find(cellfun('isclass',varargin,'char'));
+for i = i_char
+    if strcmpi(varargin{i}(1:3), 'thi')
+         method = 1;
+         disp('Use reconstruction by Thiel.')
+    end
+    if strcmpi(varargin{i}(1:3), 'hir')
+         method = 2;
+         disp('Use reconstruction by Hirata')
+    end
+    if strcmpi(varargin{i}(1:3), 'uni')
+         P = (0:1/(N-1):1)';
+         disp('Use uniform distribution.')
+    end
+end
 
 if nargin > 1 & isnumeric(varargin{2})==1
     P = varargin{2};
-elseif nargin > 1 & ischar(varargin{2})==1
-    if strcmpi(varargin{2}(1:3), 'uni')
-        P = (0:1/(N(1)-1):1)';
-        disp('uni')
-    end
 end
 P(isnan(P)) = []; % remove NaN from distribution
 
+h_waitbar = waitbar(0,'Reconstruct trajectory'); drawnow
 
-h_waitbar = waitbar(0,'Remove double columns'); drawnow
+% number of non-neighbours
 
-% remove double rows
-errcode = 2;
-[X2 uniI J]=unique(X,'rows');
-uniI = sort(uniI);
+if method == 2;
+%%%%% Method by Hirata
+    % calculate weights
+    waitbar(0, h_waitbar, 'Calculate weights'); drawnow
+    W = 10000*ones(N,N);
+    errcode = 3;
+    for i = 1:N, waitbar(i/N)
+       neighb = find(X(i,:)); % neighbours of i
+       % number of intersections
+       Ninter = sum(X(neighb,:).*repmat(X(i,:),length(neighb),1),2); 
+       % number of union
+       Nunion = sum((X(neighb,:)+repmat(X(i,:),length(neighb),1))>0,2); 
+       W(i, neighb) = (1 - Ninter./Nunion)';
+    end
 
-% this are the removed rows
-rem = setdiff(1:length(X2),uniI); 
-%X3(rem,:) = X(rem,:);
+    % shortest path
+    dists = zeros(N,N);
+    waitbar(0, h_waitbar, 'Calculate shortest paths'); drawnow
+    for i = 1:N, waitbar(i/N)
+        sds = dijkstra(W,i);
+        dists(i,:) = sds;
+    end
 
-% this are the corresponding columns
-for i = 1:length(rem), waitbar(i/length(rem))
-%    [dummy, dblI(i), a2] = intersect(X, X(rem(i),:), 'rows');
-    X2 = X;
-    X2(rem,:) = -1;
-    dblI(i) = find(all((X2 == repmat(X(rem(i),:), size(X2,2), 1))'));
-end
 
-% number of non-neighbours
-errcode = 3;
-
-waitbar(0, h_waitbar, 'Search neighbours'); drawnow
-NN = maxN*ones(N);
-%N = zeros(size(X1));
-for i = 1:size(X,1), waitbar(i/size(X,1))
-    if ~ismember(i,rem)
-        % indices of RPs at i
-        I = find(X(i,:));
-        I(ismember(I,rem)) = [];
-
-        % number of non-neighbours
-        n = sum((X(I,:) - repmat(X(i,:),length(I),1)) == 1,2);
-        NN(i,I) = n';
+    r = cmdscale(dists); 
+    xneu = r(:,1);
+
+
+else
+
+%%%%% Method by Thiel
+
+    % remove double rows
+    errcode = 2;
+    [X2 uniI J]=unique(X,'rows');
+    uniI = sort(uniI);
+
+    % this are the removed rows
+    rem = setdiff(1:length(X2),uniI); 
+    %X3(rem,:) = X(rem,:);
+
+    % this are the corresponding columns
+    for i = 1:length(rem), waitbar(i/length(rem))
+    %    [dummy, dblI(i), a2] = intersect(X, X(rem(i),:), 'rows');
+        X2 = X;
+        X2(rem,:) = -1;
+        dblI(i) = find(all((X2 == repmat(X(rem(i),:), size(X2,2), 1))'));
     end
-    
-end
-errcode = 4;
-
-% remove entries on the main diagonal in N
-i(1:size(NN,1)) = 1; NN(find(diag(i))) = maxN;
-
-% determine the both indices where N == 0 for all i
-iNaN = NN == maxN; i = find(iNaN); NN2 = NN; NN2(i) = zeros(size(i)); % remove NaNs
-n = sum(NN2,1);
-borderI=find(n==0);
-borderI(ismember(borderI,rem)) = [];
-
-k = borderI(1);
-r = k; % start point in the rank order vector
-
-errcode = 5;
-waitbar(0, h_waitbar, 'Reconstructing time series'); drawnow
-while min(NN(:)) < maxN
-    % look for the minimal N(k,:)
-    waitbar(length(r)/N(1))
-    [dummy kneu] = min(NN(k,:));
-    if ~isempty(dummy) kneu = find(NN(k,:) == dummy); else break, end
-    if length(kneu) > 1 % if a set of minimal N exist
-        [dummy kneu_index] = min(NN(kneu,k));
-        kneu=kneu(kneu_index);
+
+    % number of non-neighbours
+    errcode = 3;
+
+    waitbar(0, h_waitbar, 'Search neighbours'); drawnow
+    NN = maxN*ones(N);
+    %N = zeros(size(X1));
+    for i = 1:size(X,1), waitbar(i/size(X,1))
+        if ~ismember(i,rem)
+            % indices of RPs at i
+            I = find(X(i,:));
+            I(ismember(I,rem)) = [];
+
+            % number of non-neighbours
+            n = sum((X(I,:) - repmat(X(i,:),length(I),1)) == 1,2);
+            NN(i,I) = n';
+        end
+
     end
-    NN(:,k) = maxN;
-    kold = k;
-    k=kneu;
-    r = [r; k]; % add the new found index to the rank order vector
-    if length(r) > N(1) + 1
-        delete(h_waitbar)
-        disp(['Critical abort. Could not find enough corresponding neigbours.',10,'Perhaps the recurrence threshold is too small.'])
-        if nargout
-            varargout{1} = NaN;
+    errcode = 4;
+
+    % remove entries on the main diagonal in N
+    i(1:size(NN,1)) = 1; NN(find(diag(i))) = maxN;
+
+    % determine the both indices where N == 0 for all i
+    iNaN = NN == maxN; i = find(iNaN); NN2 = NN; NN2(i) = zeros(size(i)); % remove NaNs
+    n = sum(NN2,1);
+    borderI=find(n==0);
+    borderI(ismember(borderI,rem)) = [];
+
+    k = borderI(1);
+    r = k; % start point in the rank order vector
+
+    errcode = 5;
+    waitbar(0, h_waitbar, 'Reconstructing time series'); drawnow
+    while min(NN(:)) < maxN
+        % look for the minimal N(k,:)
+        waitbar(length(r)/N(1))
+        [dummy kneu] = min(NN(k,:));
+        if ~isempty(dummy) kneu = find(NN(k,:) == dummy); else break, end
+        if length(kneu) > 1 % if a set of minimal N exist
+            [dummy kneu_index] = min(NN(kneu,k));
+            kneu=kneu(kneu_index);
+        end
+        NN(:,k) = maxN;
+        kold = k;
+        k=kneu;
+        r = [r; k]; % add the new found index to the rank order vector
+        if length(r) > N(1) + 1
+            delete(h_waitbar)
+            disp(['Critical abort. Could not find enough corresponding neigbours.',10,'Perhaps the recurrence threshold is too small.'])
+            if nargout
+                varargout{1} = NaN;
+            end
+            return
         end
-        return
+
     end
-    
-end
 
-errcode = 6;
-% adjust length of the distribution
-P = interp1((1:length(P))/length(P), P, (1:length(r))/length(r),'cubic');
 
-% assign the distribution to the rank order series
-waitbar(0, h_waitbar, 'Assigning distribution'); drawnow
-clear xneu
-for i = 1:length(r);
-    xneu(r(i)) = P(i);
-    waitbar(i/length(r))
-end
 
-errcode = 7;
 
-% add the removed (doubled) elements
-for i =1:length(rem)
-    xneu(rem(i)) = xneu(dblI(i));
+    errcode = 6;
+    % adjust length of the distribution
+    P = interp1((1:length(P))/length(P), P, (1:length(r))/length(r),'cubic');
+
+    % assign the distribution to the rank order series
+    waitbar(0, h_waitbar, 'Assigning distribution'); drawnow
+    clear xneu
+    for i = 1:length(r);
+        xneu(r(i)) = P(i);
+        waitbar(i/length(r))
+    end
+
+    errcode = 7;
+
+    % add the removed (doubled) elements
+    for i =1:length(rem)
+        xneu(rem(i)) = xneu(dblI(i));
+    end
+
 end
 delete(h_waitbar)
 
@@ -248,3 +317,53 @@ catch
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function sds = dijkstra(dm, s)
+% DIJSKSTRA    Shortest path for weighted graph.
+%    SDS = DIJKSTRA(DM, S) calculates the shortest paths
+%    from vertex S to all other vertices of the graph DM.
+
+%DMAX = length(dm) + 1;
+DMAX = 100001;
+
+n = size(dm,1);
+visited = zeros(1,n);
+
+% allocate 1 at the starting node and DMAX for the rest.
+
+visited = zeros(n,1);
+sds = DMAX * ones(n,1);
+sds(s) = 0;
+
+c = n-1;
+
+while(c > 0)
+  % find the shortest path among not-visited nodes.
+  j = 1;
+  while (visited(j) == 1)
+    j = j+1;
+  end
+  dmin = sds(j);
+  k = j;
+  for i=j:n
+    if(visited(i) == 0)
+      if(dmin > sds(i))
+	    dmin = sds(i);
+	    k = i;
+	  end
+    end
+  end
+
+  visited(k) = 1;
+  c = c - 1;
+%    [c k]
+
+  % find nodes connected to k and increment
+  for i=1:n
+    if((visited(i) == 0)&(dm(k,i) >= 0))
+      if (sds(i) > sds(k) + dm(k,i))
+	    sds(i) = sds(k) + dm(k,i);
+	  end
+    end
+  end
+end
+