IGT_Toolbox_4_ModelRecovery.m

%%%%%%%%%%%%%%%%%%%%%%% IGTtoolbox Recovery - Modeling %%%%%%%%%%%%%%%%%%%%%
% This script performs the model recovery analysis, which is a slow and
% computationally intensive step, as it requires to fit each of the
% candidate models to the data generated by each of the candidate models.
% The results of this analysis can be analyzed using the 4bis script.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Romain Ligneul / romain.ligneul@gmail.com / August 2018.
% v1.1

% clean the environment
clear all;close all;clc;

% ask user for directory (what should be analyzed)
model_dir = [uigetdir('IGTmodelfit', 'Select folder containing the models of interest') '/'];

% get all models
dummy_list = dir(model_dir)
for d = 1:length(dummy_list); model_list{d,1} = dummy_list(d).name; end; model_list(1:2)=[];

% select models to use for the recovery analysis
model_selection = listdlg('PromptString', 'Select the models you want to compare:', ...
    'ListString', model_list, 'ListSize', [300 200], 'InitialValue', 1:length(model_list));
model_list = model_list(model_selection);

% load relevant dataset (matching its folder name in IGTdata/)
data_dir = [uigetdir('IGTdata', 'Select folder corresponding to the modeling analysis') '/'];
load([data_dir 'IGTdata.mat']);

% add to path
addpath(genpath([pwd '/Tools/VBA/']));
addpath(genpath([pwd '/Tools/MODELS/']));
addpath(genpath([pwd '/Tools/OTHERS/']));

% create the input structure based on each simulated dataset and fit each
% model (resulting number of model fit: MxM)

for m=1:length(model_list)
    
    % load model and create directory for the recovery analysis
    load([model_dir model_list{m} '/fitted_model.mat'],'R', 'A')
    A.output.dir = [pwd '/IGTmodelrecovery/' A.output_name '/' model_list{m} '/'];
    mkdir(A.output.dir);
    
    % adjust option structure
    A.simulate_and_recover = 0;
    A.complete_save = -1;
    
    for s = 1:length(data)
        
        % subject info (numeric or string)
        try
            A.fit.subject(s,1) = data{s}.id;
        catch % if id is string
            A.fit.subject{s,1} = data{s}.id;
        end
        % if available, report condition
        try
            A.fit.condition(s,1) = data{s}.cond;
        end
        
        %  inputs
        A.fit.u{s}(1,:) = 0:A.fit.ntrials-1;                                % indices of the trials where evolution function should run.
        A.fit.u{s}(2,:) = [0, R.simulation.simulated_choices(s,1:end-1)];   % previous choice
        A.fit.u{s}(3,:) = R.simulation.simulated_gains(s,:);                % previous win
        A.fit.u{s}(4,:) = R.simulation.simulated_loss(s,:);                 % previous lose
        A.fit.u{s}(5,:) = 1:A.fit.ntrials;                                  % indices of the trials where observation function should run.
        A.fit.u{s}(6,:) = R.simulation.simulated_choices(s,:);
        % multinomial observation vector
        A.fit.y{s} = zeros(4,length(A.fit.u{s}));
        for t = 1:length(A.fit.u{s})
            A.fit.y{s}(A.fit.u{s}(6,t),t) = 1;
        end
        
    end
    
    %% run the fitting procedure
    % cluster output will only works if you have the qsub functions provided by
    %fieldtrip and a torque-compatible grid computing architecture.
    cluster_ouput = 'torque_output/';
%     mkdir(cluster_ouput);
    return_dir = pwd;
    for m = 1:length(A.fit.models)
        if A.cluster_run ==1
            try cd(cluster_ouput);end
            qsubfeval(A.fit.models{m}, A, data, 'memreq', 4*(1024^3), 'timreq', 600*60);
        else
            A.fit.models{m}(A,data);
        end
    end
    cd(return_dir);
    
end