function tae_ss_pf = tae_ss_fit_pf(cfg)

addpath ../Toolboxes/Palamedes

subject     = cfg.subject;
data        = cfg.data;

stim_levels = -12:6:12;
    

for iStimLevel = 1:length(stim_levels)

    numPos(iStimLevel) = nansum(data(data(:,5) == stim_levels(iStimLevel),9));
    OutOfNum(iStimLevel) = sum(data(:,5) == stim_levels(iStimLevel)) - sum(isnan(data(data(:,5) == stim_levels(iStimLevel),9)));

    p(iStimLevel) = numPos(iStimLevel)/OutOfNum(iStimLevel);

    numPos_bias_ccw(iStimLevel) = nansum(data(data(:,5) == stim_levels(iStimLevel) & data(:,6) == -1,9));
    OutOfNum_bias_ccw(iStimLevel) = sum(data(:,5) == stim_levels(iStimLevel) & data(:,6) == -1) - sum(isnan(data(data(:,5) == stim_levels(iStimLevel) & data(:,6) == -1,9)));

    p_bias_ccw(iStimLevel) = numPos_bias_ccw(iStimLevel)/OutOfNum_bias_ccw(iStimLevel);

    numPos_bias_cw(iStimLevel) = nansum(data(data(:,5) == stim_levels(iStimLevel) & data(:,6) == 1,9));
    OutOfNum_bias_cw(iStimLevel) = sum(data(:,5) == stim_levels(iStimLevel) & data(:,6) == 1) - sum(isnan(data(data(:,5) == stim_levels(iStimLevel) & data(:,6) == 1,9)));

    p_bias_cw(iStimLevel) = numPos_bias_cw(iStimLevel)/OutOfNum_bias_cw(iStimLevel);

end

% Fit psychometric function
PF = @PAL_CumulativeNormal;

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

%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 = -10:.001:10;
searchGrid.beta = logspace(0,3,101);
searchGrid.gamma = 0.01;  %scalar here (since fixed) but may be vector
searchGrid.lambda = 0.01;  %ditto
%searchGrid.gamma = 0:.01:.15;  %scalar here (since fixed) but may be vector
%searchGrid.lambda = 0:.01:0.15;  %ditto


%Perform fit
[pf_ccw, ~, ~] = ...
    PAL_PFML_Fit(stim_levels,numPos_bias_ccw, OutOfNum_bias_ccw,searchGrid,paramsFree,PF);
    %PAL_PFML_Fit(stim_levels,nTrials_right_low_noise, outOfNum_low_noise,searchGrid,paramsFree,PF);

[pf_cw, ~, ~] = ...
    PAL_PFML_Fit(stim_levels,numPos_bias_cw, OutOfNum_bias_cw,searchGrid,paramsFree,PF);
    %PAL_PFML_Fit(stim_levels,nTrials_right_high_noise, outOfNum_high_noise,searchGrid,paramsFree,PF);


%% Write in output structure

tae_ss_pf.pf_params_cww         = pf_ccw;
tae_ss_pf.pf_params_cw          = pf_cw;
tae_ss_pf.p_bias_ccw        	= p_bias_ccw;
tae_ss_pf.p_bias_cw             = p_bias_cw;


%% Show figure if opted
if cfg.plot == true

    StimLevelsFineGrain=[min(stim_levels):max(stim_levels)./1000:max(stim_levels)];

    fit_ccw =  fit(stim_levels', p_bias_ccw', 'smoothingspline');
    fit_cw  =  fit(stim_levels', p_bias_cw', 'smoothingspline');

    figure; hold on;
    plot(stim_levels, p_bias_ccw,'+r');
    plot(stim_levels, p_bias_cw,'+b');
   
    g1 = plot(fit_ccw, '-r');
    g2 = plot(fit_cw, '-b');
    
    xlabel('Stimulus orientation')
    ylabel('p("right" response)')
    
    xlim([-14 14]);
    ylim([0 1]);
    
    set(gca,'XTick',stim_levels)
    set(gca,'XTickLabel',{'-12';'-6';'0';'6';'12'});
    
    legend([g1 g2], 'ccw inducer','cw inducer')

    % title(['Subject ' num2str(subject)]);
    title(['Subject ' num2str(subject)]);

    hold off;
    
    % saveas(gcf,['Figures/SingleSubjectFits/tae_ss_fit_pf_' num2str(subject)'.png'])

end

end

