% 02/15/2020 Makoto. Modified.
% 02/04/2020 Makoto. Further fix made. Re-calculated.
% 01/27/2020 Makoto. Bug fixed. Re-calculated.
% 01/22/2020 Makoto. Used for test.
% 01/22/2017 Makoto. Bug in choosing common channels fixed.
% 12/27/2017 Makoto. Changed to fit rigid body only.
% 12/15/2017 Makoto. Used.
% 09/06/2017 Makoto. Data point rejection based on distance from Head/Torso centroid implemented.
% 08/25/2017 Makoto. Created.

addpath /data/mobi/Audiomaze/makoto/code/phaseSpaceCorrectionTools

% 152, 153, 154 (Subj829) skipped for no reason.
% 173 suffers from 40.9% NaN--but let's give it a try anyway. 
%allSets = dir('/data/mobi/Audiomaze/makoto/p01_imported/subj829*.set');
allSets = dir('/data/mobi/Audiomaze/makoto/p01_imported/*.set');

fullChannelList = zeros(length(allSets), 5);
for setIdx = 1:length(allSets)
    
    % Load .set file.
    EEG = pop_loadset('filename', allSets(setIdx).name, 'filepath', '/data/mobi/Audiomaze/makoto/p01_imported', 'loadmode', 'info');
    
    % Obtain phaseSpace data.
    phaseSpaceData = EEG.etc.audiomaze.phaseSpaceData;
    
    % Skip if phaseSpaceData is all zeros.
    if ~any(phaseSpaceData(:))
        continue
    end
    
    % Remove the last channel. Then there are 4 x 32 = 128 time series; x, y, z, and signal quality.
    phaseSpaceData = phaseSpaceData(1:end-1,:);
    phaseSpaceData = reshape(phaseSpaceData, [4 32 size(phaseSpaceData,2)]);
    
    phaseSpaceData(phaseSpaceData==0) = NaN;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Fit rigid body for the five body parts. %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    lessThanThreeMarkerIdx       = logical(zeros(5, size(phaseSpaceData,3)));
    lessThanThreeLedMarkerLabels = {'Head' 'Torso' 'Hand' 'footR' 'footL'};
    
            %%%%%%%%%%%%%%%%%%%%%
            %%% Process head. %%%
            %%%%%%%%%%%%%%%%%%%%%
            headIdx = [1:7];
            phaseSpaceDataSelected  = squeeze(phaseSpaceData(1, headIdx,:));
            goodIdx   = find(sum(isnan(phaseSpaceDataSelected))==0);
            goodHeadRatio = length(goodIdx)/length(phaseSpaceDataSelected);
            
            %%%%%%%%%%%%%%%%%%%%%%
            %%% Process torso. %%%
            %%%%%%%%%%%%%%%%%%%%%%
            torsoIdx = [8:11 21 27];
            phaseSpaceDataSelected  = squeeze(phaseSpaceData(1, torsoIdx,:));
            goodIdx   = find(sum(isnan(phaseSpaceDataSelected))==0);
            goodTorsoRatio = length(goodIdx)/length(phaseSpaceDataSelected);

            %%%%%%%%%%%%%%%%%%%%%
            %%% Process hand. %%%
            %%%%%%%%%%%%%%%%%%%%%
            handIdx = [16:20];
            phaseSpaceDataSelected  = squeeze(phaseSpaceData(1, handIdx,:));
            goodIdx   = find(sum(isnan(phaseSpaceDataSelected))==0);
            goodHandRatio = length(goodIdx)/length(phaseSpaceDataSelected);
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%% Process right foot. %%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%
            footR_Idx = [23:26];
            phaseSpaceDataSelected  = squeeze(phaseSpaceData(1, footR_Idx,:));
            goodIdx   = find(sum(isnan(phaseSpaceDataSelected))==0);
            goodFootR_Ratio = length(goodIdx)/length(phaseSpaceDataSelected);
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%
            %%% Process left foot. %%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%
            footL_Idx = [29:32];
            phaseSpaceDataSelected  = squeeze(phaseSpaceData(1, footL_Idx,:));
            goodIdx   = find(sum(isnan(phaseSpaceDataSelected))==0);
            goodFootL_Ratio = length(goodIdx)/length(phaseSpaceDataSelected);
    
            fullChannelList(setIdx,:) = [goodHeadRatio goodTorsoRatio goodHandRatio goodFootR_Ratio goodFootL_Ratio];
end

save fullChannelList fullChannelList


%%

load('/data/common/mobi/Experiments/Audiomaze/makoto/code/fullChannelList')
load('/data/common/mobi/Experiments/Audiomaze/makoto/code/rigidBodyApplicableRatio')


figure

subplot(1,2,1)
partsMean = mean(fullChannelList);
partsStd  = std(fullChannelList);
partsMin  = min(fullChannelList);
partsMax  = max(fullChannelList);
barHandle1 = bar(partsMean);
hold on
errorbarHandle1 = errorbar(1:5, partsMean, partsStd/sqrt(17));
set(errorbarHandle1, 'linestyle', 'none', 'color', [0 0 0])
xlim([0.5 5.5])
ylim([0 1])
set(gca, 'xticklabel', {'Head' 'Torso' 'Hand' 'FootR' 'FootL'})
set(barHandle1, 'faceColor', [0.93 0.96 1])
ylabel('Probability')
title('All markers available (+/-SE)')



subplot(1,2,2)
rigidBodyApplicableRatio(10,9,:) = rigidBodyApplicableRatio(10,10,:);
rigidBodyApplicableRatio(:,10,:) = [];
mazeMean = cat(2, nanmean(rigidBodyApplicableRatio(:,1:3,:),2), ...
                  nanmean(rigidBodyApplicableRatio(:,4:6,:),2), ...
                  nanmean(rigidBodyApplicableRatio(:,7:9,:),2), ...
                  nanmean(rigidBodyApplicableRatio(:,10:12,:),2));
barHandle2 = bar(squeeze(nanmean(nanmean(mazeMean,1),2)));
hold on
permutedData = permute(mazeMean, [3 1 2]);
errorbarHandle2 = errorbar(1:5, squeeze(nanmean(nanmean(mazeMean,1),2)), nanstd(permutedData(:,:), 0, 2)/sqrt(size(permutedData,2)*size(permutedData,3)));
set(errorbarHandle2, 'linestyle', 'none', 'color', [0 0 0])
xlim([0.5 5.5])
ylim([0 1])
set(gca, 'xticklabel', {'Head' 'Torso' 'Hand' 'FootR' 'FootL'})
set(barHandle2, 'faceColor', [0.93 0.96 1])
ylabel('Probability')
title('>= 3 markers available (+/-SE)')



set(findall(gcf, '-property', 'fontsize'), 'fontsize', 12)
set(gcf, 'position', [5         884        1018         430]);
print -dsvg beforeAndAfterRigidBodyCorrection



