% Andrea Trevino
% 1/15/15

% Function to use a CLSPF catalog to estimate the CLS curve of inputted
% listener data

% data format is 2 columns of trial-by-trial data, with the first column
% being the level and the second being the CU response
% range of input levels is (0, 120) dB

% makeplot is 1 or 0 depending on whether you want to see a plot or not

function [clslevel,clsCU,clspf]=catalogCLS_1db(storelikeli,data,makeplot,storeclspf)

data(:,1)=1+data(:,1); % make into indices for matrix
data(:,2)=1+data(:,2)/5; % make into indices for matrix

tot_like=zeros(size(storelikeli,1),1); % total likelihood to update iteratively
for t=1:length(data)
    tot_like=tot_like+squeeze(storelikeli(:,data(t,1),data(t,2))); % likelihood is saved in log form, so sum is total
end
[temp,loc]=max(tot_like); % find the maximum likelihood family
ml=exp(squeeze(storelikeli(loc,:,:))); % select the family
clspf=cumsum(ml')';
clspf=clspf(:,1:end-1); % remove a category that isn't there
% done w/clspfs

% find the cls function
this_cls=storeclspf(loc,:);
fa=0.01;
D=[1-10*fa:fa:1-fa];
A=[fa:fa:fa*10];
B=1./this_cls(1:10);
C=[this_cls(11) this_cls(11)+cumsum(this_cls(12:end))]+20;
modelclspf=[];
for cat=1:10
    modelclspf=[modelclspf; D(cat)+((A(cat)-D(cat))./(1+(([0:120]+20)./C(cat)).^(B(cat))))];
end
[temp,loc]=min(abs(modelclspf'-0.5));
clslevel=loc-1;
clsCU=[5:5:50]-2.5;

if makeplot
    figure(1)
    clf(1)
    subplot(1,2,1)
    set(gca,'FontSize',14)
    % plot resulting family
    plot([0:120],clspf,':')            
    % show the raw data overlaid
    cat=unique(data(:,2));
    levels=unique(data(:,1));
    origclspf=zeros(length(levels),length(cat));
    for c=1:length(cat)
        for lev=1:length(levels)
            origclspf(lev,c)=sum(data(data(:,2)==cat(c),1)==levels(lev))/sum(data(:,1)==levels(lev));
        end
    end
    hold on
    plot(levels-1,cumsum(origclspf')','-')
    title('CLSPFs')
    ylim([0 1]) 
    xlim([0 125])
    xlabel('Level (dB SPL)')
    ylabel('Cumulative Probability')
    set(gca,'XTick',[0:20:100])
    
    subplot(1,2,2)
    set(gca,'FontSize',14)
    % plot CLS curve
    plot(clslevel,clsCU,'k.-','LineWidth',1.5)
    ylim([0 50]) 
    xlim([0 120])
    xlabel('Level (dB SPL)')
    ylabel('CU')
    title('CLS Curve')
    set(gca,'XTick',[0:20:100])
    xlim([-20 120])
end