Added paper wrapper.

Author: micahgallen

fig4_get_params.m

@@ -6,7 +6,6 @@
 
 %% get parameters for intero and extero
 
-clear all
 evparam =[0.5:.1:1];
 ivparam =[0.5:.1:1];
 
@@ -16,8 +15,8 @@
     for ev = 1:numel(evparam)
 
 
-                filestr = sprintf(['/MDPfiles/MDP_90vol_%dip_%dep.mat'],ivparam(iv)*100, evparam(ev)*100);
-                filename = [pwd filestr];
+                filestr = sprintf(['/MDP_files/MDP_90vol_%dip_%dep.mat'],ivparam(iv)*100, evparam(ev)*100);
+                filename = [datpath filestr];
                 load(filename, 'MDP')
 
                 param_mat_extero(iv, ev) = sum(calc_h_extero(MDP)); % exteroceptive confidence/entropy

figure_3_plots.m

@@ -5,10 +5,12 @@
 %% initialize 
 clear h1 h2 spiderprob1 spiderprob2 mdp1 mdp2 
 close all
-load MDP_30n_100t_healthy.mat
-load MDP_30n_100t_lesion.mat
 
-figdir = [pwd '\figures\'];
+load(fullfile(datpath, 'MDP_30n_100t_healthy.mat'), 'mdp1')
+load(fullfile(datpath, 'MDP_30n_100t_lesion.mat'), 'mdp2')
+
+
+figdir = [pwd '\Figures\'];
 
 if ~exist(figdir)
     mkdir(figdir)
@@ -61,7 +63,7 @@
 xlim([0 40])
 box off
 fname = [figdir 'FearExpectation.pdf'];
-print(fname, '-dpdf', '-r600')
+% print(fname, '-dpdf', '-r600') % uncomment to save figure to disk 
 
 %% figure 3A - binned beats
 
@@ -192,4 +194,4 @@
 set(gca, 'FontSize', 12)
 box off
 fname = [figdir 'cardiac_response.pdf'];
-print(fname, '-dpdf', '-r600')
+%print(fname, '-dpdf', '-r600') % uncomment to save figure to disk 

figure_3_simulations.m

@@ -1,19 +1,6 @@
 %% script to run the simulations in figure 3
 % runs MDP simulation with the first 15 trials as fixed inputs. 
 
-%% get canonical variables for each cardiac phase
-close all
-clear all
-
-% path to code repository - change this
-addpath('C:\Users\Micah Allen\Google Drive\FEP_Heartbeat\code\cbrewer')
-
-% path to SPM 12 - change this
-addpath('C:\Users\Micah Allen\Dropbox\toolboxes\MatlabToolboxes\SPM12_latest')
-
-% intialize and close spm - need MDP functions on path
-spm fmri -nogui
-spm quit
 
 %% run model 1 - healthy control 
 % healthy interoception
@@ -39,28 +26,37 @@
 input.T = 100; % specify number of total trials to simulate
 
 %% simulate 30 healthy subjects in a parfor loop
-nsubjects = 30;
-
-parfor n = 1:nsubjects
-
-
-
-tic
-fprintf('\nSubject %d/%d Started\n', n, nsubjects)
-mdp1(n)= run_mdp_simulation(input);
-t=toc;
-fprintf('\nSubject %d/%d Simulated in %02f seconds \n'...
-    , n,nsubjects, t)
-
-
-
+% if you do not have parallel computing toolbox, just change to a for loop
+% and it should run. It may take a while. 
+if ~exist(fullfile(datpath, 'MDP_30n_100t_healthy.mat'), 'file')    
+    nsubjects = 30;
+    
+    parfor n = 1:nsubjects
+        
+        
+        
+        tic
+        fprintf('\nSubject %d/%d Started\n', n, nsubjects)
+        mdp1(n)= run_mdp_simulation(input);
+        t=toc;
+        fprintf('\nSubject %d/%d Simulated in %02f seconds \n'...
+            , n,nsubjects, t)
+        
+        
+        
+    end
+    
+    save(fullfile(datpath, 'MDP_30n_100t_healthy.mat'), 'mdp1')
+    
+else
+    
+    load(fullfile(datpath, 'MDP_30n_100t_healthy.mat'), 'mdp1')
+    
 end
 
-save MDP_30n_100t_healthy mdp1
-
-%% optional plotting
-spm_figure('GetWin','Figure 1');
-spm_MDP_VB_trial(mdp1(end))
+%% optional plotting of state trajectories 
+% spm_figure('GetWin','Figure 1');
+% spm_MDP_VB_trial(mdp1(end))
 
 %% simulation #2 - lesioned interoceptive precision
 
@@ -80,6 +76,7 @@
 input.s = [repmat([1 2 3], 1, length(input.s)./3); input.s];
 input.T = 100;
 
+if ~exist(fullfile(datpath, 'MDP_30n_100t_lesion.mat'), 'file')   
 
 nsubjects = 30;
 
@@ -98,13 +95,21 @@
 
 end
 
+save(fullfile(datpath, 'MDP_30n_100t_lesion.mat'), 'mdp2')
+  
+else
+   
 
+load(fullfile(datpath, 'MDP_30n_100t_lesion.mat'), 'mdp2')
+    
+    
+end
 
-save MDP_30n_100t_lesion mdp2
+    
 
-%% optional plotting
-spm_figure('GetWin','Figure 1');
-spm_MDP_VB_trial(mdp2(end))
+%% optional plotting of state trajectories 
+% spm_figure('GetWin','Figure 1');
+% spm_MDP_VB_trial(mdp2(end))
 
 
 

figure_4_cardiac_cycle_plots.m

@@ -1,24 +1,7 @@
 %% wrapper script to run simulations and
 %% get canonical variables for each cardiac phase
 %% reproduces Figure 4B and 4D
-close all
-clear all
 
-% path to code repository - must be run from inside code repo dir
-addpath(genpath(pwd))
-% path to SPM 12 - change this
-addpath('C:\Users\Micah Allen\Dropbox\toolboxes\MatlabToolboxes\SPM12_latest')
-
-% intialize and close spm - need MDP functions on path
-spm fmri -nogui
-spm quit
-
-
-figdir = [pwd '\figures\'];
-
-if ~exist(figdir)
-    mkdir(figdir)
-end
 
 %% run model 1 - healthy extero 
 
@@ -72,7 +55,7 @@
 [h1_3, h2_3, h3_3] = split_h_cardiac_extero(MDP3);
 
 figure; 
-cmap = cbrewer('seq', 'PuRd', 3);
+cmap = cbrewer2('seq', 'PuRd', 3);
 y = [h1_1 h2_1 h3_1; h1_2 h2_2 h3_2; h1_3 h2_3 h3_3];%; h1_4 h2_4 h3_4];
 b = bar(y,'FaceColor','flat');
 
@@ -88,15 +71,15 @@
 title('Cardiac Cycle and Exteroceptive Uncertainty')
 
 fname = [figdir 'fig4b_barplot_cardiac_cycle_extero.pdf'];
-print(fname, '-dpdf', '-r600')
+%print(fname, '-dpdf', '-r600')
 
 %%  figure 3D
 [h1_1, h2_1, h3_1] = split_h_cardiac_intero(MDP1);
 [h1_2, h2_2, h3_2] = split_h_cardiac_intero(MDP2);
 [h1_3, h2_3, h3_3] = split_h_cardiac_intero(MDP3);
 
 figure; 
-cmap = cbrewer('seq', 'PuRd', 3);
+cmap = cbrewer2('seq', 'PuRd', 3);
 y = [h1_1 h2_1 h3_1; h1_2 h2_2 h3_2; h1_3 h2_3 h3_3];%; h1_4 h2_4 h3_4];
 b = bar(y,'FaceColor','flat');
 
@@ -112,4 +95,4 @@
 title('Cardiac Cycle and Interoceptive Uncertainty')
 
 fname = [figdir 'fig4D_barplot_cardiac_cycle_intero.pdf'];
-print(fname, '-dpdf', '-r600')
+%print(fname, '-dpdf', '-r600')

figure_4_plot_params.m

@@ -1,7 +1,7 @@
 %% initialize
 
 
-figdir = [pwd '\figures\'];
+figdir = [pwd '\Figures\'];
 
 if ~exist(figdir)
     mkdir(figdir)
@@ -23,7 +23,7 @@
 %% figure
 
 figure
-cmap = cbrewer('div', 'Spectral', 50);
+cmap = cbrewer2('div', 'Spectral', 50);
 
 s=surf(tickparam, tickparam, param_mat_in); zlabel(zlab), xlabel(xlab),ylabel(ylab); colormap(cmap)
 set(gca, 'XTickLabel', tickparam, 'XTick',tickparam, 'YTickLabel', tickparam, 'YTick', tickparam)
@@ -36,7 +36,7 @@
 
 fname = [figdir 'fig4A_extero.pdf'];
 title('Exteroceptive Uncertainty')
-print(fname, '-dpdf', '-r600');
+%print(fname, '-dpdf', '-r600');
 
 
 %% parameters - intero
@@ -53,7 +53,7 @@
 %% figure
 
 figure
-cmap = cbrewer('div', 'Spectral', 50);
+cmap = cbrewer2('div', 'Spectral', 50);
 
 s=surf(tickparam, tickparam, param_mat_in); zlabel(zlab), xlabel(xlab),ylabel(ylab); colormap(cmap)
 set(gca, 'XTickLabel', tickparam, 'XTick',tickparam, 'YTickLabel', tickparam, 'YTick', tickparam)
@@ -65,4 +65,4 @@
 pbaspect([1 1 1])
 title('Interoceptive Uncertainty')
 fname = [figdir 'fig4C_intero.pdf'];
-print(fname, '-dpdf', '-r600');
+%print(fname, '-dpdf', '-r600');

figure_4_plots_wrapper.m

@@ -1,7 +1,6 @@
 %% wrapper script to generate the plots in figure 4
 
 %% first, run the simulations for panels A & C
-% note this will take a *long* time
 
 figure_4_simulations_wrapper;
 

figure_4_simulations.m

@@ -0,0 +1,25 @@
+%%  script to run many simulations in parallel,
+% for figure 4A and 4C. 
+
+
+if ~isdir(fullfile(datpath, '/MDP_files'))
+    mkdir(fullfile(datpath, '/MDP_files'))
+    addpath(fullfile(datpath, '/MDP_files'))
+end
+
+%% this will run *many* simulations in parallel and
+% save them. Will take at least a day to run! 
+eparam =[0.5:0.1:1];
+iparam =[0.5:0.1:1];
+
+for ev = 1:numel(eparam)
+    
+    
+    parfor iv = 1:numel(iparam)
+        % run sim
+        MDP = sim_save(iparam(iv), eparam(ev),datpath);
+        
+        
+    end
+    
+end

figure_4_simulations_wrapper.m

@@ -1,37 +1,23 @@
-%% wrapper script to run many simulations in parallel,
+%%  script to run many simulations in parallel,
 % for figure 4A and 4C. 
 
-close all
-clear all
 
-% path to code repository - must be run from inside code repo dir
-addpath(genpath(pwd))
-% path to SPM 12 - change this
-addpath('C:\Users\Micah Allen\Dropbox\toolboxes\MatlabToolboxes\SPM12_latest')
-
-% intialize and close spm - need MDP functions on path
-spm fmri -nogui
-spm quit
-
-% parameters to loop over - here the full range from 
-% lesioned precision to hyper-precision in steps of 10%
-
-eparam =[0.5:0.1:1];
-iparam =[0.5:0.1:1];
-
-if ~isdir('MDP_files')
-    mkdir('MDP_files')
+if ~isdir(fullfile(datpath, '/MDP_files'))
+    mkdir(fullfile(datpath, '/MDP_files'))
+    addpath(fullfile(datpath, '/MDP_files'))
 end
 
 %% this will run *many* simulations in parallel and
 % save them. Will take at least a day to run! 
+eparam =[0.5:0.1:1];
+iparam =[0.5:0.1:1];
 
 for ev = 1:numel(eparam)
 
 
     parfor iv = 1:numel(iparam)
         % run sim
-        MDP = sim_save(iparam(iv), eparam(ev));
+        MDP = sim_save(iparam(iv), eparam(ev),datpath);
 
 
     end

figure_5_hrv_plots.m

@@ -1,5 +1,5 @@
 
-load MDP_hrv_models.mat
+
 
 
 %% healthy model
@@ -12,15 +12,15 @@
 plot_data(1).data = [summary_data.pLF summary_data.pHF]; % output data for final barplot
 
 hrv_rr_plot % rr 
-print([figdir 'rr_plot_healthy.pdf'], '-dpdf', '-r600')
+%print([figdir 'rr_plot_healthy.pdf'], '-dpdf', '-r600')
 
 hrv_trace_plot % trace
 
-print([figdir 'hrv_traceplot_healthy.pdf'], '-dpdf', '-r600')
+%print([figdir 'hrv_traceplot_healthy.pdf'], '-dpdf', '-r600')
 
 psd_plot_hrv % spectral frequency
 
-print([figdir 'hrv_spectral_healthy.pdf'], '-dpdf', '-r600')
+%print([figdir 'hrv_spectral_healthy.pdf'], '-dpdf', '-r600')
 
 
 %% visceral hyper precision model
@@ -34,7 +34,7 @@
 
 psd_plot_hrv % spectral frequency
 
-print([figdir 'hrv_spectral_hyper.pdf'], '-dpdf', '-r600')
+%print([figdir 'hrv_spectral_hyper.pdf'], '-dpdf', '-r600')
 
 %% visceral hyper arousal prior model
 
@@ -47,4 +47,4 @@
 
 psd_plot_hrv % spectral frequency
 
-print([figdir 'hrv_spectral_arousalprior.pdf'], '-dpdf', '-r600')
+%print([figdir 'hrv_spectral_arousalprior.pdf'], '-dpdf', '-r600')