diff --git a/AutomaticQC/SpikeClassifiers/SavGol.m b/AutomaticQC/SpikeClassifiers/SavGol.m new file mode 100644 index 000000000..32ec62cbc --- /dev/null +++ b/AutomaticQC/SpikeClassifiers/SavGol.m @@ -0,0 +1,42 @@ +function g = SavGol (f, nl, nr, M) + +% SAVGOL SavGol smoothes the data in the vector f by means of a +% Savitzky-Golay smoothing filter. +% +% Input: f : noisy data +% nl: number of points to the left of the reference point +% nr: number of points to the right of the reference point +% M : degree of the least squares polynomial +% +% Output: g: smoothed data +% +% W. H. Press and S. A. Teukolsky, +% Savitzky-Golay Smoothing Filters, +% Computers in Physics, 4 (1990), pp. 669-672. + +% matrix A +A = ones (nl+nr+1, M+1); +for j = M:-1:1, + A (:, j) = [-nl:nr]' .* A (:, j+1); +end + +% filter coefficients c +[Q, R] = qr (A); +c = Q (:, M+1) / R (M+1, M+1); + +% smoothing of the noisy data +% Note that there are two equivalent ways to apply the Savitzky-Golay +% filter to the vector f. In the first case we use a for-loop whereas +% in the second case we use the faster built-in function filter. +% +% g = f; +% n = size (f); +% for i = 1+nl:n-nr, +% g (i) = c' * f (i-nl:i+nr); +% end +% +n = length (f); +g = filter (c (nl+nr+1:-1:1), 1, f); +g (1:nl) = f (1:nl); +g (nl+1:n-nr) = g (nl+nr+1:n); +g (n-nr+1:n) = f (n-nr+1:n); diff --git a/AutomaticQC/SpikeClassifiers/imosSpikeClassifierHampel.m b/AutomaticQC/SpikeClassifiers/imosSpikeClassifierHampel.m new file mode 100644 index 000000000..0eea233e7 --- /dev/null +++ b/AutomaticQC/SpikeClassifiers/imosSpikeClassifierHampel.m @@ -0,0 +1,57 @@ +function [spikes, fsignal] = imosSpikeClassifierHampel(signal, window, stdfactor) +% function spikes, fsignal = imosSpikeClassifierHampel(signal, window, stdfactor) +% +% A wrapper to Hampel, a median windowed filter, for both +% Burst and Equally spaced time series. +% +% The hampel method is a robust parametric denoising filter that uses +% the MAD deviation (scaled) from the median in a pre-defined +% indexed centred window. +% +% See also hampel.m. +% +% Inputs: +% +% signal - A 1-d signal. +% window - A window index range. If signal is in burst mode, this is ignored. +% stdfactor - A multiplying scale for the standard deviation. +% +% Outputs: +% +% spikes - An array with spikes indexes. +% fsignal - A filtered signal where the spikes are subsituted +% by the median of the window. +% +% Example: +% +% +% author: hugo.oliveira@utas.edu.au +% + +% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . +% + +if nargin < 3 + window = 3; + stdfactor = 10; +elseif nargin < 2 + stdfactor = 10; +end + +[fsignal, bind, ~, ~] = hampel(signal, window, stdfactor); +spikes = find(bind); +end diff --git a/AutomaticQC/SpikeClassifiers/imosSpikeClassifierNonBurstSavGolOTSU.m b/AutomaticQC/SpikeClassifiers/imosSpikeClassifierNonBurstSavGolOTSU.m new file mode 100644 index 000000000..102946e6c --- /dev/null +++ b/AutomaticQC/SpikeClassifiers/imosSpikeClassifierNonBurstSavGolOTSU.m @@ -0,0 +1,111 @@ +function [spikes] = imosSpikeClassifierNonBurstSavGolOTSU(signal, window, pdeg, nbins, oscale) +% function spikes = imosSpikeClassifierNonBurstSavGolOTSU(signal, window, pdeg, nbins, oscale) +% +% Detect spikes in signal by using the Otsu threshold method +% applied to the noise estimated by a Savitzy Golay Filter applied on the signal. +% +% The routine uses the SavGol.m implementation provided by +% Ken Ridgway. +% +% Inputs: +% +% signal - the signal. +% window - the SavGol window argument (odd integer). +% pdeg - the SavGol lsq polynomial order (integer). +% nbins - the number of histogram bins in the OTSU threshold method (positive integer) +% oscale - a constant scale factor to weight the OTSU threshold method (positive number) +% +% Outputs: +% +% spikes - spikes indexes. +% +% Example: +% +% t = 0:3600:86400*5; +% omega = 2*pi/86400; +% signal = 10*sin(omega*t); +% signal([20]) = max(signal).^3; +% [spikes] = imosSpikeClassifierNonBurstSavGolOTSU(signal,5,2,100,5); +% assert(spikes,20) +% +% author: Unknown +% author: Ken Ridgway +% author: hugo.oliveira@utas.edu.au [rewrite/refactoring] + +% +% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . +% + +if nargin == 1 + window = 5; %minimum filtering + pdeg = 2; %quadratic polynomial + nbins = 256; % default otsu bins + oscale = 1; %arbitrary scale used to reduce the otsu threshold +elseif nargin == 2 + pdeg = 2; %quadratic polynomial + nbins = 256; % default bins in despiking1 + oscale = 1; %arbitrary scale used to reduce the otsu threshold +elseif nargin == 3 + nbins = 256; % default bins in despiking1 + oscale = 1; %arbitrary scale used to reduce the otsu threshold +elseif nargin == 4 + oscale = 1; %arbitrary scale used to reduce the otsu threshold +end + +if window <= 3 + error('Window need to be even and larger than 3') +else + + if rem(window, 2) == 0 + window = window - 1; + end + +end + +fsignal = SavGol(signal, (window - 1) / 2, (window - 1) / 2, pdeg); +noise = signal - fsignal; +abs_noise = abs(noise); +otsu = otsu_threshold(abs_noise, nbins); +threshold = otsu / oscale; +ring_spikes = find(abs_noise > threshold | abs_noise <- threshold); +spikes = centred_indexes(ring_spikes); + +end + +function [cgind] = centred_indexes(sind) +% Return the central value of a group of +% consecutive integers in a sorted array of integers. + +s_start = sind(1); +s_end = sind(end); +cgind = sind*NaN; +c=1; +for k = 1:length(sind) - 1 + adjacent = sind(k + 1) - sind(k) <= 1; + if adjacent + continue + end + + s_end = sind(k); + cgind(c) = floor((s_end + s_start) / 2); + s_start = sind(k + 1); + c = c+1; +end +cgind(end) = floor((s_end + s_start) / 2); + +cgind = cgind(~isnan(cgind)); +end diff --git a/AutomaticQC/SpikeClassifiers/imosSpikeClassifierOTSU.m b/AutomaticQC/SpikeClassifiers/imosSpikeClassifierOTSU.m new file mode 100644 index 000000000..f76746133 --- /dev/null +++ b/AutomaticQC/SpikeClassifiers/imosSpikeClassifierOTSU.m @@ -0,0 +1,114 @@ +function [spikes,threshold] = imosSpikeClassifierOTSU(signal, nbins, oscale, centralize) +% function [spikes,threshold] = imosSpikeClassifierOTSU(signal, nbins, oscale, centralize) +% +% Detect spikes in by using the Otsu threshold method. +% +% Inputs: +% +% signal - the signal. +% nbins - the number of histogram bins in the OTSU threshold method (positive integer) +% oscale - a constant scale factor to weight the OTSU threshold method (positive number) +% centralize - a boolean to average/centralize the detected indexes (default: true) +% +% Outputs: +% +% spikes - spikes indexes. +% threshold - the threshold used for spike detection. +% +% Example: +% +% signal = [0,1,0]; +% [spikes] = imosSpikeClassifierOTSU(signal,3,1,false); +% assert(spikes==2) +% +% %centralized +% signal = [0,1,1,1,0] +% [spikes] = imosSpikeClassifierOTSU(signal,3,1,false); +% assert(spikes==[4,5,6]) +% [spikes] = imosSpikeClassifierOTSU(signal,3,1,true); +% assert(spikes==[5]) +% +% %harmonic +% omega = 2*pi/86400; +% t = 0:3600:86400*5; +% signal = 10*sin(omega*t) +% signal([10,20]) = max(signal)^3; +% [spikes] = imosSpikeClassifierOTSU(signal,100,1,true); +% assert(spikes==[10,20]) +% +% % +% omega = 2*pi/86400; +% t = 0:3600:86400*5; +% signal = 10*sin(omega*t) +% signal(20) = signal(20)*2; +% signal(21:end) = signal(21:end)+signal(20); +% signal([33,55]) = sign(signal([33,55])).*signal([33,55]).^2; +% [spikes] = imosSpikeClassifierOTSU(signal,100,1,true); +% assert(spikes==[33,55]) +% +% +% author: hugo.oliveira@utas.edu.au [rewrite/refactoring] + +% +% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . +% + +if nargin == 1 + nbins = 256; + oscale = 1; + centralize = true; +elseif nargin == 2 + oscale = 1; + centralize = true; +elseif nargin == 3 + centralize = true; +end + +otsu = otsu_threshold(signal, nbins); +threshold = otsu / oscale; +spikes = find(signal > abs(threshold) | signal < -1*abs(threshold)); +if centralize + spikes = centred_indexes(spikes); +end + +end + +function [cgind] = centred_indexes(sind) +% Return the central value of a group of +% consecutive integers in a sorted array of integers. + +s_start = sind(1); +cgind = sind * NaN; +c = 1; + +for k = 1:length(sind) - 1 + adjacent = sind(k + 1) - sind(k) <= 1; + + if adjacent + continue + end + + s_end = sind(k); + cgind(c) = floor((s_end + s_start) / 2); + s_start = sind(k + 1); + c = c + 1; +end + +s_end = sind(end); +cgind(end) = floor((s_end + s_start) / 2); +cgind = cgind(~isnan(cgind)); +end diff --git a/AutomaticQC/SpikeClassifiers/otsu_threshold.m b/AutomaticQC/SpikeClassifiers/otsu_threshold.m new file mode 100644 index 000000000..fda34d551 --- /dev/null +++ b/AutomaticQC/SpikeClassifiers/otsu_threshold.m @@ -0,0 +1,79 @@ +function [threshold] = otsu_threshold(signal, nbins) +% function threshold = otsu_threshold(signal,nbins) +% +% Compute the Otsu threshold of a one dimensional signal. +% +% The Otsu threshold is a binary quantization threshold +% that minimize (maximize) the within (between) class variance. +% +% Inputs: +% +% signal - the input signal. +% nbins - the number of bins/partitions in the histogram +% Default: 256. +% +% Outputs: +% +% threshold - The threshold that maximize +% +% Example: +% +% +% [threshold] = compute_otsu_threshold(signal,nbins) +% assert() +% +% author: hugo.oliveira@utas.edu.au +% + +% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . +% + +if nargin == 1 + nbins = 256; %based on original - grayscale levels in images +end + +nsample = length(signal); +[c, amplitude] = hist(signal, nbins); +norm_counts = c ./ nsample; + +prev_variance = 0; +threshold = 0; + +for i = 1:nbins - 1 + + if norm_counts(i) == 0 + continue + end + + %left class + a1range = 1:i; + A1 = sum(norm_counts(a1range)); % A1 = P(left_class) + A1m = sum(amplitude(a1range) .* norm_counts(a1range)) / A1; + + %right class + a2range = (i + 1):nbins; + A2 = 1 - A1; % optmized since range is [0,1] + A2m = sum(amplitude(a2range) .* norm_counts(a2range)) / A2; + + variance = A1 * A2 * (A1m - A2m)^2; %variance between class + + if variance > prev_variance + threshold = amplitude(i + 1); % t = the next histogram bin limit + prev_variance = variance; + end + +end diff --git a/AutomaticQC/imosTimeSeriesSpikeQC.m b/AutomaticQC/imosTimeSeriesSpikeQC.m new file mode 100644 index 000000000..d0adcd9ef --- /dev/null +++ b/AutomaticQC/imosTimeSeriesSpikeQC.m @@ -0,0 +1,236 @@ +function [data, flags, paramsLog] = imosTimeSeriesSpikeQC(sample_data, data, k, type, auto) +%function [data, flags, paramsLog] = imosTimeSeriesSpikeQC( sample_data, data, k, type, auto ) +% +% The top-level function to initialize Spike Tests over data. +% in and out water time. +% +% Inputs: +% sample_data - struct containing the entire data set and dimension data. +% +% data - the vector of data to check. +% +% k - Index into the sample_data.variables vector. +% +% type - dimensions/variables type to check in sample_data. +% +% auto - logical, run QC in batch mode +% +% Outputs: +% data - Same as input. +% +% flags - Vector the same size as data, with before in-water samples +% flagged. +% +% paramsLog - string containing details about params' procedure to include in QC log +% +% + +% +% Copyright (C) 2017, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% 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 General Public License for more details. + +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . +% +narginchk(4, 5); +if ~isstruct(sample_data), error('sample_data must be a struct'); end +if ~isscalar(k) ||~isnumeric(k), error('k must be a numeric scalar'); end +if ~ischar(type), error('type must be a string'); end + +global apply_to_all_spike_qc +global button_action_spike_qc + +% auto logical in input to enable running under batch processing +if nargin < 5, auto = false; end + +paramsLog = []; +flags = []; + +is_dimension_processing = strcmp(type, 'dimensions'); + +if is_dimension_processing + return +else + first_call = k == 1; + if first_call + apply_to_all_spike_qc = []; + button_action_spike_qc = 'ok'; + else + abort_requested = strcmpi(button_action_spike_qc, 'abort'); + if abort_requested + return + end + end +end + + +default_opts_file = ['AutomaticQC' filesep 'imosTimeSeriesSpikeQC.txt']; +allowed_variables = strsplit(readProperty('allowed_variables', default_opts_file), ','); +allowed_variables = cellfun(@strtrim, allowed_variables, 'UniformOutput', false); + +selected_variable = sample_data.variables{k}.name; + +not_supported = ~contains(selected_variable, allowed_variables); +if not_supported + return +end + +has_burst = ~isempty(sample_data.instrument_burst_duration) ||~isempty(sample_data.instrument_burst_interval); + +if has_burst + warning('Spike burst processing not implemented yet') + return +end + +if strcmp(readProperty('toolbox.mode'), 'profile') + warning('Spike profile processing not implemented yet') + return +end + +classifiers_visible_names = containers.Map; +classifiers_visible_names('Hampel') = load_hampel_opts(default_opts_file); +classifiers_visible_names('OTSU-Thresholding/despiking1') = load_otsu_threshold_opts(default_opts_file); +classifiers_visible_names('OTSU-Savgol/despiking2') = load_otsu_savgol_opts(default_opts_file); + +%set defaults and ask user +if isempty(apply_to_all_spike_qc) + [method_selected_in_popup,button_action_spike_qc] = init_popup(classifiers_visible_names, selected_variable); + method = classifiers_visible_names(method_selected_in_popup); +else + button_action_spike_qc = 'ok'; + method = classifiers_visible_names(apply_to_all_spike_qc); +end + +bailout = ~strcmpi(button_action_spike_qc, 'ok') && ~isempty(button_action_spike_qc); +if bailout + return +end + +paramsLog = sprintf('imosTimeSeriesSpikeQC: %s method used', func2str(method.fun)); +all_args = [data, method.args]; +spikes = method.fun(all_args{:}); + +qcSet = str2double(readProperty('toolbox.qc_set')); +badFlag = imosQCFlag('bad', qcSet, 'flag'); +flags = zeros(length(data), 1, 'int8'); +flags(spikes) = badFlag; + +end + +%ok variable can be processed - ask the user what to use. +function [method_selected_in_popup,button_action_spike_qc] = init_popup(method_dict, selected_variable) + +method_selected_in_popup = 'Hampel'; % default and first to be shown +menu_options = keys(method_dict); +global apply_to_all_spike_qc + +figure_opts = { +'Name', sprintf('imosToolbox: Spike QC test for %s', selected_variable), ... + 'Visible', 'on', ... + 'MenuBar', 'none', ... + 'Resize', 'on', ... + 'WindowStyle', 'Modal', ... + 'NumberTitle', 'on', ... + }; +figwindow = figure(figure_opts{:}); + +text_tip = uicontrol('Style', 'text', 'String', sprintf('Select the spike test for %s', selected_variable)); +check_box = uicontrol('Style', 'checkbox'); +popup = uicontrol(figwindow, ... + 'Style', 'popupmenu', ... + 'String', menu_options, ... + 'Callback', @menu_selection); + +function apply_to_all(~, ~) + if check_box.Value + apply_to_all_spike_qc = method_selected_in_popup; + else + apply_to_all_spike_qc = []; + end +end + +function menu_selection(~, ~) +sval = popup.Value; +cstr = popup.String; +method_selected_in_popup = cstr{sval}; +end + +function abortTest(~, ~) +button_action_spike_qc = 'abort'; +delete(figwindow) +end + +function skipTest(~, ~) +button_action_spike_qc = 'skip'; +delete(figwindow) +end + +function confirmTest(~, ~) +button_action_spike_qc = 'ok'; +delete(figwindow) +end + +cancelButton = uicontrol(figwindow, 'Style', 'pushbutton', 'String', 'Abort', 'Callback', @abortTest, 'Position', [0, 0.0, 0.5, 0.1]); +skipButton = uicontrol(figwindow, 'Style', 'pushbutton', 'String', 'Skip', 'Callback', @skipTest, 'Position', [0., 0.0, 0.2, 0.2]); +confirmButton = uicontrol(figwindow, 'Style', 'pushbutton', 'String', 'Ok', 'Callback', @confirmTest, 'Position', [0.5, 0.0, 0.5, 0.1]); + +state_list = {figwindow, check_box, text_tip, popup, cancelButton, skipButton, confirmButton}; +normalize_units(state_list); + +set(figwindow, 'Position', [0.3, 0.45, 0.4, 0.1]); +set(check_box, 'Position', [0.7, 0.32, 0.3, 0.3],'String','Apply to all variables','Callback',@apply_to_all); %[0.25,0.2,.3,.5]); +set(text_tip, 'Position', [0.23, 0.025, 0.3, 0.6]); %[0.25,0.2,.3,.5]); +set(popup, 'Position', [0.25, .25, .4, 0.25]); +set(cancelButton, 'Position', [0, 0, 0.2, 0.2]); +set(skipButton, 'Position', [0.2, 0., 0.6, 0.2]); +set(confirmButton, 'Position', [0.8, 0., 0.2, 0.2]); + +waitfor(figwindow); + +end + +function [s] = load_hampel_opts(file) +s = struct(); +read_func = @(x) str2func(readProperty(x, file)); +read_numeric = @(x) str2double(readProperty(x, file)); +s.fun = read_func('hampel_function'); +opts = {'hampel_window', 'hampel_stdfactor'}; +s.args = cellfun(read_numeric, opts, 'UniformOutput', false); +end + +function [s] = load_otsu_savgol_opts(file) +s = struct(); +read_func = @(x) str2func(readProperty(x, file)); +read_numeric = @(x) str2double(readProperty(x, file)); +s.fun = read_func('otsu_savgol_function'); +opts = {'otsu_savgol_window', 'otsu_savgol_pdeg', 'otsu_savgol_nbins', 'otsu_savgol_oscale'}; +s.args = cellfun(read_numeric, opts, 'UniformOutput', false); +end + +function [s] = load_otsu_threshold_opts(file) +s = struct(); +read_func = @(x) str2func(readProperty(x, file)); +read_numeric = @(x) str2double(readProperty(x, file)); +s.fun = read_func('otsu_threshold_function'); +opts = {'otsu_threshold_nbins', 'otsu_threshold_oscale', 'otsu_threshold_centralize'}; +s.args = cellfun(read_numeric, opts, 'UniformOutput', false); +s.args{end} = logical(s.args{end}); +end + +function normalize_units(items) + +for k = 1:length(items) + set(items{k}, 'Units', 'normalize'); +end + +end diff --git a/AutomaticQC/imosTimeSeriesSpikeQC.txt b/AutomaticQC/imosTimeSeriesSpikeQC.txt new file mode 100644 index 000000000..b19bc2c18 --- /dev/null +++ b/AutomaticQC/imosTimeSeriesSpikeQC.txt @@ -0,0 +1,30 @@ +% comma separated variable names/subnames +allowed_variables = CDOM, CNDC, CPHL, OXY, PAR, PSAL, SST, TEMP, TURB + +hampel_function = imosSpikeClassifierHampel +%window of data at each side +hampel_window = 4 +%multiplicative factor for the standard deviation within the window +hampel_stdfactor = 5 + +%SavGol Filter with OTSU Threshold - despiking2 + +otsu_savgol_function = imosSpikeClassifierNonBurstSavGolOTSU +%odd window filtering size +otsu_savgol_window = 5 +%polynomial order +otsu_savgol_pdeg = 2 +%number of histogram bins in otsu thresholding +otsu_savgol_nbins = 256 +%scale to reduce the otsu threshold +otsu_savgol_oscale = 1 + +%OTSU Threhsold - despiking1 + +otsu_threshold_function = imosSpikeClassifierOTSU +%number of histogram bins in otsu thresholding +otsu_threshold_nbins = 256 +%scale to reduce the otsu threshold +otsu_threshold_oscale = 1 +%centralization of detections +otsu_threshold_centralize = 1 diff --git a/Java/ddb.jar b/Java/ddb.jar index 00328843c..3a149b541 100644 Binary files a/Java/ddb.jar and b/Java/ddb.jar differ diff --git a/Util/Path/setToolboxPaths.m b/Util/Path/setToolboxPaths.m index a03213d6a..17733c1e7 100644 --- a/Util/Path/setToolboxPaths.m +++ b/Util/Path/setToolboxPaths.m @@ -35,37 +35,137 @@ function setToolboxPaths(toolbox_path) % If not, see . % -% This will add all required paths for the toolbox to run -folders = { -'AutomaticQC', ... -'Seawater/TEOS10', ... -'Seawater/TEOS10/library', ... -'Seawater/EOS80', ... -'Util', ... -'Util/CellUtils', ... -'Util/StructUtils', ... -'Util/Path', ... -'Util/Schema', ... -'Util/units', ... -'IMOS', ... -'NetCDF', ... -'Parser', ... -'Parser/GenericParser', ... -'Parser/GenericParser/InstrumentRules', ... -'Parser/GenericParser/InstrumentParsers', ... -'test', ... -'test/Parser', ... -'test/Preprocessing', ... -'test/Util', ... -'test/Util/Schema', ... -'test/Util/CellUtils', ... -'DDB', ... -'Preprocessing'}; +% Toolbox root level folders and subfolders that contains matlab functions +folders = { ... + frdir(fullfile(toolbox_path, 'AutomaticQC')), ... + 'DDB', ... + 'FlowManager', ... + 'Geomag', ... + 'Graph', ... + 'GUI', ... + 'IMOS', ... + 'Java', ... + 'Java/UCanAccess-3.0.2-bin', ... + 'Java/UCanAccess-3.0.2-bin/lib', ... + frdir(fullfile(toolbox_path, 'Java/bin')), ... + 'NetCDF', ... + frdir(fullfile(toolbox_path, 'NetCDF')), ... + 'Parser', ... + frdir(fullfile(toolbox_path, 'Parser')), ... + 'Preprocessing', ... + frdir(fullfile(toolbox_path, 'Preprocessing')), ... + 'Seawater/TEOS10', ... + 'Seawater/TEOS10/library', ... + 'Seawater/EOS80', ... + 'test', ... + frdir(fullfile(toolbox_path, 'test')), ... + 'Util', ... + frdir(fullfile(toolbox_path, 'Util')), ... + }; addpath(toolbox_path) for k = 1:length(folders) - addpath(fullfile(toolbox_path,folders{k})); + + if isstring(folders{k}) || ischar(folders{k}) + addpath(fullfile(toolbox_path, folders{k})); + elseif iscell(folders{k}) + rfolders = folders{k}; + + for kk = 1:length(rfolders) + addpath(rfolders{kk}) + end + + end + +end + +end + +%functions required to be in this scope since we +%may need to call setToolboxPaths from anywhere. + +function [folders] = subFolders(path) +% function files = Folders(path) +% +% Return fullpath of sub-folders. +% +% Inputs: +% +% path - a path string +% +% Outputs: +% +% folders - a cell with fullpath strings of the subfolders. +% +% Example: +% % assumes /dev/shm got 2 folders, "a","b" +% path = '/dev/shm'; +% [subfolders] = subFolders(path); +% assert(any(contains(sf,'a'))); +% assert(any(contains(sf,'b'))); +% +% author: hugo.oliveira@utas.edu.au +% + +dobj = dir(path); +folders = cell(0, 0); +p = 1; + +for k = 1:length(dobj) + name = dobj(k).name; + isfolder = dobj(k).isdir; + + if isfolder &&~strcmp(name, '.') &&~strcmp(name, '..') + folders{p} = fullfile(path, name); + p = p + 1; + end + +end + +end + +function [folders] = frdir(path) +% function folders = frdir(path) +% +% recursive dir - List all children folders given a path +% +% Inputs: +% +% path - a path string +% +% Outputs: +% +% `folders` - a cell with fullpath strings of all subfolders at all levels. +% +% Example: +% % assumes /dev/shm/c/d/e/f/g/h/j is a empty folder +% path = '/dev/shm' +% [allfiles,allfolders] = rdir(path); +% assert(strcmp(allfolders{1},'/dev/shm/c')); +% assert(strcmp(allfolders{2},'/dev/shm/c/d')); +% assert(strcmp(allfolders{end},'/dev/shm/c/d/e/f/g/h/j')); +% +% author: hugo.oliveira@utas.edu.au +% + +[folders] = subFolders(path); + +if isempty(folders) + return +end + +folders_to_walk = string(folders); + +while ~isempty(folders_to_walk) + f = folders_to_walk{1}; + folders_to_walk(1) = []; + [newfolders] = frdir(f); + + if ~isempty(newfolders) + folders = cat(2, folders, newfolders); + end + end end diff --git a/Util/StructUtils/combineStructFields.m b/Util/StructUtils/combineStructFields.m index 02f4178e2..a7a6af33b 100644 --- a/Util/StructUtils/combineStructFields.m +++ b/Util/StructUtils/combineStructFields.m @@ -5,7 +5,7 @@ % the different fieldnames in the argument structures. % % Inputs: - + % % any number of structures % % Outputs: @@ -13,12 +13,29 @@ % combined - a structure with all fieldnames of the arguments combined % % Example: + % % basic case + % a = struct(); b=a; % a.x = 1; % b.x = 2; % b.y = 3; % combined = combineStructFields(a,b); - % assert(combined.x == 2); - % assert(combined.y == 3); + % assert(combined.x==2); + % assert(combined.y==3); + % + % % multi structure case + % a = struct(); b=a; c=a; + % a.x = 1; + % b.x = 2; + % b(2).x = 3; + % b(3).x = 4; + % c.x = 3; + % c(2).x = 4; + % c(3).x = 5; + % combined = combineStructFields(a,b,c); + % assert(isequal(size(combined),[1,3])) + % assert(combined(1).x==3) + % assert(combined(2).x==4) + % assert(combined(3).x==5) % % author: hugo.oliveira@utas.edu.au @@ -42,15 +59,15 @@ combined = struct(); - for k = 1:nargin - s = varargin{k}; - snames = fieldnames(s); - - for kk = 1:length(snames) - name = snames{kk}; - combined.(name) = s.(name); + for n = 1:nargin + arg = varargin{n}; + snames = fieldnames(arg); + for s = 1:numel(arg) + for k = 1:length(snames) + name = snames{k}; + combined(s).(name) = arg(s).(name); + end end - end end diff --git a/Util/Time/stepsInDay.m b/Util/Time/stepsInDay.m new file mode 100644 index 000000000..dc45361d8 --- /dev/null +++ b/Util/Time/stepsInDay.m @@ -0,0 +1,55 @@ +function [steps] = stepsInDay(tunit) +% function [steps] = timeQuantisationStep(prec) +% +% How many steps are defined in day, for a given time unit. +% +% +% Inputs: +% +% tunit - string representing a time/precision unit. +% +% Outputs: +% +% steps - a floating point number +% +% Example: +% +% [steps] = stepsInDay('minute'); +% assert(isequal(steps,1440)) +% +% +% author: hugo.oliveira@utas.edu.au +% + +% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . + +switch tunit + case 'microsecond' + steps = 8.64e10; + case 'milisecond' + steps = 8.64e7; + case 'second' + steps = 8.64e4; + case 'minute' + steps = 1440; + case 'hour' + steps = 24; + case 'day' + steps = 1; + otherwise + error('%s not supported',prec) +end diff --git a/Util/Time/timeQuantisation.m b/Util/Time/timeQuantisation.m new file mode 100644 index 000000000..9279730de --- /dev/null +++ b/Util/Time/timeQuantisation.m @@ -0,0 +1,48 @@ +function [tq] = timeQuantisation(time, prec) +% function [tq] = timeQuantisation(time, prec) +% +% Quantize time by a certain precision. +% +% Inputs: +% +% time - a datenum time array - i.e units in days +% prec - the precision string +% +% Outputs: +% +% tq - quantize time +% +% Example: +% +% [tq] = timeQuantisation([1,2,3,4+0.001/86400],'microsecond'); +% assert(isequal(tq,[1,2,3,4+0.001/86400])) +% [tq] = timeQuantisation([1,2,3,4+0.001/86400],'milisecond'); +% assert(isequal(tq,[1,2,3,4+0.001/86400])) +% [tq] = timeQuantisation([1,2,3,4+0.001/86400],'second'); +% assert(isequal(tq,[1,2,3,4])) +% +% +% author: hugo.oliveira@utas.edu.au +% + +% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . +% + +step = stepsInDay(prec); +tq = uniformQuantise(time, 1 ./ step); + +end diff --git a/Util/Time/timeSamplingInfo.m b/Util/Time/timeSamplingInfo.m new file mode 100644 index 000000000..bb55da280 --- /dev/null +++ b/Util/Time/timeSamplingInfo.m @@ -0,0 +1,194 @@ +function [info] = timeSamplingInfo(time, prec) + % function [info] = timeSamplingInfo(time, prec) + % + % Report information about a datenum time array + % given a time precision, including burst Information, if any. + % + % See Examples below to quantise time and to + % allow (ignore) round ups (timescales). + % + % Inputs: + % + % time - a datenum time array + % prec - a time precision string for time interval matching. + % Default to 'microsecond', i.e. + % half-microsecond or lower sampling is quantised/floored. + % + % Outputs: + % + % info - a structure with information regarding time sampling. + % + % % sampling information + % .assumed_time_precision - the precision + % .constant_vector - if time is constant + % .number_of_sampling_steps - the number of different samplings at given precision + % .sampling_steps - the sampling periods in days + % .sampling_steps_in_seconds - ditto but in seconds + % .sampling_steps_median - the most frequent value of sampling dt. + % .sampling_steps_median_in_seconds - ditto but in seconds + % .sampling_step_mode - the most ocurring value of the sampling dt + % .sampling_step_mode_in_seconds - ditto but in seconds + % .consistent_sampling - if median == mode. + % .repeated_samples - if the time variable got repeated values + % .unique_sampling - if sampling intervals are unique + % .uniform_sampling - if sampling intervals are uniform + % .monotonic_sampling - if sampling are monotonic [incr/decr] + % .progressive_sampling - values are generally increasing + % .regressive_sampling - values are generally decreasing + % .non_uniform_sampling_indexes - left-most indexes of non-uniform sampling. + % .jump_magnitude - the time distance between the non-uniform and the next uniform or non-uniform sampling + % .jump_forward_indexes - non_uniform boolean array where positive non-uniform sampling is found + % .jump_backward_indexes - as above, but for negative shifts. + % + % % The burst sampling info. + % .is_burst_sampling - true if burst were detected + % .is_burst_sampling_exact - true if bursts are exact at `.burst_sampling_units`. + % .is_burst_sampling_repeats - true if repeated sampling at burst level, such as repeated timestamps. + % .burst_sampling_units = the precision for burst sampling - may be different from `prec`. + % .burst_sampling_interval - the burst sampling interval. + % .burst_sampling_interval_in_seconds - as above, but in seconds. + % .burst_duration - the burst duration. + % .burst_duration_in_seconds - as above, but in seconds. + % + % Example: + % % basic - ignore small drifts at requested scale + % + % info = timeSamplingInfo([1,2,3,4,5.3],'day'); + % assert(info.uniform_sampling); + % assert(all([info.unique_sampling,info.monotonic_sampling,info.progressive_sampling])); + % + % % extra reports about sampling + % + % info = timeSamplingInfo([1,2,1,3,1,4,1,5,1],'day'); + % assert(~all([info.unique_sampling,info.progressive_sampling,info.regressive_sampling,info.monotonic_sampling])); + % assert(isequal(info.sampling_steps_mode,[])); + % assert(all(info.jump_forward_indexes==[1 0 1 0 1 0 1 0])); + % assert(all(info.jump_backward_indexes==[0 1 0 1 0 1 0 1])); + % + % % ignore a different scale - microsecond. + % + % info = timeSamplingInfo([1,2,3,4+0.49/86400,5],'second'); + % assert(info.repeated_samples==0); + % assert(info.sampling_steps==1); + % + % % trigger different sampling since diff is at or above half the requested scale. + % + % info = timeSamplingInfo([1,2,3,4+0.5/86400,5],'second'); + % assert(info.repeated_samples==0); + % assert(isequal(info.sampling_steps*86400,[86400,86401])); + % + % % ignore constant offsets in clock drifting + % + % drifts = [0,0,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]/86400; + % times = [1,2,3,4,5,6,7,8,9,10,11,12,13]; + % info = timeSamplingInfo(times+drifts,'second'); + % assert(info.sampling_steps_in_seconds==86400); + % + % % report different samplings since drifts are at or above half the requested scale. + % + % drifts = [0,0.5,0.5,0.5,0.75,1,1.25,1.5,3,6,9,12,24]/86400; + % times = [1,2,3,4,5,6,7,8,9,10,11,12,13]; + % info = timeSamplingInfo(times+drifts,'second'); + % assert(isequal(info.sampling_steps_in_seconds,86400+[0,1,2,3,12])); + % + % + % author: hugo.oliveira@utas.edu.au + % + + % Copyright (C) 2019, Australian Ocean Data Network (AODN) and Integrated + % Marine Observing System (IMOS). + % + % This program is free software: you can redistribute it and/or modify + % it under the terms of the GNU General Public License as published by + % the Free Software Foundation version 3 of the License. + % + % 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 General Public License for more details. + % + % You should have received a copy of the GNU General Public License + % along with this program. + % If not, see . + % + + if nargin < 2 + prec = 'microsecond'; + end + + info = struct(); + info.assumed_time_precision = prec; + info.number_of_sampling_steps = 0; + info.sampling_steps = []; + info.sampling_steps_in_seconds = []; + info.sampling_steps_median = []; + info.sampling_steps_median_in_seconds = []; + info.sampling_steps_mode = []; + info.sampling_steps_mode_in_seconds = []; + info.consistent_sampling = []; + info.constant_vector = []; + info.repeated_samples = []; + info.unique_sampling = []; + info.uniform_sampling = []; + info.monotonic_sampling = []; + info.progressive_sampling = []; + info.regressive_sampling = []; + info.non_uniform_sampling_indexes = []; + info.jump_magnitude = []; + info.jump_forward_indexes = []; + info.jump_backward_indexes = []; + + if isempty(time) || length(time) == 1 + return + end + + info.unique_sampling = isunique(time); + + tdiff = timeQuantisation(diff(time),prec); + + info.sampling_steps = unique(tdiff); + info.number_of_sampling_steps = length(info.sampling_steps); + info.sampling_steps_in_seconds = unique(timeQuantisation(tdiff*86400,'second')); + info.constant_vector = all(info.sampling_steps == 0); + info.repeated_samples = info.sampling_steps(1) == 0; + info.uniform_sampling = info.number_of_sampling_steps == 1; + + mfdt = median(tdiff); + modt = mode(tdiff); + tdiff_sign = sign(tdiff); + + info.consistent_sampling = modt == mfdt; + info.sampling_steps_median = mfdt; + info.sampling_steps_median_in_seconds = mfdt*86400; + info.sampling_steps_mode = modt; + info.sampling_steps_mode_in_seconds = modt*86400; + info.monotonic_sampling = ~all(tdiff_sign==0) && ~any(diff(tdiff_sign)); + + mfdt_sign = mode(tdiff_sign); % use median to avoid symmetric false positives + info.progressive_sampling = mfdt_sign > 0; + info.regressive_sampling = mfdt_sign < 0; + + tdiff_jumps = tdiff ~= mfdt; % uses median to avoid symmetric false positives + info.non_uniform_sampling_indexes = find(tdiff_jumps); + jval = tdiff(tdiff_jumps); + has_jumps = ~isempty(jval); + if has_jumps + info.jump_magnitude = jval; + if length(jval) > 1 + info.jump_forward_indexes = jval > 0; + info.jump_backward_indexes = jval < 0; + else + + if jval > 0 + info.jump_forward_indexes = find(tdiff == jval); + info.jump_backward_indexes = []; + else + info.jump_forward_indexes = []; + info.jump_backward_indexes = find(tdiff == jval); + end + + end + + end + +end diff --git a/Util/Time/uniformQuantise.m b/Util/Time/uniformQuantise.m new file mode 100644 index 000000000..fce1b3d09 --- /dev/null +++ b/Util/Time/uniformQuantise.m @@ -0,0 +1,46 @@ +function [aq] = uniformQuantise(a, step) +% function [aq] = uniformQuantise(a,step) +% +% A wrapper to uniform quantisation of an array by a step. +% +% Inputs: +% +% a - an array +% step - a quantise step +% +% Outputs: +% +% aq - quantise array +% +% Example: +% %discard miliseconds +% msec = 0.001/86400; % ~1e-08 +% y = uniformQuantise(msec, 1e-7); +% assert(y==0); +% +% author: hugo.oliveira@utas.edu.au +% + +% Copyright (C) 2020, Australian Ocean Data Network (AODN) and Integrated +% Marine Observing System (IMOS). +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation version 3 of the License. +% +% 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. +% If not, see . +% +if step > 1 + error('Quantization step should be less than 1.') +end + +aq = step .* floor(a ./ step + 0.5); + +end diff --git a/Util/units/toCelsius.m b/Util/Units/toCelsius.m similarity index 100% rename from Util/units/toCelsius.m rename to Util/Units/toCelsius.m diff --git a/Util/timeSamplingInfo.m b/Util/timeSamplingInfo.m deleted file mode 100644 index a02f01ac8..000000000 --- a/Util/timeSamplingInfo.m +++ /dev/null @@ -1,144 +0,0 @@ -function [info] = timeSamplingInfo(t, prec), - % function [info] = timeSamplingInfo(t) - % - % Report information about a - % time array t, usually - % the returned value of datenum. - % - % Inputs: - % - % t - an array with numerical values - % prec - a time precision string for - % time interval matching. - % Default to 'milisecond'. - % - % Outputs: - % - % info - a structure with information regarding time sampling. - % .sampling_steps - the sampling periods - % .sampling_steps_median - the median of the sampling dt - % .unique_sampling - if sampling intervals are unique - % .uniform_sampling - if sampling intervals are uniform - % .monotonic_sampling - if sampling are monotonic [incr/decr] - % .progressive_sampling - values are generally increasing - % .regressive_sampling - values are generally decreasing - % .non_uniform_sampling_indexes - left-most indexes of non-uniform sampling. - % .jump_magnitude - the index distance between the non-uniform and the next uniform or non-uniform sampling - % .jump_forward_indexes - the indexes where positive non-uniform sampling is found - % .jump_backward_indexes - as above, but for negative sampling. - % - % Example: - % - % info = timeSamplingInfo([1,2,3,4]); - % assert(info.uniform_sampling); - % assert(all([info.unique_sampling,info.monotonic_sampling,info.progressive_sampling])); - % - % info = timeSamplingInfo([1,2,1,3,1,4,1,5,1]); - % assert(~all([info.unique_sampling,info.progressive_sampling,info.regressive_sampling,info.monotonic_sampling])); - % assert(isequal(info.sampling_steps_median,[])); - % assert(all(info.jump_forward_indexes==[1 0 1 0 1 0 1 0])); - % assert(all(info.jump_backward_indexes==[0 1 0 1 0 1 0 1])); - % - % author: hugo.oliveira@utas.edu.au - % - - % Copyright (C) 2019, Australian Ocean Data Network (AODN) and Integrated - % Marine Observing System (IMOS). - % - % This program is free software: you can redistribute it and/or modify - % it under the terms of the GNU General Public License as published by - % the Free Software Foundation version 3 of the License. - % - % 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 General Public License for more details. - % - % You should have received a copy of the GNU General Public License - % along with this program. - % If not, see . - % - if nargin < 2, - prec = 'milisecond'; - end - - switch prec - case 'day' - tequalizer = 1; - case 'hour' - tequalizer = 24; - case 'second' - tequalizer = 86400; - case 'milisecond' - tequalizer = 86400 * 1e3; - case 'microsecond' - tequalizer = 86400 * 1e6; - end - - tdiff = diff(round(t * tequalizer, 0)) / tequalizer; - - sampling_steps = unique(tdiff); - constant_vector = all(sampling_steps == 0); - repeated_samples = sampling_steps(1) == 0; - - mfdt = median(tdiff); - - tdiff_sign = sign(tdiff); - mfdt_sign = median(tdiff_sign); - - tdiff_jumps = tdiff ~= mfdt; - - info = struct(); - info.constant_vector = constant_vector; - info.sampling_steps = sampling_steps; - info.repeated_samples = repeated_samples; - - if find(mfdt == info.sampling_steps), - info.sampling_steps_median = mfdt; - else - info.sampling_steps_median = []; - end - - info.unique_sampling = isunique(t); - info.uniform_sampling = length(info.sampling_steps) == 1; - info.monotonic_sampling = ~any(diff(tdiff_sign)); - info.progressive_sampling = mfdt_sign > 0; - info.regressive_sampling = mfdt_sign < 0; - - jind = find(tdiff_jumps); - - if jind == 0 - info.non_uniform_sampling_indexes = []; - else - info.non_uniform_sampling_indexes = jind; - end - - jval = tdiff(tdiff_jumps); - - no_jumps = isempty(jval); - - if no_jumps, - info.jump_magnitude = []; - info.jump_forward_indexes = []; - info.jump_backward_indexes = []; - else - info.jump_magnitude = jval; - - if length(jval) > 1, - info.jump_forward_indexes = jval > 0; - info.jump_backward_indexes = jval < 0; - else, - - if jval > 0 - info.jump_forward_indexes = find(tdiff == jval); - info.jump_backward_indexes = []; - else - info.jump_forward_indexes = []; - info.jump_backward_indexes = find(tdiff == jval); - end - - end - - end - -end diff --git a/build.py b/build.py index 0c9732569..ad9b887ad 100755 --- a/build.py +++ b/build.py @@ -1,5 +1,5 @@ #!/usr/bin/env python3 -""" Build the imos stand-alone toolbox GUI using only the Matlab compiler. +"""Build the imos stand-alone toolbox GUI using only the Matlab compiler. This script ignores untracked and testfiles within the repository and weakly attach a string version to the binary filename. @@ -15,6 +15,7 @@ --root_path= The repository root path --dist_path= Where Compiled items should be stored --arch= One of win64,glnxa64,maci64 + """ import os @@ -22,28 +23,56 @@ from getpass import getuser from pathlib import Path from shutil import copy -from typing import List +from typing import List, Tuple from docopt import docopt from git import Repo -VALID_ARCHS = ['glnxa64', 'win64', 'maci64'] +VALID_ARCHS = ["glnxa64", "win64", "maci64"] ARCH_NAMES = { - 'glnxa64': 'Linux64', - 'maci64': 'Mac64', - 'win64': 'Win64', + "glnxa64": "Linux64", + "maci64": "Mac64", + "win64": "Win64", } -VERSION_FILE = '.standalone_canonical_version' +VERSION_FILE = ".standalone_canonical_version" +MATLAB_VERSION = "R2018b" +MATLAB_TOOLBOXES_TO_INCLUDE: List[str] = ["signal"] + + +def find_matlab_rootpath(arch_str): + """Locate the matlab installation root path.""" + + if arch_str != "win64": + fname = f"*/MATLAB/{MATLAB_VERSION}" + if arch_str == "glnxa64": + lookup_folders = [ + "/usr/local/", + "/opt", + f"/home/{os.environ['USER']}", + "/usr", + ] + elif arch_str == "maci64": + lookup_folders = ["/Applications/", "/Users"] + else: + fname = f"**\\MATLAB\\{MATLAB_VERSION}" + lookup_folders = ["C:\\Program Files\\"] + + for folder in lookup_folders: + m_rootpath = list(Path(folder).glob(fname)) + if m_rootpath: + return m_rootpath[0] + raise ValueError("Unable to find matlab root path") -def run(x: str): - proc = sp.run(x, stdout=sp.PIPE, stderr=sp.PIPE, shell=True) - return proc.stdout.decode('utf-8').strip(), proc.stderr.decode( - 'utf-8').strip() +def run(command: str): + """Run a command at shell/cmdline.""" + proc = sp.run(command, stdout=sp.PIPE, stderr=sp.PIPE, shell=True, check=True) + return proc.stdout.decode("utf-8").strip(), proc.stderr.decode("utf-8").strip() def create_java_call_sig_compile(root_path: str) -> str: - java_path = os.path.join(root_path, 'Java') + """Create the call to build java dependencies.""" + java_path = os.path.join(root_path, "Java") return f"cd {java_path} && ant install && cd {root_path}" @@ -67,148 +96,181 @@ def git_info(root_path: str) -> dict: raise TypeError(f"{root_path} is not a git repository") info = {} - info['username'] = getuser() - info['is_dirty'] = repo.is_dirty() - info['branch'] = str(repo.active_branch) - info['tag'] = str(repo.git.describe('--tags')) - info['modified_files'] = [x.a_path for x in repo.index.diff(None)] - info['staged_files'] = [x.a_path for x in repo.index.diff("HEAD")] - info['untracked_files'] = repo.untracked_files - - info['is_tag_release'] = len(info['tag'].split('-')) < 3 - info['is_master'] = info['branch'] == 'master' - info['is_clean'] = not info['is_dirty'] and not info['staged_files'] - info['is_official_release'] = info['is_tag_release'] and info['is_clean'] - - version = '' - if info['is_official_release']: - version += info['tag'] + info["username"] = getuser() + info["is_dirty"] = repo.is_dirty() + info["branch"] = str(repo.active_branch) + info["tag"] = str(repo.git.describe("--tags")) + info["modified_files"] = [x.a_path for x in repo.index.diff(None)] + info["staged_files"] = [x.a_path for x in repo.index.diff("HEAD")] + info["untracked_files"] = repo.untracked_files + + info["is_tag_release"] = len(info["tag"].split("-")) < 3 + info["is_master"] = info["branch"] == "master" + info["is_clean"] = not info["is_dirty"] and not info["staged_files"] + info["is_official_release"] = info["is_tag_release"] and info["is_clean"] + + version = "" + if info["is_official_release"]: + version += info["tag"] else: - version += info['branch'] - version += '-' + info['username'] - version += '-' + info['tag'] - if info['is_dirty']: - version += '-dirty' + version += info["branch"] + version += "-" + info["username"] + version += "-" + info["tag"] + if info["is_dirty"]: + version += "-dirty" - info['version'] = version + info["version"] = version return info def find_files(folder: str, ftype: str) -> list: - """ Find all the files within the repository that are not testfiles """ - for file in Path(folder).glob(os.path.join('**/', ftype)): + """Find all the files within the repository that are not testfiles.""" + for file in Path(folder).glob(os.path.join("**/", ftype)): filepath = file.as_posix() - if 'data/testfiles' not in filepath and 'tmpbuild.m' not in filepath: + if "data/testfiles" not in filepath and "tmpbuild.m" not in filepath: yield filepath def get_required_files(root_path: str, info: dict) -> (list, list): - """ obtain the mfiles and matfiles for binary compilation """ - allmfiles = list(find_files(root_path, '*.m')) - allmatfiles = list(find_files(root_path, '*.mat')) - if info['is_dirty']: + """obtain the mfiles and matfiles for binary compilation.""" + allmfiles = list(find_files(root_path, "*.m")) + allmatfiles = list(find_files(root_path, "*.mat")) + if info["is_dirty"]: print(f"Warning: {root_path} repository is dirty") - if info['modified_files']: + if info["modified_files"]: print(f"Warning: {root_path} got modified files not staged") - if info['staged_files']: + if info["staged_files"]: print(f"Warning: {root_path} got staged files not commited") - if info['untracked_files']: - print( - f"Warning: {root_path} got untracked files - Ignoring them all...") - for ufile in info['untracked_files']: + if info["untracked_files"]: + print(f"Warning: {root_path} got untracked files - Ignoring them all...") + for ufile in info["untracked_files"]: if ufile in allmfiles: allmfiles.remove(ufile) if ufile in allmatfiles: allmatfiles.remove(ufile) - mind = [ - ind for ind, file in enumerate(allmfiles) if 'imosToolbox.m' in file - ] + mind = [ind for ind, file in enumerate(allmfiles) if "imosToolbox.m" in file] mentry = allmfiles.pop(mind[0]) allmfiles = [mentry] + allmfiles return allmfiles, allmatfiles -def create_mcc_call_sig(mcc_path: str, root_path: str, dist_path: str, - mfiles: List[str], matfiles: List[str], - output_name: str, arch: str) -> List[str]: - standalone_flags = '-m' - verbose_flag = '-v' +def create_mcc_call_sig( + mcc_path: str, + root_path: str, + dist_path: str, + mfiles: List[str], + matfiles: List[str], + output_name: str, + arch: str, +) -> List[str]: + """Create the mcc call signature based on argument options.""" + standalone_flags = "-m" + verbose_flag = "-v" outdir_flag = f"-d {dist_path}" - clearpath_flag = '-N' + clearpath_flag = "-N" outname_flag = f"-o {output_name}" - verbose_flag = '-v' - warning_flag = '-w enable' - mfiles_str = ' '.join([f"{file}" for file in mfiles]) - matfiles_str = ' '.join([f"-a {file}" for file in matfiles]) - - mcc_arguments = ' '.join([ - standalone_flags, verbose_flag, outdir_flag, clearpath_flag, - outname_flag, verbose_flag, warning_flag, mfiles_str, matfiles_str - ]) - - if arch == 'win64': - tmpscript = os.path.join(root_path, 'tmpbuild.m') - with open(tmpscript, - 'w') as tmpfile: # overcome cmd.exe limitation of characters - tmpfile.write(f"eval('{mcc_path} {mcc_arguments}')") - external_call = f"matlab -nodesktop -nosplash -wait -r \"run('{tmpscript}');exit\"" - rm_call = f"del {tmpscript}" - return [external_call, rm_call] - else: - return [ - ' '.join([ - mcc_path, standalone_flags, verbose_flag, outdir_flag, - clearpath_flag, outname_flag, verbose_flag, warning_flag, - mfiles_str, matfiles_str - ]) + verbose_flag = "-v" + warning_flag = "-w enable" + mfiles_str = " ".join([f"{file}" for file in mfiles]) + matfiles_str = " ".join([f"-a {file}" for file in matfiles]) + + if MATLAB_TOOLBOXES_TO_INCLUDE: + matlab_root_path = find_matlab_rootpath(arch) + extra_toolbox_paths = [ + os.path.join(matlab_root_path, "toolbox", x, x) + for x in MATLAB_TOOLBOXES_TO_INCLUDE ] - -def get_args(args: dict) -> (str, str, str, str): - """process the cmdline arguments to return the root_path, dist_path, mcc binary, and architecture""" - if not args['--root_path']: + mcc_arguments = " ".join( + [ + standalone_flags, + verbose_flag, + outdir_flag, + clearpath_flag, + outname_flag, + verbose_flag, + warning_flag, + mfiles_str, + matfiles_str, + ], + ) + + if arch == "win64": + tmpscript = os.path.join(root_path, "tmpbuild.m") + # overcome cmd.exe limitation of characters by creating a tmp mfile + with open(tmpscript, "w") as tmpfile: + # overcome limitation of spaces within eval. + if MATLAB_TOOLBOXES_TO_INCLUDE: + # include other paths in the build call + extra_requirements = " ".join( + ["-I " + "'" + x + "'" for x in extra_toolbox_paths], + ) + cmdline = f'eval("{mcc_path} {mcc_arguments} {extra_requirements}")' + else: + cmdline = f'eval("{mcc_path} {mcc_arguments}")' + + tmpfile.write(cmdline) + external_call = ( + f"matlab -nodesktop -nosplash -wait -r \"run('{tmpscript}');exit\"" + ) + rm_call = f"del {tmpscript}" + rlist = [external_call, rm_call] + else: + # no need for extra quotes in linux, since we are in the shell. + if MATLAB_TOOLBOXES_TO_INCLUDE: + extra_requirements = " ".join(["-I " + x for x in extra_toolbox_paths]) + rlist = [" ".join([mcc_path, mcc_arguments, extra_requirements])] + else: + rlist = [" ".join([mcc_path, mcc_arguments])] + return rlist + + +def get_args(args: dict) -> Tuple[str, str, str, str]: + """process the cmdline arguments to return the root_path, dist_path, mcc binary, and + architecture.""" + if not args["--root_path"]: root_path = os.getcwd() else: - root_path = args['--root_path'] + root_path = args["--root_path"] - if not args['--dist_path']: - dist_path = os.path.join(root_path, 'dist') + if not args["--dist_path"]: + dist_path = os.path.join(root_path, "dist") exists = os.path.lexists(dist_path) if not exists: os.mkdir(dist_path) else: - dist_path = args['--dist_path'] + dist_path = args["--dist_path"] - if not args['--mcc_path']: - mcc_path = 'mcc' + if not args["--mcc_path"]: + mcc_path = "mcc" else: - mcc_path = args['--mcc_path'] + mcc_path = args["--mcc_path"] - if args['--arch'] in VALID_ARCHS: - ind = VALID_ARCHS.index(args['--arch']) + if args["--arch"] in VALID_ARCHS: + ind = VALID_ARCHS.index(args["--arch"]) arch = VALID_ARCHS[ind] else: raise Exception( - f"{args['--arch']} is not one of the valid architectures {VALID_ARCHS}" + f"{args['--arch']} is not one of the valid architectures {VALID_ARCHS}", ) return root_path, dist_path, mcc_path, arch def write_version(afile: str, version: str) -> None: - with open(afile, 'w') as file: - file.writelines([version + '\n']) + with open(afile, "w") as file: + file.writelines([version + "\n"]) -if __name__ == '__main__': - args = docopt(__doc__, version='0.1') +if __name__ == "__main__": + args = docopt(__doc__, version="0.1") root_path, dist_path, mcc_path, arch = get_args(args) repo_info = git_info(root_path) print("Starting build process.") print(f"Marking {repo_info['version']} as the standalone version") - write_version(VERSION_FILE, repo_info['version']) + write_version(VERSION_FILE, repo_info["version"]) # output name of binary is restricted in mcc tmp_name = f"imosToolbox_{ARCH_NAMES[arch]}" @@ -217,8 +279,9 @@ def write_version(afile: str, version: str) -> None: mfiles, matfiles = get_required_files(root_path, repo_info) java_call = create_java_call_sig_compile(root_path) - mcc_call = create_mcc_call_sig(mcc_path, root_path, dist_path, mfiles, - matfiles, tmp_name, arch) + mcc_call = create_mcc_call_sig( + mcc_path, root_path, dist_path, mfiles, matfiles, tmp_name, arch + ) print("Current repo information:") print(repo_info) @@ -226,7 +289,6 @@ def write_version(afile: str, version: str) -> None: print(f"Calling {java_call}..") stdout, stderr = run(java_call) if stderr: - __import__('pdb').set_trace() raise Exception(f"{jcall} failed") print(f"Calling {mcc_call}...") @@ -241,11 +303,15 @@ def write_version(afile: str, version: str) -> None: print(f"Updating binary at {root_path}....") # mcc append .exe to end of file - if arch == 'win64': - copy(os.path.join(dist_path, tmp_name + '.exe'), - os.path.join(root_path, tmp_name + '.exe')) + if arch == "win64": + copy( + os.path.join(dist_path, tmp_name + ".exe"), + os.path.join(root_path, tmp_name + ".exe"), + ) else: - copy(os.path.join(dist_path, tmp_name), - os.path.join(root_path, tmp_name + '.bin')) + copy( + os.path.join(dist_path, tmp_name), + os.path.join(root_path, tmp_name + ".bin"), + ) print("Build Finished.") diff --git a/imosToolbox.m b/imosToolbox.m index fab2f9060..e9411d0a6 100644 --- a/imosToolbox.m +++ b/imosToolbox.m @@ -8,7 +8,7 @@ function imosToolbox(auto, varargin) % toolbox is output. if 'auto', the toolbox is executed % automatically, with no user interaction. Any other string will % result in the toolbox being executed normally. -% varargin - In 'auto' mode, any other parameters passed in are passed +% varargin - In 'auto' mode, any other parameters passed in are passed % through to the autoIMOSToolbox function - see the % documentation for that function for details. % @@ -19,7 +19,7 @@ function imosToolbox(auto, varargin) % % -% Copyright (C) 2017, Australian Ocean Data Network (AODN) and Integrated +% Copyright (C) 2017, Australian Ocean Data Network (AODN) and Integrated % Marine Observing System (IMOS). % % This program is free software: you can redistribute it and/or modify @@ -37,7 +37,7 @@ function imosToolbox(auto, varargin) % % Set current toolbox version -toolboxVersion = ['2.6.5 - ' computer]; +toolboxVersion = ['2.6.6 - ' computer]; if nargin == 0, auto = 'manual'; end @@ -49,7 +49,7 @@ function imosToolbox(auto, varargin) path = ''; if ~isdeployed [path, ~, ~] = fileparts(which('imosToolbox.m')); - + % set Matlab path for this session (add all recursive directories to Matlab % path) searchPath = textscan(genpath(path), '%s', 'Delimiter', pathsep); @@ -69,7 +69,7 @@ function imosToolbox(auto, varargin) javaaddpath(jars{j}); end -switch auto +switch auto case 'auto', autoIMOSToolbox(toolboxVersion, varargin{:}); otherwise, flowManager(toolboxVersion); end diff --git a/imosToolbox_Linux64.bin b/imosToolbox_Linux64.bin index ca733e8c3..7933c6cd4 100755 Binary files a/imosToolbox_Linux64.bin and b/imosToolbox_Linux64.bin differ diff --git a/imosToolbox_Win64.exe b/imosToolbox_Win64.exe index fc306c762..52b3444d3 100644 Binary files a/imosToolbox_Win64.exe and b/imosToolbox_Win64.exe differ diff --git a/test/AutomaticQC/test_SpikeClassifiers.m b/test/AutomaticQC/test_SpikeClassifiers.m new file mode 100644 index 000000000..5c9b68b09 --- /dev/null +++ b/test/AutomaticQC/test_SpikeClassifiers.m @@ -0,0 +1,324 @@ +classdef test_SpikeClassifiers < matlab.unittest.TestCase + + properties (TestParameter) + plot_figure = {true}; %,false}; + end + + methods (Test) + + function test_allf_for_single_large_spike_pass(~) + signal = gen_spike_signal(1:1/24:10, [0.1, 10], [5, 1], 0, 0); + signal(15) = signal(15) * std(signal) * 10; + r_hampel = imosSpikeClassifierHampel(signal); + r_otsu = imosSpikeClassifierOTSU(signal); + r_savgol_otsu = imosSpikeClassifierNonBurstSavGolOTSU(signal); + loc = unique([r_hampel, r_otsu, r_savgol_otsu]); + assert(isequal(loc, 15)); + end + + function test_allf_low_std_otsu_fails(~) + signal = gen_spike_signal(1:1/24:10, [0.1, 10], [5, 1], 0, 0); + signal(15) = signal(15) * std(signal); + r_hampel = imosSpikeClassifierHampel(signal); + r_otsu = imosSpikeClassifierOTSU(signal); + r_savgol_otsu = imosSpikeClassifierNonBurstSavGolOTSU(signal); + assert(isequal(r_hampel, 15)) + assert(~isequal(r_otsu, 15))% the computed threshold is not enough to split signal and spike + assert(isequal(r_savgol_otsu, 15)) + + %adhoc solution - look iteratively for the largest threshold + nbins = maximize_otsu_threshold(signal); + r_otsu = imosSpikeClassifierOTSU(signal, nbins); + assert(isequal(r_otsu, 15)) + + end + + function test_allf_multiple_spikes_savgol_fails(~) + signal = gen_spike_signal(1:1/24:10, [.1, 10], [5, 1], 0, 0); + signal(15) = signal(15) * std(signal) * 10; + signal(20) = signal(20) * std(signal) * 10; + signal(30) = signal(30) * std(signal) * 10; + r_hampel = imosSpikeClassifierHampel(signal); + r_otsu = imosSpikeClassifierOTSU(signal); + r_savgol_otsu = imosSpikeClassifierNonBurstSavGolOTSU(signal); + assert(isequal(r_hampel, [15, 20, 30])) + assert(isequal(r_otsu, [15, 20, 30])) + assert(~isequal(r_savgol_otsu, [15, 20, 30])) + % no easy solution to overcome savgol otsu - the exact locations never match + end + + function test_allf_multiple_spikes_similar_amplitude_savgol_fails(~) + signal = gen_spike_signal(1:1/24:10, [10.1, 10], [5, 1], 0, 0); + signal(15) = signal(15) * std(signal) * 10; + signal(20) = signal(20) * std(signal) * 10; + signal(30) = signal(30) * std(signal) * 10; + r_hampel = imosSpikeClassifierHampel(signal); + r_otsu = imosSpikeClassifierOTSU(signal); + r_savgol_otsu = imosSpikeClassifierNonBurstSavGolOTSU(signal); + assert(isequal(r_hampel, [15, 20, 30])) + assert(isequal(r_otsu, [15, 20, 30])) + assert(~isequal(r_savgol_otsu, [15, 20, 30])) + % no easy solution - one spike always leaks in savgol. + end + + function test_allf_multiple_spikes_with_different_stds_hampel_pass(~) + signal = gen_spike_signal(1:1/24:10, [10.1, 10], [5, 1], 0, 0); + signal(15) = signal(15) * std(signal) * 10; + signal(20) = signal(20) * std(signal) * 5; + signal(30) = signal(30) * std(signal) * 1; + r_hampel = imosSpikeClassifierHampel(signal); + r_otsu = imosSpikeClassifierOTSU(signal); + r_savgol_otsu = imosSpikeClassifierNonBurstSavGolOTSU(signal); + assert(isequal(r_hampel, [15, 20, 30])) + assert(~isequal(r_otsu, [15, 20, 30])) + assert(~isequal(r_savgol_otsu, [15, 20, 30])) + % there are no parameters available to make exactly the spikes in savgol. + end + + function test_allf_multiple_spikes_with_different_stds_hampel_fails(~) + signal = gen_spike_signal(1:1/24:10, [10.1, 10], [5, 1], 0, 0); + signal(15) = signal(15) * std(signal) * 1; + signal(20) = signal(20) * std(signal) * 1; + signal(30) = signal(30) * std(signal) * 1; + r_hampel = imosSpikeClassifierHampel(signal); + r_otsu = imosSpikeClassifierOTSU(signal); + r_savgol_otsu = imosSpikeClassifierNonBurstSavGolOTSU(signal); + assert(~isequal(r_hampel, [15, 20, 30])) + assert(~isequal(r_otsu, [15, 20, 30])) + assert(~isequal(r_savgol_otsu, [15, 20, 30])) + % feasible solution - reduce the standard deviation factor to detect lower + % variability spikes + r_hampel = imosSpikeClassifierHampel(signal, 3, 1); + assert(isequal(r_hampel, [15, 20, 30])); + % solution is feasible even with different windows, since signal is simple + r_hampel = imosSpikeClassifierHampel(signal, 6, 1); + assert(isequal(r_hampel, [15, 20, 30])); + % However, one needs to know the approximate scale of the spike noise + % or it will overflow return false positives. + r_hampel = imosSpikeClassifierHampel(signal, 12, 1); + assert(~isequal(r_hampel, [15, 20, 30])); + end + + function test_hampel_robust_in_multi_spike_detection(~, plot_figure) + time = 0:1/24:10; + nspikes = 25; + mfac = 5; + [osig, ssig, spikes] = gen_spike_signal(time, [0.2, 8, 2], [1/24, 5, 30], nspikes, mfac); + hs = imosSpikeClassifierHampel(ssig, nspikes / 5, mfac * 1.2); + + nbins = maximize_otsu_threshold(ssig); + os = imosSpikeClassifierOTSU(ssig, nbins); + ss = imosSpikeClassifierNonBurstSavGolOTSU(ssig, nspikes / 5, 2, nbins); + + assert(isequal(hs, spikes)) + assert(~isequal(os, spikes)) + assert(~isequal(ss, spikes)) + + if plot_figure + plot_spikes(time, osig, ssig, spikes, hs) + title('Simple spike detection - Hampel') + hold on + plot_spikes(time, osig, ssig, spikes, os) + title('Simple spike detection - OTSU') + plot_spikes(time, osig, ssig, spikes, ss) + title('Simple spike detection - Savgol-OTSU') + + end + + end + + function test_OTSU_noise_spike_detection_is_incomplete(~, plot_figure) + time = 0:1/24:10; + nspikes = 25; + mfac = 1; + [osig, ssig, spikes] = gen_spike_signal(time, [20, 5, 20, 1], [4/24, 12/24, 5, 30], nspikes, mfac); + dspikes = imosSpikeClassifierOTSU(ssig); + + if plot_figure + plot_spikes(time, osig, ssig, spikes, dspikes) + title('Incomplete Spike Detection - no False Positives - OTSU') + end + + matches = intersect(dspikes, spikes); + incomplete_detection = length(matches) < length(spikes); + ratio_detection = length(matches) / length(spikes); + no_false_positives = isempty(setdiff(dspikes, spikes)); + assert(incomplete_detection) + assert(ratio_detection > 0.7) + assert(no_false_positives) + assert(length(dspikes) == 19) + end + + function test_hampel_noise_spike_detection_incomplete(~, plot_figure) + time = 0:1/24:10; + nspikes = 25; + mfac = 1; + [osig, ssig, spikes] = gen_spike_signal(time, [20, 5, 20, 1], [4/24, 12/24, 5, 30], nspikes, mfac); + dspikes = imosSpikeClassifierHampel(ssig); + + if plot_figure + plot_spikes(time, osig, ssig, spikes, dspikes) + title('Incomplete Spike Detection - Hampel') + end + + matches = intersect(dspikes, spikes); + incomplete_detection = length(matches) < length(spikes); + ratio_detection = length(matches) / length(spikes); + no_false_positives = isempty(setdiff(dspikes, spikes)); + assert(incomplete_detection) + assert(ratio_detection < 0.5) + assert(no_false_positives) + assert(length(dspikes) == 9) + + end + + function test_savgolOTSU_noise_spike_detection_has_false_positives(~, plot_figure) + time = 0:1/24:10; + nspikes = 25; + mfac = 1; + [osig, ssig, spikes] = gen_spike_signal(time, [20, 5, 20, 1], [4/24, 12/24, 5, 30], nspikes, mfac); + dspikes = imosSpikeClassifierNonBurstSavGolOTSU(ssig); + + if plot_figure + plot_spikes(time, osig, ssig, spikes, dspikes) + title('Incomplete Spike Detection - with False positives - SavGol-OTSU') + end + + matches = intersect(dspikes, spikes); + incomplete_detection = length(matches) < length(spikes); + ratio_detection = length(matches) / length(spikes); + false_positives = ~isempty(setdiff(dspikes, spikes)); + assert(incomplete_detection) + assert(ratio_detection > 0.5) + assert(false_positives) + assert(length(dspikes) == 18) + end + + function test_hampel_under_detection_false_positives(~, plot_figure) + time = 0:1/24:10; + nspikes = 25; + mfac = 1; + [osig, ssig, spikes] = gen_spike_signal(time, [20, 5, 20, 1], [4/24, 12/24, 5, 30], nspikes, mfac); + dspikes = imosSpikeClassifierHampel(ssig, 6, mfac * 1.2); % + + if plot_figure + plot_spikes(time, osig, ssig, spikes, dspikes) + title('Incomplete Spike detection - Hampel method') + end + + matches = intersect(dspikes, spikes); + incomplete_detection = length(matches) < length(spikes); + ratio_detection = length(matches) / length(spikes); + false_positives = ~isempty(setdiff(dspikes, spikes)); + nfalse = length(setdiff(dspikes, spikes)); + assert(incomplete_detection) + assert(ratio_detection > 0.9) + assert(false_positives) + assert(length(dspikes) == 44) + assert(length(nfalse) == 1) + end + + end + +end + +function [osig, ssig, spikes] = gen_spike_signal(time, amps, periods, nspikes, stdmag) +% function e = gen_spike_signal() +% +% Generate a periodic signal with fixed interval and +% a fixed number of randomized spikes with amplitude proportional to stdmag and +% the distribution of the amplitudes provided. +% +% Inputs: +% time - a time range, time series, or sorted datenum array. +% The smallest time interval is used for generation. +% amps - a vector with normalized signal amplitudes. +% periods - a vector with the signal period (days). +% nspikes - the number of spikes to introduce. +% If missing, the number of spikes is a +% random number between 0 and half the length of the time vector. +% stdmag - a standard deviation scale for the spike. +% +% Outputs: +% +% osig - the clean signal. +% ssig - the spiked signal. +% spikes - the location of the spikes. +% +% author: hugo.oliveira@utas.edu.au +% + +rng(0, 'simdTwister')%simd is faster + +if length(amps) ~= length(periods) + error("Amplitude and Periods shou should be the same length") +end + +dts = timeQuantisation(unique(diff(time)), 'second'); +dt = dts(1); +stime = 0:dt:max(time); +osig = zeros(size(stime)); + +if nargin < 4 + nspikes = randi([0, floor(length(stime) / 2)]); + stdmag = 1; +end + +[~, time_ind] = intersect(stime, time); + +for k = 1:length(amps) + a = amps(k); + omega = 2 * pi / periods(k); + osig = osig + a * cos(omega * stime); +end + +if nspikes == 0 + return +end + +spikes = 0; + +while length(spikes) ~= nspikes + spikes = unique(randi([1, length(time_ind)], 1, nspikes)); +end + +if length(amps) > 1 + rmax = max(amps) * std(amps) * stdmag; + rmin = -1 * min(amps) * std(amps) * stdmag; +else + rmax = amps(1) * stdmag; + rmin = -1 * rmax; +end + +spike_mag = abs(rmin) - (rmax - rmin) * rand(1, nspikes); +osig = osig(time_ind); +ssig = osig; +ssig(spikes) = ssig(spikes) + spike_mag; + +end + +function plot_spikes(time, osig, ssig, spikes, dspikes) +figure(); +plot(time, ssig, 'r') +hold on +plot(time, osig, 'b') +plot(time(spikes), ssig(spikes), 'rs') +plot(time(dspikes), ssig(dspikes), 'k*') +legend('spiked signal', 'clean signal', 'real spikes', 'detected spikes') +end + +function nbins = maximize_otsu_threshold(signal) +ot = 0; +nbins = 0; + +for k = 1:512 + new_ot = abs(otsu_threshold(signal, k)); + + if new_ot > ot + ot = new_ot; + nbins = k; + end + +end + +end diff --git a/test/Util/StructUtils/test_combineStructFields.m b/test/Util/StructUtils/test_combineStructFields.m new file mode 100644 index 000000000..031c3f627 --- /dev/null +++ b/test/Util/StructUtils/test_combineStructFields.m @@ -0,0 +1,34 @@ +classdef test_combineStructFields < matlab.unittest.TestCase + + methods (Test) + + function test_basic(~) + a = struct(); b=a; + a.x = 1; + b.x = 2; + b.y = 3; + combined = combineStructFields(a,b); + assert(combined.x==2); + assert(combined.y==3); + end + + function test_multi_structure_case(~) + % multi structure case + a = struct(); b=a; c=a; + a.x = 1; + b.x = 2; + b(2).x = 3; + b(3).x = 4; + c.x = 3; + c(2).x = 4; + c(3).x = 5; + combined = combineStructFields(a,b,c); + assert(isequal(size(combined),[1,3])) + assert(combined(1).x==3) + assert(combined(2).x==4) + assert(combined(3).x==5) + end + + end + +end diff --git a/test/Util/Time/test_timeSamplingInfo.m b/test/Util/Time/test_timeSamplingInfo.m new file mode 100644 index 000000000..ab4429477 --- /dev/null +++ b/test/Util/Time/test_timeSamplingInfo.m @@ -0,0 +1,97 @@ +classdef test_timeSamplingInfo < matlab.unittest.TestCase + + methods (Test) + + function test_basic(~) + info = timeSamplingInfo([1, 2, 3, 4, 5.3], 'day'); + assert(info.uniform_sampling) + assert(all([info.unique_sampling, info.monotonic_sampling, info.progressive_sampling])) + end + + function test_symmetric_cases(~) + time = [1, 2, 1, 3, 1, 4, 1, 5, 1]; + info = timeSamplingInfo(time, 'day'); + assert(~all([info.unique_sampling, info.progressive_sampling, info.regressive_sampling, info.monotonic_sampling])) + assert(info.sampling_steps_median==0) % symmetric values + assert(info.sampling_steps_mode==-4) %sorted first occurrence + assert(all(info.jump_magnitude == [1, -1, 2, -2, 3, -3, 4, -4])) + assert(all(info.jump_forward_indexes == [1 0 1 0 1 0 1 0])) + assert(all(info.jump_backward_indexes == [0 1 0 1 0 1 0 1])) + + end + + function test_perfect_sampling(~) + time = 1:10; + info = timeSamplingInfo(time, 'day'); + assert(all([info.consistent_sampling, info.unique_sampling, info.uniform_sampling, info.monotonic_sampling, info.progressive_sampling])) + assert(isequal(info.number_of_sampling_steps,1)) + assert(isequal(info.sampling_steps, 1)) + assert(isequal(info.sampling_steps_median, 1)) + assert(isequal(info.sampling_steps_mode, 1)) + end + + function test_detect_time_jumps(~) + time = [1:10,60:80,130:140]; + info = timeSamplingInfo(time,'day'); + assert(all([info.consistent_sampling,info.unique_sampling,info.monotonic_sampling, info.progressive_sampling])) + assert(isequal(info.jump_magnitude,[50,50])) + end + + function test_null_case(~) + time = repmat(0,10,1); + info = timeSamplingInfo(time,'day'); + assert(info.number_of_sampling_steps==1); + assert(info.sampling_steps_median==0); + assert(info.sampling_steps_mode==0); + assert(all([info.consistent_sampling,info.constant_vector,info.repeated_samples,info.uniform_sampling])) + assert(isempty(info.jump_magnitude)) + assert(isempty(info.jump_forward_indexes)) + end + + function test_repeated_time_entries(~) + time = [1:10,60,60,60,60,60,65:80,130:140]; + info = timeSamplingInfo(time,'day'); + assert(isequal(info.sampling_steps,[0,1,5,50])) + assert(info.sampling_steps_mode==1) + assert(all([info.consistent_sampling,info.repeated_samples,info.progressive_sampling])) + assert(isequal(info.jump_magnitude,[50, 0, 0, 0, 0, 5, 50])) + assert(isequal(info.jump_forward_indexes,[1, 0, 0, 0, 0, 1, 1])) + end + + function test_ignore_microsecond(~) + info = timeSamplingInfo([1, 2, 3, 4 + 0.49/86400, 5], 'second'); + assert(info.repeated_samples == 0) + assert(info.sampling_steps == 1) + end + + function test_detect_sampling_above_scale(~) + info = timeSamplingInfo([1, 2, 3, 4 + 0.5/86400, 5], 'second'); + assert(info.repeated_samples == 0) + assert(isequal(info.sampling_steps * 86400, [86400, 86401])) + end + + function test_ignore_sequential_offsets_below_scale(~) + drifts = [0, 0, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1] / 86400; + times = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]; + info = timeSamplingInfo(times + drifts, 'second'); + assert(info.sampling_steps_in_seconds == 86400) + end + + function test_detect_offset_above_scale(~) + drifts = [0, 0, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.9] / 86400; + times = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]; + info = timeSamplingInfo(times + drifts, 'second'); + assert(length(info.sampling_steps) == 2) + assert(isequal(info.sampling_steps_in_seconds, [86400, 86401])) + end + + function test_several_sampling_drifts(~) + drifts = [0, 0.5, 0.5, 0.5, 0.75, 1, 1.25, 1.5, 3, 6, 9, 12, 24] / 86400; + times = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]; + info = timeSamplingInfo(times + drifts, 'second'); + assert(isequal(info.sampling_steps_in_seconds, 86400 + [0, 1, 2, 3, 12])) + end + + end + +end diff --git a/test/Util/testtimeSamplingInfo.m b/test/Util/testtimeSamplingInfo.m deleted file mode 100644 index 235af5189..000000000 --- a/test/Util/testtimeSamplingInfo.m +++ /dev/null @@ -1,42 +0,0 @@ -classdef testtimeSamplingInfo < matlab.unittest.TestCase - - methods (Test) - - function test_simple_monotonic(testCase), - info = timeSamplingInfo([1, 2, 3, 4]); - assert(info.uniform_sampling); - assert(all([info.unique_sampling, info.monotonic_sampling, info.progressive_sampling])); - end - - function test_decreasing_with_jump_at_end(testCase), - info = timeSamplingInfo([4, 3, 2, 1, -1]); - assert(~info.uniform_sampling); - assert(all(info.sampling_steps == [-2, -1])); - assert(all(info.sampling_steps_median == -1)); - assert(info.regressive_sampling); - assert(all([info.unique_sampling, info.monotonic_sampling])); - end - - function test_increasing_with_jump_at_end(testCase), - info = timeSamplingInfo([1, 2, 3, 5]); - assert(info.non_uniform_sampling_indexes == 3); - end - - function test_increasing_with_repeat(testCase), - info = timeSamplingInfo([1, 2, 3, 3]); - assert(~info.uniform_sampling); - assert(~all([info.unique_sampling, info.monotonic_sampling])); - assert(info.non_uniform_sampling_indexes == [3]); - end - - function test_jump_forward_backward(testCase), - info = timeSamplingInfo([1, 2, 1, 3, 1, 4, 1, 5, 1]); - assert(~all([info.unique_sampling, info.progressive_sampling, info.regressive_sampling, info.monotonic_sampling])); - assert(isequal(info.sampling_steps_median, [])); - assert(all(info.jump_forward_indexes == [1 0 1 0 1 0 1 0])); - assert(all(info.jump_backward_indexes == [0 1 0 1 0 1 0 1])); - end - - end - -end