%% Group Sequential Experimental Design
% This file is a Monte-Carlo examination, extending the file
% GroupSequential.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,too, is a source code for the algorithm found
% for group sequential experimental design for GLM
% It is using the example from section 4.4 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 Monter-Carlo iterations
iter=0; % current iteration

% In this example we assume that all variables are codes to the region [-1,1]
m=4; % Number of variables
p=5; % number of coefficiens

PLD=8; % Recommended number of points in a good Local D-optimal design

% Prior Distribution
betar=[ -3  3
    4  10
    5  11
    -6  0
    -2.5  3.5];
NN=10000;
NI=rand(NN,p);
beta=[];
for j=1:p
    beta(:,j)=betar(j,1)*ones(NN,1) + ( betar(j,2)-betar(j,1) ) * NI(:,j);
end;

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

BatchSize=2;
EL=48; % Length of the experiment (total number of observations)
L=round(EL/BatchSize); % Length of the experiment (total number of Batches)

% more initialization
N1=500; % N1 - number of vectors to calculate the weighted sum of determinants of information matrices
RealMin=realmin;
clo=clock;
Eff=[];
X0=[]; % This is the matrix with the experimental steps already done
RResult=[];
ss=0;

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


%% begin iterations
while iter<MaxIter
    iter=iter+1;
    %% "True" parameters
    NI=rand(1,p);
    Tbeta=[];
    for j=1:p
        Tbeta(:,j)=betar(j,1) + ( betar(j,2)-betar(j,1) ) * NI(:,j);
    end;
    X0=[]; % This is the matrix with the experimental steps already done
    RResult=[];
    ss=0; % ss=size(X0,1);
    Weights=ones(NN,1)/NN;

    for Rem=L:-1:1
        fprintf('%g ',Rem);
        % in group sequential we replace the weighted median with
        % a weighted K-means; here - simple as K-means on a weighted sample
        % from the prior
        C0 = beta( randsample(NN,NN,true,Weights) , : ); % start='Wsample'
        [idx,C]=kmeans(C0,BatchSize,'maxiter',500,'emptyaction','singleton','distance','cityblock','display','off');

        X0t=[];
        for i=1:BatchSize % Create a first guess by using the medians
            X0Ci{i}=GLMDAUG(X0,C(i,:)',m,ss+PLD); % PLD=8
            X0t=[X0t ; median(X0Ci{i})];
        end;
        X0t=[X0t ; X0];
        if (size(X0,1)+PLD)<p
            X0C=[];
            for i=1:BatchSize
                X0C=[X0C ; X0Ci{i}];
            end;
            [idx, CC]=kmeans(X0C,PLD);
            X0t=[X0t ; CC];
        end;
        % now, try to improve by replacements
        F0t=x2fx(X0t);
        ThisOne=0;
        for i2=1:N1
            ThisOnej=max(det(InfoMtrxGLM(F0t,beta(i2,:)',family,link)),RealMin);
            ThisOne=ThisOne + Weights(i2)*log(ThisOnej);
        end;
        BestSoFar= ThisOne;
        % Now we do a simplification of a mimic for an exchange algorithm
        for ii=1:2 % repeat twice - two cycles
            for i=1:BatchSize
                for j=1:PLD
                    ThisOne=0;
                    X0tR=X0t;
                    X0tR(i,:)=X0Ci{i}(j,:);
                    F0tR=x2fx(X0tR);
                    for i2=1:N1
                        ThisOnej=max(det(InfoMtrxGLM(F0tR,beta(i2,:)',family,link)),RealMin);
                        ThisOne=ThisOne + Weights(i2)*log(ThisOnej);
                    end;
                    if ThisOne>BestSoFar
                        BestSoFar=ThisOne;
                        X0t=X0tR;
                    end;
                end;
            end;
        end; % for ii

        X0=[X0 ; X0t(1:BatchSize,:)];
        ss=ss+BatchSize;
        F0=x2fx(X0);

        %% randomize results for the batch, according to the "true" model
        for i=(ss-BatchSize+1):ss
            Texpc=1/(1+exp(-(Tbeta*F0(i,:)')));
            RResult(end+1)= (rand(1)<Texpc); % Observation
        end;

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

    end; % for Rem

    %% Calculate Efficiency
    XO=DGLM(Tbeta,m,size(X0,1),model,family,link,1e-4); % optimal design from True beta
    Eff(iter)=(det(InfoMtrxGLM(x2fx(X0),Tbeta,family,link))/det(InfoMtrxGLM(x2fx(XO),Tbeta,family,link)))^(1/p);
    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]);

