Skip to content

Commit 377cf1c

Browse files
brendanjohnharrisbrendanjohnharris
committedJan 9, 2024
Start testing RAD on other systems
·
v0.2.1v0.2.0
1 parent 7e6c344 commit 377cf1c

File tree

11 files changed

+623
-466
lines changed

11 files changed

+623
-466
lines changed
 

‎.vscode/settings.json

+3
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
{
2+
"julia.environmentPath": "/home/brendan/OneDrive/PhD/Code/Features/Criticality"
3+
}

‎Additional_functions/reWriter/reWriter.m

+19-2
Original file line numberDiff line numberDiff line change
@@ -1,50 +1,67 @@
11
classdef reWriter < handle
2-
properties
2+
3+
properties
34
LastBytes = 0;
45
Active = 1;
56
thatWasMe = 0;
67
ID = [];
78
end
9+
810
methods (Static)
11+
912
function out = setgetWhichOne(x)
1013
% So, this is a vector
1114
% The first element is the currently active reWriter
1215
% The subsequent elements refer to created reWriters (whether
1316
% thay have been deleted or not), with a unique index
1417
% defined as the ID property of each instance.
1518
persistent whichOne;
19+
1620
if nargin
1721
whichOne = x;
1822
end
23+
1924
out = whichOne;
2025
end
26+
2127
end
28+
2229
methods
30+
2331
function obj = reWriter()
32+
2433
if isempty(reWriter.setgetWhichOne)
2534
obj.ID = 1;
2635
reWriter.setgetWhichOne([0, 1]); % You are the only instance, but you are not yet activated
2736
else
2837
obj.ID = max(reWriter.setgetWhichOne) + 1; % You can have the next availiable ID
2938
reWriter.setgetWhichOne([reWriter.setgetWhichOne, obj.ID]);
30-
end
39+
end
40+
3141
end
42+
3243
function obj = reWrite(obj, message, varargin)
3344
whichOne = reWriter.setgetWhichOne;
45+
3446
if obj.LastBytes > 0 && obj.Active && (obj.ID == whichOne(1))
3547
fprintf(repmat('\b', 1, obj.LastBytes));
3648
end
49+
3750
obj = Write(obj, message, varargin{:});
3851
end
52+
3953
function obj = Write(obj, message, varargin)
4054
Bytes = builtin('fprintf', message, varargin{:});
4155
obj.LastBytes = Bytes;
4256
whichOne = reWriter.setgetWhichOne;
4357
whichOne(1) = obj.ID; % Pick me
4458
reWriter.setgetWhichOne(whichOne);
4559
end
60+
4661
function obj = inactivate(obj)
4762
obj.Active = 0;
4863
end
64+
4965
end
66+
5067
end

‎Feature_Analysis/Plot_Data/plot_feature_vals_directly.m

+11-4
Original file line numberDiff line numberDiff line change
@@ -1,33 +1,40 @@
11
function plot_feature_vals_directly(ops, mops, input_file, cp_range, etarange)
2-
% Modification of 'predict_direction' that only plots the feature values
2+
% Modification of 'predict_direction' that only plots the feature values
33
if isstring(ops) || ischar(ops)
44
ops = SQL_add('ops', ops, 0, 0);
55
end
6+
67
if nargin < 2 || isempty(mops)
78
mops = SQL_add('mops', 'INP_mops.txt', 0, 0);
89
end
10+
911
[ops, mops] = TS_LinkOperationsWithMasters(ops, mops);
12+
1013
if ischar(input_file) || isstring(input_file)
1114
f = struct2cell(load(input_file));
1215
parameters = f{1};
1316
else
1417
parameters = input_file;
1518
end
19+
1620
if nargin > 3 && ~isempty(cp_range)
1721
parameters.cp_range = cp_range;
1822
end
23+
1924
if nargin > 4 && ~isempty(etarange)
2025
parameters.etarange = etarange;
2126
end
27+
2228
etarange = parameters.etarange;
29+
2330
for ind = 1:length(etarange)
2431
p = parameters;
2532
p.etarange = etarange(ind);
2633
time_series_data = time_series_generator('input_struct', p);
2734
feature_vals = generate_feature_vals(time_series_data, ops, mops, 0);
28-
% figure
35+
% figure
2936
plot(p.cp_range, feature_vals', '-', 'markersize', 20, 'Marker', '.')
30-
title(sprintf('\\eta = %g', etarange(ind)), 'interpreter', 'Tex')
37+
title(sprintf('\\eta = %g', etarange(ind)), 'interpreter', 'Tex')
3138
end
32-
end
3339

40+
end

‎Feature_Analysis/Plot_Data/scatterScript.m

+1
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
st.String = '';%['g) ', st.String];
2727
pdfprint('NewFeature_Scatter.pdf', '-dpdf')
2828

29+
%%
2930
st = plot_feature_vals(8009, time_series_data, 'noise', 1, [5, 25, 50, 75, 100], 1); ax = gca; ax.XTick = ax.XTick(1:2:end);
3031
st.String = '';%['g) ', st.String];
3132
pdfprint('fitSupercriticalHopfRadius_Scatter.pdf', '-dpdf')

‎Time_Series_Generation/time_series_generator/TSG_systems.m

+10
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,16 @@
8383
end
8484
end
8585
rout = abs(rout);
86+
87+
case 'supercritical_hopf_radius_(measurement)'
88+
for n = 2:numpoints-1
89+
r = r + (mu.*r - (r.^3)).*dt + 0.1.*sqrt(dt).*randn(Wl, 1); % Fixed dynamical eta of 0.1
90+
if n >= transient_cutoff && ~mod(n - transient_cutoff - 1, savestep)
91+
rout(fails, 1 + (n - transient_cutoff - 1)./savestep) = r;
92+
end
93+
end
94+
rout = abs(rout);
95+
rout = rout + eta.*sqrt(dt).*randn(size(rout)); % Measurement noise
8696

8797
case 'supercritical_hopf_radius_(strogatz)_pink'
8898
for n = 2:numpoints-1

‎Time_Series_Generation/time_series_generator/time_series_generator.m

+477-431
Large diffs are not rendered by default.

‎paper/measurement_noise/measurement_noise.m

+28
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
addpath(genpath('../'));
2+
load('./time_series_data.mat');
3+
num_hctsa = 7873; % The number of hctsa features. Any above this are custom
4+
5+
% * Do a similar scatter for RAD, but varying measurement noise between the scatters
6+
mops = SQL_add('mops', 'radMops.txt', 0, 0);
7+
ops = SQL_add('ops', 'radOps.txt', 0, 0);
8+
ops = ops(1, :);
9+
10+
input_file = time_series_data(1).Inputs;
11+
input_file.save_cp_split = 0;
12+
input_file.foldername = "";
13+
cp_range = -1:0.1:0.0;
14+
etarange = 0.01:0.1:1;
15+
16+
input_file.etarange = etarange;
17+
input_file.cp_range = cp_range;
18+
19+
time_series_data = time_series_generator('input_struct', input_file);
20+
21+
[ops, mops] = TS_LinkOperationsWithMasters(ops, mops);
22+
feature_vals = generate_feature_vals(time_series_data, ops, mops, 0);
23+
24+
f = figure();
25+
hold on
26+
plot_feature_vals_directly(ops, mops, input_file, cp_range, etarange)
27+
28+
% * Then do a summary figure that can go in the paper

‎paper/other_systems/gensystem.m

+20
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
% * Does RAD perform as well for other normal forms? Set up a simulation of other normal forms matching the subcritical hopf dataset, then compute RAD, AC, and SD on this dataset.
2+
% * I'm thinking transcritical, subcritical hopf/pitchfork, and sadle-node
3+
4+
% ? Start with subcritical hopf, should be easy. We will have to shift the window back a little to avoid the exact critical point, and reject any trajectories that cross the unstable threshold.
5+
system = 'supercritical_hopf_radius_(strogatz)';
6+
dryrun = false;
7+
8+
if ~dryrun
9+
testTimeseries(system, ['./Data/' system '/'], 0)
10+
cd(['./Data/' system '/results/'])
11+
TS_init('timeseries.mat', 'radMops.txt', 'radOps.txt', 0);
12+
TS_compute(0, [], [], [], [], 0);
13+
save_data('./time_series_data.mat', system, 'paper', 'HCTSA.mat', '../inputs.mat');
14+
group_by_noise('time_series_data.mat', 'time_series_data.mat')
15+
find_correlation('time_series_data.mat', 'Pearson', [-1, 0], 'time_series_data.mat');
16+
cd('../../../')
17+
end
18+
19+
% cd(['./Data/' system '/results/'])
20+
% plot_feature_vals(93, time_series_data, 'noise', true, [1, 25, 50, 75, 100], true)

‎test/radMops.txt

+6
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
RAD(x,true,1) CR_RAD_1_1
2+
RAD(x,false,1) CR_RAD_0_1
3+
CO_AutoCorr(x,1,'Fourier') AC_1
4+
DN_Spread(x,'std') DN_Spread_std
5+
DN_Moments(x,3) DN_Moments_3
6+
DN_Moments(x,4) DN_Moments_4

‎test/radOps.txt

+6
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
CR_RAD_1_1 CR_RAD_1_1 distribution,bifurcation,delay
2+
CR_RAD_0_1 CR_RAD_0_1 distribution,bifurcation,delay
3+
AC_1 AC_1 correlation
4+
DN_Spread_std standard_deviation distribution,spread,raw,spreaddep
5+
DN_Moments_3 DN_Moments_3 distribution,moment,shape
6+
DN_Moments_4 DN_Moments_4 distribution,moment,shape

‎test/testTimeseries.m

+42-29
Original file line numberDiff line numberDiff line change
@@ -1,38 +1,51 @@
1-
function testTimeseries(system, outfolder)
2-
% Generate time series for the main analysis; this should only take a few
3-
% minutes, and will produce an `inputs.mat` file, a `results` folder and 100
4-
% `time_series_data` subfolders (each with a `timeseries.mat` file and an `inputs_out.mat` file)
1+
function testTimeseries(system, outfolder, save_cp_split)
2+
% Generate time series for the main analysis; this should only take a few
3+
% minutes, and will produce an `inputs.mat` file, a `results` folder and 100
4+
% `time_series_data` subfolders (each with a `timeseries.mat` file and an `inputs_out.mat` file)
5+
if nargin < 3 || isempty(save_cp_split)
6+
save_cp_split = 100;
7+
end
8+
59
switch system
610
case 'supercritical_hopf_radius_(strogatz)'
7-
inputs = make_input_struct(0, 'cp_range', -1:0.01:0, 'etarange', 0.01:0.01:1,...
8-
'bifurcation_point', 0, 'initial_conditions', 0,...
9-
'foldername', fullfile(outfolder, 'results'), 'numpoints', 1000000,...
10-
'parameters', [], 'savelength', 5000,...
11-
'system_type', 'supercritical_hopf_radius_(strogatz)',...
12-
'tmax', 1000, 'T', 500, 'save_cp_split', 100, ...
13-
'randomise', 0, 'rngseed', 17062020, 'vocal', 1);
14-
11+
inputs = make_input_struct(0, 'cp_range', -1:0.01:0, 'etarange', 0.01:0.01:1, ...
12+
'bifurcation_point', 0, 'initial_conditions', 0, ...
13+
'foldername', fullfile(outfolder, 'results'), 'numpoints', 1000000, ...
14+
'parameters', [], 'savelength', 5000, ...
15+
'system_type', 'supercritical_hopf_radius_(strogatz)', ...
16+
'tmax', 1000, 'T', 500, 'save_cp_split', save_cp_split, ...
17+
'randomise', 0, 'rngseed', 17062020, 'vocal', 1);
18+
19+
case 'supercritical_hopf_radius_(measurement)'
20+
inputs = make_input_struct(0, 'cp_range', -1:0.01:0, 'etarange', 0.01:0.01:1, ...
21+
'bifurcation_point', 0, 'initial_conditions', 0, ...
22+
'foldername', fullfile(outfolder, 'results'), 'numpoints', 1000000, ...
23+
'parameters', [], 'savelength', 5000, ...
24+
'system_type', 'supercritical_hopf_radius_(strogatz)', ...
25+
'tmax', 1000, 'T', 500, 'save_cp_split', save_cp_split, ...
26+
'randomise', 0, 'rngseed', 17062020, 'vocal', 1);
27+
1528
case 'saddle_node'
16-
inputs = make_input_struct(0, 'cp_range', -1:0.01:-0.01, 'etarange', 0.0001:0.0001:0.01,...
17-
'bifurcation_point', 0, 'initial_conditions', 0,...
18-
'foldername', fullfile(outfolder, 'results'), 'numpoints', 1000000,...
19-
'parameters', [], 'savelength', 5000,...
20-
'system_type', 'saddle_node',...
21-
'tmax', 1000, 'T', 500, 'save_cp_split', 100, ...
22-
'randomise', 0, 'rngseed', 29062020, 'vocal', 1, 'criteria', 'max(rout, [], 2) < 3.*sqrt(-cp_range)', 'maxAttempts', 10);
23-
29+
inputs = make_input_struct(0, 'cp_range', -1.05:0.01:-0.05, 'etarange', 0.001:0.001:0.1, ...
30+
'bifurcation_point', 0, 'initial_conditions', 0, ...
31+
'foldername', fullfile(outfolder, 'results'), 'numpoints', 1000000, ...
32+
'parameters', [], 'savelength', 5000, ...
33+
'system_type', 'saddle_node', ...
34+
'tmax', 1000, 'T', 500, 'save_cp_split', save_cp_split, ...
35+
'randomise', 0, 'rngseed', 29062020, 'vocal', 1, 'criteria', "max(rout, [], 2)' < 3.*sqrt(-cp_range)", 'maxAttempts', 1000);
36+
2437
case 'lowHighNoise'
25-
inputs = make_input_struct(0, 'cp_range', -1:0.001:0, 'etarange', [0.0001, 1],...
26-
'bifurcation_point', 0, 'initial_conditions', 0,...
27-
'foldername', fullfile(outfolder, 'results'), 'numpoints', 1000000,...
28-
'parameters', [], 'savelength', 5000,...
29-
'system_type', 'supercritical_hopf_radius_(strogatz)',...
30-
'tmax', 1000, 'T', 500, ...
31-
'randomise', 0, 'rngseed', 30062020, 'vocal', 1);
32-
38+
inputs = make_input_struct(0, 'cp_range', -1:0.001:0, 'etarange', [0.0001, 1], ...
39+
'bifurcation_point', 0, 'initial_conditions', 0, ...
40+
'foldername', fullfile(outfolder, 'results'), 'numpoints', 1000000, ...
41+
'parameters', [], 'savelength', 5000, ...
42+
'system_type', 'supercritical_hopf_radius_(strogatz)', ...
43+
'tmax', 1000, 'T', 500, ...
44+
'randomise', 0, 'rngseed', 30062020, 'vocal', 1);
45+
3346
end
3447

3548
time_series_generator('input_struct', inputs);
3649

3750
save(fullfile(outfolder, 'inputs.mat'), 'inputs')
38-
end
51+
end

0 commit comments

Comments
 (0)
Please sign in to comment.