%% Group Sequential Experimental Design
%
% This file 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
%
% At each step, the algorithm will suggest a locations for the 
% next group of observations, such as (for BatchSize=2):
%           Recommendation:   
%           -1.0000    0.6055   -1.0000   -1.0000
%            0.1606    0.0702         0   -0.0001
% You can either press 'enter' to accept, or enter a different location
% such as:  "enter"=accept, or enter value:    [-1 -1 -1 -1 ; 0 1 0 1]
%
% 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

% 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;
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));


    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

        % Recommend a value for the next group of observations
        fprintf('Recommendation:\n');
        disp(X0t(1:BatchSize,:));
        XChosen=input(' "enter"=accept, or enter value:'); % input where to take the next observations
        if isempty(XChosen), XChosen=X0t(1:BatchSize,:); end;

        X0=[X0 ; XChosen];
        ss=ss+BatchSize;
        F0=x2fx(X0);
        
        for i=(ss-BatchSize+1):ss
            RResultTemp=[];
            fprintf('For observation at: ');
            disp(X0(i,:));
            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;
        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

