text:
function ImP_SQI_mimic_data_importer
% ImP_SQI_mimic_data_importer imports data from the MIMIC II database with
% which to assess the performance of the impedance pneumography SQI, and
% outputs the dataset in Matlab format.
%
%               ImP_SQI_mimic_data_importer
%
%	Inputs:
%       Just specify the relevant paths in the "universal_parameters"
%       function below. No other inputs are required as this script
%       downloads the required data files from PhysioNet automatically.
%
%	Outputs:
%       A Matlab file containing the collated data, called "mimic_imp_sqi_data.mat"
%
%   The WFDB Toolbox:
%       This requires the WFDB Toolbox, which can be downloaded from:
%           https://archive.physionet.org/physiotools/matlab/wfdb-app-matlab/
%       I installed it using the 'Quick Start' instructions at:
%           https://archive.physionet.org/physiotools/matlab/wfdb-app-matlab/#quick-start
%               (including modifying the download path to be https://archive.physionet.org ... )
%       The Toolbox is now available at:
%           https://physionet.org/content/wfdb-swig-matlab/1.0.0/
%
%   Further Information:
%       This version of the ImP_SQI_mimic_data_importer is provided to facilitate
%       reproduction of the MIMIC subset used in:
%           Charlton P. H. et al., "An impedance pneumography signal
%           quality index for respiratory rate monitoring: design,
%           assessment and application", [under review]
%       Further information on this study can be obtained at:
%           http://peterhcharlton.github.io/RRest/imp_sqi.html
%       In addition, further information on RRest, a toolbox of respiratory
%       algorithms which can be used with this dataset, can be obtained at:
%           http://peterhcharlton.github.io/RRest/index.html
%
%   Comments, Questions, Criticisms, Feedback, Contributions:
%       See: http://peterhcharlton.github.io/RRest/contributions.html
%
%   Version:
%       v.0.1 - accompanying peer review, 6th Aug 2020 by Peter H Charlton
%
%   Source:
%       This script has been adapted from 'MIMICII_data_importer.m', which
%       is one of the scripts in the RRest toolbox, available at:
%           https://github.com/peterhcharlton/RRest
%
%   Licence:
%       This program is available under the GNU public license. 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, either version 3 of the License, or
%       (at your option) any later version. This program is distributed in
%       the hope that it will be useful, but WITHOUT ANY WARRANTY; without
%       even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%       PARTICULAR PURPOSE.  See the GNU General Public License for more
%       details: <http://www.gnu.org/licenses/>.
%

%% Setup
% ~~~ This function needs adjusting by the individual user as it contains paths
% specific to their computer ~~~
up = universal_parameters;

%% Create ReadMe and License files
create_readme_and_license_files(up);

%% Download data
download_data(up);

%% Extract data from files on computer
extract_from_mimic_ii(up);

%% Add in manual annotations
add_manual_annotations(up);

% %% Save data in RRest format
% if up.output.rrest_format, save_data_in_rrest_format(up); end
%
% %% Save data in CSV format
% if up.output.csv_format, save_data_in_csv_format(up); end
%
% %% Save data in WFDB format
% if up.output.csv_format, save_data_in_wfdb_format(up); end

fprintf(['\n\n -- Finished importing the ImP SQI dataset. Data saved at:\n\n       ', up.paths.data_root, '\n\n'])

end

function up = universal_parameters

fprintf('\n -- Setting up Universal Parameters')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% PARAMETERS TO BE SPECIFIED %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Specify the root data directory (where the data will be stored)
up.paths.data_root = '/Users/petercharlton/Documents/Data/Imp_SQI_mimic_dataset_test/';
fprintf('\n   - The data will be stored at %s', up.paths.data_root);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%   OTHER PARAMETERS   %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%   (don't need editing)   %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

provide_details_of_script;

% Specify the web address of the data to be downloaded
rel_database = 'mimic3wdb/matched';
rel_database_part = 'p00';
up.paths.database_dir_root = ['https://archive.physionet.org/physiobank/database/', rel_database, '/'];
up.paths.database_dir = [up.paths.database_dir_root, rel_database_part];
up.paths.analysis_path = [up.paths.data_root, 'Analysis_files', filesep, 'Data_for_Analysis', filesep ];
up.paths.records_filepath = [up.paths.analysis_path, 'RECORDS'];
up.paths.download_details_file = [up.paths.analysis_path, 'download_details'];
up.no_pt_stays = 100;

% extraction definitions
up.extraction.rel_sigs = {'II', 'PLETH', 'RESP'};    % required signals
up.extraction.rel_nums = {'HR', 'PULSE', 'RESP', 'SpO2'};    % required numerics
up.fs = 125;  % Hz
up.req_data_time = 0*60;  % data taken from 0 minutes into record (secs)
up.req_data_duration = 60*60;  % segments of data extracted (secs)
up.dataset_name = 'imp_sqi_mimic';
up.csv.new_line = '\n';
up.units.signals = {'II', 'AVR', 'V',  'III', 'I',  'MCL1', 'MCL', 'PLETH', 'RESP', 'ABP',  'CVP',  'ART',   'P1',  'UAP',  'PAP'};
up.units.units   = {'mV', 'mV',  'mV', 'mV',  'mV', 'mV',   'mV',  'NU',    'pm'  , 'mmHg', 'mmHg', 'mmHg', 'mmHg', 'mmHg', 'mmHg'};
up.units.numerics    = {'HR',  'PULSE',  'SpO2', 'RESP', 'ABPSys', 'ABPDias', 'ABPMean', 'NBPSys', 'MBPDia', 'MBPMean'};
up.units.num_units   = {'bpm', 'bpm',    '%',    'pm',   'mmHg',   'mmHg',    'mmHg',    'mmHg',   'mmHg',   'mmHg'   };

% paths
up.paths.file_location = [up.paths.data_root, 'physionet.org', filesep, 'physiobank' ,filesep 'database', filesep];
up.paths.extracted_data_file = [up.paths.data_root, 'mimic_imp_sqi_data.mat'];

% url paths
up.paths.url.records = [up.paths.database_dir_root, 'RECORDS'];

% Add all functions within the directory
addpath(genpath(fileparts(mfilename('fullpath'))));

% check that all the required folders exist. If not, create them:
if ~exist(up.paths.data_root, 'dir')
    fprintf('\n The specified data folder does not exist \n - please either: \n 1) select a different location (speifying it in "universal_parameters"), or\n 2) create a folder at the specified location: %s', up.paths.data_root);
    error('Specified folder does not exist')
end
req_folders = {up.paths.analysis_path};
displayed_msg = 0;
for req_folder_no = 1 : length(req_folders)
    curr_folder = req_folders{req_folder_no};
    if ~exist(curr_folder, 'dir')
        if displayed_msg == 0
            fprintf('\n -- Creating folder in which to store data and analysis')
            displayed_msg = 1;
        end
        mkdir(curr_folder)
    end
end

% Type(s) of output (currently not used)
up.output.wfdb_format = true; % Save in PhysioNet's waveform database format
up.output.csv_format = true;  % Save as comma-separated-value format
up.output.rrest_format = true;% Save in a format suitable for use with the RRest Toolbox of algorithms

% Check that the WFDB Matlab Toolbox has been installed
if ~exist('getWfdbClass', 'file')
    error('Couldn''t find the WFDB Matlab Toolbox. Please install as described at the top of the file.')
end

end

function provide_details_of_script

licence_details = ['\n\n ImP_SQI_data_importer', ...
    '\n This program comes with ABSOLUTELY NO WARRANTY', ...
    '\n and is available under the GNU public license.', ...
    '\n For details see the accompanying LICENSE.txt file.\n\n'];

fprintf(licence_details)

end

function download_data(up)

fprintf('\n -- Downloading Data:')

%% set current directory to that of wget:
cd(up.paths.data_root)

%% Find out what records are available:
fprintf(' List of records, ');
% download list of records
if ~exist(up.paths.records_filepath, 'file')
    outfilename = websave(up.paths.records_filepath,up.paths.url.records);
end

% form list of patient stays from downloaded file
pt_stays = importdata(up.paths.records_filepath);

%% Download each record
file_counter = 0;

for pt_stay_no = 1 : length(pt_stays)

    % specify patient stay
    pt_stay = pt_stays{pt_stay_no};
    temp = strfind(pt_stay, '/');
    database_subset = pt_stay(1:temp(1)-1);
    pt_no = pt_stay(temp(1)+1:end-1); clear temp

    % set save location
    save_folder = [up.paths.file_location, pt_stay];
    if ~exist(save_folder, 'dir')
        mkdir(save_folder)
    end

    % download records for this pt stay
    pt_stay_records_url = [up.paths.database_dir_root, pt_stay, 'RECORDS'];
    save_path = [save_folder,'RECORDS'];
    if ~exist(save_path, 'file')
        outfilename = websave(save_path, pt_stay_records_url);
    end
    pt_stay_records = importdata(save_path);

    % download overall header file
    pt_stay_segment = pt_stay_records{1};
    curr_file_name = [pt_stay_segment, '.hea'];
    overall_header_file = [up.paths.database_dir_root, pt_stay, curr_file_name];
    overall_header_file_save_loc = [save_folder, curr_file_name];
    if ~exist(overall_header_file_save_loc, 'file')
        outfilename = websave(overall_header_file_save_loc,overall_header_file);
    end

    % find out what individual files are in this pt stay:
    fileID = fopen(overall_header_file_save_loc);
    individual_files.name = {};
    individual_files.duration = []; individual_files.cum_duration_before = []; individual_files.cum_duration_after = [];
    line_no = 0;
    while ~feof(fileID)
        curr_text_line = fgetl(fileID); line_no = line_no+1;
        if line_no == 1
            temp = strfind(curr_text_line, ':');
            individual_files.start_time = curr_text_line(temp(1)-2:temp(1)+9);
        elseif ~isempty(strfind(curr_text_line, '_')) & length(strfind(curr_text_line, ' '))==1 & isempty(strfind(curr_text_line, 'layout'))
            individual_files.name{end+1} = curr_text_line(1:strfind(curr_text_line, '_')+4);
            individual_files.duration(end+1) = str2double(curr_text_line(strfind(curr_text_line, '_')+6:end));
        end
        clear curr_text_line
    end
    clear line_no
    fclose(fileID); clear fileID

    % Find out which files contain sufficient data
    rel_file_els = find(individual_files.duration >= (up.req_data_duration*up.fs+1));
    rel_file_names = individual_files.name(rel_file_els);

    % cycle through each file containing sufficient data
    for file_no = 1 : length(rel_file_names)
        curr_rel_file = rel_file_names{file_no};

        % download data header file
        curr_file_name = [curr_rel_file, '.hea'];
        down_loc = [up.paths.database_dir_root, pt_stay, curr_file_name];
        save_loc = [save_folder, curr_file_name];
        if ~exist(save_loc, 'file')
            outfilename = websave(save_loc,down_loc);
        end
        data_header_file = save_loc;

        % identify the signals in this file
        fileID = fopen(save_loc);
        line_no = 0;
        signals = [];
        while ~feof(fileID)
            curr_text_line = fgetl(fileID); line_no = line_no+1;
            if line_no == 1
                temp = strfind(curr_text_line, ':');
                curr_file.start_time = curr_text_line(temp(1)-2:end);
            elseif contains(curr_text_line, '.dat')
                temp = strfind(curr_text_line, ' ');
                signals{end+1} = curr_text_line(temp(end)+1:end);
            end
            clear curr_text_line
        end
        clear line_no
        fclose(fileID); clear fileID

        % check to see whether this file contains the required signals
        if sum(strcmp(signals, 'II')) && sum(strcmp(signals, 'PLETH')) && sum(strcmp(signals, 'RESP'))

            % store details
            file_counter = file_counter+1;
            fprintf([' ' num2str(file_counter)]);
            download.source.waves_header_file{file_counter} = down_loc;
            download.source.overall_header_file{file_counter} = overall_header_file;
            download.saved.pt_stay{file_counter} = pt_stay;
            download.saved.pt_no{file_counter} = pt_no;
            download.saved.pt_stay_segment{file_counter} = pt_stay_segment;
            download.saved.waves_header_file{file_counter} = data_header_file;
            download.saved.root_folder{file_counter} = save_folder;
            download.saved.overall_header_file{file_counter} = overall_header_file_save_loc;

            % download data file
            curr_file_name = [curr_rel_file, '.dat'];
            download.saved.file_name{file_counter} = curr_file_name;
            down_loc = [up.paths.database_dir_root, pt_stay, curr_file_name];
            download.source.waves_data_file{file_counter} = down_loc;
            download.source.waves_data_filename{file_counter} = curr_file_name;
            save_loc = [save_folder, curr_file_name];
            if ~exist(save_loc, 'file')
                outfilename = websave(save_loc,down_loc);
            end
            download.saved.waves_data_file{file_counter} = save_loc;
            clear outfilename save_loc down_loc curr_file_name

            % download numerics header file
            curr_file_name = [pt_stay_segment, 'n.hea'];
            down_loc = [up.paths.database_dir_root, pt_stay, curr_file_name];
            save_loc = [save_folder, curr_file_name];
            if ~exist(save_loc, 'file')
                outfilename = websave(save_loc,down_loc);
                download.source.numerics_header_file{file_counter} = down_loc;
            end
            download.saved.numerics_header_file{file_counter} = save_loc;
            clear outfilename save_loc down_loc curr_file_name

            % download numerics data file
            temp = strfind(curr_rel_file, '_');
            curr_file_name = [curr_rel_file(1:temp-1), 'n.dat']; clear temp
            down_loc = [up.paths.database_dir_root, pt_stay, curr_file_name];
            download.source.numerics_data_file{file_counter} = down_loc;
            save_loc = [save_folder, curr_file_name];
            if ~exist(save_loc, 'file')
                outfilename = websave(save_loc,down_loc);
            end
            download.saved.numerics_data_file{file_counter} = save_loc;
            clear outfilename save_loc down_loc curr_file_name

            % skip any remaining files for this subject
            break

        end

        clear outfilename save_loc down_loc curr_file_name
    end

    % finish if reached the desired number of pt stays with data
    if file_counter >= up.no_pt_stays
        break
    end

end

% Save the download details
save(up.paths.download_details_file, 'download');

end

function extract_from_mimic_ii(up)
%
% Based on RDSAMP, part of the WFDB Toolbox. Further details are available
% at:
%       http://physionet.org/physiotools/matlab/wfdb-app-matlab/
%
% Reads MIMIC II data and outputs it as a structure of .t and .v values.

fprintf('\n -- Extracting Data from MIMIC files: ')

% Load download details
load(up.paths.download_details_file);

for pt_stay_no = 1 : length(download.saved.pt_stay)
    fprintf([' ' num2str(pt_stay_no)]);

    % set current dir
    dir_path = download.saved.root_folder{pt_stay_no};
    cd(dir_path)

    % setup java
    % persistent javaWfdbExec config
    [javaWfdbExec, config] = deal([]);
    if(isempty(javaWfdbExec))
        [javaWfdbExec,config]=getWfdbClass('rdsamp');
    end

    % Remove file extension
    record_str = download.saved.file_name{pt_stay_no}(1:end-4);

    % Extract file information (including which signals are present)
    siginfo = extract_file_info(record_str);

    % create wfdb argument
    N0 = (siginfo.fs * up.req_data_time);
    wfdb_argument={'-r',record_str,'-Ps','-f',['s' num2str(N0)]}; clear record_str

    % calculate start time of numerics data
    [numerics_start, waveforms_start] = calc_start_times(download, pt_stay_no);

    % Make command to only extract the required signals
    wfdb_argument{end+1}='-s ';
    %-1 is necessary because WFDB is 0 based indexed.
    for sInd=1:siginfo.no_sigs
        wfdb_argument{end+1} = num2str(sInd-1);
    end
    clear sInd

    % determine how many samples to extract from this file
    N= N0 +(siginfo.fs * up.req_data_duration);

    % additional bits of wfdb_argument
    wfdb_argument{end+1}='-t';
    wfdb_argument{end+1}=['s' num2str(N+1)];

    % extract waveform data
    extracted_data=javaWfdbExec.execToDoubleArray(wfdb_argument);
    clear wfdb_argument

    % Store waveform data
    waves.fs = siginfo.fs;
    for wave_no = 1 : length(siginfo.sig_names)
        eval(['waves.' siginfo.sig_names{wave_no} '.t = extracted_data(:, 1)-(N0/125);']);
        eval(['waves.' siginfo.sig_names{wave_no} '.v = extracted_data(:, wave_no+1);']);
    end
    clear wave_no N N0 current_waves extracted_data Fs ListCapacity signalList javaWfdbExec config

    %% Extract corresponding numerics

    % Extract file information (including which signals are present)
    header_filepath = download.saved.numerics_header_file{pt_stay_no};
    do_numerics = 1;
    if ~isempty(header_filepath) && do_numerics
        header_filepath = header_filepath(1:end-4);
        numinfo = extract_file_info(header_filepath);

        % setup java
        persistent javaWfdbExec_num config_num
        if(isempty(javaWfdbExec_num))
            [javaWfdbExec_num,config_num]=getWfdbClass('rdsamp');
        end

        % Identify relevant numerics
        num_inds = 1:numinfo.no_sigs;
        num_log = zeros(length(num_inds),1);
        for num_no = 1 : length(num_inds)
            if sum(strcmp(up.extraction.rel_nums, numinfo.sig_names(num_no)))
                num_log(num_no) = 1;
            end
        end

        numerics.fs = numinfo.fs;

        % create wfdb argument
        if length(numinfo.start_time) == 5, numinfo.start_time = ['00:', numinfo.start_time]; end
%         start_time_diff = (datenum(['00/01/0000 ' numinfo.start_time], 'dd/mm/yyyy HH:MM:SS') - datenum(['00/01/0000 ' siginfo.start_time], 'dd/mm/yyyy HH:MM:SS')) *24*60*60;  % in secs
        start_time_diff = round((waveforms_start - numerics_start) *24*60*60);  % in secs
        N0 = (numerics.fs * up.req_data_time) + (start_time_diff*numerics.fs);
        clear waveforms_start numerics_start start_time_diff
        N = N0 +(numerics.fs * up.req_data_duration);
        [~, header_filepath, ~] = fileparts(header_filepath);
        wfdb_argument={'-r',header_filepath,'-Ps','-f',['s' num2str(N0)]}; clear header_filepath

        % Make command to only extract the required signals
        numList = num_inds(logical(num_log)); clear num_log num_inds
        if(~isempty(numList))
            wfdb_argument{end+1}='-s ';
            %-1 is necessary because WFDB is 0 based indexed.
            for sInd=1:length(numList)
                wfdb_argument{end+1}=[num2str(numList(sInd)-1)];
            end
            clear sInd
        end

        % Additional wfdb_argument bits
        wfdb_argument{end+1}='-t';
        wfdb_argument{end+1}=['s' num2str(N+1)];

        % extract numerics
        extracted_data = javaWfdbExec_num.execToDoubleArray(wfdb_argument);
        clear wfdb_argument

        % store the numerics
        numNames = numinfo.sig_names(numList);
        for num_no = 1 : length(up.extraction.rel_nums)
            rel_num_ind = find(strcmp(numNames, up.extraction.rel_nums{num_no}));
            eval(['numerics.' up.extraction.rel_nums{num_no} '.t = round(extracted_data(:, 1)-N0);']);
            eval(['numerics.' up.extraction.rel_nums{num_no} '.v = extracted_data(:, rel_num_ind+1);']);
            clear rel_num_ind
        end
        data{pt_stay_no}.numerics = numerics;

        clear numerics num_no numinfo siginfo current_numerics extracted_data numNames ListCapacity numinfo numList javaWfdbExec_num config_num N N0

    end

    %% Extract fixed data

    fileID = fopen(download.saved.overall_header_file{pt_stay_no});
    while ~feof(fileID)
        curr_text_line = fgetl(fileID);
        if ~isempty(strfind(curr_text_line, 'Location'))
            fix.loc = curr_text_line(strfind(curr_text_line, ': ')+2:end);
        end
        if ~isempty(strfind(curr_text_line, 'age'))
            fix.age = curr_text_line(strfind(curr_text_line, 'age>: ')+6:strfind(curr_text_line, '<sex>: ')-2);
            if isempty(strfind(fix.age, '?')) && isempty(strfind(fix.age, '+'))
                fix.age = str2double(fix.age);
            elseif strcmp(fix.age, '??')
                fix.age = nan;
            end
        end
        if ~isempty(strfind(curr_text_line, 'sex'))
            fix.sex = curr_text_line(strfind(curr_text_line, 'sex>: ')+6:end);
        end
        clear curr_text_line
    end
    clear line_no
    fclose(fileID); clear fileID

    %% Store data for this patient stay in the output
    data{pt_stay_no}.loc = fix.loc;
    if sum(strcmp(fieldnames(fix), 'age'))
        data{pt_stay_no}.age = fix.age;
    end
    if sum(strcmp(fieldnames(fix), 'sex'))
        data{pt_stay_no}.sex = fix.sex;
    end
    data{pt_stay_no}.ID = download.saved.pt_stay{pt_stay_no};
    data{pt_stay_no}.data_segment = download.saved.pt_stay_segment{pt_stay_no};
    data{pt_stay_no}.waves = waves;
    clear waves numerics fix

end

%% Convert to required variable for analysis

old_data = data; clear data

for s = 1:length(old_data)
    data(s).ref.resp_sig.imp.v = old_data{s}.waves.RESP.v;
    data(s).ref.resp_sig.imp.fs = old_data{s}.waves.fs;
    data(s).ref.params.hr.t = old_data{s}.numerics.HR.t;
    data(s).ref.params.hr.v = old_data{s}.numerics.HR.v;
    data(s).ref.params.hr.units.t = 'secs since beginning of signal';
    data(s).ref.params.hr.units.v = 'beats per minute';
    data(s).ref.params.pr.t = old_data{s}.numerics.PULSE.t;
    data(s).ref.params.pr.v = old_data{s}.numerics.PULSE.v;
    data(s).ref.params.pr.units.t = 'secs since beginning of signal';
    data(s).ref.params.pr.units.v = 'beats per minute';
    data(s).ref.params.rr.t = old_data{s}.numerics.RESP.t;
    data(s).ref.params.rr.v = old_data{s}.numerics.RESP.v;
    data(s).ref.params.rr.units.t = 'secs since beginning of signal';
    data(s).ref.params.rr.units.v = 'breaths per minute';
    data(s).ref.params.spo2.t = old_data{s}.numerics.SpO2.t;
    data(s).ref.params.spo2.v = old_data{s}.numerics.SpO2.v;
    data(s).ref.params.spo2.units.t = 'secs since beginning of signal';
    data(s).ref.params.spo2.units.v = '%';
    data(s).ekg.v = old_data{s}.waves.II.v;
    data(s).ekg.fs = old_data{s}.waves.fs;
    data(s).ppg.v = old_data{s}.waves.PLETH.v;
    data(s).ppg.fs = old_data{s}.waves.fs;
    data(s).fix.id = old_data{s}.ID;
    data(s).fix.loc = old_data{s}.loc;
    data(s).fix.source.waves = download.source.waves_data_file{s};
    data(s).fix.source.numerics = download.source.numerics_data_file{s};

end

%% Save extracted data
save(up.paths.extracted_data_file, 'data')

end

function [numerics_start, waveforms_start] = calc_start_times(download, pt_stay_no)

% Identify relevant header files
all_waves_header_file = download.saved.overall_header_file{pt_stay_no};
specific_waves_header_file = download.saved.waves_header_file{pt_stay_no};
numerics_header_file = download.saved.numerics_header_file{pt_stay_no};

% Find start time of all waveform data
fid = fopen(all_waves_header_file);
lines = textscan(fid,'%s','Delimiter','\n'); header_line = lines{1,1}{1};
fclose(fid);
colons = strfind(header_line, ':');
timings.all_waves.start_time = header_line(colons(1)-2:colons(2)+6);
slashes = strfind(header_line, '/'); slashes = slashes(slashes>colons(2));
timings.all_waves.start_date = header_line(slashes(1)-2:slashes(2)+4);
timings.all_waves.deb = datenum([timings.all_waves.start_date, ' ', timings.all_waves.start_time], 'dd/mm/yyyy HH:MM:SS.FFF');
clear slashes colons fid lines

% Find start time of the waveform data in this file
fid = fopen(all_waves_header_file);
lines = textscan(fid,'%s','Delimiter','\n'); lines = lines{1,1};
fclose(fid);
% First line is a header line:
lines = lines(2:end);
% Exclude next line if it's a layout line
if contains(lines{1}, 'layout 0')
    lines = lines(2:end);
end
% Exclude last line if it's a location line
if strcmp(lines{end}(1), '#')
    lines = lines(1:end-1);
end
% Identify duration (and cumulative duration) of each waveform file
for line_no = 1 : length(lines)
    space_loc = strfind(lines{line_no}, ' ');
    line_data_duration(line_no,1) = str2double(lines{line_no}(space_loc+1:end));
    line_data_name{line_no,1} = lines{line_no}(1:space_loc-1);
end
clear fid lines line_no space_loc
line_cum_duration = [0;cumsum(line_data_duration)]; % in ms
% Find start time of the waveform file of interest
rel_line = find(strcmp(line_data_name, download.source.waves_data_filename{pt_stay_no}(1:end-4)));
timings.rel_waves.deb = timings.all_waves.deb + (line_cum_duration(rel_line)/(125*60*60*24));

% % Find start time of waveforms data
% fid = fopen(specific_waves_header_file);
% lines = textscan(fid,'%s','Delimiter','\n'); header_line = lines{1,1}{1};
% fclose(fid);
% colons = strfind(header_line, ':');
% timings.rel_waves.start_time = header_line(colons(1)-2:colons(2)+6);
% slashes = strfind(header_line, '/'); slashes = slashes(slashes>colons(2));
% timings.rel_waves.start_date = header_line(slashes(1)-2:slashes(2)+4);
% timings.rel_waves.deb = datenum([timings.rel_waves.start_date, ' ', timings.rel_waves.start_time], 'dd/mm/yyyy HH:MM:SS');
% clear slashes colons fid lines

% Find start time of all numerics data
fid = fopen(numerics_header_file);
lines = textscan(fid,'%s','Delimiter','\n'); header_line = lines{1,1}{1};
fclose(fid);
colons = strfind(header_line, ':');
timings.all_numerics.start_time = header_line(colons(1)-2:colons(2)+6);
slashes = strfind(header_line, '/'); slashes = slashes(slashes>colons(2));
timings.all_numerics.start_date = header_line(slashes(1)-2:slashes(2)+4);
timings.all_numerics.deb = datenum([timings.all_numerics.start_date, ' ', timings.all_numerics.start_time], 'dd/mm/yyyy HH:MM:SS');
clear slashes colons fid lines

numerics_start = timings.all_numerics.deb;
waveforms_start = timings.rel_waves.deb;

end

function headers = extract_file_info(filepath)

%% Obtain info on signals

fileID = fopen([filepath, '.hea']);
header_lines = textscan(fileID,'%s','Delimiter','\n'); header_lines = header_lines{1,1};
fclose all;

% Extract the relevant header information
temp.top_line = textscan(header_lines{1}, '%s', 'delimiter', sprintf(' ')); temp.top_line = temp.top_line{1,1};  % adapted from http://uk.mathworks.com/matlabcentral/newsreader/view_thread/249016
headers.no_sigs = str2double(temp.top_line{2});
headers.fs = str2double(temp.top_line{3});
headers.no_samps = str2double(temp.top_line{4});
rel_line = find(~cellfun(@isempty, strfind(temp.top_line, ':')));
headers.start_time = temp.top_line{rel_line};
clear rel_line

% extract signal names
headers.sig_names = cell(0);
for sig_no = 1 : headers.no_sigs
    temp.curr_line = textscan(header_lines{1+sig_no}, '%s', 'delimiter', sprintf(' ')); temp.curr_line = temp.curr_line{1,1};  % adapted from http://uk.mathworks.com/matlabcentral/newsreader/view_thread/249016
    headers.sig_names = [headers.sig_names; temp.curr_line{end}];
end

% extract location
loc_line = header_lines{end};
headers.loc = loc_line(13:end);

end

function save_data_in_rrest_format(up)

fprintf('\n -- Saving data in RRest toolbox format')

temp = dir([up.paths.individual_files_folder, '*.mat']);
created_files = extractfield(temp, 'name'); clear temp

for subj_el = 1:length(created_files)

    % Load this subject's data
    subj_data = load([up.paths.individual_files_folder, created_files{subj_el}]);

    % Insert fixed params
    data(1,subj_el).fix.id = created_files{subj_el}(10:15);
    data(1,subj_el).fix.data_segment = created_files{subj_el}(17:end);
    data(1,subj_el).fix.loc = subj_data.location;
    data(1,subj_el).fix.ventilation = 'unknown';
    data(1,subj_el).fix.recording_conditions = 'critical care';
    data(1,subj_el).fix.source = subj_data.download_source;

    % insert PPG signal
    rel_col = strcmp(subj_data.header, 'PLETH');
    data(1,subj_el).ppg.v = subj_data.data(:,rel_col);
    data(1,subj_el).ppg.fs = subj_data.fs(rel_col);
    data(1,subj_el).ppg.method = 'measured using clinical monitor';
    % insert EKG signal
    rel_col = strcmp(subj_data.header, 'II');
    data(1,subj_el).ekg.v = subj_data.data(:,rel_col);
    data(1,subj_el).ekg.fs = subj_data.fs(rel_col);
    data(1,subj_el).ekg.method = 'measured using clinical monitor';
    % insert RESP signal
    rel_col = strcmp(subj_data.header, 'RESP');
    data(1,subj_el).ref.resp_sig.imp.v = subj_data.data(:,rel_col);
    data(1,subj_el).ref.resp_sig.imp.fs = subj_data.fs(rel_col);
    data(1,subj_el).ref.resp_sig.imp.method = 'thoracic impedance measured using clinical monitor';
    % insert annotations
    data(1,subj_el).ref.breaths.ann1 = subj_data.ann.resp.breaths.x1;
    data(1,subj_el).ref.breaths.ann2 = subj_data.ann.resp.breaths.x2;
    data(1,subj_el).ref.breaths.method = 'two independent sets of breath annotations using the impedance respiratory signal';

    % insert numerics
    if sum(strcmp(fieldnames(subj_data), 'numerics'))
        data(1,subj_el).ref.params.rr.t = subj_data.numerics.RESP.t;
        data(1,subj_el).ref.params.rr.v = subj_data.numerics.RESP.v;
        data(1,subj_el).ref.params.rr.method = 'derived from thoracic impedance measured using clinical monitor';
        data(1,subj_el).ref.params.rr.units.t = 's';
        data(1,subj_el).ref.params.rr.units.v = 'breaths/min';
        data(1,subj_el).ref.params.hr.t = subj_data.numerics.HR.t;
        data(1,subj_el).ref.params.hr.v = subj_data.numerics.HR.v;
        data(1,subj_el).ref.params.hr.method = 'derived from ecg measured using clinical monitor';
        data(1,subj_el).ref.params.hr.units.t = 's';
        data(1,subj_el).ref.params.hr.units.v = 'beats/min';
        data(1,subj_el).ref.params.pr.t = subj_data.numerics.PULSE.t;
        data(1,subj_el).ref.params.pr.v = subj_data.numerics.PULSE.v;
        data(1,subj_el).ref.params.pr.method = 'derived from ppg measured using clinical monitor';
        data(1,subj_el).ref.params.pr.units.t = 's';
        data(1,subj_el).ref.params.pr.units.v = 'beats/min';
        data(1,subj_el).ref.params.spo2.t = subj_data.numerics.SpO2.t;
        data(1,subj_el).ref.params.spo2.v = subj_data.numerics.SpO2.v;
        data(1,subj_el).ref.params.spo2.method = 'measured using clinical monitor';
        data(1,subj_el).ref.params.spo2.units.t = 's';
        data(1,subj_el).ref.params.spo2.units.v = '%';
    end

end

% Save to file
savepath = [up.paths.data_root, up.dataset_name, '_data'];
save(savepath, 'data')

fprintf(['\n --- Dataset saved to: ' savepath])

end

function save_data_in_csv_format(up)

fprintf('\n -- Converting to CSV format');

% lower case of dataset name to keep to usual format of PhysioNet files:
dataset_name = up.dataset_name;

% Create folder to save data in
data_type = 'csv';
save_folder = [up.paths.data_root, dataset_name, '_', data_type, filesep];
if ~exist(save_folder, 'dir')
    mkdir(save_folder)
end
cd(save_folder)

% Load data
load(up.paths.extracted_data_file);

% cycle through each subject
no_subjs = length(data);
file_names = cell(no_subjs,1);
for subj_no = 1 : no_subjs

    % setup command
    file_names{subj_no} = [up.dataset_name, '_', num2str(subj_no, '%02.f')];
    fs = data{subj_no}.waves.fs;
    sig_names = fieldnames(data{subj_no}.waves); sig_names = sig_names(~strcmp(sig_names, 'fs'));
    num_names = fieldnames(data{subj_no}.numerics); num_names = num_names(~strcmp(num_names, 'fs'));

    % create a file of fixed variables for this subject
    curr_filename = [file_names{subj_no}, '_Fix.txt'];
    fid = fopen(curr_filename, 'wt');

    signals_line = 'Time [s],';
    data_mat = data{subj_no}.waves.II.t;
    for sig_no = 1 : length(sig_names)
        signals_line = [signals_line, ' ', sig_names{sig_no}];
        eval(['curr_sig = data{subj_no}.waves.' sig_names{sig_no} ';'])
        data_mat = [data_mat, curr_sig.v];
        % check whether time vector is consistent
        if ~isequal(curr_sig.t, data{subj_no}.waves.II.t)
            error('Check this')
        end

        if sig_no < length(sig_names)
            signals_line = [signals_line, ';'];
        end
    end
    clear sig_no curr_sig

    numerics_line = 'Time [s],';
    numerics_mat = data{subj_no}.numerics.HR.t;
    numerics_fs = data{subj_no}.numerics.fs;
    for num_no = 1 : length(num_names)
        numerics_line = [numerics_line, ' ', num_names{num_no}];
        eval(['curr_num = data{subj_no}.numerics.' num_names{num_no} ';'])
        numerics_mat = [numerics_mat, curr_num.v];
        % check whether time vector is consistent
        if ~isequal(curr_num.t, data{subj_no}.numerics.HR.t)
            error('Check this')
        end

        if num_no < length(num_names)
            numerics_line = [numerics_line, ';'];
        end
    end
    clear num_no curr_num

    % create a text file containing the fixed variables for this subject
    text_to_write = [file_names{subj_no}, up.csv.new_line, ...
        'Signals:' signals_line(10:end), up.csv.new_line, ...
        'Signals sampling frequency: ' num2str(fs) ' Hz', up.csv.new_line, ...
        'Numerics:' numerics_line(10:end), up.csv.new_line, ...
        'Numerics sampling frequency: ' num2str(numerics_fs) ' Hz', up.csv.new_line, ...
        'Age: ' num2str(data{subj_no}.age), up.csv.new_line, ...
        'Gender: ' data{subj_no}.sex, up.csv.new_line, ...
        'Location: ' data{subj_no}.loc, up.csv.new_line, ...
        'MIMIC II matched wdb ID: ' data{subj_no}.ID, up.csv.new_line, ...
        'MIMIC II data segment: ' data{subj_no}.data_segment, up.csv.new_line, ...
        'Original overall header file: ' data{subj_no}.download_source.overall_header_file, up.csv.new_line, ...
        'Original signals header file: ' data{subj_no}.download_source.waves_header_file, up.csv.new_line, ...
        'Original signals data file: ' data{subj_no}.download_source.waves_data_file, up.csv.new_line, ...
        'Original numerics header file: ' data{subj_no}.download_source.numerics_header_file, up.csv.new_line, ...
        'Original numerics data file: ' data{subj_no}.download_source.numerics_data_file, up.csv.new_line, ...
        ];
    fprintf(fid, text_to_write);
    fclose(fid);
    clear text_to_write fid

    % create a csv file of signals for this subject
    curr_filename = [file_names{subj_no}, '_Signals.csv'];
    fid = fopen(curr_filename,'w');
    signals_line = strrep(signals_line, ';', ',');
    fprintf(fid,[signals_line, up.csv.new_line]); clear signals_line
    fclose(fid); clear fid
    dlmwrite(curr_filename, data_mat, '-append'); clear data_mat curr_filename

    % create a csv file of numerics for this subject
    curr_filename = [file_names{subj_no}, '_Numerics.csv'];
    fid = fopen(curr_filename,'w');
    numerics_line = strrep(numerics_line, ';', ',');
    fprintf(fid,[numerics_line, up.csv.new_line]); clear numerics_line
    fclose(fid); clear fid
    dlmwrite(curr_filename, numerics_mat, '-append'); clear numerics_mat curr_filename

    % create a csv file of annotations for this subject
    curr_filename = [file_names{subj_no}, '_Breaths.csv'];
    fid = fopen(curr_filename,'w');
    header_line = 'breaths ann1 [signal sample no], breaths ann2 [signal sample no]';
    fprintf(fid,[header_line, up.csv.new_line]); clear header_line
    fclose(fid); clear fid
    length_diff = length(data{subj_no}.ann.resp.breaths.x1(:)) - length(data{subj_no}.ann.resp.breaths.x2(:));
    if length_diff > 0
        data{subj_no}.ann.resp.breaths.x2 = [data{subj_no}.ann.resp.breaths.x2, nan(1, length_diff)];
    elseif length_diff < 0
        data{subj_no}.ann.resp.breaths.x1 = [data{subj_no}.ann.resp.breaths.x1, nan(1, abs(length_diff))];
    end
    ann_mat = [data{subj_no}.ann.resp.breaths.x1(:), data{subj_no}.ann.resp.breaths.x2(:)];
    dlmwrite(curr_filename, ann_mat, '-append'); clear ann_mat curr_filename

    clear fs sig_names signals_line curr_filename numerics_line

end
clear subj_no

end

function save_data_in_wfdb_format(up)

fprintf('\n -- Converting to WFDB format');

% Create folder to save data in
data_type = 'wfdb';
save_folder = [up.paths.data_root, up.dataset_name, '_', data_type, filesep];
if ~exist(save_folder, 'dir')
    mkdir(save_folder)
end
cd(save_folder)

% Load data
load(up.paths.extracted_data_file);

% cycle through each subject
no_subjs = length(data);
file_names = cell(no_subjs,1);
for subj_no = 1 : no_subjs

    % setup command
    file_names{subj_no} = [up.dataset_name, sprintf('%.2d', subj_no)];
    fs = data{subj_no}.waves.fs;
    sig_names = fieldnames(data{subj_no}.waves); sig_names = sig_names(~strcmp(sig_names, 'fs'));
    num_names = fieldnames(data{subj_no}.numerics); num_names = num_names(~strcmp(num_names, 'fs'));

    % Extract signals data
    data_mat = [];
    for sig_no = 1 : length(sig_names)
        eval(['curr_sig = data{subj_no}.waves.' sig_names{sig_no} ';'])
        data_mat = [data_mat, curr_sig.v];
        % check whether time vector is consistent
        if ~isequal(curr_sig.t, data{subj_no}.waves.II.t)
            error('Check this')
        end
    end
    clear sig_no curr_sig

    % Create signals file
    descrip = ['<age>: ' num2str(data{subj_no}.age) ' <sex>: ' data{subj_no}.sex ' <location>: ' data{subj_no}.loc ' <source>: https://physionet.org/physiobank/database/mimic2wdb/matched/' data{subj_no}.ID '/ <data segment>: ' data{subj_no}.data_segment '<modifications>: this is not a verbatim copy of an original file. Please see the accompanying LICENSE.txt file for further details.'];
    units = '';
    for sig_no = 1 : length(sig_names)
        signal_names(sig_no,1) = {[sig_names{sig_no} ', ' ]};
        rel_unit = up.units.units{strcmp(up.units.signals, sig_names{sig_no})};
        units = [units, rel_unit, '/']; clear rel_unit
    end
    clear sig_no
    units = units(1:end-1);
    % convert this subject's data into WFDB format
    mat2wfdb(data_mat, file_names{subj_no}, fs, [], units, descrip, [], signal_names);
    clear data_mat fs units descrip

    % Extract numerics data
    numerics_mat = [];
    for num_no = 1 : length(num_names)
        eval(['curr_num = data{subj_no}.numerics.' num_names{num_no} ';'])
        numerics_mat = [numerics_mat, curr_num.v];
        % check whether time vector is consistent
        if ~isequal(curr_num.t, data{subj_no}.numerics.HR.t)
            error('Check this')
        end
    end
    clear num_no curr_num

    % Create numerics file
    descrip = ['<age>: ' num2str(data{subj_no}.age) ' <sex>: ' data{subj_no}.sex ' <location>: ' data{subj_no}.loc ' <source>: ' data{subj_no}.download_source.numerics_data_file ' <modifications>: this is not a verbatim copy of an original file. Please see the accompanying LICENSE.txt file for further details.'];
    units = '';
    for num_no = 1 : length(num_names)
        numeric_names(num_no,1) = {[num_names{num_no} ', ' ]};
        rel_unit = up.units.num_units{strcmp(up.units.numerics, num_names{num_no})};
        units = [units, rel_unit, '/']; clear rel_unit
    end
    clear sig_no
    units = units(1:end-1);
    % convert this subject's data into WFDB format
    fs = data{subj_no}.numerics.fs;
    mat2wfdb(numerics_mat, [file_names{subj_no}, 'n'], fs, [], units, descrip, [], numeric_names);
    clear numerics_mat fs units descrip

    % create annotation files
    ann1 = data{subj_no}.ann.resp.breaths.x1;
    ann2 = data{subj_no}.ann.resp.breaths.x2;
    ann = [ann1(:); ann2(:)]; % col vector of samples of anns
    chan = find(~cellfun(@isempty, strfind(sig_names, 'RESP')))*ones(size(ann))-1;
    num = zeros(size(ann));
    comments = [repmat({'ann1'}, [length(ann1),1]); repmat({'ann2'}, [length(ann2),1])]; % ann1 or ann2 as char
    type = repmat('"', size(ann));
    subtype = zeros(size(num));
    wrann(file_names{subj_no},'breath',ann,type,subtype,chan,num,comments);

    clear sig_names ann1 ann2 ann chan num type subtype

end
clear no_subjs subj_no

% create a RECORDS file
file_name = 'RECORDS';
file_path = [save_folder, file_name];
fid = fopen(file_path, 'w');
formatSpec = '%s\r\n';
for line_no = 1 : length(file_names)
    fprintf(fid, formatSpec, file_names{line_no});
end

% create an ANNOTATORS file
file_name = 'ANNOTATORS';
file_path = [save_folder, file_name];
fid = fopen(file_path, 'w');
fprintf(fid, 'ann1\tbreath annotations by annotator 1\n');
fprintf(fid, 'ann2\tbreath annotations by annotator 2');

% create a DBS file
file_name = 'DBS';
file_path = [save_folder, file_name];
fid = fopen(file_path, 'w');
fprintf(fid, [up.dataset_name, '\t', 'BIDMC ECG, PPG and Respiration Dataset']);

fclose all;

end

function create_readme_and_license_files(up)

fprintf('\n -- Saving ReadMe and License files')

%% License file
licence_file_path = [up.paths.data_root, 'LICENSE'];
fid = fopen(licence_file_path, 'wt');
fprintf(fid, ['The mimic_imp_sqi_data file produced by this script contains modified versions of PhysioNet files. The modified files are distributed under the terms for distribution of modified data files listed at: https://archive.physionet.org/copying.shtml#data .',...
    '\n\n The following modifications have been made:', ...
    '\n\n - Waveform, numerics and fixed data (such as age) were downloaded from the PhysioNet MIMIC III matched waveform database. 1-hour segments of data were extracted for each of the 100 ICU stays included in the dataset.',...
    '\n\n - The data were imported into Matlab. In addition, manually annotated breaths are provided for some 32 second windows of most subjects'' impedance pneumography signals, and estimated respiratory rates (calculated from these manual annotations) are provided.', ...
    '\n\n Further details of the dataset are provided in the accompanying ReadMe file.']);
fclose all;

readme_file_path = [up.paths.data_root, 'README'];
fid = fopen(readme_file_path, 'wt');

%% ReadMe file

readme_text = [...
'Impedance Pneumography Signal Quality Index Dataset',...
'\n=============',...
'\n\nTable of Contents',...
'\n=================',...
'\n\n  * [The Dataset](#the-dataset)',...
'\n  * [Revisions](#revisions)',...
'\n  * [Description](#description)',...
'\n\n# The Dataset',...
'\n\nThe Impedance Pneumography Signal Quality Index Dataset was first reported in the following publication:',...
'\n\nCharlton, P.H. et al. An impedance pneumography signal quality index: design, assessment and application to respiratory rate monitoring, [under review], 2020',...
'\n\nIn this publication it was used to assess the performance of a signal quality index for the impedance pneumography signal. The dataset was extracted from the much larger [MIMIC III matched waveform Database](https://archive.physionet.org/physiobank/database/mimic3wdb/matched/), which was acquired from critically-ill patients during hospital care at the Beth Israel Deaconess Medical Centre (Boston, MA, USA). The 100 recordings within the dataset, each of 1-hour duration, each contain:',...
'\n  * Physiological signals: the impedance pneumography respiratory signal (ImP), photoplethysmogram (PPG) and electrocardiogram (ECG). These are sampled at 125 Hz.',...
'\n  * Physiological parameters, such as the heart rate (HR), pulse rate (PR), respiratory rate (RR), and blood oxygen saturation level (SpO2). These are sampled at 1 Hz.',...
'\n  * Fixed parameters, such as age, sex, ward and patient ID',...
'\n  * Manual annotations of breaths, and RRs calculated from these: A single annotator manually annotated individual breaths in approximately five 32 sec segments of each recording using the impedance respiratory signal.',...
'\n\nThe dataset is distributed in Matlab (r) format.',...
'\n\nFor more information about the dataset, please contact author at:',...
'\npc657@medschl.cam.ac.uk .',...
'\n\n# Revisions',...
'\n\nR1: 	2017-Aug-06 		initial release',...
'\n\n## Description',...
'\n### Matlab (r) Format',...
'\nThe *mimic_imp_sqi_data.mat* file contains the following subset of the dataset in a single Matlab (r) variable named *data*. The following are provided for each of the 100 recordings:',...
'\n* *ekg*:   Lead II ECG signal. Each signal is provided in a structure, where the *v* field denotes the signal values, and *fs* is the sampling frequency.',...
'\n* *ppg*:   Photoplethysmogram signal',...
'\n* *ref.resp_sig.imp*:  Impedance respiratory signal',...
'\n* *ref.ann*:  Manual annotations of breaths provided by a single annotator. A vector of sample numbers is provided (*breath_els*), which correspond to the signal sample numbers. Breaths were only annotated in certain segments of the signal, rather than the whole signal: the start and end samples of each segment in which breaths were annotated are provided in *win_start_els* and *win_end_els*. RRs calculated from the manual annotations are provided (*rr*, one per window).',...
'\n* *ref.params*:  Physiological parameters: *rr* (respiratory rate, derived by the monitor from the impedance signal, breaths per minute), *hr* (heart rate, derived from the ECG, beats per minute), *pr* (pulse rate, derived from the PPG, beats per minute), *spo2* (blood oxygen saturation level, %%).',...
'\n* *fix*: A structure of fixed variables: *id* (the MIMIC III matched waveform database subject ID and recording identifier), *loc* (the ward location), and download source (specified for both the waves and numerics).',...
'\n\nFurther details on the contents of each file are provided [here](https://archive.physionet.org/tutorials/creating-records.shtml).'];

fprintf(fid, readme_text);
fclose all;

end

function add_manual_annotations(up)

add_manual_annotations = 1;

if ~add_manual_annotations, return, else clear add_manual_annotations, end

fprintf('\n -- Adding in Manual Annotations')

% load extracted data file
load(up.paths.extracted_data_file)

% specify folder containing manual annotations files
ann_folder = '/Users/petercharlton/Documents/Data/Imp_SQI_mimic_dataset_test/2019_annotations/';
manual_annotator = 'PC';
temp = dir([ann_folder, '*', manual_annotator, '_an.mat']);
ann_files = extractfield(temp, 'name'); clear temp

do_check = 0;
% import annotations for each manually annotated subject
for subj_no = 1 : length(data)

    % create fields in which to store annotations
    data(subj_no).ref.ann.win_start_els = [];
    data(subj_no).ref.ann.win_end_els = [];
    data(subj_no).ref.ann.breath_els = [];
    data(subj_no).ref.ann.rr = [];

    % identify candidate annotations file for this subject
    annotations_filename = ['ImP_SQI' num2str(subj_no) '_' manual_annotator '_an.mat'];
    annotations_filepath = [ann_folder, annotations_filename];

    % skip if no annotations for this subject
    if ~exist(annotations_filepath)
        continue
    end

    % load annotations (which are stored in 'pk_anns' and 'qual')
    temp = load(annotations_filepath);

    % identify windows
    imp.v = data(subj_no).ref.resp_sig.imp.v; imp.fs = 125;
    imp.t = [0:(length(imp.v)-1)]/imp.fs;
    up.paramSet.winLeng = 32;     % the duration of each window in secs (no overlap)
    win_starts = imp.t(1):up.paramSet.winLeng:imp.t(end);
    win_ends = win_starts+up.paramSet.winLeng;

    for win_no = 1 : length(win_starts)

        % Identify annotations within this window
        if sum(strcmp(fieldnames(temp), 'pk_anns'))
            win_start = win_starts(win_no);
            win_end = win_ends(win_no);
            rel_ann_els = temp.pk_anns.t >= win_start & temp.pk_anns.t <= win_end;
            if sum(rel_ann_els) > 1
                first_breath = min(temp.pk_anns.t(rel_ann_els));
                last_breath = max(temp.pk_anns.t(rel_ann_els));
                no_breaths = sum(rel_ann_els)-1;
                rr_ann = 60/((last_breath-first_breath)/no_breaths);

                data(subj_no).ref.ann.win_start_els = [data(subj_no).ref.ann.win_start_els; win_start*125+1];
                data(subj_no).ref.ann.win_end_els = [data(subj_no).ref.ann.win_end_els; win_end*125+1];
                data(subj_no).ref.ann.breath_els = [data(subj_no).ref.ann.breath_els; round(temp.pk_anns.t(rel_ann_els)*125+1)];
                data(subj_no).ref.ann.rr = [data(subj_no).ref.ann.rr; rr_ann];

                % to check:
                if do_check
                    sig_els = (win_starts(win_no)*125+1):(win_ends(win_no)*125+1); plot(imp.t(sig_els), imp.v(sig_els)), hold on, rel_breath_els = data(subj_no).ref.ann.breath_els(data(subj_no).ref.ann.breath_els >= sig_els(1) & data(subj_no).ref.ann.breath_els <= sig_els(end)); plot(imp.t(rel_breath_els), imp.v(rel_breath_els), '*r')
                    close all
                end

                clear rr_ann no_breaths first_breath last_breath
            end
            clear win_start win_end rel_ann_els
        end

    end
    % to check
    if do_check
        plot(imp.t,imp.v), hold on, plot(imp.t(data(subj_no).ref.ann.breath_els), imp.v(data(subj_no).ref.ann.breath_els), '*r')
        close all
    end

    clear win_no temp imp win_starts win_ends

end

% Save new data file with annotations added
save(up.paths.extracted_data_file, 'data')

end

author: peterhcharlton

repo: RRest


text:
function n = norm(T)
%NORM Frobenius norm of a tenmat.
%
%   NORM(X) returns the Frobenius norm of a tenmat.
%
%   See also TENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

v = reshape(T.data, numel(T.data), 1);
n = norm(v);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%% 2dmatrix hard_threshold(2dmatrix)
% S - sparse 2dmatrix
% O - outliers 2dmatrix
%
function O = hard_threshold(S)
  displog('Applying hard threshold...');

  % beta = 0.5*(3*std(S(:))/20)^2; % min beta, lower bound: suppose SNR <= 20
  beta = 0.5*(std(S(:)))^2; % begin beta, start from a big value

  % direct hard thresholding if no smoothness
  O = double(0.5*S.^2 > beta);
end

repo: lrslibrary


text:
function varargout = GUIEditSegyHeader(varargin)
% GUIEditSegyHeader : GUI for editing the SGY header
%
% Call :
%   SH=GUIEditSegyHeader(SH);
%
%
% Example :
%   [Data,STH,SH]=ReadSegy('841_m.sgy');
%   SH=GUIEditSegyHeader(SH);
%
%

% GUIEDITSEGYHEADER Application M-file for GUIEditSegyTraceHeader.fig
%    FIG = GUIEDITSEGYHEADER launch GUIEditSegyTraceHeader GUI.
%    GUIEDITSEGYHEADER('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.5 03-Mar-2003 15:33:07

if (nargin == 0)|(isstruct(varargin{1}))  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');

	% Use system color scheme for figure:
	set(fig,'Color',get(0,'defaultUicontrolBackgroundColor'));

	% Generate a structure of handles to pass to callbacks, and store it.
	handles = guihandles(fig);
	guidata(fig, handles);

    if nargin>0,
       data = guihandles(fig); % initialize it to contain handles
       data.SegyHeader=varargin{1};
       data.SegyHeaderOrig=varargin{1};
       guidata(fig,data);

       set(fig,'HandleVisibility','On')
       GUIEditSegyHeader('actionSetHeader',fig,handles)
       GUIEditSegyHeader('actionUpdateHeader',fig,handles)
       set(fig,'HandleVisibility','CallBack')

    end

	% Wait for callbacks to run and window to be dismissed:
	uiwait(fig);

	if nargout > 0
	  data=guidata(fig);
	  if isstruct(data.SegyHeader)
	    varargout{1}=data.SegyHeader;
	  else
	    varargout{1} = fig;
	  end
    end

    close(handles.figure1);
	set(handles.figure1,'Visible','off')

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end

%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and
%| sets objects' callback properties to call them through the FEVAL
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the
%| callback type separated by '_', eSweepFrequencyStart.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, eSweepFrequencyStart.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.

% --------------------------------------------------------------------

% --------------------------------------------------------------------
function varargout = figure1_ResizeFcn(h, eventdata, handles, varargin)

% --------------------------------------------------------------------
function varargout = pbDone_Callback(h, eventdata, handles, varargin)
uiresume;

% --------------------------------------------------------------------
function varargout = pbCancel_Callback(h, eventdata, handles, varargin)
data=guidata(h);
data.SegyHeader=data.SegyHeaderOrig;
guidata(h,data);
uiresume;

clear name;
function actionSetHeader(h,handles)
data=guidata(h);
SH=data.SegyHeader;

% Revision
revision=SH.SegyFormatRevisionNumber;
r=1; % Default old revisoin
if revision>50, r=2;end
if revision==0, r=1;end
for i=1:length(SH.Rev), Revision{i}=SH.Rev(i).name;end
set(handles.popSegyFormatRevision,'String',Revision);

clear name;
% DataSampleFormat
DSF=SH.Rev(r).DataSampleFormat;
for i=1:length(DSF), name{i}=DSF(i).name;end
set(handles.popDataSampleFormat,'String',name);

clear name;
% TraceSorting
TS=SH.Rev(r).TraceSorting;
for i=1:length(TS), name{i}=TS(i).name;end
set(handles.popTraceSorting,'String',name);

clear name;
% Sweep Type
ST=SH.Rev(r).SweepType;
for i=1:length(ST), name{i}=ST(i).name;end
set(handles.popSweepType,'String',name);

clear name;
% BinaryGain
BG=SH.Rev(r).BinaryGain;
for i=1:length(BG), name{i}=BG(i).name;end
set(handles.popBinaryGain,'String',name);

clear name;
% TaperType
TT=SH.Rev(r).TaperType;
for i=1:length(TT), name{i}=TT(i).name;end
set(handles.popTaperType,'String',name);

clear name;
% CorrelatedDataTraces
CD=SH.Rev(r).CorrelatedDataTraces;
for i=1:length(CD), name{i}=CD(i).name;end
set(handles.popCorrelatedDataTraces,'String',name);

clear name;
% AmplitudeRecoveryMethod
ARM=SH.Rev(r).AmplitudeRecoveryMethod;
for i=1:length(ARM), name{i}=ARM(i).name;end
set(handles.popAmplitudeRecoveryMethod,'String',name);

clear name;
% MeasurementSystem
M=SH.Rev(r).MeasurementSystem;
for i=1:length(M), name{i}=M(i).name;end
set(handles.popMeasurementSystem,'String',name);

clear name;
% ImpulseSignalPolarity
ISR=SH.Rev(r).ImpulseSignalPolarity;
for i=1:length(ISR), name{i}=ISR(i).name;end
set(handles.popImpulseSignalPolarity,'String',name);

clear name;
% popVibratoryPolarityCode
VPC=SH.Rev(r).VibratoryPolarityCode;
for i=1:length(VPC), name{i}=VPC(i).name;end
set(handles.popVibratoryPolarityCode,'String',name);

clear name;
% popFixedLengthTraceFlag
FLT=SH.Rev(r).FixedLengthTraceFlag;
for i=1:length(FLT), name{i}=FLT(i).name;end
set(handles.popFixedLengthTraceFlag,'String',name);

% EDIT BOXES
set(handles.eJob,'String',SH.Job);
set(handles.eLine,'String',SH.Line);
set(handles.eReel,'String',SH.Reel);
set(handles.eDataTracePerEnsemble,'String',SH.DataTracePerEnsemble);
set(handles.eAuxiliaryTracePerEnsemble,'String',SH.AuxiliaryTracePerEnsemble);
set(handles.edt,'String',SH.dt);
set(handles.edtOrig,'String',SH.dtOrig);
set(handles.ens,'String',SH.ns);
set(handles.ensOrig,'String',SH.nsOrig);
set(handles.eEnsembleFold,'String',SH.EnsembleFold);
set(handles.eVerticalSumCode,'String',SH.VerticalSumCode);

set(handles.eSweepFrequencyStart,'String',SH.SweepFrequencyStart);
set(handles.eSweepFrequencyEnd,'String',SH.SweepFrequencyEnd);
set(handles.eSweepLength,'String',SH.SweepLength);
set(handles.eSweepChannel,'String',SH.SweepChannel);
set(handles.eSweepTaperlengthStart,'String',SH.SweepTaperlengthStart);
set(handles.eSweepTaperLengthEnd,'String',SH.SweepTaperLengthEnd);

%---------------------------------------------------------
function r=actionGetrevision(h);
data=guidata(h);
SH=data.SegyHeader;
% Revision
revision=SH.SegyFormatRevisionNumber;
r=1; % Default old revisoin
if revision>50, r=2;end
if revision==0, r=1;end

%---------------------------------------------------------
function actionUpdateHeader(h,handles)
data=guidata(h);

SH=data.SegyHeader;

% Revision
revision=SH.SegyFormatRevisionNumber;
r=1; % Default old revisoin
if revision>50, r=2;end
if revision==0, r=1;end
set(handles.popSegyFormatRevision,'value',r);

% DataSampleFormat
DSF=SH.Rev(r).DataSampleFormat;
if length(get(handles.popDataSampleFormat,'String'))<SH.DataSampleFormat,
    SH.DataSampleFormat=1;
    disp([mfilename,' : Using DataSampleFormat ',num2str(SH.DataSampleFormat)])
end
set(handles.popDataSampleFormat,'value',SH.DataSampleFormat);

% MAYBE THE NEXT LINES COULD BE USED FOR MORE HEADERS
% TraceSorting
TraceSortingSelect=1;
for i=1:length(SH.Rev(r).TraceSorting)
    if SH.TraceSorting==SH.Rev(r).TraceSorting(i).value;
        TraceSortingSelect=i;
        SH.TraceSorting;
    end
end

if TraceSortingSelect>length(get(handles.popTraceSorting,'String'));
    TraceSortingSelect=1;
end
set(handles.popTraceSorting,'value',TraceSortingSelect);

SweepTypeSelect=1;
for i=1:length(SH.Rev(r).SweepType)
    if SH.SweepType==SH.Rev(r).SweepType(i).value;
        SweepTypeSelect=i;
        SH.SweepType;
    end
end
if SweepTypeSelect>length(get(handles.popSweepType,'String'));
    SweepTypeSelect=1;
end
set(handles.popSweepType,'value',SweepTypeSelect);

data.SegyHeader=SH;
guidata(h,data)

% --------------------------------------------------------------------
function varargout = popSegyFormatRevision_Callback(h, eventdata, handles, varargin)
data=guidata(h);
sel=get(h,'value');
data.SegyHeader.SegyFormatRevisionNumber =  data.SegyHeader.Rev(sel).SegyFormatRevisionNumber;

guidata(h,data);
actionSetHeader(h,handles)
actionUpdateHeader(h,handles)

% --------------------------------------------------------------------
function varargout = popDataSampleFormat_Callback(h, eventdata, handles, varargin)
data=guidata(h);
sel=get(h,'value');
data.SegyHeader.DataSampleFormat=get(h,'value');
guidata(h,data);

% --------------------------------------------------------------------
function varargout = popTraceSorting_Callback(h, eventdata, handles, varargin)
data=guidata(h);
TraceSortingSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.TraceSorting=data.SegyHeader.Rev(r).TraceSorting(TraceSortingSelect).value;
guidata(h,data);

function popSweepType_Callback(h, eventdata, handles)
data=guidata(h);
SweepTypeSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.SweepType=data.SegyHeader.Rev(r).SweepType(SweepTypeSelect).value;
guidata(h,data);

function popTaperType_Callback(h, eventdata, handles)
data=guidata(h);
TaperTypeSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.TaperType=data.SegyHeader.Rev(r).TaperType(TaperTypeSelect).value;
guidata(h,data);
actionUpdateHeader(h,handles)

function popCorrelatedDataTraces_Callback(h, eventdata, handles)
data=guidata(h);
CorrelatedDataTracesSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.CorrelatedDataTraces=data.SegyHeader.Rev(r).CorrelatedDataTraces(CorrelatedDataTracesSelect).value;
guidata(h,data);

function popBinaryGain_Callback(h, eventdata, handles)
data=guidata(h);
BinaryGainSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.BinaryGain=data.SegyHeader.Rev(r).BinaryGain(BinaryGainSelect).value;
guidata(h,data);

function popAmplitudeRecoveryMethod_Callback(h, eventdata, handles)
data=guidata(h);
AmplitudeRecoveryMethodSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.AmplitudeRecoveryMethod=data.SegyHeader.Rev(r).AmplitudeRecoveryMethod(AmplitudeRecoveryMethodSelect).value;
guidata(h,data);

function varargout = popMeasurementSystem_Callback(h, eventdata, handles, varargin)
data=guidata(h);
MSSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.MeasurementSystem=data.SegyHeader.Rev(r).MeasurementSystem(MSSelect).value;
guidata(h,data);

function popImpulseSignalPolarity_Callback(h, eventdata, handles)
data=guidata(h);
ImpulseSignalPolaritySelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.ImpulseSignalPolarity=data.SegyHeader.Rev(r).ImpulseSignalPolarity(ImpulseSignalPolaritySelect).value;
guidata(h,data);

function popVibratoryPolarityCode_Callback(h, eventdata, handles)
data=guidata(h);
VibratoryPolarityCodeSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.VibratoryPolarityCode=data.SegyHeader.Rev(r).VibratoryPolarityCode(VibratoryPolarityCodeSelect).value;
guidata(h,data);

% --- Executes on selection change in popFixedLengthTraceFlag.
function popFixedLengthTraceFlag_Callback(h, eventdata, handles)
data=guidata(h);
FixedLengthTraceFlagSelect=get(h,'value');
r=actionGetrevision(h);
data.SegyHeader.FixedLengthTraceFlag=data.SegyHeader.Rev(r).FixedLengthTraceFlag(FixedLengthTraceFlagSelect).value;
guidata(h,data);

% --- Executes on button press in pbEditTextualFileHeader.
function pbEditTextualFileHeader_Callback(h, eventdata, handles)

data=guidata(h);
data.SegyHeader=GUIEditTextualFileHeader(data.SegyHeader);
guidata(h,data);

function eJob_Callback(hObject, eventdata, handles)

function eLine_Callback(hObject, eventdata, handles)

function eReel_Callback(hObject, eventdata, handles)

function eDataTracePerEnsemble_Callback(hObject, eventdata, handles)

function eAuxiliaryTracePerEnsemble_Callback(hObject, eventdata, handles)

function edt_Callback(hObject, eventdata, handles)

function edtOrig_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function ens_CreateFcn(hObject, eventdata, handles)

% --- Executes on button press in ens.
function ens_Callback(hObject, eventdata, handles)

% --- Executes on button press in ensOrig.
function ensOrig_Callback(hObject, eventdata, handles)

% --- Executes on button press in eEnsembleFold.
function eEnsembleFold_Callback(hObject, eventdata, handles)

% --- Executes on button press in eVerticalSumCode.
function eVerticalSumCode_Callback(hObject, eventdata, handles)

% --- Executes on button press in eSweepFrequencyStart.
function eSweepFrequencyStart_Callback(hObject, eventdata, handles)

function eSweepFrequencyEnd_Callback(hObject, eventdata, handles)

function eSweepLength_Callback(hObject, eventdata, handles)

function eSweepChannel_Callback(hObject, eventdata, handles)

function eSweepTaperlengthStart_Callback(hObject, eventdata, handles)

function eSweepTaperLengthEnd_Callback(hObject, eventdata, handles)

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
%% Grain Exchange Symmetry
%
%%
% %%
% MTEX allows to identify antipodal directions to model axes and to
% identify misorientations with opposite rotational angle. The later is
% required when working with misorientations between grains of the same
% phase and the order of the grains is arbitrary.

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%% KUKA sunrise toolbox example.
% moving end-effector of the robot on an ellipse.

% Copy right: Mohammad SAFEEA
% 10th-Nov-2017
clear all;
close all;
clc;
% add the path of the KST to the working folder of Matlab
addpath getTheKSTDirectory(pwd)
% Initial configuration
jPos={0., pi / 180 * 30, 0, -pi / 180 * 60, 0,...
                        pi / 180 * 90, 0};
% Start the KST, move robot to initial configuration, start directServo function
[tKuka,flag]=connectToKuka(jPos );
% return;
if flag==false
    fprintf('Can not connect to KST \n');
    fprintf('Program terminated \n');
    return;
end

    % move a little bit back on the X direction
    deltaX=-60;deltaY=0;deltaZ=0.;
    Pos{1}=deltaX;
    Pos{2}=deltaY;
    Pos{3}=deltaZ;
    vel=50;
    movePTPLineEefRelBase( tKuka , Pos, vel);
    % put the pen on the level of the page
    deltaX=0;deltaY=0;deltaZ=-85.;
    Pos{1}=deltaX;
    Pos{2}=deltaY;
    Pos{3}=deltaZ;
    vel=50;
    movePTPLineEefRelBase( tKuka , Pos, vel);
    pause(1);
    % get joints angles of robot
    jPos  = getJointsPos( tKuka );

    % start the direct servo
     realTime_startDirectServoJoints(tKuka);

    % calculate current position of flange point of the robot
    qs=zeros(7,1);
for i=1:7
    qs(i)=jPos{i};
end
        TefTool=eye(4);
        T0=directKinematics(qs,TefTool); % EEF frame transformation matrix
        p0=T0(1:3,4);
        Tt=T0;

%% Define the ellipse, dimentsions are in (meter)
p=p0; % this is the starting point of the ellipse
c=p+([0; 50; 0]/1000); % this is the center of the ellipse
dir=[1 0 0]; % the direction vector of the (a) axis of the ellipse
ratio=0.5; % the radious ratio (a/b) of the ellipse
%% Claulcate ellipse parameters
[R,theta0,a,b,c,errorFlag]=getEllipseParameters(p,c,ratio,dir);
if errorFlag==true % this happens when the start point is on the (dir) vector direction.
    % in such case the ellipse plane can not be specified in space.
    % Stop the realtime motion
    realTime_stopDirectServoJoints(tKuka);
    % Turn off the server
    net_turnOffServer( tKuka );
    return;
end
%% Draw the length of the ellipse arc as a function of the parametric angle (theta)
theta=2*pi;

theta1=theta0+theta;

[ thetaVec,sVec ] = getEllipseLengthVector( a,b, theta0,theta1 );
sizesVec=max(size(sVec));
L=sVec(end); % L is the length of the curve

%% Calculate the times:
velocity=40/1000; % stable velcotiy, m/sec;
accel=25/1000; % linear acceleration, m/sec2
[t0,t1,t2]=calculateInterpolationTimes(L,velocity,accel);

%% Joint space control

    [Ttemp,J]=directKinematics(qs,TefTool);
    vec=Ttemp(1:3,4);

% dls solver parameters
        numberOfIterationForSolver=100;
        lambda=0.5;
        TefTool=eye(4);

    sCoordinate=0;
    interpolationCounter=1;
    thetaVar=theta0;

% start of time
        dateVector0=datevec(now);
        time0=dateVector0(6)+dateVector0(5)*60+dateVector0(4)*60*60; % calculate time at this instant
while true

    % Calculate the elapsed time
        dateVector=datevec(now);
        timeNow=dateVector(6)+dateVector(5)*60+dateVector(4)*60*60; % calculate time at this instant
        deltaT=timeNow-time0; % elapsed is zero at first excution

    % calculate position of servo point
    if deltaT<t0
        sCoordinate=0.5*accel*deltaT*deltaT;
    elseif(deltaT<t1)
        s0=0.5*accel*t0*t0;
        sCoordinate=velocity*(deltaT-t0)+s0;
    elseif(deltaT<t2)
           s0=0.5*accel*t0*t0;
           s1=velocity*(t1-t0)+s0;
           sCoordinate=velocity*(deltaT-t1)-0.5*accel*(deltaT-t1)*(deltaT-t1)+s1;
    end
 %% when time ends break the loop
    if deltaT>t2
        break;
    end
%% Interpolate theta from the sVec,thetaVec
    for counter=interpolationCounter:(sizesVec-1)
        if(sCoordinate>sVec(counter))
            tetaRange=thetaVec(counter+1)-thetaVec(counter);
            sRange=sVec(counter+1)-sVec(counter);
            thetaVar=(sCoordinate-sVec(counter))...
                *tetaRange/sRange+thetaVec(counter);
            interpolationCounter=counter;
        end
    end

    pPrime=R*ellipseParametricFunction( a,b,thetaVar );
    p=c+pPrime;

    % calculate target transform
    Tt(1:3,4)=p;

    [ qs ] = kukaDLSSolver( qs, Tt, TefTool,numberOfIterationForSolver,lambda );

    for i=1:7
        jPos{i}=qs(i);
    end

	%% Send joint positions to robot
	sendJointsPositions( tKuka,jPos);

end

    % Stop the realtime motion
    realTime_stopDirectServoJoints(tKuka);

    % put the pen up the page
    deltaX=0;deltaY=0;deltaZ=+85.;
    Pos{1}=deltaX;
    Pos{2}=deltaY;
    Pos{3}=deltaZ;
    vel=50;
    movePTPLineEefRelBase( tKuka , Pos, vel);

% Turn off the server
    net_turnOffServer( tKuka );
    fclose(tKuka);

repo: KST-Kuka-Sunrise-Toolbox


text:
function [val, pos] = min3(f)
%MIN3   Global minimum of a CHEBFUN3.
%   [VAL, POS] = MIN3(F) returns the global minimum VAL and the position
%   POS of the global minimum of a CHEBFUN3 object F.
%
% See also CHEBFUN3/MINANDMAX3.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Call MINANDMAX3:
[val, pos] = minandmax3(f);

% Extract minimum:
val = val(1);
pos = pos(1, :);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function airwayAxisMenu(command, varargin)
%"airwayAxisMenu"
%   Handles callbacks from the right click menus for all CERR axes.  Also
%   creates new right click menus in passed axes, and interfaces with the
%   ud.axisInfo field that all CERR viewer axes have.
%
%APA 5/11/2021
%
%Usage:
%   function airwayAxisMenu(command, varargin)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global airwayStateS stateS

if ishandle(command) & strcmpi(get(command, 'type'), 'axes')
    varargin{1} = command;
    command     = 'init';
elseif ~ischar(command)
    error('Invalid call to airwayAxisMenu.');
end

switch upper(command)

    case 'UPDATE_MENU'
        hMenu = gcbo;
        ud = get(hMenu, 'userdata');
        hpV = ud{1};
        elemDoseV = ud{2};
        elemRadV = ud{3};
        elemDiffDistV = ud{4};

        %Wipe out old submenus.
        kids = get(hMenu, 'children');
        delete(kids);

        %Create top level menus.
        hDose        = uimenu(hMenu, 'Label', 'dose - Gy','Callback',...
            'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hpV, elemDoseV});
        hSize        = uimenu(hMenu, 'Label', 'min surface dist - cm', 'Callback',...
            'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hpV, elemRadV});
        %hArea        = uimenu(hMenu, 'Label', 'Area', 'Callback',...
        %    'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hpV, elemAreaV});
        hSizeChange  = uimenu(hMenu, 'Label', 'change in min surface distance - %','Callback',...
            'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hpV, elemDiffDistV});
        nodeStart  = uimenu(hMenu, 'Label', 'Select starting node','Callback',...
            'airwayAxisMenu(''SELECT_START_NODE'')', 'userdata', {hpV});
        nodeStop  = uimenu(hMenu, 'Label', 'Select ending node','Callback',...
            'airwayAxisMenu(''SELECT_END_NODE'')', 'userdata', {hpV});

    case 'INIT_START_STOP_VIEW_NODES'
        hMenu = gcbo;
        %Wipe out old submenus.
        kids = get(hMenu, 'children');
        delete(kids);
        ud = get(hMenu, 'userdata');
        axisType = ud{1};
        hAxis = ud{2};
        hPlot = ud{3};
        basePt = ud{4};
%         hStartStop  = uimenu(hMenu, 'Label', 'Add nodes','Callback',...
%             'airwayAxisMenu(''ADD_START_STOP_NODES'')', 'userdata', {hAxis,hPlot,axisType});
%         hRemoveNodes  = uimenu(hMenu, 'Label', 'Remove nodes','Callback',...
%             'airwayAxisMenu(''REMOVE_START_STOP_NODES'')', 'userdata', {hAxis,hPlot,axisType});
        if strcmpi(axisType,'base')
            %             vf = ud{5};
            %             xFieldV = ud{6};
            %             yFieldV = ud{7};
            %             zUnifV = ud{8};
            %             followupPt = ud{9};
            %             minDistBaseV = ud{10};
            %             radiusDiffV = ud{11};
            %             nodeOrigDoseV = ud{12};
            nodeXyzInterpM = ud{5};
            followupPt = ud{6};
            minDistBaseV = ud{7};
            radiusDiffV = ud{8};
            nodeOrigDoseV = ud{9};
%             hCerrViewer  = uimenu(hMenu, 'Label', 'View location in viewer','Callback',...
%                 'airwayAxisMenu(''SHOW_IN_CERR_VIEWER'')', 'userdata', ...
%                 {hPlot,axisType,basePt,vf,xFieldV,yFieldV,zUnifV,followupPt});
%             hCerrViewer  = uimenu(hMenu, 'Label', 'View location in viewer','Callback',...
%                 'airwayAxisMenu(''SHOW_IN_CERR_VIEWER'')', 'userdata', ...
%                 {hPlot,axisType,basePt,nodeXyzInterpM,followupPt});
            hDose        = uimenu(hMenu, 'Label', 'dose - Gy','Callback',...
                'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hPlot, nodeOrigDoseV});
            hSize        = uimenu(hMenu, 'Label', 'min surface dist - cm', 'Callback',...
                'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hPlot, minDistBaseV});
            %hArea        = uimenu(hMenu, 'Label', 'Area', 'Callback',...
            %    'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hpV, elemAreaV});
            hSizeChange  = uimenu(hMenu, 'Label', 'change in min surface distance - %','Callback',...
                'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hPlot, radiusDiffV});
            hNoColor  = uimenu(hMenu, 'Label', 'No color','Callback',...
                'airwayAxisMenu(''CHANGE_MAP'')', 'userdata', {hPlot, [0.5,0.5,0.5]});

%             if airwayStateS.baseAddNodes
%                 set(hStartStop,'checked', 'on')
%             end
%             if airwayStateS.baseShowInCerrViewer
%                 set(hCerrViewer,'checked', 'on')
%             end
%             if airwayStateS.baseRemoveNodes
%                 set(hRemoveNodes,'checked', 'on')
%             end
        else
%             hCerrViewer  = uimenu(hMenu, 'Label', 'View location in viewer','Callback',...
%                 'airwayAxisMenu(''SHOW_IN_CERR_VIEWER'')', 'userdata', {hPlot,axisType,basePt});
%             if airwayStateS.followupAddNodes
%                 set(hStartStop,'checked', 'on')
%             end
%             if airwayStateS.followupShowInCerrViewer
%                 set(hCerrViewer,'checked', 'on')
%             end
%             if airwayStateS.followupRemoveNodes
%                 set(hRemoveNodes,'checked', 'on')
%             end
        end

    case 'ADD_START_STOP_NODES'
        disp('starting node')
        menuLabel   = get(gcbo, 'Label');
        ud       = get(gcbo, 'userdata');
        hAxis = ud{1};
        hPlot = ud{2};
        axisType = ud{3};
        if strcmpi(axisType,'base')
            airwayStateS.baseAddNodes = 1;
            airwayStateS.baseShowInCerrViewer = 0;
            airwayStateS.baseRemoveNodes = 0;
        else
            airwayStateS.followupAddNodes = 1;
            airwayStateS.followupShowInCerrViewer = 0;
            airwayStateS.followupRemoveNodes = 0;
        end
        % Callback to add nodes
        buttonDownFcn = @pickStatrStopNodes;
        set(hPlot,'ButtonDownFcn',{buttonDownFcn,hAxis,axisType});
        %hp = ud{1};
        %set(hp,'ButtonDownFcn',{@pickStatrStopNodes});

    case 'REMOVE_START_STOP_NODES'
        disp('ending node')
        ud       = get(gcbo, 'userdata');
        hAxis = ud{1};
        hPlot = ud{2};
        axisType = ud{3};
        if strcmpi(axisType,'base')
            airwayStateS.baseAddNodes = 0;
            airwayStateS.baseShowInCerrViewer = 0;
            airwayStateS.baseRemoveNodes = 1;
        else
            airwayStateS.followupAddNodes = 0;
            airwayStateS.followupShowInCerrViewer = 0;
            airwayStateS.followupRemoveNodes = 1;
        end
        % Callback to add nodes
        buttonDownFcn = @pickStatrStopNodes;
        set(hPlot,'ButtonDownFcn',{buttonDownFcn,hAxis,axisType});
        %hp = ud{1};
        %set(hp,'ButtonDownFcn',{@pickStatrStopNodes});

    case 'SHOW_IN_CERR_VIEWER'
        airwayStateS.addNodes = 0;
        airwayStateS.showInCerrViewer = 1;
        airwayStateS.removeNodes = 0;
        disp('location in viewer')
        ud       = get(gcbo, 'userdata');
        %axisType = ud{1};
        hPlot = ud{1};
        % Callback to show nodes in Viewer
        buttonDownFcn = @showNodeInCerrViewer;
        %set(hPlotBase,'ButtonDownFcn',{buttonDownFcn,axisType,vf,xFieldV,yFieldV,zUnifV,...
        %    basePt,followupPt})
        set(hPlot,'ButtonDownFcn',{buttonDownFcn,ud{2:end}})

    case 'CHANGE_MAP'
        %         stateS.viewChanged = 1;
        if isfield(airwayStateS,'hCbar') && ishandle(airwayStateS.hCbar)
            delete(airwayStateS.hCbar)
        end
        menuLabel   = get(gcbo, 'Label');
        ud       = get(gcbo, 'userdata');
        hp = ud{1};
        elemValV = ud{2};
        hAxis = get(hp,'parent');
        set(hp,'cData',elemValV)
        if length(elemValV)==3 %assume single rgb color
            title(hAxis,'Baseline')
            return;
        end
        minDist = min(elemValV);
        maxDist = max(elemValV);
        %         maxMinusMinDist = maxDist - minDist;
        cmapM = CERRColorMap('jetmod');
        %         cmapSiz = size(cmapM,1)-1;
        %         cmapIndV = round((elemValV - minDist) / maxMinusMinDist * cmapSiz) + 1;
        %         for i = 1:length(hp)
        %             if ~isnan(cmapIndV(i))
        %             set(hp,'MarkerFaceColor',[cmapM(cmapIndV(i),:)])
        %             end
        %         end
        set(hAxis,'cLim',[minDist,maxDist])
        colormap(cmapM)
        ticksV = linspace(minDist,maxDist,5);
        tickC = {};
        for i = 1:length(ticksV)
            tickC{i} = num2str(ticksV(i));
        end
        airwayStateS.hCbar = colorbar(hAxis,'Ticks',ticksV,'TickLabels',tickC,...
            'units','normalized','position',[0.45,0.7,0.02,0.2]);
        titleStr = ['Baseline (',menuLabel,')'];
        title(hAxis,titleStr)
        %hCbarLine = line('parent',hCbar,'xdata',[-5 5],...
        %    'ydata',[level level],'color','k','LineWidth',3);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%%%%%%%VERSION 1

%The Gabor filter is basically a Gaussian (with variances sx and sy along x and y-axes respectively)
%modulated by a complex sinusoid (with centre frequencies U and V along x and y-axes respectively)
%described by the following equation
%%
%               1                -1     x  ^    y  ^
%%% G(x,y) = ---------- * exp ([----{(----) 2+(----) 2}+2*pi*i*(Ux+Vy)])
%            2*pi*sx*sy           2    sx       sy

%% Describtion :

%% I : Input image
%% Sx & Sy : Variances along x and y-axes respectively
%% U & V : Centre frequencies  along x and y-axes respectively

%% G : The output filter as described above
%% gabout : The output filtered image

%%  Author : Ahmad poursaberi  e-mail : a.poursaberi@ece.ut.ac.ir
%%          Faulty of Engineering, Electrical&Computer Department,Tehran
%%          University,Iran,June 2004

function [G,gabout] = gaborfilter(I,Sx,Sy,U,V);
I=imread('G:\eye.jpg');

if isa(I,'double')~=1
    I = double(I);
end

for x = -fix(Sx):fix(Sx)
    for y = -fix(Sy):fix(Sy)
        G(fix(Sx)+x+1,fix(Sy)+y+1) = (1/(2*pi*Sx*Sy))*exp(-.5*((x/Sx)^2+(y/Sy)^2)+2*pi*i*(U*x+V*y));
    end
end

Imgabout = conv2(I,double(imag(G)),'same');
Regabout = conv2(I,double(real(G)),'same');

gabout = uint8(sqrt(Imgabout.*Imgabout + Regabout.*Regabout));

author: Sable

repo: mcbench-benchmarks


text:
host = 'https://hedtools.ucsd.edu/hed';
%host = 'http://127.0.0.1:5000/';
errorMap = containers.Map('KeyType', 'char', 'ValueType', 'any');
errorMap('testGetServices') = testGetServices(host);
errorMap('testEventServices') = testEventServices(host);
errorMap('testEventSearchServices') = testEventSearchServices(host);
errorMap('testSidecarServices') = testSidecarServices(host);
errorMap('testSpreadsheetServices') = testSpreadsheetServices(host);
errorMap('testStringServices') = testStringServices(host);

%% Output the errors
fprintf('\n\nOverall error report:\n');
keys = errorMap.keys();
for k = 1:length(keys)
    errors = errorMap(keys{k});
    if isempty(errors)
        fprintf('\t%s: no errors\n', keys{k});
    else
        fprintf('\t%s:\n', keys{k});
        for n = 1:length(errors)
           fprintf('\t\t%s\n', keys{k}, errors{n});
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% len is frame size, len1 is overlap
function wav = complexSpec2wav(complex_x, framelen, overlap)
nFFT = size(complex_x,2)*2-2;
img = sqrt(-1);

complex_x(:,nFFT/2+2:nFFT) = conj(complex_x(:,nFFT/2:-1:2));
xi = ifft(complex_x');
xi = real(xi);
wav = my_ola(xi, framelen, overlap);

%A = [1 -0.97];
%wav = filter(1, A, wav);
end

author: singaxiong

repo: SignalGraph


text:
function s = ymdf_to_s_common ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_S_COMMON writes a Common YMDF date into a string.
%
%  Format:
%
%    CE YYYY/MM/DD.FF
%    BCE YYYY/MM/DD.FF
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, real F, the YMDF date.
%
%    Output, string S, a representation of the date.
%

%
%  Copy the input.
%
  y2 = y;
  m2 = m;
  d2 = d;
  f2 = f;
%
%  Check the input.
%
  [ y2, m2, d2, f2, ierror ] = ymdf_check_common ( y2, m2, d2, f2 );

  if ( ierror ~= 0 )
    s = sprintf ( '?' );
    return
  end

  if ( 0 <= y2 )
    s = sprintf ( 'CE %d/%d/%d/%f', y2, m2, d2, f2 );
  else
    s = sprintf ( 'BCE %d/%d/%d/%f', -y2, m2, d2, f2 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rts_clearInputBuffer(obj, event)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

    flushinput(obj);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% RangeFrame Demo
%
%  James Houghton
%  James.P.Houghton@gmail.com
%  December 4 2009
%
% Requires rangeframe.m
%
% In this demo we use the range frame to display the range of the first and last quartile,
% and the mean of a bivariate scatter plot. The same range frame could display other information,
% such as median, standard deviation, etc. with appropriate explanation.
%

%% Generate data to plot
x = exp(0 + .5*randn(1,100));                 % the log-normal is an interesting distribution
y = x + .15*randn(1,100);

%% determine relevant range frame points
x_sort = sort(x);
x_Q1 = x_sort(round(length(x_sort)/4));       %approximate 1st quartile
x_Q3 = x_sort(round(length(x_sort)*3/4));     %approximate 3rd quartile

y_sort = sort(y);
y_Q1 = y_sort(round(length(y_sort)/4));       %approximate 1st quartile
y_Q3 = y_sort(round(length(y_sort)*3/4));     %approximate 3rd quartile

x_frame = [min(x), x_Q1, mean(x), x_Q3, max(x)];
y_frame = [min(y), y_Q1, mean(y), y_Q3, max(y)];

%% plot the data
% we can use a simple plot command, as the changes are made at the axis level
% instead of at the figure level.
plot(x, y, '.');

% apply a range frame to the data we are interested in
ax1 = rangeframe(gca, x_frame, y_frame);

% we can make changes to the axis after it is returned from the rangeframe function
set(ax1, 'DefaultTextInterpreter', 'Latex');    %use the latex font interpreter
set(ax1, 'XTick', [0, 1, 2, 3, 4], ...          %set the axis values to display
         'Ytick', [0, 1, 2, 3, 4]);

xlabel(ax1, 'X data');                          %label everything
ylabel(ax1, 'Y data');
title(ax1, 'Statistical data with quartiles and mean');


text:
  function ss = outer_sum(xx,yy)
%|function ss = outer_sum(xx,yy)
%|
%| compute an "outer sum" x + y'
%| that is analogous to the "outer product" x * y'
%|
%| in
%|	xx	[nx 1]
%|	yy	[1 ny]
%|		more generally: xx [(dim)] + yy [L,1] -> xx [(dim) LL]
%| out
%|	ss [nx ny]	ss(i,j) = xx(i) + yy(j)
%|
%| Copyright 2001, Jeff Fessler, University of Michigan

if ~nargin, ir_usage, end
if streq(xx, 'test'), outer_sum_test, return, end

% for 1D vectors, allow rows or cols for backward compatibility
if ndims(xx) == 2 && min(size(xx)) == 1 && ndims(yy) == 2 && min(size(yy)) == 1
	nx = length(xx);
	ny = length(yy);
	xx = repmat(xx(:), [1 ny]);
	yy = repmat(yy(:)', [nx 1]);
	ss = xx + yy;
return
end

%if size(xx,1) == 1
%	warn 'xx is a row vector? are you sure?'
%end

% otherwise, xx is not a vector, but yy must be

xdim = size(xx);
ydim = size(yy);
if ndims(yy) ~= 2 || min(size(yy)) ~= 1
	error 'yy must be a vector'
end
sdim = [xdim length(yy)];
xo = repmat(xx, [ones(1, length(xdim)) length(yy)]); % [xdim] -> [xdim ny]
yo = repmat(yy(:), [1 xdim]); yo = permute(yo, [2 3 1]);
% yo = repmat(yy(:)', [xdim 1 1]); % not sure how to make this work.
ss = xo + yo;

function outer_sum_test
%xx = [1:4];
%yy = [0:10:50];
%pr outer_sum(xx,yy)
pr 'outer_sum([1:4], [0:10:50])'
xx = outer_sum([1:4], [0:10:50]);
pr 'outer_sum(xx, [100 200])'

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Visualize validation results of system identification

%% Show training and prediction
clear ph
h = figure;
ccolors = get(gca,'colororder');
subplot(3,1,1), box on, hold on
plot([tA(1),tA(1)],[-15 260],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
text(tA(1)-30,230,'Training', 'FontSize',12)
text(5+tA(1),230,'Validation', 'FontSize',12)
plot([t;tA],[x(:,1);xA(:,1)],'Color',[.4 .4 .4],'LineWidth',1.5);
plot(tA,xA(:,1),'k','LineWidth',1.5);
plot(tB,xB(:,1),'--','Color','r','LineWidth',1.5);
grid on
ylim([-25 25])
ylabel('x_1','FontSize',13)
set(gca,'FontSize',13)
subplot(3,1,2), box on, hold on
plot([tA(1) tA(1)],[0 60],':')
plot([t;tA],[x(:,2);xA(:,2)],'Color',[.4 .4 .4],'LineWidth',1.5);
plot(tA,xA(:,2),'k','LineWidth',1.5);
plot(tB,xB(:,2),'--','Color','r','LineWidth',1.5);
grid on
ylim([-25 25])
ylabel('x_3','FontSize',13)
set(gca,'FontSize',13)

subplot(3,1,3), box on, hold on
plot([tA(1) tA(1)],[0 60],':')
ph(1) = plot([t;tA],[x(:,3);xA(:,3)],'Color',[.4 .4 .4],'LineWidth',1.5);
ph(2) = plot(tA,xA(:,3),'k','LineWidth',1.5);
ph(3) = plot(tB,xB(:,3),'--','Color','r','LineWidth',1.5);
l1=legend(ph,'Training','Validation',ModelName);
set(l1,'Location','NorthWest')
grid on
ylim([0 50])
ylabel('x_3','FontSize',13)
set(gca,'FontSize',13)
xlabel('Time','FontSize',13)
set(gca,'FontSize',13)

set(h,'Units','Inches');
set(gcf,'Position',[1 1 6. 5.5])
pos = get(h,'Position');
set(h,'PaperPositionMode','Auto','PaperSize',[pos(3), pos(4)])
print(h,'-painters','-depsc2', '-loose','-cmyk', [figpath,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation.eps']],'-r0');

%%

clear ph
figure,box on, hold on,
ccolors = get(gca,'colororder');
plot([tA(1),tA(1)],[min(xB(:)) max([xA(:);xB(:)])],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
ylim([min([xA(:);xB(:)]) max([xA(:);xB(:)])])
plot(t,x(:,1),'Color',ccolors(1,:),'LineWidth',1);
plot(t,x(:,2),'Color',ccolors(2,:),'LineWidth',1);
ph(1) = plot([t;tA],[x(:,1);xA(:,1)],'Color',ccolors(1,:),'LineWidth',1);
ph(2) = plot([t;tA],[x(:,2);xA(:,2)],'Color',ccolors(2,:),'LineWidth',1);
ph(3) = plot([t;tA],[x(:,3);xA(:,3)],'Color',ccolors(3,:),'LineWidth',1);
ph(4) = plot(tB,xB(:,1),'-.','Color',ccolors(1,:)-[0 0.2 0.2],'LineWidth',2);
ph(5) = plot(tB,xB(:,2),'-.','Color',ccolors(2,:)-[0.1 0.2 0.09],'LineWidth',2);
ph(6) = plot(tB,xB(:,3),'-.','Color',ccolors(3,:)-[0.1 0.2 0.09],'LineWidth',2);
grid off
xlim([0 tB(end)])
xlabel('Time')
ylabel('xi')
t1 = text(tA(1)-5,40,'Training', 'FontSize',12);
t2 = text(1+tA(1),40,'Validation', 'FontSize',12);
lh = legend(ph([1,4]),'True',ModelName);
lh.Position = [lh.Position(1)-0.5,lh.Position(2)-0.15,lh.Position(3:4)];
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_OneFig.eps']);

delete(lh); delete(t1), delete(t2);
print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Validation_OneFig_noleg.eps']);
%% Actuation signal
clear ph
figure,box on, hold on,
ccolors = get(gca,'colororder');
plot([tA(1),tA(1)],[min([u;uv]) max([u;uv])],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
plot(t,u,'-k','LineWidth',1);
plot(tv,uv,'-k','LineWidth',1);
grid off
ylim([min([u;uv])+0.05*min([u;uv]) max([u;uv])+0.05*max([u;uv])])
xlim([0 tB(end)])
xlabel('Time')
ylabel('Input')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_Actuation_OneFig.eps']);

author: eurika-kaiser

repo: SINDY-MPC


text:
function [cl,cl2] = hewma_plot_bivariate(varargin)
% Visualize a change_point map stored in hewma_cp.img
% and a run-length map stored in hewma_runlen.img
% (output from hewma2)
%
% :Usage:
% ::
%
%     [cl,cl2] = hewma_plot_bivariate([Method],[optional overlay image])
%
% and classify voxels into groupings based on locations in the bivariate
% space

Method = 'group';
if length(varargin), Method = varargin{1};, end
if length(varargin) > 1, ovl = varargin{2}; else, ovl = which('scalped_single_subj_T1.img');, end

doplot = 1;

% ----------------------------------------------------
% *
% *
% * get data, x, and list of all voxels to go with each data point, CLU
% *
% *
% ----------------------------------------------------
switch Method

    case 'group'

        mask = 'hewma_sig.img';
        t = 'hewma_t.img';

        name1 = 'hewma_cp.img';
        name2 = 'hewma_runlen.img';
        name2 = 'gausslongest.img';

    case 'difference'
        mask = 'hewma_sigdiff.img';
        t = 'hewma_tdiff.img';

        name1 = 'hewma_cpdiff.img';
        name2 = 'hewma_runlendiff.img';

    otherwise
        error('Method (first arg) must be group or difference.');
end

% ----------------------------------------------------
% initial viewing of sig map
% ----------------------------------------------------

cl = mask2clusters(mask,t);

if doplot, cluster_orthviews(cl,'overlay',ovl,'bivalent'); %,{[1 0 0]});,
    %cluster_orthviews(cl,{[1 0 0]},'add');,
    %cluster_orthviews(cl,{[1 0 0]},'add');,
end

CLU = clusters2CLU(cl);

% ----------------------------------------------------
% get mask for pos and neg activation
% ----------------------------------------------------
v = spm_read_vols(spm_vol(t));

% ----------------------------------------------------
% load data, get data for each voxel in CLU
% ----------------------------------------------------

v1 = spm_read_vols(spm_vol(name1));
v2 = spm_read_vols(spm_vol(name2));

for i = 1:size(CLU.XYZ,2)
    % data in voxel list
    x(i,1) = v1(CLU.XYZ(1,i),CLU.XYZ(2,i),CLU.XYZ(3,i));
    x(i,2) = v2(CLU.XYZ(1,i),CLU.XYZ(2,i),CLU.XYZ(3,i));
end

x = round(x);

go = input('Make tables?');
if go, mincs = input('Min cluster size?  (in voxels): ');, end

% ----------------------------------------------------
% POSITIVE
% ----------------------------------------------------

% select positive points and CLU
wh = find(CLU.Z > 0);
dat = x; dat = dat(wh,:);
CLU2 = CLU; CLU2.XYZ = CLU2.XYZ(:,wh); CLU2.XYZmm = CLU2.XYZmm(:,wh); CLU2.Z = CLU2.Z(:,wh);

[cl2,nclasses,colors] = cluster_kmeans_parcel(dat,CLU2,1);

% -------------------------------------------------------------------
% * tables
% -------------------------------------------------------------------

if go

    fprintf(1,'\nPositive deviations\n-----------------------------------\n');

    for i = 1:length(cl2);

        whomit = cat(1,cl2{i}.numVox); whomit(whomit >= mincs) = 0; cl2{i}(find(whomit)) = [];

        fprintf(1,'Class %3.0f, color = %s\n-----------------------------------\n',i,num2str(colors{i}));
        if isempty(cl2{i})
            fprintf(1,'No clusters large enough.\n');
        else
            cluster_table(cl2{i});
        end
        fprintf(1,'\n-----------------------------------\n')
    end
end

cl = cl2;

% ----------------------------------------------------
% NEGATIVE
% ----------------------------------------------------

dodecrease = input('Look at deactivations?');

if dodecrease

    % select negative points and CLU
    wh = find(CLU.Z < 0);
    dat = x; dat = dat(wh,:);
    CLU2 = CLU; CLU2.XYZ = CLU2.XYZ(:,wh); CLU2.XYZmm = CLU2.XYZmm(:,wh); CLU2.Z = CLU2.Z(:,wh);

    [cl2,nclasses,colors] = cluster_kmeans_parcel(dat,CLU2,1);

    % -------------------------------------------------------------------
    % * tables
    % -------------------------------------------------------------------

    if go

        fprintf(1,'\nNegative deviations\n-----------------------------------\n');

        for i = 1:length(cl2);

           whomit = cat(1,cl2{i}.numVox); whomit(whomit >= mincs) = 0; cl2{i}(find(whomit)) = [];

            fprintf(1,'Class %3.0f, color = %s\n-----------------------------------\n',i,num2str(colors{i}));
            if isempty(cl2{i})
                fprintf(1,'No clusters large enough.\n');
            else
                cluster_table(cl2{i});
            end
            fprintf(1,'\n-----------------------------------\n')
        end
    end

    cl = [cl cl2];

end

cl2 = cl; % cl2 has output in cell array, sep. by class
cl = cat(2,cl{:}); % cl is single vector

% -------------------------------------------------------------------
% * cluster extraction
% -------------------------------------------------------------------

go = input('Save timeseries in clusters?');

if go
    nm= 'hewma_timeseries.mat'; ctr = 1;

    while exist(nm) == 2
        nm = [nm(1:end-4) num2str(ctr) '.mat'];
        ctr = ctr+1;
    end
    fprintf(1,'Extracting cluster data and saving in %s\n',nm);
    mincs = input('Min cluster size?  (in voxels): ');
    whomit = cat(1,cl.numVox); whomit(whomit >= mincs) = 0; cl(find(whomit>0)) = [];

    cl = hewma_save_timeseries(cl,mincs);
    eval(['save ' nm ' cl']);

end

% -------------------------------------------------------------------
% * timeseries plotting
% -------------------------------------------------------------------

go = input('Plot timeseries interactively?');

if go
    if ~(exist('EXPT') == 1)
        file = spm_get(1,'*mat','Select EXPT.mat',pwd);
        [dd,ff,ee] = fileparts(file);
        cd(dd)
        load(file)
    end

    % used in button-up fcn callback
    E2 = EXPT;
    clear EXPT

    global VOL
    global f
    global f2
    global EXPT
    EXPT = E2;

    set(gcf,'WindowButtonUpFcn','[dat,files,stats,mycov] = hewma_plot_coord_btnupfcn;')

    % get coordinate mapping matrix
    VOL = struct('M',cl(1).M);

    % prepare figure
    f1 = figure('Color','w','Name','Hewma plots');
    f = f1;     % button callback uses figure f

    stop = input('Click on figure to plot.  Press return to exit.');

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore