%% Script to import and analyse data from jsPsych online experiment

clc; clear;
set(0,'DefaultFigureWindowStyle','docked');

% subno = input('Subject No');

%% Import data from text file
% Script for importing data from the following text file:
%
%    filename: /Users/riccihannah/Google Drive/Online_Stop_Signal/jspsych-psychophysics-master/Data/All_Stop_Data_xp8q_smalltarg_slowest.csv
%
% Auto-generated by MATLAB on 16-Mar-2021 14:33:25

%% Setup the Import Options
opts = delimitedTextImportOptions("NumVariables", 40);

% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";

% Specify column names and types
opts.VariableNames = ["rt", "response_type", "key_press", "avg_frame_time", "center_x", "center_y", "CueType", "trial_type", "trial_index", "time_elapsed", "internal_node_id", "subject", "browser_name", "browser_version", "os_name", "os_version", "screen_resolution", "window_resolution", "click_x", "click_y", "CueDirection", "TrialType", "Staircase", "BlockName", "BlockNo", "Xpos", "Ypos", "IsHome", "ClickX", "ClickY", "TooSoon", "Jitter", "ReactTime", "MoveTime", "TargetHit", "MouseX", "MouseY", "MouseTime", "Outcome", "StopSignalDelay"];
opts.VariableTypes = ["categorical", "categorical", "categorical", "double", "double", "double", "categorical", "categorical", "double", "double", "string", "double", "categorical", "double", "categorical", "categorical", "categorical", "categorical", "double", "double", "categorical", "categorical", "categorical", "categorical", "double", "double", "double", "categorical", "double", "double", "categorical", "double", "double", "double", "categorical", "double", "double", "double", "categorical", "double"];
opts = setvaropts(opts, 11, "WhitespaceRule", "preserve");
opts = setvaropts(opts, [12, 36, 37, 38], "TrimNonNumeric", true);
opts = setvaropts(opts, [12, 36, 37, 38], "ThousandsSeparator", ",");
opts = setvaropts(opts, [1, 2, 3, 7, 8, 11, 13, 15, 16, 17, 18, 21, 22, 23, 24, 28, 31, 35, 39], "EmptyFieldRule", "auto");
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% Import the data
dataTable = readtable("/Users/riccihannah/Google Drive/Online_Stop_Signal/jspsych-psychophysics-master/Data/All_Stop_Data_xp8q_smalltarg_slowest.csv", opts);


%% Clear temporary variables
clear opts

%% Index trials
% Go trials
trl_GoAll = dataTable.BlockName == 'MainStopExp' & dataTable.TrialType == 'Go' & dataTable.TargetHit == 'Left' | ...
    dataTable.BlockName == 'MainStopExp' & dataTable.TrialType == 'Go' & dataTable.TargetHit == 'Right';
trl_GoCorr =  dataTable.BlockName == 'MainStopExp' & dataTable.Outcome == 'Correct_Go' & dataTable.TooSoon == 'false';
trl_GoErr =  dataTable.BlockName == 'MainStopExp' & dataTable.TrialType == 'Go' & dataTable.Outcome == 'Undefined'; %% Check this works - didn't have any to check

% Stop trials with the Homepad staircase
trl_StopAll_HomeStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'HomeStair' & ~isnan(dataTable.StopSignalDelay);
trl_StopCorr_HomeStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'HomeStair' & dataTable.Outcome == 'Correct_Stop_Home';
trl_StopFail_HomeStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'HomeStair' & dataTable.Outcome == 'Failed_Stop_Home' | ...
    dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'HomeStair' & dataTable.Outcome == 'Failed_Stop_Target';

% Stop trials with the Target staircase
trl_StopAll_TargetStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'TargetStair' & ~isnan(dataTable.StopSignalDelay);
trl_StopCorr_TargetStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Correct_Stop_Home' | ...
    dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Failed_Stop_Home';
trl_StopFail_TargetStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Failed_Stop_Target' ;
trl_StopFail_TargetStr_Home = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Failed_Stop_Home'; 

% How many Target staircase trials did they initiate a movement
numTrls_moved_TargetStr = sum(dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'TargetStair' & ~isnan(dataTable.ReactTime)); % looks like I initiated a movement every time!!

%% Get basic data
% All RTs and MTs
RT_GoAll = dataTable.ReactTime(trl_GoAll);
MT_GoAll = dataTable.MoveTime(trl_GoAll);

% Remove outliers from RT and MT
prctile_RT          = prctile(RT_GoAll,[25 75]);
prctile_MT          = prctile(MT_GoAll,[25 75]);

dataTable.ReactTime(dataTable.ReactTime < (prctile_RT(1) - 1.5 * iqr(RT_GoAll)) | ...
    dataTable.ReactTime >(prctile_RT(2)+ 1.5 * iqr(RT_GoAll)) | ...
    dataTable.ReactTime < 100) = NaN;

dataTable.MoveTime(dataTable.MoveTime < (prctile_MT(1) - 1.5 * iqr(RT_GoAll)) | ...
    dataTable.MoveTime>(prctile_MT(2) + 1.5 * iqr(RT_GoAll)) | ...
    dataTable.MoveTime < 100) = NaN;

RT_GoCorr = dataTable.ReactTime(trl_GoCorr);
RT_GoErr = dataTable.ReactTime(trl_GoErr);
MT_GoCorr = dataTable.MoveTime(trl_GoCorr);
MT_GoErr = dataTable.MoveTime(trl_GoErr);

RT_StopFail_HomeStair = dataTable.ReactTime(trl_StopFail_HomeStr);
MT_StopFail_HomeStair = dataTable.MoveTime(trl_StopFail_HomeStr);

RT_StopFail_TargetStair = dataTable.ReactTime(trl_StopFail_TargetStr);
MT_StopFail_TargetStair = dataTable.MoveTime(trl_StopFail_TargetStr);

% Combine RTs and MTs into one thing
TotRespT_GoCorr          = RT_GoCorr + MT_GoCorr;
TotRespT_StopFail_HomeStair  = RT_StopFail_HomeStair + MT_StopFail_HomeStair;
TotRespT_StopFail_TargetStair  = RT_StopFail_TargetStair + MT_StopFail_TargetStair;

% Mean RT and MT
RT_GoCorr_mn = nanmean(RT_GoCorr);
MT_GoCorr_mn = nanmean(MT_GoCorr);
RT_StopFail_HomeStair_mn = nanmean(RT_StopFail_HomeStair);
MT_StopFail_HomeStair_mn = nanmean(MT_StopFail_HomeStair(MT_StopFail_HomeStair > 0));
RT_StopFail_TargetStair_mn = nanmean(RT_StopFail_TargetStair);
MT_StopFail_TargetStair_mn = nanmean(MT_StopFail_TargetStair);

% All SSDs
SSD_Stop_HomeStair = dataTable.StopSignalDelay(trl_StopAll_HomeStr);
SSD_Stop_TargetStair = dataTable.StopSignalDelay(trl_StopAll_TargetStr);

% Mean SSD
SSD_Stop_HomeStair_mn = nanmean(SSD_Stop_HomeStair);
SSD_Stop_TargetStair_mn = nanmean(SSD_Stop_TargetStair);

% Stop success
pStop_HomeStair = (sum(trl_StopCorr_HomeStr) / sum(trl_StopAll_HomeStr)) *100;
pStop_TargetStair = (sum(trl_StopCorr_TargetStr) / sum(trl_StopAll_TargetStr)) *100;

%% Calculate SSRT for RTs and MTs
% RTs
if ~isempty(RT_GoErr)
    dataTable.ReactTime(trl_GoErr) = nanmax(RT_GoCorr);
    dataTable.MoveTime(trl_GoErr) = nanmax(MT_GoCorr);
    RT_GoAll = dataTable.ReactTime(trl_GoAll);% Max Go RT
    MT_GoAll = dataTable.MoveTime(trl_GoAll);
    TotRespT_GoCorr = RT_GoAll + MT_GoAll;
end
sort_RT = sort(RT_GoCorr(~isnan(RT_GoCorr)));
nth_RT_index_HomeStair  = round((pStop_HomeStair/100).*length(RT_GoCorr));
nth_RT_HomeStair        = sort_RT(nth_RT_index_HomeStair);
SSRT_HomeStair          = nth_RT_HomeStair - SSD_Stop_HomeStair_mn

% MTs
sort_TotRespT            = sort(TotRespT_GoCorr(~isnan(TotRespT_GoCorr)));
nth_RT_index_TargetStair = round((pStop_TargetStair/100).*length(sort_TotRespT));
nth_RT_TargeStair        = sort_TotRespT(nth_RT_index_TargetStair);
SSRT_TargetStair         = nth_RT_TargeStair - SSD_Stop_TargetStair_mn


%% Plot RTs and MTs
figure()

subplot(3,2,1);
histogram(RT_GoCorr,'BinWidth',25, 'Normalization','probability');
hold on
histogram(RT_StopFail_HomeStair,'BinWidth',25, 'Normalization','probability');
title('RTs: Home staircase');
xlabel('Reaction Times (ms)');
ylabel('Count');
xlim([0,700]);
% ylim([0,80]);
tx1 = (['Go RT: ' num2str(round(RT_GoCorr_mn)) 'ms']);
tx2 = (['FS RT: ' num2str(round(RT_StopFail_HomeStair_mn)) 'ms']);
text(40, 0.28,tx1);
text(40, 0.2,tx2);

subplot(3,2,2)
histogram(RT_GoCorr,'BinWidth',25, 'Normalization','probability');
hold on
histogram(RT_StopFail_TargetStair,'BinWidth',25, 'Normalization','probability');
title('RTs: Target staircase');
xlabel('Reaction Times (ms)');
ylabel('Count');
xlim([0,700]);
% ylim([0,80]);
tx1 = (['Go RT: ' num2str(round(RT_GoCorr_mn)) 'ms']);
tx2 = (['FS RT: ' num2str(round(RT_StopFail_TargetStair_mn)) 'ms']);
text(40, 0.28,tx1);
text(40, 0.2,tx2);;

subplot(3,2,3)
histogram(MT_GoCorr,'BinWidth',25, 'Normalization','probability')
hold on
histogram(MT_StopFail_HomeStair,'BinWidth',25, 'Normalization','probability');
title('MTs: Home staircase');
xlabel('Movement Times (ms)');
ylabel('Count');
xlim([0,500]);
% ylim([0,100]);
tx1 = (['Go MT: ' num2str(round(MT_GoCorr_mn)) 'ms']);
tx2 = (['FS MT: ' num2str(round(MT_StopFail_HomeStair_mn)) 'ms']);
text(20, 0.2,tx1);
text(20, 0.1,tx2);


subplot(3,2,4);
histogram(MT_GoCorr,'BinWidth',25, 'Normalization','probability');
hold on
histogram(MT_StopFail_TargetStair,'BinWidth',25, 'Normalization','probability');
title('MTs: Target staircase');
xlabel('Movement Times (ms)');
ylabel('Count');
xlim([0,500]);
% ylim([0,100]);
tx1 = (['Go MT: ' num2str(round(MT_GoCorr_mn)) 'ms']);
tx2 = (['FS MT: ' num2str(round(MT_StopFail_TargetStair_mn)) 'ms']);
text(20, 0.2,tx1);
text(20, 0.1,tx2);

subplot(3,2,5);
histogram(TotRespT_GoCorr,'BinWidth',25, 'Normalization','probability')
hold on
histogram(TotRespT_StopFail_HomeStair,'BinWidth',25, 'Normalization','probability')
title('Tot Time: Home staircase');
xlabel('Total Times (ms)');
ylabel('Count');
xlim([0,1000]);
% ylim([0,65]);
tx1 = (['Go TotT: ' num2str(round(nanmean(TotRespT_GoCorr))) 'ms']);
tx2 = (['FS TotT: ' num2str(round(nanmean(TotRespT_StopFail_HomeStair))) 'ms']);
text(40, 0.2,tx1);
text(40, 0.15,tx2);

subplot(3,2,6);
histogram(TotRespT_GoCorr,'BinWidth',25, 'Normalization','probability')
hold on
histogram(TotRespT_StopFail_TargetStair,'BinWidth',25, 'Normalization','probability')
% histogram(RT_StopFail_TargetStair,'BinWidth',25, 'Normalization','probability')
title('Tot Time: Target staircase');
xlabel('Total Times (ms)');
ylabel('Count');
xlim([0,1000]);
% ylim([0,65]);
tx1 = (['Go TotT: ' num2str(round(nanmean(TotRespT_GoCorr))) 'ms']);
tx2 = (['FS TotT: ' num2str(round(nanmean(TotRespT_StopFail_TargetStair))) 'ms']);
% tx3 = (['FS RT: ' num2str(round(RT_StopFail_TargetStair_mn)) 'ms']);
text(40, 0.28,tx1);
text(40, 0.2,tx2);
% text(40, 0.12,tx2);

print(gcf,strcat('Test3_smallT_slowest_RT_MT_Stop_Stop.pdf'),'-dpdf','-r400');

% play
% Do failed Target stops have shorter RT than successful target stops? 
% Yes, the success seems to be determined by the RTs 
% Made clearer by the fact that the MTs for Go and Failed stop are
% basically the same. That is not the decisive factor. 
% Why??? Harder to stop the further along in movement you are? 
% See that old paper - they said it might be easier, if you're already
% breaking. But I'm not breaking - shooting through target. 
trl_StopFail_TargetStr_Home = dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Failed_Stop_Home';
trl_StopFail_TargetStr_Target = dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Failed_Stop_Target';

RT_StopFail_TargetStr_Home = dataTable.ReactTime(trl_StopFail_TargetStr_Home);
RT_StopFail_TargetStair_Home_mn = nanmean(RT_StopFail_TargetStr_Home)

RT_StopFail_TargetStr_Target = dataTable.ReactTime(trl_StopFail_TargetStr_Target);
RT_StopFail_TargetStair_Target_mn = nanmean(RT_StopFail_TargetStr_Target)

figure()
bar([RT_StopFail_TargetStair_Home_mn,RT_StopFail_TargetStair_Target_mn] )


%% Plot stopping staircases
% Need to first calculate pStop at each SSD
% how many SSDs?
SSD_uniq_HomeStair = unique(SSD_Stop_HomeStair);
SSD_uniq_TargetStair = unique(SSD_Stop_TargetStair);
max_ntrls = max(length(SSD_uniq_HomeStair), length(SSD_uniq_TargetStair));

pStop_SSD = NaN(max_ntrls,2);

% Home stair case
for i = 1:length(SSD_uniq_HomeStair)
    currentSSD = SSD_uniq_HomeStair(i);
    %     count how many trials for that SSD in that stair case
    idx_all = dataTable.StopSignalDelay(trl_StopAll_HomeStr) == currentSSD;
    nTrls_all = sum(idx_all);
    idx_success = dataTable.Staircase == 'HomeStair' & dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Correct_Stop_Home';
    nTrls_correct = sum(idx_success);
    pStop_SSD(i,1 ) = (nTrls_correct / nTrls_all) * 100;
end

% Target stair case
for i = 1:length(SSD_uniq_TargetStair)
    currentSSD = SSD_uniq_TargetStair(i);
    %     count how many trials for that SSD in that stair case
    idx_all = dataTable.StopSignalDelay(trl_StopAll_TargetStr) == currentSSD;
    nTrls_all = sum(idx_all);
    idx_success = dataTable.Staircase == 'TargetStair' & dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Correct_Stop_Home' | ...
        dataTable.Staircase == 'TargetStair' & dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Failed_Stop_Home';
    nTrls_correct = sum(idx_success);
    pStop_SSD(i,2) = (nTrls_correct / nTrls_all) * 100;
end


idx_success_H = dataTable.Staircase == 'HomeStair' & ~isnan(dataTable.StopSignalDelay);
idx_success_T = dataTable.Staircase == 'TargetStair' & ~isnan(dataTable.StopSignalDelay);

figure()
plot(dataTable.StopSignalDelay(idx_success_H),  'b');
hold on
plot(dataTable.StopSignalDelay(idx_success_T),  'r');
lsline;
title('Staircasing','FontSize',14);
xlabel('Trial','FontSize',14);
ylabel('Stop signal delay (ms)','FontSize',14);
% xlim([0,max(SSD_uniq_TargetStair)]);
ylim([0,max(SSD_uniq_TargetStair)]);
legend('Home','Target','FontSize',14)
legend boxoff
print(gcf,strcat('Test3_smallT_slowest_StaircaseTracking_Stop.pdf'),'-dpdf','-r400');

figure()
scatter(SSD_uniq_HomeStair, pStop_SSD(1:length(SSD_uniq_HomeStair),1), 50, 'filled', 'b');
hold on
scatter(SSD_uniq_TargetStair, pStop_SSD(1:length(SSD_uniq_TargetStair),2), 50, 'filled', 'r');
plot([0, max(SSD_uniq_TargetStair)],[50,50],'--k');
lsline;
title('Stop functions','FontSize',14);
xlabel('Stop signal delay (ms)','FontSize',14);
ylabel('Probability of stopping (%)','FontSize',14);
xlim([0,max(SSD_uniq_TargetStair)]);
ylim([0,100]);
legend('Home','Target','50%','Hfit','Tfit','FontSize',14)
legend boxoff
print(gcf,strcat('Test3_smallT_slowestStopFunction_Stop.pdf'),'-dpdf','-r400');


%% Plot SSRT and mean SSD
figure()
subplot(3,2,1)
bar([RT_GoCorr_mn,0, 0], 'k');
hold on
bar([0,RT_StopFail_HomeStair_mn, 0],'b');
bar([0,0, RT_StopFail_TargetStair_mn], 'r');
er = errorbar([1,2,3],[RT_GoCorr_mn,RT_StopFail_HomeStair_mn,RT_StopFail_TargetStair_mn],...
    [nanstd(RT_GoCorr),nanstd(RT_StopFail_HomeStair), nanstd(RT_StopFail_TargetStair)],'k');
er.LineStyle = 'none';
title('React. Time','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2 3])
xticklabels({'Go','FS-H','FS-T'})
ylabel('RT (ms)','FontSize',14);
ylim([0,round(RT_GoCorr_mn+100)]);

subplot(3,2,2)
bar([MT_GoCorr_mn,0, 0], 'k');
hold on
bar([0,MT_StopFail_HomeStair_mn, 0],'b');
bar([0,0, MT_StopFail_TargetStair_mn], 'r');
er = errorbar([1,2,3],[MT_GoCorr_mn,MT_StopFail_HomeStair_mn,MT_StopFail_TargetStair_mn],...
    [nanstd(MT_GoCorr),nanstd(MT_StopFail_HomeStair), nanstd(MT_StopFail_TargetStair)],'k');
er.LineStyle = 'none';
title('Move. Time','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2 3])
xticklabels({'Go','FS-H','FS-T'})
ylabel('MT (ms)','FontSize',14);
ylim([0,round(MT_GoCorr_mn+50)]);

subplot(3,2,3)
bar([nanmean(TotRespT_GoCorr),0, 0], 'k');
hold on
bar([0,nanmean(TotRespT_StopFail_HomeStair), 0],'b');
bar([0,0, nanmean(TotRespT_StopFail_TargetStair)], 'r');
er = errorbar([1,2,3],[nanmean(TotRespT_GoCorr),nanmean(TotRespT_StopFail_HomeStair),nanmean(TotRespT_StopFail_TargetStair)],...
    [nanstd(TotRespT_GoCorr),nanstd(TotRespT_StopFail_HomeStair), nanstd(TotRespT_StopFail_TargetStair)],'k');
er.LineStyle = 'none';
title('Total Time','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2 3])
xticklabels({'Go','FS-H','FS-T'})
ylabel('Total Time (ms)','FontSize',14);
ylim([0,round(nanmean(TotRespT_GoCorr)+100)]);

subplot(3,2,4)
bar([pStop_HomeStair,0], 'b');
hold on
bar([0,pStop_TargetStair], 'r');
title('pStop','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2])
xticklabels({'Home','Target'})
ylabel('pStop (%)','FontSize',14);
ylim([0,100]);

subplot(3,2,5)
bar([SSD_Stop_HomeStair_mn,0], 'b');
hold on
bar([0,SSD_Stop_TargetStair_mn], 'r');
title('SSD','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2])
xticklabels({'Home','Target'})
ylabel('SSD (ms)','FontSize',14);
ylim([0,round(SSD_Stop_TargetStair_mn+20)]);

subplot(3,2,6)
bar([SSRT_HomeStair,0], 'b');
hold on
bar([0,SSRT_TargetStair], 'r');
title('SSRT','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2])
xticklabels({'Home','Target'})
ylabel('SSRT (ms)','FontSize',14);
ylim([0,round(SSRT_TargetStair+20)]);

print(gcf,strcat('Test3_smallT_slowestMeanBehav_Stop.pdf'),'-dpdf','-r400');





