clear all; close all;

addpath ../../Toolboxes/Palamedes


%% Figure 0

% experimental paradigm

%% Figure 1

% Which subjects?
subjects = [1:34];

% Which data?
whichdata = 'clean';  % alternative: 'raw'

if whichdata == 'clean'
    load('P:/3018029.09/data_clean.mat')
    folder_str = 'Clean data';
elseif whichdata == 'raw'
    load('P:/3018029.09/data.mat')
    folder_str = 'Raw data';
end

%% Figure 2


% Which stimulus levels?
stim_levels_all = table2array(unique(data(:,'relativeCoherence')));

%% a

% Calculate the proportion "more coherent" responses per stimulus level

for i = 1:length(stim_levels_all)
    
    data_now            = data(data.relativeCoherence == stim_levels_all(i), :);
    
    prop_more_coh(i)    = mean(data_now.responseDir);
    prop_more_coh_se(i) = std(data_now.responseDir)/sqrt(length(data_now.responseDir));
    
    responses(i) = sum(...
        data_now.responseDir == 1 & ...
        data_now.relativeCoherence == stim_levels_all(i));
    num(i) = size(...
        data_now.relativeCoherence(data_now.relativeCoherence == stim_levels_all(i)),1);
       
end

p = responses ./ num;

% psychometric function
PF = @PAL_CumulativeNormal;

%Threshold and Slope are free parameters, guess and lapse rate are fixed
paramsFree = [1 1 1 1];  %1: free parameter, 0: fixed parameter

%Parameter grid defining parameter space through which to perform a
%brute-force search for values to be used as initial guesses in iterative
%parameter search.
searchGrid.alpha = -0.2:.001:0.2;
searchGrid.beta = logspace(0,3,101);
searchGrid.gamma = [0:.01:.06];
searchGrid.lambda = [0:.01:.06];

%Perform fit
[fit, ~,~] = ...
    PAL_PFML_Fit(stim_levels_all',responses,num,searchGrid,paramsFree,PF,'gammaEQlambda',1,'lapseLimits',[0 0.1]);

% Save PSE into structure
PSE.prevchoice(i) = fit(1);

% Plot data and model
stimLevelsFineGrain = [min(stim_levels_all):max(stim_levels_all)./1000:max(stim_levels_all)];
y_model = PF(fit,stimLevelsFineGrain);

figure; 
subplot(1,3,1); hold on;
%    plot(stim_levels_all, prop_more_coh, (prop_more_coh-prop_more_coh_se), 'Color', [.6 .6 .6])
%    plot(stim_levels_all, prop_more_coh, (prop_more_coh+prop_more_coh_se), 'Color', [.6 .6 .6])
plot(stim_levels_all, prop_more_coh, '.k')
        plot(stimLevelsFineGrain,y_model,'k')
xlabel('Evidence')
ylabel('Proportion "more coherent" responses')

%% b & c

for iSubject = 1:length(subjects)
    
    % Select the subject's data
    
    subject_data = data(data.subject == subjects(iSubject), :);
    
    stim_levels_all = table2array(unique(subject_data(:,'relativeCoherence')));
    stim_levels     = unique(abs(stim_levels_all));
    
    for iStimulus = 1:length(stim_levels)
        
        
        % Select data for stimulus levels
        subject_data_stim = subject_data(subject_data.relativeCoherenceAbs == stim_levels(iStimulus), :);
        
        % Mean RT and confidence split for incorrect (i=1)/correct (i=2) trials
        for i = 1:2 
            responsetime_corr_split{i}(iSubject,iStimulus) = ...
                1000* nanmean(subject_data_stim.rt(~isnan(subject_data_stim.responseKey) & subject_data_stim.correct == i-1)); 
            confidence_corr_split{i}(iSubject,iStimulus) = ...
                nanmean(subject_data_stim.confidence(subject_data_stim.correct == i-1));  
        end 
        % Mean RT and confidence for all trials (i=3)
        responsetime_corr_split{3}(iSubject,iStimulus) = ...
            1000* nanmean(subject_data_stim.rt(~isnan(subject_data_stim.responseKey))); 
        confidence_corr_split{3}(iSubject,iStimulus) = ...
            nanmean(subject_data_stim.confidence);  
        
        % Calculate the number of trials for each participant for each data point
        for i = 1:2 
            responsetime_corr_split_n{i}(iSubject,iStimulus) = ...
                length(subject_data_stim.rt(~isnan(subject_data_stim.responseKey) & subject_data_stim.correct == i-1)); 
            confidence_corr_split_n{i}(iSubject,iStimulus) = ...
                length(subject_data_stim.confidence(subject_data_stim.correct == i-1));  
        end 
        responsetime_corr_split_n{3}(iSubject,iStimulus) = ...
            length(subject_data_stim.rt(~isnan(subject_data_stim.responseKey))); 
        confidence_corr_split_n{3}(iSubject,iStimulus) = ...
            length(subject_data_stim.confidence);  
               
        % Remove data points with fewer than 10 trials
        for i = 1:3
            if responsetime_corr_split_n{i}(iSubject,iStimulus) < 10
                responsetime_corr_split{i}(iSubject,iStimulus) = NaN;
            end
            if confidence_corr_split_n{i}(iSubject,iStimulus) < 10
                confidence_corr_split{i}(iSubject,iStimulus) = NaN;
            end
        end
        
    end
  
end

for i = 1:3
    
    % If there are fewer than 5 remaining participants in a column, remove
    % so that the point will not be plotted
    for iStimulus = 1:length(stim_levels)
        if sum(~isnan(responsetime_corr_split{i}(:,iStimulus))) < 10
            responsetime_corr_split{i}(:,iStimulus) = NaN;
        end
        if sum(~isnan(confidence_corr_split{i}(:,iStimulus))) < 10
            confidence_corr_split{i}(:,iStimulus) = NaN;
        end
    end
        
    % Calculate means and SEs    
    responsetime_corr_split_mean(i,:)   = nanmean(responsetime_corr_split{i});
        responsetime_corr_split_se(i,:) = nanstd(responsetime_corr_split{i})/sqrt(length(responsetime_corr_split{i}(~isnan(responsetime_corr_split{i}))));
    confidence_corr_split_mean(i,:)     = nanmean(confidence_corr_split{i});
        confidence_corr_split_se(i,:)   = nanstd(confidence_corr_split{i}/sqrt(length(confidence_corr_split{i}(~isnan(confidence_corr_split{i})))));
end

% Mean response time per stimulus level, with standard errors, split for incorrect/correct trials
subplot(1,3,2); hold on;
plot(responsetime_corr_split_mean');
    errorbar(responsetime_corr_split_mean', responsetime_corr_split_se', '.k')
legend('Incorrect responses', 'Correct responses', 'All responses')
legend('Location','northwest'); legend boxoff;
ylim([850 1600]); 
xlim([0.25 6.75]); 
set(gca,'XTick',1:6,'XTickLabel',{'smallest','','','','','largest'});
xlabel('|Evidence| strength')
ylabel('Response time (ms)')
%title('Mean response time per stimulus level');

% Mean confidence per stimulus level, split for incorrect/correct trials
subplot(1,3,3); hold on;
plot(confidence_corr_split_mean');
    errorbar(confidence_corr_split_mean', confidence_corr_split_se', '.k')
legend('Incorrect responses', 'Correct responses', 'All responses')
legend('Location','northwest'); legend boxoff;
ylim([1 4]); 
xlim([0.25 6.75]); 
set(gca,'XTick',1:6,'XTickLabel',{'smallest','','','','','largest'});
xlabel('|Evidence| strength')
ylabel('Confidence (1-4)')
%title('Mean confidence per stimulus level');

% Change figure size and save
set(gcf, 'Position',  [100, 100, 1300, 400])
saveas(gcf,['P:/3018029.09/Draft/Figure 2 after revision.png'])
saveas(gcf,['P:/3018029.09/Draft/Figure 2 after revision.eps'])

%% Figure 3

% Which subjects?
subjects = [1:34];
example_subject = 3;

% Which data?
whichdata = 'clean';  % alternative: 'raw'

if whichdata == 'clean'
    load('P:/3018029.09/data_sd_clean.mat')
    folder_str = 'Clean data';
elseif whichdata == 'raw'
    load('P:/3018029.09/data_sd.mat')
    folder_str = 'Raw data';
end
alldata = data_sd; clear data_sd

% Which stimulus levels?
stim_levels_all = table2array(unique(alldata(:,'currRelativeCoherence')));

% Append stimulus difficulty to data
difficulty_levels = [6 5 4 3 2 1 1 2 3 4 5 6];

prevDifficulty = difficulty_levels(alldata.prevStimLevel)';
currDifficulty = difficulty_levels(alldata.currStimLevel)';

alldata = [alldata, table(prevDifficulty), table(currDifficulty)];
    
%% a

% Select the subject's data
    
data = alldata(alldata.subject == example_subject, :);

stim_levels_all = table2array(unique(data(:,'currRelativeCoherence')));
stim_levels     = stim_levels_all((size(stim_levels_all)/2+1):size(stim_levels_all),:);

stimLevelsFineGrain=[min(stim_levels_all):max(stim_levels_all)./1000:max(stim_levels_all)];
       
% p_higher for prev_higher and prev_lower

figure; 
subplot(2,3,1); hold on;
title('Subject 3')

for i = 1:2
    for iCoherence = 1:length(stim_levels_all)
        responses_higher_split(i,iCoherence) = sum(...
            data.currResponseDir == 1 & ...
            data.currRelativeCoherence == stim_levels_all(iCoherence) & ...
            data.prevResponseDir == i-1);
        num_higher_split(i,iCoherence) = size(...
            data.currRelativeCoherence(data.currRelativeCoherence == stim_levels_all(iCoherence) & data.prevResponseDir == i-1),1);
    end
    p_higher_split(i,:) = responses_higher_split(i,:) ./ num_higher_split(i,:);

    % psychometric function
    PF = @PAL_CumulativeNormal;

    %Threshold and Slope are free parameters, guess and lapse rate are fixed
    paramsFree = [1 1 1 1];  %1: free parameter, 0: fixed parameter

    %Parameter grid defining parameter space through which to perform a
    %brute-force search for values to be used as initial guesses in iterative
    %parameter search.
    searchGrid.alpha = -0.2:.001:0.2;
    searchGrid.beta = logspace(0,3,101);
    searchGrid.gamma = [0:.01:.06];
    searchGrid.lambda = [0:.01:.06];

    %Perform fit
    [fit, ~,~] = ...
        PAL_PFML_Fit(stim_levels_all',responses_higher_split(i,:),num_higher_split(i,:),searchGrid,paramsFree,PF,'gammaEQlambda',1,'lapseLimits',[0 0.1]);

    % Save PSE into structure
    PSE.prevchoice(iSubject,i) = fit(1);

    % Plot data and model
    y_model = PF(fit,stimLevelsFineGrain);

    if i == 1 
        x = plot(stim_levels_all,p_higher_split(i,:),'.','Color',[0, 0.4470, 0.7410]);
        plot(stimLevelsFineGrain,y_model,'-','Color',[0, 0.4470, 0.7410])
        %text(stim_levels_all,p_higher_split(i,:),strsplit(num2str(num_higher_split(i,:))),'Fontsize',7)
    elseif i == 2
        y = plot(stim_levels_all,p_higher_split(i,:),'.','Color',[0.3010, 0.7450, 0.9330]);
        plot(stimLevelsFineGrain,y_model,'-','Color',[0.3010, 0.7450, 0.9330])
        %text(stim_levels_all,p_higher_split(i,:),strsplit(num2str(num_higher_split(i,:))),'Fontsize',7)
    end

end

ylim([0 1]); 
yticks([0 .25 .5 .75 1])
xlim([-0.30 0.30]); 
xlabel('Signed evidence strength')
ylabel('% "higher coherence" responses')
leg = legend([x,y],'"lower coherence"','"higher coherence"');   
title(leg,'prev. response')
legend('Location','southeast'); legend boxoff;


%% b & c

% RT

for rt = 1:2

for sj = unique(subjects)
    
    % individual sj data
    data = alldata(alldata.subject == subjects(sj), :);
        
    nbins = 6;
      
    for u = 1:nbins
                
        % first select trials for this bin
        trls = find(data.prevDifficulty == u);
                
        % find median rt
        datartsorted    = sortrows(data.prevRt(trls));
        medianrt        = datartsorted(floor(length(datartsorted)/2));
        
      % select low or high rts
        if rt == 1 
            trlssplit = trls(data.prevRt(trls) < medianrt);
        elseif rt == 2
            trlssplit = trls(data.prevRt(trls) >= medianrt);
        end
                
        % =========================================== %
        % history-dependent logistic
        % =========================================== %
        
        % previous response
        resps = [0 1];
        for r = 1:2
            
            % take a further subset of trials with previous response 1 or 0
            thisdat = data(trlssplit,:);
            thisdat = thisdat(thisdat.prevResponseDir == resps(r),:);
            
            [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                thisdat.currRelativeCoherence, thisdat.currResponseDir);
            
            % save betas
            grandavg{rt}.logisticHistory(sj, r, u, :) = [bias slope];
        end    
           
    end
          
end

% ========================================================= %
% transform history-dependent bias into repetition bias
% ========================================================= %

resp1 = -1 * squeeze(grandavg{rt}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
resp2 =      squeeze(grandavg{rt}.logisticHistory(:, 2, :, 1));

% since this is centred at 0.5, treat it that way
grandavg{rt}.repetition = (resp1 + resp2) ./ 2;

% ========================================================= %
% transform from log(odds) to probability
% ========================================================= %

logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
grandavg{rt}.repetition = logOdds2Probability(grandavg{rt}.repetition);

end

% save means and SEs for prev low and high rt

mean_repetition_lowrt = mean(grandavg{1}.repetition);
sem_repetition_lowrt = std(grandavg{1}.repetition)/sqrt(length(subjects));

mean_repetition_highrt = mean(grandavg{2}.repetition);
sem_repetition_highrt = std(grandavg{2}.repetition)/sqrt(length(subjects));

subplot(2,3,2); hold on;
%title('All subjects')
plot(mean_repetition_lowrt,'.','Color',[0.4660, 0.6740, 0.1880]);
    v = errorbar(mean_repetition_lowrt,sem_repetition_lowrt,'-','Color',[0.4660, 0.6740, 0.1880]);
plot(mean_repetition_highrt,'.','Color',[0.8500, 0.3250, 0.0980]);
    w = errorbar(mean_repetition_highrt,sem_repetition_highrt,'-','Color',[0.8500, 0.3250, 0.0980]);
xlim([.5 nbins+.5])
ylim([.3 .7]); 
set(gca,'XTick',1:7,'XTickLabel',{'smallest','','','','','largest'});
ylabel('P(repeat)')
xlabel('Previous |evidence| strength')

legend([v,w],'prev. fast response/low rt','prev. slow response/high rt','Location','southwest')
legend boxoff

% Add to p(repeat) matrix
p_repeat = [];
p_repeat = [p_repeat, grandavg{1,1}.repetition, grandavg{1,2}.repetition];

clear grandavg

% Confidence 

for conf = 1:2

for sj = unique(subjects)
    
    % individual sj data
    data = alldata(alldata.subject == subjects(sj), :);
      
    nbins = 6;
      
    for u = 1:nbins
         
        if conf == 1 
            trls = find(data.prevDifficulty == u & data.prevConfidence <= 2);
        elseif conf == 2
            trls = find(data.prevDifficulty == u & data.prevConfidence >= 3);
        end
        
        % =========================================== %
        % history-dependent logistic
        % =========================================== %
        
        % previous response
        resps = [0 1];
        for r = 1:2
                        
            % take a further subset of trials with previous response 1 or 0
            thisdat         = data(trls,:);
            thisdat         = thisdat(thisdat.prevResponseDir == resps(r),:);
            
            [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                thisdat.currRelativeCoherence, thisdat.currResponseDir);
              
            % save betas
            grandavg{conf}.logisticHistory(sj, r, u, :) = [bias slope];
        end
               
    end
          
end

% ========================================================= %
% transform history-dependent bias into repetition bias
% ========================================================= %

resp1 = -1 * squeeze(grandavg{conf}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
resp2 =      squeeze(grandavg{conf}.logisticHistory(:, 2, :, 1));

% since this is centred at 0.5, treat it that way
grandavg{conf}.repetition = (resp1 + resp2) ./ 2;

% ========================================================= %
% transform from log(odds) to probability
% ========================================================= %

logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
grandavg{conf}.repetition = logOdds2Probability(grandavg{conf}.repetition);

end

% save means and SEs for prev low and high conf

mean_repetition_lowconfidence = mean(grandavg{1}.repetition);
sem_repetition_lowconfidence = std(grandavg{1}.repetition)/sqrt(length(subjects));

mean_repetition_highconfidence = mean(grandavg{2}.repetition);
sem_repetition_highconfidence = std(grandavg{2}.repetition)/sqrt(length(subjects));

subplot(2,3,3); hold on;
plot(mean_repetition_highconfidence, '.','Color',[0.4660, 0.6740, 0.1880]);
    v = errorbar(mean_repetition_highconfidence,sem_repetition_highconfidence,'-','Color',[0.4660, 0.6740, 0.1880]);
plot(mean_repetition_lowconfidence,'.','Color',[0.8500, 0.3250, 0.0980]);
    w = errorbar(mean_repetition_lowconfidence,sem_repetition_lowconfidence,'-','Color',[0.8500, 0.3250, 0.0980]);
ylim([.3 .7]); 
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'smallest','','','','','largest'});
ylabel('P(repeat)')
xlabel('Previous |evidence| strength')

legend([v,w],'prev. high confidence (3 or 4)', 'prev. low confidence (1 or 2)','Location','southwest')
legend boxoff;

% Add to p(repeat) matrix
p_repeat = [p_repeat, grandavg{1,1}.repetition, grandavg{1,2}.repetition];

clear grandavg

%% e & f

% Correctness 

for corr = 1:2

for sj = unique(subjects)
    
    % individual sj data
    data = alldata(alldata.subject == subjects(sj), :);
      
    nbins = 6;
      
    for u = 1:nbins
         
        if corr == 1 
            trls = find(data.prevDifficulty == u & data.prevCorrect == 0);
        elseif corr == 2
            trls = find(data.prevDifficulty == u & data.prevCorrect == 1);
        end
        
        % =========================================== %
        % history-dependent logistic
        % =========================================== %
        
        % previous response
        resps = [0 1];
        for r = 1:2
                        
            % take a further subset of trials with previous response 1 or 0
            thisdat         = data(trls,:);
            thisdat         = thisdat(thisdat.prevResponseDir == resps(r),:);
            
            [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                thisdat.currRelativeCoherence, thisdat.currResponseDir);
              
            % save betas
            grandavg{corr}.logisticHistory(sj, r, u, :) = [bias slope];
        end
               
    end
          
end

% ========================================================= %
% transform history-dependent bias into repetition bias
% ========================================================= %

resp1 = -1 * squeeze(grandavg{corr}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
resp2 =      squeeze(grandavg{corr}.logisticHistory(:, 2, :, 1));

% since this is centred at 0.5, treat it that way
grandavg{corr}.repetition = (resp1 + resp2) ./ 2;

% ========================================================= %
% transform from log(odds) to probability
% ========================================================= %

logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
grandavg{corr}.repetition = logOdds2Probability(grandavg{corr}.repetition);

end

% save means and SEs for prev incorrect and correct

mean_repetition_incorrect = mean(grandavg{1}.repetition);
sem_repetition_incorrect = std(grandavg{1}.repetition)/sqrt(length(subjects));

mean_repetition_correct = mean(grandavg{2}.repetition);
sem_repetition_correct = std(grandavg{2}.repetition)/sqrt(length(subjects));

subplot(2,3,5); hold on;
plot(mean_repetition_correct, '.','Color',[0.4660, 0.6740, 0.1880]);
    v = errorbar(mean_repetition_correct,sem_repetition_correct,'-','Color',[0.4660, 0.6740, 0.1880]);
plot(mean_repetition_incorrect,'.','Color',[0.8500, 0.3250, 0.0980]);
    w = errorbar(mean_repetition_incorrect,sem_repetition_incorrect,'-','Color',[0.8500, 0.3250, 0.0980]);
ylim([.3 .7]); 
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'smallest','','','','','largest'});
ylabel('P(repeat)')
xlabel('Previous |evidence| strength')

legend([v,w],'prev. correct', 'prev. incorrect','Location','southwest')
legend boxoff;

% Add to p(repeat) matrix
p_repeat = [p_repeat, grandavg{1,1}.repetition, grandavg{1,2}.repetition];

clear grandavg


% Correctness - P(repeat STIMULUS)

for corr = 1:2

for sj = unique(subjects)
    
    % individual sj data
    data = alldata(alldata.subject == subjects(sj), :);
      
    nbins = 6;
      
    for u = 1:nbins
         
        if corr == 1 
            trls = find(data.prevDifficulty == u & data.prevCorrect == 0);
        elseif corr == 2
            trls = find(data.prevDifficulty == u & data.prevCorrect == 1);
        end
        
        % =========================================== %
        % history-dependent logistic
        % =========================================== %
        
        % previous response
        resps = [-1 1];
        for r = 1:2
                        
            % take a further subset of trials with previous STIMULUS direction -1 or 1
            thisdat         = data(trls,:);
            thisdat         = thisdat(sign(thisdat.prevRelativeCoherence) == resps(r),:);
            
            [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                thisdat.currRelativeCoherence, thisdat.currResponseDir);
              
            % save betas
            grandavg{corr}.logisticHistory(sj, r, u, :) = [bias slope];
        end
               
    end
          
end

% ========================================================= %
% transform history-dependent bias into repetition bias
% ========================================================= %

resp1 = -1 * squeeze(grandavg{corr}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
resp2 =      squeeze(grandavg{corr}.logisticHistory(:, 2, :, 1));

% since this is centred at 0.5, treat it that way
grandavg{corr}.repetition = (resp1 + resp2) ./ 2;

% ========================================================= %
% transform from log(odds) to probability
% ========================================================= %

logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
grandavg{corr}.repetition = logOdds2Probability(grandavg{corr}.repetition);

end

% save means and SEs for prev incorrect and correct

mean_repetition_stim_incorrect = mean(grandavg{1}.repetition);
sem_repetition_stim_incorrect = std(grandavg{1}.repetition)/sqrt(length(subjects));

mean_repetition_stim_correct = mean(grandavg{2}.repetition);
sem_repetition_stim_correct = std(grandavg{2}.repetition)/sqrt(length(subjects));

subplot(2,3,6); hold on;
plot(mean_repetition_stim_correct, '.','Color',[0.4660, 0.6740, 0.1880]);
    v = errorbar(mean_repetition_stim_correct,sem_repetition_stim_correct,'-','Color',[0.4660, 0.6740, 0.1880]);
plot(mean_repetition_stim_incorrect,'.','Color',[0.8500, 0.3250, 0.0980]);
    w = errorbar(mean_repetition_stim_incorrect,sem_repetition_stim_incorrect,'-','Color',[0.8500, 0.3250, 0.0980]);
ylim([.3 .7]); 
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'smallest','','','','','largest'});
ylabel('P(repeat STIMULUS)')
xlabel('Previous |evidence| strength')

legend([v,w],'prev. correct', 'prev. incorrect','Location','southwest')
legend boxoff;

% Change figure size and save
set(gcf, 'Position',  [0, 0, 1300, 800])
saveas(gcf,['P:/3018029.09/Draft/Figure 3 after revision.png'])
saveas(gcf,['P:/3018029.09/Draft/Figure 3 after revision.eps'])

% Save p(repeat) matrix
p_repeat = [p_repeat, grandavg{1,1}.repetition, grandavg{1,2}.repetition];
csvwrite('P:/3018029.09/p_repeat.csv', p_repeat)

clear grandavg


