% plot rms_err as a function of number of trials for MEI tracking
function adtrk2
global all_tbt all_ent
load catalog2
nce=size(storeclspf,1);
fa=0.005;
te=6;
nt=15;
nr=2;
lev=(0:120)';
[lis,med]=simlis_prepare('coeffdata',lev,fa);
ns=size(lis,1);
nf=size(lis,2);
re=zeros(nt,1);
all_tbt=zeros(nr,ns,nf,nt+2,2);
all_ent=zeros(nr,ns,nf,nt);
ent=2*log2(nce);
fprintf('catalog: size=%d, entropy=%.1f bits\n',nce,ent)
get_err(lev,lis,med,nr,0,storelikeli,storeclspf);
for t=1:nt
    re(t)=get_err(lev,lis,med,nr,t,storelikeli,storeclspf);
end
tt=find(re<te,1)*nf;
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:nt)*nf)';
ent=squeeze(sum(mean(mean(all_ent,1),2),3));
semilogx(x,re,x,ent)
axis([xlim 0 20])
legend('dB rms','entropy')
write_data('adtrk2.txt',[x re ent])
return

function rmserr=get_err(lev,lis,med,nr,t,storelikeli,storeclspf)
global all_tbt all_ent
ns=size(lis,1);
nf=size(lis,2);
mkplt=0;
sserr=0;
count=0;
for r=1:nr
   for j=1:ns
      for i=1:nf
         tbt=squeeze(all_tbt(r,j,i,:,:));
         lisji=squeeze(lis(j,i,:,:));
         [tbt,ent]=track(t,lev,lisji,storelikeli,tbt);
         tbtt=tbt(1:(t+2),:);
         lv=catalogCLS_1db(storelikeli,tbtt,mkplt,storeclspf);
         md=squeeze(med(j,i,:))';
         sserr=sserr+sum((lv-md).^2);
         count=count+length(lv);
         all_tbt(r,j,i,:,:)=tbt;
         if (t>0) all_ent(r,j,i,t)=ent; end
      end
   end
end
rmserr=sqrt(sserr/count);
fprintf('%3d: rms_err = %6.2f\n',t,rmserr)
return

function [tbt,ent]=track(t,lev,lis,storelikeli,tbt)
if (t<1)
   range=get_range(lev,lis);
   tbt(1:2,:)=[range(1),0;range(2),50];
   loglikeli=zeros(size(storelikeli,1),1);
   for k=1:2
      loglikeli=postlike(tbt(k,:),loglikeli,storelikeli);
   end
   ent=get_ntrpy(loglikeli);
else
   % initialize likelihood based on previous trials
   loglikeli=zeros(size(storelikeli,1),1);
   for k=1:(t+1)
      loglikeli=postlike(tbt(k,:),loglikeli,storelikeli);
   end
   % adaptive trials
   stm=mei_lev(loglikeli,storelikeli);
   rsp=get_resp(lis(lev==stm,:));
   tbt(t+2,:)=[stm,rsp];
   loglikeli=postlike(tbt(t+2,:),loglikeli,storelikeli);
   ent=get_ntrpy(loglikeli);
end
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 [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);
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 resp=get_resp(lisstm)
r=rand(1);
if (r>lisstm(10))
    resp=50;
else
    resp=(find(r<lisstm,1)-1)*5;
end
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
