diff --git a/sami_examples/DemoOne_BehavEmotionRating.txt b/sami_examples/DemoOne_BehavEmotionRating.txt new file mode 100644 index 0000000..10c0fdf --- /dev/null +++ b/sami_examples/DemoOne_BehavEmotionRating.txt @@ -0,0 +1,24 @@ +s07t21v1.c3d,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1 +s07t20v2.c3d,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1 +s03t26v1.c3d,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 +s08t12v1.c3d,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 +s05t12v2.c3d,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1 +s11t07v1.c3d,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,1 +s08t17v2.c3d,2,2,2,2,2,2,2,1,3,2,2,3,2,2,3,2,2,2,2 +s10t20v2.c3d,2,2,2,1,2,2,2,2,2,2,2,1,2,2,1,1,2,2,2 +s03t22v2.c3d,2,2,2,2,2,2,2,2,1,2,2,1,2,2,2,2,3,2,2 +s03t28v1.c3d,1,1,2,1,2,2,2,1,2,2,1,1,2,2,1,1,1,4,2 +s03t19v2.c3d,1,1,2,2,2,2,2,2,2,1,2,1,2,2,1,2,3,2,2 +s07t26v3.c3d,3,2,2,1,2,2,2,2,3,1,3,2,2,1,1,1,2,2,2 +s11t20v2.c3d,3,3,3,3,3,3,3,3,3,3,3,3,3,1,3,3,3,3,3 +s03t32v1.c3d,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,3,3 +s10t25v2.c3d,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3 +s07t34v2.c3d,3,3,3,3,3,3,3,3,3,3,3,1,3,3,3,3,3,3,3 +s07t31v2.c3d,3,3,3,3,3,3,3,3,3,3,3,1,3,3,3,3,3,3,3 +s07t35v1.c3d,3,3,1,3,3,3,3,3,3,3,3,3,3,3,3,3,2,3,3 +s10t06v1.c3d,4,4,4,4,4,4,4,3,4,4,4,4,4,4,4,4,4,4,4 +s10t08v1.c3d,4,4,4,4,4,4,1,4,4,4,4,4,4,4,4,4,4,4,4 +s07t03v1.c3d,4,4,4,4,4,4,1,4,4,4,4,4,1,4,4,4,4,4,4 +s07t32v1.c3d,4,3,4,4,3,4,4,3,4,4,4,4,4,4,4,4,4,4,4 +s05t05v2.c3d,4,4,1,4,4,4,3,4,4,4,4,4,4,4,4,4,3,4,4 +s10t04v2.c3d,4,3,4,4,3,4,4,3,4,4,3,4,4,4,4,4,1,4,3 diff --git a/sami_examples/DemoOne_BehavValenceRating.txt b/sami_examples/DemoOne_BehavValenceRating.txt new file mode 100644 index 0000000..2e72169 --- /dev/null +++ b/sami_examples/DemoOne_BehavValenceRating.txt @@ -0,0 +1,24 @@ +s07t21v1.c3d,2,4,5,3,5,4,5,5,3,5,5,5,4,4,4,2,3,3,3 +s07t20v2.c3d,2,5,5,3,3,3,5,4,4,4,4,5,4,2,5,2,3,3,3 +s03t26v1.c3d,2,2,4,3,2,3,5,4,3,3,4,3,4,4,3,4,3,3,2 +s08t12v1.c3d,3,5,5,2,4,4,5,5,4,3,4,2,3,3,4,1,2,4,1 +s05t12v2.c3d,3,2,2,3,2,3,5,4,3,4,4,2,3,3,1,1,3,3,2 +s11t07v1.c3d,2,2,3,2,3,3,5,4,2,4,4,1,3,3,3,3,2,3,1 +s08t17v2.c3d,2,3,2,2,2,5,5,2,-1,2,4,-1,2,4,1,5,3,2,3 +s10t20v2.c3d,2,2,4,2,2,4,4,4,3,3,3,2,3,2,3,1,3,3,2 +s03t22v2.c3d,1,3,3,2,2,2,5,2,1,3,4,3,4,3,1,1,-3,4,2 +s03t28v1.c3d,2,2,3,1,5,3,5,2,3,4,3,2,3,5,3,2,3,-2,2 +s03t19v2.c3d,2,1,3,3,1,2,5,2,2,2,2,1,1,4,0,2,2,3,4 +s07t26v3.c3d,-2,1,3,1,1,1,4,1,-2,0,-2,2,1,2,1,2,3,1,1 +s11t20v2.c3d,-2,-3,-5,-2,-5,-1,-5,-5,-4,-3,-4,-5,-4,2,-5,-2,-3,-5,-4 +s03t32v1.c3d,-2,-2,-5,-4,-5,-2,-5,-5,-3,-4,-4,-5,-4,-5,-5,-5,2,-3,-2 +s10t25v2.c3d,-3,-3,-5,-3,-4,-1,-5,-3,-2,-2,-3,-3,-3,-2,-5,-3,-2,-3,-3 +s07t34v2.c3d,-3,-2,-2,-2,-4,-3,-5,-2,-3,-3,-2,0,-3,-2,-3,-3,-3,-2,-2 +s07t31v2.c3d,-2,-3,-2,-2,-4,-4,-5,-2,-1,-2,-4,0,-2,-4,-3,-1,-3,-1,-2 +s07t35v1.c3d,-4,-1,-3,-2,-3,-2,-4,-1,-2,-4,-3,-4,-2,-3,-2,-2,1,-1,-1 +s10t06v1.c3d,-2,-2,-2,-3,-3,-2,-5,-1,-2,-5,-2,-2,-3,-2,-3,-1,-3,-2,-2 +s10t08v1.c3d,-1,-4,-5,-3,-4,-1,3,-2,-4,-3,-4,-5,-1,-4,-1,-1,-3,-2,-2 +s07t03v1.c3d,-2,-3,-4,-1,-2,-3,2,-3,-3,-3,-3,-4,1,-3,-1,-2,-3,-3,-2 +s07t32v1.c3d,-2,-2,-3,-2,-3,-2,-3,-2,-2,-2,-2,-4,-2,-3,-2,-3,-2,-2,-1 +s05t05v2.c3d,-1,-3,0,-3,-3,-4,-2,-2,-2,-2,-2,-3,-2,-3,-1,-3,-2,-2,-1 +s10t04v2.c3d,-1,-1,-3,-2,-3,-2,-3,-3,-2,-2,0,-2,-4,-3,-2,-2,2,-2,-1 diff --git a/sami_examples/DemoOne_SAMI.m b/sami_examples/DemoOne_SAMI.m new file mode 100644 index 0000000..1cf99f7 --- /dev/null +++ b/sami_examples/DemoOne_SAMI.m @@ -0,0 +1,92 @@ +%% example script +clear all; close all; + +% add toolbox to path +returnHere = pwd; +cd ../sami_toolbox; +addpath(genpath(pwd)); +cd(returnHere) + +% +++ obligatory: loading userOptions / initializing SAMItoolbox +++++++++++++++++++++++++ +userOptions = DemoOne_defineUserOptions(); +userOptions = sami.initSAMI(userOptions,'c'); + +redo = 0; +if redo == 1 + % +++ loading c3d-files and checking them ++++++++++++++++++++++++++++++++++++++++++++ + c3dData = sami.c3d.importFiles(userOptions); + + % +++ calculate features +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + feat_idv = sami.calcFeatures(c3dData, 'idv', userOptions); + feat_itx = sami.calcFeatures(c3dData, 'itx', userOptions); +else + load(fullfile(userOptions.rootPath,'c3dData.mat')); + + load(fullfile(userOptions.rootPath,'features_idv.mat')); + load(fullfile(userOptions.rootPath,'features_itx.mat')); +end + +% +++ create feature_RDMs ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +RDMs_idv = sami.createFeatureRDMs(feat_idv, userOptions); +RDMs_itx = sami.createFeatureRDMs(feat_itx, userOptions); + +% +++ create model RDMs ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +modelThis = {'emotion','valence'}; +RDMs_stim_models_categorical = sami.createCategoryRDMs(modelThis, 'binary', userOptions); + +% +++ create behavioral RDMs +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +[RDMs_behav_val_subj, behav_val_data] = sami.createBehavioralRDMs('DemoOne_BehavValenceRating.txt', 'Valence', 'distance', userOptions); +[RDMs_behav_emo_subj, behav_emo_data] = sami.createBehavioralRDMs('DemoOne_BehavEmotionRating.txt', 'Emotion', 'binary', userOptions); + +% average across subjects +RDMs_behav_val = sami.rdm.averageRDMs(RDMs_behav_val_subj,'mean_val_RDM',[1 0 0]); +RDMs_behav_emo = sami.rdm.averageRDMs(RDMs_behav_emo_subj,'mean_emo_RDM',[1 0 0]); + +%% display RDMs/MDS plots +% +++ show RDMs ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +sami.plotRDMs([RDMs_stim_models_categorical RDMs_idv],userOptions); +sami.plotRDMs([RDMs_stim_models_categorical RDMs_itx],userOptions); +sami.plotRDMs([RDMs_behav_val_subj RDMs_behav_val],userOptions); +sami.plotRDMs([RDMs_behav_emo_subj RDMs_behav_emo],userOptions); + +% +++ MDS of Stimuli +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +sami.MDSofStimuli(RDMs_idv(1), 'emotion', userOptions); +sami.MDSofStimuli(RDMs_idv(1), 'valence', userOptions); + +% +++ 2nd order - RDM correlation matrix and MDS +++++++++++++++++++++++++++++++++++++++++ +sami.RDMsPairwiseCorrelations([RDMs_stim_models_categorical RDMs_idv], userOptions, 'CAT and IDV'); +sami.RDMsPairwiseCorrelations([RDMs_stim_models_categorical RDMs_itx], userOptions, 'CAT and ITX'); +sami.RDMsPairwiseCorrelations([RDMs_behav_val RDMs_behav_emo RDMs_idv], userOptions, 'BEHAV and IDV'); +sami.RDMsPairwiseCorrelations([RDMs_behav_val RDMs_behav_emo RDMs_itx], userOptions, 'BEHAV and ITX'); + +sami.MDSofRDMs([RDMs_stim_models_categorical RDMs_idv], userOptions,'CategoricalModels and IDV'); +sami.MDSofRDMs([RDMs_stim_models_categorical RDMs_itx], userOptions,'CategoricalModels and ITX'); +sami.MDSofRDMs([RDMs_behav_val RDMs_behav_emo RDMs_stim_models_categorical RDMs_idv], userOptions,'Behavior and CategoricalModels and IDV'); +sami.MDSofRDMs([RDMs_behav_val RDMs_behav_emo RDMs_stim_models_categorical RDMs_itx], userOptions,'Behavior and CategoricalModels and ITX'); + +%% do some statistical analyses, like RSA or ANOVA... +% +++ compare METRIC behavioral data +++++++++++++++++++++++++++++++++++++++++++++++++++++ +sami.compareBehavData(behav_val_data, 'emotion', 'Valence', userOptions); + +% +++ compare CATEGORICAL behavioral data ++++++++++++++++++++++++++++++++++++++++++++++++ +sami.compareBehavCategoricalData(behav_emo_data, 'emotion', 'Emotion', userOptions); + +% +++ compare average feature-values between specific categories +++++++++++++++++++++++++ +sami.compareFeatValues(feat_idv, 'emotion', 'IDV', userOptions); +sami.compareFeatValues(feat_idv, 'valence', 'IDV', userOptions); +sami.compareFeatValues(feat_itx, 'emotion', 'ITX', userOptions); +sami.compareFeatValues(feat_itx, 'valence', 'ITX', userOptions); + +% +++ compare categoryRDMs with movement RDMs ++++++++++++++++++++++++++++++++++++++++++++ +sami.compareCatRDMs2FeatRDMs(RDMs_stim_models_categorical,RDMs_idv,'Categorical vs IDV',userOptions); +sami.compareCatRDMs2FeatRDMs(RDMs_stim_models_categorical,RDMs_itx,'Categorical vs ITX',userOptions); + +% +++ compare behavioralRDMs with movement RDMs ++++++++++++++++++++++++++++++++++++++++++ +userOptions.rdms_pairWiseCorr = 'Pearson'; +statsA = sami.compareBehavRDMs2FeatRDMs(RDMs_behav_val_subj, RDMs_idv, 'Valence', 'IDV', userOptions); +statsB = sami.compareBehavRDMs2FeatRDMs(RDMs_behav_val_subj, RDMs_itx, 'Valence', 'ITX', userOptions); + +userOptions.rdms_pairWiseCorr = 'Kendall_taua'; +statsC = sami.compareBehavRDMs2FeatRDMs(RDMs_behav_emo_subj, RDMs_idv, 'Emotion', 'IDV', userOptions); +statsD = sami.compareBehavRDMs2FeatRDMs(RDMs_behav_emo_subj, RDMs_itx, 'Emotion', 'ITX', userOptions); + diff --git a/sami_examples/DemoOne_c3d/s03t19v2.c3d b/sami_examples/DemoOne_c3d/s03t19v2.c3d new file mode 100644 index 0000000..4ed51f1 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s03t19v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s03t22v2.c3d b/sami_examples/DemoOne_c3d/s03t22v2.c3d new file mode 100644 index 0000000..e50d5b9 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s03t22v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s03t26v1.c3d b/sami_examples/DemoOne_c3d/s03t26v1.c3d new file mode 100644 index 0000000..2f32653 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s03t26v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s03t28v1.c3d b/sami_examples/DemoOne_c3d/s03t28v1.c3d new file mode 100644 index 0000000..339d366 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s03t28v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s03t32v1.c3d b/sami_examples/DemoOne_c3d/s03t32v1.c3d new file mode 100644 index 0000000..2fee0b7 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s03t32v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s05t05v2.c3d b/sami_examples/DemoOne_c3d/s05t05v2.c3d new file mode 100644 index 0000000..b785ed4 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s05t05v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s05t12v2.c3d b/sami_examples/DemoOne_c3d/s05t12v2.c3d new file mode 100644 index 0000000..dde5109 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s05t12v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s07t03v1.c3d b/sami_examples/DemoOne_c3d/s07t03v1.c3d new file mode 100644 index 0000000..0bff7fa Binary files /dev/null and b/sami_examples/DemoOne_c3d/s07t03v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s07t20v2.c3d b/sami_examples/DemoOne_c3d/s07t20v2.c3d new file mode 100644 index 0000000..7648cd0 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s07t20v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s07t21v1.c3d b/sami_examples/DemoOne_c3d/s07t21v1.c3d new file mode 100644 index 0000000..c25da4c Binary files /dev/null and b/sami_examples/DemoOne_c3d/s07t21v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s07t26v3.c3d b/sami_examples/DemoOne_c3d/s07t26v3.c3d new file mode 100644 index 0000000..0c2f0ed Binary files /dev/null and b/sami_examples/DemoOne_c3d/s07t26v3.c3d differ diff --git a/sami_examples/DemoOne_c3d/s07t31v2.c3d b/sami_examples/DemoOne_c3d/s07t31v2.c3d new file mode 100644 index 0000000..d385405 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s07t31v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s07t32v1.c3d b/sami_examples/DemoOne_c3d/s07t32v1.c3d new file mode 100644 index 0000000..cac855b Binary files /dev/null and b/sami_examples/DemoOne_c3d/s07t32v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s07t34v2.c3d b/sami_examples/DemoOne_c3d/s07t34v2.c3d new file mode 100644 index 0000000..f41cdc4 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s07t34v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s07t35v1.c3d b/sami_examples/DemoOne_c3d/s07t35v1.c3d new file mode 100644 index 0000000..b2ea783 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s07t35v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s08t12v1.c3d b/sami_examples/DemoOne_c3d/s08t12v1.c3d new file mode 100644 index 0000000..451b836 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s08t12v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s08t17v2.c3d b/sami_examples/DemoOne_c3d/s08t17v2.c3d new file mode 100644 index 0000000..760f0c9 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s08t17v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s10t04v2.c3d b/sami_examples/DemoOne_c3d/s10t04v2.c3d new file mode 100644 index 0000000..50ca591 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s10t04v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s10t06v1.c3d b/sami_examples/DemoOne_c3d/s10t06v1.c3d new file mode 100644 index 0000000..6093b5f Binary files /dev/null and b/sami_examples/DemoOne_c3d/s10t06v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s10t08v1.c3d b/sami_examples/DemoOne_c3d/s10t08v1.c3d new file mode 100644 index 0000000..b58c38c Binary files /dev/null and b/sami_examples/DemoOne_c3d/s10t08v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s10t20v2.c3d b/sami_examples/DemoOne_c3d/s10t20v2.c3d new file mode 100644 index 0000000..e10c39e Binary files /dev/null and b/sami_examples/DemoOne_c3d/s10t20v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s10t25v2.c3d b/sami_examples/DemoOne_c3d/s10t25v2.c3d new file mode 100644 index 0000000..574dd7a Binary files /dev/null and b/sami_examples/DemoOne_c3d/s10t25v2.c3d differ diff --git a/sami_examples/DemoOne_c3d/s11t07v1.c3d b/sami_examples/DemoOne_c3d/s11t07v1.c3d new file mode 100644 index 0000000..2948217 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s11t07v1.c3d differ diff --git a/sami_examples/DemoOne_c3d/s11t20v2.c3d b/sami_examples/DemoOne_c3d/s11t20v2.c3d new file mode 100644 index 0000000..8f3b3d0 Binary files /dev/null and b/sami_examples/DemoOne_c3d/s11t20v2.c3d differ diff --git a/sami_examples/DemoOne_defineUserOptions.m b/sami_examples/DemoOne_defineUserOptions.m new file mode 100644 index 0000000..b395839 --- /dev/null +++ b/sami_examples/DemoOne_defineUserOptions.m @@ -0,0 +1,168 @@ +function userOptions = DemoOne_defineUserOptions() +% defineUserOptions is a nullary function which initialises a struct +% containing the preferences and details for a particular project. +% It should be edited to taste before a project is run, and a new +% one created for each substantially different project. +% +% For a guide to how to fill out the fields in this file, consult +% the documentation folder (particularly the userOptions_guide.m) +% +% A. Zabicki 09-2020 +%__________________________________________________________________________ + +userOptions.debug = true; % will save and/or display some more information during execution of several functions + +%% ********************************************************** +% Project details +% ********************************************************** +% *** edit how needed ********************* + +% This name identifies a collection of files which all belong to the same run of a project. +userOptions.analysisName = 'DemoOne_sami'; + +% The path leading to where the c3d files are stored. +userOptions.c3dPath = fullfile(pwd,'DemoOne_c3d'); + +% *** no need to change here anything ******************************* +% This is the root directory of the project. +userOptions.rootPath = fullfile(pwd,userOptions.analysisName); +% ******************************************************************* + + +%% ********************************************************** +% c3d file and label settings +% ********************************************************** + +% if "Vicon Plug-in Gait" Modell is used and markers are named according to Plug-in Gait +% and _c3d_personIdentifier, sami_toolbox will transform data automatically +userOptions.c3d_ViconPluginGait = false; + +% if markers are available, but named differently, use this to rename labels of "ownMarker" +% into "samiMarker" [Head, LSHO, LELB, LWRI, LHIP, LKNE, LANK, RSHO, RELB, RWRI, RHIP, RKNE, RANK] +userOptions.c3d_OwnMarker = true; +userOptions.c3d_MarkerMatching = {... + 'ownMarker','samiMarker';... + 'HAND','WRI';... + }; +% if specified, user is able to keep own indiviuum identifications, +% otherwise markers will be renamed, e.g. 'p1HEAD' and 'p2HEAD", automatically +% userOptions.c3d_personIdentifier = '[[marker]][[person]]'; + + +%% ********************************************************** +% descriptions of stimuli: regarding their categories, and how to sort them +% ********************************************************** +% set category which is used to sort the stimuli according to +userOptions.stimuli_sorting = 'emotion'; + +% if filename is specified: "stimulus_settings" will be loaded from this file +userOptions.stimuli_settings_filename = 'DemoOne_stimuli_settings.txt'; +% else: stimulus_settings have to be defined here +userOptions.stimuli_settings = {}; + +% providing labels for each category describing the stimuli +userOptions.stimuli_naming_key(1).name = 'emotion'; +userOptions.stimuli_naming_key(1).condition = {'happiness','affection','sadness','anger'}; +userOptions.stimuli_naming_key(1).color = {[0 .5 1],[.1 .8 .1],[1 .5 0],[1 0 0]}; +userOptions.stimuli_naming_key(2).name = 'valence'; +userOptions.stimuli_naming_key(2).condition = {'positive','negative'}; +userOptions.stimuli_naming_key(2).color = {[0 .5 1],[1 .5 0]}; + +% stimuli_labels for MDS plots +% ---> !!! same order as in first column in userOptions.stimuli_settings_filename !!! +% [userOptions.stimuli_MDS_labels{1:size(userOptions.stimuli_settings,1)-1}] = deal(' '); + +%% ********************************************************** +% calculating individual/interaction movement-features +% -> comparing movement-features between stimuli-categories +% -> creating movement/feature-RDMs for stimulus-set +% ********************************************************** +% set alpha and correction method which will be applied in post-hoc multiple comparisons of feature values +userOptions.feat_threshold = 0.05; % default: 0.05 +userOptions.feat_multipleTesting = 'bonferroni'; % ['bonferroni'] | 'tukey-kramer' | 'hsd' | 'lsd' | 'dunn-sidak' | 'scheffe' + +% which distance measure to use when calculating feature-RDMs. +userOptions.feat_distance = 'euclidean'; % input into pdist function, default: 'euclidean' + +% should feature-RDM-entries be rank transformed into [0,1] before they're displayed? +userOptions.feat_rankTransform = false; % !!!!!!!! NOT YET IMPLEMENTED !!!!!!!!!!!! + +%% ********************************************************** +% behavioral data: +% -> comparing average behavioral ratings between stimuli-categories +% -> creating behavioral-RDMs for each subject +% ********************************************************** +% set alpha and correction method which will be applied in post-hoc multiple comparisons of behavioral ratings +userOptions.behav_threshold = 0.05; % default: 0.05 +userOptions.behav_multipleTesting = 'bonferroni'; % ['bonferroni'] | 'tukey-kramer' | 'hsd' | 'lsd' | 'dunn-sidak' | 'scheffe' + +% which distance measure to use between stimuli-rating-values when calculating subjects behavioral-RDMs +userOptions.behav_distanceMeasure = 'euclidean'; % input into pdist function, default: 'euclidean' + +%% ********************************************************** +% MDS plots +% ********************************************************** +% what criterion to be minimised in MDS calculation? +userOptions.MDS_criterion = 'metricstress'; % default: 'metricstress' + +% style +userOptions.MDS_plotLabels = true; % show labels in MDS-Plot +userOptions.MDS_plotLegend = true; % show color-legend in MDS-plot -> useful for "stimuliMDS()" +userOptions.MDS_dotSize = 20; % default: 20 +userOptions.MDS_fontSize = 9; % default: 9 + +%% ********************************************************** +% second-order-analysis of RDMs +% ********************************************************** +% which similarity-measure is used for the pair-wise comparison of RDMs +userOptions.rdms_pairWiseCorr = 'Kendall_taua'; + +% for 'compareCatRDMs2FeatRDMs' function +% set alpha and correction method to be applied in pairwise-RDM-correlation analysis +userOptions.rdms_pairWiseCorrThreshold = 0.05; % default: 0.05 +userOptions.rdms_pairWiseCorrMultipleTesting = 'holm'; % ['holm'] | 'bonferroni' | 'FDR' + +% for 'compareBehavRDMs2FeatRDMs' function +% set test, alpha and correction method for analysing relatedness of features and behavioral RDMs +userOptions.rdms_relatednessTest = 'signedRank'; % ['signedRank'] | 'randomisation' +userOptions.rdms_relatednessThreshold = 0.05; % default: 0.05 +userOptions.rdms_relatednessMultipleTesting = 'FDR'; % ['FDR'] | 'bonferroni' | 'holm' + +% set test, alpha and correction method for analysing relatedness of features and behavioral RDMs +userOptions.rdms_differencesTest = 'signedRank'; % ['signedRank'] | 'conditionBootstrap' +userOptions.rdms_differencesThreshold = 0.05; % default: 0.05 +userOptions.rdms_differencesMultipleTesting = 'FDR'; % ['FDR'] | 'bonferroni' | 'holm' + +% some set test, alpha and correction method for analysing relatedness of features and behavioral RDMs +userOptions.rdms_orderByCorr = true; % sort features by height of relatedness? default: true +userOptions.rdms_nRandomisations = 50000; % default: 50000 (min. 10,000 highly recommended) +userOptions.rdms_nBootstrap = 1000; % default: 1000 + +%% ********************************************************** +% everything else, my be set or not +% (here, defaults are used to show what can be edited) +% ********************************************************** +userOptions.default_modelRDMcolor = [.5 .5 .5]; +userOptions.default_behavRDMcolor = [.9 .2 .1]; + +% threshold: minimal number of datapoints for parametric group-comparisons (ANOVA, t-test) +userOptions.stats_minNforSubjectRFXtests = 12; +% threshold: minimal number of datapoints for NON-parametric group-comparisons (Wilcoxon, Kruskall-Walis) +userOptions.stats_minNforNonParamTests = 5; + +%% ********************************************************** +% handling figures +% ********************************************************** +% generall displaying all figures? +userOptions.fig_display = true; % default: true + +% How should figures be outputted? +userOptions.fig_saveFIG = true; % default: false +userOptions.fig_savePDF = true; % default: false +userOptions.fig_saveSVG = true; % default: false +userOptions.fig_saveTIF = true; % default: false + +% Which dots per inch resolution do we output? +userOptions.fig_dpi = 300; + +end diff --git a/sami_examples/DemoOne_stimuli_settings.txt b/sami_examples/DemoOne_stimuli_settings.txt new file mode 100644 index 0000000..0e259cd --- /dev/null +++ b/sami_examples/DemoOne_stimuli_settings.txt @@ -0,0 +1,25 @@ +filename,emotion,valence +s07t21v1.c3d,1,1 +s07t20v2.c3d,1,1 +s03t26v1.c3d,1,1 +s08t12v1.c3d,1,1 +s05t12v2.c3d,1,1 +s11t07v1.c3d,1,1 +s08t17v2.c3d,2,1 +s10t20v2.c3d,2,1 +s03t22v2.c3d,2,1 +s03t28v1.c3d,2,1 +s03t19v2.c3d,2,1 +s07t26v3.c3d,2,1 +s11t20v2.c3d,3,2 +s03t32v1.c3d,3,2 +s10t25v2.c3d,3,2 +s07t34v2.c3d,3,2 +s07t31v2.c3d,3,2 +s07t35v1.c3d,3,2 +s10t06v1.c3d,4,2 +s10t08v1.c3d,4,2 +s07t03v1.c3d,4,2 +s07t32v1.c3d,4,2 +s05t05v2.c3d,4,2 +s10t04v2.c3d,4,2 diff --git a/sami_toolbox/+sami/+c3d/getPersonIdentifier.m b/sami_toolbox/+sami/+c3d/getPersonIdentifier.m new file mode 100644 index 0000000..454c92b --- /dev/null +++ b/sami_toolbox/+sami/+c3d/getPersonIdentifier.m @@ -0,0 +1,114 @@ +function [personIdentifier,personIndex,nPersons] = getPersonIdentifier(c3d, samiOptions) +% [personIdentifier,personIndex,nPersons] = getPersonIdentifier(c3d, samiOptions) +% +% Identifies individuals from given marker labels. and returns their indi +% Detects which marker belongs to which person, and labels the given marker according to +% sami_format for each person (i.e. p1HEAD). +% +% input: +% - c3d: +% struct. Containing c3d data of a stimulus, given by importFiles function. +% +% - samiOptions.labels: +% cell array. Defines sami_standard names and order for labels. Needed to +% identify individuals by searching for substring unequal to sami_labels. +% +% output: +% - personIdentifier: +% cell array. contains individual substrings within labels for each identified person. +% +% - personIndex: +% array. for each label, personIndex states the corresponding person. +% +% - nPersons: +% double. how many persons are identified in given labels. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 10/2020 + +% init vars +nLabelsUnsorted = numel(c3d.LabelsUnsorted); +personIndex = nan(1,nLabelsUnsorted); + +% various versions for identifying individuals from labels implemented, #3 seems to catch +% all possible naming_conventions +use_version = 3; + +switch use_version + case 1 + % detect which marker belongs to an individuum, by erase 'marker-name' +++++++++++ + % from labels, and then finding unique values + personIdentifierTMP = cell(nLabelsUnsorted,1); + for i = 1:nLabelsUnsorted + fullLabel = c3d.LabelsUnsorted{i}; + if any(cellfun(@(s) contains(fullLabel, s), samiOptions.labels)) + markerLabel = samiOptions.labels{cellfun(@(s) contains(fullLabel, s), samiOptions.labels)}; + personIdentifierTMP{i} = erase(fullLabel,markerLabel); + end + end + personIdentifierTMP = personIdentifierTMP(~cellfun('isempty',personIdentifierTMP)); + personIdentifier = unique(personIdentifierTMP); + + case 2 + % v2: ausgehend von samiLabels + personIdentifierTMP = cell(0); + for i = 1:numel(samiOptions.labels) + samiLabel = samiOptions.labels{i}; + labelIdx = find(contains(c3d.LabelsUnsorted, samiLabel)); + + for L = 1:numel(labelIdx) + fullLabel = c3d.LabelsUnsorted{labelIdx(L)}; + personIdentifierTMP = [personIdentifierTMP; {erase(fullLabel,samiLabel)}]; %#ok + end + end + personIdentifier = unique(personIdentifierTMP); + + case 3 + % first: find individuel 'personIdentifiers' in labelsUnsorted by looking into + % marker which contain "samiLabels" + personIdentifierTMP = cell(nLabelsUnsorted,1); + for i = 1:nLabelsUnsorted + fullLabel = c3d.LabelsUnsorted{i}; + tmp = cellfun(@(s) contains(fullLabel, s), samiOptions.labels); + if any(tmp) + markerLabel = samiOptions.labels{tmp}; + personIdentifierTMP{i} = erase(fullLabel,markerLabel); + end + end + personIdentifierTMP = personIdentifierTMP(cellfun(@(s) ischar(s),personIdentifierTMP)); + personIdentifier = unique(personIdentifierTMP); + + % second: assign to each unsortedLabel person# (i.e. index in 'personIdentifier') + emptyIdentity = cellfun('isempty',personIdentifier); + if ~any(emptyIdentity) + for i = 1:numel(emptyIdentity) + % loop this, bc there might be more than 2 persons in the future + personIndex(contains(c3d.LabelsUnsorted,personIdentifier(i))) = i; + end + else + % case that one person is identified by 'pure labels', i.e. no suffix to the + % marker was added, e.g. subject 1 -> 'Head' / subject 2 -> 'Head2' + + % do this: loop all non-empty-identifiers, set 'personIndex' and remember which index + % was changed, use this knowledge to identify emptyIdentity-labels + nonEmptyIdx = zeros(1,nLabelsUnsorted); + for i = find(~emptyIdentity) + thisNonEmptyIdx = contains(c3d.LabelsUnsorted,personIdentifier(i)); + nonEmptyIdx = nonEmptyIdx + double(thisNonEmptyIdx); + personIndex(thisNonEmptyIdx) = i; + end + + % set emptyIdentifier there, where no previous non-empty-identifier was found + emptyIdx = nonEmptyIdx==0; + personIndex(emptyIdx) = find(emptyIdentity); + end +end + +% find numer of persons ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +nPersons = numel(personIdentifier); +if nPersons < 1 || nPersons > 2 + fprintf('\n'); error('*** ERROR *** there are less then 1 or more then 2 persons identified in c3d-file.'); +end + +end + diff --git a/sami_toolbox/+sami/+c3d/importFiles.m b/sami_toolbox/+sami/+c3d/importFiles.m new file mode 100644 index 0000000..add516b --- /dev/null +++ b/sami_toolbox/+sami/+c3d/importFiles.m @@ -0,0 +1,151 @@ +function c3dData = importFiles(userOptions) +% c3dData = importFiles(userOptions) +% +% Imports c3d-files from a specific folder, set in userOptions, and returns struct +% containing all relevant 3D-MoCap-data for further analysis. If set, stimuli will be +% sorted according to 'userOptions.stimuli_sorting', otherwise alphabetically. +% +% input: +% - userOptions.c3dPath: +% string. Path to folder which contains all c3d files to be loaded. +% +% - userOptions.c3d_ViconPluginGait: +% boolean. In case c3d-MoCap-Data was recorded using Vicon "Plug-in Gait Full Body" +% markerset, the data will be transformed automaticaly into sami-format. +% Defaults to false. +% +% (https://docs.vicon.com/display/Nexus25/Plug-in+Gait+models+and+templates) +% +% - userOptions.c3d_OwnMarker: +% boolean. If all sami-markers are available, but named differently, this function +% can rename parts of given labels, if userOptions.c3d_OwnMarker is true. +% Defaults to false. +% +% - userOptions.c3d_MarkerMatching: +% cell array. Definies rules for renaming labels. In each row, the substring +% from the 'ownMarker# column will be replaced the string in the 'samiMarker' column. +% Example: +% userOptions.c3d_MarkerMatching = {'ownMarker','samiMarker';... +% 'Head','HEAD';... +% 'LASI','LHIP';... +% 'RASI','RHIP';... +% };% +% +% output: +% - c3dData: +% struct. Containing for each loaded c3dd-stimulis-file relevant data: +% c3dData.file: filename +% c3dData.marker: 3d data in [time*marker*dimension] format +% c3dData.labels: labels of [marker] +% c3dData.framerate: framerate in Hz +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 10/2020 + +%% define default behavior +if ~isfield(userOptions, 'c3dPath'), error('*** sami:error *** path to folder containing c3d data mus be set in userOptions.c3dPath.'); end%if +userOptions = sami.util.setIfUnset(userOptions, 'c3d_ViconPluginGait', false); +userOptions = sami.util.setIfUnset(userOptions, 'c3d_OwnMarker', false); +if userOptions.c3d_OwnMarker==true && (~isfield(userOptions, 'c3d_MarkerMatching') || isempty(userOptions.c3d_MarkerMatching) || size(userOptions.c3d_MarkerMatching,1)<2), error('*** sami:error *** please specifiy ''userOptions.c3d_MarkerMatching''.'); end%if + +disp('*** importing c3d files and checking them ***'); + +%% load samiOptions +samiOptions = sami.loadSamiOptions(); +labelOrderWanted = [strcat('p1',samiOptions.labels), strcat('p2',samiOptions.labels)]; +nLabels = numel(labelOrderWanted); + +%% read in data folder +files = dir(fullfile(userOptions.c3dPath,'*.c3d')); + +%% get order for Stimuli +[~,fileOrderNames] = sami.util.getStimOrder(userOptions); + +%% loop files + load c3d + checking +c3dData = struct(); +nFiles = numel(files); +fprintf(' ... import file (of %d): 1',nFiles); +for f = 1:nFiles + % display progress + fprintf([repmat('\b',1,numel(num2str(f-1))) '%d'],f); + + %% load c3d -> according to sort_order defined by userOptions.stimuli_sorting + thisFileIdx = find(contains({files(:).name},fileOrderNames{f})); + thisFileName = fullfile(files(thisFileIdx).folder,files(thisFileIdx).name); + [c3d.MarkersUnsorted, c3d.VideoFrameRate,~,~,~, c3d.ParameterGroup,~,~,c3d.MissingMarker] = sami.c3d.loadc3d(thisFileName); + + % find labels + idxPoint = find(contains([c3d.ParameterGroup.name],'POINT')); + idxLabels = contains([c3d.ParameterGroup(idxPoint).Parameter.name],'LABELS'); + c3d.LabelsUnsorted = c3d.ParameterGroup(idxPoint).Parameter(idxLabels).data; + + % check for missing data (residuals==-1 in c3d-file) + if any(c3d.MissingMarker(:)) + fprintf('\n'); error('*** sami:ERROR *** missing data. please label all markers without any gaps.'); + end + + %% if Vicon Plug-in Gait Modell is used, call convert_function + if userOptions.c3d_ViconPluginGait == true + c3d = sami.c3d.plugingait(c3d, samiOptions); + else + %% else, rename/sort into target sami_labels/sami_order + nLabelsUnsorted = numel(c3d.LabelsUnsorted); + + % if ownMarker names were used, rename them to sami-standard +++++++++++++++++++++ + if userOptions.c3d_OwnMarker == true + for i = 1:nLabelsUnsorted + oldLabel = c3d.LabelsUnsorted{i}; + if any( cellfun(@(s) contains(oldLabel, s), userOptions.c3d_MarkerMatching(:,1)) ) + idx = cellfun(@(s) contains(oldLabel, s), userOptions.c3d_MarkerMatching(:,1)); + c3d.LabelsUnsorted{i} = strrep(oldLabel, userOptions.c3d_MarkerMatching{idx,1} , userOptions.c3d_MarkerMatching{idx,2}); + end + end + end + + % rename labels by adding person# to sami-standard + c3d = sami.c3d.renameLabels2samiStandard(c3d,samiOptions); + + % sort labels into order defined by sami-standard ++++++++++++++++++++++++++++++++ + c3d.Labels = cell(1,nLabels); + c3d.Markers = nan(size(c3d.MarkersUnsorted,1),nLabels,size(c3d.MarkersUnsorted,3)); + for i = 1:nLabels + tmp_idx = find(strcmp(c3d.LabelsUnsorted , labelOrderWanted{i})); + if ~isempty(tmp_idx) + c3d.Labels{i} = c3d.LabelsUnsorted{tmp_idx}; + c3d.Markers(:,i,:) = c3d.MarkersUnsorted(:,tmp_idx,:); + end + end + + % check for any missing lables from config_labels ++++++++++++++++++++++++++++++++ + if any(cellfun(@isempty,c3d.Labels)) + fprintf('\n'); error('*** ERROR *** not all needed labels are found in c3d-file.'); + end + end + + %% Quality-Checks + % NANs? + if any(isnan(c3d.Markers(:))) + fprintf('\n'); error('*** ERROR *** NANs still remain in 3d-MoCap-Data. Please check input data! abort.'); + end + + %% saving structure + c3dData(f).file = files(thisFileIdx).name; + c3dData(f).marker = c3d.Markers; + c3dData(f).labels = c3d.Labels; + c3dData(f).frameRate = c3d.VideoFrameRate; + c3dData(f).nPersons = c3d.nPersons; + c3dData(f).ignore_labels = c3d.LabelsUnsorted; + c3dData(f).ignore_marker = c3d.MarkersUnsorted; +end +fprintf(' ... ok\n'); + +%% saving c3dData struct +save(fullfile(userOptions.rootPath,'c3dData.mat'),'c3dData'); + +%% finishing +fprintf(' ... c3d_data saved in rootPath as "c3dData.mat" \n'); +fprintf(' ... DONE loading c3d_data\n\n'); + +end + diff --git a/sami_toolbox/+sami/+c3d/loadc3d.m b/sami_toolbox/+sami/+c3d/loadc3d.m new file mode 100644 index 0000000..1027cce --- /dev/null +++ b/sami_toolbox/+sami/+c3d/loadc3d.m @@ -0,0 +1,301 @@ +function [Markers, VideoFrameRate, AnalogSignals, AnalogFrameRate, Event, ParameterGroup, CameraInfo, ResidualError, MissingMarker] = loadc3d(FullFileName) +% LOADC3D - Retrieves 3D coordinate/analog data from a C3D file +% +% Input: +% FullFileName File name of trial +% showWaitBar showing progress bar? (true/[false]) +% +% Outputs: +% Markers 3D-marker data [Nmarkers x NvideoFrames x Ndim(=3)] +% VideoFrameRate Frames/sec +% AnalogSignals Analog signals [Nsignals x NanalogSamples] +% AnalogFrameRate Samples/sec +% Event Event(Nevents).time ..value ..name +% ParameterGroup ParameterGroup(Ngroups).Parameters(Nparameters).data ..etc. +% CameraInfo MarkerRelated CameraInfo [Nmarkers x NvideoFrames] +% ResidualError MarkerRelated ErrorInfo [Nmarkers x NvideoFrames] +% +% AUTHOR(S) AND VERSION-HISTORY +% Ver. 1.0 Creation (Alan Morris, Toronto, October 1998) [originally named "getc3d.m"] +% Ver. 2.0 Revision (Jaap Harlaar, Amsterdam, April 2002) [renamed to "readc3d.m"] +% Ver. 2.1 Revision (Jasper Menger, Enschede, February 2005) [adapted for files with Float32-data, speed improved, renamed to "loadc3d"] +% Ver. 2.2 Revision (Adam Zabicki, Gießen, June 2020) [get rid of "progress_bar", it slowed down performance dramatically! +% also, add new *MissingMarker" variable, indicating "Residual = -1" timepoints ] +% ________________________________________________________________________________________ + +%% init vars +Markers = []; +VideoFrameRate = 0; +AnalogSignals = []; +AnalogFrameRate = 0; +Event = []; +ParameterGroup = []; +CameraInfo = []; +ResidualError = []; +MissingMarker = []; + +% ############################################### +% ## ## +% ## open the file ## +% ## ## +% ############################################### + +ind = strfind(FullFileName,'\'); +if ind>0 + FileName = FullFileName(ind(length(ind)) + 1:length(FullFileName)); +else + FileName = FullFileName; +end + +fid = fopen(FullFileName, 'r', 'n'); % native format (PC-intel) +if fid == -1 + h = errordlg(['File: ', FileName, ' could not be opened'], 'application error'); + uiwait(h) + return +end + +NrecordFirstParameterblock = fread(fid, 1, 'int8'); % Reading record number of parameter section +key = fread(fid, 1, 'int8'); % key = 80; +if key ~= 80 + h = errordlg(['File: ',FileName,' does not comply to the C3D format'], 'application error'); + uiwait(h) + fclose(fid); + return +end + +fseek(fid, 512 * (NrecordFirstParameterblock - 1) + 3, 'bof'); % jump to processortype - field +proctype = fread(fid, 1, 'int8') - 83; % proctype: 1(INTEL-PC); 2(DEC-VAX); 3(MIPS-SUN/SGI) +if proctype == 2 + fclose(fid); + fid = fopen(FullFileName, 'r', 'd'); % DEC VAX D floating point and VAX ordering +end + + +% ############################################### +% ## ## +% ## read header ## +% ## ## +% ############################################### + +% NrecordFirstParameterblock = fread(fid, 1, 'int8'); % Reading record number of parameter section +% key1 = fread(fid, 1, 'int8'); % key = 80; + +fseek(fid, 2, 'bof'); +Nmarkers = fread(fid, 1, 'int16'); % number of markers +NanalogSamplesPerVideoFrame = fread(fid, 1, 'int16'); % number of analog channels x #analog frames per video frame +StartFrame = fread(fid, 1, 'int16'); % # of first video frame +EndFrame = fread(fid, 1, 'int16'); % # of last video frame +MaxInterpolationGap = fread(fid,1, 'int16'); % maximum interpolation gap allowed (in frame) +Scale = fread(fid, 1, 'float32'); % floating-point scale factor to convert 3D-integers to ref system units +NrecordDataBlock = fread(fid, 1, 'int16'); % starting record number for 3D point and analog data + +NanalogFramesPerVideoFrame = fread(fid, 1, 'int16'); +if NanalogFramesPerVideoFrame > 0 + NanalogChannels = NanalogSamplesPerVideoFrame / NanalogFramesPerVideoFrame; +else + NanalogChannels = 0; +end +VideoFrameRate = fread(fid, 1, 'float32'); +AnalogFrameRate = VideoFrameRate * NanalogFramesPerVideoFrame; + + +% ############################################### +% ## ## +% ## read events ## +% ## ## +% ############################################### + +fseek(fid, 298, 'bof'); +EventIndicator = fread(fid, 1, 'int16'); +if EventIndicator == 12345 + Nevents = fread(fid, 1, 'int16'); + fseek(fid, 2, 'cof'); % skip one position/2 bytes + if Nevents > 0 + for i = 1:Nevents + Event(i).time = fread(fid, 1, 'float'); + end + fseek(fid, 188 * 2, 'bof'); + for i = 1:Nevents + Event(i).value = fread(fid, 1, 'int8'); + end + fseek(fid, 198 * 2, 'bof'); + for i = 1:Nevents + Event(i).name = cellstr(char(fread(fid, 4, 'char')')); + end + end +end + + +% ############################################### +% ## ## +% ## read 1st parameter block ## +% ## ## +% ############################################### + +fseek(fid, 512 * (NrecordFirstParameterblock - 1), 'bof'); +dat1 = fread(fid, 1, 'int8'); +key2 = fread(fid, 1, 'int8'); % key = 80; +NparameterRecords = fread(fid, 1, 'int8'); +proctype = fread(fid, 1, 'int8')-83; % proctype: 1(INTEL-PC); 2(DEC-VAX); 3(MIPS-SUN/SGI) +Ncharacters = fread(fid, 1, 'int8'); % characters in group/parameter name +GroupNumber = fread(fid, 1, 'int8'); % id number -ve=group / +ve=parameter + +while Ncharacters > 0 % The end of the parameter record is indicated by <0 characters for group/parameter name + + if GroupNumber < 0 % Group data + GroupNumber = abs(GroupNumber); + GroupName = fread(fid, [1, Ncharacters], 'char'); + ParameterGroup(GroupNumber).name = cellstr(char(GroupName)); %group name + offset = fread(fid, 1, 'int16'); %offset in bytes + deschars = fread(fid, 1, 'int8'); %description characters + GroupDescription = fread(fid, [1, deschars], 'char'); + ParameterGroup(GroupNumber).description = cellstr(char(GroupDescription)); %group description + ParameterNumberIndex(GroupNumber) = 0; + fseek(fid, offset - 3 - deschars, 'cof'); + + else % Parameter data + clear dimension; + ParameterNumberIndex(GroupNumber) = ParameterNumberIndex(GroupNumber) + 1; + ParameterNumber = ParameterNumberIndex(GroupNumber); % index all parameters within a group + ParameterName = fread(fid, [1, Ncharacters], 'char'); % name of parameter + + % read parameter name + if size(ParameterName) > 0 + ParameterGroup(GroupNumber).Parameter(ParameterNumber).name = cellstr(char(ParameterName)); %save parameter name + end + + % read offset + offset = fread(fid, 1, 'int16'); % offset of parameters in bytes + filepos = ftell(fid); % present file position + nextrec = filepos + offset(1) - 2; % position of beginning of next record + + % read type + type = fread(fid, 1, 'int8'); % type of data: -1=char/1=byte/2=integer*2/4=real*4 + ParameterGroup(GroupNumber).Parameter(ParameterNumber).datatype = type; + + % read number of dimensions + dimnum = fread(fid, 1, 'int8'); + if dimnum == 0 + datalength = abs(type); % length of data record + else + mult = 1; + for j = 1:dimnum + dimension(j) = fread(fid, 1, 'int8'); + mult = mult * dimension(j); + ParameterGroup(GroupNumber).Parameter(ParameterNumber).dim(j) = dimension(j); %save parameter dimension data + end + datalength = abs(type) * mult; %length of data record for multi-dimensional array + end + + if type == -1 % datatype=='char' + + wordlength=dimension(1); % length of character word + if dimnum == 2 && datalength > 0 % & parameter(idnumber, index,2).dim > 0 + for j= 1 : dimension(2) + data = fread(fid, [1, wordlength], 'char'); %character word data record for 2-D array + ParameterGroup(GroupNumber).Parameter(ParameterNumber).data(j) = cellstr(char(data)); + end + elseif dimnum == 1 && datalength > 0 + data = fread(fid, [1, wordlength], 'char'); % numerical data record of 1-D array + ParameterGroup(GroupNumber).Parameter(ParameterNumber).data = cellstr(char(data)); + end + + elseif type == 1 %1-byte for boolean + + Nparameters=datalength/abs(type); + data = fread(fid,Nparameters,'int8'); + ParameterGroup(GroupNumber).Parameter(ParameterNumber).data = data; + + elseif type == 2 && datalength > 0 %integer + + Nparameters = datalength / abs(type); + data = fread(fid, Nparameters, 'int16'); + if dimnum > 1 + ParameterGroup(GroupNumber).Parameter(ParameterNumber).data = reshape(data, dimension); + else + ParameterGroup(GroupNumber).Parameter(ParameterNumber).data = data; + end + + elseif type == 4 && datalength > 0 + + Nparameters = datalength / abs(type); + data = fread(fid, Nparameters, 'float'); + if dimnum > 1 + ParameterGroup(GroupNumber).Parameter(ParameterNumber).data = reshape(data,dimension); + else + ParameterGroup(GroupNumber).Parameter(ParameterNumber).data = data; + end + else + % error + end + + deschars = fread(fid, 1, 'int8'); % description characters + if deschars > 0 + description = fread(fid, [1, deschars], 'char'); + ParameterGroup(GroupNumber).Parameter(ParameterNumber).description = cellstr(char(description)); + end + %moving ahead to next record + fseek(fid, nextrec, 'bof'); + end + + % check group/parameter characters and idnumber to see if more records present + Ncharacters = fread(fid, 1, 'int8'); % characters in next group/parameter name + GroupNumber = fread(fid, 1, 'int8'); % id number -ve=group / +ve=parameter +end + + +% ############################################### +% ## ## +% ## read data block ## +% ## ## +% ############################################### +% Get the coordinate and analog data + +fseek(fid, (NrecordDataBlock - 1) * 512,'bof'); +NvideoFrames = EndFrame - StartFrame + 1; +Markers = zeros(NvideoFrames, Nmarkers, 3); +AnalogSignals = zeros((NanalogFramesPerVideoFrame * NvideoFrames), NanalogChannels); +CameraInfo = zeros(NvideoFrames, Nmarkers); +ResidualError = zeros(NvideoFrames, Nmarkers); +MissingMarker = nan(NvideoFrames, Nmarkers); + +if Scale < 0 + % Negative scale factor: Float32-data + for i = 1:NvideoFrames %XXXXXXXXXXXXXX hier flaschenhals XXXXXXXXX + % Marker and AD data + DataM = fread(fid, [4, Nmarkers], 'float32')'; + DataA = fread(fid, [NanalogChannels, NanalogFramesPerVideoFrame], 'float32')'; + Markers(i, 1:Nmarkers, 1:3) = DataM(1:Nmarkers, 1:3); + i_start = NanalogFramesPerVideoFrame * (i - 1) + 1; + i_stop = NanalogFramesPerVideoFrame * i; + AnalogSignals(i_start:i_stop, 1:NanalogChannels) = DataA; + % Camera Info and Residual Error + Byte(1:Nmarkers) = fix(squeeze(DataM(1:Nmarkers, 4))); + HighByte = fix(Byte / 256); + LowByte = Byte - HighByte * 256; + CameraInfo(i, 1:Nmarkers) = HighByte(1:Nmarkers); + ResidualError(i, 1:Nmarkers) = LowByte(1:Nmarkers) * abs(Scale); + % mark missing marker if residual == -1 + missing = squeeze(DataM(1:Nmarkers, 4)) == -1; + MissingMarker(i, missing) = 1; + MissingMarker(i, ~missing) = 0; + end +else + % Positive scale factor: Int-data + for i = 1:NvideoFrames + for j = 1:Nmarkers + Markers(i , j, 1:3) = fread(fid , 3, 'int16')' .* Scale; + ResidualError(i,j) = fread(fid, 1, 'int8'); + CameraInfo(i,j) = fread(fid, 1, 'int8'); + end + DataA = fread(fid, [NanalogFramesPerVideoFrame, NanalogChannels], 'int16'); + i_start = NanalogFramesPerVideoFrame * (i - 1) + 1; + i_stop = NanalogFramesPerVideoFrame * i; + AnalogSignals(i_start:i_stop, 1:NanalogChannels) = DataA; + end +end + +fclose(fid); + +end + diff --git a/sami_toolbox/+sami/+c3d/plugingait.m b/sami_toolbox/+sami/+c3d/plugingait.m new file mode 100644 index 0000000..1cf289e --- /dev/null +++ b/sami_toolbox/+sami/+c3d/plugingait.m @@ -0,0 +1,90 @@ +function c3d = plugingait(c3d, samiOptions) +% c3d = plugingait(c3d, samiOptions) +% +% This function receives 3d-MoCap as defined by the VICON NEXUS "full body plugin-gait" +% model. It expects two posterior superior iliac spine (PSIS) markers for the pelvis (not +% the single sacral % marker), and does not use KAD. It then transforms the current +% 3d-MoCap data into the SAMI-13-Markerset. +% +% input: +% - c3d: +% struct. Containing c3d data of a stimulus, given by importFiles function. +% +% - samiOptions.labels: +% cell array. Defines sami_standard names and order for labels. Needed to +% identify individuals by searching for substring unequal to sami_labels. +% string. Path to folder which contains all c3d files to be loaded. +% +% output: +% - c3d: +% struct. same as input, but automatically preprocessed from Plug-in-Gait into +% sami-data-format (e.g. averaging 4-head-markers from plug-in-gait model,...) +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 10/2020 + +%% init varis +labelOrderWanted = [strcat('p1',samiOptions.labels), strcat('p2',samiOptions.labels)]; +nLabels = numel(labelOrderWanted); + +%% rename labels by adding person# to sami-standard +c3d = sami.c3d.renameLabels2samiStandard(c3d,samiOptions); + +%% sort labels into order defined by sami-standard +c3d.Labels = cell(1,nLabels); +c3d.Markers = nan(size(c3d.MarkersUnsorted,1),nLabels,size(c3d.MarkersUnsorted,3)); +for i = 1:nLabels + actualLabel = labelOrderWanted{i}; + actualPerson = actualLabel(1:2); + + % depending on sami_label, do different things + switch actualLabel(3:end) + case 'HEAD' + c3d.Labels{i} = actualLabel; + searchThis = {[actualPerson 'LFHD'],[actualPerson 'LBHD'],[actualPerson 'RFHD'],[actualPerson 'RBHD']}; + idx = contains(c3d.LabelsUnsorted,searchThis); + c3d.Markers(:,i,:) = mean(c3d.MarkersUnsorted(:,idx,:),2); + + case 'LWRI' + c3d.Labels{i} = actualLabel; + searchThis = {[actualPerson 'LWRA'],[actualPerson 'LWRB']}; + idx = contains(c3d.LabelsUnsorted,searchThis); + c3d.Markers(:,i,:) = mean(c3d.MarkersUnsorted(:,idx,:),2); + + case 'RWRI' + c3d.Labels{i} = actualLabel; + searchThis = {[actualPerson 'RWRA'],[actualPerson 'RWRB']}; + idx = contains(c3d.LabelsUnsorted,searchThis); + c3d.Markers(:,i,:) = mean(c3d.MarkersUnsorted(:,idx,:),2); + + case 'LHIP' + c3d.Labels{i} = actualLabel; + searchThis = {[actualPerson 'LASI'],[actualPerson 'LPSI']}; + idx = contains(c3d.LabelsUnsorted,searchThis); + c3d.Markers(:,i,1:2) = mean(c3d.MarkersUnsorted(:,idx,1:2),2); + + searchThis = [actualPerson 'LASI']; + idx = contains(c3d.LabelsUnsorted,searchThis); + c3d.Markers(:,i,3) = c3d.MarkersUnsorted(:,idx,3); + + case 'RHIP' + c3d.Labels{i} = actualLabel; + searchThis = {[actualPerson 'RASI'],[actualPerson 'RPSI']}; + idx = contains(c3d.LabelsUnsorted,searchThis); + c3d.Markers(:,i,1:2) = mean(c3d.MarkersUnsorted(:,idx,1:2),2); + + searchThis = [actualPerson 'RASI']; + idx = contains(c3d.LabelsUnsorted,searchThis); + c3d.Markers(:,i,3) = c3d.MarkersUnsorted(:,idx,3); + + otherwise + tmp_idx = find(strcmp(c3d.LabelsUnsorted , actualLabel)); + if ~isempty(tmp_idx) + c3d.Labels{i} = c3d.LabelsUnsorted{tmp_idx}; + c3d.Markers(:,i,:) = c3d.MarkersUnsorted(:,tmp_idx,:); + end + end +end + +end + diff --git a/sami_toolbox/+sami/+c3d/renameLabels2samiStandard.m b/sami_toolbox/+sami/+c3d/renameLabels2samiStandard.m new file mode 100644 index 0000000..629181a --- /dev/null +++ b/sami_toolbox/+sami/+c3d/renameLabels2samiStandard.m @@ -0,0 +1,39 @@ +function c3d = renameLabels2samiStandard(c3d,samiOptions) +% c3d = renameLabels2samiStandard(c3d,samiOptions) +% +% Detects which marker belongs to which person, and labels the given marker according to +% sami_format for each person (i.e. p1HEAD). +% +% input: +% - c3d: +% struct. Containing c3d data of a stimulus, given by importFiles function. +% +% - samiOptions.labels: +% cell array. Defines sami_standard names and order for labels. Needed to +% identify individuals by searching for substring unequal to sami_labels. +% +% output: +% - c3d: +% struct. same as input, but now with renamed labels in sami_format. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 10/2020 + +%% init vars +nLabelsUnsorted = numel(c3d.LabelsUnsorted); + +%% find individual persons and their 'identifiers' in each given labels ++++++++++++++++++ +[personIdentifier,personIndex,c3d.nPersons] = sami.c3d.getPersonIdentifier(c3d, samiOptions); + +%% rename labels by adding person# to sami-standard ++++++++++++++++++++++++++++++++++++++ +for i = 1:nLabelsUnsorted + oldLabel = c3d.LabelsUnsorted{i}; + % if oldLabel is asigned to a person: erasing old personIdentifier and adding p# + if ~isnan(personIndex(i)) + personString = personIdentifier{personIndex(i)}; + c3d.LabelsUnsorted{i} = ['p' num2str(personIndex(i)) erase(oldLabel,personString)]; + end +end + +end + diff --git a/sami_toolbox/+sami/+feat/+idv/limb_angles.m b/sami_toolbox/+sami/+feat/+idv/limb_angles.m new file mode 100644 index 0000000..c84614b --- /dev/null +++ b/sami_toolbox/+sami/+feat/+idv/limb_angles.m @@ -0,0 +1,61 @@ +function results = limb_angles(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% use this markers +calc_angles = {{'RELB', 'RSHO', 'RHIP'};... + {'LELB', 'LSHO', 'LHIP'};... + {'RSHO', 'RELB', 'RWRI'};... + {'LSHO', 'LELB', 'LWRI'};... + {'RSHO', 'RHIP', 'RKNE'};... + {'LSHO', 'LHIP', 'LKNE'};... + {'RHIP', 'RKNE', 'RANK'};... + {'LHIP', 'LKNE', 'LANK'}}; + +%% init vars +timesteps = size(data.marker,1); +nP = 2; % soften ?? +nAngles = size(calc_angles,1); + +%% get marker indices in dataset in a [angle*marker*person] matrix +mIdx = nan(nAngles,3,nP); +for p = 1:nP + for a = 1:nAngles + mIdx(a,:,p) = [find(strcmp(data.labels, ['p' num2str(p) calc_angles{a}{1}])), find(strcmp(data.labels, ['p' num2str(p) calc_angles{a}{2}])), find(strcmp(data.labels, ['p' num2str(p) calc_angles{a}{3}]))]; + end +end + +%% loop: time * angles * persons +angles_ts = nan(timesteps,nAngles,nP); +for p = 1:nP + for t = 1:timesteps + for a = 1:nAngles + % extract specific marker-information for each angle to calculate + % into [marker*xyz] matrix containing 3d information of each marker + m3d = squeeze(data.marker(t,mIdx(a,:,p),:)); + + % define vector in space and calculate angle between them + vec.B2A = m3d(1,:) - m3d(2,:); + vec.B2C = m3d(3,:) - m3d(2,:); + angles_ts(t,a,p) = acosd( dot ( vec.B2A, vec.B2C ) / ( norm( vec.B2A ) * norm( vec.B2C ) ) ); + end + end +end + +% average timeseries +angle_mean = squeeze(mean(angles_ts)); +% average across persons +limbAngle_avg = mean(angle_mean,2); + +%% Results +results.feat = limbAngle_avg; +results.name = 'limb_angles'; +results.color = [.2 .5 0]; +results.unit = 'deg'; + +end diff --git a/sami_toolbox/+sami/+feat/+idv/limb_contraction.m b/sami_toolbox/+sami/+feat/+idv/limb_contraction.m new file mode 100644 index 0000000..7b42cf3 --- /dev/null +++ b/sami_toolbox/+sami/+feat/+idv/limb_contraction.m @@ -0,0 +1,50 @@ +function results = limb_contraction(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% use this markers +limbs = {'RWRI', 'LWRI','RANK', 'LANK'}; + +%% init vars +head = {'HEAD'}; +timesteps = size(data.marker,1); +nP = 2; % soften? +nLimbs = numel(limbs); +dist_ts = nan(timesteps,nLimbs,nP); + +%% calculate distances for each person +for p = 1:nP + % get head index + idxH = strcmp(data.labels, ['p' num2str(p) head{1}]); + + % loop limbs + for limb = 1:nLimbs + idxL = strcmp(data.labels, ['p' num2str(p) limbs{limb}]); + + xyz.head = squeeze(data.marker(:,idxH,:)); + xyz.limb = squeeze(data.marker(:,idxL,:)); + + tmpData = [reshape(xyz.head,1,3*timesteps);reshape(xyz.limb,1,3*timesteps)]; + tmpData = reshape(tmpData,2*timesteps,3); + tmpData = squareform(pdist(tmpData,'euclidean')); + + dist_ts(:,limb,p) = diag(tmpData(1:2:2*timesteps,2:2:2*timesteps)); + end +end + +% calculating means +limbContraction_avg = mean(mean(dist_ts),3); + +%% Results +results.feat = limbContraction_avg; +results.name = 'limb_contraction'; +results.color = [.2 .5 0]; +results.unit = 'mm'; + +end + diff --git a/sami_toolbox/+sami/+feat/+idv/symmetry.m b/sami_toolbox/+sami/+feat/+idv/symmetry.m new file mode 100644 index 0000000..cc5407e --- /dev/null +++ b/sami_toolbox/+sami/+feat/+idv/symmetry.m @@ -0,0 +1,139 @@ +function results = symmetry(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init vars +timesteps = size(data.marker,1); +nP = 2; + +%% define markers +marker_shoulder = {'LSHO', 'RSHO'}; +marker_hip = {'LHIP', 'RHIP'}; + +%% loop persons +for p = 1:nP + %% find marker indices + % shoulder + hip + marker_idx.shoulder = [find(strcmp(data.labels, ['p' num2str(p) marker_shoulder{1}])), find(strcmp(data.labels, ['p' num2str(p) marker_shoulder{2}]))]; + marker_idx.hip = [find(strcmp(data.labels, ['p' num2str(p) marker_hip{1}])), find(strcmp(data.labels, ['p' num2str(p) marker_hip{2}]))]; + + % all left and right markers + marker_idx.left = []; + marker_idx.right = []; + for i = 1:numel(data.labels) + if strcmp(data.labels{1,i}(1:3), ['p' num2str(p) 'L']) + marker_idx.left = [marker_idx.left i]; + % find coresponding right marker (should be in corresponding order, but u never know...) + tmpLabel = [data.labels{1,i}(1:2) 'R' data.labels{1,i}(4:end)]; + marker_idx.right = [marker_idx.right find(strcmp(data.labels, tmpLabel))]; + end + end + nMarker = numel(marker_idx.right); + + %% calculate difference in height (z-axes) between body-sides, for each marker in each timestep + symmetry.height.diff(:,:,p) = abs( ( data.marker(:,marker_idx.left,3) - data.marker(:,marker_idx.right,3) ) ); + + %% calculate deviations (angle and distance) from symmetry-line + for t = 1:timesteps + % calculate midline [location and pointing_vector] for each person + [ml_loc, ml_vec] = calc_midline(data,t,marker_idx); + + % loop marker and get distance/angle with respect to midline + for m = 1:nMarker + % distance to symmetry_line + tmpL = squeeze(data.marker(t,marker_idx.left(m),1:2)); + tmpR = squeeze(data.marker(t,marker_idx.right(m),1:2)); + symmetry.ml_dist.left(t,m,p) = dist_to_symmetry_line(tmpL, ml_loc, ml_loc + ml_vec); + symmetry.ml_dist.right(t,m,p) = dist_to_symmetry_line(tmpR, ml_loc, ml_loc + ml_vec); + + % circular segment + tmpLRnorm = mean([norm(tmpL - ml_loc) norm(tmpR - ml_loc)]); + tmpLalpha = asin(symmetry.ml_dist.left(t,m,p) / norm(ml_loc - tmpL)); + tmpRalpha = asin(symmetry.ml_dist.right(t,m,p) / norm(ml_loc - tmpR)); + symmetry.ml_circSeg.left(t,m,p) = tmpLRnorm * tmpLalpha; + symmetry.ml_circSeg.right(t,m,p) = tmpLRnorm * tmpRalpha; + end %marker loop + end %time loop +end %person loop + +% average anatomical marker across persons +symmetry.height.avg_per_marker = mean(symmetry.height.diff); +symmetry.height.avg_diff = mean(symmetry.height.avg_per_marker, 3)'; + +symmetry.ml_dist.diff = abs(symmetry.ml_dist.right - symmetry.ml_dist.left); +symmetry.ml_dist.avg_per_marker = mean(symmetry.ml_dist.diff); +symmetry.ml_dist.avg_diff = mean(symmetry.ml_dist.avg_per_marker, 3)'; + +symmetry.ml_circSeg.diff = abs(symmetry.ml_circSeg.right - symmetry.ml_circSeg.left); +symmetry.ml_circSeg.avg_per_marker = mean(symmetry.ml_circSeg.diff); +symmetry.ml_circSeg.avg_diff = mean(symmetry.ml_circSeg.avg_per_marker, 3)'; + +symmetry.results = [symmetry.height.avg_diff; symmetry.ml_dist.avg_diff; symmetry.ml_circSeg.avg_diff]; + +%% Results +results.feat = symmetry.results; +results.name = 'symmetry'; +results.color = [.2 .5 0]; +results.unit = 'delta(mm)'; + +end + +%% ******** SUB FUNCTIONS ****************************************************************** +function [r_mid, v_f] = calc_midline(data,t,marker_idx) +% shoulder ********************************************************************************* +% get coordinates of each marker +r_ls = permute( data.marker(t,marker_idx.shoulder(1),1:2) ,[3 1 2]); +r_rs = permute( data.marker(t,marker_idx.shoulder(2),1:2) ,[3 1 2]); + +% calc pointing vector as normal_vector between shoulders +v_lr_shoulder = r_rs - r_ls; + +v_lr_k_shoulder = norm(v_lr_shoulder); % k ~ k times n_vector is lenght of shoulder +v_lr_n_shoulder = v_lr_shoulder ./ v_lr_k_shoulder; % n ~ norm (einheitsvektor) + +% orthogonal vector, pointing forward +v_f_n_shoulder = (v_lr_n_shoulder' * [cosd(-90) -sind(-90); sind(-90) cosd(-90)])'; +v_f_shoulder = v_f_n_shoulder * 3000; + +% middle between shoulders +r_mid_shoulder = r_ls + (v_lr_k_shoulder/2 * v_lr_n_shoulder); + +% hip ************************************************************************************* +% get coordinates of each marker +r_lh = permute( data.marker(t,marker_idx.hip(1),1:2), [3,1,2]); +r_rh = permute( data.marker(t,marker_idx.hip(2),1:2), [3,1,2]); + +% calc pointing vector as normal_vector between shoulders +v_lr_hip = r_rh - r_lh; + +v_lr_k_hip = norm( v_lr_hip); % k ~ k times n_vector is lenght of shoulder +v_lr_n_hip = v_lr_hip ./ v_lr_k_hip; % n ~ norm (einheitsvektor) + +% orthogonal vector, pointing forward +v_f_n_hip= (v_lr_n_hip' * [cosd(-90) -sind(-90); sind(-90) cosd(-90)])'; +v_f_hip = v_f_n_hip * 3000; + +% middle between shoulders +r_mid_hip = r_lh + (v_lr_k_hip/2 * v_lr_n_hip); + +% mean(shoulder+hip) LOCATION ************************************************************* +r_mid = mean([r_mid_shoulder,r_mid_hip],2); + +% mean(shoulder+hip) NORMAL_VECTOR ********************************************************* +v_f = mean([v_f_shoulder,v_f_hip],2); +end + +function d = dist_to_symmetry_line(pt, v1, v2) +pt(3) = 0; +v1(3) = 0; +v2(3) = 0; + +a = v1 - v2; +b = pt - v2; +d = norm(cross(a,b)) / norm(a); +end diff --git a/sami_toolbox/+sami/+feat/+idv/vel_acc.m b/sami_toolbox/+sami/+feat/+idv/vel_acc.m new file mode 100644 index 0000000..339c152 --- /dev/null +++ b/sami_toolbox/+sami/+feat/+idv/vel_acc.m @@ -0,0 +1,42 @@ +function results = vel_acc(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init vars +dt = 1/data.frameRate; % set time interval +nMarkerPerPerson = numel(data.labels)/2; + +%% calculations +% pathlength of each marker in 3d space +path_diff = diff(data.marker); +path = sqrt(sum(path_diff.*path_diff,3)); + +% velocity +vel = path ./ dt; +vel_avg_marker = mean(vel); + +% acceleration +acc = diff(vel) ./ dt; +acc_avg_marker = mean(acc); + +% average across individuals +vel_avg = mean(reshape(vel_avg_marker,nMarkerPerPerson,2),2); +acc_avg = mean(reshape(acc_avg_marker,nMarkerPerPerson,2),2); + +%% Results +results(1).feat = vel_avg; +results(1).name = 'velocity'; +results(1).color = [.8 .8 .0]; +results(1).unit = 'mm*s^-1'; + +results(2).feat = acc_avg; +results(2).name = 'acceleration'; +results(2).color = [.8 .8 .0]; +results(2).unit = 'mm*s^-2'; + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+feat/+idv/vertical_movement.m b/sami_toolbox/+sami/+feat/+idv/vertical_movement.m new file mode 100644 index 0000000..f94af72 --- /dev/null +++ b/sami_toolbox/+sami/+feat/+idv/vertical_movement.m @@ -0,0 +1,27 @@ +function results = vertical_movement(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init vars +nMarkerPerPerson = numel(data.labels)/2; + +% calculations +vPath = abs(diff(data.marker(:,:,3))); +vPath_sum = sum(vPath); + +% average across individuals +vPath_avg = mean(reshape(vPath_sum,nMarkerPerPerson,2),2); + +%% Results +results.feat = vPath_avg; +results.name = 'vertical_movement'; +results.color = [.8 .8 0]; +results.unit = 'mm'; + + +end diff --git a/sami_toolbox/+sami/+feat/+idv/volume.m b/sami_toolbox/+sami/+feat/+idv/volume.m new file mode 100644 index 0000000..b67f3ae --- /dev/null +++ b/sami_toolbox/+sami/+feat/+idv/volume.m @@ -0,0 +1,43 @@ +function results = volume(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init vars +timesteps = size(data.marker,1); +nP = 2; + +%% get range of marker-3d-values for each timestep, each person and in each dimension +marker_range = nan(timesteps,nP,3); +for p = 1:nP + for a = 1:3 + marker_range(:,p,a) = range( data.marker(:,contains(data.labels, ['p' num2str(p)]),a)' ) ./ 1000; % convert into [m] + end +end + + +%% calculate vol for each person seperately +vol_ts = prod(marker_range,3); +vol_MN = mean(vol_ts); +vol_STD = std(vol_ts); + +% average across persons +volume_avg = mean(vol_MN); +volume_std = mean(vol_STD); + +%% Results +results(1).feat = volume_avg; +results(1).name = 'volume_average'; +results(1).color = [.2 .5 0]; +results(1).unit = 'm^3'; + +results(2).feat = volume_std; +results(2).name = 'volume_std'; +results(2).color = [.2 .5 0]; +results(2).unit = 'm^3'; + +end diff --git a/sami_toolbox/+sami/+feat/+itx/NOT_hand_movements.m b/sami_toolbox/+sami/+feat/+itx/NOT_hand_movements.m new file mode 100644 index 0000000..5e7e055 --- /dev/null +++ b/sami_toolbox/+sami/+feat/+itx/NOT_hand_movements.m @@ -0,0 +1,229 @@ +function results = hand_movements(data) +% +% EXAMPLE HOW TO *SKIP* FUNCTIONS +% just rename it and write "NOT" somewhere into the filename +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init vars +timesteps = size(data.marker,1); +dt = 1/data.frameRate; + +p1_lh = nan(timesteps,3); +p1_rh = nan(timesteps,3); +p2_lh = nan(timesteps,3); +p2_rh = nan(timesteps,3); + +%% use this markers +marker_lh = {'LWRA', 'LWRB', 'LFIN'}; +marker_rh = {'RWRA', 'RWRB', 'RFIN'}; +marker_idx_p1_lh = [find(strcmp(data.labels, marker_lh{1})),find(strcmp(data.labels.p1, marker_lh{2})),find(strcmp(data.labels.p1, marker_lh{3}))]; +marker_idx_p1_rh = [find(strcmp(data.labels, marker_rh{1})),find(strcmp(data.labels.p1, marker_rh{2})),find(strcmp(data.labels.p1, marker_rh{3}))]; +marker_idx_p2_lh = [find(strcmp(data.labels, marker_lh{1})),find(strcmp(data.labels.p2, marker_lh{2})),find(strcmp(data.labels.p2, marker_lh{3}))]; +marker_idx_p2_rh = [find(strcmp(data.labels, marker_rh{1})),find(strcmp(data.labels.p2, marker_rh{2})),find(strcmp(data.labels.p2, marker_rh{3}))]; + +%% TIME loop +for t = 1:timesteps + % calc midpoint of centroid for each Person and each Hand + m_p1_lA = permute( data.marker.p1(t,marker_idx_p1_lh(1),:) ,[3 1 2]); + m_p1_lB = permute( data.marker.p1(t,marker_idx_p1_lh(2),:) ,[3 1 2]); + m_p1_lF = permute( data.marker.p1(t,marker_idx_p1_lh(3),:) ,[3 1 2]); + p1_lh(t,1:3) = (m_p1_lA + m_p1_lB + m_p1_lF) / 3; + + m_p1_rA = permute( data.marker.p1(t,marker_idx_p1_rh(1),:) ,[3 1 2]); + m_p1_rB = permute( data.marker.p1(t,marker_idx_p1_rh(2),:) ,[3 1 2]); + m_p1_rF = permute( data.marker.p1(t,marker_idx_p1_rh(3),:) ,[3 1 2]); + p1_rh(t,1:3) = (m_p1_rA + m_p1_rB + m_p1_rF) / 3; + + m_p2_lA = permute( data.marker.p2(t,marker_idx_p2_lh(1),:) ,[3 1 2]); + m_p2_lB = permute( data.marker.p2(t,marker_idx_p2_lh(2),:) ,[3 1 2]); + m_p2_lF = permute( data.marker.p2(t,marker_idx_p2_lh(3),:) ,[3 1 2]); + p2_lh(t,1:3) = (m_p2_lA + m_p2_lB + m_p2_lF) / 3; + + m_p2_rA = permute( data.marker.p2(t,marker_idx_p2_rh(1),:) ,[3 1 2]); + m_p2_rB = permute( data.marker.p2(t,marker_idx_p2_rh(2),:) ,[3 1 2]); + m_p2_rF = permute( data.marker.p2(t,marker_idx_p2_rh(3),:) ,[3 1 2]); + p2_rh(t,1:3) = (m_p2_rA + m_p2_rB + m_p2_rF) / 3; + + %% calc parameters + if t > 1 + % way from last timestep to this timestep + handway.timeseries_p1_lh(t-1,1) = pdist([p1_lh(t,:); p1_lh(t-1,:)],'euclidean'); + handway.timeseries_p1_rh(t-1,1) = pdist([p1_rh(t,:); p1_rh(t-1,:)],'euclidean'); + handway.timeseries_p2_lh(t-1,1) = pdist([p2_lh(t,:); p2_lh(t-1,:)],'euclidean'); + handway.timeseries_p2_rh(t-1,1) = pdist([p2_rh(t,:); p2_rh(t-1,:)],'euclidean'); + end +end % end TIMESTEPS loop + +%% calc stuff in LAST TIMESTEP +% ****** save coordinates ****** +results.coordinates.coord_p1_lh = p1_lh; +results.coordinates.coord_p1_rh = p1_rh; +results.coordinates.coord_p2_lh = p2_lh; +results.coordinates.coord_p2_rh = p2_rh; + +% ****** handway ****** +results.way = handway; +results.way.p1_lh = sum(results.way.timeseries_p1_lh); +results.way.p1_rh = sum(results.way.timeseries_p1_rh); +results.way.p2_lh = sum(results.way.timeseries_p2_lh); +results.way.p2_rh = sum(results.way.timeseries_p2_rh); + +results.way.mean = mean([results.way.p1_lh; results.way.p1_rh; results.way.p2_lh; results.way.p2_rh]); +results.way.mean_p1 = mean([results.way.p1_lh; results.way.p1_rh]); +results.way.mean_p2 = mean([results.way.p2_lh; results.way.p2_rh]); + +% ****** velocity ****** +results.velo.timeseries_p1_lh = (results.way.timeseries_p1_lh) / dt; +results.velo.timeseries_p1_rh = (results.way.timeseries_p1_rh) / dt; +results.velo.timeseries_p2_lh = (results.way.timeseries_p2_lh) / dt; +results.velo.timeseries_p2_rh = (results.way.timeseries_p2_rh) / dt; + +results.velo.p1_lh = mean(results.velo.timeseries_p1_lh); +results.velo.p1_rh = mean(results.velo.timeseries_p1_rh); +results.velo.p2_lh = mean(results.velo.timeseries_p2_lh); +results.velo.p2_rh = mean(results.velo.timeseries_p2_rh); + +results.velo.mean = mean([results.velo.p1_lh; results.velo.p1_rh; results.velo.p2_lh; results.velo.p2_rh]); +results.velo.mean_p1 = mean([results.velo.p1_lh; results.velo.p1_rh]); +results.velo.mean_p2 = mean([results.velo.p2_lh; results.velo.p2_rh]); + +% ****** acceleration ****** +results.acc.timeseries_p1_lh = diff(results.velo.timeseries_p1_lh) / dt; +results.acc.timeseries_p1_rh = diff(results.velo.timeseries_p1_rh) / dt; +results.acc.timeseries_p2_lh = diff(results.velo.timeseries_p2_lh) / dt; +results.acc.timeseries_p2_rh = diff(results.velo.timeseries_p2_rh) / dt; + +results.acc.p1_lh = mean(results.acc.timeseries_p1_lh); +results.acc.p1_rh = mean(results.acc.timeseries_p1_rh); +results.acc.p2_lh = mean(results.acc.timeseries_p2_lh); +results.acc.p2_rh = mean(results.acc.timeseries_p2_rh); + +results.acc.mean = mean([results.acc.p1_lh; results.acc.p1_rh; results.acc.p2_lh; results.acc.p2_rh]); +results.acc.mean_p1 = mean([results.acc.p1_lh; results.acc.p1_rh]); +results.acc.mean_p2 = mean([results.acc.p2_lh; results.acc.p2_rh]); + + +% ****** handvol of ellipsoid ****** +[ellipsoid_p1l] = MinVolEllipse(p1_lh'./1000,0.01); [~,Q,~] = svd(ellipsoid_p1l); +handvol.p1_lh = 4/3 * pi * ( 1/sqrt(Q(1,1)) ) * ( 1/sqrt(Q(2,2)) ) * ( 1/sqrt(Q(3,3)) ); +[ellipsoid_p1r] = MinVolEllipse(p1_rh'./1000,0.01); [~,Q,~] = svd(ellipsoid_p1r); +handvol.p1_rh = 4/3 * pi * ( 1/sqrt(Q(1,1)) ) * ( 1/sqrt(Q(2,2)) ) * ( 1/sqrt(Q(3,3)) ); +[ellipsoid_p2l] = MinVolEllipse(p2_lh'./1000,0.01); [~,Q,~] = svd(ellipsoid_p2l); +handvol.p2_lh = 4/3 * pi * ( 1/sqrt(Q(1,1)) ) * ( 1/sqrt(Q(2,2)) ) * ( 1/sqrt(Q(3,3)) ); +[ellipsoid_p2r] = MinVolEllipse(p2_rh'./1000,0.01); [~,Q,~] = svd(ellipsoid_p2r); +handvol.p2_rh = 4/3 * pi * ( 1/sqrt(Q(1,1)) ) * ( 1/sqrt(Q(2,2)) ) * ( 1/sqrt(Q(3,3)) ); + +results.vol.mean = mean([handvol.p1_lh; handvol.p1_rh; handvol.p2_lh; handvol.p2_rh]); +results.vol.mean_p1 = mean([handvol.p1_lh; handvol.p1_rh]); +results.vol.mean_p2 = mean([handvol.p2_lh; handvol.p2_rh]); + +%% Results +results.feat = 123; +results.name = 'avg_distance'; +results.color = [0 .8 .2]; +results.unit = 'idk'; + +end + + + +function [A , c] = MinVolEllipse(P, tolerance) +% [A , c] = MinVolEllipse(P, tolerance) +% Finds the minimum volume enclsing ellipsoid (MVEE) of a set of data +% points stored in matrix P. The following optimization problem is solved: +% +% minimize log(det(A)) +% subject to (P_i - c)' * A * (P_i - c) <= 1 +% +% in variables A and c, where P_i is the i-th column of the matrix P. +% The solver is based on Khachiyan Algorithm, and the final solution +% is different from the optimal value by the pre-spesified amount of 'tolerance'. +% +% inputs: +%--------- +% P : (d x N) dimnesional matrix containing N points in R^d. +% tolerance : error in the solution with respect to the optimal value. +% +% outputs: +%--------- +% A : (d x d) matrix of the ellipse equation in the 'center form': +% (x-c)' * A * (x-c) = 1 +% c : 'd' dimensional vector as the center of the ellipse. +% +% example: +% -------- +% P = rand(5,100); +% [A, c] = MinVolEllipse(P, .01) +% +% To reduce the computation time, work with the boundary points only: +% +% K = convhulln(P'); +% K = unique(K(:)); +% Q = P(:,K); +% [A, c] = MinVolEllipse(Q, .01) +% +% +% Nima Moshtagh (nima@seas.upenn.edu) +% University of Pennsylvania +% +% December 2005 +% UPDATE: Jan 2009 + + + +%%%%%%%%%%%%%%%%%%%%% Solving the Dual problem%%%%%%%%%%%%%%%%%%%%%%%%%%%5 +% --------------------------------- +% data points +% ----------------------------------- +[d, N] = size(P); + +Q = zeros(d+1,N); +Q(1:d,:) = P(1:d,1:N); +Q(d+1,:) = ones(1,N); + + +% initializations +% ----------------------------------- +count = 1; +err = 1; +u = (1/N) * ones(N,1); % 1st iteration + + +% Khachiyan Algorithm +% ----------------------------------- +while err > tolerance + X = Q * diag(u) * Q'; % X = \sum_i ( u_i * q_i * q_i') is a (d+1)x(d+1) matrix +% M = diag(Q' * inv(X) * Q); % M the diagonal vector of an NxN matrix + M = diag(Q' * (X\Q)); % M the diagonal vector of an NxN matrix + [maximum, j] = max(M); + step_size = (maximum - d -1)/((d+1)*(maximum-1)); + new_u = (1 - step_size)*u ; + new_u(j) = new_u(j) + step_size; + count = count + 1; + err = norm(new_u - u); + u = new_u; +end + + + +%%%%%%%%%%%%%%%%%%% Computing the Ellipse parameters%%%%%%%%%%%%%%%%%%%%%% +% Finds the ellipse equation in the 'center form': +% (x-c)' * A * (x-c) = 1 +% It computes a dxd matrix 'A' and a d dimensional vector 'c' as the center +% of the ellipse. + +U = diag(u); + +% the A matrix for the ellipse +% -------------------------------------------- +% A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' ); +A = (1/d) / (P * U * P' - (P * u)*(P*u)' ); + + +% center of the ellipse +% -------------------------------------------- +c = P * u; + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+feat/+itx/distance_correlation.m b/sami_toolbox/+sami/+feat/+itx/distance_correlation.m new file mode 100644 index 0000000..a2e0e9e --- /dev/null +++ b/sami_toolbox/+sami/+feat/+itx/distance_correlation.m @@ -0,0 +1,148 @@ +function results = distance_correlation(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% options +mw_size = .1; % size of moving window [in seconds] + +%% init vars +nP = 2; +nMarker = size(data.labels,2); +timesteps = size(data.marker,1); +IPD = nan(timesteps,1); +dt = 1/data.frameRate; % set time interval +m_w = mw_size * data.frameRate; % moving average window + +%% use this markers to calculate virtual marker at assumed "center" of each person +marker = {'LHIP', 'RHIP'}; +marker_idx_p1 = [find(strcmp(data.labels, ['p1' marker{1}])), find(strcmp(data.labels, ['p1' marker{2}]))]; +marker_idx_p2 = [find(strcmp(data.labels, ['p2' marker{1}])), find(strcmp(data.labels, ['p2' marker{2}]))]; + +% for limb contraction +head = {'HEAD'}; +limbs = {'RWRI', 'LWRI','RANK', 'LANK'}; +nLimbs = numel(limbs); + +%% calculate timeseries of interpersonal distance +for t = 1:timesteps + % calc virtual center for each person + r_p1_m1 = data.marker(t,marker_idx_p1(1),:); + r_p1_m2 = data.marker(t,marker_idx_p1(2),:); + mdp_p1 = squeeze((r_p1_m1 + r_p1_m2) / 2); + + r_p2_m1 = data.marker(t,marker_idx_p2(1),:); + r_p2_m2 = data.marker(t,marker_idx_p2(2),:); + mdp_p2 = squeeze((r_p2_m1 + r_p2_m2) / 2); + + % calc interpersonal distance between virtual centers + IPD(t) = pdist([mdp_p1, mdp_p2]','euclidean'); +end % end TIME loop + +%% correlation: IPD vs. vel/acc +IPD_vel = IPD(2:end); +IPD_acc = IPD(3:end); +corr_vel = nan(1,numel(data.labels)); +corr_acc = nan(1,numel(data.labels)); + +% calculate velocity and acceleration of each marker +path_diff = diff(data.marker); +path = sqrt(sum(path_diff.*path_diff,3)); + +vel_raw = path ./ dt; +vel = movmean(vel_raw, m_w); +acc_raw = diff(vel)./ dt; +acc = movmean(acc_raw, m_w); + +% correlations +for m = 1:nMarker + corr_vel(m) = corr(IPD_vel, vel(:,m)); + corr_acc(m) = corr(IPD_acc, acc(:,m)); +end + +% post-processing +corr_vel = abs(reshape(corr_vel, [nMarker/2, 2])); % abs +mean_corr_vel = tanh(mean(atanh(corr_vel),2)); % fisher z + +corr_acc = abs(reshape(corr_acc, [nMarker/2, 2])); % abs +mean_corr_acc = tanh(mean(atanh(corr_acc),2)); % fisher z + +%% correlation: IPD vs. limb contraction +corr_lc = nan(4,nP); + +for p = 1:nP + % get head index + idxH = strcmp(data.labels, ['p' num2str(p) head{1}]); + + % loop limbs + for limb = 1:nLimbs + % get index of limb and 3d_data + idxL = strcmp(data.labels, ['p' num2str(p) limbs{limb}]); + + xyz.head = squeeze(data.marker(:,idxH,:)); + xyz.limb = squeeze(data.marker(:,idxL,:)); + + % timeseries of distance between head and limb + tmpData = [reshape(xyz.head,1,3*timesteps);reshape(xyz.limb,1,3*timesteps)]; + tmpData = reshape(tmpData,2*timesteps,3); + tmpData = squareform(pdist(tmpData,'euclidean')); + + lc_ts = diag(tmpData(1:2:2*timesteps,2:2:2*timesteps)); + + % correlation betwenn IPD vs LC + corr_lc(limb,p) = corr(IPD, lc_ts); + end +end + +% post-processing +corr_lc = abs(corr_lc); % abs +mean_corr_lc = tanh(mean(atanh(corr_lc),2)); % fisher z + +%% correlation: IPD vs. volume +marker_range = nan(timesteps,nP,3); + +% calculate vol +for p = 1:nP + for a = 1:3 + marker_range(:,p,a) = range( data.marker(:,contains(data.labels, ['p' num2str(p)]),a)' ); + end +end +vol_ts = prod(marker_range,3); + +% correlation +corr_vol = corr(vol_ts, IPD); + +% post-processing +corr_vol = abs(corr_vol); % abs +mean_corr_vol = tanh(mean(atanh(corr_vol))); % fisher z + +%% Results +results(1).feat = mean_corr_vel; +results(1).name = 'corr_dist_velocity'; +results(1).color = [0 .8 .6]; +results(1).unit = '|Pearson r|'; +results(1).fisherZ4paramTesting = true; + +results(2).feat = mean_corr_acc; +results(2).name = 'corr_dist_acceleration'; +results(2).color = [0 .8 .6]; +results(2).unit = '|Pearson r|'; +results(2).fisherZ4paramTesting = true; + +results(3).feat = mean_corr_lc; +results(3).name = 'corr_dist_limb_contraction'; +results(3).color = [0 .8 .6]; +results(3).unit = '|Pearson r|'; +results(3).fisherZ4paramTesting = true; + +results(4).feat = mean_corr_vol; +results(4).name = 'corr_dist_volume'; +results(4).color = [0 .8 .6]; +results(4).unit = '|Pearson r|'; +results(4).fisherZ4paramTesting = true; + +end diff --git a/sami_toolbox/+sami/+feat/+itx/interpersonal_distance.m b/sami_toolbox/+sami/+feat/+itx/interpersonal_distance.m new file mode 100644 index 0000000..97e1896 --- /dev/null +++ b/sami_toolbox/+sami/+feat/+itx/interpersonal_distance.m @@ -0,0 +1,49 @@ +function results = interpersonal_distance(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init vars +timesteps = size(data.marker,1); +IPD_ts = nan(timesteps,1); + +%% use this markers to calculate virtual marker at assumed "center" of each person +marker = {'LHIP', 'RHIP'}; +marker_idx_p1 = [find(strcmp(data.labels, ['p1' marker{1}])), find(strcmp(data.labels, ['p1' marker{2}]))]; +marker_idx_p2 = [find(strcmp(data.labels, ['p2' marker{1}])), find(strcmp(data.labels, ['p2' marker{2}]))]; + +%% loop TIME +for t = 1:timesteps + % calc virtual center for each person + r_p1_m1 = data.marker(t,marker_idx_p1(1),:); + r_p1_m2 = data.marker(t,marker_idx_p1(2),:); + mdp_p1 = squeeze((r_p1_m1 + r_p1_m2) / 2); + + r_p2_m1 = data.marker(t,marker_idx_p2(1),:); + r_p2_m2 = data.marker(t,marker_idx_p2(2),:); + mdp_p2 = squeeze((r_p2_m1 + r_p2_m2) / 2); + + % calc interpersonal distance between virtual centers + IPD_ts(t) = pdist([mdp_p1, mdp_p2]','euclidean'); +end % end TIME loop + +%% calc mean interpersonal distance +IPD_mean = mean(IPD_ts); +IPD_std = std(IPD_ts); + +%% Results +results(1).feat = IPD_mean; +results(1).name = 'IP_distance_avg'; +results(1).color = [0 .8 .6]; +results(1).unit = 'mm'; + +results(2).feat = IPD_std; +results(2).name = 'IP_distance_std'; +results(2).color = [0 .8 .6]; +results(2).unit = 'mm'; + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+feat/+itx/interpersonal_orientation.m b/sami_toolbox/+sami/+feat/+itx/interpersonal_orientation.m new file mode 100644 index 0000000..98b317d --- /dev/null +++ b/sami_toolbox/+sami/+feat/+itx/interpersonal_orientation.m @@ -0,0 +1,107 @@ +function results = interpersonal_orientation(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init vars +timesteps = size(data.marker,1); +IPO_p1_p2 = zeros(timesteps,1); % if p1 toward p2 +IPO_p2_p1 = zeros(timesteps,1); % if p2 toward p1 + +%% use this markers +marker = {'LSHO', 'RSHO'}; +marker_idx_p1 = [find(strcmp(data.labels, ['p1' marker{1}])), find(strcmp(data.labels, ['p1' marker{2}]))]; +marker_idx_p2 = [find(strcmp(data.labels, ['p2' marker{1}])), find(strcmp(data.labels, ['p2' marker{2}]))]; + +%% TIME loop +for t = 1:timesteps + % get coordinates of each marker + r_p1_ls = permute( data.marker(t,marker_idx_p1(1),1:2) ,[3 1 2]); + r_p1_rs = permute( data.marker(t,marker_idx_p1(2),1:2) ,[3 1 2]); + r_p2_ls = permute( data.marker(t,marker_idx_p2(1),1:2) ,[3 1 2]); + r_p2_rs = permute( data.marker(t,marker_idx_p2(2),1:2) ,[3 1 2]); + + r_p1_mdp = [mean([r_p1_ls(1) r_p1_rs(1)]); mean([r_p1_ls(2) r_p1_rs(2)])]; + r_p2_mdp = [mean([r_p2_ls(1) r_p2_rs(1)]); mean([r_p2_ls(2) r_p2_rs(2)])]; + + % get max dist between all shoulder markers + dist(1) = pdist([r_p1_ls, r_p2_ls]','euclidean'); + dist(2) = pdist([r_p1_ls, r_p2_rs]','euclidean'); + dist(3) = pdist([r_p1_rs, r_p2_ls]','euclidean'); + dist(4) = pdist([r_p1_rs, r_p2_rs]','euclidean'); + dist = max(dist); + + % calc pointing vector as normal_vector between shoulders + v_p1_lr = r_p1_rs - r_p1_ls; + v_p2_lr = r_p2_rs - r_p2_ls; + + v_p1_lr_k = norm( v_p1_lr); % k ~ k times n_vector is lenght of shoulder + v_p1_lr_n = v_p1_lr ./ v_p1_lr_k; % n ~ norm (einheitsvektor) + v_p2_lr_k = norm( v_p2_lr); + v_p2_lr_n = v_p2_lr ./ v_p2_lr_k; + + % orthogonal vector, pointing forward + v_p1_f_n = (v_p1_lr_n' * [cosd(-90) -sind(-90); sind(-90) cosd(-90)])'; + v_p2_f_n = (v_p2_lr_n' * [cosd(-90) -sind(-90); sind(-90) cosd(-90)])'; + v_p1_f = v_p1_f_n * dist * 1.1; + v_p2_f = v_p2_f_n * dist * 1.1; + + % check if persons are oriented towards each other + d_p1_ls_p2 = point_to_line(r_p2_mdp, r_p1_ls, r_p1_ls+v_p1_f); + d_p1_rs_p2 = point_to_line(r_p2_mdp, r_p1_rs, r_p1_rs+v_p1_f); + d_p2_ls_p1 = point_to_line(r_p1_mdp, r_p2_ls, r_p2_ls+v_p2_f); + d_p2_rs_p1 = point_to_line(r_p1_mdp, r_p2_rs, r_p2_rs+v_p2_f); + + % if p1 toward p2 + if (d_p1_ls_p2 <= v_p1_lr_k) && (d_p1_rs_p2 <= v_p1_lr_k) + IPO_p1_p2(t) = 1; + end + + % if p2 toward p1 + if (d_p2_ls_p1 <= v_p2_lr_k) && (d_p2_rs_p1 <= v_p2_lr_k) + IPO_p2_p1(t) = 1; + end +end % end TIMESTEPS loop + +% average across time +IPO_pp(1) = nanmean(IPO_p1_p2) * 100; +IPO_pp(2) = nanmean(IPO_p2_p1) * 100; + +% average IPO across persons +IPO_mean = mean(IPO_pp); + +% balance of IPO between persons +IPO_balance = 1 - abs(diff(IPO_pp)) / sum(IPO_pp); + +% catch case when both persons never look at each other +if isnan(IPO_balance) + IPO_balance = 0; +end + +%% Results +results(1).feat = IPO_mean; +results(1).name = 'IP_orientation_average'; +results(1).color = [0 .8 .6]; +results(1).unit = '%'; + +results(2).feat = IPO_balance; +results(2).name = 'IP_orientation_balance'; +results(2).color = [0 .8 .6]; +results(2).unit = 'AU'; + +end + + +function d = point_to_line(pt, v1, v2) +pt(3) = 0; +v1(3) = 0; +v2(3) = 0; + +a = v1 - v2; +b = pt - v2; +d = norm(cross(a,b)) / norm(a); +end diff --git a/sami_toolbox/+sami/+feat/+itx/motion_energy.m b/sami_toolbox/+sami/+feat/+itx/motion_energy.m new file mode 100644 index 0000000..0f0d14b --- /dev/null +++ b/sami_toolbox/+sami/+feat/+itx/motion_energy.m @@ -0,0 +1,35 @@ +function results = motion_energy(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% calculations +% indices of person specific markers +marker_idx(1,:) = find(contains(data.labels, 'p1')); +marker_idx(2,:)= find(contains(data.labels, 'p2')); + +% difference between frames for each dimension +diff_p1 = diff(data.marker(:,marker_idx(1,:),:)); +diff_p2 = diff(data.marker(:,marker_idx(2,:),:)); + +% length of 3d-interframe-difference (pythagoras) +disp_p1 = sqrt(sum(diff_p1.*diff_p1,3)); +disp_p2 = sqrt(sum(diff_p2.*diff_p2,3)); + +% balance of "motion_energy" between persons +MoEn(1) = mean(mean(disp_p1)); +MoEn(2) = mean(mean(disp_p2)); +MoEn_balance = 1 - abs(diff(MoEn)) / sum(MoEn); + +%% results +results.feat = MoEn_balance; +results.name = 'motion_energy_balance'; +results.color = [0 .8 .6]; +results.unit = 'AU'; + + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+feat/+itx/personalzone_time.m b/sami_toolbox/+sami/+feat/+itx/personalzone_time.m new file mode 100644 index 0000000..4ab587b --- /dev/null +++ b/sami_toolbox/+sami/+feat/+itx/personalzone_time.m @@ -0,0 +1,337 @@ +function results = personalzone_time(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init vars +timesteps = size(data.marker,1); +np = 2; + +%% use this markers +marker_idx.p1 = find(contains(data.labels, 'p1')); +marker_idx.p2 = find(contains(data.labels, 'p2')); +LabelsArms = {'LSHO', 'LELB', 'LWRI', 'RSHO', 'RELB', 'RWRI'}; +marker_idx.arms(1,:) = [find(strcmp(data.labels, ['p1' LabelsArms{1}])), find(strcmp(data.labels, ['p1' LabelsArms{2}])), find(strcmp(data.labels, ['p1' LabelsArms{3}])), find(strcmp(data.labels, ['p1' LabelsArms{4}])), find(strcmp(data.labels, ['p1' LabelsArms{5}])), find(strcmp(data.labels, ['p1' LabelsArms{6}]))]; +marker_idx.arms(2,:) = [find(strcmp(data.labels, ['p2' LabelsArms{1}])), find(strcmp(data.labels, ['p2' LabelsArms{2}])), find(strcmp(data.labels, ['p2' LabelsArms{3}])), find(strcmp(data.labels, ['p2' LabelsArms{4}])), find(strcmp(data.labels, ['p2' LabelsArms{5}])), find(strcmp(data.labels, ['p2' LabelsArms{6}]))]; +LabelsCore = {'LSHO', 'RSHO', 'LHIP', 'RHIP'}; +marker_idx.core(1,:) = [find(strcmp(data.labels, ['p1' LabelsCore{1}])), find(strcmp(data.labels, ['p1' LabelsCore{2}])), find(strcmp(data.labels, ['p1' LabelsCore{3}])), find(strcmp(data.labels, ['p1' LabelsCore{4}]))]; +marker_idx.core(2,:) = [find(strcmp(data.labels, ['p2' LabelsCore{1}])), find(strcmp(data.labels, ['p2' LabelsCore{2}])), find(strcmp(data.labels, ['p2' LabelsCore{3}])), find(strcmp(data.labels, ['p2' LabelsCore{4}]))]; + +%% get length arm +for p = 1:np + % all coords for each joint + tmp.lS = squeeze(data.marker(:,marker_idx.arms(p,1),:)); + tmp.lE = squeeze(data.marker(:,marker_idx.arms(p,2),:)); + tmp.lW = squeeze(data.marker(:,marker_idx.arms(p,3),:)); + tmp.rS = squeeze(data.marker(:,marker_idx.arms(p,4),:)); + tmp.rE = squeeze(data.marker(:,marker_idx.arms(p,5),:)); + tmp.rW = squeeze(data.marker(:,marker_idx.arms(p,6),:)); + + % left Shoulder_Elbow + tmp.lSE = [reshape(tmp.lS,1,3*timesteps);reshape(tmp.lE,1,3*timesteps)]; + tmp.lSE = reshape(tmp.lSE,2*timesteps,3); + tmp.lSE = squareform(pdist(tmp.lSE,'euclidean')); + % left Elbow_Wrist + tmp.lEW = [reshape(tmp.lE,1,3*timesteps);reshape(tmp.lW,1,3*timesteps)]; + tmp.lEW = reshape(tmp.lEW,2*timesteps,3); + tmp.lEW = squareform(pdist(tmp.lEW,'euclidean')); + % right Shoulder_Elbow + tmp.rSE = [reshape(tmp.rS,1,3*timesteps);reshape(tmp.rE,1,3*timesteps)]; + tmp.rSE = reshape(tmp.rSE,2*timesteps,3); + tmp.rSE = squareform(pdist(tmp.rSE,'euclidean')); + % right Elbow_Wrist + tmp.rEW = [reshape(tmp.rE,1,3*timesteps);reshape(tmp.rW,1,3*timesteps)]; + tmp.rEW = reshape(tmp.rEW,2*timesteps,3); + tmp.rEW = squareform(pdist(tmp.rEW,'euclidean')); + + % save timeseries of distances + tmp.lSE_ts(:,p) = diag(tmp.lSE(1:2:2*timesteps,2:2:2*timesteps)); + tmp.lEW_ts(:,p) = diag(tmp.lEW(1:2:2*timesteps,2:2:2*timesteps)); + tmp.rSE_ts(:,p) = diag(tmp.rSE(1:2:2*timesteps,2:2:2*timesteps)); + tmp.rEW_ts(:,p) = diag(tmp.rEW(1:2:2*timesteps,2:2:2*timesteps)); +end + +% calculating armLengths +tmp.lArm = mean(tmp.lSE_ts) + mean(tmp.lEW_ts); +tmp.rArm = mean(tmp.rSE_ts) + mean(tmp.rEW_ts); +armLength = mean(mean([tmp.lArm; tmp.rArm])); + +clear tmp; + +%% loop timesteps and check for "in-zone" +personalDist.p1in2 = zeros(timesteps,1); +personalDist.p2in1 = zeros(timesteps,1); + +for t = 1:timesteps + % calculate polygon for shoulder and hips + % get coordinates of each marker + r_p1_ls = permute( data.marker(t,marker_idx.core(1,1),:) ,[1 3 2]); + r_p1_rs = permute( data.marker(t,marker_idx.core(1,2),:) ,[1 3 2]); + r_p1_lh = permute( data.marker(t,marker_idx.core(1,3),:) ,[1 3 2]); + r_p1_rh = permute( data.marker(t,marker_idx.core(1,4),:) ,[1 3 2]); + r_p2_ls = permute( data.marker(t,marker_idx.core(2,1),:) ,[1 3 2]); + r_p2_rs = permute( data.marker(t,marker_idx.core(2,2),:) ,[1 3 2]); + r_p2_lh = permute( data.marker(t,marker_idx.core(2,3),:) ,[1 3 2]); + r_p2_rh = permute( data.marker(t,marker_idx.core(2,4),:) ,[1 3 2]); + + % get coordinates of orthogonal vec, pointing forward and backward + [p1f_ls, p1f_rs, p1b_ls, p1b_rs] = calc_point_front_back(r_p1_ls, r_p1_rs, armLength); + [p1f_lh, p1f_rh, p1b_lh, p1b_rh] = calc_point_front_back(r_p1_lh, r_p1_rh, armLength); + [p2f_ls, p2f_rs, p2b_ls, p2b_rs] = calc_point_front_back(r_p2_ls, r_p2_rs, armLength); + [p2f_lh, p2f_rh, p2b_lh, p2b_rh] = calc_point_front_back(r_p2_lh, r_p2_rh, armLength); + + % define convex hull + c_h_p1 = [p1f_ls; p1f_rs; p1b_ls; p1b_rs; p1f_lh; p1f_rh; p1b_lh; p1b_rh]; + c_h_p2 = [p2f_ls; p2f_rs; p2b_ls; p2b_rs; p2f_lh; p2f_rh; p2b_lh; p2b_rh]; + + %% check distance between each marker + % in p1 in p2 ? + for m = 1:numel(marker_idx.p1) + % check this marker + markerToCheckP1 = permute(data.marker(t, marker_idx.p1(m), :), [3,1,2]); + % inhull + inHull_p1in2_AZ = inhull(markerToCheckP1', c_h_p2); + % sphere + coreP2 = squeeze(data.marker(t, marker_idx.core(2,:), :)); + tmp_dist = pdist([markerToCheckP1'; coreP2],'euclidean'); + tmp_dist = tmp_dist(1:numel(LabelsCore)); + inSphere_p1in2 = tmp_dist <= armLength; + % check + if inHull_p1in2_AZ || any(inSphere_p1in2) + personalDist.p1in2(t) = 1; + break + end + end + % is p2 in p1 ? + for m = 1:numel(marker_idx.p2) + % check this marker + markerToCheckP2 = permute(data.marker(t, marker_idx.p2(m),:), [3,1,2]); + % inhull + inHull_p2in1_AZ = inhull(markerToCheckP2', c_h_p1); + % sphere + coreP1 = squeeze(data.marker(t, marker_idx.core(1,:), :)); + tmp_dist = pdist([markerToCheckP2'; coreP1],'euclidean'); + tmp_dist = tmp_dist(1:numel(LabelsCore)); + inSphere_p2in1 = tmp_dist <= armLength; + % check + if inHull_p2in1_AZ || any(inSphere_p2in1) + personalDist.p2in1(t) = 1; + break + end + end +end + + +%% results in % +p1in2_perc = sum(personalDist.p1in2) / timesteps * 100; +p2in1_perc = sum(personalDist.p2in1) / timesteps * 100; +result_PZ(1,1) = (p1in2_perc + p2in1_perc) / 2; + +%% results +results.feat = result_PZ; +results.name = 'personalspace'; +results.color = [0 .8 .6]; +results.unit = '%'; + +end + +%% +++++++++++ sub functions ++++++++++++++++++++++++++++++++++++++++++++++ +function [l_f, r_f, l_b, r_b] = calc_point_front_back(l_p, r_p, armLength) + % calc pointing vector as normal_vec + v_p_lr = r_p(1:2) - l_p(1:2); + v_p_lr_n = v_p_lr ./ norm( v_p_lr); % n ~ norm (einheitsvektor) + + % orthogonal vector, pointing forward and backward + v_p_f_n = v_p_lr_n * [cosd(-90) -sind(-90); sind(-90) cosd(-90)]; + v_p_f_s = v_p_f_n * armLength; + + % create corners by adding vector to shoulder/hips + l_f = l_p + [v_p_f_s, 0]; + r_f = r_p + [v_p_f_s, 0]; + l_b = l_p - [v_p_f_s, 0]; + r_b = r_p - [v_p_f_s, 0]; +end + +function in = inhull(testpts,xyz,tess,tol) +% inhull: tests if a set of points are inside a convex hull +% usage: in = inhull(testpts,xyz) +% usage: in = inhull(testpts,xyz,tess) +% usage: in = inhull(testpts,xyz,tess,tol) +% +% arguments: (input) +% testpts - nxp array to test, n data points, in p dimensions +% If you have many points to test, it is most efficient to +% call this function once with the entire set. +% +% xyz - mxp array of vertices of the convex hull, as used by +% convhulln. +% +% tess - tessellation (or triangulation) generated by convhulln +% If tess is left empty or not supplied, then it will be +% generated. +% +% tol - (OPTIONAL) tolerance on the tests for inclusion in the +% convex hull. You can think of tol as the distance a point +% may possibly lie outside the hull, and still be perceived +% as on the surface of the hull. Because of numerical slop +% nothing can ever be done exactly here. I might guess a +% semi-intelligent value of tol to be +% +% tol = 1.e-13*mean(abs(xyz(:))) +% +% In higher dimensions, the numerical issues of floating +% point arithmetic will probably suggest a larger value +% of tol. +% +% DEFAULT: tol = 0 +% +% arguments: (output) +% in - nx1 logical vector +% in(i) == 1 --> the i'th point was inside the convex hull. +% +% Example usage: The first point should be inside, the second out +% +% xy = randn(20,2); +% tess = convhulln(xy); +% testpoints = [ 0 0; 10 10]; +% in = inhull(testpoints,xy,tess) +% +% in = +% 1 +% 0 +% +% A non-zero count of the number of degenerate simplexes in the hull +% will generate a warning (in 4 or more dimensions.) This warning +% may be disabled off with the command: +% +% warning('off','inhull:degeneracy') +% +% See also: convhull, convhulln, delaunay, delaunayn, tsearch, tsearchn +% +% Author: John D'Errico +% e-mail: woodchips@rochester.rr.com +% Release: 3.0 +% Release date: 10/26/06 + +% get array sizes +% m points, p dimensions +p = size(xyz,2); +[n,c] = size(testpts); +if p ~= c + error 'testpts and xyz must have the same number of columns' +end +if p < 2 + error 'Points must lie in at least a 2-d space.' +end + +% was the convex hull supplied? +if (nargin<3) || isempty(tess) + tess = convhulln(xyz); +end +[nt,c] = size(tess); +if c ~= p + error 'tess array is incompatible with a dimension p space' +end + +% was tol supplied? +if (nargin<4) || isempty(tol) + tol = 0; +end + +% build normal vectors +switch p + case 2 + % really simple for 2-d + nrmls = (xyz(tess(:,1),:) - xyz(tess(:,2),:)) * [0 1;-1 0]; + + % Any degenerate edges? + del = sqrt(sum(nrmls.^2,2)); + degenflag = (del<(max(del)*10*eps)); + if sum(degenflag)>0 + warning('inhull:degeneracy',[num2str(sum(degenflag)), ... + ' degenerate edges identified in the convex hull']) + + % we need to delete those degenerate normal vectors + nrmls(degenflag,:) = []; + nt = size(nrmls,1); + end + case 3 + % use vectorized cross product for 3-d + ab = xyz(tess(:,1),:) - xyz(tess(:,2),:); + ac = xyz(tess(:,1),:) - xyz(tess(:,3),:); + nrmls = cross(ab,ac,2); + degenflag = false(nt,1); + otherwise + % slightly more work in higher dimensions, + nrmls = zeros(nt,p); + degenflag = false(nt,1); + for i = 1:nt + % just in case of a degeneracy + % Note that bsxfun COULD be used in this line, but I have chosen to + % not do so to maintain compatibility. This code is still used by + % users of older releases. + % nullsp = null(bsxfun(@minus,xyz(tess(i,2:end),:),xyz(tess(i,1),:)))'; + nullsp = null(xyz(tess(i,2:end),:) - repmat(xyz(tess(i,1),:),p-1,1))'; + if size(nullsp,1)>1 + degenflag(i) = true; + nrmls(i,:) = NaN; + else + nrmls(i,:) = nullsp; + end + end + if sum(degenflag)>0 + warning('inhull:degeneracy',[num2str(sum(degenflag)), ... + ' degenerate simplexes identified in the convex hull']) + + % we need to delete those degenerate normal vectors + nrmls(degenflag,:) = []; + nt = size(nrmls,1); + end +end + +% scale normal vectors to unit length +nrmllen = sqrt(sum(nrmls.^2,2)); +% again, bsxfun COULD be employed here... +% nrmls = bsxfun(@times,nrmls,1./nrmllen); +nrmls = nrmls.*repmat(1./nrmllen,1,p); + +% center point in the hull +center = mean(xyz,1); + +% any point in the plane of each simplex in the convex hull +a = xyz(tess(~degenflag,1),:); + +% ensure the normals are pointing inwards +% this line too could employ bsxfun... +% dp = sum(bsxfun(@minus,center,a).*nrmls,2); +dp = sum((repmat(center,nt,1) - a).*nrmls,2); +k = dp<0; +nrmls(k,:) = -nrmls(k,:); + +% We want to test if: dot((x - a),N) >= 0 +% If so for all faces of the hull, then x is inside +% the hull. Change this to dot(x,N) >= dot(a,N) +aN = sum(nrmls.*a,2); + +% test, be careful in case there are many points +in = false(n,1); + +% if n is too large, we need to worry about the +% dot product grabbing huge chunks of memory. +memblock = 1e6; +blocks = max(1,floor(n/(memblock/nt))); +aNr = repmat(aN,1,length(1:blocks:n)); +for i = 1:blocks + j = i:blocks:n; + if size(aNr,2) ~= length(j) + aNr = repmat(aN,1,length(j)); + end + in(j) = all((nrmls*testpts(j,:)' - aNr) >= -tol,1)'; +end + +end + diff --git a/sami_toolbox/+sami/+feat/+itx/synchronization.m b/sami_toolbox/+sami/+feat/+itx/synchronization.m new file mode 100644 index 0000000..907e5fc --- /dev/null +++ b/sami_toolbox/+sami/+feat/+itx/synchronization.m @@ -0,0 +1,126 @@ +function results = synchronization(data) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% options +step_length = .02; % in seconds +max_lag = 1; % in seconds +mw_size = .1; % size of moving window [in seconds] + +%% init vars +dt = 1/data.frameRate; % set time interval + +time_lag_int = step_length * data.frameRate; +time_lag_border = max_lag * data.frameRate; +time_lag = round(time_lag_int:time_lag_int:time_lag_border); + +marker_idx(1,:) = find(contains(data.labels, 'p1')); +marker_idx(2,:)= find(contains(data.labels, 'p2')); + +lbl_arms = {'SHO','ELB','WRI'}; +lbl_trunc = {'HEAD','HIP','KNE','ANK'}; +marker_idx_arms = contains(data.labels, lbl_arms); +marker_idx_arms = marker_idx_arms(contains(data.labels, 'p1')); +marker_idx_trunc = contains(data.labels, lbl_trunc); +marker_idx_trunc = marker_idx_trunc(contains(data.labels, 'p1')); + +corrMat_vel = nan( numel(time_lag)*2+1 , size(marker_idx,2) ); +corrMat_acc = nan( numel(time_lag)*2+1 , size(marker_idx,2) ); + +%% calculate velocity and acceleration of each marker +path_diff_3d = diff(data.marker); +path_diff = sqrt(sum(path_diff_3d.*path_diff_3d,3)); + +vel_raw = path_diff ./ dt; +vel = movmean(vel_raw,mw_size*data.frameRate); +acc_raw = diff(vel)./ dt; +acc = movmean(acc_raw,mw_size*data.frameRate); + +%% calculate correlations +% without time_shift +for m = 1:size(marker_idx,2) + corrMat_vel(numel(time_lag)+1,m) = corr(vel(:,marker_idx(1,m)), vel(:,marker_idx(2,m))); + corrMat_acc(numel(time_lag)+1,m) = corr(acc(:,marker_idx(1,m)), acc(:,marker_idx(2,m))); +end + +% correlations with time-shift +for m = 1:size(marker_idx,2) + for lag = 1:size(time_lag,2) + idx_plus = numel(time_lag) +1 + lag; + idx_minus = numel(time_lag) +1 - lag; + + % shift p1 -> left & p2 -> right + % velocity + tmp.time_shift_p1 = vel(time_lag(lag)+1:end,marker_idx(1,m)); + tmp.time_shift_p2 = vel(1:end-time_lag(lag),marker_idx(2,m)); + corrMat_vel(idx_plus,m) = corr(tmp.time_shift_p1, tmp.time_shift_p2); + + % acceleration + tmp.time_shift_p1 = acc(time_lag(lag)+1:end,marker_idx(1,m)); + tmp.time_shift_p2 = acc(1:end-time_lag(lag),marker_idx(2,m)); + corrMat_acc(idx_plus,m) = corr(tmp.time_shift_p1, tmp.time_shift_p2); + + % shift p1 -> right & p2 -> left + % velocity + tmp.time_shift_p1 = vel(1:end-time_lag(lag),marker_idx(1,m)); + tmp.time_shift_p2 = vel(time_lag(lag)+1:end,marker_idx(2,m)); + corrMat_vel(idx_minus,m) = corr(tmp.time_shift_p1, tmp.time_shift_p2); + + % acceleration + tmp.time_shift_p1 = acc(1:end-time_lag(lag),marker_idx(1,m)); + tmp.time_shift_p2 = acc(time_lag(lag)+1:end,marker_idx(2,m)); + corrMat_acc(idx_minus,m) = corr(tmp.time_shift_p1, tmp.time_shift_p2); + end +end + +% absolute correlations bc negative are also indicator for "*inter*action" +corrMat_vel = abs(corrMat_vel); +corrMat_acc = abs(corrMat_acc); + +% calc max +corr_vel = max(corrMat_vel); +corr_acc = max(corrMat_acc); + +%% results +results(1).feat = corr_vel; +results(1).name = 'sync_vel'; +results(1).color = [0 .8 .6]; +results(1).unit = '|Pearson r|'; +results(1).fisherZ4paramTesting = true; + +results(2).feat = corr_vel(marker_idx_arms); +results(2).name = 'sync_vel_arms'; +results(2).color = [0 .8 .6]; +results(2).unit = '|Pearson r|'; +results(2).fisherZ4paramTesting = true; + +results(3).feat = corr_vel(marker_idx_trunc); +results(3).name = 'sync_vel_trunc'; +results(3).color = [0 .8 .6]; +results(3).unit = '|Pearson r|'; +results(3).fisherZ4paramTesting = true; + +results(4).feat = corr_acc; +results(4).name = 'sync_acc'; +results(4).color = [0 .8 .6]; +results(4).unit = '|Pearson r|'; +results(4).fisherZ4paramTesting = true; + +results(5).feat = corr_acc(marker_idx_arms); +results(5).name = 'sync_acc_arms'; +results(5).color = [0 .8 .6]; +results(5).unit = '|Pearson r|'; +results(5).fisherZ4paramTesting = true; + +results(6).feat = corr_acc(marker_idx_trunc); +results(6).name = 'sync_acc_trunc'; +results(6).color = [0 .8 .6]; +results(6).unit = '|Pearson r|'; +results(6).fisherZ4paramTesting = true; + +end diff --git a/sami_toolbox/+sami/+fig/MDSArrangement.m b/sami_toolbox/+sami/+fig/MDSArrangement.m new file mode 100644 index 0000000..9981ecd --- /dev/null +++ b/sami_toolbox/+sami/+fig/MDSArrangement.m @@ -0,0 +1,130 @@ +function MDSArrangement(RDM, MDSOptions) +% MDSArrangement(RDM, MDSOptions) +% +% draws a multidimensional scaling (MDS) arrangement of colored dots reflecting the +% dissimilarity structure of the items whose representational dissimilarity matrix is +% passed in argument RDM. +% +% input: +% - RDM: +% struct. dissimilarity matrix as a struct RDM, containg 'RDM' and 'name' field. +% +% - MDSOpions.fig_display: +% boolean. +% defaults to true. +% +% - MDSOpions.figI: +% double. number of current figure handle. +% defaults to "next available number". +% +% - MDSOpions._____: +% boolean. +% defaults to true. +% +% - MDSOpions._____: +% boolean. +% defaults to true. +% +% - MDSOpions._____: +% boolean. +% defaults to true. +% +% - MDSOpions._____: +% boolean. +% defaults to true. +% +% - MDSOpions._____: +% boolean. +% defaults to true. +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 03/2021 + + +%% define defaults +if ~exist('MDSOptions','var') || isempty(MDSOptions), MDSOptions = struct; end +MDSOptions = sami.util.setIfUnset(MDSOptions, 'fig_display',true); +MDSOptions = sami.util.setIfUnset(MDSOptions, 'figI',sami.util.getFigI()); +MDSOptions = sami.util.setIfUnset(MDSOptions, 'MDSCriterion', 'metricstress'); +MDSOptions = sami.util.setIfUnset(MDSOptions, 'plotLabels', false); +MDSOptions = sami.util.setIfUnset(MDSOptions, 'plotLegend', false); +MDSOptions = sami.util.setIfUnset(MDSOptions, 'dotSize', 20); +MDSOptions = sami.util.setIfUnset(MDSOptions, 'dotColors', repmat([0 0 0],size(RDM.RDM,1),1)); +MDSOptions = sami.util.setIfUnset(MDSOptions, 'fontSize', 9); +MDSOptions = sami.util.setIfUnset(MDSOptions, 'titleString', sami.util.deunderscore(RDM.name)); + +%% perform multidimensional scaling (MDS) +D = sami.rdm.unwrapRDMs(RDM); +try + xyCoords = mdscale(D, 2,'criterion',MDSOptions.MDSCriterion, 'options', struct('MaxIter', 100000)); +catch + try + xyCoords = mdscale(D, 2,'criterion','stress'); + disp([' -> MDS Info: ' RDM.name ' -> reverted to stress (' MDSOptions.MDSCriterion ' failed)']); + catch + try + D2 = D + 0.2; + D2(logical(eye(length(D)))) = 0; + xyCoords = mdscale(D2, 2,'criterion',MDSOptions.MDSCriterion); + disp([' -> MDS Info: ' RDM.name ', ' MDSOptions.MDSCriterion ' , added 0.2 to distances to avoid colocalization']); + catch + disp([' -> MDS Info: ' RDM.name ', MDS failed...']); + return + end + end +end + +%% plot MDS arrangement using text labels +sami.util.selectPlot(MDSOptions.figI,MDSOptions.fig_display); +set(gcf,'Position',[200 50 750 580]); +hold on; + +% plot DOTS +for itemI = 1:size(D,1) + plot(xyCoords(itemI,1), xyCoords(itemI,2),'o',... + 'MarkerFaceColor',MDSOptions.dotColors(itemI, :),... + 'MarkerEdgeColor','none',... + 'MarkerSize', MDSOptions.dotSize); +end + +% fix axis limits for 1-dimensional MDS-plot +axLim = [min(xyCoords(:,1)) 1.1*max(xyCoords(:,1)) min(xyCoords(:,2)) 1.1*max(xyCoords(:,2))]; +if isequal(axLim([3 4]),[0 0]) + axLim([3 4]) = [-1 3]; +end + +% edit axis +axis square equal off; +axis(axLim); +title({['\fontsize{11}' MDSOptions.titleString],... + ['\fontsize{10}distance measure:\rm ' MDSOptions.MDSDistance ' | \bfcriterion:\rm ' MDSOptions.MDSCriterion]}); + +% plot legend +if MDSOptions.plotLegend == 1 + h = flipud(get(gca,'Children')); % grab all the axes handles at once + l = legend(h(MDSOptions.legend.DotIdx),MDSOptions.legend.Labels,'Location','EastOutside'); + if ~isempty(MDSOptions.legend.title) + title(l,MDSOptions.legend.title); + end +end + +% plot text labels in black +% --> shift y-coordinate of text -> "MDSOptions.dotSize" in point-unit +df = .03; +dx = range(xlim) * df * 0; +dy = range(ylim) * df * -1; + +if MDSOptions.plotLabels == 1 + for itemI = 1:size(D,1) + text(xyCoords(itemI,1) + dx, xyCoords(itemI,2) + dy,... + sami.util.deunderscore(MDSOptions.dotLabels{itemI}),... + 'Color','k',... + 'FontSize',MDSOptions.fontSize,... + 'VerticalAlignment','top',... + 'HorizontalAlignment','center'); + end +end + +end + diff --git a/sami_toolbox/+sami/+fig/createColormap.m b/sami_toolbox/+sami/+fig/createColormap.m new file mode 100644 index 0000000..3162abe --- /dev/null +++ b/sami_toolbox/+sami/+fig/createColormap.m @@ -0,0 +1,50 @@ +function cols = createColormap(style) +% cols = createColormap(style) +% +% calculates a customized colormap given the requested 'style', which is defined in +% the samiOptions.m file +% +% input: +% - style: +% string. Defining which field within samiOtions.cmap.'style' is used to +% create colormap +% +% output: +% - cols: +% colormap in RGB values (dimensions [m x 3]) +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 03/2021 + +%% loading samiOptions +samiOptions = sami.loadSamiOptions(); +nCols = samiOptions.cmap.nCols; + +%% get positions+colors for requested style +positions = samiOptions.cmap.(style).pos; +colors = samiOptions.cmap.(style).col; + +%% computations +% compute positions along the samples +colSamp = round((nCols-1)*positions)+1; + +% make the gradients among colors +cols = zeros(nCols,3); +cols(colSamp,:) = colors; +diffS = diff(colSamp)-1; +for d = 1:1:length(diffS) + if diffS(d)~=0 + col1 = colors(d,:); + col2 = colors(d+1,:); + G = zeros(diffS(d),3); + for idx = 1:3 + g = linspace(col1(idx), col2(idx), diffS(d)+2); + g([1, length(g)]) = []; + G(:,idx) = g'; + end + cols(colSamp(d)+1:colSamp(d+1)-1,:) = G; + end +end +end + + diff --git a/sami_toolbox/+sami/+fig/drawBars.m b/sami_toolbox/+sami/+fig/drawBars.m new file mode 100644 index 0000000..4d896e7 --- /dev/null +++ b/sami_toolbox/+sami/+fig/drawBars.m @@ -0,0 +1,51 @@ +function drawBars(d,category,userOptions) +% drawBars(d,category,userOptions) +% +% draws a bar plot (bar + errorbar), overlayed with dots representing each datapoint. +% +% input: +% - d: +% double. [d*g]-matrix containg 'n' datapoints for each of the 'g' groups +% +% - category: +% string. the category by which the groups were formed +% +% - userOptions.stimuli_naming_key: +% struct. contains labels for each group, according to the respective category +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 10/2020 + +% plotting options +dot_spread = .3; +dot_col = [.4 .4 .4]; +dot_size = 16; +bar_col = [.8 .8 .8]; +xAxisWhiteSpace = .6; + +% preparing +dMN = mean(d); +dSEM = std(d) ./ sqrt(size(d,1)); +[~, cI] = ismember(category,{userOptions.stimuli_naming_key.name}); + +% bar plot +bar(dMN,'FaceColor',bar_col); hold on; + +% dots +for i = 1:numel(dMN) + tmpXplus1 = mod(size(d,1),2); + pointXs = [linspace(-dot_spread,-0.1,floor(size(d,1)/2)) ,... + linspace(0.1,dot_spread,floor(size(d,1)/2)+tmpXplus1)]; + plot(pointXs+i,d(:,i),'.','Color',dot_col,'MarkerSize',dot_size); +end + +% errorbars +errorbar(1:numel(dMN),dMN,dSEM,'LineStyle','none','CapSize',0,'LineWidth',3,'Color','k') + +% edit axis +set(gca,'Xtick',1:numel(dMN),'XTickLabel',userOptions.stimuli_naming_key(cI).condition); +h=gca; h.XAxis.TickLength = [0 0]; +set(gca, 'YGrid', 'on', 'XGrid', 'off') +xlim([1-xAxisWhiteSpace numel(dMN)+xAxisWhiteSpace]); + +end diff --git a/sami_toolbox/+sami/+fig/drawBoxPlot.m b/sami_toolbox/+sami/+fig/drawBoxPlot.m new file mode 100644 index 0000000..1fa94c8 --- /dev/null +++ b/sami_toolbox/+sami/+fig/drawBoxPlot.m @@ -0,0 +1,53 @@ +function drawBoxPlot(d,category,userOptions) +% drawBoxPlot(d,category,userOptions) +% +% draws a boxplots, overlayed with dots representing each datapoint. +% +% input: +% - d: +% double. [d*g]-matrix containg 'n' datapoints for each of the 'g' groups +% +% - category: +% string. the category by which the groups were formed +% +% - userOptions.stimuli_naming_key: +% struct. contains labels for each group, according to the respective category +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 10/2020 + +% plotting options +dot_spread = .2; +dot_col = [.4 .4 .4]; +dot_size = 20; +xAxisWhiteSpace = .6; +box_col = [.8 .8 .8]; + +% preparing +[~, cI] = ismember(category,{userOptions.stimuli_naming_key.name}); + +% first box plot to obtain x-and-y-data +boxplot(d,'colors',[0 0 0],'Labels',userOptions.stimuli_naming_key(cI).condition,'Notch','off','Widths',0.6,'OutlierSize',0.001,'Symbol','w.'); +hold on; + +% patches +h = findobj(gca,'Tag','Box'); +for j=1:length(h) + patch(get(h(j),'XData'),get(h(j),'YData'),box_col,'FaceAlpha',1); +end + +% dots +for i = 1:size(d,2) + pointXs = linspace(-dot_spread,dot_spread,size(d,1)); + plot(pointXs+i,d(:,i),'.','Color',dot_col,'MarkerSize',dot_size); +end + +% second boxplot on the top +boxplot(d,'colors',[0 0 0],'Labels',userOptions.stimuli_naming_key(cI).condition,'Notch','off','Widths',0.6,'OutlierSize',0.001,'Symbol','w.'); + +% edit stuff +h=gca; h.XAxis.TickLength = [0 0]; +set(gca, 'YGrid', 'on', 'XGrid', 'off') +xlim([1-xAxisWhiteSpace size(d,2)+xAxisWhiteSpace]); + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+fig/drawRefenceTestLine.m b/sami_toolbox/+sami/+fig/drawRefenceTestLine.m new file mode 100644 index 0000000..61efa0e --- /dev/null +++ b/sami_toolbox/+sami/+fig/drawRefenceTestLine.m @@ -0,0 +1,34 @@ +function drawRefenceTestLine(m) +% drawRefenceTestLine(m) +% +% draws dashed line for testValue, against groups are tested for a significant difference +% +% input: +% - m: +% double. test value, against tests are performed +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 03/2021 + +line_col = [0 0 0]; +line_width = 2; +line_style = '--'; + +if m ~= 0 + curXLim = xlim; + curYLim = ylim; + curYTicks = get(gca,'ytick'); + line(curXLim,[m m],'Color',line_col,'LineWidth',line_width,'LineStyle',line_style); + + % adjust y-axis-limits and yTicks + if m < curYLim(1) + curYLim(1) = m * (4/5); + ylim([curYLim(1) curYTicks(end)]); + set(gca, 'YTickMode', 'auto', 'YTickLabelMode', 'auto'); + curYTicks = get(gca,'YTick'); + ylim(curYLim); + set(gca,'YTick',curYTicks); + end +end +end + diff --git a/sami_toolbox/+sami/+fig/drawSigAsterisks.m b/sami_toolbox/+sami/+fig/drawSigAsterisks.m new file mode 100644 index 0000000..1ca3ee2 --- /dev/null +++ b/sami_toolbox/+sami/+fig/drawSigAsterisks.m @@ -0,0 +1,52 @@ +function drawSigAsterisks(inpX,inpP) +% drawSigAsterisks(inpX,inpP) +% +% draws asterisks at the top of the axis, indicating significance for one-sample +% comparisons agains a testValue +% +% input: +% - inpX: +% double. index (= location on x-axis) of group which is compared against a +% specific value +% +% - inpP: +% double. pValues for each group, order according to inpX +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 03/2021 + +ast_col = [0 0 0]; +ast_size = 18; +line_col = [0 0 0]; +line_width = 1; +line_style = ':'; + +yAxisFactorLimit = .05; +yAxisFactorAsterisk = .055; + +if any(inpP < 0.05) + % get yCoords of sig. lines + curYTicks = get(gca,'YTick'); + curYLim = ylim; + curXLim = xlim; + + yMarks = curYLim(2) + range(curYLim)*yAxisFactorLimit .* [1 2]; + ylim([curYLim(1) yMarks(2)]); + set(gca,'YTick',curYTicks); + + astH = curYLim(2) + range(curYLim)*yAxisFactorAsterisk; + + % plot line and asterisks + for i = 1:numel(inpP) + if inpP(i) < 0.05 + text(inpX(i),astH,'*',... + 'FontSize',ast_size,'Color',ast_col,... + 'VerticalAlignment','middle',... + 'HorizontalAlignment','center'); + end + end + line(curXLim,[yMarks(1) yMarks(1)],... + 'Color',line_col,'LineWidth',line_width,'LineStyle',line_style); +end +end + diff --git a/sami_toolbox/+sami/+fig/drawSigLines.m b/sami_toolbox/+sami/+fig/drawSigLines.m new file mode 100644 index 0000000..ff352fc --- /dev/null +++ b/sami_toolbox/+sami/+fig/drawSigLines.m @@ -0,0 +1,48 @@ +function drawSigLines(inpX,inpP) +% drawSigLines(inpX,inpP) +% +% draws lines above bar/box-plots for significant post-hoc pairwise comparisons +% +% input: +% - inpX: +% double. a [p x 2]-matrix, indicating for p pairwise comparisons the group +% which are compared against each other (corresponding location on x-axis) +% +% - inpP: +% double. pValues for p pairwise comparisons, order according to rows in inpX +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 03/2021 + + +line_col = [0 0 0]; +line_width = 1; +line_compactness = 2; % integer (!) values from intervall [0,9] makes sense + +addYAxisWhiteSpaceFactor = .3; + +if any(inpP < 0.05) + % get yCoords of sig. lines + curYTicks = get(gca,'YTick'); + curYLim = ylim; + ylim([curYLim(1) curYLim(2)+range(curYLim)*addYAxisWhiteSpaceFactor]); + set(gca,'YTick',curYTicks); + + linH2 = linspace(curYLim(2),curYLim(2)+range(curYLim)*addYAxisWhiteSpaceFactor,size(inpP,1)+2+(2*line_compactness)); + linSpread = range(linH2(1:2)); + linH = linH2(2+line_compactness:end-1-line_compactness); + + % show pairwise sign. lines + tmpI = 0; + for i = 1:size(inpP,1) + if inpP(i) < 0.05 + tmpI = tmpI + 1; + line([inpX(i,1) inpX(i,2)],[linH(tmpI) linH(tmpI)],'Color',line_col,'LineWidth',line_width); + end + end + + % update yLim and get rid of white_space above lines + ylim([curYLim(1) linH(tmpI)+2*linSpread]); +end +end + diff --git a/sami_toolbox/+sami/+fig/handleFigure.m b/sami_toolbox/+sami/+fig/handleFigure.m new file mode 100644 index 0000000..d58666a --- /dev/null +++ b/sami_toolbox/+sami/+fig/handleFigure.m @@ -0,0 +1,78 @@ +function handleFigure(figI, fileName, userOptions) +% handleFigure(figI, fileName, userOptions) +% +% handleFigure is a function which will, based on the preferences set in userOptions, +% save the current figure as a .pdf, a .tif, a .svg, or a .fig; and also either leaving +% it open or closing it. +% +% input: +% - figI: +% figure handle. The figure handle of the figure to be saved. +% +% - fileName: +% string. The name of the file to be saved. +% +% - userOptions.fig_display +% boolean. If true, the figure remains open after it is handles. +% Defaults to true. +% +% - userOptions.fig_savePDF +% boolean. If true, the figure is saved as a PDF. +% Defaults to false. +% +% - userOptions.fig_saveTIF +% A boolean value. If true, the figure is saved as a TIF. +% Defaults to false. +% +% - userOptions.fig_saveSVG +% A boolean value. If true, the figure is saved as a SVG. +% Defaults to false. +% +% - userOptions.fig_saveFIG +% A boolean value. If true, the figure is saved as a MATLAB .fig file. +% Defaults to false. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% Set defaults and check options struct +userOptions = sami.util.setIfUnset(userOptions, 'fig_display', true); +userOptions = sami.util.setIfUnset(userOptions, 'fig_savePDF', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_saveSVG', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_saveTIF', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_saveFIG', false); + +%% set PaperSize to FigSize +figHandles = get(groot, 'Children'); +fh = figHandles([figHandles(:).Number] == figI); +fh.PaperPositionMode = 'auto'; +fig_pos = fh.PaperPosition; +fh.PaperSize = [fig_pos(3) fig_pos(4)]; + +%% saving +% saving PDF +if userOptions.fig_savePDF + print('-dpdf',fileName,'-bestfit'); +end + +% saving SVG +if userOptions.fig_saveSVG + print('-dsvg',sprintf('-r%d',userOptions.fig_dpi),fileName); +end + +% saving TIF +if userOptions.fig_saveTIF + print('-dtiff',sprintf('-r%d',userOptions.fig_dpi),fileName); +end + +% saving Matlab figure +if userOptions.fig_saveFIG + savefig(figI,fileName); +end + +% closing figure +if ~userOptions.fig_display + close(figI); +end + +end diff --git a/sami_toolbox/+sami/+rdm/averageRDMs.m b/sami_toolbox/+sami/+rdm/averageRDMs.m new file mode 100644 index 0000000..555cddd --- /dev/null +++ b/sami_toolbox/+sami/+rdm/averageRDMs.m @@ -0,0 +1,38 @@ +function avgRDM = averageRDMs(RDMs,new_name,new_color) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% preparations +if ~exist('new_name','var') || isempty(new_name), new_name = 'averageRDM'; end +if ~exist('new_color','var') || isempty(new_color), new_color = RDMs(1).color; end + +%% init +n = numel(RDMs); +suSum = nan(size(RDMs(1).RDM,1),size(RDMs(1).RDM,1),n); + +%% average +% for su = 1:n +% if su == 1 +% suSum = RDMs(su).RDM; +% else +% suSum = suSum + RDMs(su).RDM; +% end%if:su==1 +% end%for:su +% +% avgRDM.RDM = suSum ./ n; + +for su = 1:n + suSum(:,:,su) = RDMs(su).RDM; +end +avgRDM.RDM = nanmean(suSum,3); + +avgRDM.name = new_name; +avgRDM.color = new_color; + +end + diff --git a/sami_toolbox/+sami/+rdm/categoricalRDM.m b/sami_toolbox/+sami/+rdm/categoricalRDM.m new file mode 100644 index 0000000..a5b747e --- /dev/null +++ b/sami_toolbox/+sami/+rdm/categoricalRDM.m @@ -0,0 +1,38 @@ +function [binRDM, nCatCrossingsRDM] = categoricalRDM(categoryVectors, figI, monitor) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +import sami.* +import sami.fig.* +import sami.rdm.* +import sami.util.* + +%% preparations +if ~exist('monitor','var'), monitor = false; end +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +if min(size(categoryVectors)) == 1, categoryVectors = categoryVectors(:); end + +[nCond, nCats] = size(categoryVectors); + +%% count category crossings for each pair of conditions +nCatCrossingsRDM = zeros(nCond,nCond); + +for catI=1:nCats + cCatBinRDM = repmat(categoryVectors(:,catI),[1 nCond]) ~= repmat(categoryVectors(:,catI)',[nCond 1]); + nCatCrossingsRDM = nCatCrossingsRDM + cCatBinRDM; +end + +binRDM = double(logical(nCatCrossingsRDM)); + +%% visualise +if monitor + plotRDMs(concatRDMs(binRDM,nCatCrossingsRDM),figI); +end + + +end%function diff --git a/sami_toolbox/+sami/+rdm/concatRDMs.m b/sami_toolbox/+sami/+rdm/concatRDMs.m new file mode 100644 index 0000000..ec759ca --- /dev/null +++ b/sami_toolbox/+sami/+rdm/concatRDMs.m @@ -0,0 +1,20 @@ +function RDMs=concatRDMs(varargin) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +RDMs = []; +for RDMI = 1:nargin + + if ~isstruct(varargin{RDMI}) + varargin{RDMI} = sami.rdm.wrapRDMs(varargin{RDMI}); + end + + RDMs = [RDMs, varargin{RDMI}]; +end + +end%function \ No newline at end of file diff --git a/sami_toolbox/+sami/+rdm/nameRDMs.m b/sami_toolbox/+sami/+rdm/nameRDMs.m new file mode 100644 index 0000000..afd1f73 --- /dev/null +++ b/sami_toolbox/+sami/+rdm/nameRDMs.m @@ -0,0 +1,14 @@ +function RDMs = nameRDMs(RDMs,names) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +for nameI = 1:numel(names) + RDMs(nameI).name = names{nameI}; +end + +end diff --git a/sami_toolbox/+sami/+rdm/unwrapRDMs.m b/sami_toolbox/+sami/+rdm/unwrapRDMs.m new file mode 100644 index 0000000..e72e4f6 --- /dev/null +++ b/sami_toolbox/+sami/+rdm/unwrapRDMs.m @@ -0,0 +1,37 @@ +function [RDMs, nRDMs, RDMNames] = unwrapRDMs(RDMs_struct) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +nonameIdx = 1; + +if isstruct(RDMs_struct) + % in struct form + nRDMs = size(RDMs_struct,2); + + RDMNames = cell(nRDMs,1); + RDMs = nan(size(RDMs_struct(1).RDM,1),size(RDMs_struct(1).RDM,2),nRDMs); + + for i = 1:nRDMs + RDMs(:,:,i) = double(RDMs_struct(i).RDM); + try + RDMNames{i} = RDMs_struct(i).name; + catch + RDMNames{i} = ['unnamed_RDM',num2str(nonameIdx)]; + nonameIdx = nonameIdx + 1; + end + end +else + % bare already + RDMs = RDMs_struct; + nRDMs = size(RDMs,3); + + RDMNames = cell(nRDMs,1); + [RDMNames{:}] = deal('unnamed_RDM'); +end + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+rdm/vectorizeRDM.m b/sami_toolbox/+sami/+rdm/vectorizeRDM.m new file mode 100644 index 0000000..4362305 --- /dev/null +++ b/sami_toolbox/+sami/+rdm/vectorizeRDM.m @@ -0,0 +1,16 @@ +function RDM = vectorizeRDM(RDM) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + + +if size(RDM,1)==size(RDM,2) + RDM(logical(eye(size(RDM))))=0; % fix diagonal: zero by definition + RDM = squareform(RDM); +end + +end%function diff --git a/sami_toolbox/+sami/+rdm/vectorizeRDMs.m b/sami_toolbox/+sami/+rdm/vectorizeRDMs.m new file mode 100644 index 0000000..e8f5693 --- /dev/null +++ b/sami_toolbox/+sami/+rdm/vectorizeRDMs.m @@ -0,0 +1,32 @@ +function ltv_RDMs = vectorizeRDMs(RDMs) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +if isstruct(RDMs) + % wrapped + nRDMs = size(RDMs,2); + ltv_RDMs = []; + + for iRDM = 1:nRDMs + thisRDM = RDMs(iRDM).RDM; + ltv_RDMs = cat(3,ltv_RDMs,sami.rdm.vectorizeRDM(thisRDM)); + end + + ltv_RDMs = sami.rdm.wrapRDMs(ltv_RDMs,RDMs); +else + % bare + nRDMs = size(RDMs,3); + ltv_RDMs = []; + + for iRDM = 1:nRDMs + ltv_RDMs = cat(3,ltv_RDMs,sami.rdm.vectorizeRDM(RDMs(:,:,iRDM))); + end +end + + +end%function diff --git a/sami_toolbox/+sami/+rdm/wrapAndNameRDMs.m b/sami_toolbox/+sami/+rdm/wrapAndNameRDMs.m new file mode 100644 index 0000000..84cd435 --- /dev/null +++ b/sami_toolbox/+sami/+rdm/wrapAndNameRDMs.m @@ -0,0 +1,13 @@ +function RDMs = wrapAndNameRDMs(RDMs,names) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +RDMs = sami.rdm.wrapRDMs(RDMs); +RDMs = sami.rdm.nameRDMs(RDMs,names); + +end diff --git a/sami_toolbox/+sami/+rdm/wrapRDMs.m b/sami_toolbox/+sami/+rdm/wrapRDMs.m new file mode 100644 index 0000000..ebfcc21 --- /dev/null +++ b/sami_toolbox/+sami/+rdm/wrapRDMs.m @@ -0,0 +1,37 @@ +function RDMs_struct = wrapRDMs(RDMs,refRDMs_struct) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +if isstruct(RDMs) + % wrapped already, but replace the wrapping + nRDMs = numel(RDMs); +else + nRDMs = size(RDMs,3); +end + +if ~exist('refRDMs_struct','var') + for iRDM = 1:nRDMs + refRDMs_struct(iRDM).name = '[unnamed RDM]'; + refRDMs_struct(iRDM).color = [0 0 0]; + end +end + +RDMs_struct = refRDMs_struct; +if isstruct(RDMs) + % wrapped already, but replace the wrapping + for iRDM = 1:nRDMs + RDMs_struct(iRDM).RDM = RDMs(iRDM).RDM; + end +else + % RDMs need wrapping + for iRDM = 1:nRDMs + RDMs_struct(iRDM).RDM = RDMs(:,:,iRDM); + end +end + +end%function diff --git a/sami_toolbox/+sami/+stat/ANOVA.m b/sami_toolbox/+sami/+stat/ANOVA.m new file mode 100644 index 0000000..0dca72a --- /dev/null +++ b/sami_toolbox/+sami/+stat/ANOVA.m @@ -0,0 +1,33 @@ +function results = ANOVA(d,CorrMethod,alpha) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +% ANOVA +[aR.p,aR.tbl,aR.stats] = anova1(d,[],'off'); +aR.comp = multcompare(aR.stats,'display','off','ctype',CorrMethod,'alpha',alpha); + +% plot horizontal sign. lines if sign. pairwise comparisons exist +sami.fig.drawSigLines(aR.comp(:,1:2),aR.comp(:,6)) + +% calculate eta^2 +etaSq = aR.tbl{2,2} / aR.tbl{4, 2}; + +% sorting ANOVA results +results.test = 'ANOVA'; +results.S = 'F'; +results.V = aR.tbl{2,5}; +results.df1 = aR.tbl{2,3}; +results.df2 = aR.tbl{3,3}; +results.p = aR.tbl{2,6}; +results.etaSq = etaSq; + +% store anova- and multCompare-tables +results.tbl = aR.tbl; +results.multComp = aR.comp; + +end diff --git a/sami_toolbox/+sami/+stat/FDR.m b/sami_toolbox/+sami/+stat/FDR.m new file mode 100644 index 0000000..b57fa71 --- /dev/null +++ b/sami_toolbox/+sami/+stat/FDR.m @@ -0,0 +1,48 @@ +function [pCorr, pCrit] = FDR(pVals,alpha) +% +% +% +% References: +% Benjamini, Y. & Hochberg, Y. (1995) Controlling the false discovery rate: A practical +% and powerful approach to multiple testing. Journal of the Royal Statistical +% Society, Series B (Methodological). 57(1), 289-300. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +if nargin < 2 + alpha = .05; +end + +% sort p-values +[pSorted, sortIdx] = sort(pVals); +n = length(pSorted); %number of tests + +% BH procedure +thresh = (1:n)'*alpha/n; +pW = n*pSorted./(1:n)'; + + +% compute adjusted p-values +pCorr = nan(n,1); +[pWsorted, pWindex] = sort(pW); +runI = 1; +for k = 1 : n + if pWindex(k) >= runI + pCorr(runI:pWindex(k)) = pWsorted(k); + runI = pWindex(k)+1; + if runI > n + break; + end + end +end +pCorr(sortIdx,1) = pCorr; + +% find greatest significant pvalue +r = pSorted <= thresh; +maxID = find(r,1,'last'); +if isempty(maxID) + pCrit = 0; +else + pCrit = pSorted(maxID); +end diff --git a/sami_toolbox/+sami/+stat/Friedman.m b/sami_toolbox/+sami/+stat/Friedman.m new file mode 100644 index 0000000..192c93f --- /dev/null +++ b/sami_toolbox/+sami/+stat/Friedman.m @@ -0,0 +1,31 @@ +function results = Friedman(d,CorrMethod,alpha) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +% KruskalWallis +[F.p,F.tbl,F.stats] = friedman(d,1,'off'); +F.comp = multcompare(F.stats,'display','off','ctype',CorrMethod,'alpha',alpha); + +% plot horizontal sign. lines if sign. pairwise comparisons exist +sami.fig.drawSigLines(F.comp(:,1:2),F.comp(:,6)) + +% calculate effect size +KendallW = F.tbl{2,5} / (size(d,1) * (size(d,2)-1)); + +% sorting ANOVA results +results.test = 'Friedman'; +results.S = 'Chi^2'; +results.V = F.tbl{2,5}; +results.df = F.tbl{2,3}; +results.p = F.tbl{2,6}; +results.KendallW = KendallW; + +results.tbl = F.tbl; +results.multComp = F.comp; + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+stat/KruskalWallis.m b/sami_toolbox/+sami/+stat/KruskalWallis.m new file mode 100644 index 0000000..b7577d6 --- /dev/null +++ b/sami_toolbox/+sami/+stat/KruskalWallis.m @@ -0,0 +1,31 @@ +function results = KruskalWallis(d,CorrMethod,alpha) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +% KruskalWallis +[KW.p,KW.tbl,KW.stats] = kruskalwallis(d,[],'off'); +KW.comp = multcompare(KW.stats,'display','off','ctype',CorrMethod,'alpha',alpha); + +% plot horizontal sign. lines if sign. pairwise comparisons exist +sami.fig.drawSigLines(KW.comp(:,1:2),KW.comp(:,6)) + +% calculate effect size +etaSq = (KW.tbl{2,5} - size(d,2) + 1) / (numel(d) - size(d,2)); + +% sorting ANOVA results +results.test = 'Kruskal-Wallis'; +results.S = 'H'; +results.V = KW.tbl{2,5}; +results.df = KW.tbl{2,3}; +results.p = KW.tbl{2,6}; +results.etaSq = etaSq; + +results.tbl = KW.tbl; +results.multComp = KW.comp; + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+stat/RDMCorrMat.m b/sami_toolbox/+sami/+stat/RDMCorrMat.m new file mode 100644 index 0000000..af4f778 --- /dev/null +++ b/sami_toolbox/+sami/+stat/RDMCorrMat.m @@ -0,0 +1,34 @@ +function [corrMat, pValMat] = RDMCorrMat(RDMs,type) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + + +if ~exist('type','var') || isempty(type), type = 'Spearman'; end + +nRDMs = size(RDMs,2); +RDMs_cols = sami.rdm.unwrapRDMs(sami.rdm.vectorizeRDMs(RDMs)); +RDMs_cols = permute(RDMs_cols,[2 3 1]); + +% For each pair of RDMs, ignore missing data only for this pair of RDMs +% (unlike just using corr, which would ignore it if ANY RDM had missing +% data at this point). +corrMat = nan(nRDMs); +pValMat = nan(nRDMs); +for iRDM = 1:nRDMs + for jRDM = 1 : nRDMs + if isequal(type,'Kendall_taua') + [corrMat(iRDM,jRDM),pValMat(iRDM,jRDM)] = sami.stat.rankCorr_Kendall_taua(RDMs_cols(:,iRDM), RDMs_cols(:,jRDM)); + else + [corrMat(iRDM,jRDM),pValMat(iRDM,jRDM)] = corr(RDMs_cols(:,iRDM), RDMs_cols(:,jRDM), 'type', type, 'rows', 'complete'); + end + end +end + +corrMat(logical(eye(nRDMs))) = 1; % make the diagonal artificially one + +end diff --git a/sami_toolbox/+sami/+stat/Wilcoxon.m b/sami_toolbox/+sami/+stat/Wilcoxon.m new file mode 100644 index 0000000..b47a049 --- /dev/null +++ b/sami_toolbox/+sami/+stat/Wilcoxon.m @@ -0,0 +1,118 @@ +function results = Wilcoxon(d,type,m) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 03/2021 + +%% define default behavior +if ~exist('m','var') || isempty(m), m = 0; end + +%% perform appropriate Wilcoxon, as specified by 'type' +switch type + case 'one_sample' + % naming t-Test results + results.test = 'one-sample Wilcoxon signed rank test'; + + % Wilcoxon test + for i=1:size(d,2) + [WX.p(i),~,WX.stats(i)] = signrank(d(:,i),m); + + % calculate effect size as matched rank-biserial correlation + dDiff = d(:,i) - m; + + epsdiff = eps(d(:,i)); + t = (abs(dDiff) <= epsdiff); + dDiff(t) = []; + epsdiff(t) = []; + + tr = tiedrank(abs(dDiff),0,0,epsdiff); + W(1) = sum(tr(dDiff>0)); + W(2) = sum(tr(dDiff<0)); + + n = size(dDiff,1); + WX.r(i) = 4 * abs( min(W) - mean(W) ) / (n*(n+1)); + end + + % bonferonni correction of p-values + if size(d,2) > 1 + results.info = ['p values bonferonni corrected by ' num2str(size(d,2)) ' indepedent tests']; + WX.p = WX.p .* size(d,2); + end + + % draw line, if tested against ~=0 + sami.fig.drawRefenceTestLine(m); + + % plot asterisks if sign. difference from 'm' + sami.fig.drawSigAsterisks(1:size(d,2),WX.p); + + case 'two_sample' + % naming t-Test results + results.test = 'Wilcoxon rank sum test'; + + % two-sided Wilcoxon rank sum test for independant samples + [WX.p,~,WX.stats] = ranksum(d(:,1),d(:,2)); + + % calculate effect size as rank biserial correlation + n = size(d,1); + tdrnk = tiedrank(d(:)); + WX.r = 2*(mean(tdrnk(1:n))-mean(tdrnk(n+1:end))) / numel(d); + + % plot horizontal sign. lines if sign. pairwise comparisons exist + sami.fig.drawSigLines([1 2],WX.p); + + case 'paired' + % naming t-Test results + results.test = 'paired two-sample Wilcoxon signed rank test'; + + % Wilcoxon test + [WX.p,~,WX.stats] = signrank(d(:,1),d(:,2)); + + % calculate effect size as matched rank-biserial correlation + dDiff = d(:,1)-d(:,2); + + epsdiff = eps(d(:,1)) + eps(d(:,2)); + t = (abs(dDiff) <= epsdiff); + dDiff(t) = []; + epsdiff(t) = []; + + tr = tiedrank(abs(dDiff),0,0,epsdiff); + W(1) = sum(tr(dDiff>0)); + W(2) = sum(tr(dDiff<0)); + + n = size(dDiff,1); + WX.r = 4 * (min(W) - mean(W)) / (n*(n+1)); + + % correct sign of r + WX.r = WX.r .* sign(diff(W)); + + %-------- to compare with JASP results, which makes some rounding errors --------- +% dDiff = diff(round(d,3),[],2); +% +% t = (abs(dDiff) == 0); +% dDiff(t) = []; +% +% tr = tiedrank(abs(dDiff)); +% W(1) = sum(tr(dDiff<0)); +% W(2) = sum(tr(dDiff>0)); +% +% n = size(dDiff,1); +% WX.r_JASP = 4 * ( min(W) - mean(W) ) / (n*(n+1)); + %--------------------------------------------------------------------------------- + + % plot horizontal sign. lines if sign. pairwise comparisons exist + sami.fig.drawSigLines([1 2],WX.p) + + otherwise + error('*** sami:ERROR *** Wilcoxon test ''type'' is not specified correctly. Please check input. returning.'); +end + + +%% sorting/saving results into output variable +results.p = WX.p; +results.r = WX.r; +results.stats = WX.stats; + +end diff --git a/sami_toolbox/+sami/+stat/bonf.m b/sami_toolbox/+sami/+stat/bonf.m new file mode 100644 index 0000000..2ffcaac --- /dev/null +++ b/sami_toolbox/+sami/+stat/bonf.m @@ -0,0 +1,14 @@ +function pCorr = bonf(p) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +nTests = size(p,1); +pCorr = max(p .* nTests,1); + +end + diff --git a/sami_toolbox/+sami/+stat/bootConfint.m b/sami_toolbox/+sami/+stat/bootConfint.m new file mode 100644 index 0000000..71b1d5a --- /dev/null +++ b/sami_toolbox/+sami/+stat/bootConfint.m @@ -0,0 +1,172 @@ +function [ci_lo, ci_up, p] = bootConfint(x, B, alternative, user_opts) +% bootConfint +% Compute confidence intervals based on bootstrap distribution using +% various methods and use them to derive p-values for hypothesis testing. +% +% [ci_lo, ci_up, p] = bootConfint(x, B, alternative, user_opts) +% +% Currently two methods are available, "percentile" and "bca". P-values +% can be computed for both one-tailed and two-tailed testing. +% +% INPUT +% ----- +% x : vector, original sample of length #subjects +% B : vector, bootstrap distribution of the statistic of +% interest of length #bootstrap_iterations +% alternative : string, either 'two-tailed' or 'greater' +% user_opts : struct, user options used by RSA Toolbox. Options are ... +% user_opts.bootCImethod: +% - "percentile": bootstrap percentile method (default) +% - "bca": bias-corrected and accelerated CI +% user_opts.bootCIalpha: type 1 error probability to use +% for CI definition, default is .05. +% +% OUTPUT +% ------ +% ci_lo, ci_up: lower and upper bounds of CI, respectively +% p : p value of hypothesis test +% +% Literature +% ---------- +% "An Introduction to the Bootstrap" (Efron & Tibshirani, 1993) +% +% 2018-01-27 michael.bannert@tuebingen.mpg.de + +user_opts = sami.util.setIfUnset(user_opts, 'bootCImethod', 'percentile'); +user_opts = sami.util.setIfUnset(user_opts, 'bootCIalpha', .05); + +method = user_opts.bootCImethod; +ci_alpha = user_opts.bootCIalpha; + +n_B = length(B); +n_samples = length(x); +B_sorted = sort(B); + +switch lower(method) + case 'bca' + % Bias-correction, eq. 14.14 + z_hat0 = norminv(sum(B < mean(x)) / n_B); + + % Compute acceleration using jackknife + theta_i = nan(n_samples, 1); + for n = 1:n_samples + x_i = x(setdiff(1:n_samples, n)); + theta_i(n) = mean(x_i); + end + theta_pt = mean(theta_i); + + % Eq. 14.15 + a_hat = sum((theta_pt * ones(n_samples, 1) - theta_i).^3) / ... + 6 * (sum((theta_pt * ones(n_samples, 1) - theta_i).^2)) ^(3/2); + + z_alpha = norminv(ci_alpha / 2); + z_one_minus_alpha = norminv(1 - ci_alpha / 2); + + % Eq. 14.10 + alpha1 = normcdf( ... + z_hat0 + (z_hat0 + z_alpha) / ... + (1 - a_hat * (z_hat0 + z_alpha)) ... + ); + + alpha2 = normcdf( ... + z_hat0 + (z_hat0 + z_one_minus_alpha) / ... + (1 - a_hat * (z_hat0 + z_one_minus_alpha)) ... + ); + + % Find indices of bootstrap distribution corresponding to desired + % alpha + alpha1_n = max(floor(alpha1 * n_B), 1); % At least 1 + alpha2_n = ceil(alpha2 * n_B); + + % Eq. 14.9 + ci_lo = B_sorted(alpha1_n); + ci_up = B_sorted(alpha2_n); + + % Calculate p-value by solving eq. 14.10 for alpha. P-value should + % be one-tailed for the null hypothesis that the true value is <= 0. + % Find the z_alpha for which the threshold 0 equals the lower CI + % bound. + + % 1. Precompute inverse cumulative Gaussian of lower confidence + % bound that would just include 0. (Used several times in equation + % below.) + + % Index must be at least 1. Increase number of bootstrap iterations + % for more precise p-values + alpha0_n = find(B_sorted > 0, 1) - 1; + alpha1_norminv = norminv(alpha0_n / n_B); + + % 2. Find CI alpha for which the lower bound would just enclose + % 0. + % This is eq. 14.10 solved for z_alpha + z_alpha0 = ... + (alpha1_norminv * (1 - a_hat * z_hat0) + ... + z_hat0 * (z_hat0 * a_hat - 2)) / ... + (a_hat * (alpha1_norminv - z_hat0) - 1); + + % Super non-significant, i.e., no positive values in bootstrap + % distribution + if isempty(z_alpha0) + z_alpha0 = -inf; + end + + % Super significant, i.e., only positive values in bootstrap + % distribution. P-value only limited by number of bootstrap + % iterations. Increasing number of iterations would possibly help. + if isnan(z_alpha0) + z_alpha0 = norminv(1 - 1/n_B); + end + + switch lower(alternative) + case 'greater' + % Convert z_alpha to p-value + p = 1 - normcdf(z_alpha0); + + case 'two-tailed' + % For which z_alpha does the upper CI bound equal 0? See + % one-tailed case for explanation + + % Index must be at least n_B - 1. Increase number of + % bootstraps for more precise p-values + alpha0_n_up = min(n_B-1, find(B_sorted < 0, 1, 'last') + 1); + if isempty(alpha0_n_up) + alpha0_n_up = n_B-1; + end + alpha1_norminv_up = norminv(alpha0_n_up / n_B); + + z_alpha0_up = ... + (alpha1_norminv_up * (1 - a_hat * z_hat0) + ... + z_hat0 * (z_hat0 * a_hat - 2)) / ... + (a_hat * (alpha1_norminv_up - z_hat0) - 1); + + % Compare the two one-tailed p-values and return the + % smaller one multiplied by two. + p = min(1 - normcdf(z_alpha0), normcdf(z_alpha0_up)) * 2; + if isempty(p) || p > 1 + keyboard; + end + end + + case 'percentile' + % Efron & Tibshirani (1993, p. 170) + ci_lo = B_sorted(max(1, floor((ci_alpha / 2) * n_B))); + ci_up = B_sorted(ceil((1 - ci_alpha / 2) * n_B)); + + switch lower(alternative) + case 'greater' + % One-tailed p-value for null hypothesis that true value is + % below or equal to 0. + p = max(sum(B < 0) / n_B, 1/n_B); + + case 'two-tailed' + % Two-tailed p-value for null hypothesis that true value + % lies within confidence interval + p = max((min(sum(B < 0), sum(B >= 0)) / n_B) * 2, 1/n_B); + end + + otherwise + error(['Valid method ''%s'' to compute bootstrap CIs\nValid ' ... + 'methods are: ''bca'', ''percentile'''], method); +end + +assert(p <= 1 && p >= 0, '''p'' is not a probability.'); diff --git a/sami_toolbox/+sami/+stat/bootstrapRDMs.m b/sami_toolbox/+sami/+stat/bootstrapRDMs.m new file mode 100644 index 0000000..7450ec6 --- /dev/null +++ b/sami_toolbox/+sami/+stat/bootstrapRDMs.m @@ -0,0 +1,72 @@ +function bootstrapRs = bootstrapRDMs(bootstrappableRefRDMs, featRDMs, userOptions) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +% Sort out defaults +userOptions = sami.util.setIfUnset(userOptions, 'rdms_nBootstrap', 1000); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_pairWiseCorr', 'Kendall_taua'); + +% Constants +nBS = userOptions.rdms_nBootstrap; +nConditions = size(bootstrappableRefRDMs, 1); +nSubjects = size(bootstrappableRefRDMs, 3); +nFeatRDMs = size(featRDMs, 3); + +if ~(size(bootstrappableRefRDMs, 1) == size(featRDMs, 1)) + error('bootstrapRDMComparison:DifferentSizedRDMs', 'Two RDMs being compared are of different sizes. This is incompatible with bootstrap methods!'); +end + +fprintf([' -> resampling ''conditions'' ' num2str(nBS) ' times: 0%%']); +tic; + +% Come up with the random samples (with replacement) +resampledConditionIs = ceil(nConditions * rand(nBS, nConditions)); + +% Preallocation +bootstrapRs = nan(nFeatRDMs, nBS); + +% replace the diagonals for each instance of the candidate RDMs with NaN entries +for subI = 1:nSubjects + temp = bootstrappableRefRDMs(:,:,subI); + temp(logical(eye(size(temp,1)))) = nan; + bootstrappableRefRDMs(:,:,subI) = temp; +end + +% Need to create one candidate RDM for each reference because we don't want +% to correlate averaged RDMs. +candRDMs3rd = cell(nFeatRDMs, 1); +for candRDMI = 1:nFeatRDMs + candRDMs3rd{candRDMI} = repmat(featRDMs(:,:,candRDMI), ... + [1, 1, nSubjects]); +end + +% Bootstrap +for b = 1:nBS + for candRDMI = 1:nFeatRDMs + localReferenceRDMs = bootstrappableRefRDMs(resampledConditionIs(b,:),resampledConditionIs(b,:),:); + localTestRDM = candRDMs3rd{candRDMI}(resampledConditionIs(b,:), resampledConditionIs(b,:), :); + + if isequal(userOptions.rdms_pairWiseCorr,'Kendall_taua') + bootstrapRs(candRDMI, b) = mean(diag(sami.stat.rankCorr_Kendall_taua(squeeze(sami.rdm.vectorizeRDMs(localReferenceRDMs)),squeeze(sami.rdm.vectorizeRDMs(localTestRDM))))); + else + c1 = squeeze(sami.rdm.vectorizeRDMs(localReferenceRDMs)); + c2 = squeeze(sami.rdm.vectorizeRDMs(localTestRDM)); + tmp = corr(c1, c2, 'type',userOptions.rdms_pairWiseCorr,'rows','pairwise'); + bootstrapRs(candRDMI, b) = mean(diag( tmp )); + end + end + + if mod(b,nBS/100) == 0 + perc = round(b / nBS * 100); + fprintf([repmat('\b',1,numel(num2str(perc-1)) + 1) '%d%%'],perc); + end +end + +t = toc; +fprintf([' ... DONE [in ' num2str(ceil(t)) 's]\n']); +end diff --git a/sami_toolbox/+sami/+stat/holm.m b/sami_toolbox/+sami/+stat/holm.m new file mode 100644 index 0000000..2701ecb --- /dev/null +++ b/sami_toolbox/+sami/+stat/holm.m @@ -0,0 +1,33 @@ +function pCorr = holm(p) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +[numtests, nS] = size(p); +pCorr = nan(size(p)); + +% loop samples +for i = 1:nS + % next we will calculate the p value using the homls method. + [psorted, idx_p] = sort(p(:,i)); + + p_holm_sorted = zeros(numtests,1); + p_holm_sorted(1,1) = numtests*psorted(1,1); + for j=2:numtests + p_holm_sorted(j,1) = (numtests-j+1)*psorted(j,1); +% p_holm_sorted(j,1) = max(p_holm_sorted(j-1,1), (numtests-j+1)*psorted(j,1)); + end + + p_holm_sorted(:,1) = min(p_holm_sorted(:,1),1); + for j=1:numtests + pCorr(j,i) = p_holm_sorted(idx_p==j,1); + end +end + +% p_holm_sorted(idx_p) + +end diff --git a/sami_toolbox/+sami/+stat/iterativeUpperBound.m b/sami_toolbox/+sami/+stat/iterativeUpperBound.m new file mode 100644 index 0000000..c1980d4 --- /dev/null +++ b/sami_toolbox/+sami/+stat/iterativeUpperBound.m @@ -0,0 +1,135 @@ +function ceiling_upperBound = iterativeUpperBound(refRDMs, bestFitRDM, meanRDM_avgCorr, monitor) +% Iterative approach to optimise the estimate by searching an RDM that yields a higher +% average tau a correlation with the single-subject RDMs. +% +% Iterations stops and returns best estimate of the upper bound, when A) convergence criteria +% are fulfilled, or B) maxIter are reached +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + + +%% preparations +if ~exist('monitor','var'), monitor = false; end +[~, nDiss, nSubj] = size(refRDMs); + +%% find best-fit RDM by gradient descent to estimate upper bound for ceiling +fprintf(' -> iterative optimisation of the upper bound... of max. iterations: 0%%'); +tic; +nExploredDirections = 50; +maxIter = 100; + +stepSize = 1/nSubj/10; +bestFitRDMcorrHistory = nan(maxIter,1); + +if monitor + propOfSuccessSamplesHistory = nan(maxIter,1); + propOfStationSamplesHistory = nan(maxIter,1); + stepSizeHistory = nan(maxIter,1); +end + +convergenceZoneLength = 10; +convergenceStdThreshold = 0.000001; +converged = false; + +nImproved = 0; +iter = 0; + +while true + iter = iter + 1; + if mod(iter,maxIter/100) == 0 + perc = round(iter / maxIter * 100); + fprintf([repmat('\b',1,numel(num2str(perc-1)) + 1) '%d%%'], perc); + end + + % create RANDOM candidate RDMs + candRDMs = repmat(bestFitRDM,[1 1 nExploredDirections])+stepSize*randn(1,nDiss,nExploredDirections); + candRDMs = cat(3,candRDMs,bestFitRDM); % the last one is the current best-fit RDM + + % calculate correlations for each candidate RDM + RDMcorrs = nan(nExploredDirections,nSubj); + for iSubj = 1:nSubj + for directionI = 1:nExploredDirections+1 + RDMcorrs(directionI,iSubj) = sami.stat.rankCorr_Kendall_taua(candRDMs(1,:,directionI),refRDMs(1,:,iSubj)); + end + end + currBestFitRDMcorr = mean(RDMcorrs(end,:),2); + candRDMcorrs = mean(RDMcorrs(1:end-1,:),2); + bestFitRDMcorrHistory(iter) = currBestFitRDMcorr; + + % test for convergence + if iter >= convergenceZoneLength && std(bestFitRDMcorrHistory(iter-convergenceZoneLength+1:iter)) < convergenceStdThreshold + converged = true; + end + + % visualise + if monitor && (mod(iter,2) == 0 || converged || iter >= maxIter) + h=figure(400); set(h,'Color','w'); clf; + + subplot(3,1,1); hold on; + plot([1 iter],[meanRDM_avgCorr meanRDM_avgCorr],'-k','LineWidth',5); + plot(bestFitRDMcorrHistory(1:iter),'r','LineWidth',2); + xlabel('time step'); ylabel('current best-fit RDM''s average corr. to ref. RDM estimates','LineWidth',2); + legend('mean-RDM avg. corr.','optimised-RDM avg. corr.','Location','SouthEast'); + title(['current best-fit RDM''s average Kendall_taua correlation to the RDM estimates: ',num2str(currBestFitRDMcorr)]); + + subplot(3,2,3); plot(stepSizeHistory(1:iter),'LineWidth',3); xlabel('iteration'); ylabel('step size'); + subplot(3,2,4); hold on; + plot(propOfSuccessSamplesHistory(1:iter),'r','LineWidth',4); xlabel('iteration'); + plot(propOfStationSamplesHistory(1:iter),'b','LineWidth',2); + xlabel('iteration'); ylabel('proportions of samples'); + legend('successful','stationary','Location','NorthEast'); + + subplot(3,2,5); hold on; + image(squareform(bestFitRDM),'CDataMapping','scaled','AlphaData',~isnan(squareform(bestFitRDM))); + colormap(gca, sami.fig.createColormap('RDMs')); + set(gca,'XTick',[],'YTick',[]); + title('current best-fit RDM'); + axis square off; + drawnow; + end + + % break if converged or time out or subject pressed 'p' for 'proceed' + if converged + t = toc; + fprintf('\n -> converged + upper bound improved %d times... DONE [in %ds]\n', nImproved, ceil(t)) + break; + end + + if iter >= maxIter + t = toc; + warning([' >>> WARNING <<< noiseCeilingOfAvgRDMcorr: iterative ceiling estimation procedure did not converge [within ' num2str(t) 's].']); + break; + end + + % check if any of the explored direction led to an improvement + avgRDMcorrDiffs = candRDMcorrs - currBestFitRDMcorr; + + [mxImprovement, mxImprovementI] = max(avgRDMcorrDiffs); + if mxImprovement > 0 + bestFitRDM = candRDMs(1,:,mxImprovementI); + nImproved = nImproved + 1; + end + + if monitor + propOfSuccessSamplesHistory(iter) = sum(avgRDMcorrDiffs > 0) / nExploredDirections; + propOfStationSamplesHistory(iter) = sum(avgRDMcorrDiffs == 0) / nExploredDirections; + stepSizeHistory(iter) = stepSize; + end + + if 0.3 < sum(avgRDMcorrDiffs > 0) / nExploredDirections || 0.1 < sum(avgRDMcorrDiffs == 0) / nExploredDirections + % more than 30% of our samples improve the fit or more than 10% of our samples are + % on a local plateau, either way: take bigger steps + stepSize = stepSize*(2^(1/2)); % move faster + elseif sum(avgRDMcorrDiffs > 0) / nExploredDirections < 0.1 + % less than 10% of the explored directions improve the fit (and we're not on a plateau). + % might be close to the peak, where every way is down, so look more locally. + stepSize = stepSize/(2^(1/2)); % move slower + end +end % while true + +% define ceiling upper bound +ceiling_upperBound = currBestFitRDMcorr; + +end + diff --git a/sami_toolbox/+sami/+stat/multipleTesting.m b/sami_toolbox/+sami/+stat/multipleTesting.m new file mode 100644 index 0000000..b37a885 --- /dev/null +++ b/sami_toolbox/+sami/+stat/multipleTesting.m @@ -0,0 +1,35 @@ +function pCorr = multipleTesting(input,mtMethod,alpha) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +if ~exist('alpha','var') || isempty(alpha), alpha = 0.05; end + +pCorr = nan(size(input)); + +% loop columns +[~, nCols] = size(input); +for c = 1:nCols + p = input(:,c); + + switch mtMethod + case 'none' + pCorr(:,c) = p; + case 'bonf' + pCorr(:,c) = sami.stat.bonf(p); + case 'holm' + pCorr(:,c) = sami.stat.holm(p); + case 'FDR' + pCorr(:,c) = sami.stat.FDR(p,alpha); + + otherwise + pCorr = input; + warning('multipleTesting:unknownMethod', 'Please set multiple testing method to "bonf", "holm", "FDR" or "none". Now, "none" correction will be done!!!'); + return; + end +end +end diff --git a/sami_toolbox/+sami/+stat/noiseCeilingOfAvgRDMcorr.m b/sami_toolbox/+sami/+stat/noiseCeilingOfAvgRDMcorr.m new file mode 100644 index 0000000..ef0a0eb --- /dev/null +++ b/sami_toolbox/+sami/+stat/noiseCeilingOfAvgRDMcorr.m @@ -0,0 +1,88 @@ +function [ceiling_upperBound, ceiling_lowerBound, bestFitRDM] = noiseCeilingOfAvgRDMcorr(refRDMs, corrType) +% This function estimates the lower and upper bounds of the noise ceilling, i.e. the +% heighest average RDM correlation, across a given set of reference RDMs [refRDMs], a true +% model's RDM prediction could achieve, depending on the variability of the reference RDMs. +% +% The lower bound is estimated via a leave-one-subject-out crossvalidation procedure. +% +% The upper bound is estimated, depending on the type of correlation, by finding the +% hypothetical model RDM, that maximises the average correlation to the reference RDMs. +% +% For further information and a more detailed explanation please see: +% Nili H, et al., (2014) A Toolbox for Representational Similarity Analysis. +% PLoS Comput Biol 10(4): e1003553. https://doi.org/10.1371/journal.pcbi.1003553 +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + + +%% preparations +if ~exist('corrType','var'), corrType = 'Kendall_taua'; end +monitor_itUpperBound = false; + +disp(' -> Estimating the ceiling for the average RDM correlation...'); + +refRDMs = sami.rdm.vectorizeRDMs(sami.rdm.unwrapRDMs(refRDMs)); +[~, nDiss, nSubj] = size(refRDMs); + +%% compute the upper bound for the ceiling +switch corrType + case 'Pearson' + % Pearson correlation distance = const * Euclidean dist. squared in z-transformed + % RDM space. Thus, z-transform RDMs to make the mean RDM minimise the average + % correlation distance to the single-subject RDMs. + refRDMs = refRDMs - repmat(mean(refRDMs,2),[1 nDiss 1]); + refRDMs = refRDMs ./ repmat(std(refRDMs,[],2),[1 nDiss 1]); + bestFitRDM = mean(refRDMs,3); + + case 'Spearman' + % Spearman correlation distance = const * Euclidean dist squares in rank-transformed + % RDM space. Thus, rank-transform RDMs to make the mean RDM minimise the average + % correlation distance to the single-subject RDMs. + refRDMs = reshape( tiedrank( reshape(refRDMs,[nDiss nSubj]) ) ,[1 nDiss nSubj]); + bestFitRDM = mean(refRDMs,3); + + case 'Kendall_taua' + % No closed-form solution providing a tight upper bound for the ceiling (to our + % knowledge), so initialise with the mean of the rank-transformed RDMs and + % optimise iteratively. + refRDMs = reshape( tiedrank( reshape(refRDMs,[nDiss nSubj]) ) ,[1 nDiss nSubj]); + bestFitRDM = mean(refRDMs,3); +end + +%% estimate lower bound on the ceiling +LOSOcorrs = nan(nSubj,1); +for iSubj = 1:nSubj + currSubjRDM = refRDMs(:,:,iSubj); + + LOSOrefRDMs = refRDMs; + LOSOrefRDMs(:,:,iSubj) = []; + avgLOSORefRDM = nanmean(LOSOrefRDMs,3); + + LOSOcorrs(iSubj) = correlation(currSubjRDM,avgLOSORefRDM,corrType); +end +ceiling_lowerBound = mean(LOSOcorrs); + +%% estimate the upper bound on the ceiling + avgRDM_corrs = nan(nSubj,1); + for iSubj = 1:nSubj + avgRDM_corrs(iSubj) = correlation(bestFitRDM,refRDMs(1,:,iSubj),corrType); + end +if isequal(corrType,'Pearson') || isequal(corrType,'Spearman') + ceiling_upperBound = mean(avgRDM_corrs); + +elseif isequal(corrType,'Kendall_taua') + ceiling_upperBound = sami.stat.iterativeUpperBound(refRDMs, bestFitRDM, mean(avgRDM_corrs), monitor_itUpperBound); +end + +end % noise_ceiling_function + +%% compute correlations of all types +function r = correlation(a,b,corrType) + switch corrType + case 'Kendall_taua' + r = sami.stat.rankCorr_Kendall_taua(a(:),b(:)); + otherwise + r = corr(a(:),b(:),'type',corrType); + end +end \ No newline at end of file diff --git a/sami_toolbox/+sami/+stat/rankCorr_Kendall_taua.m b/sami_toolbox/+sami/+stat/rankCorr_Kendall_taua.m new file mode 100644 index 0000000..a87cc4d --- /dev/null +++ b/sami_toolbox/+sami/+stat/rankCorr_Kendall_taua.m @@ -0,0 +1,41 @@ +function [taua,pVal] = rankCorr_Kendall_taua(a,b) +% computes the Kendall's tau a correlation coefficient between the input vectors (a and b). +% NaN entries would be removed from both. +% Also, p-Value is calculated to check for significance. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + + +%% preparations +a = a(:); +b = b(:); +validEntryIs = ~isnan(a) & ~isnan(b); +a = a(validEntryIs); +b = b(validEntryIs); +n = size(a,1); + +%% compute Kendall rank correlation coefficient tau-a +K = 0; +for k = 1:n-1 + pairRelations_a = sign(a(k)-a(k+1:n)); + pairRelations_b = sign(b(k)-b(k+1:n)); + K = K + sum(pairRelations_a .* pairRelations_b); +end +taua = K / (n*(n-1)/2 ); % normalise by the total number of pairs + + +%% compute significance (https://www.real-statistics.com/correlation/kendalls-tau-correlation/kendalls-tau-normal-approximation/) +C = nchoosek(n,2); + +se = sqrt( (2*n+5)/ C ) / 3; +Z = taua / se; + +pVal = normcdf(-abs(Z)) * 2; + +%% in case statistics toolbox is missing +% F = @(x)(exp (-0.5*(x.^2))./sqrt (2*pi)); +% p = integral (F, abs(Z), 100); +% pVal = 2*p; + +end diff --git a/sami_toolbox/+sami/+stat/relRank.m b/sami_toolbox/+sami/+stat/relRank.m new file mode 100644 index 0000000..39a1142 --- /dev/null +++ b/sami_toolbox/+sami/+stat/relRank.m @@ -0,0 +1,13 @@ +function p = relRank(null,value) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +null = [null(:); value]; +p = sum( null(:) < value ) / numel(null); + +end diff --git a/sami_toolbox/+sami/+stat/rmANOVA.m b/sami_toolbox/+sami/+stat/rmANOVA.m new file mode 100644 index 0000000..56d7b02 --- /dev/null +++ b/sami_toolbox/+sami/+stat/rmANOVA.m @@ -0,0 +1,100 @@ +function results = rmANOVA(d,CorrMethod,alpha) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +% prepare rmANOVA +design = ''; +for i = 1:size(d,2) + design = [design 'd' num2str(i) ',']; %#ok +end +design = [design(1:end-1) '~1']; +condT = table((1:size(d,2))','VariableNames',{'Conditions'}); + +% fit the repeated measures model +rmModel = fitrm(array2table(d),design,'WithinDesign',condT); + +% get results for rmANOVA + post-hoc pairwise comparisons +[aR,~,cMat] = ranova(rmModel); +posthoc = multcompare(rmModel,'Conditions','ComparisonType',CorrMethod,'alpha',alpha); + +% calculate partial eta^2 +aR.pEtaSq(1) = aR.SumSq(1)/sum(aR.SumSq(1:2)); + +% calculat mauchly test + epsilon and add to anova table +res_mauchly = mauchly(rmModel,cMat); +res_epsilon = epsilon(rmModel,cMat); + +aR.mauchlyChiSq(1) = res_mauchly.ChiStat; +aR.mauchlyDF(1) = res_mauchly.DF; +aR.mauchlyPvalue(1) = res_mauchly.pValue; +aR.epsilonGG = reshape(repmat(res_epsilon.GreenhouseGeisser',[2 1]),numel(res_epsilon.GreenhouseGeisser)*2,1); +aR.epsilonHF = reshape(repmat(res_epsilon.HuynhFeldt',[2 1]),numel(res_epsilon.HuynhFeldt)*2,1); + +% adjust degrees of freedom +aR.dfGG = aR.DF .* aR.epsilonGG; +aR.dfHF = aR.DF .* aR.epsilonHF; + +% decide if, and which, correction will be used, according to Girden (1992) +if aR.mauchlyPvalue(1) > 0.05 + aR.correction{1} = 'none'; +elseif aR.epsilonGG(1) > .75 + aR.correction{1} = 'HF'; +else + aR.correction{1} = 'GG'; +end + +% plot horizontal sign. lines if sign. pairwise comparisons exist +ph_index = diff([posthoc.Conditions_1,posthoc.Conditions_2],[],2) > 0; +sami.fig.drawSigLines([posthoc.Conditions_1(ph_index),posthoc.Conditions_2(ph_index)],posthoc.pValue(ph_index)); + +% sorting rmANOVA results, depending on sphericity +results.test = 'rmANOVA'; +switch aR.correction{1} + case 'none' + results.info = sprintf('no correction <- Mauchly: chi^2(%d) = %.2f, %s',... + aR.mauchlyDF(1),... + aR.mauchlyChiSq(1),... + sami.util.getPString(aR.mauchlyPvalue(1))); + results.S = 'F'; + results.V = aR.F(1); + results.df1 = aR.DF(1); + results.df2 = aR.DF(2); + results.p = aR.pValue(1); + + case 'GG' + results.info = sprintf('Greenhouse-Geisser correction <- Mauchly: chi^2(%d) = %.2f, %s, epsilon = %.3f',... + aR.mauchlyDF(1),... + aR.mauchlyChiSq(1),... + sami.util.getPString(aR.mauchlyPvalue(1)),... + aR.epsilonGG(1)); + results.S = 'F'; + results.V = aR.F(1); + results.df1 = aR.dfGG(1); + results.df2 = aR.dfGG(2); + results.p = aR.pValueGG(1); + + case 'HF' + results.info = sprintf('Huyn-Feldt correction <- Mauchly: chi^2(%d) = %.2f, %s, epsilon = %.3f',... + aR.mauchlyDF(1),... + aR.mauchlyChiSq(1),... + sami.util.getPString(aR.mauchlyPvalue(1)),... + aR.epsilonGG(1)); + results.S = 'F'; + results.V = aR.F(1); + results.df1 = aR.dfHF(1); + results.df2 = aR.dfHF(2); + results.p = aR.pValueHF(1); + +end +results.pEtaSq = aR.pEtaSq(1); + +% store anova- and multCompare-tables +results.tbl = aR; +results.multComp = posthoc; + +end diff --git a/sami_toolbox/+sami/+stat/signrank_onesided.m b/sami_toolbox/+sami/+stat/signrank_onesided.m new file mode 100644 index 0000000..d87d01a --- /dev/null +++ b/sami_toolbox/+sami/+stat/signrank_onesided.m @@ -0,0 +1,21 @@ +function p = signrank_onesided(x) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +p_signrank = signrank(x,[],'alpha',0.05,'method','exact'); + +% correct p value for testing values > 0 +if median(x) > 0 + p_signrank = p_signrank/2; +else + p_signrank = 1-p_signrank/2; +end + +p = p_signrank; + +end%function diff --git a/sami_toolbox/+sami/+stat/tTest.m b/sami_toolbox/+sami/+stat/tTest.m new file mode 100644 index 0000000..c254e1a --- /dev/null +++ b/sami_toolbox/+sami/+stat/tTest.m @@ -0,0 +1,78 @@ +function results = tTest(d,type,m) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 03/2021 + +%% define default behavior +if ~exist('m','var') || isempty(m), m = 0; end + +%% perform appropriate t-test, as specified by 'type' +switch type + case 'one_sample' + % naming t-Test results + results.test = 'one-sample t-Test'; + + % t-Test + [~,tt.p,tt.ci,tt.stats] = ttest(d,m); + + % bonferonni correction of p-values + if size(d,2) > 1 + results.info = ['p values bonferonni corrected by ' num2str(size(d,2)) ' indepedent tests']; + tt.p = tt.p .* size(d,2); + end + + % calculate cohen's d + cd = abs( mean(d)-m ) ./ std(d); + + % plot asterisks if sign. difference from 'm' + sami.fig.drawSigAsterisks(1:size(d,2),tt.p); + + % draw line, if tested against ~=0 + sami.fig.drawRefenceTestLine(m); + + case 'two_sample' + % naming t-Test results + results.test = 'two-sample t-Test'; + + % t-Test + [~,tt.p,tt.ci,tt.stats] = ttest2(d(:,1),d(:,2)); + + % calculate cohen's d + SDp = sqrt((std(d(:,1))^2+std(d(:,2))^2)/2); + cd = abs(diff(mean(d))) / SDp; + + % plot horizontal sign. lines if sign. pairwise comparisons exist + sami.fig.drawSigLines([1 2],tt.p); + + case 'paired' + % naming t-Test results + results.test = 'paired t-Test'; + + % t-Test + [~,tt.p,tt.ci,tt.stats] = ttest(d(:,1),d(:,2)); + + % calculate cohen's d + D = d(:,1) - d(:,2); + cd = mean(D) / std(D); + + % plot horizontal sign. lines if sign. pairwise comparisons exist + sami.fig.drawSigLines([1 2],tt.p) + + otherwise + error('*** sami:ERROR *** t-test ''type'' is not specified correctly. Please check input. returning.'); +end + +%% sorting/saving results into output variable +results.S = 't'; +results.V = tt.stats.tstat; +results.df = tt.stats.df; +results.p = tt.p; +results.d = cd; +results.ci = tt.ci; +results.stats = tt.stats; + +end diff --git a/sami_toolbox/+sami/+util/deblank.m b/sami_toolbox/+sami/+util/deblank.m new file mode 100644 index 0000000..0e47581 --- /dev/null +++ b/sami_toolbox/+sami/+util/deblank.m @@ -0,0 +1,20 @@ +function txt = deblank(txt) +% gets an input string and replaces "spaces" (' ') with underscores ('_'). +% This avoids the problems when saving files, this way with displaying the names in which +% the first character after the underscore would be displayed as an index. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +% replace spaces +if iscell(txt) + for i = 1:numel(txt) + line = txt{i}; + line(line==32) = '_'; + txt{i} = line; + end +else + txt(txt==32)='_'; +end + +end diff --git a/sami_toolbox/+sami/+util/deunderscore.m b/sami_toolbox/+sami/+util/deunderscore.m new file mode 100644 index 0000000..eef5135 --- /dev/null +++ b/sami_toolbox/+sami/+util/deunderscore.m @@ -0,0 +1,20 @@ +function txt = deunderscore(txt) +% gets an input string and adds to an underscores ('_') a backslash +% ('\_').This avoids the problems with displaying the names in which the +% first character after the underscore would be displayed as an index. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +% replace underscores +if iscell(txt) + for i = 1:numel(txt) + line = txt{i}; + line = strrep(line,'_','\_'); + txt{i} = line; + end +else + txt = strrep(txt,'_','\_'); +end + +end diff --git a/sami_toolbox/+sami/+util/getFigI.m b/sami_toolbox/+sami/+util/getFigI.m new file mode 100644 index 0000000..b330289 --- /dev/null +++ b/sami_toolbox/+sami/+util/getFigI.m @@ -0,0 +1,32 @@ +function figI = getFigI(varargin) +% find figure number (or numbers), which is (are) not in use by any open figure +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% init +if isempty(varargin) + n = 1; +else + n = varargin{1}; +end + +%% get open figs +figHandles = get(groot, 'Children'); +if ~isempty(figHandles) + inUse = [figHandles(:).Number]; +else + inUse = []; +end + +%% get fist n possible figNumbers +figI = nan(1,n); +for f = 1:n + allUsed = [inUse figI(~isnan(figI))]; + + figI(f) = 1; + while ~isempty(allUsed) && any(ismember(allUsed,figI(f))) + figI(f) = figI(f) + 1; + end +end +end diff --git a/sami_toolbox/+sami/+util/getPString.m b/sami_toolbox/+sami/+util/getPString.m new file mode 100644 index 0000000..7d0bee6 --- /dev/null +++ b/sami_toolbox/+sami/+util/getPString.m @@ -0,0 +1,30 @@ +function [pStr, stars] = getPString(p,depth) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +if ~exist('depth','var') || isempty(depth), depth = 3; end + +if p < 0.0001 && depth >= 4 + pStr = 'p < .0001'; + stars = '****'; +elseif p < 0.001 && depth >= 3 + pStr = 'p < .001'; + stars = '***'; +elseif p < 0.01 && depth >= 2 + pStr = 'p < .01'; + stars = '**'; +elseif p < 0.05 + pStr = 'p < .05'; + stars = '*'; +else + pStr = sprintf('p = %0.2f', p); + stars = 'n.s.'; +end + +end + diff --git a/sami_toolbox/+sami/+util/getStars.m b/sami_toolbox/+sami/+util/getStars.m new file mode 100644 index 0000000..45fe83f --- /dev/null +++ b/sami_toolbox/+sami/+util/getStars.m @@ -0,0 +1,25 @@ +function stars = getStars(p,depth) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +if ~exist('depth','var') || isempty(depth), depth = 3; end + +stars = ''; + +if p < 0.0001 && depth >= 4 + stars = '****'; +elseif p < 0.001 && depth >= 3 + stars = '***'; +elseif p < 0.01 && depth >= 2 + stars = '**'; +elseif p < 0.05 + stars = '*'; +end + +end + diff --git a/sami_toolbox/+sami/+util/getStimColors.m b/sami_toolbox/+sami/+util/getStimColors.m new file mode 100644 index 0000000..b0c3dfb --- /dev/null +++ b/sami_toolbox/+sami/+util/getStimColors.m @@ -0,0 +1,30 @@ +function stimColors = getStimColors(category,userOptions) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +stimOrder = sami.util.getStimOrder(userOptions); + +if ismember(category , userOptions.stimuli_settings(1,2:end)) + [~, idx_Set] = ismember(category, userOptions.stimuli_settings(1,2:end)); + [~, idx_Key] = ismember(category, {userOptions.stimuli_naming_key.name}); +else + idx_Set = 1; + [~, idx_Key] = ismember(userOptions.stimuli_settings(1,2), {userOptions.stimuli_naming_key.name}); +end + +stimCategories = cell2mat(userOptions.stimuli_settings(2:end,idx_Set + 1)); +categoryColors = userOptions.stimuli_naming_key(idx_Key).color; + +stimColors = nan(size(stimOrder,1),3); +for i = 1:size(stimColors,1) + thisCat = stimCategories(stimOrder(i)); + stimColors(i,:) = categoryColors{thisCat}; +end + +end + diff --git a/sami_toolbox/+sami/+util/getStimOrder.m b/sami_toolbox/+sami/+util/getStimOrder.m new file mode 100644 index 0000000..a35d937 --- /dev/null +++ b/sami_toolbox/+sami/+util/getStimOrder.m @@ -0,0 +1,48 @@ +function [StimOrder,sortedFileNames] = getStimOrder(userOptions) +% [StimOrder,sortedFileNames] = getStimOrder(userOptions) +% +% Sorts stimuli with respect to a specified stimulus-category and returns order, as well +% as a sorted cell array containing filenames. +% +% input: +% - userOptions.stimuli_sorting: +% string. 'name_of_category' which is used to sort the stimuli according to. Must +% be present in the first row of 'userOptions.stimuli_settings'. +% Defaults to 1 +% +% - userOptions.stimuli_settings: +% cell array. Each stimulus (column 1) is described by various categories +% (column 2-...). Each category-column starts with 'name_of_category', followed +% by numbers coding for category-groups. +% +% output: +% - StimOrder: +% array. Order of stimuli in first column in 'userOptions.stimuli_settings', +% grouped according to category set in 'userOptions.stimuli_sorting' +% +% - sortedFileNames: +% cell array. Filenames of stimuli sorted as 'StimOrder'. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 10/2020 + +if isfield(userOptions,'stimuli_sorting') && isfield(userOptions,'stimuli_settings') + % sort according to specified condition in userOptions.stimuli_sorting + if ismember(userOptions.stimuli_sorting , userOptions.stimuli_settings(1,2:end)) + [~, sort_cat] = ismember(userOptions.stimuli_sorting, userOptions.stimuli_settings(1,2:end)); + else + sort_cat = 1; + end + + [~, StimOrder] = sort(cell2mat(userOptions.stimuli_settings(2:end,sort_cat + 1))); + sortedFileNames = userOptions.stimuli_settings(StimOrder+1,1); +else + % in case user does not specified stimuli_settings OR stimuli_sorting, stimuli will + % be sorted according to 'dir' function (this is usefull when someone wants to + % analyse behavioral data only... + files = dir(fullfile(userOptions.c3dPath,'*.c3d')); + StimOrder = 1:numel(files); + sortedFileNames = {files(StimOrder).name}'; +end + + diff --git a/sami_toolbox/+sami/+util/gotoDir.m b/sami_toolbox/+sami/+util/gotoDir.m new file mode 100644 index 0000000..cfaf991 --- /dev/null +++ b/sami_toolbox/+sami/+util/gotoDir.m @@ -0,0 +1,35 @@ +function gotoDir(path) +% gotoDir(path) +% Goes to path, making all required directories on the way. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +sIndices = strfind(path, filesep); + +% test all directories within path +for i = 1:numel(sIndices) + if i == 1 && sIndices(i) == 1 + continue; + end + + try + cd(path(1:sIndices(i)-1)); + catch +% fprintf(['The directory "' path(1:sIndices(i)-1) '" doesn''t exist... making it.\n']); + mkdir(path(1:sIndices(i)-1)); + cd(path(1:sIndices(i)-1)); + end +end + +% cleanup final directory! +try + cd(path); +catch +% fprintf(['The directory "' path '" doesn''t exist... making it.\n']); + mkdir(path); +end + +cd(path); + +end diff --git a/sami_toolbox/+sami/+util/selectPlot.m b/sami_toolbox/+sami/+util/selectPlot.m new file mode 100644 index 0000000..0dd73e8 --- /dev/null +++ b/sami_toolbox/+sami/+util/selectPlot.m @@ -0,0 +1,43 @@ +function [hf,ha,figI] = selectPlot(figI,visibleFig) +% - activates or creates figure figI(1) +% - activates or creates plot [figI(2:4)], +% or [1 1 1] if figI(2:4) are missing +% - sets background color to white +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +if ~exist('figI','var'); figI = 0; end +if ~exist('visibleFig','var'); visibleFig = true; end + +% default values +background_color = 'w'; + +% switch visivility from BOOL to on/off +if visibleFig, visibleFig = 'on'; else, visibleFig = 'off'; end + +% handle figure +if figI + if ishghandle(figI(1)) + set(0, 'currentfigure', figI(1)); + hf = gcf; + else + hf = figure(figI(1)); + set(hf,'Visible',visibleFig); + set(hf,'Color',background_color); + end + + if numel(figI) >= 4 + ha = subplot(figI(2),figI(3),figI(4:end)); + else + ha = subplot(1,1,1); + end +else + hf = figure; + set(hf,'Visible',visibleFig); + set(hf,'Color',background_color); + figI(1) = hf; + ha = subplot(1,1,1); +end + +end diff --git a/sami_toolbox/+sami/+util/setIfUnset.m b/sami_toolbox/+sami/+util/setIfUnset.m new file mode 100644 index 0000000..44645b2 --- /dev/null +++ b/sami_toolbox/+sami/+util/setIfUnset.m @@ -0,0 +1,11 @@ +function options = setIfUnset(options,field,value) +% if options.(field) is empty or doesn't exist, this function sets options.(field) to value. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +if ~isfield(options, field) || isempty(options.(field)) + options.(field) = value; +end + +end diff --git a/sami_toolbox/+sami/MDSofRDMs.m b/sami_toolbox/+sami/MDSofRDMs.m new file mode 100644 index 0000000..532a85a --- /dev/null +++ b/sami_toolbox/+sami/MDSofRDMs.m @@ -0,0 +1,83 @@ +function MDSofRDMs(RDMs, userOptions, title_suffix, figI) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% Set defaults and check options struct +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +if ~exist('title_suffix','var') || isempty(title_suffix), title_suffix = ''; end +userOptions = sami.util.setIfUnset(userOptions, 'rdms_pairWiseCorr', 'Spearman'); +userOptions = sami.util.setIfUnset(userOptions, 'MDS_criterion', 'metricstress'); +userOptions = sami.util.setIfUnset(userOptions, 'fig_display', true); +userOptions = sami.util.setIfUnset(userOptions, 'fig_savePDF', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_saveFig', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_savePS', false); + +if numel(RDMs) < 3, warning('RDMsPLotMDS:NotEnoughRDMs', ['Only ' num2str(numel(RDMs)) ' RDMs is not enough. 3 is a minimum for MDS to work; skipping.']); return; end + +%% init +if isempty(title_suffix) + txt1 = ''; txt2 = ''; txt3 = ''; +else + txt1 = '_'; txt2 = '['; txt3 = ']'; +end + +nRDMs = numel(RDMs); + +%% Set MDSOptions +MDSOptions.fig_display = userOptions.fig_display; +MDSOptions.figI = figI; +MDSOptions.fileName = ['2ndOrderMDSofRDMs' txt1 sami.util.deblank(title_suffix)]; +MDSOptions.titleString = ['2nd order MDS of RDMs ' txt2 sami.util.deunderscore(title_suffix) txt3]; + +MDSOptions.MDSCriterion = userOptions.MDS_criterion; +MDSOptions.MDSDistance = ['1 - \it' sami.util.deunderscore(userOptions.rdms_pairWiseCorr) '\rm correaltion']; +MDSOptions.dotLabels = {RDMs(:).name}; +MDSOptions.dotColors = reshape([RDMs(:).color],[3 nRDMs])'; + + +if isfield(userOptions, 'MDS_plotLabels') + MDSOptions.plotLabels = userOptions.MDS_plotLabels; +end +if isfield(userOptions, 'MDS_dotSize') + MDSOptions.dotSize = userOptions.MDS_dotSize; +end +if isfield(userOptions, 'MDS_fontSize') + MDSOptions.fontSize = userOptions.MDS_fontSize; +end + +% prepare legend +if isfield(userOptions, 'MDS_plotLegend') + MDSOptions.plotLegend = userOptions.MDS_plotLegend; + if MDSOptions.plotLegend == 1 + leg.title = []; + leg.Labels = sami.util.deunderscore(MDSOptions.dotLabels); + leg.Colors = MDSOptions.dotColors; + leg.DotIdx = 1:numel(leg.Labels); + MDSOptions.legend = leg; + end +end + +%% calculate RDM of RDMs +distanceMatrix.RDM = 1 - sami.stat.RDMCorrMat(RDMs,userOptions.rdms_pairWiseCorr); +distanceMatrix.name = ['Pairwise RDM correlations: ' title_suffix]; + +%% plotting +disp(['*** Drawing MDS: "' MDSOptions.titleString '" [fig. ' num2str(figI) ']']); +sami.fig.MDSArrangement(distanceMatrix, MDSOptions); + +%% saving +returnHere = pwd; +thisFileName = MDSOptions.fileName; + +% figure +sami.util.gotoDir(fullfile(userOptions.rootPath, 'MDSplots')); +disp([' -> saving MDS_plot FIGURE to ' fullfile(pwd, thisFileName)]); +sami.fig.handleFigure(figI, thisFileName, userOptions); + +cd(returnHere); + diff --git a/sami_toolbox/+sami/MDSofStimuli.m b/sami_toolbox/+sami/MDSofStimuli.m new file mode 100644 index 0000000..59de57b --- /dev/null +++ b/sami_toolbox/+sami/MDSofStimuli.m @@ -0,0 +1,98 @@ +function MDSofStimuli(RDMs, colorCategory, userOptions, figI) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + + +%% Set defaults and check options struct +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(numel(RDMs)); end +if numel(figI) ~= numel(RDMs), error('*** sami:ERROR *** amount of Figure-Numbers and RDMs are not equal.'); end +userOptions = sami.util.setIfUnset(userOptions, 'MDS_criterion', 'metricstress'); +userOptions = sami.util.setIfUnset(userOptions, 'feat_distance', 'Euclidean'); +userOptions = sami.util.setIfUnset(userOptions, 'fig_display', true); +userOptions = sami.util.setIfUnset(userOptions, 'fig_savePDF', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_saveFig', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_savePS', false); + +if userOptions.MDS_plotLabels == 1 + if (~isfield(userOptions, 'stimuli_MDS_labels') || isempty(userOptions.stimuli_MDS_labels)) + disp('*** sami:INFO *** MDS_plotLabels is ON, but there are no stimuli_MDS_Labels provided. Turning it OFF.'); + userOptions.MDS_plotLabels = 0; + elseif size(RDMs(1).RDM,1) ~= numel(userOptions.stimuli_MDS_labels) + error('*** sami:ERROR *** amount of "Stimuli Labels" and "size of RDMs" are not equal.'); + end +end + +%% Set MDSOptions +MDSOptions.fig_display = userOptions.fig_display; +MDSOptions.MDSCriterion = userOptions.MDS_criterion; +if strcmpi(userOptions.feat_distance,'euclidean') + MDSOptions.MDSDistance = ['\it' sami.util.deunderscore(userOptions.feat_distance) '\rm']; +else + MDSOptions.MDSDistance = ['1 - \it' sami.util.deunderscore(userOptions.feat_distance) '\rm correaltion']; +end +MDSOptions.dotColors = sami.util.getStimColors(colorCategory,userOptions); + +if userOptions.MDS_plotLabels == 1 + stimOrder = sami.util.getStimOrder(userOptions); + MDSOptions.dotLabels = userOptions.stimuli_MDS_labels(stimOrder); +end + +if isfield(userOptions, 'MDS_plotLabels') + MDSOptions.plotLabels = userOptions.MDS_plotLabels; +end +if isfield(userOptions, 'MDS_dotSize') + MDSOptions.dotSize = userOptions.MDS_dotSize; +end +if isfield(userOptions, 'MDS_fontSize') + MDSOptions.fontSize = userOptions.MDS_fontSize; +end + +% prepare legend +if isfield(userOptions, 'MDS_plotLegend') + MDSOptions.plotLegend = userOptions.MDS_plotLegend; + if MDSOptions.plotLegend == 1 + [~, leg_idx_Key] = ismember(colorCategory, {userOptions.stimuli_naming_key.name}); + leg.title = userOptions.stimuli_naming_key(leg_idx_Key).name; + leg.Labels = userOptions.stimuli_naming_key(leg_idx_Key).condition; + leg.Colors = userOptions.stimuli_naming_key(leg_idx_Key).color; + for i = 1:numel(leg.Colors) + j = 1; + while any(MDSOptions.dotColors(j,:) ~= leg.Colors{i}) + j = j+1; + end + leg.DotIdx(i) = j; + end + MDSOptions.legend = leg; + end +end + +%% loop RDMs +for iRDM = 1:numel(RDMs) + RDMName = RDMs(iRDM).name; + thisFigI = figI(iRDM); + + localOptions = MDSOptions; + localOptions.figI = thisFigI; + localOptions.fileName = ['Stimuli_MDS_' sami.util.deblank(RDMName) '_by_' colorCategory]; + localOptions.titleString = ['Stimuli MDS for RDM: "' sami.util.deunderscore(RDMName) '"']; + + disp(['*** Drawing MDS: "' localOptions.titleString '" [fig. ' num2str(thisFigI) ']']); + sami.fig.MDSArrangement(RDMs(iRDM), localOptions); + + %% saving + returnHere = pwd; + thisFileName = localOptions.fileName; + + % figure + sami.util.gotoDir(fullfile(userOptions.rootPath, 'MDSplots')); + disp([' -> saving MDS_plot FIGURE to ' fullfile(pwd, thisFileName)]); + sami.fig.handleFigure(figI, thisFileName, userOptions); + + cd(returnHere); +end + diff --git a/sami_toolbox/+sami/RDMsPairwiseCorrelations.m b/sami_toolbox/+sami/RDMsPairwiseCorrelations.m new file mode 100644 index 0000000..5448236 --- /dev/null +++ b/sami_toolbox/+sami/RDMsPairwiseCorrelations.m @@ -0,0 +1,70 @@ +function RDMsPairwiseCorrelations(RDMs, userOptions, infoTXT, figI) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + + +%% Set defaults and check options struct +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +userOptions = sami.util.setIfUnset(userOptions, 'rdms_pairWiseCorr', 'Spearman'); +userOptions = sami.util.setIfUnset(userOptions, 'fig_display', true); +userOptions = sami.util.setIfUnset(userOptions, 'fig_savePDF', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_saveFig', false); +userOptions = sami.util.setIfUnset(userOptions, 'fig_savePS', false); + +%% init +if ~exist('infoTXT','var') || isempty(infoTXT) + titleTxt = ''; + fileNameSufix = ''; +else + titleTxt = [' - ' sami.util.deunderscore(infoTXT)]; + fileNameSufix = ['_' sami.util.deblank(infoTXT)]; +end + +cmp = sami.fig.createColormap('RDMofRDMs'); +nRDMs = size(RDMs,2); + +disp(['*** Drawing RDMs [fig. ' num2str(figI) ']']); + +%% calc correlation matrix +corrMat = sami.stat.RDMCorrMat(RDMs, userOptions.rdms_pairWiseCorr); + +%% figure +sami.util.selectPlot(figI,userOptions.fig_display); +set(gcf,'Position',[200 50 500 500]); +imagesc(corrMat,[-1 1]); +axis square on; +box on; + +RDMnames = cellfun(@sami.util.deunderscore,{RDMs(:).name},'UniformOutput',false); + +set(gca,'xTick',1:nRDMs,'xTickLabel',RDMnames,'XTickLabelRotation', 90,'TickLength',[0 0],'fontsize',10); +set(gca,'yTick',1:nRDMs,'yTickLabel',RDMnames); +title(['RDM correlation matrix' titleTxt],'FontSize',12); + +colormap(cmp); +cb = colorbar; +cb.Label.String = [sami.util.deunderscore(userOptions.rdms_pairWiseCorr) ' correlation coefficient']; + + +%% saving +returnHere = pwd; +thisFileName = ['RDMsCorrMat' fileNameSufix]; + +% correlation matrix +sami.util.gotoDir(fullfile(userOptions.rootPath, 'RDMsCorrMat')); +disp([' -> saving 2nd-order correlation MATRIX to ' fullfile(pwd, [thisFileName '.mat'])]); +save([thisFileName '.mat'], 'corrMat'); + +% figure +sami.util.gotoDir(fullfile(userOptions.rootPath, 'RDMsCorrMat', 'figs')); +disp([' -> saving 2nd-order correlation FIGURE to ' fullfile(pwd, thisFileName)]); +sami.fig.handleFigure(figI, thisFileName, userOptions); + +cd(returnHere); + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/calcFeatures.m b/sami_toolbox/+sami/calcFeatures.m new file mode 100644 index 0000000..67df528 --- /dev/null +++ b/sami_toolbox/+sami/calcFeatures.m @@ -0,0 +1,91 @@ +function feat = calcFeatures(c3dData, featureType, userOptions, fileNameSufix) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +disp(['*** calculate _' featureType '_ features ***']); + +%% define default behavior +if ~exist('featureType','var') || isempty(featureType), featureType = 'idv'; end +if ~exist('fileNameSufix','var') || isempty(fileNameSufix), fileNameSufix = ''; end + +%% init vars +nFiles = size(c3dData,2); +stb = what('sami_toolbox'); +feat = struct('fSet',[]); +nActFeat = 0; % need to know how many features are already calculated + +%% fetch all functions for feature calculation +switch featureType + case 'idv' + import sami.feat.idv.* + fcn = dir(fullfile(stb.path,'+sami','+feat','+idv','*.m')); + case 'itx' + import sami.feat.itx.* + fcn = dir(fullfile(stb.path,'+sami','+feat','+itx','*.m')); + otherwise + error('*** sami:error *** "featureType" unknown.'); +end +fcn = {fcn(:).name}'; + +% ignore all functions containing "NOT" in filename +fcn(contains(fcn,'NOT')) = []; + +%% loop feature-functions +dur = nan(nFiles,numel(fcn)); +% fprintf(' ... computing features (of %d): 1',numel(fcn)); +for fcnI = 1:numel(fcn) + % display progress + fprintf(' ... loading function #%d (of %d) ... processing stimulus (of %d): 0',fcnI,numel(fcn),nFiles); + + % loop stimulus_files + for fileI = 1:nFiles + fprintf([repmat('\b',1,numel(num2str(fileI-1))) '%d'],fileI); + + tic; + tmp_feat = feval(fcn{fcnI}(1:end-2), c3dData(fileI)); + dur(fileI,fcnI) = toc; + + % save each calculated feature_sets + for i = 1:numel(tmp_feat) + feat(nActFeat+i).fSet(:,fileI) = tmp_feat(i).feat; + end + end + + % get meta_data + for i = 1:numel(tmp_feat) + feat(nActFeat+i).name = tmp_feat(i).name; + feat(nActFeat+i).color = tmp_feat(i).color; + feat(nActFeat+i).unit = tmp_feat(i).unit; + if isfield(tmp_feat,'fisherZ4paramTesting') && ~isempty(tmp_feat(i).fisherZ4paramTesting) + feat(nActFeat+i).fisherZ4paramTesting = tmp_feat(i).fisherZ4paramTesting; + else + feat(nActFeat+i).fisherZ4paramTesting = false; + end + end + + % updating variable: how many features are already calculated + nActFeat = nActFeat + numel(tmp_feat); + fprintf(' ... ok\n'); +end + +%% rename feature variable and saving it +feat_filename = ['features_' featureType fileNameSufix '.mat']; +feat_varname = ['feat_' featureType fileNameSufix]; +eval([feat_varname '= feat;']); +save(fullfile(userOptions.rootPath,feat_filename),feat_varname); + +%% finishing +fprintf([' ... features saved in rootPath as "' feat_filename '" \n']); +fprintf(' ... DONE calculating features\n\n'); + +%% some debug/information +if isfield(userOptions,'debug') && userOptions.debug == true + dur_info = [fcn, num2cell(sum(dur)')]; + save(fullfile(userOptions.rootPath,['dur_calc_feat_' featureType fileNameSufix '.mat']),'dur','dur_info'); +end + diff --git a/sami_toolbox/+sami/compareBehavCategoricalData.m b/sami_toolbox/+sami/compareBehavCategoricalData.m new file mode 100644 index 0000000..6829300 --- /dev/null +++ b/sami_toolbox/+sami/compareBehavCategoricalData.m @@ -0,0 +1,293 @@ +function stats = compareBehavCategoricalData(data, category, depVar, userOptions, figI) +% stats = compareBehavCategoricalData(data, category, depVar, userOptions, figI) +% Detailed explanation goes here +% +% IMPORTANT NOTE: this function assumes an "within-factorial" design, i.e. participants +% are rating each category, leading to a data s*p-data_matrix, where +% s=stimulus and p=participant +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +%% define default behavior +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +if ~exist('category','var') || isempty(category), category = userOptions.stimuli_sorting; end +if ~exist('depVar','var') || isempty(depVar), depVar = 'NOT_SPECIFIED'; end +userOptions = sami.util.setIfUnset(userOptions, 'behav_threshold', 0.05); +userOptions = sami.util.setIfUnset(userOptions, 'behav_multipleTesting', 'bonferroni'); +userOptions = sami.util.setIfUnset(userOptions, 'stats_minNforSubjectRFXtests', 12); +userOptions = sami.util.setIfUnset(userOptions, 'stats_minNforNonParamTests', 5); + +% check +if ismember(category, userOptions.stimuli_settings(1,2:end)) == 0, error('*** sami:ERROR *** input for ''category'' is not specified in ''userOptions.stimuli_settings''. Please check input. returning.'); end + +% print progress +disp(['*** compare ''categorical'' behavioral data [fig. ' num2str(figI) '] ***']); + +%% init vars +fileNameSufix = ['_' sami.util.deblank(depVar) '_by_' sami.util.deblank(category)]; +cmp = sami.fig.createColormap('confMatrix'); + +%% get category for 'data'-input according to 'category'-input +fileOrder = sami.util.getStimOrder(userOptions); +[~, catIdx] = ismember(category, userOptions.stimuli_settings(1,2:end)); +catOfStim = cell2mat(userOptions.stimuli_settings(fileOrder+1,catIdx + 1)); +catNums = unique(catOfStim); +nCat = numel(catNums); +nStim = numel(catOfStim); +nSubj = size(data,2); +nRatings = (nStim/nCat) * nSubj; + +%% data preprocessing +% create confusionMatrix by counting every single ratings and comparing to 'category'-input +% in 'absolute' as well as 'percentage' values +confMat_abs = zeros(nCat); +for i = 1:nStim + thisRef = catOfStim(i); + for j = 1:nSubj + thisRat = data(i,j); + confMat_abs(thisRef,thisRat) = confMat_abs(thisRef,thisRat) + 1; + end +end +confMat_perc = confMat_abs ./ nRatings .* 100; +% calc total values +confMat_abs_totalTarget = sum(confMat_abs,2); +confMat_abs_totalResponse = sum(confMat_abs,1); +confMat_perc_totalTarget = confMat_abs_totalTarget ./ sum(confMat_abs_totalTarget) .* 100; +confMat_perc_totalResponse = confMat_abs_totalResponse ./ sum(confMat_abs_totalResponse) .* 100; +% complete matrices +confMat_abs_total = [confMat_abs,confMat_abs_totalTarget;confMat_abs_totalResponse,nan]; +confMat_perc_total = [confMat_perc,confMat_perc_totalTarget;confMat_perc_totalResponse,nan]; +% save stats +stats.confMat_abs = confMat_abs; +stats.confMat_perc = confMat_perc; +stats.confMat_abs_total = confMat_abs_total; +stats.confMat_perc_total = confMat_perc_total; + +% for subject RFX tests: get accuracy and d-prime for each subject +data_acc = nan(nSubj,nCat); +data_dp = nan(nSubj,nCat); +for s = 1:nSubj + for c = catNums' + thisCatIdx = catOfStim == catNums(c); + + % accuracy per category + data_acc(s,c) = mean(data(thisCatIdx,s) == c) * 100; + + % hit + false alarm rate to calculate d' + catRated = data(:,s) == c; + catPresent = thisCatIdx; + catNotPresent = ~thisCatIdx; + + hit = sum(catRated & catPresent) / sum(catPresent); + fa = sum(catRated & catNotPresent) / sum(catNotPresent); + + % note: proportions of 0 are replaced with 0.5/N, proportions of 1 are replaced with (N-0.5)/N + % where N is the number of 'category'-trials, or 'non-category#-trials respectively + if hit == 1, hit = (sum(catPresent)-0.5)./sum(catPresent); end + if hit == 0, hit = 0.5/sum(catPresent); end + if fa == 1, fa = (sum(catNotPresent)-0.5)./sum(catNotPresent); end + if fa == 0, fa = 0.5/sum(catNotPresent); end + + data_dp(s,c) = norminv(hit)-norminv(fa); + end +end + +%% plotting : confusion matrix +sami.util.selectPlot([figI 2 2 1],userOptions.fig_display); cla; +% set(gcf,'Position',[2600 50 850 650]); +set(gcf,'Position',[200 10 830 830]); + + +image(confMat_perc,'CDataMapping','scaled') +hold on; box on; axis equal; +xlim([.5 nCat+.5]); ylim([.5 nCat+.5]); + +title(['\fontsize{11}confusion matrix: ''' depVar ''' responses']); + +xlabel('\fontsize{10}\bfbehavioral response'); +ylabel('\fontsize{10}\bftarget'); +set(gca,... + 'XTick',1:nCat,... + 'XTickLabel',userOptions.stimuli_naming_key(catIdx).condition,... + 'XTickLabelRotation',90,... + 'YTick',1:nCat,... + 'YTickLabel',userOptions.stimuli_naming_key(catIdx).condition,... + 'XAxisLocation','top',... + 'FontSize',9); + +caxis([0 100]); +colormap(cmp) +c = colorbar; +c.Label.String = '[%]'; +c.Label.FontSize = 10; + +%% plotting : confusion matrix - TABLE +sami.util.selectPlot([figI 2 2 3],userOptions.fig_display); cla; +box on; axis equal + +xlabel('\fontsize{10}\bfbehavioral response'); +ylabel('\fontsize{10}\bftarget'); +xlim([0.5 nCat+1.5]); ylim([0.5 nCat+1.5]); +set(gca,... + 'XTick',1:nCat+1,... + 'XTickLabel',[userOptions.stimuli_naming_key(catIdx).condition,{'\Sigma'}],... + 'XTickLabelRotation',90,... + 'YTick',1:nCat+1,... + 'YTickLabel',[userOptions.stimuli_naming_key(catIdx).condition,{'\Sigma'}],... + 'XAxisLocation','top',... + 'FontSize',9,... + 'TickLength',[0 0],... + 'Ydir','reverse'); + +line([0.5 nCat+1.5],[nCat+.5 nCat+.5],'Color',[.7 .7 .7],'LineWidth',0.5); +line([nCat+.5 nCat+.5],[0.5 nCat+1.5],'Color',[.7 .7 .7],'LineWidth',0.5); + +for t = 1:nCat+1 + for r = 1:nCat+1 + if ~isnan(confMat_abs_total(t,r)) + if r == t + txtB = 'bold'; + else + txtB = 'normal'; + end + str = {sprintf('%.2f%%',confMat_perc_total(t,r)),' '}; + text(r,t,str,'HorizontalAlignment','center','Fontsize',9,'fontweight',txtB); + str = {' ',['(' num2str(confMat_abs_total(t,r)) ')']}; + text(r,t,str,'HorizontalAlignment','center','Fontsize',7,'fontweight',txtB); + end + end +end + +%% plotting : accuracy + d' barplots in a loop +for plt = 1:2 + switch plt + case 1 + thisData = data_acc; + thisM = 100/nCat; + + thisFig = 2; + statName = 'accuracy'; + + thisTitleP = {['accuracy of behavioral responses (' depVar ')'],... + ['averaged across subjects (n=' num2str(nSubj) ')']}; + thisYlabelP = 'accuracy (\pm SEM) [%]'; + + thisTitleNP = {['behavioral response (' depVar ') of subjects (n=' num2str(nSubj) ')'],... + ['(subjects: averaged across ' num2str(nStim) ' stimuli)']}; + thisYlabelNP = 'accuracy [%]'; + case 2 + thisData = data_dp; + thisM = 0; + + thisFig = 4; + statName = 'dprime'; + + thisTitleP = {['sensitivity of stimulus categories (' depVar ') '],... + ['averaged across subjects (n=' num2str(nSubj) ')']}; + thisYlabelP = 'd'' (\pm SEM)'; + + thisTitleNP = {['sensitivity for stimulus categories (' depVar ') of subjects (n=' num2str(nSubj) ')'],... + ['(subjects: averaged across ' num2str(nStim) ' stimuli)']}; + thisYlabelNP = 'd'''; + end + + sami.util.selectPlot([figI 2 2 thisFig],userOptions.fig_display); cla; + + % decide what kind of statistical analysis and figures to plot with current data + if nSubj >= userOptions.stats_minNforSubjectRFXtests + % ---------- parametric --- show barplot ----------------------------------------- + sami.fig.drawBars(thisData,category,userOptions); + + % calculate one-way-ANOVA t-test + if nCat > 2 + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.rmANOVA(thisData,userOptions.behav_multipleTesting,userOptions.behav_threshold); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + if isinteger([stats.(statName).groupDiffs.df1 stats.(statName).groupDiffs.df2]) + statStr = sprintf(['%s(%d,%d) = %.2f , ' pStr ', pEta^2 = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df1, stats.(statName).groupDiffs.df2, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.pEtaSq); + else + statStr = sprintf(['%s(%.2f,%.2f) = %.2f , ' pStr ', pEta^2 = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df1, stats.(statName).groupDiffs.df2, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.pEtaSq); + end + else + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.tTest(thisData,'paired'); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s(%d) = %.2f , ' pStr ', d = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.d); + end + + % calculate one-sample t-test comparing against chance/zero + stats.(statName).testVsValue = sami.stat.tTest(thisData,'one_sample',thisM); + + % show title/ylabel + title([thisTitleP,statStr], 'fontsize',10) + ylabel(thisYlabelP); + else + % ---------- non-parameteric --- boxplot ----------------------------------------- + sami.fig.drawBoxPlot(thisData,category,userOptions); + + % calculate Kruskal-Wallis Mann-Whitney-U-Test *IF* StimPerCat > minNforNPtests + if nSubj >= userOptions.stats_minNforNonParamTests + if nCat > 2 + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.Friedman(thisData,userOptions.behav_multipleTesting,userOptions.behav_threshold); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s(%d) = %.2f , ' pStr ', W = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.KendallW); + else + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.Wilcoxon(thisData,'paired'); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s = %.2f , ' pStr ', r = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.V, stats.(statName).groupDiffs.r); + end + else + statStr = ''; + end + + % calculate one-sided Wilcoxon signed rank test comparing against chance/zero + stats.(statName).testVsValue = sami.stat.Wilcoxon(thisData,'one_sample',thisM); + + % show title/ylabel + title([thisTitleNP,statStr], 'fontsize',10) + ylabel(thisYlabelNP); + end + + % edit title and add statistics + curAx = gca; + curTitle = curAx.Title.String; + if stats.(statName).groupDiffs.p < 0.05 && nCat > 2 + multTestStr = ['\rm(threshold: ' num2str(userOptions.behav_threshold) ', multiple testing: ' userOptions.behav_multipleTesting ')']; + title([curTitle; multTestStr]); + end + + % save stats + stats.(statName).dataPoints = thisData; + stats.(statName).dataMEAN = mean(thisData); + stats.(statName).dataSEM = std(thisData) ./ sqrt(size(thisData,1)); +end + +%% save +returnHere = pwd; +thisFileName = ['compBehavCategoricalData' fileNameSufix]; + +% ANOVA results +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compBehavCategoricalData')); +disp([' -> saving STATISTIC results to ' fullfile(pwd, thisFileName)]); +save([thisFileName '.mat'], 'stats'); + +% figure +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compBehavCategoricalData','figs')); +disp([' -> saving behavioral data FIGURE to ' fullfile(pwd, thisFileName)]); +sami.fig.handleFigure(figI, thisFileName, userOptions); + +cd(returnHere); + +end + diff --git a/sami_toolbox/+sami/compareBehavData.m b/sami_toolbox/+sami/compareBehavData.m new file mode 100644 index 0000000..311c799 --- /dev/null +++ b/sami_toolbox/+sami/compareBehavData.m @@ -0,0 +1,161 @@ +function stats = compareBehavData(data, category, depVar, userOptions, figI) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +%% define default behavior +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +if ~exist('category','var') || isempty(category), category = userOptions.stimuli_sorting; end +if ~exist('depVar','var') || isempty(depVar), depVar = 'NOT_SPECIFIED'; end +userOptions = sami.util.setIfUnset(userOptions, 'behav_threshold', 0.05); +userOptions = sami.util.setIfUnset(userOptions, 'behav_multipleTesting', 'bonferroni'); +userOptions = sami.util.setIfUnset(userOptions, 'stats_minNforSubjectRFXtests', 12); +userOptions = sami.util.setIfUnset(userOptions, 'stats_minNforNonParamTests', 5); + +% check +if ismember(category, userOptions.stimuli_settings(1,2:end)) == 0, error('*** sami:ERROR *** input for ''category'' is not specified in ''userOptions.stimuli_settings''. Please check input. returning.'); end + +% print progress +disp(['*** compare behavioral data [fig. ' num2str(figI) '] ***']); + +%% init vars +fileNameSufix = ['_' sami.util.deblank(depVar) '_by_' sami.util.deblank(category)]; + +%% get category for 'data'-input according to 'category'-input +fileOrder = sami.util.getStimOrder(userOptions); +[~, catIdx] = ismember(category, userOptions.stimuli_settings(1,2:end)); +catOfStim = cell2mat(userOptions.stimuli_settings(fileOrder+1,catIdx + 1)); +catNums = unique(catOfStim); +nCat = numel(catNums); + +%% data preprocessing +% average across stimuli in category for each subject +dataSubj = nan(size(data,2),numel(catNums)); +for i = 1:nCat + dataSubj(:,i) = mean(data(catOfStim == catNums(i),:)); +end + +% average across subjects in category for each stimulus +dataStim = nan(size(data,1)/numel(catNums),numel(catNums)); +for i = 1:nCat + dataStim(:,i) = mean(data(catOfStim == catNums(i),:),2); +end + +%% loop both data_set and plot them, and do parametric/non_parametric group comparisons, if applicable +for i = 1:2 + % fetch data + switch i + case 1 + % fetch data + thisData = dataSubj; + title_info = {'subjects' 'stimuli'}; + statName = 'compSubjects'; + case 2 + % fetch data + thisData = dataStim; + title_info = {'stimuli' 'subjects'}; + statName = 'compStimuli'; + end + nThisData = size(thisData,1); + + % selectPlot + sami.util.selectPlot([figI 1 2 i],userOptions.fig_display); cla; + set(gcf,'Position',[200 50 800 400]); + + % decide what kind of plot to plot with current data + if nThisData >= userOptions.stats_minNforSubjectRFXtests + % ---------- show barplot -------------------------------------------------- + sami.fig.drawBars(thisData,category,userOptions); + + % calculate one-way-ANOVA t-test + if nCat > 2 + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.ANOVA(thisData,userOptions.behav_multipleTesting,userOptions.behav_threshold); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s(%d,%d) = %.2f , ' pStr ', eta^2 = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df1, stats.(statName).groupDiffs.df2, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.etaSq); + else + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.tTest(thisData,'paired'); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s(%d) = %.2f , ' pStr ', d = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.d); + end + + % show title/ylabel + title({['behavioral rating (' depVar ')'],... + ['averaged across ' title_info{1} ' (n=' num2str(nThisData) ')'],... + statStr},... + 'fontsize',10) + ylabel('average rating (\pm SEM)'); + else + % ---------- boxplot -------------------------------------------------- + sami.fig.drawBoxPlot(thisData,category,userOptions); + + % calculate Kruskal-Wallis Mann-Whitney-U-Test *IF* StimPerCat > minNforKW + if nThisData >= userOptions.stats_minNforNonParamTests + if nCat > 2 + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.KruskalWallis(thisData,userOptions.behav_multipleTesting,userOptions.behav_threshold); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s(%d) = %.2f , ' pStr ', eta^2 = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.etaSq); + else + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.Wilcoxon(thisData,'paired'); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s = %.2f , ' pStr ', r = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.V, stats.(statName).groupDiffs.r); + end + else + statStr = ''; + end + + % show title/ylabel + title({['behavioral rating (' depVar ') of ' title_info{1} ' (n=' num2str(nThisData) ')'],... + ['(' title_info{1} ', averaged across ' title_info{2} ')'],... + statStr},... + 'fontsize',10) + ylabel('behavioral rating values'); + end + + % edit title: add statistics and add multCompare info in case of significant differences + curAx = gca; + curTitle = curAx.Title.String; + if stats.(statName).groupDiffs.p < 0.05 && nCat > 2 + multTestStr = ['\rm(threshold: ' num2str(userOptions.behav_threshold) ', multiple testing: ' userOptions.behav_multipleTesting ')']; + title([curTitle; multTestStr]); + end + + % save stats + stats.(statName).dataPoints = thisData; + stats.(statName).dataMEAN = mean(thisData); + stats.(statName).dataSEM = std(thisData) ./ sqrt(size(thisData,1)); +end + +%% save +returnHere = pwd; +thisFileName = ['compBehavioralData' fileNameSufix]; + +% ANOVA results +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compBehavData')); +disp([' -> saving ANOVA results to ' fullfile(pwd, thisFileName)]); +save([thisFileName '.mat'], 'stats'); + +% figure +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compBehavData','figs')); +disp([' -> saving behavioral data FIGURE to ' fullfile(pwd, thisFileName)]); +sami.fig.handleFigure(figI, thisFileName, userOptions); + +cd(returnHere); + +end + diff --git a/sami_toolbox/+sami/compareBehavRDMs2FeatRDMs.m b/sami_toolbox/+sami/compareBehavRDMs2FeatRDMs.m new file mode 100644 index 0000000..cb41532 --- /dev/null +++ b/sami_toolbox/+sami/compareBehavRDMs2FeatRDMs.m @@ -0,0 +1,431 @@ +function stats = compareBehavRDMs2FeatRDMs(refRDMs_input, featRDMs_input, refName_input, featName_input, userOptions, figI) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +%% define default behavior +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +userOptions = sami.util.setIfUnset(userOptions, 'rdms_pairWiseCorr', 'Spearman'); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_relatednessTest', 'signedRank'); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_relatednessThreshold', 0.05); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_relatednessMultipleTesting', 'FDR'); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_differencesTest','signedRank'); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_differencesThreshold',0.05); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_differencesMultipleTesting','FDR'); +userOptions = sami.util.setIfUnset(userOptions, 'stats_minNforSubjectRFXtests', 12); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_nRandomisations',50000); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_orderByCorr',true); + +%% check input +if ~any(ismember({'signedRank','randomisation','none'}, userOptions.rdms_relatednessTest)) + error('*** sami:ERROR *** wrong specification in ''rdms_relatednessTest''. Please check. returning.') +end +if ~any(ismember({'signedRank','conditionBootstrap','none'}, userOptions.rdms_relatednessTest)) + error('*** sami:ERROR *** wrong specification in ''rdms_relatednessTest''. Please check. returning.') +end + +disp(['*** compare behavioralRDMs with featuresRDMs [fig. ' num2str(figI) '] ***']); + +%% init vars +fileNameSufix = ['_' refName_input '_with_' featName_input]; + +if ~exist('refName_input','var') || isempty(refName_input) + refName = 'reference'; +else + refName = [refName_input ' -']; +end +if ~exist('featName_input','var') || isempty(featName_input) + featName = 'feature'; +else + featName = featName_input; +end + +% figure style +style.barFace = [.4 .4 .4]; +style.barEdge = 'none'; + +style.ceilingFace = [.7 .7 .7]; +style.ceilingEdge = 'none'; +style.ceilingAlpha = 0.5; + +style.errorBarColor = [0 0 0]; +style.errorBarWidth = 4; +style.errorBarCapSize = 0; +style.errorBarLineStyle = 'none'; + +style.starsColor = [0 0 0]; +style.starsSize = 20; + +style.featRDMsColor = [0 0 0]; +style.featRDMsRotation = 45; + +style.vertAxisWidth = 1; +style.vertAxisColor = [0 0 0]; + +style.pwCompColors = [0 0 0;... % black + 1 .4667 0;... % orange + .7725 0 0]; % red + +%% prepare ref/feat RDMs +[refRDMs,nRefRDMs] = sami.rdm.unwrapRDMs(refRDMs_input); +meanRefRDM = mean(refRDMs,3); +[featRDMs,nFeatRDMs,featNames] = sami.rdm.unwrapRDMs(featRDMs_input); +[nStim,~,nSubj] = size(refRDMs); + +%% check if all entries in any RDM are valid +if any(isnan(refRDMs(:))) || any(isnan(featRDMs(:))) + error('*** sami:ERROR *** NANs found in one of the reference or feature RDMs. Cant''t deal with that. Please check and remove NANs from RDMs. returning.') +end + +%% estimate the ceiling for the refRDM +[ceilingUpperBound, ceilingLowerBound] = sami.stat.noiseCeilingOfAvgRDMcorr(refRDMs,userOptions.rdms_pairWiseCorr); +stats.ceiling = [ceilingLowerBound,ceilingUpperBound]; + +%% calculate the average correlation +% correlate the average of the featRDMs with each instance of the refRDM +feat2refCorrs = nan(nSubj,nFeatRDMs); +for iFeat = 1:nFeatRDMs + for iSubj = 1:nSubj + if isequal(userOptions.rdms_pairWiseCorr,'Kendall_taua') + feat2refCorrs(iSubj,iFeat) = sami.stat.rankCorr_Kendall_taua(sami.rdm.vectorizeRDMs(featRDMs(:,:,iFeat))',sami.rdm.vectorizeRDMs(refRDMs(:,:,iSubj))'); + else + feat2refCorrs(iSubj,iFeat) = corr(sami.rdm.vectorizeRDMs(featRDMs(:,:,iFeat))',sami.rdm.vectorizeRDMs(refRDMs(:,:,iSubj))','type',userOptions.rdms_pairWiseCorr,'rows','pairwise'); + end + end +end +y = mean(feat2refCorrs,1); + +% check if bars will be sorted +if userOptions.rdms_orderByCorr + [~,sortedIs] = sort(y,'descend'); +else + sortedIs = 1:nFeatRDMs; +end +featNames = featNames(sortedIs); +y_sorted = y(sortedIs); + +% save stats +stats.featNames = featNames; +stats.featRelatedness_r = y_sorted; +stats.featRelatedness_subjects_r = feat2refCorrs(:,sortedIs); + +%% decide inference procedures and test for RDM relatedness +% check if there are enough subjects for subject RFX tests +if nRefRDMs >= userOptions.stats_minNforSubjectRFXtests + subjectRFX = true; + fprintf(' -> (info) Found %d instances of the reference RDMs.\n',nRefRDMs); +else + subjectRFX = false; + fprintf(' -> (info) Found less than 12 of reference RDMs instances. Cannot do subject random-effects inference.\n'); +end + +% decide RDM-relatedness test +if isequal(userOptions.rdms_relatednessTest,'signedRank') && ~subjectRFX + userOptions.rdms_relatednessTest = 'randomisation'; +end + +% perform choosen test +switch userOptions.rdms_relatednessTest + case 'signedRank' + fprintf(' -> Performing signed-rank test of RDM relatedness (subject as random effect).\n'); + pVals = nan(nFeatRDMs,1); + for iFeat = 1:nFeatRDMs + [pVals(iFeat)] = sami.stat.signrank_onesided(feat2refCorrs(:,iFeat)); + end + pVals = pVals(sortedIs); + + % correct pVals for MC + stats.featRelatedness_p = sami.stat.multipleTesting(pVals,userOptions.rdms_relatednessMultipleTesting,userOptions.rdms_relatednessThreshold); + + case 'randomisation' + fprintf(' -> Performing condition-label randomisation test of RDM relatedness (fixed effects).\n'); + fprintf(' -> of %d randomisations: 0%%', userOptions.rdms_nRandomisations); + tic; + nRand = userOptions.rdms_nRandomisations; + rdms = nan(nFeatRDMs,nStim*(nStim-1)/2); + for iRDM = 1:nFeatRDMs + rdms(iRDM,:) = sami.rdm.vectorizeRDM(featRDMs(:,:,iRDM)); + end + + % do exhaustive permutations if number of stimuli < 8 + exPerm = false; + if nStim < 8 + allPerm = perms(1:nStim); + nRand = size(allPerm, 1); + exPerm = true; + end + + % loop randomisations + nullCorr = nan(nRand,nFeatRDMs); + for iRand = 1:nRand + if exPerm + rndIdx = allPerm(iRand, :); + else + rndIdx = randperm(nStim); + end + + rdmA_rand_vec = sami.rdm.vectorizeRDM(meanRefRDM(rndIdx,rndIdx)); + + if isequal(userOptions.rdms_pairWiseCorr,'Kendall_taua') + for iFeat = 1:nFeatRDMs + nullCorr(iRand,iFeat) = sami.stat.rankCorr_Kendall_taua(rdmA_rand_vec',rdms(iFeat,:)'); + end + else + nullCorr(iRand,:) = corr(rdmA_rand_vec',rdms','type',userOptions.rdms_pairWiseCorr,'rows','pairwise'); + end + + if mod(iRand,nRand/100) == 0 + fprintf([repmat('\b',1,numel(num2str(round(100*iRand/nRand)-1))+1) '%d%%'],round(100*iRand/nRand)) + end + end + t = toc; + fprintf([' ... DONE [in ' num2str(ceil(t)) 's]\n']); + + % calculate p-values from the randomisation test + pVals = nan(nFeatRDMs,1); + for iFeat = 1:nFeatRDMs + pVals(iFeat) = 1 - sami.stat.relRank(nullCorr(:,iFeat), y(iFeat)); + end + pVals = pVals(sortedIs); + + % correct pVals for MC + stats.featRelatedness_p = sami.stat.multipleTesting(pVals,userOptions.rdms_relatednessMultipleTesting,userOptions.rdms_relatednessThreshold); + + % display warning text + if exPerm + fprintf(' >>> WARNING <<< Comparing RDMs with fewer than 8 conditions (per conditions set) will produce unrealiable results!\n'); + fprintf(' I''ll partially compensate by using exhaustive instead of random permutations...\n'); + end + + otherwise + fprintf(' -> Not performing any test of RDM relatedness.\n'); +end + +%% decide inference type and perform feature-RDM-comparison test +% pairwise difference in average feat2refCorrs +for i = 1:nFeatRDMs + for j = 1:nFeatRDMs + stats.featDifferences_r(i,j) = y_sorted(i)-y_sorted(j); + end +end + +% decide feature-RDM-comparison test +if isequal(userOptions.rdms_differencesTest,'signedRank') + if subjectRFX + fprintf(' -> (info) using %d instances of the reference RDM for random-effects tests comparing pairs of feature RDMs.\n', nRefRDMs); + else + fprintf(' -> (info) attempting to revert to condition-bootstrap tests comparing pairs of feature RDMs.\n'); + if nStim >= 20 + fprintf(' -> (info) reverting to condition bootstrap tests for comparing pairs of feature RDMs.\n'); + userOptions.rdms_differencesTest = 'conditionBootstrap'; + else + fprintf(' -> (info) there are less than 20 conditions. can not do tests for comparing pairs of feature RDMs.\n'); + userOptions.rdms_differencesTest = 'none'; + end + end +end + +% perform choosen test +switch userOptions.rdms_differencesTest + case 'signedRank' + fprintf(' -> Performing signed-rank test for feature RDM comparisons (subject as random effect).\n'); + % do a one sided signrank test on the similarity of the each feature RDM,averaged + % over all instances with the different instances of the reference RDM + pairWisePs = nan(nFeatRDMs); + for iFeat = 1:nFeatRDMs + for jFeat = 1:nFeatRDMs + [pairWisePs(iFeat,jFeat)] = signrank(feat2refCorrs(:,iFeat),feat2refCorrs(:,jFeat),'alpha',0.05,'method','exact'); + end + end + stats.featDifferences_pUncorr = pairWisePs(sortedIs,sortedIs); + pairWisePs = pairWisePs(sortedIs,sortedIs); + + % correct pVals for MC + pUncorr = sami.rdm.vectorizeRDM(pairWisePs)'; + pCorr = sami.stat.multipleTesting(pUncorr,userOptions.rdms_differencesMultipleTesting,userOptions.rdms_differencesThreshold); + pCorr = squareform(pCorr); + pCorr(logical(eye(size(pCorr,1)))) = 1; + stats.featDifferences_p = pCorr; + stats.featDifferences_pUncorr = pairWisePs; + + % standard error + stats.featRelatedness_SEs = std(feat2refCorrs(:,sortedIs))/sqrt(nSubj); + + case 'conditionBootstrap' + fprintf(' -> Performing condition bootstrap test comparing feature RDMs (subject as random effect).\n'); + % bootstrap correlation values + bootstrapRs = sami.stat.bootstrapRDMs(refRDMs, featRDMs, userOptions); + + % find p-values based on bootstraped distributions of correlation values + pairWisePs = nan(nFeatRDMs); + for iFeat = 1:(nFeatRDMs-1) + bsI = bootstrapRs(iFeat,:); + feat2refSimsI = feat2refCorrs(:,iFeat); + + for jFeat = (iFeat+1):nFeatRDMs + bsJ = bootstrapRs(jFeat,:); + feat2refSimsJ = feat2refCorrs(:,jFeat); + + [~, ~, pairWisePs(iFeat, jFeat)] = sami.stat.bootConfint(feat2refSimsI - feat2refSimsJ, bsI - bsJ, 'two-tailed', userOptions); + pairWisePs(jFeat, iFeat) = pairWisePs(iFeat, jFeat); + end + end + stats.featDifferences_pUncorr = pairWisePs(sortedIs,sortedIs); + pairWisePs = pairWisePs(sortedIs,sortedIs); + + % correct pVals for MC + pUncorr = sami.rdm.vectorizeRDM(pairWisePs)'; + pCorr = sami.stat.multipleTesting(pUncorr,userOptions.rdms_differencesMultipleTesting,userOptions.rdms_differencesThreshold); + pCorr = squareform(pCorr); + pCorr(logical(eye(size(pCorr,1)))) = 1; + stats.featDifferences_p = pCorr; + stats.featDifferences_pUncorr = pairWisePs; + + % standard error + bootstrapEs = std(bootstrapRs, 0, 2); + stats.featRelatedness_SEs = bootstrapEs(sortedIs); + + otherwise + fprintf(' -> Not performing any test for comparing pairs of feature RDMs.\n'); +end + +%% creating figure >>> bar-plot +sami.util.selectPlot([figI 1 2 1],userOptions.fig_display); +set(gcf,'Position',[130 150 950 650]); +hold on; + +% average RDM-relatedness bars +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +for barI = 1:nFeatRDMs + patch([barI-0.4 barI-0.4 barI+0.4 barI+0.4],[0 y_sorted(barI) y_sorted(barI) 0],[-0.01 -0.01 -0.01 -0.01],... + style.barFace,... + 'edgecolor',style.barEdge); +end + +% ceiling ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +h = patch([0.1 0.1 nFeatRDMs+0.9 nFeatRDMs+0.9],... + [ceilingLowerBound ceilingUpperBound ceilingUpperBound ceilingLowerBound],... + [0.1 0.1 0.1 0.1],... + style.ceilingFace,... + 'edgecolor',style.ceilingEdge); +alpha(h,style.ceilingAlpha); + +% errorbars ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +errorbar(1:nFeatRDMs,y_sorted,stats.featRelatedness_SEs,... + 'Color',style.errorBarColor,... + 'CapSize',style.errorBarCapSize,... + 'LineWidth',style.errorBarWidth,... + 'LineStyle',style.errorBarLineStyle); + +% p-values from RDM relatedness tests ++++++++++++++++++++++++++++++++++++++++++++++++++++ +if ~isequal(userOptions.rdms_relatednessTest,'none') + pVals = stats.featRelatedness_p; + for test = 1:nFeatRDMs + if pVals(test) <= userOptions.rdms_relatednessThreshold + H = stats.featRelatedness_r(test) + stats.featRelatedness_SEs(test) + max(stats.featRelatedness_SEs); + text(test, H, ['\bf',sami.util.getStars(pVals(test))],... + 'Rotation', 0,... + 'Color', style.starsColor,... + 'FontSize',style.starsSize,... + 'HorizontalAlignment','center'); + end + end +end + +% label the bars with the names of the feature RDMs ++++++++++++++++++++++++++++++++++++++ +Ymin = min(min(y)-max(stats.featRelatedness_SEs),0); +for test = 1:nFeatRDMs + text(test, Ymin, ['\bf',sami.util.deunderscore(featNames{test})],... + 'Rotation', style.featRDMsRotation,... + 'Color', style.featRDMsColor,... + 'HorizontalAlignment','right'); +end + +% plot pretty vertical axis ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +axis off; +maxYTickI = ceil(max( [y ceilingUpperBound] ) * 10); +for YTickI = 0:maxYTickI + plot([0 0.2],[YTickI YTickI]./10,'Color',style.vertAxisColor,'LineWidth',style.vertAxisWidth); + text(0,double(YTickI/10),num2str(YTickI/10,1),'HorizontalAlignment','right'); +end +plot([0.1 0.1],[0 YTickI]./10,'k','LineWidth',style.vertAxisWidth); +% text(-1.2,double(maxYTickI/10/2),... +text(-.15*range(xlim),double(maxYTickI/10/2),... + {'\bf RDM correlation ',['\rm[',sami.util.deunderscore(userOptions.rdms_pairWiseCorr),', averaged across ',num2str(nRefRDMs),' subjects]']},... + 'HorizontalAlignment','center',... + 'Rotation',90); + +% add title/description text above +++++++++++++++++++++++++++++++++++++++++++++++++++++++ +ylim(ylim + [0 .01]); +title({['\bf\fontsize{10} Relatedness between ' refName '-RDM and each of the ' featName '-RDMs'],... + ['\fontsize{8}RDM relatedness tests: \rm',userOptions.rdms_relatednessTest,... + '\rm (threshold: ',num2str(userOptions.rdms_relatednessThreshold),', multiple testing: ',userOptions.rdms_relatednessMultipleTesting,')'],... + ['\bfsignificance:\rm *: p < ' num2str(userOptions.rdms_relatednessThreshold, '%0.2f') ', **: p < ',num2str(userOptions.rdms_relatednessThreshold/5, '%0.2f') ', ***: p < ',num2str(userOptions.rdms_relatednessThreshold/50, '%0.3f')]}); + +%% creating figure >>> matrices for pairwise feature RDM comparisons +sami.util.selectPlot([figI 1 2 2],userOptions.fig_display); + +% prepare matrix: get data and color_code according to p-values ++++++++++++++++++++++++++ +plt_thresh = userOptions.rdms_differencesThreshold; +plt_invisible = true(nFeatRDMs); +plt_invisible(stats.featDifferences_r > 0) = false; + +plt_pMat = stats.featDifferences_p; +plt_sigLevel = plt_pMat; +plt_sigLevel(plt_pMat > plt_thresh) = 1; % black +plt_sigLevel(plt_pMat <= plt_thresh & plt_pMat > plt_thresh/5) = 2; % intermediate color +plt_sigLevel(plt_pMat <= plt_thresh/5) = 3; % red + +plt_image = nan(size(plt_sigLevel,1),size(plt_sigLevel,2),3); +for k1 = 1:nFeatRDMs + for k2 = 1:nFeatRDMs + if plt_invisible(k1,k2) + plt_image(k1,k2,:) = [1 1 1]; % set every "invisible cell = white" + else + plt_image(k1,k2,:) = style.pwCompColors(plt_sigLevel(k1,k2),:); + end + end +end + +% display matrix and add tick_labels +++++++++++++++++++++++++++++++++++++++++++++++++++++ +imagesc(plt_image(1:nFeatRDMs-1,2:nFeatRDMs,:)); +axis square; +set(gca,'xTick',1:nFeatRDMs-1,... + 'xTickLabel',sami.util.deunderscore(stats.featNames(2:nFeatRDMs)),... + 'XTickLabelRotation', 90,... + 'yTick',1:nFeatRDMs-1,... + 'yTickLabel',sami.util.deunderscore(stats.featNames(1:nFeatRDMs-1)),... + 'TickLength',[0 0],... + 'fontsize',10); + +% add title/description text above +++++++++++++++++++++++++++++++++++++++++++++++++++++++ +title({['\bf\fontsize{10}Pairwise comparison of ''relatedness'' between ' featName ' RMDs'],... + ['\fontsize{8}pairwise comparison tests: \rm',userOptions.rdms_differencesTest,... + '\rm (threshold: ',num2str(userOptions.rdms_differencesThreshold),', multiple testing: ',userOptions.rdms_differencesMultipleTesting,')'],... + ['\bfsignificance:\rm black: n.s., orange: p < ' num2str(userOptions.rdms_differencesThreshold, '%0.2f') ', red: p < ',num2str(userOptions.rdms_differencesThreshold/5, '%0.2f')]}); + +% edit size of axes ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +set(gca,'Position',get(gca,'Position') .* [1.1 2 .5 1]); + +%% save +returnHere = pwd; +thisFileName = ['compare' fileNameSufix '_RDMs']; + +% correlation matrix +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compBehav2Feat')); +disp([' -> saving comparing STATISTICS to ' fullfile(pwd, thisFileName)]); +save([thisFileName '.mat'], 'stats'); + +% figure +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compBehav2Feat','figs')); +disp([' -> saving 2nd-order FIGURE to ' fullfile(pwd, thisFileName)]); +sami.fig.handleFigure(figI, thisFileName, userOptions); + +cd(returnHere); +end + + + diff --git a/sami_toolbox/+sami/compareCatRDMs2FeatRDMs.m b/sami_toolbox/+sami/compareCatRDMs2FeatRDMs.m new file mode 100644 index 0000000..81c1951 --- /dev/null +++ b/sami_toolbox/+sami/compareCatRDMs2FeatRDMs.m @@ -0,0 +1,139 @@ +function stats = compareCatRDMs2FeatRDMs(refRDMs, featRDMs, infoTXT, userOptions, figI) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% define default behavior +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +userOptions = sami.util.setIfUnset(userOptions, 'rdms_pairWiseCorr', 'Spearman'); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_pairWiseCorrMultipleTesting', 'holm'); +userOptions = sami.util.setIfUnset(userOptions, 'rdms_pairWiseCorrThreshold', 0.05); + +% print progress +disp(['*** compare categoryRDMs with featuresRDMs [fig. ' num2str(figI) '] ***']); + +%% init vars +if ~exist('infoTXT','var') || isempty(infoTXT) + titleTxt = ['compare Category and Feature RDMs (multiple testing: ' userOptions.rdms_pairWiseCorrMultipleTesting ')']; + + fNr = 1; + while exist(fullfile(userOptions.rootPath, 'compCat2Feat', 'figs', [fprintf('compCat2FeatRDMs_%04d',fNr) '.*']),'file') + fNr = fNr + 1; + end + fileNameSufix = sprintf('%04d',fNr); +else + titleTxt = [sami.util.deunderscore(infoTXT) ' (multiple testing: ' userOptions.rdms_pairWiseCorrMultipleTesting ')']; + fileNameSufix = sami.util.deblank(infoTXT); +end + +cmp = sami.fig.createColormap('RDMofRDMs'); +nRefRDMs = size(refRDMs,2); +nFeatRDMs = size(featRDMs,2); + +%% calculate correlation +[corrMat, pValMat] = sami.stat.RDMCorrMat([refRDMs featRDMs], userOptions.rdms_pairWiseCorr); + +% get values interested in +rVals = corrMat(1:nRefRDMs,nRefRDMs+1:end); +pVals = pValMat(1:nRefRDMs,nRefRDMs+1:end); + +% correct pVals for multiple testing +pValsCorr = sami.stat.multipleTesting(pVals',userOptions.rdms_pairWiseCorrMultipleTesting,userOptions.rdms_pairWiseCorrThreshold)'; + +%% ploting +% +++++++++++++++++ corr. Matrix +++++++++++++++++ +sami.util.selectPlot([figI 2 1 1],userOptions.fig_display); +set(gcf,'Position',[120 20 800 800]); + +imagesc(rVals,[-1 1]); +axis equal; box on; + +set(gca,'xTick',1:nFeatRDMs,'xTickLabel',sami.util.deunderscore({featRDMs(:).name}),'XTickLabelRotation', 45,'TickLength',[0 0],'fontsize',10); +set(gca,'yTick',1:nRefRDMs,'yTickLabel',sami.util.deunderscore({refRDMs(:).name})); + +xlim([0.5 nFeatRDMs + 0.5]); +ylim([0.5 nRefRDMs + 0.5]); + +for i = 1:nRefRDMs-1 + line([0.5 0.5+nFeatRDMs],[i i]+0.5,'Color','w','LineWidth',2) +end + +% plot sig. stars +for i = 1:nRefRDMs + for j = 1:nFeatRDMs + % text(j,i,['i' num2str(i) 'j' num2str(j)],'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',12); + if pValsCorr(i,j) < 0.05 + text(j,i,sami.util.getStars(pValsCorr(i,j)),'HorizontalAlignment','center','VerticalAlignment','middle','FontSize',22,'FontWeight','normal'); + end + end +end + +% plot colorbar +colormap(gca, cmp); +cb = colorbar; +cb.Label.String = {' correlation coefficient',['(' sami.util.deunderscore(userOptions.rdms_pairWiseCorr) ')']}; + +% plot optional title +title(titleTxt); + +% +++++++++++++++++ stats table +++++++++++++++++ +sami.util.selectPlot([figI 2 1 2],userOptions.fig_display); +hold on; axis equal; box off; + +title(['r/p - Values (multiple testing: ' userOptions.rdms_pairWiseCorrMultipleTesting ')']); +set(gca,'TickLength',[0 0],'XTickLabelRotation',45,'Ydir','reverse',... + 'XTick',1:nFeatRDMs,'XTicklabel',sami.util.deunderscore({featRDMs(:).name}),... + 'YTick',1:nRefRDMs,'YTicklabel',sami.util.deunderscore({refRDMs(:).name})); +xlim([0.5 nFeatRDMs + 0.5]); +ylim([0.5 nRefRDMs + 0.5]); + +% plot lines +for i = 1:nRefRDMs+1 + line([0.5 0.5+nFeatRDMs],[0.5 0.5]+-1+i,'Color','k'); +end +for i = 1:nFeatRDMs+1 + line([0.5 0.5]+-1+i,[0.5 0.5+nRefRDMs],'Color','k'); +end + +% text +for i = 1:nRefRDMs + for j = 1:nFeatRDMs + if pValsCorr(i,j) > 0.05 + thisFW = 'normal'; + else + thisFW = 'bold'; + end + thisTXT = {sprintf('r = %.2f',rVals(i,j)),sami.util.getPString(pValsCorr(i,j))}; + text(j,i,thisTXT,'HorizontalAlignment','center','FontSize',7,'Fontweight',thisFW); + end +end + +%% return stats +stats.r = rVals; +stats.p = pValsCorr; +stats.p_uncorr = pVals; + +%% save +returnHere = pwd; +thisFileName = ['compCat2FeatRDMs_' fileNameSufix]; + +% correlation matrix +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compCat2Feat')); +disp([' -> saving 2nd-order STATISTICS to ' fullfile(pwd, thisFileName)]); +save([thisFileName '.mat'], 'stats'); + +% figure +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compCat2Feat', 'figs')); +disp([' -> saving 2nd-order FIGURE to ' fullfile(pwd, thisFileName) '.*']); +sami.fig.handleFigure(figI, thisFileName, userOptions); + +cd(returnHere); + +end + + + diff --git a/sami_toolbox/+sami/compareFeatValues.m b/sami_toolbox/+sami/compareFeatValues.m new file mode 100644 index 0000000..c561d72 --- /dev/null +++ b/sami_toolbox/+sami/compareFeatValues.m @@ -0,0 +1,178 @@ +function stats = compareFeatValues(dataInput, category, fSetInfo, userOptions, figI) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +%% define default behavior +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +if ~exist('category','var') || isempty(category), category = userOptions.stimuli_sorting; end +if ~exist('fSetInfo','var') || isempty(fSetInfo), fSetInfo = 'NOT_SPECIFIED'; end +userOptions = sami.util.setIfUnset(userOptions, 'feat_threshold', 0.05); +userOptions = sami.util.setIfUnset(userOptions, 'feat_multipleTesting', 'bonferroni'); +userOptions = sami.util.setIfUnset(userOptions, 'stats_minNforSubjectRFXtests', 12); +userOptions = sami.util.setIfUnset(userOptions, 'stats_minNforNonParamTests', 5); + +% check +if ismember(category, userOptions.stimuli_settings(1,2:end)) == 0, error('*** sami:ERROR *** input for ''category'' is not specified in ''userOptions.stimuli_settings''. Please check input. returning.'); end + +% print progress +disp(['*** compare features [fig. ' num2str(figI) '] ***']); + +%% init vars +fileNameSufix = ['_' sami.util.deblank(fSetInfo) '_by_' sami.util.deblank(category)]; +nFeat = numel(dataInput); + +%% get category for 'data'-input according to 'category'-input +fileOrder = sami.util.getStimOrder(userOptions); +[~, catIdx] = ismember(category, userOptions.stimuli_settings(1,2:end)); +catOfStim = cell2mat(userOptions.stimuli_settings(fileOrder+1,catIdx + 1)); +catNums = unique(catOfStim); + +nCat = numel(catNums); + +% works ONLY if there is an equal amount of stimuli per category !!!!!!! +catCounts = accumarray(catOfStim,1); +if all(catCounts == catCounts(1)) + nStimPerCat = catCounts(1); +else + error('*** sami:ERROR *** toolbox is (currently) unable to deal with unequally distributed stimuli across categories. please provide equal n for each category. returning.'); +end + +%% data preprocessing +% loop each feature-set +data = nan(nStimPerCat,nCat,nFeat); +units = cell(nFeat,1); +fisherZ = false(nFeat,1); +for f = 1:nFeat + % average multivariate data to obtain one data_point per stimulus + if dataInput(f).fisherZ4paramTesting + tmpData = tanh(mean(atanh(dataInput(f).fSet),1))'; + else + tmpData = mean(dataInput(f).fSet,1)'; + end + % sort according to categories + for c = 1:nCat + data(:,c,f) = tmpData(catOfStim == catNums(c)); + end + % get meta info + units{f} = dataInput(f).unit; + fisherZ(f) = dataInput(f).fisherZ4paramTesting; +end + +%% loop feature-sets -> plot data -> do ANOVA if applicable +% prepare figure +nHorPan = ceil(sqrt(2/3 * nFeat)); +nVerPan = ceil(nFeat/nHorPan); +sami.util.selectPlot(figI,userOptions.fig_display); +set(gcf,'Position',[200 10 1500 1000]); + +for f = 1:nFeat + % select panel + sami.util.selectPlot([figI nVerPan nHorPan f],userOptions.fig_display); cla; + + % select data + statName = dataInput(f).name; + thisData = data(:,:,f); + + % decide what kind of plot to plot with current data + if nStimPerCat >= userOptions.stats_minNforSubjectRFXtests + % ---------- show barplot -------------------------------------------------- + sami.fig.drawBars(thisData,category,userOptions); + + % apply Fisher-z-transformation for 'correlational' data + if fisherZ(f) + thisData = atanh(thisData); + end + + % calculate one-way-ANOVA OR t-test + if nCat > 2 + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.ANOVA(thisData,userOptions.feat_multipleTesting,userOptions.feat_threshold); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s(%d,%d) = %.2f , ' pStr ', eta^2 = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df1, stats.(statName).groupDiffs.df2, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.etaSq); + else + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.tTest(thisData,'two_sample'); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s(%d) = %.2f , ' pStr ', d = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.d); + end + + % show title/ylabel + title({['"' sami.util.deunderscore(statName) '"'],... + ['averaged across stimuli (n=' num2str(size(thisData,1)) ')'],... + statStr},... + 'fontsize',10) + ylabel({['feature values [' units{f} ']'],'(MEAN \pm SEM)'}); + else + % ---------- boxplot -------------------------------------------------- + sami.fig.drawBoxPlot(thisData,category,userOptions); + + % calculate Kruskal-Wallis Mann-Whitney-U-Test *IF* StimPerCat > minNforKW + if nStimPerCat >= userOptions.stats_minNforNonParamTests + if nCat > 2 + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.KruskalWallis(thisData,userOptions.feat_multipleTesting,userOptions.feat_threshold); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s(%d) = %.2f , ' pStr ', eta^2 = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.df, stats.(statName).groupDiffs.V,stats.(statName).groupDiffs.etaSq); + else + % test for group difference + add pairwise significance lines + stats.(statName).groupDiffs = sami.stat.Wilcoxon(thisData,'two_sample'); + + % prepare string for title, depending on statistical test + pStr = sami.util.getPString(stats.(statName).groupDiffs.p); + statStr = sprintf(['%s = %.2f , ' pStr ', r = %.2f'], stats.(statName).groupDiffs.S, stats.(statName).groupDiffs.V, stats.(statName).groupDiffs.r); + end + else + statStr = ''; + end + + % show title/ylabel + title({['"' sami.util.deunderscore(statName) '"'],... + ['distribution of stimuli (n=' num2str(size(thisData,1)) ')'],... + statStr},... + 'fontsize',10) + ylabel(['feature values [' units{f} ']']); + end + + % edit title: add statistics and add multCompare info in case of significant differences + curAx = gca; + curTitle = curAx.Title.String; + if stats.(statName).groupDiffs.p < 0.05 && nCat > 2 + multTestStr = ['\rm(threshold: ' num2str(userOptions.feat_threshold) ', multiple testing: ' userOptions.feat_multipleTesting ')']; + title([curTitle; multTestStr]); + end + + % save stats + stats.(statName).dataPoints = thisData; + stats.(statName).dataMEAN = mean(thisData); + stats.(statName).dataSEM = std(thisData) ./ sqrt(size(thisData,1)); + stats.(statName).unit = units{f}; +end + +%% save +returnHere = pwd; +thisFileName = ['compFeatureValues' fileNameSufix]; + +% ANOVA results +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compFeatValues')); +disp([' -> saving ANOVA results to ' fullfile(pwd, thisFileName)]); +save([thisFileName '.mat'], 'stats'); + +% figure +sami.util.gotoDir(fullfile(userOptions.rootPath, 'compFeatValues','figs')); +disp([' -> saving feature data FIGURE to ' fullfile(pwd, thisFileName)]); +sami.fig.handleFigure(figI, thisFileName, userOptions); + +cd(returnHere); + +end diff --git a/sami_toolbox/+sami/createBehavioralRDMs.m b/sami_toolbox/+sami/createBehavioralRDMs.m new file mode 100644 index 0000000..d20f263 --- /dev/null +++ b/sami_toolbox/+sami/createBehavioralRDMs.m @@ -0,0 +1,67 @@ +function [RDMs, behavData] = createBehavioralRDMs(input, nameSuffix, type, userOptions, makeFig) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 09/2020 + +%% preparations +if ~exist('makeFig','var') || isempty(makeFig), makeFig = false; end +userOptions = sami.util.setIfUnset(userOptions, 'default_behavRDMcolor', [.9 .2 .1]); +userOptions = sami.util.setIfUnset(userOptions, 'behav_distanceMeasure', 'euclidean'); + + +%% load data from file +try + inputData = table2cell(readtable(input)); +catch + fprintf('\n'); error('*** sami:ERROR *** check input for creatBehavioralRDMs function. Could not read in data. returning.'); +end + +%% sort behavioral input according to userOptions +% get fileOrder for 'input' as well as 'global template' +fileOrderInput = inputData(:,1); +[~,fileOrderTemplate] = sami.util.getStimOrder(userOptions); + +% find sorting indices from input -> template order +[isMem, fileSortIdx] = ismember(fileOrderTemplate,fileOrderInput); + +% return if input-stimuli are not equal to c3d-files +if any(~isMem) + fprintf('\n'); error('*** sami:ERROR *** behavioral stimuli-names does not correspond to available c3d files. please check. returning.'); +end + +% sort input data +behavData = cell2mat(inputData(fileSortIdx,2:end)); + +%% create RDM for each subject +for s = 1:size(behavData,2) + % obtain differently calculated RDMs + [tmp_binRDM,tmp_crossRDM] = sami.rdm.categoricalRDM( behavData(:,s),[],false ); + tmp_distRDM = squareform( pdist( behavData(:,s), userOptions.behav_distanceMeasure ) ); + + % Store the RDM in a struct with the right names and things! + switch type + case 'binary' + RDMs(s).RDM = tmp_binRDM; + case 'catDiffs' + RDMs(s).RDM = tmp_crossRDM; + case 'distance' + RDMs(s).RDM = tmp_distRDM; + otherwise + fprintf('\n'); error('*** sami:ERROR *** please specify type of calculation wanted.'); + end + + RDMs(s).name = ['subj' num2str(s) '_' nameSuffix]; + RDMs(s).color = userOptions.default_behavRDMcolor; +end + +%% visualise +if makeFig + sami.plotRDMs(RDMs); +end + + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/createCategoryRDMs.m b/sami_toolbox/+sami/createCategoryRDMs.m new file mode 100644 index 0000000..26e14b9 --- /dev/null +++ b/sami_toolbox/+sami/createCategoryRDMs.m @@ -0,0 +1,71 @@ +function RDMs = createCategoryRDMs(categories, type, userOptions, categoryShift, makeFig) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% preparations +if ~exist('makeFig','var') || isempty(makeFig), makeFig = false; end +if ~exist('categoryShift','var') || isempty(categoryShift), categoryShift = []; end +userOptions = sami.util.setIfUnset(userOptions, 'default_modelRDMcolor', [.5 .5 .5]); + +%% loop defined category_combinations +nRDMs = numel(categories); +StimOrder = sami.util.getStimOrder(userOptions); + +for n = 1:nRDMs + % find indizes of categories in userOptions.stimuli_settings + this_cat_cut = regexp(categories{n},'*','split'); + + % check if categories are existing + if any(~ismember(this_cat_cut, userOptions.stimuli_settings(1,2:end))) + fprintf('\n'); error('*** sami:ERROR *** definition of categories failed. check names and if there are equal to first row in "userOptions.stimuli_settings".'); + end + + % create category vectors as input for sami.rdm.categoricalRDM + cat_data = nan(numel(StimOrder),numel(this_cat_cut)); + for c = 1:numel(this_cat_cut) + [~, sort_cat] = ismember(this_cat_cut{c}, userOptions.stimuli_settings(1,2:end)); + cat_data_raw = cell2mat(userOptions.stimuli_settings(2:end,sort_cat + 1)); + cat_data(:,c) = cat_data_raw(StimOrder); + end + + % shift categories, if wanted + shifted_title = ''; + if ~isempty(categoryShift) + for cs = 1:size(categoryShift,1) + cat_data(cat_data==categoryShift(cs,1)) = categoryShift(cs,2); + end + shifted_title = '_shifted'; + end + + % analyse categorical vectors and obtain RDMs + [tmp_binRDM,tmp_crossRDM] = sami.rdm.categoricalRDM(cat_data); + tmp_distRDM = squareform(pdist(cat_data,'euclidean')); + + % Store the RDM in a struct with the right names and things! + switch type + case 'binary' + RDMs(n).RDM = tmp_binRDM; + case 'catDiffs' + RDMs(n).RDM = tmp_crossRDM; + case 'distance' + RDMs(n).RDM = tmp_distRDM; + otherwise + fprintf('\n'); error('*** sami:ERROR *** please specify type of calculation wanted.'); + end + + RDMs(n).name = [categories{n} shifted_title]; + RDMs(n).color = userOptions.default_modelRDMcolor; +end + +%% visualise +if makeFig + sami.plotRDMs(RDMs,userOptions); +end + + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/createFeatureRDMs.m b/sami_toolbox/+sami/createFeatureRDMs.m new file mode 100644 index 0000000..4e186f6 --- /dev/null +++ b/sami_toolbox/+sami/createFeatureRDMs.m @@ -0,0 +1,31 @@ +function RDMs = createFeatureRDMs(feat, userOptions) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + + +%% define default behavior +userOptions = sami.util.setIfUnset(userOptions, 'feat_distance', 'Euclidean'); + +%% looping feature-sets +nFeat = numel(feat); +for f = 1:nFeat + % Get features for this stimulus + thisFeatures = feat(f).fSet'; + + % Calculate the RDM + localRDM = squareform( pdist( thisFeatures, userOptions.feat_distance ) ); + + % Store the RDM in a struct with the right names and things! + RDMs(f).RDM = localRDM; + RDMs(f).name = feat(f).name; + RDMs(f).color = feat(f).color; + + clear localRDM; +end + +end \ No newline at end of file diff --git a/sami_toolbox/+sami/initSAMI.m b/sami_toolbox/+sami/initSAMI.m new file mode 100644 index 0000000..30fb394 --- /dev/null +++ b/sami_toolbox/+sami/initSAMI.m @@ -0,0 +1,80 @@ +function userOptions = initSAMI(userOptions,do) +% userOptions = initSAMI(userOptions[,do]) +% +% Needs to be called at the very beginning. This function: +% - checks if necessary 'userOptions' are set +% - checks if rootPath exists and asks what the user want to do if it does +% - loads 'stimuli_settings' file, if provided by user and set in userOptions +% - creates folder 'rootPath', where everything will be saved to +% +% input: +% - userOptions: struct, containing all the options set by user beforehand +% +% optional input: +% - do: a string, specifying behavior in case of an existing rootPath folder: +% 'a' = abort +% 'd' = delete everything and start over +% 'c' = continue and use existing data/files, which can (and will) be overwritten +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 10/2020 + +%% check for needed userOptions to be set +if ~isfield(userOptions, 'analysisName'), error('*** sami:error *** analysisName must be set. See help'); end%if + +%% define default behavior +if ~exist('do','var') || isempty(do), do = 'ask'; end + +disp('*** initializing sami_toolbox ... '); + +%% set root directory of the project, if not set yet +if ~isfield(userOptions, 'rootPath') + userOptions.rootPath = fullfile(pwd,userOptions.analysisName); +end + +%% check for existing folder/data +if exist(userOptions.rootPath,'dir') + % ask what to do + while strcmp(do,'ask') || ~(strcmp(do,'a') || strcmp(do,'d') || strcmp(do,'c')) + do = input(['??? "rootPath" folder already exists ??? Please choose: \n',... + ' [a]: abort \n',... + ' d : delete everything and start over \n',... + ' c : continue and use existing data/files \n',... + ' >> '],'s'); + if isempty(do), do = 'a'; end + end + + % do the right thing + switch do + case 'a' + disp(' '); + error('... abort. doing nothing.'); + case 'd' + disp(' -> deleting existing data/files.'); + rmdir(fullfile(userOptions.rootPath), 's'); + init_folder(userOptions); + case 'c' + disp(' -> continue using existing files. Be aware: data/files will be overwritten!!!'); + end + +else + disp(' -> creating "rootPath" folder.'); + init_folder(userOptions); +end + +%% if 'userOptions.stimuli_settings_filename' is set, load this file into 'userOptions.stimuli_settings' +if ~isempty(userOptions.stimuli_settings_filename) + fprintf(' -> loading ''stimuli_settings'' file...'); + T = readtable(userOptions.stimuli_settings_filename); + userOptions.stimuli_settings = [T.Properties.VariableNames; table2cell(T)]; + fprintf(' done\n'); +end + +end + + +%% sub_function -> create folder and save userOptions +function init_folder(userOptions) + mkdir(userOptions.rootPath); + save(fullfile(userOptions.rootPath,'userOptions.mat'),'userOptions'); +end \ No newline at end of file diff --git a/sami_toolbox/+sami/loadSamiOptions.m b/sami_toolbox/+sami/loadSamiOptions.m new file mode 100644 index 0000000..4736e4b --- /dev/null +++ b/sami_toolbox/+sami/loadSamiOptions.m @@ -0,0 +1,52 @@ +function samiOptions = loadSamiOptions() +% samiOptions = loadSamiOptions() +% this function serves as some kind of 'config' file and contains variables, which are +% used at different locations across the toolbox. these variables are thought to be +% rather static, and does not to be edited by the user usually. +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + +%% definition of labels used +% this cell array also defines the order in which the markers will be sorted during the +% reading process of raw c3d-files. +samiOptions.labels = {'HEAD','LSHO','LELB','LWRI','LHIP','LKNE','LANK','RSHO','RELB','RWRI','RHIP','RKNE','RANK'}; + +%% definition of colormaps for different purposes +% each 'style' is defined by 'positions' and corresponding 'colors'. then, +% sami.fig.createColormap('style') uses this information to generate color-gradients +% between these position. +% +% - samiOptions.cmap.'style'.pos: +% double. Array defining color positions between 0 and 1. Note that the first +% position must be 0, and the last one must be 1 +% +% - samiOptions.cmap.'style'.col: +% double. Colors to place in each position. Must be a RGB matrix [n_colors x 3] + +% resolution (# of columns) of colormap +samiOptions.cmap.nCols = 256; + +% 1st order RDMs showing dissimilarities between features/particpants +samiOptions.cmap.RDMs.pos = [0 .25 .5 .75 1]; +samiOptions.cmap.RDMs.col = [ 0 0 0.5;... + 0 0.625 0.625;... + 0.75 0.75 0.75;... + 0.875 0 0;... + 1 1 0]; + +% 2nd order RDMs showing dissimilarities between 1st order RDMs +samiOptions.cmap.RDMofRDMs.pos = [ 0 .5 1]; +samiOptions.cmap.RDMofRDMs.col = [ 1 0 0;... + .8 .8 .8;... + 0 0 1]; + +% confusion matrix for categorical behavioral rating analysis +samiOptions.cmap.confMatrix.pos = [0 .01 .1 1]; +samiOptions.cmap.confMatrix.col = [ 1 1 1;... + .95 .95 .95;... + .8 .8 .8;... + 0 .27 .53]; + +end + diff --git a/sami_toolbox/+sami/plotRDMs.m b/sami_toolbox/+sami/plotRDMs.m new file mode 100644 index 0000000..4e39be5 --- /dev/null +++ b/sami_toolbox/+sami/plotRDMs.m @@ -0,0 +1,61 @@ +function RDMs = plotRDMs(RDMs, userOptions, fileNameSufix, figI, aspect) +% +% +% +% +%__________________________________________________________________________ +% A. Zabicki (azabicki@posteo.de) +% v1: 08/2020 + + +%% define default behavior +if ~exist('figI','var') || isempty(figI), figI = sami.util.getFigI(); end +if ~exist('aspect', 'var') || isempty(aspect), aspect = 2/3; end + +if ~exist('fileNameSufix','var') || isempty(fileNameSufix) + fNr = 1; + while ~isempty(dir(fullfile(userOptions.rootPath, 'RDMplots',[sprintf('RDMplot_%04d',fNr) '.*']))) + fNr = fNr + 1; + end + fileNameSufix = sprintf('%04d',fNr); +end + +%% init vars +nRDMs = numel(RDMs); +cmp = sami.fig.createColormap('RDMs'); + +disp(['*** Drawing RDMs [fig. ' num2str(figI) ']']); + +% figure +% sami.util.selectPlot(figI,userOptions.fig_display); +nHorPan = ceil(sqrt(aspect * nRDMs)); +nVerPan = ceil(nRDMs/nHorPan); + +%% display dissimilarity matrices +for iRDM = 1:nRDMs + sami.util.selectPlot([figI nVerPan nHorPan iRDM],userOptions.fig_display); cla; + set(gcf,'Position',[230 70 800 800*aspect]); + + thisRDM = RDMs(iRDM).RDM; + + % display data + image(thisRDM,'CDataMapping','scaled','AlphaData',~isnan(thisRDM)); + colormap(gca, cmp); + set(gca,'XTick',[],'YTick',[]); + title(['\bf' sami.util.deunderscore(RDMs(iRDM).name)]); + axis square off; + colorbar; +end + +%% handle figure: export and/or close it appropriately +returnHere = pwd; +thisFileName = ['RDMplot_' fileNameSufix]; + +% figure +sami.util.gotoDir(fullfile(userOptions.rootPath, 'RDMplots')); +disp([' -> saving RDMs FIGURE to ' fullfile(pwd, thisFileName) '.*']); +sami.fig.handleFigure(figI, thisFileName, userOptions); + +cd(returnHere); + +end diff --git a/sami_toolbox/documentation/will_follow.txt b/sami_toolbox/documentation/will_follow.txt new file mode 100644 index 0000000..e69de29