%count_repaired
%%This script count the number of repaired volumes in each run. written by
%%Jin Wang 5/20/2019. It will print out the name of the run (the first column) the number of volumes being
%%replaced (the second column) and how many chunks of more than or equal to 6 consecutive volumes being
%%replaced (the third column) for each subject.

global CCN;
subjects = {};
root='/gpfs51/dors2/gpc/JamesBooth/JBooth-Lab/BDL/ELP/Jin_code/bids';
addpath(genpath('/gpfs51/dors2/gpc/JamesBooth/JBooth-Lab/BDL/ELP/Jin_code/2preprocessing'));
CCN.func_n='sub*.tsv';
CCN.ses='ses-5';
n=6; %number of consecutive volumes being replaced. no more than 6 consecutive volumes being repaired.
writefile='acc_file_ses-5.txt';


%%%%%%%%%%%%%%should not edit below%%%%%%%%%%%%%%%%%%%%%%
if isempty(subjects) %if user did not specify the subjects list, then it will read in all the data in bids
    cd(root);
    listing=dir(root);
    all_list=extractfield(listing,'name');
    index=strfind(all_list,'sub');
    idx=find(not(cellfun('isempty',index)));
    subjects=all_list(idx);
end

cd(root);
if exist(writefile)
    delete(writefile);
end
fid=fopen(writefile,'w');
hdr='subjects run_name num_repaired chunks cond1 acc1 rt1 cond2 acc2 rt2 cond3 acc3 rt3 cond4 acc4 rt4';
fprintf(fid, '%s', hdr);
fprintf(fid, '\n');
for i=1:length(subjects)
    func_p=[root '/' subjects{i}];
    func_f=expand_path([func_p '/[ses]/func/[func_n]']);
    for j=1:length(func_f)
        run_n=func_f{j};
        [run_p, run_name]=fileparts(run_n);
        %get the movement data from art_repair
        cd(run_n);
        fileid=fopen('art_repaired.txt');
        m=fscanf(fileid, '%f');
        [num_repaired, col]=size(m);
        N=n; %no more than 6 consecutive volumes repaired
        x=diff(m')==1;
        ii=strfind([0 x 0], [0 1]);
        jj=strfind([0 x 0], [1 0]);
        out=((jj-ii)>=N);
        if out==0
            chunks=0;
        else
            chunks=sum(out(:)==1);
        end
        
        %get accuracies and rt from events.tsv
        [~, run_name]=fileparts(func_f{j});
        run_name=[run_name(1:end-7), '_bold'];
        data=tdfread(func_f{j});
        conditions=unique(data.trial_type,'row');
        acc_by_condition=[];
        rt_by_condition=[];
        rt_allvalue=[];
        rt_all_correcttrials=data.calculated_RT(data.calculated_accuracy==1, :);
        for jj=1: size(rt_all_correcttrials,1)
            if ~contains(num2str(rt_all_correcttrials(jj,:)), 'n/a')
                if isa(rt_all_correcttrials(jj,:),'double')
                    rt_allvalue=[rt_allvalue; rt_all_correcttrials(jj,:)];
                else
                    rt_allvalue=[rt_allvalue; str2double(rt_all_correcttrials(jj,:))];
                end
            end
        end
        M=mean(rt_allvalue);
        SD3=3*std(rt_allvalue);
        min1=M-SD3;
        min2=0.25;
        mini=max([min1;min2]);
        maxi=M+SD3;
        for ii=1:size(conditions,1)
            thiscondition=conditions(ii,:);
            [~, lencond]=size(conditions); 
            acc_thiscond=data.calculated_accuracy((sum((data.trial_type==thiscondition)')==lencond)');
            rt_thiscond=data.calculated_RT(((sum((data.trial_type==thiscondition)')==lencond)'),:);
            rt_thiscond_new_count=0;
            rt_thiscond_new_value=[];
            for mm=1: size(rt_thiscond,1)
                if ischar(rt_thiscond(mm,:)) && ~contains(rt_thiscond(mm,:),'n/a')
                    cur_rt=str2double(rt_thiscond(mm,:));
                    if acc_thiscond(mm)==1 && cur_rt>mini && cur_rt<maxi
                        rt_thiscond_new_count=rt_thiscond_new_count+1;
                        rt_thiscond_new_value=[rt_thiscond_new_value; cur_rt];
                    end
                elseif ~ischar(rt_thiscond(mm,:))
                    cur_rt=rt_thiscond(mm,:);
                    if acc_thiscond(mm)==1 && cur_rt>mini && cur_rt<maxi
                        rt_thiscond_new_count=rt_thiscond_new_count+1;
                        rt_thiscond_new_value=[rt_thiscond_new_value; cur_rt];
                    end
                end
                
            end
            average_rt=sum(rt_thiscond_new_value)/rt_thiscond_new_count;
            average_acc=sum(acc_thiscond)/size(acc_thiscond,1);
            acc_by_condition=[acc_by_condition; average_acc];
            rt_by_condition=[rt_by_condition; average_rt];
        end
        
        %get the shifted dates of acquisition date for each run
        fname=[run_name '.json'];
        if exist(fname)
            val=jsondecode(fileread(fname));
            shifted_date_acq=val.ShiftedAquisitionDate;
        else
            shifted_data_acq='NaN';
        end
        
        %save all the values to txt
        fprintf(fid,'%s %s %d %d %s %.4f %.6f %s %.4f %.6f %s %.4f %.6f %s %.4f %.6f\n', ...
            subjects{i}, run_name, num_repaired, chunks, conditions(1,:), acc_by_condition(1,:), rt_by_condition(1,:), ...
            conditions(2,:), acc_by_condition(2,:), rt_by_condition(2,:), conditions(3,:), acc_by_condition(3,:), rt_by_condition(3,:),...
            conditions(4,:), acc_by_condition(4,:), rt_by_condition(4,:));
    end
end