%% Experimental Design for a Sensitivity Test
% This file is a source code for the algorithm found
% for sequential experimental design for GLM
% and specifically for the simple case of "sensitivity tests"
% using the example from section 4.1 of our paper,
% available at http:/www.math.tau.ac.il/~dms/GLM_Design
%
% This script was written by Hovav Dror, Tel-Aviv University, 2006
%
% If you make changes to the script, or if you implement it
% to a (to be) published problem, please inform us:
% dms@post.tau.ac.il ; hovavdror@gmail.com
%
% further details, including papers describing the algorithm, can be
% found at: http://www.math.tau.ac.il/~dms/GLM_Design


%% Input Parameters

% ScaleMin & ScaleMax designate the physically allowed values
% for the experiment (the design region)
ScaleMin=1;
ScaleMax=50;

ScaleX=[ScaleMin ScaleMax]; % Scaling from experiment grid to [-1,1]
RangeX=ScaleX(2)-ScaleX(1);

% Prior Distribution
NN=10000; % Number of parameter vectors to represent the prior using a discretization
MedianMean=17;
MuS=lognrnd(log(MedianMean),0.5,NN,1);
MedianSigma=0.7;
SigmaS=lognrnd(log(MedianSigma),1,NN,1);

family='binomial';
link='logit';
model='linear';

L=40; % Length of the experiment (total number of observations)

% more initialization
N1=500; % N1 - number of vectors to calculate the weighted sum of determinants of information matrices
RealMin=realmin;
clo=clock;
beta=[-MuS./SigmaS  1./SigmaS];
X0=[]; % This is the matrix with the experimental steps already done
RResult=[];
ss=0;
xrange0=quantile(MuS,[0.1 0.9]);% Basic range for plotting the figures


fprintf('Process began at %g:%g:%1.0f \n',clo(4:6));


%% plot prepartion
%---------------------------------
figure(1)
plot([xrange0(1) xrange0(1) xrange0(2)],[1 0 0],'w'); % empty plot - preparation only
scrsz = get(0,'ScreenSize');
set(gcf,'Position',[2*scrsz(3)/3 2*scrsz(4)/3 scrsz(3)/3 scrsz(4)/3]); % small figure upper right of the screen

figure(2)
plot([xrange0(1) xrange0(1) xrange0(2)],[0 0 1],'w'); % empty plot - preparation only
scrsz = get(0,'ScreenSize');
set(gcf,'Position',[2*scrsz(3)/3 scrsz(4)/6 scrsz(3)/3 scrsz(4)/3]); % small figure upper right of the screen

figure(3)
plot([xrange0(1) xrange0(1) xrange0(2)],[1 0 0],'w'); % empty plot - preparation only
set(gcf,'Position',[0 2*scrsz(4)/3 scrsz(3)/3 scrsz(4)/3]); % small figure upper right of the screen

figure(1);

PMSS=500; % Parameter Space plot Sampe Size - how many vectors of coefficient values to plot
CSS=150; % Curves plot Sample Size - How many curves to plot



Weights=ones(NN,1)/NN;

for Rem=L:-1:1
    fprintf('%g ',ss+1); % the number of the current observation

    % find the Weighted Median (the "best guess" for the parameter vector)
    Mu1=sortrows([MuS Weights]);
    Mu1(:,2)=cumsum(Mu1(:,2)); % Cumulative sum of weightes based on likelihood ratios
    Mu1a=abs( Mu1(:,2) - 0.5 );
    MedMu=Mu1(Mu1a==min(Mu1a),1);
    MedMu=median(MedMu); % in case there was more than one Mu that met the criterion

    Sigma1=sortrows([SigmaS Weights]);
    Sigma1(:,2)=cumsum(Sigma1(:,2)); % Cumulative sum of weightes based on likelihood ratios
    Sigma1a=abs( Sigma1(:,2) - 0.5 );
    MedSigma=Sigma1(Sigma1a==min(Sigma1a),1);
    MedSigma=median(MedSigma); % in case there was more than one Mu that met the criterion

    betaMed=[-MedMu/MedSigma 1/MedSigma];

    %% Now, use this "best guess" to choose the next point
    betaMedX(1,1)=betaMed(1)+betaMed(2)*(ScaleX(1)+0.5*(RangeX)); % Scale
    betaMedX(1,2)=betaMed(2)*ScaleX(2)-betaMed(2)*(ScaleX(1)+0.5*RangeX);
    if betaMedX(2)<25 % if not very sharp curve
        Y0=-1+(X0-ScaleX(1))./RangeX.*2;
        X0t=GLMDAUG(Y0,betaMedX',1,ss+2); % Augment 2 points, D-optimaly
        X0t=ScaleX(1)+(X0t+1)./2.*RangeX; % Scale Back
    else % if very sharp curve - rescale
        ScaleX2=[max([MedMu-4*MedSigma ScaleX(1)]) min([MedMu+4*MedSigma ScaleX(2)])];
        RangeX2=ScaleX2(2)-ScaleX2(1);
        betaMedX2(1,1)=betaMed(1)+betaMed(2)*(ScaleX2(1)+0.5*(RangeX2));
        betaMedX2(1,2)=betaMed(2)*ScaleX2(2)-betaMed(2)*(ScaleX2(1)+0.5*RangeX2);
        Y0=-1+(X0-ScaleX2(1))./RangeX2.*2;
        X0t=GLMDAUG(Y0,betaMedX2',1,ss+2);% Augment 2 points, D-optimaly
        X0t=ScaleX2(1)+(X0t+1)./2.*RangeX2;
    end;

    % X0t has the candidates points for the choice of the next point
    X0t=[median(X0t) ;X0t ]; % Add the median of the 2 points augmentation
    BestSoFar=-inf;
    for i=1:3 % there are 3 candidates
        ThisOne=0;
        X0ti=[X0; X0t(i)];
        if ss==0, X0ti=[X0ti ; X0t(2:3,:)]; end;
        F0ti=x2fx(X0ti);
        for j=1:N1   % Calculate each candidates contribution to the posterior
            ThisOnej=max(det(InfoMtrxGLM(F0ti,beta(j,:)',family,link) ) ,RealMin);
            ThisOne=ThisOne + Weights(j)*log( ThisOnej );
        end;
        if ThisOne>BestSoFar
            BestSoFar=ThisOne;
            X0Best=X0t(i);
        end;
    end;

    RResultTemp=[];
    %% **********************plot***************************

    figure(1) % upper right figure, prior vs. posterior distribution
    plot(MuS(1:PMSS),SigmaS(1:PMSS),'k.','MarkerSize',0.5);% Original Sample is small dots
    hold on
    SampleVectors=randsample(NN,PMSS,true,Weights); % Choose a weighted random sample
                                                                                                                                                        % representing the Posterior
    plot(MuS(SampleVectors),SigmaS(SampleVectors),'k.'); % Plot Weighted sample (large dots)
    plot(MedMu,MedSigma,'o','MarkerFaceColor','g'); % Estimated point ("best guess") in green
    hold off
    xlabel('Mu'); ylabel('Sigma');

    figure(3) % Same as figure(1) but on a log scale (prior vs. posterior), upper left figure
    loglog(MuS(1:PMSS),SigmaS(1:PMSS),'k.','MarkerSize',0.5);% Original Sample
    hold on
    plot(MuS(SampleVectors),SigmaS(SampleVectors),'k.'); % Weighted sample
    plot(MedMu,MedSigma,'o','MarkerFaceColor','g'); % Estimated point
    text(min(MuS(1:PMSS)),max(SigmaS(1:PMSS)),sprintf('     est: \\mu = %2.2f, \\sigma = %2.2f',[MedMu MedSigma]))
    hold off
    xlabel('Mu'); ylabel('Sigma');

    figure(2) % lower right, a sample of possible models, from the posterior
    plot([xrange0(1) xrange0(1) xrange0(2)],[0 0 1],'w'); % empty figure
    hold on
    P=[];
    for ii=1:CSS % Plot a sampe of possible models, from the Weighted Sample
        j=0;
        xrange=min([min(X0) xrange0(1)]):.01:max([max(X0) xrange0(2)]); % range for plotting figures
        for x=xrange
            j=j+1;
            P(j)=1 / (1 + exp(-beta(SampleVectors(ii),1)-beta(SampleVectors(ii),2)*x) );
        end;
        plot(xrange,P,'k');
    end; %for
    j=0; % and plot the current guess setup, in green
    P=[];
    for x=xrange
        j=j+1;
        P(j)=1/(1+exp(MedMu/MedSigma-(1/MedSigma)*x));
    end;
    plot(xrange,P,'g','LineWidth',4);

    for jj=1:(L-Rem) % plot the location of observations, and their outcome
        if RResult(jj)==1
            plot(X0(jj),jj/L,'^','MarkerSize',8,'MarkerFaceColor','m'); % upper triangle - "response"
        else
            plot(X0(jj),jj/L,'V','MarkerSize',8,'MarkerFaceColor','c'); % lower triange - "no response"
        end;
    end;
    % and emphasize the next chosen point, in yelow
    plot(X0Best,0.5,'s','MarkerSize',8,'MarkerFaceColor','y');
    set(gca,'XLim',[min(xrange) max(xrange)]);
    xlabel('Dose Value')
    ylabel('Probability of positive Response');

    % if we can get MLE by regression, plot it, too:
    if ss>2 
        if max(X0(RResult==0))>min(X0(RResult==1))
                [b,dev,stats] = glmfit(X0,[RResult' ones(ss,1)],'binomial',link);
            hold off
            figure(3);
            hold on;
            text(min(MuS(1:PMSS)),max(SigmaS(1:PMSS)),sprintf('\n     reg: \\mu = %2.2f, \\sigma = %2.2f',[-b(1)/b(2) 1/b(2)]));
            hold off;

            figure(2);
            hold on;
            [yfit,dlo,dhi] = glmval(b,xrange',link,stats,0.95);
            plot(xrange,yfit-dlo,'r:','LineWidth',3);
            plot(xrange,yfit+dhi,'r:','LineWidth',3);
        end;
    end;
    hold off

    drawnow;
    %**********************endplot**********************

    % Recommend a value for the next observations
    text1=sprintf('Recommendation: %g  ; "enter"=accept, or enter value: ',X0Best);
    XChosen=input(text1); % input where to take the next observations
    if isempty(XChosen), XChosen=X0Best; end;

    X0=[X0 ; XChosen];
    ss=ss+1;
    F0=x2fx(X0);

    while isempty(RResultTemp) || ~isnumeric(RResultTemp) || (RResultTemp~=0 & RResultTemp~=1)
        RResultTemp=input('Enter the result (0 /1): ');  % Input the result of the observation (0=no response, 1=response)
    end;
    RResult(end+1)=RResultTemp;

    %% Now modify the Posterior distribution:
    RResultM=ones(NN,1)*RResult;
    ex=exp(-F0*beta');
    bprob=(1 ./ (1+ex))';
    mprob=1-bprob;
    LRi = prod( bprob.^RResultM ,2 ) .* prod( mprob.^(1-RResultM) ,2);
    Weights=LRi/sum(LRi);

end; % for Rem






