% plot rms_err and entropy as functions versus trials for MEI tracking
function adtrksim
global all_tbt all_ent sse cnt
load catalog2
nce=size(storeclspf,1);
fa=0.005;
te=6;
nt=45;
nr=2;
lev=(0:120)';
[lis,med]=simlis_prepare('coeffdata',lev,fa);
ns=size(lis,1);
nf=size(lis,2);
all_tbt=zeros(nr,ns,nf,nt+2,2);
all_ent=zeros(nr,ns,nf,nt+1);
sse=zeros(nt+1,1);
cnt=zeros(nt+1,1);
ent=2*log2(nce);
fprintf('catalog: size=%d, entropy=%.1f bits\n',nce,ent)
get_err(lev,lis,med,nr,ns,nf,nt,storelikeli,storeclspf);
re=zeros(nt+1,1);
for t=0:nt
   rmserr=sqrt(sse(t+1)/cnt(t+1));
   avgent=sum(mean(mean(all_ent(:,:,:,t+1),1),2),3);
   fprintf('%3d: rms_err = %6.2f, entropy = %6.2f\n',t,rmserr,avgent)
   re(t+1)=rmserr;
end
tt=find(re<te,1)-1;
if (isempty(tt))
    ts=sprintf('> %d',nt);
else
    ts=sprintf('= %d',tt);
end
fprintf('target: err = %.2f, trials %s\n',te,ts)
%
figure(1);clf
xlim=[1 100];
x=1+(0:nt)';
ent=squeeze(sum(mean(mean(all_ent,1),2),3));
semilogx(x,re,x,ent)
axis([xlim 0 20])
legend('dB rms','entropy')
xlabel('trial number')
write_data('adtrk1.txt',[x re ent])
%
plot_subj( 1,2,2);
plot_subj(25,2,3);
%
return

function get_err(lev,lis,med,nr,ns,nf,nt,storelikeli,storeclspf)
tbt=zeros(nt+2,2);
ent=zeros(nt+1,1);
for j=1:ns                      % loop over subject
   for i=1:nf                   % loop over frequency
      ls=squeeze(lis(j,i,:,:)); % listener probabilty
      md=squeeze(med(j,i,:))';  % median SPL
      for r=1:nr                % loop over repeats
         % initialize range
         lv_range=get_range(lev,ls);
         cu_range=[5,50];
         loglikeli=zeros(size(storelikeli,1),1);
         for k=1:2
            trial=[lv_range(k),cu_range(k)];
            loglikeli=postlike(trial,loglikeli,storelikeli);
            tbt(k,:)=trial;
         end
         ent(1)=get_ntrpy(loglikeli);
         dyn=abs(lv_range(1)-lv_range(2)); % dynamic range
         prv=tbt(1,1);                 % set previous level to threshold
         bot=max(  0,lv_range(1)-10);  % minimum level
         top=min(110,lv_range(2)+10);  % maximum level
         mnch=round(dyn/5);            % minimum level change
         mxch=round(dyn/3);            % maximum level change
         % perform trials
         for t=1:nt
            [stm,~,pcat]=mei_lev(loglikeli,storelikeli);
            stm=limit_level(stm,prv,bot,top,mnch,mxch,pcat);
            rsp=get_resp(ls(lev==stm,:));         % simulate listener
            trial=[stm,rsp];
            loglikeli=postlike(trial,loglikeli,storelikeli);
            % store tbt and ent
            tbt(t+2,:)=trial;
            ent(t+1)=get_ntrpy(loglikeli);
            prv=stm;
         end
         % store rms error and entropy
         store_sse_cnt(nt,tbt,md,storelikeli,storeclspf);
         store_tbt_ent(r,j,i,tbt,ent);
      end
   end
end
return

function stm=limit_level(stm,prv,bot,top,mnch,mxch,pcat)
% convert level arguments to indices
istm = stm + 1;
iprv = prv + 1;
itop = min(top,prv+mxch) + 1;
ibot = max(bot,prv-mxch) + 1;
jtop = min(itop,iprv+mnch);
jbot = max(ibot,iprv-mnch);
% apply "max change" restriction
if (istm>itop) istm=itop; end
if (istm<ibot) istm=ibot; end
% apply "min change" restriction
if ((istm<jtop)&&(istm>jbot))
   imid=(itop+ibot)/2;
   if (abs(jtop-imid)<abs(jbot-imid))
      istm = jtop;
   else
      istm = jbot;
   end
end
% calculate expected loudness
[~,jstm]=max(pcat(:,istm));
% apply "expected inaudible" restriction
while ((jstm<1)&&(istm<itop))
   istm = istm + 1;
   [~,jstm]=max(pcat(:,istm));
end
% apply "expected too-loud" restriction
while ((jstm>10)&&(istm>ibot))
   istm = istm - 1;
   [~,jstm]=max(pcat(:,istm));
end
% convert stimulus index to stimulus level
stm = istm - 1;
return

function [range,trial]=get_range(lev,lis)
thr=min(lev);
% 8 fixed stim, so they get a feel for the range
stim=[50 35 20 5 15 0 25 45];
trial=0;
for k=1:length(stim)
    stm=stim(k);
    if (stm<(thr-30)) 
        continue; 
    end
    rsp=get_resp(lis(lev==stm,:));
    trial=trial+1;
    if ((rsp == 00) && (thr<stm)) 
        thr=stm; 
    elseif (thr>stm)
        thr=stm - 5; 
    end
end
% 4 stim re threshold
sl=[-10 0 10 20];
for k=1:length(sl)
    stm=thr+sl(k);
    if (stm<0)
        continue; 
    end
    rsp=get_resp(lis(lev==stm,:));
    trial=trial+1;
    if ((rsp == 00) && (thr<stm)) 
        thr=stm; 
    elseif (thr>stm)
        thr=stm - 5; 
    end
end
%range(1)=max(0,thr-10);
range(1)=max(0,thr); % lower limit is audible
tool=0;
L105=0;
ucl=110;
mcl=max(65,thr);
stm=mcl;
while ((tool<2)&&(L105<3))
    rsp=get_resp(lis(lev==stm,:));
    trial=trial+1;
    if (rsp==50)
        tool=tool+1;
        if (ucl>stm)
            ucl=stm;
        end
    elseif (ucl<stm)
        ucl=stm+5;
    end
    if (stm==105)
        L105=L105+1;
        stm=mcl;
    elseif (stm<90)
        stm=stm+10;
    else
        stm=stm+5;
    end
end
range(2)=ucl;
return

function loglikeli=postlike(trial,loglikeli,storelikeli)
ilv=1+trial(1);   % index level
icu=1+trial(2)/5; % index CU
% likelihood is saved in log form, so sum is total
loglikeli=loglikeli+squeeze(storelikeli(:,ilv,icu)); 
return

%========================================

function [lis,med]=simlis_prepare(dfn,lev,fa)
load(dfn);
nlv=length(lev);
nnh=size(all_coeff1,1);
nhi=size(all_coeff1_HL,1);
ns=nnh+nhi;
coeff=zeros(2,ns,20);
coeff(1,:,:)=[all_coeff1;all_coeff1_HL];
coeff(2,:,:)=[all_coeff4;all_coeff4_HL];
lis=zeros(ns,2,nlv,10);
med=zeros(ns,2,10);
for i=1:2
    for j=1:ns
        C=0;
        for k=1:10
            A=k*fa;
            B=coeff(i,j,k);
            C=coeff(i,j,k+10)+C;
            D=1+(k-10)*fa;
            y=D+(A-D)./(1+(lev./C).^B);
            lis(j,i,:,k)=y;
            med(j,i,k)=C;
        end
    end
end
return

function resp=get_resp(lisstm)
r=rand(1);
if (r>lisstm(10))
    resp=50;
else
    resp=(find(r<lisstm,1)-1)*5;
end
return

function store_sse_cnt(nt,tbt,md,storelikeli,storeclspf)
global sse cnt
for t=0:nt
   tbtt=tbt(1:(t+2),:);
   lv=catalogCLS_1db(storelikeli,tbtt,0,storeclspf);
   sse(t+1)=sse(t+1)+sum((lv-md).^2);
   cnt(t+1)=cnt(t+1)+length(lv);
end
return

function store_tbt_ent(r,j,i,tbt,ent)
global all_tbt all_ent
all_tbt(r,j,i,:,:)=tbt;
all_ent(r,j,i,:)=ent;
return

%========================================

function plot_subj(sbj,frq,fig)
global all_tbt
figure(fig);clf
tbt=squeeze(all_tbt(1,sbj,frq,:,:));
lv=tbt(:,1);
cu=tbt(:,2);
nt=size(tbt,1)-2;
t=(-1:nt)';
plot(t,lv,'b',t,cu,'r--');
axis([-1 max(t)+5 -5 115])
f=[1 4];
title(sprintf('sbj = %d, f = %d kHz',sbj,f(frq)))
legend('level','CU')
xlabel('trial number')
write_data(sprintf('adtrk%d.txt',fig),[t lv cu])
return

%========================================

function write_data(fn,data)
[nr,nc] = size(data);
fp=fopen(fn,'wt');
fprintf(fp,'; %s\n', fn);
for i=1:nr
    for j=1:nc;
        fprintf(fp,' %14.5g',data(i,j));
    end
    fprintf(fp,'\n');
end
fclose(fp);
return
