%% Sequential Experimental Design:
%% Two competing models + Poisson response
%
% This file is a Monte-Carlo examination, extending the file
% TwoModelsAndPoissonSequential.m
% In each iteration, a "true" setup is chosen at random
% from the prior distribution,
% The algorithm runs automatically,
% and the Local D-Efficiency is recorded.
%
% This file is a source code for the algorithm found
% for multivariate sequential experimental design for GLM
% using the example from section 4.5 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

MaxIter=1000; % Number of quasi-monter carlo iterations
iter=0;

PLD=8; % Recommended number of points in a good Local D-optimal design
L=48; % Length of the experiment (total number of observations)
m=5; %number of variables
family='poisson';
link='log';
model1=[zeros(1,m) ; eye(m)  ;  1 1 zeros(1,m-2) ;  1 0 1 zeros(1,m-3)];
model2='linear';
        if isnumeric(model1)
            p1=length(model1);
        else
            p1=length(x2fx(ones(1,m),model1));
        end;
        if isnumeric(model2)
            p2=length(model2);
        else
            p2=length(x2fx(ones(1,m),model2));
        end;
% Prior distribution        
BetaMeans1=[   -2.3479   -5.5296   -2.9929   -3.9494   -0.8555    0.4139   -2.0666   -1.1263];
BetaSE1=[ 0.6863    0.9408    0.8164    0.5876    0.5389    0.3649    1.3195    0.9848];
BetaMeans2=[   -1.5165   -4.2972   -1.7859   -3.3852   -0.2806    0.2282];
BetaSE2=[  0.2058    0.2011    0.1581    0.2367    0.3185    0.3026];
NN=5000;
NI2=rand(NN,max([p1 p2]));
NI1=NI2(1:NN,1:p1);
NI2=NI2(1:NN,1:p2);
% Now transform it according to distribution
beta1 = norminv(NI1,ones(NN,1)*BetaMeans1,ones(NN,1)*BetaSE1);
beta2 = norminv(NI2,ones(NN,1)*BetaMeans2,ones(NN,1)*BetaSE2);

N=500; % N = number of vectors in the low-descrepancy sequence representing the parameter space
N1=500; % N1 - number of vectors to calculate the weighted sum of determinants of information matrices

RealMin=realmin;
RealMax=realmax;
LRealMin=log(RealMin);
LRealMax=log(RealMax);

Eff=[];
clo=clock;

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

%% begin iterations
while iter<MaxIter
    iter=iter+1;
    %% "True" parameters
    Tmodel=1+round(rand(1)); % toss a coin which of the two models is true
    if Tmodel==1  % and choose coefficient values
        NI=rand(1,p1);
        Tbeta=norminv(NI,BetaMeans1,BetaSE1);
        CBeta=Tbeta;
    else
        NI=rand(1,p2);
        Tbeta=norminv(NI,BetaMeans2,BetaSE2);
        CBeta=Tbeta;
    end;
    X0=[]; % This is the matrix with the experimental steps already done
    RResult=[];
    ResultFactorial=[];
    SumModelWeights=[];
    ss=0; % ss=size(X0,1);

    Weights1=0.5*ones(NN,1)/NN;
    Weights2=0.5*ones(NN,1)/NN;
    
    for Rem=L:-1:1

        fprintf('%g ',Rem);
        for modelN=1:2
            if modelN==1
                p=p1;
                beta=beta1;
                Weights=Weights1;
            else
                p=p2;
                beta=beta2;
                Weights=Weights2;
            end;
            betaMed=[];
            AA=(1:NN)';
            for i=1:p
                betai=sortrows([beta(:,i) Weights]);
                betai=betai(betai(:,2)>0,:);
                if ~isempty(betai)
                    betai(:,3)=cumsum(betai(:,2));
                    betai(:,3)=abs(betai(:,3)-betai(end,3)/2);
                    betaMedi=betai(betai(:,3)==min(betai(:,3)),1);
                    betaMed(1,i)=median(betaMedi);% in case there was more than one Mu that met the criterion
                else
                    betaMed(1,i)=beta(1,i); % just a vector, zero weight anyhow
                end;
            end;
            if modelN==1, betaMed1=betaMed; else betaMed2=betaMed; end;
        end;
         
        % plots of the estimated coefficient values, vs. the "true" model
        % (extreme right point)
        if Tmodel==1  % just to help see the convergence of the algorithm
            CBeta(end,:)=betaMed1;
            CBeta(end+1,:)=Tbeta;
            SumModelWeights=[SumModelWeights ; sum(Weights1)];
        else
            CBeta(end,:)=betaMed2;
            CBeta(end+1,:)=Tbeta;
            SumModelWeights=[SumModelWeights ; sum(Weights2)];
        end;
        figure(1);
        plot(CBeta,'.-');
        hold on
        plot(SumModelWeights,'k-');
        text(ss+1,SumModelWeights(end),sprintf('%1.2f',SumModelWeights(end)));
        set(gca,'XLim',[0 50]);
        hold off
        xlabel('Observation number');
        ylabel('coefficient values');
        drawnow
            
        X0t1=GLMDAUG(X0,betaMed1',m,ss+PLD,model1,family,link); 
        X0t2=GLMDAUG(X0,betaMed2',m,ss+PLD,model2,family,link); 
        X0t1=[median(X0t1) ;X0t1 ];
        X0t2=[median(X0t2) ;X0t2 ];
        X0t=[X0t1 ; X0t2];
        BestSoFar=-inf;
        if max(Weights1)>0
            SampleVectors1=randsample(NN,N1,true,Weights1);
        else
            SampleVectors1=1:N1;
        end;
        if max(Weights2)>0
            SampleVectors2=randsample(NN,N1,true,Weights2);
        else
            SampleVectors2=1:N1;
        end;

        for i=1:size(X0t,1)
            ThisOne=0;
            X0ti=[X0; X0t(i,:)];
            if ss<(PLD-1), X0ti=[X0ti ; X0t1(2:end,:)]; end;
            F0ti1=x2fx(X0ti,model1);
            F0ti2=x2fx(X0ti,model2);
            
            for j=1:N1
                 ThisOnej=max(det(InfoMtrxGLM(F0ti1,beta1(SampleVectors1(j),1:p1)',family,link) ) ,RealMin);
                 ThisOne=ThisOne + Weights1(SampleVectors1(j))*log( ThisOnej );
                 ThisOnej=max(det(InfoMtrxGLM(F0ti2,beta2(SampleVectors2(j),1:p2)',family,link) ) ,RealMin);
                 ThisOne=ThisOne + Weights2(SampleVectors2(j))*log( ThisOnej );
            end;
            if ThisOne>BestSoFar
                BestSoFar=ThisOne;
                X0Best=X0t(i,:);
            end;
        end;

        X0=[X0 ; X0Best];
        ss=ss+1;
        F01=x2fx(X0,model1);
        F02=x2fx(X0,model2);
        
        %% randomize result
        if Tmodel==1
            lambda=exp(Tbeta*F01(end,:)');
        else
            lambda=exp(Tbeta*F02(end,:)');
        end;
        RResult(end+1)=poissrnd(lambda);
        if RResult(end)<21
            ResultFactorial(end+1)=factorial(RResult(end));
        else
            ResultFactorial(end+1)=-1; % else - use normal approximation and not the exact poisson distribution
        end;
        
        
%% modify the Posterior distribution

        RealMinSS=RealMin^(1/ss);
        RResultM=ones(NN,1)*RResult;
        bx1=(F01*beta1')';
        lambda1=exp(bx1);
        LRi1=[]; LRi2=[];
        for i=1:NN
            LRi1(i)=0;
            for j=1:ss
                if lambda1(i,j)>20 | RResult(j)>20  % if it is better to use normal approximation
                    LRi1(i)=LRi1(i)  -0.5*log(2*pi*lambda1(i,j)) - (RResult(j) - lambda1(i,j))^2 / (2*lambda1(i,j)) ;
                else % for small lambda & RResult use exact poisson distribution
                    LRi1(i)=LRi1(i) +log( max([lambda1(i,j)^RResult(j) * exp(-lambda1(i,j)) / ResultFactorial(j) RealMinSS]));
                end; % if lambda
            end; % for j
        end; % for i
        bx2=(F02*beta2')';
        lambda2=exp(bx2);
        for i=1:NN
            LRi2(i)=0;
            for j=1:ss
                if lambda2(i,j)>20 | RResult(j)>20  % if it is better to use normal approximation
                    LRi2(i)=LRi2(i)  -0.5*log(2*pi*lambda2(i,j)) - (RResult(j) - lambda2(i,j))^2 / (2*lambda2(i,j)) ;
                else % for small lambda & RResult use exact poisson distribution
                    LRi2(i)=LRi2(i) +log (max([lambda2(i,j)^RResult(j) * exp(-lambda2(i,j)) / ResultFactorial(j)  RealMinSS]) );
                end; % if lambda
            end; % for j
        end; % for i

        MaxLogLike=max([LRi1 LRi2]);
        MinLogLike=min([LRi1 LRi2]);
        if MinLogLike<LRealMin
            ScaleAddition=LRealMin-(MinLogLike-MaxLogLike);
            LRi1=LRi1+ScaleAddition;
            LRi2=LRi2+ScaleAddition;
            MaxLogLike=max([LRi1 LRi2]);
        end;
            
        if exp(MaxLogLike)>0 % if max likelihood is computable
            LRi21=exp(LRi1-MaxLogLike);
            LRi22=exp(LRi2-MaxLogLike);
            Sum12=sum([LRi21 LRi22]);
            LRi21=LRi21/Sum12;
            LRi22=LRi22/Sum12;
        else
            fprintf('\nMax Likelihood =0 ! \n');
        end;
        
        
        %MuSigmaWeights= [ MuS0 SigmaS0 LRi2 ];
        Weights1(:,end)=LRi21;
        Weights2(:,end)=LRi22;
        
        
    end; % for Rem
    
%% Calculate Efficiency
    if Tmodel==1
        XO0=DGLM(Tbeta,m,2*PLD,model1,family,link,1e-4); % optimal design from True beta
        Eff(iter)=2*PLD/size(X0,1)*(det(InfoMtrxGLM(x2fx(X0,model1),Tbeta,family,link))/det(InfoMtrxGLM(x2fx(XO0,model1),Tbeta,family,link)))^(1/p1);
    else
        XO0=DGLM(Tbeta,m,2*PLD,model2,family,link,1e-4); % optimal design from True beta
        Eff(iter)=2*PLD/size(X0,1)*(det(InfoMtrxGLM(x2fx(X0,model2),Tbeta,family,link))/det(InfoMtrxGLM(x2fx(XO0,model2),Tbeta,family,link)))^(1/p2);
    end;
    MinEff=min(Eff);
    Q5Eff=quantile(Eff,0.05);
    MedEff=median(Eff);
    meanEff=mean(Eff);
    if iter>1, stdEff=std(Eff); else stdEff=0; end;
    fprintf('\n\n');
    fprintf('---------------------------------------------------\n');
    fprintf('Iter Last  Min    5%%     Med   Mean  std \n');
    fprintf(' %2.0f        %2.2f      %2.2f     %2.2f     %2.2f    %2.2f      %2.2f \n',[iter Eff(iter) MinEff Q5Eff MedEff meanEff stdEff]);
    fprintf('---------------------------------------------------\n');
    
        
end; % while iter<MaxIter

figure()
hist(Eff);
x=axis;
axis([0 1 x(3:4)]);

ElapsedTime=etime(clock,clo);
fprintf('Total: %2.1f hours, or %2.1f minutes per run \n\n',[ElapsedTime/3600 ElapsedTime/iter/60]);

