-
Notifications
You must be signed in to change notification settings - Fork 6
/
rp.m
202 lines (179 loc) · 7.12 KB
/
rp.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
function [r,d] = rp(varargin)
%RP Calculates a recurrence plot
% R=RP(X,E,THRESH,NORM,ALG) calculates the recurrence plot R
% from an embedding vector X and using the threshold E.
% X is a N-by-M matrix corresponding to N time points
% and M embedding dimensions.
%
% [R,D]=RP(...) outputs the recurrence plot R and the
% underlying distance matrix D.
%
% Optional arguments:
% NORM - is a string setting the norm for distance
% calculation in phasespace. Can be 'euc'
% for euclidian norm (default) or 'max'
% for maximum norm.
% ALG - is a string specifying the algorithm of
% calculating the distance matrix. Can be
% 'loops', 'vector' (default), or
% 'matlabvector'.
% THRESH - is a string that specifies how the threshold
% epsilon will be calculated. With 'fix' (default)
% the RP is computed with a fixed threshold
% epsilon specified by the input parameter E.
% With 'var' the RP is computed with a fixed
% threshold epsilon, which corresponds to the
% lower 'E'-quantile (specified by E) of the
% distance distribution of all points in
% phasespace. With 'fan' the RP is computed with
% a variable threshold resulting in a fixed amount
% of nearest neighbours in phasespace, specified
% by the fraction E of recurrence points.
%
% Reference:
% Marwan, N., Romano, M. C., Thiel, M., Kurths, J. (2007).
% Recurrence plots for the analysis of complex systems.
% Physics Reports, 438, 237-329.
% Kraemer, K. H., Donner, R. V., Heitzig, J., & Marwan, N.
% (2018). Recurrence threshold selection for obtaining robust
% recurrence characteristics in different embedding dimensions.
% Chaos, 28, 085720.
%
% Example:
% N = 300; % length of time series
% x = .9*sin((1:N)*2*pi/70); % exemplary time series
% xVec = embed(x,2,17); % embed into 2 dimensions using delay 17
% R = rp(xVec,.1,'fix','max'); % calculate RP using maximum norm and fixed threshold
% imagesc(R)
% Copyright (c) 2016-2019
% Potsdam Institute for Climate Impact Research, Germany
% Institute of Geosciences, University of Potsdam, Germany
% Norbert Marwan, K. Hauke Kraemer
% http://www.pik-potsdam.de
%
% This program is free software: you can redistribute it and/or modify it under the terms of the
% GNU Affero General Public License as published by the Free Software Foundation, either
% version 3 of the License, or (at your option) any later version.
% This program is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more
% details.
% You should have received a copy of the GNU Affero General Public License along with this
% program. If not, see <http://www.gnu.org/licenses/>.
%% check input
narginchk(1,5)
nargoutchk(0,2)
%% set default values for input parameters
algoLib={'loops','vector','matlabvector'}; % the possible algorithms
methLib={'euc','max'}; % the possible norms
thresLib={'fix','var','fan'}; % the possible ways of threshold computation
algorithm = 'vector';
meth = 'euc'; % norm
thres = 'fix'; % threshold algorithm
e = 1; % recurrence threshold
%% get input arguments
% set the algorithm of computing the RP
if nargin > 4
if isa(varargin{5},'char') & ismember(varargin{5},algoLib)
algorithm = varargin{5};
else
warning(['Specified algorithm should be one of the following possible values:',...
10,sprintf('''%s'' ',algoLib{:})])
end
end
% set norm
if nargin > 3
if isa(varargin{4},'char') & ismember(varargin{4},methLib)
meth = lower(varargin{4});
else
warning(['Specified norm should be one of the following possible values:',...
10,sprintf('''%s'' ',methLib{:})])
end
end
% set method for threshold calculation
if nargin > 2
if isa(varargin{3},'char') & ismember(varargin{3},thresLib)
thres = lower(varargin{3});
else
warning(['Specified way of calculating threshold should be one of the following possible values:',...
10,sprintf('''%s'' ',thresLib{:})])
end
end
% set threshold value
if nargin > 1
if isa(varargin{2},'double')
e = varargin{2};
else
warning('Threshold has to be numeric.')
end
end
% embedding vector
x = varargin{1};
N = size(x); % size of the embedding vector
if N(1) < N(2)
error('Embedding dimension is larger than the length of the vector. Please check!')
end
%% calculate distance matrix D
switch algorithm
% calculation using loops
case 'loops'
d = zeros(N(1),N(1)); % allocate matrix
for i = 1:N(1)
for j = 1:N(1)
switch meth % select norm
case 'euc'
% euclidean distance between two phase space vectors
s = (x(i,:) - x(j,:)).^2;
d(i,j) = sqrt(sum(s));
otherwise
% max-norm distance between two phase space vectors
s = abs(x(i,:) - x(j,:));
d(i,j) = max(s);
end
end
end
% calculation using full vectorisation
case 'vector'
% repeat vector entries to have all pair-wise combinations
x1 = repmat(x,N(1),1);
x2 = reshape(repmat(reshape(x,N(1)*N(2),1),1,N(1))',N(1)*N(1),N(2));
switch meth
case 'euc'
% euclidean distance between two phase space vectors
s = (x1 - x2).^2;
d = sqrt(sum(s,2));
otherwise
% max-norm distance between two phase space vectors
s = abs(x1 - x2);
d = max(s,[],2);
end
d = reshape(d,N(1), N(1)); % reshape to have a square matrix
% calculation using MATLAB's pdist function
case 'matlabvector'
switch meth
case 'euc'
% euclidean distance between two phase space vectors
d = pdist2(x,x);
otherwise
% max-norm distance between two phase space vectors
d = pdist2(x,x,'chebychev');
end
end
%% apply threshold to get the final RP
% calculate the threshold
if strcmp(thres,'fix')
% simply apply threshold to the distance matrix
r = double(d < e);
elseif strcmp(thres,'var')
% threshold based on lower (e*100)%-quantile of distance-distribution
% get the lower e%-quantile
epsilon = quantile(d(:),e);
r = double(d < epsilon);
elseif strcmp(thres,'fan')
% variable threshold for each point in order to get fixed
% number of nearest neighbours
q = quantile(d,e); % distance that corresponds to the fraction e of rec. points per column
thresholds = repmat(q,N(1),1); % q has to be applied for each row in d
% apply individual thresholds
r = double(d<thresholds);
end