Skip to content

Commit cb2c528

Browse files
pucicupucicu
committed
iso-directional RP and perpendicular RP added
1 parent eda7eee commit cb2c528

File tree

5 files changed

+251
-0
lines changed

5 files changed

+251
-0
lines changed

‎example_isoRP.m

+45
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
%% MATLAB Example for calculating perpendicular and iso-directional RP
2+
3+
%% calculate example system Lorenz
4+
[t x] = ode45('lorenz',[0 200],rand(1,3));
5+
%x = resample(x(10300:11000,:),4,1);
6+
x = x(10300:11000,:);
7+
8+
% perpendicular RP
9+
[R1, SP, R0] = rp_perp(x,5,.5);
10+
11+
% iso-directional RP
12+
R2 = rp_iso(x,5,.01);
13+
14+
subplot(131)
15+
imagesc(R0)
16+
title(sprintf('RR=%.3f',mean(R0(:))))
17+
18+
subplot(132)
19+
imagesc(R1)
20+
title(sprintf('RR=%.3f',mean(R1(:))))
21+
22+
subplot(133)
23+
imagesc(R2)
24+
title(sprintf('RR=%.3f',mean(R2(:))))
25+
26+
27+
%% calculate example sine
28+
xSin = sin(linspace(0,2*pi*7,1000));
29+
x = embed(xSin,2,36);
30+
31+
% perpendicular RP (should be empty)
32+
[R1, SP, R0] = rp_perp(x,.5,.1);
33+
34+
% iso-directional RP
35+
R2 = rp_iso(x,.5,.001);
36+
37+
38+
subplot(131)
39+
imagesc(R0)
40+
41+
subplot(132)
42+
imagesc(R1)
43+
44+
subplot(133)
45+
imagesc(R2)

‎icon.png

47.9 KB
Loading

‎lorenz.m

+22
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
function dy = lorenz(t,y,p1,p2)
2+
% LORENZ ODE for the Lorenz system
3+
% Usage with ode45:
4+
% [t x] = ode45('lorenz',[0 200],rand(1,3));
5+
% Changing default parameter r
6+
% r = 210;
7+
% [t x] = ode45('lorenz',[0 200],rand(1,3),[],r);
8+
9+
% parameter r
10+
if nargin<4
11+
r=28;
12+
else
13+
r=p2;
14+
end
15+
16+
dy = zeros(3,1);
17+
18+
% ODE of Lorenz 63
19+
dy(1) = 10 * (y(2) - y(1));
20+
dy(2) = r*y(1) - y(2) - y(1)*y(3);
21+
dy(3) = y(1)*y(2) - (8/3)*y(3);
22+

‎rp_iso.m

+92
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
function [iso_r,sp,r] = rp_iso(varargin)
2+
% RP_ISO Calculates the isodirectional recurrence plot
3+
% R=RP_ISO(X,E,W) calculates the isodirectional recurrence plot R
4+
% from an embedding vector X and using the threshold E for the
5+
% vector distances and threshold W for the angle to be
6+
% considered as isodirectional.
7+
%
8+
% R=RP_ISO(X,E,W,TAU) estimate tangential vector using time delay TAU.
9+
10+
11+
%% check input
12+
narginchk(1,5)
13+
nargoutchk(0,3)
14+
15+
%% set default values for input parameters
16+
e = 1; % recurrence threshold
17+
w = .001; % angle threshold
18+
tau = 1; % distance between points at phase space vector for constructing the tangential trajectory
19+
20+
%% get input arguments
21+
% embedding vector
22+
x = varargin{1};
23+
N = size(x); % size of the embedding vector
24+
if N(1) < N(2)
25+
error('Embedding dimension is larger than the length of the vector. Please check!')
26+
end
27+
28+
% set threshold value
29+
if nargin > 1
30+
if isa(varargin{2},'double')
31+
e = varargin{2};
32+
else
33+
warning('Threshold has to be numeric.')
34+
end
35+
end
36+
37+
% set angle threshold value
38+
if nargin > 2
39+
if isa(varargin{3},'double')
40+
w = varargin{3};
41+
else
42+
warning('Threshold has to be numeric.')
43+
end
44+
end
45+
46+
% set time delay
47+
if nargin > 3
48+
if isa(varargin{4},'double')
49+
tau = varargin{4};
50+
else
51+
warning('Time delay has to be numeric.')
52+
end
53+
end
54+
55+
56+
%% calculate distance matrix D and RP
57+
% distance matrix using MATLAB's pdist function
58+
d = squareform(pdist(x));
59+
60+
% apply threshold and get the RP
61+
r = d < e;
62+
63+
64+
%% estimate the tangential vector
65+
% two estimation variants:
66+
% (1) use pre- and successor points
67+
% (2) use reference point and another point in the future, tau time steps ahead
68+
if 0
69+
% use pre- and successor points
70+
tangential_vec = zeros(size(x));
71+
tangential_vec(2:end-1,:) = x(3:end,:) - x(1:end-2,:);
72+
else
73+
% use some delay (due to Horai et al, 2002)
74+
tangential_vec = zeros(size(x));
75+
tangential_vec(1:end-tau,:) = x(1+tau:end,:) - x(1:end-tau,:);
76+
end
77+
78+
%% calculate dot product between the tangential vectors, if ~1 then parallel
79+
sp = zeros(N(1),N(1));
80+
for i = 1:N
81+
% dot product
82+
sp(i,:) = dot(tangential_vec, repmat(tangential_vec(i,:),N(1),1),2)';
83+
% normalize
84+
sp(i,:) = sp(i,:) ./ (vecnorm(tangential_vec,2,2) .* vecnorm(repmat(tangential_vec(i,:),N(1),1),2,2))';
85+
end
86+
87+
% apply threshold to dot product matrix to check whether parallel (when dot product ~ 1)
88+
iso_r = r .* (abs(sp-1) < w);
89+
90+
91+
92+

‎rp_perp.m

+92
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
function [perp_r,sp,r] = rp_perp(varargin)
2+
% RP_PERP Calculates the perpendicular recurrence plot
3+
% R=RP_PERP(X,E,W) calculates the perpendicular recurrence plot R
4+
% from an embedding vector X and using the threshold E for the
5+
% vector distances and threshold W for the angle to be
6+
% considered as perpendicular.
7+
%
8+
% R=RP_PERP(X,E,W,TAU) estimate tangential vector using time delay TAU.
9+
10+
%% check input
11+
narginchk(1,5)
12+
nargoutchk(0,3)
13+
14+
%% set default values for input parameters
15+
e = 1; % recurrence threshold
16+
w = .001; % angle threshold
17+
tau = 1; % distance between points at phase space vector for constructing the tangential trajectory
18+
19+
%% get input arguments
20+
% embedding vector
21+
x = varargin{1};
22+
N = size(x); % size of the embedding vector
23+
if N(1) < N(2)
24+
error('Embedding dimension is larger than the length of the vector. Please check!')
25+
end
26+
27+
% set threshold value
28+
if nargin > 1
29+
if isa(varargin{2},'double')
30+
e = varargin{2};
31+
else
32+
warning('Threshold has to be numeric.')
33+
end
34+
end
35+
36+
% set angle threshold value
37+
if nargin > 2
38+
if isa(varargin{3},'double')
39+
w = varargin{3};
40+
else
41+
warning('Threshold has to be numeric.')
42+
end
43+
end
44+
45+
% set time delay
46+
if nargin > 3
47+
if isa(varargin{4},'double')
48+
tau = varargin{4};
49+
else
50+
warning('Time delay has to be numeric.')
51+
end
52+
end
53+
54+
55+
%% calculate distance matrix D and RP
56+
% distance matrix using MATLAB's pdist function
57+
d = squareform(pdist(x));
58+
59+
% apply threshold and get the RP
60+
r = d < e;
61+
62+
%% estimate the tangential vector
63+
% two estimation variants:
64+
% (1) use pre- and successor points
65+
% (2) use reference point and another point in the future, tau time steps ahead
66+
if 1
67+
% use pre- and successor points
68+
tangential_vec = zeros(size(x));
69+
tangential_vec(2:end-1,:) = x(3:end,:) - x(1:end-2,:);
70+
else
71+
% use some delay (due to Horai et al, 2002)
72+
tangential_vec = zeros(size(x));
73+
tangential_vec(1:end-tau,:) = x(1+tau:end,:) - x(1:end-tau,:);
74+
end
75+
76+
%% calculate dot product
77+
sp = zeros(N(1),N(1));
78+
for i = 1:N
79+
% create distance vector between reference point and potential neighbours
80+
dist_vect = x - repmat(x(i,:),N(1),1);
81+
% dot product
82+
sp(i,:) = abs(dot(dist_vect, repmat(tangential_vec(i,:),N(1),1),2))';
83+
% normalize
84+
sp(i,:) = sp(i,:) ./ (vecnorm(dist_vect,2,2) .* vecnorm(repmat(tangential_vec(i,:),N(1),1),2,2))';
85+
end
86+
87+
% apply threshold to dot product matrix to check whether perpendicular
88+
perp_r = r .* (sp < w);
89+
90+
91+
92+

0 commit comments

Comments
 (0)