% mei_lev - implements mei algorithm and selects next stimulus level
% [adapted from % Andrea Trevino - 4/23/15]

function [stimlevel,xlikeli,pcat]=mei_lev(loglikeli,storelikeli) 
% response has been recorded in CLS.m under listener_data=[lev loudness_resp]
% ID is the listener name
% ear is R or L
% listener data is [freq stimlevel catresp phase]
% iter is stimuli_number

nce=size(storelikeli,1); % number of catelog entries
nlv=size(storelikeli,2); % number of stimlus levels
ncu=size(storelikeli,3); % number of loudness categories
likeli=exp(storelikeli); % nce x nlv x ncu, general likelihoods
xlikeli=exp(loglikeli); % nce x 1, posterior likelihoods
% posterior probability
post=xlikeli/sum(xlikeli); % 1460 x 1, the posterior for each theta_i
% current entropy H_t
%ht=-1*sum(post.*log2(post));  % won't actually be used in alg, but good reference
%meient(mc,(rep-1)*numreps+s)=ht;
% FOR EACH LEVEL, find x+1=L likelihood, posterior, and entropy
hnext=zeros(ncu,nlv); % cat x levels, h_{t+1}
pcat=zeros(ncu,nlv); % cat x levels, P(L|X_t)
for cat=1:ncu
    lognextlikeli=repmat(log(xlikeli),1,nlv)+squeeze(storelikeli(:,:,cat));
    nextlikeli=exp(lognextlikeli);
    sumnextlikeli=repmat(sum(nextlikeli),nce,1);
    nextpost=nextlikeli./sumnextlikeli;
    log2nextpost=(lognextlikeli-log(sumnextlikeli))/log(2);
    hnext(cat,:)=-1*sum(nextpost.*log2nextpost);
    pcat(cat,:)=sum(squeeze(likeli(:,:,cat)).*repmat(post,1,121));
end
exph=sum(pcat.*hnext); % expected entropy, sum over loudness categories
[~,ilv]=min(exph); % find lev that minimizes exp entropy
stimlevel=ilv-1; % fix due to indices
end