% This file refers to section 4 of:
% Technical Report RP-SOR-0601, Robust Experimental Design for Multivariate Generalized Linear Models, Hovav A. Dror and David M. Steinberg, January 2006.
% Available at http://www.math.tau.ac.il/~dms/GLM_Design

% It is a work through for finding a robust design using clustering
% to a Logisitic model example, 4-dimensional, with a uniform prior
% The model is the crystallography experiment described in section 6
% (example 2) of
% Woods, D.C., Lewis, S.M., Eccleston, J.A. and Russell, K.G. (2005) Designs for generalized linear models with several variables and model uncertainty, Technometrics, to appear.

% You can run the script and see its stages

% The Parameter Space
BetaV=[];
betar=[ -3  3
                  4  10
                  5  11
                  -6  0
                  -2.5  3.5];
 p=size(betar,1);
 m=4; % number of variables
 model='linear';
 family='binomial';
 link='logit';
 
 plot([0 0 1],[1 0 0],'k'); % 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

 % ------------------------------------------------------------------------
 % This step is optional.
 % Its purpose is to choose the number of support points
 % in each Local D-optimal design
 fprintf('----------------------------------------------------------------------------------------------------\n');
 fprintf('\nThis is an optional step \n');
fprintf('It helps choosing the number of support \n');
fprintf('points for a Local D-optimal Design \n');
t1=0; % time for first step
A=input('Choose "1" to use the option, \nor press "Enter" to skip: ');
if ~isempty(A)
    t1=clock;
    fprintf('The number of support points tested now is:');
     BetaMean= mean(betar,2); % The centroid of the beta space
    NSP=p:2^m;   % Sequence of support points values to check
    SPE=[];          % Vector with the Efficiency for each NSP value
    i=0;
    for SP=NSP
        i=i+1;
        fprintf(' %g',SP);
        [X,IM,SPE(i)]=DGLM(BetaMean,m,SP,model,family,link,10^-4);
    end; % for SP
    fprintf(' \n \n');
    plot(NSP,SPE/max(SPE),'.-k');
    axis([min(NSP) max(NSP) 0 1]);
    grid on;
    t1=etime(clock,t1); % Time to compute
    fprintf('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n');
    fprintf('(Finished optional part in %1.0f seconds)\n',t1);
    xlabel('Number of Support Points in Beta''s Centroid Local D-Optimal Design');
    ylabel('Relative Efficiency');
    title('Choose the proper number of support points for a Local Design');

    LSP=input('Using this graph, choose the number of \nSupport Points for Local Designs (default=8): ');
    if isempty(LSP),  LSP=8; end;
else % if skipped
    LSP=input('Choose the number of Support Points \nfor Local Designs (default=8): ');
    if isempty(LSP),  LSP=8; end;
end; % if use the option to choose LSP by graph

% -------------------------------------------------------------------------
% First Step: Create Local Designs all over the Parameter Space
N=1000; % Number of vectors from the parameter sapce
% Use a Niederreiter sequence of proper size
%NI=Nied(N,p); % If function exist: Niederreiter sequence of the proper size
%otherwise:
t2=clock; % timing the Local Database creation
NI=load('Nied1000.m');
NI=NI(1:N,1:p);
% Now transform it according to the values of betar:
for j=1:p        
    BetaV(:,j)=betar(j,1)*ones(N,1) + ( betar(j,2)-betar(j,1) ) * NI(:,j);
end;

% Now Search for the local Designs, and put them in the matrix X
X=[]; % Designs
IM=[]; %Information Matrices
IMD=[]; % D-Criteria value of the regression matrices
fprintf('\n----------------------------------------------------------------------------------------------------\n');
fprintf('Next: Finding 100 Local D-Optimal designs \n');
N2=input('press "enter" to continue, \n(or choose a different number of Local Designs): ');
if ~isempty(N2), N=N2; else N=100; end;
fprintf('Now Finding Design #:');
for i=1:N
    fprintf(' %g',i);
    [Xi,IMi, Di]=DGLM(BetaV(i,:),m,LSP,model,family,link);
    X=[ X ; Xi ] ;
    %IM{i}=IMi;
    IMD=[IMD ; Di];
end;
X=X-sign(X).*rand(size(X)).*0.0001; % Jittering
t2=etime(clock,t2);
fprintf('\n- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -');
fprintf('\nDone! Database Creation time %2.0f seconds.\n \n',t2);

% -----------------------------------------------------------------------------------------
% Next Step is choosing the number of Support points for the compromise
% design (What should be the K of the K-means clustering?)
t3=clock; %timing this step
fprintf('----------------------------------------------------------------------------------------------------\n');
fprintf('Time to choose "K" for K-means Clustering. \n');
fprintf('Repeateing the procedure only once for each K \n');
fprintf('(may cause the graph to be a bit jumpy) \n');
input('press "enter" to continue');
fprintf('Now Trying K= ');
KS=[p:2^m 2^m:2:2^(m+1) 2^(m+1):4:2^(m+2)]; % sequence of tested K values
KX=[]; % Design matrix by clustering (for different K's);
Eff=[];
MedEff=[];
MinEff=[];
i=0;
for K=KS
    i=i+1;
    fprintf(' %g',K);
    [IDX,KX,s] =kmeans(X,K,'maxiter',500,'emptyaction','singleton','distance','cityblock');
    for j=1:N
        Eff(i,j)=( det( InfoMtrxGLM(x2fx(KX,model),BetaV(j,:)',family,link) )^(1/p) ) / K / IMD(j);
    end; % for j
    MedEff(i)=median(Eff(i,:));
    MinEff(i)=min(Eff(i,:));
end; % for K
t3=etime(clock,t3);
fprintf('\n- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -');
fprintf('\nDone in %2.0f seconds (All phases togeter %2.0f seconds)\n',[t3 t1+t2+t3]);
plot(KS,MedEff,'.-k');
hold on
plot(KS,MinEff,'s:k','markersize',3);
hold off
xlabel('Number of Support Points');
ylabel('Approximate Efficiency');
legend('Median Efficiency','Minimum Efficiency','Location','SouthEast')
title('Choose "K" - the number of support points for the robust design');

K=input('Using this graph, choose the number of Support Points, K, \nfor The Robust Design (default=48): ');
if isempty(K),  K=48; end;

% -----------------------------------------------------------------------------------------
% Now use K-means clustering to produce the Robust Design
fprintf('\n----------------------------------------------------------------------------------------------------\n');
fprintf('Last step - repeat the clustering with the chosen K (%g), \nand choose the Robust Design\n',K);
%Nrepeat=N; % Number of times the clustering will be performed
Nrepeat=40; % Number of times the clustering will be performed
Nrepeats=input(sprintf('Choose the number of clustering repetitions \n(or press "enter" to use %g repetitions): ',Nrepeat));
if isempty(Nrepeats), Nrepeats=Nrepeat; end;
fprintf('Number of Clustering repetition processed: ');
t4=clock; % timing the step
BDV = -inf ; % Best D-criteria value found
for i=1:Nrepeats
    fprintf(' %g',i);
    [IDX,KX,s] =kmeans(X,K,'maxiter',500,'emptyaction','singleton','distance','cityblock');
    DV=1; % D-criteria representing value for this clustering
    for j=1:N
        DV=DV + log( det(InfoMtrxGLM(x2fx(KX,model),BetaV(j,:)',family,link) ) );
    end; % for j
    if DV>BDV
        BX=KX;
        BDV=DV;
    end; % if improved the design update BX - the best design found so far
end; % for i (Nrpeats)
t4=etime(clock,t4);
fprintf('\n- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n');
fprintf('Done! This step took %2.0f seconds. Overall the process took %2.0f seconds',[t4 t1+t2+t3+t4]);
fprintf('\nThe best design found, the matrix BX, is ready to be used\n');
fprintf('----------------------------------------------------------------------------------------------------\n');


% -------------------------------------------------------------------------
% Last, without computing "many" local designs, we can't produce
% a good histogram of efficiencies.
% The following is a rough approximation, based on the N designs created
% before
input('Press "enter" for a plot representing a rough estimation \nof the efficiency histogram\n');
HEff=[];
for i=1:N
     HEff = [ HEff  ( det( InfoMtrxGLM(x2fx(BX,model),BetaV(i,:)',family,link) )^(1/p) ) / K / IMD(i) ];
end;
hist(HEff);
x=axis;
axis([0 1 x(3:4)]);
xlabel('Approximate Efficiency');
title(sprintf('Efficiency Estimation based on the same %g local designs',N));
