text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% BREATHING FILTER DETECTION TASK ANALYSIS %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% -------------------------------------------------------------------------
% Author: Olivia Harrison
% Created: 14/08/2018
%
% This software is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation, either version 3 of the License, or (at your
% option) any later version. This software is distributed in the hope that
% it will be useful, but WITHOUT ANY WARRANTY; without even the implied
% warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details: http://www.gnu.org/licenses/
% -------------------------------------------------------------------------

% This analysis script uses either Brian Maniscalco's single subect
% analysis or Steve Fleming's Hierarchical Bayesian toolbox (for group
% analysis) to calculate perceptual decision and metacognitive metrics,
% specific to data produced by running the tapas_filter_detection_task.

% IMPORTANT ANALYSIS NOTES:
%   1) The code requires JAGS software to run, which can be found here:
%      http://mcmc-jags.sourceforge.net/
%      The code has been tested on JAGS 3.4.0; there are compatibility
%      issues between matjags and JAGS 4.X.
%   2) The code requires the HMeta-d toolbox (provided as a submodule
%      from https://github.com/metacoglab/HMeta-d) to be on your path.
%   3) This script is looking for specific results files that were output
%      from the tapas_filter_detection_task, which are located in the
%      FDT/results/ folder. The results are expected to contain a
%      results.filterThreshold.xx structure, where values for xx:
%         - results.filterThreshold.filterNum: a value to determine the
%           number of filters where trials were performed
%         - results.filterThreshold.filters: a vector with 1 meaning
%           filters were presented, 0 when filters were absent (dummy)
%         - results.filterThreshold.response: a vector with 1 meaning
%           response was 'yes', 0 when response was 'no'
%         - results.filterThreshold.confidence: a vector containing
%           confidence score (1-10) on each trial
%      If your results are formatted differently, refer directly to the
%      original scripts from the HMeta-d toolbox
%      (https://github.com/metacoglab/HMeta-d).

% TO RUN THIS SCRIPT:
% type tapas_filter_detection_analysis into the MATLAB terminal from the
% main FDT folder, and follow the prompts.

% ANALYSIS OPTIONS:
% This script allows you to run either:
%    - A single subject analysis: A non-Bayesian estimation of meta-d',
%      fitting one subject at a time using a single point maximum
%      likelihood estimation (from the original Maniscalco & Lau 2012
%      paper, more info can be found here:
%      http://www.columbia.edu/~bsm2105/type2sdt/archive/index.html).
%      HOWEVER: Simulations have shown this analysis to be VERY unstable
%      when estimating meta-d' using 60 trials. It is STRONGLY encouraged
%      to utilise the hierarchical models, or collect many more trials
%      (200+ trials) for each subject to have a more reliable measure of
%      meta-d'.
%    - A group mean analysis: A hierarchical Bayesian analysis, whereby all
%      subjects are fit together and information from the group mean is
%      used as prior information for each subject. This hierarchical model
%      helps with much more accurate estimations of meta-d' with small
%      trial numbers, such as using 60 trials per subject.
%    - A group regression analysis: A hierarchical Bayesian regression
%      analysis that is an extension on the group mean analysis, whereby
%      subjects are fit together in one model, and the relationship between
%      single subject log(meta-d'/d') values and the covariate(s) is
%      estimated within the hierarchical model structure.
%    - A group difference analysis: A hierarchical Bayesian analysis,
%      whereby each group of subjects is fitted separately, and the groups
%      are then compared. Frequentist statistics (such as parametric
%      unpaired T-tests, or non-parametric Wilcoxon signed-rank tests) can
%      be used for all values that are not fitted using hierarchical
%      information, such as d', c, filter number, accuracy and average
%      confidence. As the group values for log(meta-d'/d') are calculated
%      using two separate hierarchical models, group differences are then
%      inferred via comparison of the resulting log(meta-d'/d')
%      distributions, and whether the 95% highest-density interval (HDI) of
%      the difference between the distributions spans zero. The HDI is the
%      shortest possible interval containing 95% of the MCMC samples, and
%      may not be symmetric. NOTE: If single subject values for
%      log(meta-d'/d') (or any related meta-d' metric) are required for
%      further analyses that span both groups (such as entry into general
%      linear models), it is recommended to fit all subjects together in
%      one regression model with a regressor denoting group identity.
%    - A session difference analysis (paired group difference): A
%      hierarchical Bayesian analysis, whereby two sessions / measures from
%      the same participants are fitted in a single model using a
%      multivariate normal distribution. This distribution allows for the
%      non-independence between sessions for each participant. NOTE:
%      Participants must be listed in the same order for analysis, and must
%      have data for both sessions / each measure.

% NOTES ON NON-BAYESIAN SINGLE SUBJECT ANALYSIS FROM AUTHORS' WEBSITE
% (columbia.edu/~bsm2105/type2sdt/archive/index.html):
% This analysis is intended to quantify metacognitive sensitivity (i.e. the
% efficacy with which confidence ratings discriminate between correct and
% incorrect judgments) in a signal detection theory framework. A central
% idea is that primary task performance can influence metacognitive
% sensitivity, and it is informative to take this influence into account.
% Description of the methodology can be found here: Maniscalco, B., & Lau,
% H. (2012). A signal detection theoretic approach for estimating
% metacognitive sensitivity from confidence ratings. Consciousness and
% Cognition, 21(1), 422–430. doi:10.1016/j.concog.2011.09.021
% If you use these analysis files, please reference the Consciousness &
% Cognition paper and website. Brian Maniscalco: brian@psych.columbia.edu

% NOTES ON HIERARCHICAL TOOLBOX FROM ORIGINAL SCRIPTS (FOR GROUP ANALYSES):
% This MATLAB toolbox implements the meta-d’ model (Maniscalco & Lau, 2012)
% in a hierarchical Bayesian framework using Matlab and JAGS, a program for
% conducting MCMC inference on arbitrary Bayesian models. A paper with more
% details on the method and the advantages of estimating meta-d’ in a
% hierarchal Bayesian framework is available here https://academic.oup.com/
% nc/article/doi/10.1093/nc/nix007/3748261/HMeta-d-hierarchical-Bayesian-
% estimation-of. For a more general introduction to Bayesian models of
% cognition see Lee & Wagenmakers, Bayesian Cognitive Modeling: A Practical
% Course http://bayesmodels.com/. This code is being released with a
% permissive open-source license. You should feel free to use or adapt the
% utility code as long as you follow the terms of the license provided. If
% you use the toolbox in a publication we ask that you cite the following
% paper: Fleming, S.M. (2017) HMeta-d: hierarchical Bayesian estimation of
% metacognitive efficiency from confidence ratings, Neuroscience of
% Consciousness, 3(1) nix007, https://doi.org/10.1093/nc/nix007. Copyright
% (c) 2017, Stephen Fleming. For more information and/or licences, please
% see the original code: https://github.com/metacoglab/HMeta-d.

% KEY ANALYSIS OUTPUT VALUES:
% The primary measures to come out of this analysis can be found in
% analysis.{single/groupMean/groupDiff}.xx, and they are:
%   1) xx = filterNum: Number of filters. Less filters means more sensitive
%           to changes in breathing.
%   2) xx = d1: d prime, discriminability between filter and dummy trials.
%           Larger d' means more discriminability at specific filter number
%   3) xx = c1: Decision criterion, or where the decision boundary exists.
%           Negative criterion values mean bias towards 'yes' (lower, more
%           liberal criterion), while positive values mean bias towards
%           'no' response (higher, more strngent criterion).
%   4) xx = meta_d: A measure of metacognition, or 'type 2' sensitivity.
%           This reflects how much information, in signal-to-noise units,
%           is available for metacognition (see Maniscalco & Lau, 2012).
%           NB: ISSUES WITH THIS ESTIMATION (AND ALL RELATED ESTIMATIONS)
%           IF USING SINGLE SUBJECT ANALYSES, AS DESCRIBED ABOVE.
%   5) xx = Mratio: Meta-d' to d' ratio, as a measure of metacognitive
%           'efficiency' (see Rouault et al., 2018).
%   6) xx = log_Mratio: log(meta_d/d1), reported in papers to help with
%           normalisation of data (see Rouault et al., 2018).
%   7) xx = avgConfidence: A second measure of metacognition, thought to be
%           independent from meta-d'. NB: Average confidence should only be
%           compared between two groups if there is no difference in d'
%           between the groups (i.e. task difficulty was comparable).

% ADDITIONAL GROUP DIFFERENCE OUTPUT VALUES:
% If the analysis is a two-group or two-session (paired) difference, the
% following will also be calculated for the non-hierarchical measures
% (test = groupDiff or sessionDiff; xx = filterNum, d1, c1 and
% avgConfidence):
%   1) analysis.test.xx.h = Results of null hypothesis test, where h =
%      1 for rejection of the null, and h = 0 for no rejection.
%   2) analysis.test.xx.p = The p-value for the statistical test.
%   3) analysis.test.xx.stats = Further statistics (such as tstats,
%      number of samples, degrees of freedom) associated with the test.
%   4) analysis.test.xx.ci = Confidence interval of the difference,
%      calculated only when T-test is specified.
% The highest density interval (HDI) for the difference in log(meta-d'/d')
% between the groups / sessions will be calculated and recorded, as
% frequentist statistics cannot be used here. A summary Figure for each of
% these metrics will be created and saved in the FDT/analysis/ folder.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET UP THE ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function tapas_filter_detection_analysis()

% Check folder location is main FDT folder
[~,dir_name] = fileparts(pwd);
if ~strcmp(dir_name,'FDT')
   error('Not currently in main FDT folder. Please move to FDT folder and try again.');
end

% Add relevant paths
addpath('analysis');

% Display setup on screen
fprintf('\n________________________________________\n\n    SET UP ANALYSIS FOR FILTER TASK\n________________________________________\n');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CHOOSE THE TYPE OF ANALYSIS TO RUN AND SPECIFY FILES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

try
    analysis.type = input('Type of analysis (single or group) = ', 's'); % Ask for type of analysis to be run
catch
    fprintf('\nOption does not exist, please try again.\n');
    return
end

try
    if strcmp(analysis.type,'single') == 1 % If single subject analysis specified
        analysis.PPIDs = input('PPID = ', 's'); % Ask for PPID
        try
            analysis.data = load(fullfile('results', ['filter_task_results_', analysis.PPIDs, '.mat'])); % Load data
            analysis.ratings = analysis.data.results.setup.confidenceUpper - analysis.data.results.setup.confidenceLower + 1; % Calculate number of confidence rating bins
        catch
            fprintf('\nInvalid PPID.\n');
            return
        end
    elseif strcmp(analysis.type,'group') == 1 % If group analysis specified
        analysis.type = input('Type of analysis (mean, diff, paired or regress) = ', 's'); % Ask for type of group analysis to be run
        if strcmp(analysis.type,'mean') == 1 % If group mean analysis specified
            analysis.PPIDs = input('Input group PPIDs (e.g. {''001'', ''002'', ''003'',...}) = '); % Ask for PPIDs
            analysis.groupsize(1) = length(analysis.PPIDs);
            for n = 1:length(analysis.PPIDs)
                PPID = char(analysis.PPIDs(n));
                try
                    analysis.data(n) = load(fullfile('results', ['filter_task_results_', PPID, '.mat'])); % Load data
                catch
                    fprintf('\nInvalid PPIDs.\n');
                    return
                end
            end
            analysis.ratings = analysis.data(1).results.setup.confidenceUpper - analysis.data(1).results.setup.confidenceLower + 1; % Calculate number of confidence rating bins
        elseif strcmp(analysis.type,'diff') == 1 % If two group difference analysis specified
            analysis.PPIDs.group1 = input('Input group 1 PPIDs (e.g. {''001'', ''002'', ''003'',...}) = '); % Ask for PPIDs
            analysis.PPIDs.group2 = input('Input group 2 PPIDs (e.g. {''004'', ''005'', ''006'',...}) = '); % Ask for PPIDs
            analysis.groupDiff.test = input('Type 1 variable tests (ttest or wilcoxon) = ', 's'); % Ask for parametric or non-parmetric test options for type 1 variables (d', c etc.)
            analysis.groupsize(1) = length(analysis.PPIDs.group1);
            analysis.groupsize(2) = length(analysis.PPIDs.group2);
            for n = 1:analysis.groupsize(1)
                PPID = char(analysis.PPIDs.group1(n));
                try
                    analysis.group1.data(n) = load(fullfile('results', ['filter_task_results_', PPID, '.mat'])); % Load data
                catch
                    fprintf('\nInvalid PPIDs for group 1.\n');
                    return
                end
            end
            for n = 1:analysis.groupsize(2)
                PPID = char(analysis.PPIDs.group2(n));
                try
                    analysis.group2.data(n) = load(fullfile('results', ['filter_task_results_', PPID, '.mat'])); % Load data
                catch
                    fprintf('\nInvalid PPIDs for group 2.\n');
                    return
                end
            end
            analysis.ratings = analysis.group1.data(1).results.setup.confidenceUpper - analysis.group1.data(1).results.setup.confidenceLower + 1; % Calculate number of confidence rating bins
        elseif strcmp(analysis.type,'paired') == 1 % If paired difference analysis specified
            fprintf('\nNOTE: PPIDs from each session need to be paired and in the same order\n');
            analysis.PPIDs.session1 = input('Input session 1 PPIDs (e.g. {''001a'', ''002a'', ''003a'',...}) = '); % Ask for PPIDs
            analysis.PPIDs.session2 = input('Input session 2 PPIDs (e.g. {''001b'', ''002b'', ''003b'',...}) = '); % Ask for PPIDs
            analysis.sessionDiff.test = input('Type 1 variable tests (ttest or wilcoxon) = ', 's'); % Ask for parametric or non-parmetric test options for type 1 variables (d', c etc.)
            analysis.sessionSize(1) = length(analysis.PPIDs.session1);
            analysis.sessionSize(2) = length(analysis.PPIDs.session2);
            if analysis.sessionSize(1) ~= analysis.sessionSize(2)
                error('Number of PPIDs in each session does not match!');
            end
            for n = 1:analysis.sessionSize(1)
                PPID = char(analysis.PPIDs.session1(n));
                try
                    analysis.session1.data(n) = load(fullfile('results', ['filter_task_results_', PPID, '.mat'])); % Load data
                catch
                    fprintf('\nInvalid PPIDs for session 1.\n');
                    return
                end
            end
            for n = 1:analysis.sessionSize(2)
                PPID = char(analysis.PPIDs.session2(n));
                try
                    analysis.session2.data(n) = load(fullfile('results', ['filter_task_results_', PPID, '.mat'])); % Load data
                catch
                    fprintf('\nInvalid PPIDs for session 2.\n');
                    return
                end
            end
            analysis.ratings = analysis.session1.data(1).results.setup.confidenceUpper - analysis.session1.data(1).results.setup.confidenceLower + 1; % Calculate number of confidence rating bins
        elseif strcmp(analysis.type,'regress') == 1 % If regression analysis specified
            analysis.PPIDs = input('Input group PPIDs (e.g. {''001'', ''002'', ''003'',...}) = '); % Ask for PPIDs
            analysis.groupsize(1) = length(analysis.PPIDs);
            for n = 1:length(analysis.PPIDs)
                PPID = char(analysis.PPIDs(n));
                try
                    analysis.data(n) = load(fullfile('results', ['filter_task_results_', PPID, '.mat'])); % Load data
                catch
                    fprintf('\nInvalid PPIDs.\n');
                    return
                end
            end
            analysis.ratings = analysis.data(1).results.setup.confidenceUpper - analysis.data(1).results.setup.confidenceLower + 1; % Calculate number of confidence rating bins
            fprintf('\nNote: Covariate text file required to be placed in results folder\n--> Scores need to be in the SAME ORDER as PPID input order.\n');
            analysis.covariate.fileName = input('Input covariate file name (e.g. covariate_example.txt) = ', 's'); % Ask for covariate file name
            try
                analysis.covariate.data = load(fullfile('results', analysis.covariate.fileName)); % Load data
            catch
                fprintf('\nCovariate file cannot be found in results folder.\n');
                return
            end
            [a,b] = size(analysis.covariate.data);
            if a > b % Re-shape vector if needed
                analysis.covariate.data = analysis.covariate.data';
            end
        end
    end
catch
    fprintf('\nInvalid.\n');
    return
end

% Save results
if strcmp(analysis.type,'single') == 1
    resultsFile = fullfile('analysis', ['filter_task_analysis_', analysis.type, analysis.PPIDs]); % Create figure file name
else
    resultsFile = fullfile('analysis', ['filter_task_analysis_', analysis.type]); % Create figure file name
end
save(resultsFile, 'analysis'); % Save results

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN TRIALS2COUNTS SCRIPT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(analysis.type,'single') == 1 || strcmp(analysis.type,'mean') == 1 || strcmp(analysis.type,'regress') == 1

    % Run trials2counts script for single subject or whole group
    for n = 1:length(analysis.data)
        stimID = analysis.data(n).results.thresholdTrials.filters;
        response = analysis.data(n).results.thresholdTrials.response;
        rating = analysis.data(n).results.thresholdTrials.confidence;
        [analysis.trials2counts.nR_S1{n}, analysis.trials2counts.nR_S2{n}] = trials2counts(stimID,response,rating,analysis.ratings,0); % Run for number of ratings and no cell padding
    end

elseif strcmp(analysis.type,'diff') == 1

    % Run trials2counts script for group 1
    for n = 1:length(analysis.group1.data)
        stimID = analysis.group1.data(n).results.thresholdTrials.filters;
        response = analysis.group1.data(n).results.thresholdTrials.response;
        rating = analysis.group1.data(n).results.thresholdTrials.confidence;
        [analysis.group1.trials2counts.nR_S1{n}, analysis.group1.trials2counts.nR_S2{n}] = trials2counts(stimID,response,rating,analysis.ratings,0); % Run for number of ratings and no cell padding
    end
    % Run trials2counts script for group 2
    for n = 1:length(analysis.group2.data)
        stimID = analysis.group2.data(n).results.thresholdTrials.filters;
        response = analysis.group2.data(n).results.thresholdTrials.response;
        rating = analysis.group2.data(n).results.thresholdTrials.confidence;
        [analysis.group2.trials2counts.nR_S1{n}, analysis.group2.trials2counts.nR_S2{n}] = trials2counts(stimID,response,rating,analysis.ratings,0); % Run for number of ratings and no cell padding
    end

elseif strcmp(analysis.type,'paired') == 1

    % Run trials2counts script for session 1
    for n = 1:length(analysis.session1.data)
        stimID = analysis.session1.data(n).results.thresholdTrials.filters;
        response = analysis.session1.data(n).results.thresholdTrials.response;
        rating = analysis.session1.data(n).results.thresholdTrials.confidence;
        [analysis.trials2counts.nR_S1(1).counts{n}, analysis.trials2counts.nR_S2(1).counts{n}] = trials2counts(stimID,response,rating,analysis.ratings,0); % Run for number of ratings and no cell padding
    end
    % Run trials2counts script for session 2
    for n = 1:length(analysis.session2.data)
        stimID = analysis.session2.data(n).results.thresholdTrials.filters;
        response = analysis.session2.data(n).results.thresholdTrials.response;
        rating = analysis.session2.data(n).results.thresholdTrials.confidence;
        [analysis.trials2counts.nR_S1(2).counts{n}, analysis.trials2counts.nR_S2(2).counts{n}] = trials2counts(stimID,response,rating,analysis.ratings,0); % Run for number of ratings and no cell padding
    end

end

% Save results
save(resultsFile, 'analysis');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN SINGLE SUBJECT ANALYSIS IF SPECIFIED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(analysis.type,'single') == 1 % If single subject analysis specified

    % Pad ratings so that zero counts do not interfere with fit
    analysis.trials2counts.nR_S1_padded = analysis.trials2counts.nR_S1{1} + 1/(2*10);
    analysis.trials2counts.nR_S2_padded = analysis.trials2counts.nR_S2{1} + 1/(2*10);

    % Fit data
    analysis.single.fit = fit_meta_d_MLE(analysis.trials2counts.nR_S1_padded, analysis.trials2counts.nR_S2_padded);

    % Pull out key variables
    analysis.single.filterNum = mean(analysis.data.results.thresholdTrials.filterNum);
    analysis.single.accuracy = analysis.data.results.thresholdTrials.accuracyTotal;
    analysis.single.d1 = analysis.single.fit.d1;
    analysis.single.c1 = analysis.single.fit.c1;
    analysis.single.meta_d = analysis.single.fit.meta_d;
    analysis.single.Mratio = analysis.single.fit.M_ratio;
    analysis.single.fit.log_Mratio = log(analysis.single.fit.M_ratio);
    analysis.single.avgConfidence = mean(analysis.data(1).results.thresholdTrials.confidence);

    % Save results
    save(resultsFile, 'analysis'); % Save results

    % Display completion message on screen
    fprintf('\n________________________________________\n\n   COMPLETED SINGLE SUBJECT ANALYSIS\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    fprintf('\n             MRATIO = %.2f\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n', analysis.single.Mratio);
    fprintf('\n        RESULTS CAN BE FOUND IN: \n             FDT/analysis/\n________________________________________\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN GROUP MEAN ANALYSIS IF SPECIFIED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(analysis.type,'mean') == 1 % If group mean analysis specified

    % Specify parameters
    analysis.groupMean.mcmc_params = fit_meta_d_params;

    % Fit group data all at once
    analysis.groupMean.fit = fit_meta_d_mcmc_group(analysis.trials2counts.nR_S1, analysis.trials2counts.nR_S2, analysis.groupMean.mcmc_params);

    % Calculate log of Mratio for single subjects values
    for n = 1:length(analysis.groupMean.fit.Mratio)
        analysis.groupMean.log_Mratio.singleSubject(n) = log(analysis.groupMean.fit.Mratio(n));
    end

    % Pull out group mean values
    analysis.groupMean.d1.groupMean = mean(analysis.groupMean.fit.d1);
    analysis.groupMean.d1.singleSubject = analysis.groupMean.fit.d1;
    analysis.groupMean.c1.groupMean = mean(analysis.groupMean.fit.c1);
    analysis.groupMean.c1.singleSubject = analysis.groupMean.fit.c1;
    analysis.groupMean.meta_d.groupMean = mean(analysis.groupMean.fit.meta_d);
    analysis.groupMean.meta_d.singleSubject = analysis.groupMean.fit.meta_d;
    analysis.groupMean.Mratio.groupMean = exp(analysis.groupMean.fit.mu_logMratio);
    analysis.groupMean.Mratio.singleSubject = analysis.groupMean.fit.Mratio;
    analysis.groupMean.Mratio.hdi = calc_HDI(exp(analysis.groupMean.fit.mcmc.samples.mu_logMratio(:)));
    analysis.groupMean.log_Mratio.groupMean = analysis.groupMean.fit.mu_logMratio;
    analysis.groupMean.log_Mratio.hdi = calc_HDI(analysis.groupMean.fit.mcmc.samples.mu_logMratio(:));

    % Calculate average confidence, add filter number and accuracy to analysis results
    for n = 1:length(analysis.data)
        analysis.groupMean.avgConfidence.singleSubject(n) = mean(analysis.data(n).results.thresholdTrials.confidence);
        analysis.groupMean.filterNum.singleSubject(n) = mean(analysis.data(n).results.thresholdTrials.filterNum);
        analysis.groupMean.accuracy.singleSubject(n) = analysis.data(n).results.thresholdTrials.accuracyTotal;
    end
    analysis.groupMean.avgConfidence.groupMean = mean(analysis.groupMean.avgConfidence.singleSubject);
    analysis.groupMean.filterNum.groupMean = mean(analysis.groupMean.filterNum.singleSubject);
    analysis.groupMean.accuracy.groupMean = mean(analysis.groupMean.accuracy.singleSubject);

    % Save results
    save(resultsFile, 'analysis'); % Save results

    % Display completion message on screen
    fprintf('\n________________________________________\n\n     COMPLETED GROUP MEAN ANALYSIS\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    fprintf('\n GROUP MRATIO = %.2f (HDI: %.2f to %.2f)\n', analysis.groupMean.Mratio.groupMean, analysis.groupMean.Mratio.hdi(1), analysis.groupMean.Mratio.hdi(2));
    if (analysis.groupMean.Mratio.hdi(1) > 0 && analysis.groupMean.Mratio.hdi(2) > 0) || (analysis.groupMean.Mratio.hdi(1) < 0 && analysis.groupMean.Mratio.hdi(2) < 0)
        fprintf('   HDI DOES NOT SPAN ZERO: MRATIO SIG.\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    else
        fprintf('   HDI SPANS ZERO: MRATIO NOT SIG.\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    end
    fprintf('\n        RESULTS CAN BE FOUND IN: \n             FDT/analysis/\n________________________________________\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN GROUP DIFFERENCE ANALYSIS (UNPAIRED) IF SPECIFIED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(analysis.type,'diff') == 1 % If group difference analysis specified

    % For two-group difference analysis: Print message about groups being fit separately --> no t-tests possible on meta-d' metrics
    fprintf('\nNOTE: Groups are fit using separate hierarchical models.\nFrequentist statistics (such as t-tests) are not possible on any metrics involving meta-d''.\nFrequentist statistics will still be applied to all type 1 metrics and average confidence values.\n');

    % Specify parameters
    analysis.groupDiff.mcmc_params = fit_meta_d_params;

    % Fit each group using a separate hierarchical Bayesian model
    analysis.group1.fit = fit_meta_d_mcmc_group(analysis.group1.trials2counts.nR_S1, analysis.group1.trials2counts.nR_S2, analysis.groupDiff.mcmc_params);
    analysis.group2.fit = fit_meta_d_mcmc_group(analysis.group2.trials2counts.nR_S1, analysis.group2.trials2counts.nR_S2, analysis.groupDiff.mcmc_params);

    % Compute HDI of difference for log(meta-d'/d')
    analysis.group1.logMratio.groupMean = analysis.group1.fit.mu_logMratio;
    analysis.group1.logMratio.singleSubject = log(analysis.group1.fit.Mratio);
    analysis.group1.Mratio.groupMean = exp(analysis.group1.fit.mu_logMratio);
    analysis.group1.Mratio.singleSubject = analysis.group1.fit.Mratio;
    analysis.group2.logMratio.groupMean = analysis.group2.fit.mu_logMratio;
    analysis.group2.logMratio.singleSubject = log(analysis.group2.fit.Mratio);
    analysis.group2.Mratio.groupMean = exp(analysis.group2.fit.mu_logMratio);
    analysis.group2.Mratio.singleSubject = analysis.group2.fit.Mratio;
    analysis.groupDiff.log_Mratio.sampleDiff = analysis.group1.fit.mcmc.samples.mu_logMratio - analysis.group2.fit.mcmc.samples.mu_logMratio;
    analysis.groupDiff.log_Mratio.mean = analysis.group1.fit.mu_logMratio - analysis.group2.fit.mu_logMratio;
    analysis.groupDiff.log_Mratio.hdi = calc_HDI(analysis.groupDiff.log_Mratio.sampleDiff(:));
    analysis.groupDiff.Mratio.mean = exp(analysis.group1.fit.mu_logMratio) - exp(analysis.group2.fit.mu_logMratio);
    analysis.groupDiff.Mratio.hdi = calc_HDI((exp(analysis.group1.fit.mcmc.samples.mu_logMratio(:)) - exp(analysis.group2.fit.mcmc.samples.mu_logMratio(:))));

    % Pull out group mean values for all non-hierarchical metrics
    analysis.group1.d1.groupMean = mean(analysis.group1.fit.d1);
    analysis.group1.d1.singleSubject = analysis.group1.fit.d1;
    analysis.group1.c1.groupMean = mean(analysis.group1.fit.c1);
    analysis.group1.c1.singleSubject = analysis.group1.fit.c1;
    analysis.group2.d1.groupMean = mean(analysis.group2.fit.d1);
    analysis.group2.d1.singleSubject = analysis.group2.fit.d1;
    analysis.group2.c1.groupMean = mean(analysis.group2.fit.c1);
    analysis.group2.c1.singleSubject = analysis.group2.fit.c1;

    % Calculate average confidence, add filter number and accuracy to results for group 1
    for n = 1:length(analysis.group1.data)
        analysis.group1.avgConfidence.singleSubject(n) = mean(analysis.group1.data(n).results.thresholdTrials.confidence);
        analysis.group1.filterNum.singleSubject(n) = mean(analysis.group1.data(n).results.thresholdTrials.filterNum);
        analysis.group1.accuracy.singleSubject(n) = analysis.group1.data(n).results.thresholdTrials.accuracyTotal;
    end
    analysis.group1.avgConfidence.groupMean = mean(analysis.group1.avgConfidence.singleSubject);
    analysis.group1.filterNum.groupMean = mean(analysis.group1.filterNum.singleSubject);
    analysis.group1.accuracy.groupMean = mean(analysis.group1.accuracy.singleSubject);

    % Calculate average confidence, add filter number and accuracy to results for group 2
    for n = 1:length(analysis.group2.data)
        analysis.group2.avgConfidence.singleSubject(n) = mean(analysis.group2.data(n).results.thresholdTrials.confidence);
        analysis.group2.filterNum.singleSubject(n) = mean(analysis.group2.data(n).results.thresholdTrials.filterNum);
        analysis.group2.accuracy.singleSubject(n) = analysis.group2.data(n).results.thresholdTrials.accuracyTotal;
    end
    analysis.group2.avgConfidence.groupMean = mean(analysis.group2.avgConfidence.singleSubject);
    analysis.group2.filterNum.groupMean = mean(analysis.group2.filterNum.singleSubject);
    analysis.group2.accuracy.groupMean = mean(analysis.group2.accuracy.singleSubject);

    % Calculate group difference values for all non-hierarchical metrics
    analysis.groupDiff.filterNum.meanDiff = analysis.group1.filterNum.groupMean - analysis.group2.filterNum.groupMean;
    analysis.groupDiff.accuracy.meanDiff = analysis.group1.accuracy.groupMean - analysis.group2.accuracy.groupMean;
    analysis.groupDiff.d1.meanDiff = analysis.group1.d1.groupMean - analysis.group2.d1.groupMean;
    analysis.groupDiff.c1.meanDiff = analysis.group1.c1.groupMean - analysis.group2.c1.groupMean;
    analysis.groupDiff.avgConfidence.meanDiff = analysis.group1.avgConfidence.groupMean - analysis.group2.avgConfidence.groupMean;

    % Run specified statistical test on all non-hierarchical metrics
    if strcmp(analysis.groupDiff.test,'wilcoxon') == 1
        analysis.groupDiff.testType = 'Wilcoxon rank sum test';
        [analysis.groupDiff.filterNum.p,analysis.groupDiff.filterNum.h,analysis.groupDiff.filterNum.stats] = ranksum(analysis.group1.filterNum.singleSubject,analysis.group2.filterNum.singleSubject);
        [analysis.groupDiff.accuracy.p,analysis.groupDiff.accuracy.h,analysis.groupDiff.accuracy.stats] = ranksum(analysis.group1.accuracy.singleSubject,analysis.group2.accuracy.singleSubject);
        [analysis.groupDiff.d1.p,analysis.groupDiff.d1.h,analysis.groupDiff.d1.stats] = ranksum(analysis.group1.d1.singleSubject,analysis.group2.d1.singleSubject);
        [analysis.groupDiff.c1.p,analysis.groupDiff.c1.h,analysis.groupDiff.c1.stats] = ranksum(analysis.group1.c1.singleSubject,analysis.group2.c1.singleSubject);
        [analysis.groupDiff.avgConfidence.p,analysis.groupDiff.avgConfidence.h,analysis.groupDiff.avgConfidence.stats] = ranksum(analysis.group1.avgConfidence.singleSubject,analysis.group2.avgConfidence.singleSubject);
    elseif strcmp(analysis.groupDiff.test,'ttest') == 1
        analysis.groupDiff.testType = 'Unpaired T-test';
        [analysis.groupDiff.filterNum.h,analysis.groupDiff.filterNum.p,analysis.groupDiff.filterNum.ci,analysis.groupDiff.filterNum.stats] = ttest2(analysis.group1.filterNum.singleSubject,analysis.group2.filterNum.singleSubject);
        [analysis.groupDiff.accuracy.h,analysis.groupDiff.accuracy.p,analysis.groupDiff.accuracy.ci,analysis.groupDiff.accuracy.stats] = ttest2(analysis.group1.accuracy.singleSubject,analysis.group2.accuracy.singleSubject);
        [analysis.groupDiff.d1.h,analysis.groupDiff.d1.p,analysis.groupDiff.d1.ci,analysis.groupDiff.d1.stats] = ttest2(analysis.group1.d1.singleSubject,analysis.group2.d1.singleSubject);
        [analysis.groupDiff.c1.h,analysis.groupDiff.c1.p,analysis.groupDiff.c1.ci,analysis.groupDiff.c1.stats] = ttest2(analysis.group1.c1.singleSubject,analysis.group2.c1.singleSubject);
        [analysis.groupDiff.avgConfidence.h,analysis.groupDiff.avgConfidence.p,analysis.groupDiff.avgConfidence.ci,analysis.groupDiff.avgConfidence.stats] = ttest2(analysis.group1.avgConfidence.singleSubject,analysis.group2.avgConfidence.singleSubject);
    end

    % Save results
    save(resultsFile, 'analysis'); % Save results

    % Create Figure to display results
    figure
    set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 0.6 0.6]);
    labels = ['Group 1'; 'Group 2'];
    % Plot Type 1 metrics on top line
    ax1 = subplot(2,4,1);
    bar(ax1, [analysis.group1.filterNum.groupMean; analysis.group2.filterNum.groupMean]);
    hold on
    errorbar(ax1, [analysis.group1.filterNum.groupMean; analysis.group2.filterNum.groupMean], [std(analysis.group1.filterNum.singleSubject), std(analysis.group2.filterNum.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.groupDiff.filterNum.h == 0
        title('FILTER NUMBER: NON-SIG DIFF')
    elseif analysis.groupDiff.filterNum.h > 0
        title('FILTER NUMBER: SIG DIFF')
    end
    hold off
    ax2 = subplot(2,4,2);
    bar(ax2, [analysis.group1.accuracy.groupMean; analysis.group2.accuracy.groupMean]);
    hold on
    errorbar(ax2, [analysis.group1.accuracy.groupMean; analysis.group2.accuracy.groupMean], [std(analysis.group1.accuracy.singleSubject), std(analysis.group2.accuracy.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.groupDiff.accuracy.h == 0
        title('ACCURACY: NON-SIG DIFF')
    elseif analysis.groupDiff.accuracy.h > 0
        title('ACCURACY: SIG DIFF')
    end
    ax3 = subplot(2,4,3);
    bar(ax3, [analysis.group1.d1.groupMean; analysis.group2.d1.groupMean]);
    hold on
    errorbar(ax3, [analysis.group1.d1.groupMean; analysis.group2.d1.groupMean], [std(analysis.group1.d1.singleSubject), std(analysis.group2.d1.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.groupDiff.d1.h == 0
        title('d'': NON-SIG DIFF')
    elseif analysis.groupDiff.d1.h > 0
        title('d'': SIG DIFF')
    end
    hold off
    ax4 = subplot(2,4,4);
    bar(ax4, [analysis.group1.c1.groupMean; analysis.group2.c1.groupMean]);
    hold on
    errorbar(ax4, [analysis.group1.c1.groupMean; analysis.group2.c1.groupMean], [std(analysis.group1.c1.singleSubject), std(analysis.group2.c1.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.groupDiff.c1.h == 0
        title('CRITERION: NON-SIG DIFF')
    elseif analysis.groupDiff.c1.h > 0
        title('CRITERION: SIG DIFF')
    end
    hold off
    % Plot group Mratio and the difference
    ax5 = subplot(2,4,5);
    histogram(ax5, exp(analysis.group1.fit.mcmc.samples.mu_logMratio))
    xlim([0 2])
    xlabel('MRatio');
    ylabel('Sample count');
    hold on
    histogram(ax5, exp(analysis.group2.fit.mcmc.samples.mu_logMratio))
    hold off
    lgd = legend('Group 1', 'Group 2', 'Location', 'northeast');
    lgd.FontSize = 10;
    if (analysis.groupDiff.Mratio.hdi(1) > 0 && analysis.groupDiff.Mratio.hdi(2) > 0) || (analysis.groupDiff.Mratio.hdi(1) < 0 && analysis.groupDiff.Mratio.hdi(2) < 0)
        title('MRATIO: SIG DIFF')
    else
        title('MRATIO: NON-SIG DIFF')
    end
    ax6 = subplot(2,4,6);
    histogram(ax6, analysis.group1.fit.mcmc.samples.mu_logMratio)
    xlabel('log(MRatio)');
    ylabel('Sample count');
    hold on
    histogram(ax6, analysis.group2.fit.mcmc.samples.mu_logMratio)
    hold off
    lgd = legend('Group 1', 'Group 2', 'Location', 'northwest');
    lgd.FontSize = 10;
    if (analysis.groupDiff.log_Mratio.hdi(1) > 0 && analysis.groupDiff.log_Mratio.hdi(2) > 0) || (analysis.groupDiff.log_Mratio.hdi(1) < 0 && analysis.groupDiff.log_Mratio.hdi(2) < 0)
        title('LOG(MRATIO): SIG DIFF')
    else
        title('LOG(MRATIO): NON-SIG DIFF')
    end
    ax7 = subplot(2,4,7);
    histogram(ax7, analysis.groupDiff.log_Mratio.sampleDiff)
    xlabel('log(MRatio)');
    ylabel('Sample count');
    if (analysis.groupDiff.log_Mratio.hdi(1) > 0 && analysis.groupDiff.log_Mratio.hdi(2) > 0) || (analysis.groupDiff.log_Mratio.hdi(1) < 0 && analysis.groupDiff.log_Mratio.hdi(2) < 0)
        title('LOG(MRATIO) DIFF: SIG DIFF')
    else
        title('LOG(MRATIO) DIFF: NON-SIG DIFF')
    end
    hold on
    ln2 = line([analysis.groupDiff.log_Mratio.hdi(1) analysis.groupDiff.log_Mratio.hdi(1)], [0 1800]);
    ln2.Color = 'r';
    ln2.LineWidth = 1.5;
    ln2.LineStyle = '--';
    ln2 = line([analysis.groupDiff.log_Mratio.hdi(2) analysis.groupDiff.log_Mratio.hdi(2)], [0 1800]);
    ln2.Color = 'r';
    ln2.LineWidth = 1.5;
    ln2.LineStyle = '--';
    hold off
    % Plot average confidence
    ax8 = subplot(2,4,8);
    bar(ax8, [analysis.group1.avgConfidence.groupMean; analysis.group2.avgConfidence.groupMean]);
    hold on
    errorbar(ax8, [analysis.group1.avgConfidence.groupMean; analysis.group2.avgConfidence.groupMean], [std(analysis.group1.avgConfidence.singleSubject), std(analysis.group2.avgConfidence.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.groupDiff.avgConfidence.h == 0
        title('CONFIDENCE: NON-SIG DIFF')
    elseif analysis.groupDiff.avgConfidence.h > 0
        title('CONFIDENCE: SIG DIFF')
    end
    hold off

    % Print figure
    figureFile = fullfile('analysis', ['filter_task_analysis_', analysis.type]); % Create figure file name
    print(figureFile, '-dtiff');

    % Display completion message on screen
    fprintf('\n________________________________________\n\n  COMPLETED GROUP DIFFERENCE ANALYSIS\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    fprintf('\n DIFF MRATIO = %.2f (HDI: %.2f to %.2f)\n', analysis.groupDiff.Mratio.mean, analysis.groupDiff.Mratio.hdi(1), analysis.groupDiff.Mratio.hdi(2));
    if (analysis.groupDiff.Mratio.hdi(1) > 0 && analysis.groupDiff.Mratio.hdi(2) > 0) || (analysis.groupDiff.Mratio.hdi(1) < 0 && analysis.groupDiff.Mratio.hdi(2) < 0)
        fprintf('HDI DOES NOT SPAN ZERO: MRATIO DIFF SIG.\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    else
        fprintf('  HDI SPANS ZERO: MRATIO DIFF NOT SIG.\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    end
    fprintf('\n        RESULTS CAN BE FOUND IN: \n             FDT/analysis/\n________________________________________\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN SESSION DIFFERENCE ANALYSIS (PAIRED) IF SPECIFIED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(analysis.type,'paired') == 1 % If session difference analysis specified

    % For two-session difference analysis: Print message about groups being fit together --> t-tests possible on meta-d' metrics
    fprintf('\nNOTE: Groups are fit using one hierarchical model.\nFrequentist statistics (such as t-tests) are possible on any metrics involving meta-d''.\nFrequentist statistics are applied to all type 1 metrics and average confidence values.\n');

    % Specify parameters
    analysis.sessionDiff.mcmc_params = fit_meta_d_params;

    % Fit each group using a separate hierarchical Bayesian model
    analysis.sessionDiff.fit = fit_meta_d_mcmc_groupCorr(analysis.trials2counts.nR_S1, analysis.trials2counts.nR_S2, analysis.sessionDiff.mcmc_params);

    % Compute HDI of difference for log(meta-d'/d')
    analysis.session1.logMratio.sessionMean = analysis.sessionDiff.fit.mu_logMratio(1);
    analysis.session1.logMratio.singleSubject = log(analysis.sessionDiff.fit.Mratio(:,1));
    analysis.session1.Mratio.sessionMean = exp(analysis.sessionDiff.fit.mu_logMratio(1));
    analysis.session1.Mratio.singleSubject = analysis.sessionDiff.fit.Mratio(:,1);
    analysis.session2.logMratio.sessionMean = analysis.sessionDiff.fit.mu_logMratio(2);
    analysis.session2.logMratio.singleSubject = log(analysis.sessionDiff.fit.Mratio(:,2));
    analysis.session2.Mratio.sessionMean = exp(analysis.sessionDiff.fit.mu_logMratio(2));
    analysis.session2.Mratio.singleSubject = analysis.sessionDiff.fit.Mratio(:,2);
    analysis.sessionDiff.log_Mratio.sampleDiff = analysis.sessionDiff.fit.mcmc.samples.mu_logMratio(:,:,1) - analysis.sessionDiff.fit.mcmc.samples.mu_logMratio(:,:,2);
    analysis.sessionDiff.log_Mratio.mean = analysis.session1.logMratio.sessionMean - analysis.session2.logMratio.sessionMean;
    analysis.sessionDiff.log_Mratio.hdi = calc_HDI(analysis.sessionDiff.log_Mratio.sampleDiff(:));
    analysis.sessionDiff.Mratio.mean = exp(analysis.session1.logMratio.sessionMean) - exp(analysis.session2.logMratio.sessionMean);
    analysis.sessionDiff.Mratio.hdi = calc_HDI((exp(analysis.sessionDiff.fit.mcmc.samples.mu_logMratio(:,:,1)) - exp(analysis.sessionDiff.fit.mcmc.samples.mu_logMratio(:,:,2))));

    % Pull out session mean values for all non-hierarchical metrics
    analysis.session1.d1.sessionMean = mean(analysis.sessionDiff.fit.d1(:,1));
    analysis.session1.d1.singleSubject = analysis.sessionDiff.fit.d1(:,1);
    analysis.session1.c1.sessionMean = mean(analysis.sessionDiff.fit.c1(:,1));
    analysis.session1.c1.singleSubject = analysis.sessionDiff.fit.c1(:,1);
    analysis.session2.d1.sessionMean = mean(analysis.sessionDiff.fit.d1(:,2));
    analysis.session2.d1.singleSubject = analysis.sessionDiff.fit.d1(:,2);
    analysis.session2.c1.sessionMean = mean(analysis.sessionDiff.fit.c1(:,2));
    analysis.session2.c1.singleSubject = analysis.sessionDiff.fit.c1(:,2);

    % Calculate average confidence, add filter number and accuracy to results for session 1
    for n = 1:length(analysis.session1.data)
        analysis.session1.avgConfidence.singleSubject(n) = mean(analysis.session1.data(n).results.thresholdTrials.confidence);
        analysis.session1.filterNum.singleSubject(n) = mean(analysis.session1.data(n).results.thresholdTrials.filterNum);
        analysis.session1.accuracy.singleSubject(n) = analysis.session1.data(n).results.thresholdTrials.accuracyTotal;
    end
    analysis.session1.avgConfidence.sessionMean = mean(analysis.session1.avgConfidence.singleSubject);
    analysis.session1.filterNum.sessionMean = mean(analysis.session1.filterNum.singleSubject);
    analysis.session1.accuracy.sessionMean = mean(analysis.session1.accuracy.singleSubject);

    % Calculate average confidence, add filter number and accuracy to results for session 2
    for n = 1:length(analysis.session2.data)
        analysis.session2.avgConfidence.singleSubject(n) = mean(analysis.session2.data(n).results.thresholdTrials.confidence);
        analysis.session2.filterNum.singleSubject(n) = mean(analysis.session2.data(n).results.thresholdTrials.filterNum);
        analysis.session2.accuracy.singleSubject(n) = analysis.session2.data(n).results.thresholdTrials.accuracyTotal;
    end
    analysis.session2.avgConfidence.sessionMean = mean(analysis.session2.avgConfidence.singleSubject);
    analysis.session2.filterNum.sessionMean = mean(analysis.session2.filterNum.singleSubject);
    analysis.session2.accuracy.sessionMean = mean(analysis.session2.accuracy.singleSubject);

    % Calculate session difference values for all non-hierarchical metrics
    analysis.sessionDiff.filterNum.meanDiff = analysis.session1.filterNum.sessionMean - analysis.session2.filterNum.sessionMean;
    analysis.sessionDiff.accuracy.meanDiff = analysis.session1.accuracy.sessionMean - analysis.session2.accuracy.sessionMean;
    analysis.sessionDiff.d1.meanDiff = analysis.session1.d1.sessionMean - analysis.session2.d1.sessionMean;
    analysis.sessionDiff.c1.meanDiff = analysis.session1.c1.sessionMean - analysis.session2.c1.sessionMean;
    analysis.sessionDiff.avgConfidence.meanDiff = analysis.session1.avgConfidence.sessionMean - analysis.session2.avgConfidence.sessionMean;

    % Run specified statistical test on all non-hierarchical metrics
    if strcmp(analysis.sessionDiff.test,'wilcoxon') == 1
        analysis.sessionDiff.testType = 'Wilcoxon signed rank test';
        [analysis.sessionDiff.filterNum.p,analysis.sessionDiff.filterNum.h,analysis.sessionDiff.filterNum.stats] = signrank(analysis.session1.filterNum.singleSubject,analysis.session2.filterNum.singleSubject);
        [analysis.sessionDiff.accuracy.p,analysis.sessionDiff.accuracy.h,analysis.sessionDiff.accuracy.stats] = signrank(analysis.session1.accuracy.singleSubject,analysis.session2.accuracy.singleSubject);
        [analysis.sessionDiff.d1.p,analysis.sessionDiff.d1.h,analysis.sessionDiff.d1.stats] = signrank(analysis.session1.d1.singleSubject,analysis.session2.d1.singleSubject);
        [analysis.sessionDiff.c1.p,analysis.sessionDiff.c1.h,analysis.sessionDiff.c1.stats] = signrank(analysis.session1.c1.singleSubject,analysis.session2.c1.singleSubject);
        [analysis.sessionDiff.avgConfidence.p,analysis.sessionDiff.avgConfidence.h,analysis.sessionDiff.avgConfidence.stats] = signrank(analysis.session1.avgConfidence.singleSubject,analysis.session2.avgConfidence.singleSubject);
    elseif strcmp(analysis.sessionDiff.test,'ttest') == 1
        analysis.sessionDiff.testType = 'Paired T-test';
        [analysis.sessionDiff.filterNum.h,analysis.sessionDiff.filterNum.p,analysis.sessionDiff.filterNum.ci,analysis.sessionDiff.filterNum.stats] = ttest(analysis.session1.filterNum.singleSubject,analysis.session2.filterNum.singleSubject);
        [analysis.sessionDiff.accuracy.h,analysis.sessionDiff.accuracy.p,analysis.sessionDiff.accuracy.ci,analysis.sessionDiff.accuracy.stats] = ttest(analysis.session1.accuracy.singleSubject,analysis.session2.accuracy.singleSubject);
        [analysis.sessionDiff.d1.h,analysis.sessionDiff.d1.p,analysis.sessionDiff.d1.ci,analysis.sessionDiff.d1.stats] = ttest(analysis.session1.d1.singleSubject,analysis.session2.d1.singleSubject);
        [analysis.sessionDiff.c1.h,analysis.sessionDiff.c1.p,analysis.sessionDiff.c1.ci,analysis.sessionDiff.c1.stats] = ttest(analysis.session1.c1.singleSubject,analysis.session2.c1.singleSubject);
        [analysis.sessionDiff.avgConfidence.h,analysis.sessionDiff.avgConfidence.p,analysis.sessionDiff.avgConfidence.ci,analysis.sessionDiff.avgConfidence.stats] = ttest(analysis.session1.avgConfidence.singleSubject,analysis.session2.avgConfidence.singleSubject);
    end

    % Save results
    save(resultsFile, 'analysis'); % Save results

    % Create Figure to display results
    figure
    set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 0.6 0.6]);
    labels = ['Session 1'; 'Session 2'];
    % Plot Type 1 metrics on top line
    ax1 = subplot(2,4,1);
    bar(ax1, [analysis.session1.filterNum.sessionMean; analysis.session2.filterNum.sessionMean]);
    hold on
    errorbar(ax1, [analysis.session1.filterNum.sessionMean; analysis.session2.filterNum.sessionMean], [std(analysis.session1.filterNum.singleSubject), std(analysis.session2.filterNum.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.sessionDiff.filterNum.h == 0
        title('FILTER NUMBER: NON-SIG DIFF')
    elseif analysis.sessionDiff.filterNum.h > 0
        title('FILTER NUMBER: SIG DIFF')
    end
    hold off
    ax2 = subplot(2,4,2);
    bar(ax2, [analysis.session1.accuracy.sessionMean; analysis.session2.accuracy.sessionMean]);
    hold on
    errorbar(ax2, [analysis.session1.accuracy.sessionMean; analysis.session2.accuracy.sessionMean], [std(analysis.session1.accuracy.singleSubject), std(analysis.session2.accuracy.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.sessionDiff.accuracy.h == 0
        title('ACCURACY: NON-SIG DIFF')
    elseif analysis.sessionDiff.accuracy.h > 0
        title('ACCURACY: SIG DIFF')
    end
    ax3 = subplot(2,4,3);
    bar(ax3, [analysis.session1.d1.sessionMean; analysis.session2.d1.sessionMean]);
    hold on
    errorbar(ax3, [analysis.session1.d1.sessionMean; analysis.session2.d1.sessionMean], [std(analysis.session1.d1.singleSubject), std(analysis.session2.d1.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.sessionDiff.d1.h == 0
        title('d'': NON-SIG DIFF')
    elseif analysis.sessionDiff.d1.h > 0
        title('d'': SIG DIFF')
    end
    hold off
    ax4 = subplot(2,4,4);
    bar(ax4, [analysis.session1.c1.sessionMean; analysis.session2.c1.sessionMean]);
    hold on
    errorbar(ax4, [analysis.session1.c1.sessionMean; analysis.session2.c1.sessionMean], [std(analysis.session1.c1.singleSubject), std(analysis.session2.c1.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.sessionDiff.c1.h == 0
        title('CRITERION: NON-SIG DIFF')
    elseif analysis.sessionDiff.c1.h > 0
        title('CRITERION: SIG DIFF')
    end
    hold off
    % Plot session Mratio and the difference
    ax5 = subplot(2,4,5);
    histogram(ax5, exp(analysis.sessionDiff.fit.mcmc.samples.mu_logMratio(:,:,1)))
    xlim([0 2])
    xlabel('MRatio');
    ylabel('Sample count');
    hold on
    histogram(ax5, exp(analysis.sessionDiff.fit.mcmc.samples.mu_logMratio(:,:,2)))
    hold off
    lgd = legend('Session 1', 'Session 2', 'Location', 'northeast');
    lgd.FontSize = 10;
    if (analysis.sessionDiff.Mratio.hdi(1) > 0 && analysis.sessionDiff.Mratio.hdi(2) > 0) || (analysis.sessionDiff.Mratio.hdi(1) < 0 && analysis.sessionDiff.Mratio.hdi(2) < 0)
        title('MRATIO: SIG DIFF')
    else
        title('MRATIO: NON-SIG DIFF')
    end
    ax6 = subplot(2,4,6);
    histogram(ax6, analysis.sessionDiff.fit.mcmc.samples.mu_logMratio(:,:,1))
    xlabel('log(MRatio)');
    ylabel('Sample count');
    hold on
    histogram(ax6, analysis.sessionDiff.fit.mcmc.samples.mu_logMratio(:,:,2))
    hold off
    lgd = legend('Session 1', 'Session 2', 'Location', 'northwest');
    lgd.FontSize = 10;
    if (analysis.sessionDiff.log_Mratio.hdi(1) > 0 && analysis.sessionDiff.log_Mratio.hdi(2) > 0) || (analysis.sessionDiff.log_Mratio.hdi(1) < 0 && analysis.sessionDiff.log_Mratio.hdi(2) < 0)
        title('LOG(MRATIO): SIG DIFF')
    else
        title('LOG(MRATIO): NON-SIG DIFF')
    end
    ax7 = subplot(2,4,7);
    histogram(ax7, analysis.sessionDiff.log_Mratio.sampleDiff)
    xlabel('log(MRatio)');
    ylabel('Sample count');
    if (analysis.sessionDiff.log_Mratio.hdi(1) > 0 && analysis.sessionDiff.log_Mratio.hdi(2) > 0) || (analysis.sessionDiff.log_Mratio.hdi(1) < 0 && analysis.sessionDiff.log_Mratio.hdi(2) < 0)
        title('LOG(MRATIO) DIFF: SIG DIFF')
    else
        title('LOG(MRATIO) DIFF: NON-SIG DIFF')
    end
    hold on
    ln2 = line([analysis.sessionDiff.log_Mratio.hdi(1) analysis.sessionDiff.log_Mratio.hdi(1)], [0 1800]);
    ln2.Color = 'r';
    ln2.LineWidth = 1.5;
    ln2.LineStyle = '--';
    ln2 = line([analysis.sessionDiff.log_Mratio.hdi(2) analysis.sessionDiff.log_Mratio.hdi(2)], [0 1800]);
    ln2.Color = 'r';
    ln2.LineWidth = 1.5;
    ln2.LineStyle = '--';
    hold off
    % Plot average confidence
    ax8 = subplot(2,4,8);
    bar(ax8, [analysis.session1.avgConfidence.sessionMean; analysis.session2.avgConfidence.sessionMean]);
    hold on
    errorbar(ax8, [analysis.session1.avgConfidence.sessionMean; analysis.session2.avgConfidence.sessionMean], [std(analysis.session1.avgConfidence.singleSubject), std(analysis.session2.avgConfidence.singleSubject)], 'o', 'marker', 'none', 'linewidth', 1, 'Color','k');
    set(gca,'XTickLabel',labels);
    set(gca,'XTickLabelRotation',90)
    if analysis.sessionDiff.avgConfidence.h == 0
        title('CONFIDENCE: NON-SIG DIFF')
    elseif analysis.sessionDiff.avgConfidence.h > 0
        title('CONFIDENCE: SIG DIFF')
    end
    hold off

    % Print figure
    figureFile = fullfile('analysis', ['filter_task_analysis_', analysis.type]); % Create figure file name
    print(figureFile, '-dtiff');

    % Display completion message on screen
    fprintf('\n________________________________________\n\n  COMPLETED SESSION DIFFERENCE ANALYSIS\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    fprintf('\n DIFF MRATIO = %.2f (HDI: %.2f to %.2f)\n', analysis.sessionDiff.Mratio.mean, analysis.sessionDiff.Mratio.hdi(1), analysis.sessionDiff.Mratio.hdi(2));
    if (analysis.sessionDiff.Mratio.hdi(1) > 0 && analysis.sessionDiff.Mratio.hdi(2) > 0) || (analysis.sessionDiff.Mratio.hdi(1) < 0 && analysis.sessionDiff.Mratio.hdi(2) < 0)
        fprintf('HDI DOES NOT SPAN ZERO: MRATIO DIFF SIG.\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    else
        fprintf('  HDI SPANS ZERO: MRATIO DIFF NOT SIG.\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    end
    fprintf('\n        RESULTS CAN BE FOUND IN: \n             FDT/analysis/\n________________________________________\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RUN GROUP REGRESSION ANALYSIS IF SPECIFIED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(analysis.type,'regress') == 1 % If group regression analysis specified

    % Specify parameters
    analysis.groupRegression.mcmc_params = fit_meta_d_params;

    % Standardise covariate
    analysis.groupRegression.cov = zscore(analysis.covariate.data);

    % Fit group data all at once
    analysis.groupRegression.fit = fit_meta_d_mcmc_regression(analysis.trials2counts.nR_S1, analysis.trials2counts.nR_S2, analysis.groupRegression.cov, analysis.groupRegression.mcmc_params);

    % Calculate log of Mratio for single subjects values
    for n = 1:length(analysis.groupRegression.fit.Mratio)
        analysis.groupRegression.log_Mratio.singleSubject(n) = log(analysis.groupRegression.fit.Mratio(n));
    end

    % Pull out group and single subject values
    analysis.groupRegression.d1.groupMean = mean(analysis.groupRegression.fit.d1);
    analysis.groupRegression.d1.singleSubject = analysis.groupRegression.fit.d1;
    analysis.groupRegression.c1.groupMean = mean(analysis.groupRegression.fit.c1);
    analysis.groupRegression.c1.singleSubject = analysis.groupRegression.fit.c1;
    analysis.groupRegression.meta_d.groupMean = mean(analysis.groupRegression.fit.meta_d);
    analysis.groupRegression.meta_d.singleSubject = analysis.groupRegression.fit.meta_d;
    analysis.groupRegression.Mratio.groupMean = exp(analysis.groupRegression.fit.mu_logMratio);
    analysis.groupRegression.Mratio.singleSubject = analysis.groupRegression.fit.Mratio;
    analysis.groupRegression.Mratio.hdi = calc_HDI(exp(analysis.groupRegression.fit.mcmc.samples.mu_logMratio(:)));
    analysis.groupRegression.log_Mratio.groupMean = analysis.groupRegression.fit.mu_logMratio;
    analysis.groupRegression.log_Mratio.hdi = calc_HDI(analysis.groupRegression.fit.mcmc.samples.mu_logMratio(:));
    analysis.groupRegression.beta1.groupMean = analysis.groupRegression.fit.mu_beta1;
    analysis.groupRegression.beta1.hdi = calc_HDI(analysis.groupRegression.fit.mcmc.samples.mu_beta1(:));

    % Calculate average confidence, add filter number and accuracy to analysis results
    for n = 1:length(analysis.data)
        analysis.groupRegression.avgConfidence.singleSubject(n) = mean(analysis.data(n).results.thresholdTrials.confidence);
        analysis.groupRegression.filterNum.singleSubject(n) = mean(analysis.data(n).results.thresholdTrials.filterNum);
        analysis.groupRegression.accuracy.singleSubject(n) = analysis.data(n).results.thresholdTrials.accuracyTotal;
    end
    analysis.groupRegression.avgConfidence.groupMean = mean(analysis.groupRegression.avgConfidence.singleSubject);
    analysis.groupRegression.filterNum.groupMean = mean(analysis.groupRegression.filterNum.singleSubject);
    analysis.groupRegression.accuracy.groupMean = mean(analysis.groupRegression.accuracy.singleSubject);

    % Save results
    save(resultsFile, 'analysis'); % Save results

    % Display completion message on screen plus results
    fprintf('\n________________________________________\n\n  COMPLETED GROUP REGRESSION ANALYSIS\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    fprintf('\n GROUP MRATIO = %.2f (HDI: %.2f to %.2f)\n', analysis.groupRegression.Mratio.groupMean, analysis.groupRegression.Mratio.hdi(1), analysis.groupRegression.Mratio.hdi(2));
    if analysis.groupRegression.Mratio.hdi(1) > 0
        fprintf('   HDI DOES NOT SPAN ZERO: MRATIO SIG.\n');
    else
        fprintf('   HDI SPANS ZERO: MRATIO NOT SIG.\n');
    end
    fprintf('\n GROUP BETA = %.2f (HDI: %.2f to %.2f)\n', analysis.groupRegression.beta1.groupMean, analysis.groupRegression.beta1.hdi(1), analysis.groupRegression.beta1.hdi(2));
    if analysis.groupRegression.beta1.hdi(1) > 0
        fprintf('    HDI DOES NOT SPAN ZERO: BETA SIG.\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    else
        fprintf('      HDI SPANS ZERO: BETA NOT SIG.\n _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _\n');
    end
    fprintf('\n        RESULTS CAN BE FOUND IN: \n             FDT/analysis/\n________________________________________\n');

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: translationalneuromodeling

repo: tapas


text:
% importevent() - Import experimental events from data file or Matlab
%                 array into a structure.
%
% Usage: >> eventstruct = importevent( event, oldevent, srate);
%        >> eventstruct = importevent( event, oldevent, srate, 'key1', 'value1', ...);
%
% Input:
%   event     - [ 'filename'|array ] Filename of a text file, or name of
%               Matlab array in the global workspace containing an
%               array of events in the folowing format: The first column
%               is the type of the event, the second the latency.
%               The others are user-defined. The function can read
%               either numeric or text entries in ascii files.
%   oldevent  - Old event structure. Used for aligning new events with old
%               ones. Enter [] is no such structure exist.
%   srate     - Sampling rate of the underlying data. Event latencies are
%               expressed in terms of latency in sample point with respect
%               to the data.
%
% Optional file or array input:
%  'fields'   - [Cell array] List of the name of each user-defined column,
%               optionally followed by a description. Ex: { 'type', 'latency' }
%  'skipline' - [Interger] Number of header rows to skip in the text file
%  'timeunit' - [ latency unit rel. to seconds ]. Default unit is 1 = seconds.
%               NaN indicates that the latencies are given in time points.
%  'delim'    - [string] String of delimiting characters in the input file.
%               Default is tab|space.
%
% Optional oldevent input:
%  'append'   - ['yes'|'no'] 'yes' = Append events to the current events in
%               the EEG dataset {default}: 'no' = Erase the previous events.
%  'indices'  - {integer vector] Vector indicating the indices of the events to
%               modify.
%  'align'    - [num] Align the first event latency to the latency of existing
%               event number (num), and check latency consistency.
%               A value of 0 indicates that the first events of the pre-defined
%               events and imported events will be aligned. A positive value (num)
%               aligns the first event to the num-th pre-existing event.
%               A negative value can also be used; then event number (-num)
%               is aligned to the first pre-existing event.  Default is 0.
%               (NaN-> no alignment).
%  'optimalign' - ['on'|'off'] Optimize the sampling rate of the new events so
%               they best align with old events. Default is 'on'.
%
% Outputs:
%   eventstruct - Event structure containing imported events
%
% Example: >> eventstruct = importevent( 'event_values.txt', [], 256, ...
%         'fields', {'type', 'latency','condition' }, ...
%         'append', 'no', 'align', 0, 'timeunit', 1E-3 );
%
%         This loads the ascii file 'event_values.txt' containing 3 columns
%         (event_type, latency, and condition). Latencies in the file are
%         in ms (1E-3). The first event latency is re-aligned with the
%         beginning of the dataset ('align', 0). Any previous events
%         are erased ('append', 'no').
%
% Author: Arnaud Delorme & Scott Makeig, CNL / Salk Institute, 2004
%
% See also: pop_importevent()

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, 2004, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%% REVISION HISTORY
%% function definition, input check
function event = importevent(event, oldevent, srate, varargin)

if nargin < 1
   help importevent;
   return;
end;

I = [];

% remove the event field
% ----------------------
if ~isempty(oldevent), allfields = fieldnames(oldevent);
else                   allfields = {}; end;

g = finputcheck( varargin, { 'fields'    'cell'     []         {};
                         'skipline'  'integer'  [0 Inf]    0;
                         'indices'   'integer'  [1 Inf]    [];
                         'append'    'string'   {'yes' 'no' '''yes''' '''no''' }         'yes';
                         'timeunit'  'real'     []         1;
                         'event'     { 'cell' 'real' 'string' }     []    [];
                         'align'     'integer'  []         NaN;
                         'optimalign' 'string'  { 'on' 'off' }         'on';
                         'delim'     {'integer' 'string'}   []         char([9 32 44])}, 'importevent');
if isstr(g), error(g); end;
if ~isempty(g.indices), g.append = 'yes'; end;
g.delim = char(g.delim);

% call from pop_importevent
% -------------------------
if ~isempty(g.event)
    event = g.event;
    if ~isstr(event)
        if size(event,2) > size(event,1), event = event'; end;
    end;
end;
g.event = event;

% determine latency for old event alignment
% -----------------------------------------
tmpalign    = g.align;
g.align     = [];
g.align.val = tmpalign;
if ~isnan(g.align.val)
    if isempty(oldevent)
        error('Setevent: no pre-existing event, cannot perform alignment');
    end;
    if ~isfield(oldevent, 'latency')
        error('Setevent: pre-existing events do not have a latency field for re-alignment');
    end;
    switch g.append
        % What is the reasoning behind this warning? It seems misleading. TobyF
        % 3/28/2007
        case {'yes' '''yes'''}, disp('Setevent warning: using align, events should not be appended but erased');
    end;
    if g.align.val < 0
        g.align.event = oldevent(1).latency;
    else
        g.align.event = oldevent(g.align.val+1).latency;
    end
    g.align.nbevent = length(oldevent);
    g.oldevents = oldevent;
    g.align.txt = sprintf([ 'Check alignment between pre-existing (old) and loaded event' ...
                          ' latencies:\nOld event latencies (10 first): %s ...\n' ], ...
                          int2str([ oldevent(1:min(10, length(oldevent))).latency ]));
else
    g.oldevents = [];
end;

tmpfields = fieldnames(g);
event = oldevent;
%% scan all the fields of g
% ------------------------
for curfield = tmpfields'
    if ~isempty(event), allfields = fieldnames(event);
    else                    allfields = {}; end;
    switch lower(curfield{1})
        case {'append', 'fields', 'skipline', 'indices', 'timeunit', 'align', 'delim' }, ; % do nothing now
        case 'event', % load an ascii file
            switch g.append
                case { '''no''' 'no' } % ''no'' for backward compatibility
                      event = load_file_or_array( g.event, g.skipline, g.delim );
                      allfields = g.fields(1:min(length(g.fields), size(event,2)));
                      event = eeg_eventformat(event, 'struct', allfields);
					  % generate ori fields
					  % -------------------
                      if ~isnan(g.timeunit)
                          for index = 1:length(event)
                              event(index).init_index = index;
                              event(index).init_time  = event(index).latency*g.timeunit;
                          end;
                      end;
                      event = recomputelatency( event, 1:length(event), srate, ...
                                                    g.timeunit, g.align, g.oldevents, g.optimalign);
                case { '''yes''' 'yes' }
                      % match existing fields
                      % ---------------------
                      tmparray = load_file_or_array( g.event, g.skipline, g.delim );
                      if isempty(g.indices) g.indices = [1:size(tmparray,1)] + length(event); end;
                      if length(g.indices) ~= size(tmparray,1)
                            error('Set error: number of row in file does not match the number of event given as input');
                      end;

                      % add field
                      % ---------
                      g.fields = getnewfields( g.fields, size(tmparray,2)-length(g.fields));

                      % add new values
                      % ---------------------
                      for eventfield = 1:size(tmparray,2)
                          event = setstruct( event, g.fields{eventfield}, g.indices, { tmparray{:,eventfield} });
                      end;
					  % generate ori fields
					  % -------------------
                      offset = length(event)-size(tmparray,1);
                      for index = 1:size(tmparray,1)
                          event(index+offset).init_index = index;
                          event(index+offset).init_time  = event(index+offset).latency*g.timeunit;
                      end;
                      event = recomputelatency( event, g.indices, srate, g.timeunit, ...
                                                    g.align, g.oldevents, g.optimalign);
            end;
      end;
end;

if isempty(event) % usefull 0xNB empty structure
    event = [];
end;

%% remove the events wit out-of-bound latencies
% --------------------------------------------
if isfield(event, 'latency')
    try
        res = cellfun('isempty', { event.latency });
        res = find(res);
        if ~isempty(res)
            fprintf( 'importevent warning: %d/%d event(s) have invalid latencies and were removed\n', ...
                     length(res), length(event));
            event( res ) = [];
        end;
    end;
end;

%% interpret the variable name
% ---------------------------
function array = load_file_or_array( varname, skipline, delim );
    if isstr(varname) & exist(varname) == 2  % mean that it is a filename
                                             % --------------------------
        array = loadtxt( varname, 'skipline', skipline, 'delim', delim, 'blankcell', 'off' );

    else % variable in the global workspace
         % --------------------------
         if isstr(varname)
             array = evalin('base', varname);
             if ~iscell(array)
                 array = mattocell(array, ones(1, size(array,1)), ones(1, size(array,2)));
             end;
         else
             array = varname;
         end;
    end;
return;

%% update latency values
% ---------------------
function event = recomputelatency( event, indices, srate, timeunit, align, oldevents, optimalign);

    % update time unit
    % ----------------
    if ~isfield(event, 'latency'),
        if isfield(event, 'duration')
            error('A duration field cannot be defined if a latency field has not been defined');
        end;
        return;
    end;
    if ~isnan(timeunit)
        for index = indices
            event(index).latency  = event(index).latency*srate*timeunit;
            if isfield(event, 'duration')
                event(index).duration = event(index).duration*srate*timeunit;
            end;
        end;
    end;

    % alignment with old events
    % -------------------------
    if ~isnan( align.val )
        if align.val >= 0, alignlatency = event(1).latency;
        else               alignlatency = event(-align.val+1).latency;
        end;
        for index = indices
             event(index).latency = event(index).latency-alignlatency+align.event;
        end;
        if length(event) ~= align.nbevent
            disp('Setevent warning: the number of pre-existing events do not correspond to the number of events');
            disp('                  that were read, so their latencies may have been wrongly re-aligned');
        end;
        fprintf(align.txt);
        fprintf('New event latencies (10 first): %s ...\n', int2str(round([ event(1:min(10, length(event))).latency ])));
    end;
    if strcmpi(optimalign, 'on') & ~isempty(oldevents)
        newlat = [ event.latency     ];
        oldlat = [ oldevents.latency ];

        newlat = repmat(newlat, [length(oldlat) 1]);
        oldlat = repmat(oldlat', [1 size(newlat,2)]);
        if align.val >= 0
            newlat = newlat-newlat(1);
            oldlat = oldlat-oldlat(1+align.val);
        else
            newlat = newlat-newlat(1-align.val);
            oldlat = oldlat-oldlat(1);
        end;

        newfactor = fminsearch('eventalign',1,[],newlat, oldlat);
        fprintf('Best sampling rate ratio found is %1.7f. Below latencies after adjustment\n', newfactor);
        if newfactor > 1.01 | newfactor < 0.99
            disp('Difference is more than 1%, something is wrong; ignoring ratio');
            newfactor = 1;
        else
            difference1 = eventalign(1        ,newlat,oldlat);
            difference2 = eventalign(newfactor,newlat,oldlat);
            fprintf('The average difference before correction was %f sample points\n', difference1);
            fprintf('The average difference after correction is %f sample points\n', difference2);
        end;

        %diffarray = abs(newfactor*newlat-oldlat)';
        %[allmins poss] = min(diffarray);
        %figure; hist(allmins);
    else
        newfactor = 1;
    end;
    if ~isnan( align.val ) & newfactor ~= 1
        if align.val >= 0
            latfirstevent = event(1).latency;
        else
            latfirstevent = event(-align.val+1).latency;
        end;
        for index = setdiff(indices, 1)
            event(index).latency = round(event(index).latency-latfirstevent)*newfactor+latfirstevent;
        end;
        if ~isempty(oldevents)
            fprintf('Old event latencies (10 first): %s ...\n', int2str(round([ oldevents(1:min(10, length(oldevents))).latency ])));
            fprintf('New event latencies (10 first): %s ...\n', int2str(round([ event(1:min(10, length(event))).latency ])));
        end;
    else
        % must add one (because first sample point has latency 0
        % ------------------------------------------------------
        for index = indices
            event(index).latency = round((event(index).latency+1)*1000*newfactor)/1000;
        end;
    end;

%% create new field names
% ----------------------
function epochfield = getnewfields( epochfield, nbfields )
   count = 1;
   while nbfields > 0
       if isempty( strmatch([ 'var' int2str(count) ], epochfield ) )
               epochfield =  { epochfield{:} [ 'var' int2str(count) ] };
               nbfields = nbfields-1;
       else    count = count+1;
       end;
   end
return;

%%
% ----------------------
function var = setstruct( var, fieldname, indices, values )
    if exist('indices') ~= 1, indices = 1:length(var); end
    if ~isempty(values)
        for index = 1:length(indices)
            var = setfield(var, {indices(index)}, fieldname, values{index});
        end
    else
        for index = 1:length(indices)
            var = setfield(var, {indices(index)}, fieldname, '');
        end
    end
return;

author: PatternRecognition

repo: OpenBMI


text:
function Sk = mriphantom(kx,ky)
%========================================================================
% Sk = function mriphantom(kx,ky)
%
% Creates simulated raw MRI data of a Shepp Logan head phantom
% for given k space points
% input is the kx and ky coordinates in the image frequency space
% default values for kx and ky are for a Carthesian sampling.
% Input
%   kx, ky  Matrices with kx and ky sampling coordinates
%   E       Optional analytical phantom data matrix as produced by
%           image toolbox function phantom.m (see comments below).
% Data are from an analytical expression for a continuous phantom.
% Core equations for this are taken form :
%  Rik van der Walle et al IEEE Trans Med Imag 19(12) p1160.
%
% Note that the standard matlab function just gives the phantom
% sampled in image space on a carthesian grid.
% This function delivers the simulated raw MRI data sampled through
% a user defined path in k-space (image frequency space)
%
% Author : Ronald Ouwekerk Johns Hopkins University Dept. Radiology
% 601 N. Caroline Street Room 4250 Baltimore, MD 21287-0845 USA
% rouwerke@mri.jhu.edu
% Written :May 2002
% Modified June 2002 :  Revised to give a non-NaN result for k=0
%                       Suppress output if nargout == 0
%                       Added comments to explain the difference with phantom.m
%
% Whenever this software is used for publications an acknowledgment
% would be much appreciated.
%                           RO
%========================================================================
FOV = 2;
N = 64;

if nargin < 2
    %=====================================================================
    % Set the field of view to 24 cm (see image toolbox phantom.m)
    % resolution 128x128 points
    % create k-space coordinates for Carthesion sampling
    % Should be equivalent to IFFT of the phantom
    %=====================================================================
    FOV = 0.24;
    res = FOV/N;
    Kmax = 1/(2*res);
    ky = linspace(-Kmax,Kmax, N);
    ky = ky(ones(N,1),:);
    kx = permute(ky, [2,1]);
end;
if ( (nargin < 3 ) | isempty(E) )  & (exist('phantom.m')==2)
    %========================================================================
    % Create 2D phantom to get defining matrix E
    %========================================================================
    [P,E] = phantom;
elseif ~(exist('phantom.m')==2)
    P = [];
    E = modified_shepp_logan;
end;

%========================================================================
% copy the elements of E and scale to field of view
% Matlab phantom.m default yields FOV =[-1,1]= 2
% Scale all dimensions to desired FOV (*FOV/2).
%   Column 1:  rho  the additive intensity value of the ellipse
%   Column 2:  a    the length of the horizontal semi-axis of the ellipse
%   Column 3:  b    the length of the vertical semi-axis of the ellipse
%   Column 4:  x0   the x-coordinate of the center of the ellipse
%   Column 5:  y0   the y-coordinate of the center of the ellipse
%   Column 6:  alpha  the angle (in degrees) between the horizontal semi-axis
%                   of the ellipse and the x-axis of the image
%========================================================================
[me,ne] = size(E);
rho = E(:,1)';
nellipses = length(rho);
a   = FOV/2 * E(:,2)';
b   = FOV/2 * E(:,3)';
x0  = FOV/2 * E(:,4)';
y0  = FOV/2 * E(:,5)';
alpha = pi*E(:,6)'/180;

%========================================================================
% Stretch kspace coordinates to column vector (retain original size info)
%========================================================================
[mk,nk] = size(kx);
kx = kx(:);
ky = ky(:);
k = kx + j*ky;
klen = length(kx);

%========================================================================
% Calculate angle of the vector to the ellipse center in real space.
%========================================================================
gamma = atan2(y0, x0);
te = sqrt(x0.^2 + y0.^2);

%========================================================================
% Replicate all the ellipse defining vectors to the number of k space values
%========================================================================
gamma = gamma(ones(klen,1),:);
te = te(ones(klen,1),:);
alpha = alpha(ones(klen,1),:);
rho = rho(ones(klen,1),:);
a = a(ones(klen,1),:);
b = b(ones(klen,1),:);

%========================================================================
% Calculate k space vectors angle and magnitude
%========================================================================
kphi = angle(k);
kr = abs(k);

%========================================================================
% Replicate k space and time vectors to the number of ellipses
%========================================================================
kr = kr(:,ones(1,nellipses));
kphi = kphi(:,ones(1,nellipses));

%========================================================================
% Precalculate some values
%========================================================================
pi2 = 2*pi;
cose = cos(gamma - kphi);

%========================================================================
% Projection angles:
% Difference of the k space vector angle and the angles of the ellipses
%========================================================================
sind = sin(kphi - alpha);
cosd = cos(kphi - alpha);

akphi = sqrt(a.^2 .* cosd.^2 + b.^2 .* sind.^2);
%========================================================================
% Calculate the signals for each ellips and each k space point kx,ky
% avoid division by zero    besselj(1,0) = 0 so set 0/0 = 1
%========================================================================
knz = (akphi.*kr ~=0);
bess = ones(size(akphi.*kr));
bess(knz) = besselj(1,pi2.*akphi(knz).*kr(knz)) ./(akphi(knz).*kr(knz));
Sk = exp(-j*pi2.*kr.*te.*cose) .*rho.*a.*b.* bess;
%========================================================================
% Sum the signals for all ellipses and reshape the data back to
% the matrix size of the kx and ky input data
%========================================================================
Sk = sum(Sk,2);
Sk = reshape(Sk, mk,nk);
%========================================================================
% Show the result, including the reconstruction if default Carthesian
% sampling was used
%========================================================================
if  nargin >= 2
    fh = figure;
    subplot(1,2,1);
    if ~isempty(P)
        imagesc(abs(P));
    end;
    title('Software head phantom')
    axis image

    subplot(1,2,2);
    imagesc(abs(Sk));
    title('Simulated raw data')
    return
else
    fh = figure;
    subplot(1,3,1);
    if ~isempty(P)
        imagesc(abs(P));
    end;
    axis image

    subplot(1,3,2);
    imagesc(abs(Sk));
    title('Simulated raw data')
    axis image

    subplot(1,3,3);
    I = fftshift(fft2(Sk));
    imagesc(abs(I));
    axis image
    title('2DFFT Reconstructed image')
end;
%========================================================================
% Suppress output if not asked for any
%========================================================================
if nargout < 1
    Sk = [];
end;

%========================================================================
% END
%========================================================================

%========================================================================
% local function Copied from function phantom.m
% to produce a result without the image toolbox
%========================================================================

%========================================================================
function toft=modified_shepp_logan
%
%   This head phantom is the same as the Shepp-Logan except
%   the intensities are changed to yield higher contrast in
%   the image.  Taken from Toft, 199-200.
%
%         A    a     b    x0    y0    phi
%        ---------------------------------
toft = [  1   .69   .92    0     0     0
        -.8  .6624 .8740   0  -.0184   0
        -.2  .1100 .3100  .22    0    -18
        -.2  .1600 .4100 -.22    0     18
         .1  .2100 .2500   0    .35    0
         .1  .0460 .0460   0    .1     0
         .1  .0460 .0460   0   -.1     0
         .1  .0460 .0230 -.08  -.605   0
         .1  .0230 .0230   0   -.606   0
         .1  .0230 .0460  .06  -.605   0   ];
%========================================================================
% end local function toft=modified_shepp_logan
%========================================================================

author: Sable

repo: mcbench-benchmarks


text:
function x = fgplvmOptimisePoint(model, x, y, display, iters)

% FGPLVMOPTIMISEPOINT Optimise the postion of a latent point.
% FORMAT
% DESC optimises the location of a single point in latent space
% given an initialisation and an observed data point. Note that it
% ignores any dynamics present in the model.
% ARG model : the model for which the point will be optimised.
% ARG x : the initialisation of the point in the latent space.
% ARG y : the observed data point for which the latent point is to
% be optimised.
% ARG display : whether or not to display the iterations of the
% optimisation (default: true)
% ARG iters : maximum number of iterations for the optimisation
% (default 2000).
% RETURN x : the optimised location in the latent space.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
%
% SEEALSO : fgplvmCreate, fgplvmOptimiseSequence, fgplvmPointObjective, fgplvmPointGradient

% FGPLVM

if nargin < 5
  iters = 2000;
  if nargin < 4
    display = true;
  end
end

options = optOptions;
if display
  options(1) = 1;
  options(9) = 1;
end
options(14) = iters;

if isfield(model, 'optimiser')
  optim = str2func(model.optimiser);
else
  optim = str2func('scg');
end

if strcmp(func2str(optim), 'optimiMinimize')
  % Carl Rasmussen's minimize function
  x = optim('fgplvmPointObjectiveGradient', x, options, model, y);
else
  % NETLAB style optimization.
  x = optim('fgplvmPointObjective', x,  options, ...
            'fgplvmPointGradient', model, y);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [matchedPoints1,matchedPoints2]=matchFeaturePoints(I1,I2,thresh)
imagePoints1 = detectMinEigenFeatures(rgb2gray(I1), 'MinQuality', thresh);
tracker = vision.PointTracker('MaxBidirectionalError', 1, 'NumPyramidLevels', 5);
imagePoints1 = imagePoints1.Location;
initialize(tracker, imagePoints1, I1);
[imagePoints2, validIdx] = step(tracker, I2);
matchedPoints1 = imagePoints1(validIdx, :);
matchedPoints2 = imagePoints2(validIdx, :);
end

author: yihui-he


text:
function [avg] = read_ns_avg(filename)

% READ_NS_AVG read a NeuroScan 3.x or 4.x AVG File
%
% [avg] = read_ns_avg(filename)
%
%  The output data structure avg has the fields:
%   avg.data        - ERP signal in uV (Nchan x Npnt)
%   avg.nsweeps     - number of accepted trials/sweeps in avg
%   avg.variance    - variance of the signal (Nchan x Npnt)
%   avg.label       - electrode labels
%   avg.nchan       - number of channels
%   avg.npnt        - number of samplepoints in ERP waveform
%   avg.rate        - sample rate (Hz)
%   avg.time        - time for each sample OR
%   avg.frequency   - frequency for each sample
%   hdr.domain      - flag indicating time (0) or frequency (1) domain
%   avg.xmin        - prestimulus epoch start (e.g., -100 msec)
%   avg.xmax        - poststimulus epoch end (e.g., 900 msec)

% Copyright (C) 2002, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_ns_avg.m 945 2010-04-21 17:41:20Z roboos $

% read the neuroscan header
avg = read_ns_hdr(filename);

% create a time or frequency axis
if avg.domain==1
  avg.frequency = linspace(avg.xmin, avg.xmax, avg.npnt) / 1000;    % in Hz instead of mili-Hz
else
  avg.time = linspace(avg.xmin, avg.xmax, avg.npnt);            % in ms
end

% open the file and seek towards the place where the raw data is
fid = fopen(filename,'r','ieee-le');
if fid<0
  error(['cannot open ', filename]);
else
  fseek(fid, 900, 'cof');       % skip general header
  fseek(fid, 75*avg.nchan, 'cof');  % skip channel headers
end;

% read raw signal data and convert to uV
avg.data = zeros(avg.nchan, avg.npnt);
for elec = 1:avg.nchan
    fseek(fid, 5, 'cof'); % skip sweeps header
    raw = fread(fid, avg.npnt, 'float32');
    avg.data(elec,:) = (raw' - avg.baseline(elec)) * avg.calib(elec) / avg.nsweeps;
end;

% read signal variance if present
if avg.variance
    variance = zeros(avg.npnt, avg.nchan);
    for elec = 1:avg.nchan,
        variance(:, elec) = fread(fid, avg.npnt, 'float32');
    end;
    avg.variance = variance';
else
    avg.variance = [];
end;

fclose(fid);


text:
filename='Cantilever_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function plotNonlinCon(prob,data,linear)
%PLOTNONLINCON Plot Nonlinear Constraints on the current figure
%   plotNonlinCon(prob)

%   Copyright (C) 2013 Jonathan Currie (I2C2)

if(nargin < 3), linear = false; end

xl = xlim; yl = ylim;
hold on;

%Colour
dkg = [0.2 0.2 0.2];
data.linear = linear;

%Ensure we have row format
if(isempty(prob.cl))
    prob = nmix2row(prob);
end

%Determine number of nonlinear constraints
no = length(prob.cl);
%Generate Constraint Surface Points
[x1,x2] = meshgrid(linspace(xl(1),xl(2),data.npts),linspace(yl(1),yl(2),data.npts));
nox = size(x1); noy = size(x2);
%Form Constraint Function
if(prob.sizes.ndec==1)
    con = @(x) prob.nlcon(x(1));
else
    con = prob.nlcon;
end

%Create surface of all constraints
obj = zeros(nox(1),noy(2),no);
for n = 1:nox(1)
    for m = 1:noy(2)
        x = data.fixval;
        if(data.ndec > 1)
            x(data.idx) = [x1(n,m) x2(n,m)]';
        else
            x(data.idx) = x1(n,m);
        end
        obj(n,m,:) = con(x);
    end
end

%Plot All General Nonlinear Constraints
plotNLCon(con,obj,prob.cl,prob.cu,x1,x2,dkg,data);

hold off;

% OLD CODE
% if(isempty(prob.nle))
%     prob = nrow2mix(prob,0);
% end
%
% %Generate Nonlinear Constraint Contours & Plot
% nocon = length(prob.nle);
% if(nocon)
%     [x1,x2] = meshgrid(linspace(xl(1),xl(2),npts),linspace(yl(1),yl(2),npts));
%     nox = size(x1);
%     noy = size(x2);
%     obj = zeros(nox(1),noy(2),nocon);
%     % create surface of all constraints
%     for n = 1:nox(1)
%         for m = 1:noy(2)
%             x = [x1(n,m) x2(n,m)]';
%             obj(n,m,:) = prob.nlcon(x);
%         end
%     end
%     %Plot them
%     for i = 1:nocon
%         %Inequality
%         if(prob.nle(i))
%             c = contour(x1,x2,obj(:,:,i),'color',dkg,'levellist',prob.nlrhs(i));
%             %Plot Hatch
%             if(~isempty(c))
%                 %See if we have multiple contours
%                 len = size(c,2)-1;
%                 if(c(2,1) ~= len)
%                     %Build contour array
%                     cstrt = 2; cend = []; n = 2; ind = 1;
%                     while(ind <= len)
%                         ind = ind + c(2,ind) + 1;
%                         cend(n-1) = ind-1; %#ok<AGROW>
%                         cstrt(n) = ind+1;  %#ok<AGROW>
%                         n = n + 1;
%                     end
%                 else
%                     cstrt = 2;
%                     cend = len;
%                 end
%                 %Plot each contour hatch
%                 for n = 1:length(cend)
%                     %Get contour vectors
%                     vecx = diff(c(1,cstrt(n):cend(n)));
%                     vecy = diff(c(2,cstrt(n):cend(n)));
%                     if(isempty(vecx) || isempty(vecy))
%                         continue;
%                     end
%                     %Rotate hatch lines based on infeasible region
%                     xt = [c(1,cstrt(n))+vecy(1) c(2,cstrt(n))-vecx(1)]'; %check rotated -90
%                     fval = prob.nlcon(xt);
%                     if((fval(i) <= prob.nlrhs(i) && prob.nle(i) == -1) || (fval(i) >= prob.nlrhs(i) && prob.nle(i) == 1)) %rotate 90
%                         hvecx = -vecy;
%                         hvecy = vecx;
%                     else %rotate -90
%                         hvecx = vecy;
%                         hvecy = -vecx;
%                     end
%                     %Normalize
%                     av = mean(sqrt(hvecx.^2 + hvecy.^2));
%                     dirs = atan2(hvecy,hvecx);
%                     hvecx = av*cos(dirs);
%                     hvecy = av*sin(dirs);
%                     %Shift origin
%                     hvecx = c(1,cstrt(n):cend(n)-1) + hvecx;
%                     hvecy = c(2,cstrt(n):cend(n)-1) + hvecy;
%                     %Plot
%                     line([c(1,cstrt(n):cend(n)-1)' hvecx']',[c(2,cstrt(n):cend(n)-1)' hvecy']','Color','k')
%                 end
%             else
%                 optiwarn('opti:plot','Cannot plot inequality constraint as contour data is empty!');
%             end
%         else
%             contour(x1,x2,obj(:,:,i),'color',[0 0 1],'levellist',prob.nlrhs(i));
%         end
%     end
% end
%
% hold off;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function  p2 = project_3D_to_2D(p3)
d = distmatrix(p3);
p2_0 = p3(:,1:2);
fun = @(x)myfunc(x,d);

warnid = {'optim:fsolve:NonSquareSystem'};
for ii = 1:length(warnid)
    warning('OFF', warnid{ii});
end

p2 = fsolve(fun, p2_0, optimset('Display','off'));

for ii = 1:length(warnid)
    warning('ON', warnid{ii});
end

% -------------------------------------------------------
function  F = myfunc(p, d)
F = [];
kk = 1;
c = 1.2;
for ii = 1:size(d,1)
    for jj = 1:size(d,2)
        if d(ii,jj) == 0
            continue
        end
        F(kk) = sqrt((p(ii,1) - p(jj,1))^2 + (p(ii,2) - p(jj,2))^2) - d(ii,jj)*c;
        kk = kk+1;
    end
end

% -------------------------------------------------------------------
function d = distmatrix(varargin)

d=[];

v1 = double(varargin{1});
if size(v1,2) == 2
    v1 = [v1, zeros(size(v1,1),1)];
end

if length(varargin)==1

    n = size(v1,1);
    d = zeros(n);
    for i = 1:n
        for j = i+1:n
            d(i,j) = ((v1(i,1) - v1(j,1))^2 + ...
                      (v1(i,2) - v1(j,2))^2 + ...
                      (v1(i,3) - v1(j,3))^2)^0.5;
        end
    end

elseif length(varargin)==2

    v2 = double(varargin{2});
    m = size(v1,1);
    n = size(v2,1);
    d = zeros(m,n);
    for i = 1:m
        for j=1:n
            d(i,j) = ((v1(i,1) - v2(j,1))^2 + ...
                      (v1(i,2) - v2(j,2))^2 + ...
                      (v1(i,3) - v2(j,3))^2)^0.5;
        end
    end

end

author: BUNPC


text:
function lambda = parlett_eigenvalues ( )

%*****************************************************************************80
%
%% PARLETT_EIGENVALUES returns the eigenvalues of the PARLETT matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real LAMBDA(100,1), the eigenvalues.
%
  n = 100;
  lambda = zeros ( n, 1 );

  for i = 1 : n
    lambda(i,1) = i;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Slopes of temperature distribution with altitude.
m12=-6.5e-3;
m34=9.92e-4;
m45=2.78e-3;
m67=-1.96e-3;
m78=-3.94e-3;


text:
classdef PTKImageFileFormat
    % PTKImageFileFormat. An enumeration used to specify how a medical image is
    % stored
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    enumeration
        Dicom       (MimImageFileFormat.Dicom),      % DICOM format
        Metaheader  (MimImageFileFormat.Metaheader), % Metaheader (mha/mhd) plus raw data
        Matlab      (MimImageFileFormat.Matlab),     % Matlab matrix
        Analyze     (MimImageFileFormat.Analyze)     % Analyze format
    end

    properties (SetAccess = immutable)
        MimImageFileFormat
    end

    methods

        function obj = PTKImageFileFormat(mim_format)
            obj.MimImageFileFormat = mim_format;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
addpath(genpath('D:\GitHub\KiloSort2')) % path to kilosort folder
addpath('D:\GitHub\npy-matlab')

pathToYourConfigFile = 'D:\GitHub\KiloSort2\configFiles'; % take from Github folder and put it somewhere else (together with the master_file)
run(fullfile(pathToYourConfigFile, 'configFile384.m'))

ops.chanMap  = 'D:\GitHub\KiloSort2\configFiles\neuropixPhase3B2_kilosortChanMap.mat';

% common options for every probe
ops.NchanTOT    = 384;% total number of channels in your recording
ops.trange      =[0 Inf];% TIME RANGE IN SECONDS TO PROCESS

% find the binary file in this folder
rootZ = 'H:\DATA\Spikes\Allen\recording1';
fs = [dir(fullfile(rootZ, '*.dat')) dir(fullfile(rootZ, '*.bin'))];
fname = fs(1).name;
ops.fbinary     = fullfile(rootZ,  fname);

% path to whitened, filtered proc file (on a fast SSD)
rootH = 'H:\DATA\Spikes\temp\';
ops.fproc       = fullfile(rootH, 'temp_wh.dat'); % proc file on a fast SSD

% preprocess data to create temp_wh.dat
rez = preprocessDataSub(ops);

%%
% pre-clustering to re-order batches by depth
% rez = clusterSingleBatches(rez);
clusterSingleBatches2;

%%
% main optimization
rez = learnAndSolve8(rez);

% this does splits
rez    = splitAllClusters(rez);

% this saves to Phy
rezToPhy(rez, rootZ);

fname = fullfile(rootZ, 'rez.mat');
save(fname, 'rez', '-v7.3');

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function display(r)
% called to display a robot object

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

	disp(' ');
	disp([inputname(1), ' = '])
	disp(' ');
	disp(char(r))
	disp(' ');

author: petercorke

repo: robotics-toolbox-matlab