exp_SonEMS_3.m

Go to the documentation of this file.

% =========================================================================== %
%> @brief final analysis script for the SonEMS data
%>
%> this script filters the data of the feet and compares certain gait
%> parameters with the VICON data. the data is from the SonEMS experiments.
%>
%> @file exp_SonEMS_3.m
%>
%> @note this script is written to use the parallel processing toolbox. it
%> works only with the "Processes" pool. The script is very memory intensive,
%> therefore the number of workers should be set accordingly.
%>
%> @copyright see the file @ref LICENSE in the root directory of the repository
% =========================================================================== %
 
clear;
close all;
clc;
 
pause(2); % give the system some time to close the figures
 
%%> TASKS TO DO
%> here the tasks that need to be done again can be marked. if a task is not done
%> yet OR is marked as true, it will be executed. if a task is marked as false,
%> it will be skipped (if it has be done before).
 
%% base folder, participant list and ekf tuning parameters
basepath = "D:/ETH_Data/";
info_sheet = readtable(fullfile(basepath, "participant_info.xlsx"), ...
    FileType = "spreadsheet", ...
    TextType = "string", ...
    DatetimeType = "datetime", ...
    VariableNamingRule = "modify", ...
    ReadVariableNames = true, ...
    Range = "A1:I38");
 
% get the participants
current_run = "Test"; % "Tuning" or "Test"
ptps = unique(info_sheet.ID(info_sheet.(current_run) == 1));
 
% ekf tuning parameters file
if strcmp(current_run, "Test")
    tuning_files = "C:\Users\moss\Code\MT_moss_IMU_sensor_fusion\experiments\ekf_tuning.json"; % Test
elseif strcmp(current_run, "Tuning")
    tuning_files = dir(fullfile(basepath, "Tuning", "Trial_14", "*.json"));
    tuning_files = string(fullfile({tuning_files.folder}, {tuning_files.name}))';
    tuning_files = "C:\Users\moss\Code\MT_moss_IMU_sensor_fusion\experiments\ekf_tuning.json"; % Test
else
    error("Unknown Dataset.")
end
 
% log file
diary(fullfile(basepath, "exp_SonEMS_3.log"))
fprintf("\n\nThis is the beginning of a new log file.\n\nThe column headers are:\n")
fprintf("Time; Log Level; Participant; Message\n\n")
 
log_entry("Start processing SonEMS data")
log_entry("Running " + current_run + " set", 0, 4)
 
% Select the processes pool to use (not the thread pool)
num_workers = 16; % change this value, if you want to use more or less workers
num_workers = min(num_workers, numel(ptps)); % limit the number of workers to the number of participants
currentPool = gcp("nocreate");
if ~isempty(currentPool)
    if ~strcmp(class(currentPool), "parallel.ProcessPool")
        delete(currentPool);
        parpool("Processes", num_workers);
        log_entry("Changed the pool to Processes", 0, 5)
    elseif currentPool.NumWorkers ~= num_workers
        delete(currentPool);
        parpool("Processes", num_workers);
        log_entry("Changed the pool to Processes with " + num_workers + " workers", 0, 5)
    else
        log_entry("Using the existing pool of Processes", 0, 5)
    end
else
    parpool("Processes", num_workers);
    log_entry("Created a new pool of Processes", 0, 5)
end
 
% define a vartiabel to store the time of the last log entry of the RAM usage
ram_last_logged = datetime(1970, 1, 1); % preinitialize to a long time ago
 
% loop over the tuning files
wb = waitbar(0, "Processing participants", "Name", "Processing participants");
 
% allocate the grid search results table
grid_search_results = table('Size', [numel(tuning_files), 20], ...
    'VariableTypes', ["string", "datetime", repmat("double", 1, 18)], ...
    'VariableNames', [...
    "TuningFile", ...
    "time", ...
    "StrideLengthMu", ...
    "StrideLengthDiffMu", ...
    "StepLengthMu", ...
    "StepLengthDiffMu", ...
    "StepWidthMu", ...
    "StepWidthDiffMu", ...
    "StrideLengthStd", ...
    "StrideLengthDiffStd", ...
    "StepLengthStd", ...
    "StepLengthDiffStd", ...
    "StepWidthStd", ...
    "StepWidthDiffStd", ...
    "StrideLengthMSE", ...
    "StrideLengthDiffMSE", ...
    "StepLengthMSE", ...
    "StepLengthDiffMSE", ...
    "StepWidthMSE", ...
    "StepWidthDiffMSE"
    ]);
 
for i_tf = 1:numel(tuning_files)
    
    % get the tuning file
    tuning_file = tuning_files(i_tf);
    
    % save the tuning file
    grid_search_results.TuningFile(i_tf) = tuning_file;
    
    if numel(tuning_files) > 1
        log_entry("Start processing with tuning file " + tuning_file, 0, 4)
    end
    
    % preallocate the futures
    f_filter(1:numel(ptps)) = parallel.FevalFuture;
    f_gait_parameters(1:numel(ptps)) = parallel.FevalFuture;
    
    % queue all participants for filter processing
    for ptp = 1:numel(ptps)
        
        f_filter(ptp) = parfeval(...
            @process_participant, ... % function to execute
            0, ... % number of output arguments
            basepath, ...
            ptps(ptp), ...
            ekf_tuning_file    = tuning_file, ...
            GT_POS             = false, ...
            FIND_WALKING_AXES  = false, ...
            FIND_WALKING_SPEED = false, ...
            ACCUMULATE_DIST    = false, ...
            ALIGN_IMU_VICON    = false, ...
            FIND_YAW           = true, ...
            FILTER             = true);
        
        log_entry("Participant " + ptps(ptp) + " queued")
        
    end
    
    % wait for all futures to finish, show intermediate results in waitbar
    filter_finished = false;
    error_occured = 0;
    while ~filter_finished
        finished = sum(string({f_filter.State}) == "finished");
        queued = sum(string({f_filter.State}) == "queued");
        running = sum(string({f_filter.State}) == "running");
        failed = sum(string({f_filter.State}) == "failed") + error_occured;
        
        waitbar(finished / numel(f_filter), wb, sprintf("Filtering (%i/%i) ...\nFinished: %d, Queued: %d, Running: %d, Failed: %d", i_tf, numel(tuning_files), finished, queued, running, failed))
        
        if finished + failed >= numel(f_filter)
            filter_finished = true;
        end
        
        % log RAM usage every 1 minute
        if seconds(datetime("now") - ram_last_logged) > 60
            [~, curr_ram] = memory;
            
            log_entry(sprintf("RAM: %.1f GB of %.1f GB remaining.", curr_ram.PhysicalMemory.Available / 1e9, curr_ram.PhysicalMemory.Total / 1e9), 0, 5)
            
            ram_last_logged = datetime("now");
        end
        
        pause(1)
        
        % check the state of the futures
        for ptp = 1:numel(f_filter)
            
            % check if future is already read
            if f_filter(ptp).Read == 1
                continue
            end
            
            % check for errors
            if ~isempty(f_filter(ptp).Error)
                
                % get the SonEMS String
                err_ptp = f_filter(ptp).InputArguments{2};
                
                % stop the future
                f_filter(ptp).cancel
                
                % log the error
                log_entry("Unhandled error in participant " + err_ptp, 0, 1)
                log_entry(f_filter(ptp).Error.message, 0, 2)
                log_entry("Stack trace:", 0, 2)
                for i = 1:numel(f_filter(ptp).Error.stack)
                    log_entry(f_filter(ptp).Error.stack(i).name + " (Line: " + string(f_filter(ptp).Error.stack(i).line) + ")", 0, 2)
                end
                
                % read the future (this will throw an exeption but will mark it as read)
                try
                    fetchOutputs(f_filter(ptp));
                catch
                    % do nothing
                end
                
                error_occured = error_occured + 1;
                
                continue
            end
            
            % check if the future is finished
            if strcmp(f_filter(ptp).State, "finished")
                
                % read the future (this can throw an exeption but will mark it as read)
                try
                    fetchOutputs(f_filter(ptp));
                catch
                    % do nothing
                end
                
                % log the finished participant
                log_entry("Participant " + ptps(ptp) + " finished with filter processing", 0, 4)
                
                % start the future for the gait parameters
                f_gait_parameters(ptp) = parfeval(...
                    @calculate_gait_parameters, ... % function to execute
                    0, ... % number of output arguments
                    basepath, ...
                    ptps(ptp), ...
                    RECALCULATE = true, ...
                    create_plots = true);
                
                % log the queued participant
                log_entry("Participant " + ptps(ptp) + " queued for gait parameter calculation", 0, 4)
                
            end
        end
    end
    
    log_entry("All participants processed (filtered)")
    
    % wait for all gait parameters to finish, show intermediate results in waitbar
    gait_parameters_finished = false;
    while ~gait_parameters_finished
        finished = sum(string({f_gait_parameters.State}) == "finished");
        queued = sum(string({f_gait_parameters.State}) == "queued");
        running = sum(string({f_gait_parameters.State}) == "running");
        failed = sum(string({f_gait_parameters.State}) == "failed");
        unavailable = sum(string({f_gait_parameters.State}) == "unavailable");
        
        waitbar(finished / numel(f_filter), wb, sprintf("Gait Parameters (%i/%i) ...\nFinished: %d, Queued: %d, Running: %d, Failed: %d", i_tf, numel(tuning_files), finished, queued, running, failed))
        
        if finished + failed + unavailable >= numel(f_gait_parameters)
            gait_parameters_finished = true;
        end
        
        pause(1)
        
        % check the state of the futures
        for ptp = 1:numel(f_gait_parameters)
            
            % check if future is already read
            if f_gait_parameters(ptp).Read == 1
                continue
            end
            
            % check if the future is finished
            if strcmp(f_gait_parameters(ptp).State, "finished")
                
                % read the future (this can throw an exeption but will mark it as read)
                try
                    fetchOutputs(f_gait_parameters(ptp));
                catch
                    % do nothing
                end
                
                % log the finished participant
                log_entry("Participant " + ptps(ptp) + " finished with gait parameter calculation")
                
            end
            
            % check for errors
            if ~isempty(f_gait_parameters(ptp).Error)
                
                % get the SonEMS String
                err_ptp = f_gait_parameters(ptp).InputArguments{2};
                
                % stop the future
                f_gait_parameters(ptp).cancel
                
                % log the error
                log_entry("Unhandled eror in participant " + err_ptp, 0, 1)
                log_entry(f_gait_parameters(ptp).Error.message, 0, 2)
                log_entry("Stack trace:", 0, 2)
                for i = 1:numel(f_gait_parameters(ptp).Error.stack)
                    log_entry(f_gait_parameters(ptp).Error.stack(i).name + " (Line: " + string(f_gait_parameters(ptp).Error.stack(i).line) + ")", 0, 2)
                end
                
                % read the future (this will throw an exeption but will mark it as read)
                try
                    fetchOutputs(f_gait_parameters(ptp));
                catch
                    % do nothing
                end
            end
            
        end
        
    end
    
    %% create a summary of the gait parameters
    log_entry("Start creating gait parameter summary")
    results = create_gait_parameter_summary(basepath, ptps, tuning_file, current_run);
    
    % save the results
    grid_search_results.time(i_tf) = datetime("now");
    grid_search_results.StrideLengthMu(i_tf) = results.mu_stride;
    grid_search_results.StepLengthMu(i_tf) = results.mu_step;
    grid_search_results.StepWidthMu(i_tf) = results.mu_width;
    grid_search_results.StrideLengthDiffMu(i_tf) = results.mu_stride_d;
    grid_search_results.StepLengthDiffMu(i_tf) = results.mu_step_d;
    grid_search_results.StepWidthDiffMu(i_tf) = results.mu_width_d;
    grid_search_results.StrideLengthStd(i_tf) = results.std_stride;
    grid_search_results.StepLengthStd(i_tf) = results.std_step;
    grid_search_results.StepWidthStd(i_tf) = results.std_width;
    grid_search_results.StrideLengthDiffStd(i_tf) = results.std_stride_d;
    grid_search_results.StepLengthDiffStd(i_tf) = results.std_step_d;
    grid_search_results.StepWidthDiffStd(i_tf) = results.std_width_d;
    grid_search_results.StrideLengthMSE(i_tf) = results.mse_stride;
    grid_search_results.StepLengthMSE(i_tf) = results.mse_step;
    grid_search_results.StepWidthMSE(i_tf) = results.mse_width;
    grid_search_results.StrideLengthDiffMSE(i_tf) = results.mse_stride_d;
    grid_search_results.StepLengthDiffMSE(i_tf) = results.mse_step_d;
    grid_search_results.StepWidthDiffMSE(i_tf) = results.mse_width_d;
    
    if current_run == "Tuning"
        save(fullfile(basepath, "Tuning", "grid_search_results.mat"), "grid_search_results");
        copyfile(fullfile(basepath, "Tuning", "grid_search_results.mat"), fullfile(basepath, "Tuning", string(datetime("now", "Format", "yy-MM-dd_HH-mm-ss")) + "_grid_search_results.mat"));
    end
    
    % log the end of the processing
    log_entry("All participants processed (gait parameters calculated)")
    
end
 
% close the waitbar
close(wb)