%% 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_StopChange_Data_yoay.csv
%
% Auto-generated by MATLAB on 12-Mar-2021 09:23:40

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

% 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", "CueHeight", "ChangeHeight", "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", "categorical", "categorical", "double", "categorical", "categorical", "categorical", "categorical", "double", "double", "categorical", "categorical", "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, [38, 39, 40], "TrimNonNumeric", true);
opts = setvaropts(opts, [38, 39, 40], "ThousandsSeparator", ",");
opts = setvaropts(opts, [1, 2, 3, 7, 8, 11, 12, 13, 15, 16, 17, 18, 21, 22, 23, 24, 25, 26, 30, 33, 37, 41], "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_StopChange_Data_5h6q.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.TargetHit == 'Undefined'; %% Check this works - didn't have any to check

% Stop trials with the Nearpad staircase
trl_StopAll_NearStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'NearJump' & ~isnan(dataTable.StopSignalDelay);
trl_StopCorr_NearStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'NearJump' & dataTable.Outcome == 'Successful_Stop';
trl_StopFail_NearStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'NearJump' & dataTable.Outcome == 'Failed_Stop';

% Stop trials with the Far staircase
trl_StopAll_FarStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'FarJump' & ~isnan(dataTable.StopSignalDelay);
trl_StopCorr_FarStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'FarJump' & dataTable.Outcome == 'Successful_Stop';
trl_StopFail_FarStr = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'FarJump' & dataTable.Outcome == 'Failed_Stop';

% How many Far staircase trials did they initiate a movement
numTrls_moved_FarStr = sum(dataTable.Staircase == 'FarJump' & ~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);    
% Mean RT and MT
RT_GoCorr_mn = nanmean(RT_GoCorr);
MT_GoCorr_mn = nanmean(MT_GoCorr);
RT_StopFail_NearJump_mn = nanmean(RT_StopFail_NearJump);
MT_StopFail_NearJump_mn = nanmean(MT_StopFail_NearJump(MT_StopFail_NearJump > 0));
RT_StopFail_FarJump_mn = nanmean(RT_StopFail_FarJump);
MT_StopFail_FarJump_mn = nanmean(MT_StopFail_FarJump);

% All SSDs
SSD_Stop_NearJump = dataTable.StopSignalDelay(trl_StopAll_NearStr);
SSD_Stop_FarJump = dataTable.StopSignalDelay(trl_StopAll_FarStr);

% Mean SSD
SSD_Stop_NearJump_mn = nanmean(SSD_Stop_NearJump);
SSD_Stop_FarJump_mn = nanmean(SSD_Stop_FarJump);

% Stop success
pStop_NearJump = (sum(trl_StopCorr_NearStr) / sum(trl_StopAll_NearStr)) *100;
pStop_FarJump = (sum(trl_StopCorr_FarStr) / sum(trl_StopAll_FarStr)) *100;

%% Calculate SSRT for RTs and MTs
% Near
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_NearJump,'BinWidth',25, 'Normalization','probability');
title('RTs: Near 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_NearJump_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_FarJump,'BinWidth',25, 'Normalization','probability');
title('RTs: Far 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_FarJump_mn)) 'ms']);
text(40, 0.28,tx1);
text(40, 0.23,tx2);

subplot(3,2,3)
histogram(MT_GoCorr,'BinWidth',25, 'Normalization','probability')
hold on
histogram(MT_StopFail_NearJump,'BinWidth',25, 'Normalization','probability');
title('MTs: Near staircase');
xlabel('Movement Times (ms)');
ylabel('Count');
xlim([0,400]);
% ylim([0,100]);
tx1 = (['Go MT: ' num2str(round(MT_GoCorr_mn)) 'ms']);
tx2 = (['FS MT: ' num2str(round(MT_StopFail_NearJump_mn)) 'ms']);
text(20, 0.28,tx1);
text(20, 0.23,tx2);


subplot(3,2,4);
histogram(MT_GoCorr,'BinWidth',25, 'Normalization','probability');
hold on
histogram(MT_StopFail_FarJump,'BinWidth',25, 'Normalization','probability');
title('MTs: Far staircase');
xlabel('Movement Times (ms)');
ylabel('Count');
xlim([0,400]);
% ylim([0,100]);
tx1 = (['Go MT: ' num2str(round(MT_GoCorr_mn)) 'ms']);
tx2 = (['FS MT: ' num2str(round(MT_StopFail_FarJump_mn)) 'ms']);
text(20, 0.28,tx1);
text(20, 0.23,tx2);

subplot(3,2,5);
histogram(TotRespT_GoCorr,'BinWidth',20, 'Normalization','probability')
hold on
histogram(TotRespT_StopFail_NearJump,'BinWidth',20, 'Normalization','probability')
title('Tot Time: Near 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_NearJump))) 'ms']);
text(40, 0.28,tx1);
text(40, 0.23,tx2);

subplot(3,2,6);
histogram(TotRespT_GoCorr,'BinWidth',20, 'Normalization','probability')
hold on
histogram(TotRespT_StopFail_FarJump,'BinWidth',20, 'Normalization','probability')
histogram(RT_StopFail_FarJump,'BinWidth',20, 'Normalization','probability')
title('Tot Time: Far 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_FarJump))) 'ms']);
tx3 = (['FS RT: ' num2str(round(RT_StopFail_FarJump_mn)) 'ms']);
text(40, 0.28,tx1);
text(40, 0.23,tx2);
text(40, 0.18,tx3);

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

% play
% Do failed Far 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_FarStr_Near = dataTable.Staircase == 'FarJump' & dataTable.Outcome == 'Failed_Stop_Near';
trl_StopFail_FarStr_Far = dataTable.Staircase == 'FarJump' & dataTable.Outcome == 'Failed_Stop_Far';

RT_StopFail_FarStr_Near = dataTable.ReactTime(trl_StopFail_FarStr_Near);
RT_StopFail_FarJump_Near_mn = nanmean(RT_StopFail_FarStr_Near)

RT_StopFail_FarStr_Far = dataTable.ReactTime(trl_StopFail_FarStr_Far);
RT_StopFail_FarJump_Far_mn = nanmean(RT_StopFail_FarStr_Far)

% figure()
% bar([RT_StopFail_FarJump_Near_mn,RT_StopFail_FarJump_Far_mn] )


%% Plot stopping staircases
% Need to first calculate pStop at each SSD
% how many SSDs?
SSD_uniq_NearJump = unique(SSD_Stop_NearJump);
SSD_uniq_FarJump = unique(SSD_Stop_FarJump);
max_ntrls = max(length(SSD_uniq_NearJump), length(SSD_uniq_FarJump));

pStop_SSD = NaN(max_ntrls,2);

% Near stair case
for i = 1:length(SSD_uniq_NearJump)
    currentSSD = SSD_uniq_NearJump(i);
    %     count how many trials for that SSD in that stair case
    idx_all = dataTable.StopSignalDelay(trl_StopAll_NearStr) == currentSSD;
    nTrls_all = sum(idx_all);
    idx_success = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'NearJump' & ... 
    dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Successful_Stop';
    nTrls_correct = sum(idx_success);
    pStop_SSD(i,1 ) = (nTrls_correct / nTrls_all) * 100;
end

% Far stair case
for i = 1:length(SSD_uniq_FarJump)
    currentSSD = SSD_uniq_FarJump(i);
    %     count how many trials for that SSD in that stair case
    idx_all = dataTable.StopSignalDelay(trl_StopAll_FarStr) == currentSSD;
    nTrls_all = sum(idx_all);
    idx_success = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'FarJump' & ...
        dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Successful_Stop';
    nTrls_correct = sum(idx_success);
    pStop_SSD(i,2) = (nTrls_correct / nTrls_all) * 100;
end


idx_success_N = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'NearJump' & ~isnan(dataTable.StopSignalDelay);
idx_success_F = dataTable.BlockName == 'MainStopExp' & dataTable.Staircase == 'FarJump' & ~isnan(dataTable.StopSignalDelay);

figure()
plot(dataTable.StopSignalDelay(idx_success_N),  'b');
hold on
plot(dataTable.StopSignalDelay(idx_success_F),  'r');
lsline;
title('Staircasing','FontSize',14);
xlabel('Trial','FontSize',14);
ylabel('Stop signal delay (ms)','FontSize',14);
% xlim([0,max(SSD_uniq_FarJump)]);
ylim([0,max(SSD_uniq_FarJump)]);
legend('Near','Far','FontSize',14)
legend boxoff
print(gcf,strcat('Test2_StaircaseTracking_StopChange.pdf'),'-dpdf','-r400');

figure()
scatter(SSD_uniq_NearJump, pStop_SSD(1:length(SSD_uniq_NearJump),1), 50, 'filled', 'b');
hold on
scatter(SSD_uniq_FarJump, pStop_SSD(1:length(SSD_uniq_FarJump),2), 50, 'filled', 'r');
plot([0, max(SSD_uniq_FarJump)],[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_FarJump)]);
ylim([0,100]);
legend('Near','Far','50%','Hfit','Tfit','FontSize',14)
legend boxoff
print(gcf,strcat('Test2_StopFunction_StopChange.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_NearJump_mn, 0],'b');
% bar([0,0, RT_StopFail_FarJump_mn], 'r');
% er = errorbar([1,2,3],[RT_GoCorr_mn,RT_StopFail_NearJump_mn,RT_StopFail_FarJump_mn],...
%     [nanstd(RT_GoCorr),nanstd(RT_StopFail_NearJump), nanstd(RT_StopFail_FarJump)],'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_NearJump_mn, 0],'b');
% bar([0,0, MT_StopFail_FarJump_mn], 'r');
% er = errorbar([1,2,3],[MT_GoCorr_mn,MT_StopFail_NearJump_mn,MT_StopFail_FarJump_mn],...
%     [nanstd(MT_GoCorr),nanstd(MT_StopFail_NearJump), nanstd(MT_StopFail_FarJump)],'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_NearJump), 0],'b');
% bar([0,0, nanmean(TotRespT_StopFail_FarJump)], 'r');
% er = errorbar([1,2,3],[nanmean(TotRespT_GoCorr),nanmean(TotRespT_StopFail_NearJump),nanmean(TotRespT_StopFail_FarJump)],...
%     [nanstd(TotRespT_GoCorr),nanstd(TotRespT_StopFail_NearJump), nanstd(TotRespT_StopFail_FarJump)],'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,1,1)
bar([pStop_NearJump,0], 'b');
hold on
bar([0,pStop_FarJump], 'r');
title('pStop','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2])
xticklabels({'Near','Far'})
ylabel('pStop (%)','FontSize',14);
ylim([0,100]);

subplot(3,1,2)
bar([SSD_Stop_NearJump_mn,0], 'b');
hold on
bar([0,SSD_Stop_FarJump_mn], 'r');
title('SSD','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2])
xticklabels({'Near','Far'})
ylabel('SSD (ms)','FontSize',14);
ylim([0,round(SSD_Stop_FarJump_mn+20)]);

subplot(3,1,3)
bar([SSRT_NearJump,0], 'b');
hold on
bar([0,SSRT_FarJump], 'r');
title('SSRT','FontSize',14);
xlabel('Staircase','FontSize',14);
xticks([1 2])
xticklabels({'Near','Far'})
ylabel('SSRT (ms)','FontSize',14);
ylim([0,round(SSRT_FarJump+20)]);

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