function [X,IM,D]=DGLM(beta,m,nruns,model,family,link, RP,NRP)
%
% DGLM Finds Local D-optimal Experimental Design for 
% multivariate Generalized Linear Models (GLM)
% using sequential Monte Carlo 
% and Federov's Exchange Algorithm as implemented in MATLAB statistical toolbox (required)
%
% Syntax:
% X=DGLM(beta,m)
% [X,IM,D]=DGLM(beta,m,nruns,model,family,link, RP,NRP)
%
% Description:
% X=DGLM(beta,m) is the minimal form, which finds local D-optimal design
%                                   for a binomial First-Order model with the logit link
%                                   m - a scalar value representing the number of variables,
%                                         all assumed to be constrained to the region [-1,1]
%                                   beta - a vector containing the  coefficient values
%                                   X is a matrix containing the output m  variables
%                                        values at the suggested experimental design
%
% Optional input variables are:
% nruns: number of runs in the final design; a scalar. If [], default=2^m
% model: a model in the form legitimate for the function x2fx
%                 options are 'linear' (default), 'interaction',
%                 'quadratic', 'purequadratic', or a MATRIX of model terms
%                  A matrix example: model=[0 0 0 ; 1 0 0; 0 1 0 ; 0 0 1] is 
%                      equivalent to model='linear', for m=3
% family: implemented families are: 'binomial' (default), and 'poisson'
% link: 'logit' (default for binomial), 'probit', 'cll', 'log' (default for poisson)
% RP: Required Precision (default 10^-2) of X matrix elements
% NRP: Number of random points around each base (best guess) of X from the
%           previous iteration step (default = 50)
%
% Values not supplied must be replaced by []
% for example X=DGLM([0 2 2 0],2,[],'interaction',[],'cll',10^-4)
%
% [X,IM]=DGLM(...)   returns, in the matrix IM, the information matrix of X
% [X, IM,D]=DGLM(...)   returns in the scalar D the D-criteria value for
%                                       the model; that is the p-th root of the determinant of 
%                                       the information matrix IM, divided  by the number of runs (nruns).
%                                       p is the length of the vector beta
%                                       (the number of model coefficients)
%
% This script was written by Hovav Dror, Tel-Aviv University, 2005
%
% Information can be found at:
% 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
%
%
% If you make changes to the script, or if you implement it
% to a (to be) published problem, please inform us:
% hovavdror@gmail.com   ;  dms@post.tau.ac.il 
 

t=clock;
 
if nargin<8, NRP=[];
    if nargin<7, RP=[];
        if nargin<6, link=[];
            if nargin<5, family=[];
                if nargin<4, model=[];
                    if nargin<3, nruns=[];
                    end, end, end, end, end, end
 
if nargin==0
        fprintf('----------------------------------------------------------------------------------\n');
        fprintf('[X,IM,D]=DGLM(beta,m,nruns,model,family,link, RP,NRP) \n');
        fprintf('Type "help DGLM" or see the function''s code\n');
        fprintf('for help and details \n');
%        fprintf('use: help DGLM \n');
        fprintf('----------------------------------------------------------------------------------\n');
elseif nargin<2
    fprintf('Error: Not enough arguments. You must supply at least values for beta, and m  \n');
else % ----------- if there are enough parameters, we can begin
    
    % ---------- substitute defaults ------------------
    if isempty(nruns), nruns=2^m;end;
    if isempty(model), model='linear'; end;  
    if isempty(family), family='binomial'; end;
    if isempty(link)
        if strcmp(family,'poisson'), link='log';
        else link='logit'; end, end
    if isempty(RP), RP=10^-2; end;
    if isempty(NRP), NRP=50; end;
        
    if isnumeric(model) % Decide what is p
            p=length(model);
     else
            p=length(x2fx(ones(1,m),model));
     end; % if isnumeric(model)
     if size(beta,1)==1, beta=beta'; end;
     if p > length(beta)
            fprintf('---------------------------------------------------------------------------------------------------\n');
            fprintf('Error: Size of Beta is shorter than the requested model ! \n');
            fprintf('---------------------------------------------------------------------------------------------------\n');
     else
            if p<length(beta)
                beta=beta(1:p);
                fprintf('Attention: beta is too long. Truncating it. \n');
           end;
    end;  % if p~=length(beta)
 
    
    % ----------- More initialization ---------------
        
    SearchR=1; % Initial Search Radius
    
    X=rand(nruns,m)*2-1; % Initial random design
    X=sortrows(X);
    if (nruns*NRP) < 1000  % ensure initial grid is big enough
        CS=rand(1000-nruns*NRP,m)*2-1;    % CS: Candidate Set
    else
        CS=[];
    end;
    % ---------------- Begin Search Loop------------------------
    while SearchR>RP % while the Required Precision not achieved
        CS=[CS; X];
        for i=1:nruns
            CS=[CS; ones(NRP,1)*X(i,:) + randn(NRP,m)*SearchR];% Add random points around X
        end;
        CS(CS < -1) = -1;% ensure the variables limits are [-1,1]
        CS(CS > 1) = 1 ;
        
        CSF=x2fx(CS,model); % Turn into a design matrix (without the weights, yet)
        WCSF=diag( sw(CSF,beta,family,link) ) * CSF; % Putting the (square root of ) weights
        
        % -------- Using the "regular" D-optimal algorithm ---------
        rows=candexch(WCSF,nruns,'display','off','maxiter',50,'init',WCSF(1:nruns,:)); %Choose the D-optimal nruns from C. return their row numbers.
        % WCSF is the Candidate Set, its first nruns rows are the initial design
 
        % ----------- Updating the Search Radius
        XN=sortrows(CS(rows,:));
        G=XN - X;
        DiffG=max(max(abs(G)));
        SearchR1=min( [DiffG   SearchR]);
        if (SearchR1==SearchR)  % if no improvement in Search radius, check if new design really better
            XN2=x2fx(XN,model);
            X2=x2fx(X,model);
            XN2=diag( sw(XN2,beta,family,link) )*XN2;
            X2=diag( sw(X2,beta,family,link) )*X2;
            if (det(XN2'*XN2) / det(X2'*X2)) < 1 % if no improvement
                XN=X; % retain old design
            end;
        end;
        SearchR=max([SearchR1 0.3*SearchR]);% Ensuring no RMradius=0 or too small change in RMradius
        CS=[]; % Empty the Candidate Set
       
        X=XN; % Updating the best design found so far
        
        if etime(clock,t)>10 % every 10 seconds show a sign of life with accuracy report
            t=clock;
            fprintf('Current Search Radius %g \n', SearchR);
        end;
        
    end; % While SearchR>RP
    
end; % if nargin==0
 
% -------------------- Calculate Results for IM and D -----------------------------------
if nargout>1
    FX=x2fx(X,model);
    IM=diag( sw(FX,beta,family,link) ) * FX;
    IM = IM' * IM;
    if nargout>2
        D=det(IM)^(1/p) / nruns; 
    end;
end; % if nargout>1
 
    
    % --------------- nested functions ----------------------------
    function y=sw(F,b,family,link) % Decision on (Sqrt of) Weights
        if strcmp(family, 'binomial')
            if strcmp(link,'logit')
                y=bsw(F,b);
            elseif strcmp(link,'cll') | strcmp(link,'comploglog');
                y=cllsw(F,b);
            elseif strcmp(link,'probit')
                y=probitsw(F,b);
            else
                fprintf('----------------------------------------------------\n');
                fprintf('Error: Unkown link function \n');
                fprintf('----------------------------------------------------\n');
           end; % if link
        elseif strcmp(family,'poisson')
            y=sqrt(exf(F,b)); % assuming log link automaticaly
        else
                fprintf('----------------------------------------------------\n');
                fprintf('Error: Unkown family function \n');
                fprintf('----------------------------------------------------\n');          
        end;
    end % sw
    
    % ---------------  Additional Support Functions -----------------------
    function y=bsw(F,b)
        ex=exf(F,b);
        y= sqrt (ex ./ (1+ex).^2 );
    end
    function y=cllsw(F,b)
            ex=exf(F,b);
            y=sqrt(ex.*ex.*exp(-ex) ./ (1.-exp(-ex)+realmin) );
    end
    function y=probitsw(F,b)
            P=normcdf((F*b)')';
            y= sqrt( 1./((P.*(1-P)+realmin)) .* exp(-0.5*(F*b).^2).^2 ./ (2*pi) );
            y(P.*(1-P)<realmin)=0;
    end
        
    function y=exf(F,b)
        y=exp(F*b);
    end
 % ------------------------------------------ end of support nested functions ------------------------------

 
end % DGLM function
 

 