function df = multistat(input_files_Y,input_files_sd,input_files_df, ...
     input_files_fwhm,X,contrast,output_file_base,which_stats, ...
     fwhm_varatio,niter)

%MULTISTAT 
%
% Combines effects (E) and their standard errors (S) using a linear mixed 
% effects model:     E = X b + e_fixed + e_random,     where
% b is a vector of unknown coefficients,
% e_fixed  is normal with mean zero, standard deviation S,
% e_random is normal with mean zero, standard deviation sigma (unknown).
% The model is fitted by REML using the EM algorithm with NITER iterations. 
% WARNING: 
% 1. Multistat is very slow if X is a matrix with more than one column, 
% since it loops over voxels, rather than doing calculations in parallel.
% 2. There may be some spurious extreme values in the results where the
% algorithm has not converged, often in voxels outside the brain.
% 
% DF = MULTISTAT( INPUT_FILES_EFFECT , INPUT_FILES_SDEFFECT , 
%           INPUT_FILES_DF [, INPUT_FILES_FWHM [, X [, 
%           [, CONTRAST [, OUTPUT_FILE_BASE [, WHICH_STATS 
%           [, FWHM_VARATIO [, NITER ]]]]]]]] )
% 
% INPUT_FILES_EFFECT is the input fmri effect files, the dependent variables,
% usually the _effect.mnc files, padded with extra blanks if necessary; 
% they will be removed before use.
% 
% INPUT_FILES_SDEFFECT is the input fmri sdeffect files, the standard
% deviations of the dependent variables, usually the _sdeffect.mnc files,
% padded with extra blanks if necessary; they will be removed before use.
% If INPUT_FILES_SDEFFECT=[], then INPUT_FILES_SDEFFECT is assumed to be
% zero for all voxels, INPUT_FILES_DF is set to Inf, and FWHM_VARATIO now
% smoothes the voxel sd. This allows multistat to duplicate DOT for 
% analysing PET data (with smoothed voxel sd, rather than pooled sd).
% 
% INPUT_FILES_DF is the row vector of degrees of freedom of the input files
% as printed out by fmrilm.m. If it is a scalar, it is repeated to fill.
%
% INPUT_FILES_FWHM is the fwhm in mm of the original fMRI data. It is only  
% used to calculate the degrees of freedom, printed out at the beginning. 
% Default is 6.
%
% X is the design matrix, whose rows are the files, and columns
% are the explanatory (independent) variables of interest. 
% Default is X=[1; 1; 1; ..1] which just averages the files. If the
% rank of X equals the number of files, e.g. if X is square, then 
% the random effects cannot be estinmated, but the fixed effects
% sd's can be used for the standard error. This is done very quickly.
% 
% CONTRAST is a row vector of contrasts for the statistic images.
% Default is [1 0 ... 0], i.e. it picks out the first column of X.
% 
% OUTPUT_FILE_BASE: Base for output statistics. Default is the first 
% INPUT_FILES_EFFECT name minus the .mnc extension, or minus the .mnc.gz
% extension if it is compressed.
%
% WHICH_STATS: logical row vector indicating output statistics by 1.
% Columns correspond to:
%   1: T statistic image, OUTPUT_FILE_BASE_tstat.mnc. The degrees 
%      of freedom is DF. Note that tstat=effect/sdeffect. 
%   2: effect image, OUTPUT_FILE_BASE_effect.mnc
%   3: standard deviation of the effect, OUTPUT_FILE_BASE_sdeffect.mnc. 
%   4: ratio of random effects standard deviation to fixed effects 
%      standard deviation, OUTPUT_FILE_BASE_sdratio.mnc.
% If the number of columns is less than 4 it is padded with zeros. 
% Default is 1. 
% 
% FWHM_VARATIO is the fwhm in mm of the Gaussian filter used to smooth the
% ratio of the random effects variance divided by the fixed effects variance.
%  - 0 will do no smoothing, and give a purely random effects analysis;
%  - Inf will do complete smoothing to a global ratio of one, giving a 
%    purely fixed effects analysis. 
% The higher the FWHM_VARATIO, the higher the ultimate degrees of
% freedom DF of the tstat image (printed out at the beginning of the run), 
% and the more sensitive the test. However too much smoothing will
% bias the results. The program prints and returns DF as its value; 
% if it is too low, press control c to cancel, and try increasing 
% FWHM_VARATIO. Default is 15. 
%
% NITER is the number of iterations of the EM algorithm. Default is 10.
%
% DF is the approximate (residual) degrees of freedom of the T statistics.

%############################################################################
% COPYRIGHT:   Copyright 2000 K.J. Worsley and C. Liao, 
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@math.mcgill.ca, liao@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that the above copyright
%              notice appear in all copies.  The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
%############################################################################

% Defaults:

if nargin < 4
   input_files_fwhm=6
end
if nargin < 5
   X=ones(size(input_files_Y,1),1)
end
if nargin < 6
   contrast=[1 zeros(1,size(X,2)-1)]
end
if nargin < 7 
   first_file=deblank(input_files_Y(1,:))
   sinf=size(first_file,2)
   if input_files_Y(1,sinf) == 'c'
      output_file_base=first_file(1,1:(sinf-4))
   else
      output_file_base=first_file(1,1:(sinf-7)) 
   end
end
if nargin < 8
   which_stats=1
end
if nargin < 9
   fwhm_varatio=15
end
if nargin < 10
   niter=10
end

% Open images:

n=size(input_files_Y,1)
in_handle_Y=zeros(1,n);
in_handle_sd=zeros(1,n);
for ifile=1:n
   ifile
   in_handle_Y(ifile)=openimage(deblank(input_files_Y(ifile,:)));
   if ~isempty(input_files_sd)
      in_handle_sd(ifile)=openimage(deblank(input_files_sd(ifile,:)));
   end
end

numslices=getimageinfo(in_handle_Y(1), 'NumSlices')
numys=getimageinfo(in_handle_Y(1), 'ImageHeight')
numxs=getimageinfo(in_handle_Y(1), 'ImageWidth')
numpix=numys*numxs;
Steps=getimageinfo(in_handle_Y(1), 'Steps')
numcolX=size(contrast,2)

% Open files for writing:

which_stats=[which_stats zeros(1,4-size(which_stats,2))]
if which_stats(1)  
   out_handle_tstat=newimage( ... 
      [deblank(output_file_base) '_tstat.mnc'], ...
      [0 numslices],input_files_Y(1,:),'float');
end
if which_stats(2)
   out_handle_effect=newimage( ... 
      [deblank(output_file_base) '_effect.mnc'], ...      
      [0 numslices],input_files_Y(1,:),'float');
end
if which_stats(3)
   out_handle_sdeffect=newimage( ...
      [deblank(output_file_base) '_sdeffect.mnc'], ...
      [0 numslices],input_files_Y(1,:),'float');
end  

X2=X'.^2;
p=rank(X)
Df=n-p
if Df>0
   % Degrees of freedom is greater than zero, so do mixed effects analysis:
   
   if which_stats(4)
      out_handle_sdratio=newimage( ...
         [deblank(output_file_base) '_sdratio.mnc'], ...
         [0 numslices],input_files_Y(1,:),'float');
   end   
   if isempty(input_files_sd)
      input_files_df=Inf
   end
   if length(input_files_df)==1
      input_files_df=ones(1,n)*input_files_df
   end
   df_fixed=sum(input_files_df)
   
   % Set-up for loop over voxels:
   
   Y=zeros(n,numpix);
   S=ones(n,numpix);
   varatio_vol=zeros(numpix, numslices);
   
   if (fwhm_varatio<Inf)
      
      % Smoothing varatio in slice is done using conv2 with a kernel ker_xy.
      % Smoothing betwen slices is done by straight matrix multiplication
      % by a toeplitz matrix K normalized so that the column sums are 1.
      
      if (fwhm_varatio>0)
         fwhm_x=fwhm_varatio/abs(Steps(1));
         ker_x=exp(-(-ceil(fwhm_x):ceil(fwhm_x)).^2*4*log(2)/fwhm_x^2);
         ker_x=ker_x/sum(ker_x);
         fwhm_y=fwhm_varatio/abs(Steps(2));
         ker_y=exp(-(-ceil(fwhm_y):ceil(fwhm_y)).^2*4*log(2)/fwhm_y^2);
         ker_y=ker_y/sum(ker_y);
         ker_xy=ker_x'*ker_y;
         fwhm_z=fwhm_varatio/abs(Steps(3));
         ker_z=exp(-(0:(numslices-1)).^2*4*log(2)/fwhm_z^2);
         K=toeplitz(ker_z);
         K=K./(ones(numslices)*K);
      else
         ker_xy=1;
         K=eye(numslices);
      end
      
      % Calculate degrees of freedom:
      
      df_indep=Df/(sum(sum(ker_xy.^2))*sum(sum(K.^2))/numslices)
      df_correl=Df*(2*(fwhm_varatio/input_files_fwhm)^2+1)^(3/2)
      df=round(1/(1/min(df_indep,df_correl)+1/df_fixed))
      
      % Now loop over voxels to get variance ratio, varatio:
      
      Sreduction=0.99
      for slice=1:numslices
         slice
         for ifile=1:n
            Y(ifile,:)=getimages(in_handle_Y(ifile),slice)';
         end
         sigma2=sum((Y-X*pinv(X)*Y).^2,1)/Df;
         if ~isempty(input_files_sd)
            for ifile=1:n
               S(ifile,:)=getimages(in_handle_sd(ifile),slice)'.^2;
            end
            minS=min(S)*Sreduction;
            Sm=S-ones(n,1)*minS;
            if size(X,2)==1
               % When X is a vector, calculations can be done in parallel:
               for iter=1:niter
                  W=1./(Sm+ones(n,1)*sigma2);
                  XWXinv=1./(X2*W);
                  betahat=XWXinv.*(X'*(W.*Y));
                  R=W.*(Y-X*betahat);
                  ptrS=p+sum(Sm.*W,1)-(X2*(Sm.*W.^2)).*XWXinv;
                  sigma2=(sigma2.*ptrS+(sigma2.^2).*sum(R.^2,1))/n; 
               end
               sigma2=sigma2-minS;
            else
               % Otherwise when X is a matrix, we have to loop over voxels:
               for pix=1:numpix
                  sigma2_pix=sigma2(pix);
                  Sm_pix=Sm(:,pix);
                  Y_pix=Y(:,pix);
                  for iter=1:niter
                     W=1./(Sm_pix+sigma2_pix);
                     Whalf=diag(sqrt(W));
                     WhalfX=Whalf*X;
                     pinvX=pinv(WhalfX);
                     WhalfY=Whalf*Y_pix;
                     betahat=pinvX*WhalfY;
                     R=WhalfY-WhalfX*betahat;
                     SW=diag(Sm_pix.*W);
                     ptrS=p+sum(Sm_pix.*W,1)-sum(sum((SW*WhalfX).*pinvX'));   
                     sigma2_pix=(sigma2_pix.*ptrS+ ...
                        (sigma2_pix.^2).*sum(W.*(R.^2),1))/n; 
                  end
                  sigma2(pix)=sigma2_pix-minS(pix);
               end         
            end
            
            varfix=input_files_df*S/df_fixed;   
            varatio=sigma2./(varfix+(varfix<=0)).*(varfix>0);
            varatio_slice=reshape(varatio,numxs,numys);
            varatio_slice=conv2(varatio_slice,ker_xy,'same');
            varatio_vol(:,slice)=reshape(varatio_slice,numpix,1);        
         end
      end
      
      % Smooth varatio in z direction:
      
      if fwhm_varatio>0 
         varatio_vol=varatio_vol*K;
      end
      
   else
      df=round(df_fixed)
   end
   
   effect_slice=zeros(numpix,1);
   sdeffect_slice=zeros(numpix,1);
   
   % Second loop over slices to get statistics:
   
   for slice=1:numslices
      slice            
      for ifile=1:n
         Y(ifile,:)=getimages(in_handle_Y(ifile),slice)';
      end
      if ~isempty(input_files_sd)
         for ifile=1:n
            S(ifile,:)=getimages(in_handle_sd(ifile),slice)'.^2;
         end
         varfix=input_files_df*S/df_fixed;   
         sigma2=varatio_vol(:,slice)'.*varfix;
         Sigma=S+ones(n,1)*sigma2;
         W=(Sigma>0)./(Sigma+(Sigma<=0));
      else
         sigma2=varatio_vol(:,slice)';
         W=1./(ones(n,1)*sigma2);
      end
      if size(X,2)==1
         % Do calculations in partallel:
         XWXinv=1./(X2*W);
         betahat=XWXinv.*(X'*(W.*Y));
         effect_slice=contrast*betahat';
         sdeffect_slice=contrast*sqrt(XWXinv)';
      else
         % Loop over pixels:
         for pix=1:numpix
            Whalf=diag(sqrt(W(:,pix)));
            pinvX=pinv(Whalf*X);
            betahat=pinvX*Whalf*Y(:,pix);
            effect_slice(pix)=contrast*betahat;
            sdeffect_slice(pix)=norm(contrast*pinvX);
         end  
      end
      tstat_slice=effect_slice./(sdeffect_slice+(sdeffect_slice<=0)) ...
         .*(sdeffect_slice>0);
      sdratio_slice=sqrt(varatio_vol(:,slice)+1);
      
      if which_stats(1)
         putimages(out_handle_tstat,tstat_slice,slice);
      end
      if which_stats(2)
         putimages(out_handle_effect,effect_slice,slice);
      end
      if which_stats(3)
         putimages(out_handle_sdeffect,sdeffect_slice,slice);
      end
      if which_stats(4)
         putimages(out_handle_sdratio,sdratio_slice,slice);
      end
   end
   if which_stats(4)
      closeimage(out_handle_sdratio);
   end
   
else
   % If degrees of freedom is zero, estimate effects by least squares,
   % and use the standard errors to estimate the sdeffect.
  
   Y=zeros(n,numpix);
   S=ones(n,numpix);
   cXinv=contrast*pinv(X);
   cXinv2=cXinv.^2;    
   for slice=1:numslices
      slice            
      for ifile=1:n
         Y(ifile,:)=getimages(in_handle_Y(ifile),slice)';
      end
      if ~isempty(input_files_sd)
         for ifile=1:n
            S(ifile,:)=getimages(in_handle_sd(ifile),slice)'.^2;
         end
      end
      effect_slice=(cXinv*Y)';
      sdeffect_slice=sqrt(cXinv2*S)';
      tstat_slice=effect_slice./(sdeffect_slice+(sdeffect_slice<=0)) ...
         .*(sdeffect_slice>0);
      if which_stats(1)
         putimages(out_handle_tstat,tstat_slice,slice);
      end
      if which_stats(2)
         putimages(out_handle_effect,effect_slice,slice);
      end
      if which_stats(3)
         putimages(out_handle_sdeffect,sdeffect_slice,slice);
      end
   end
   
end

% Tidy up:

if which_stats(1)
   closeimage(out_handle_tstat);
end
if which_stats(2)
   closeimage(out_handle_effect);
end
if which_stats(3)
   closeimage(out_handle_sdeffect);
end
for ifile=1:n
   closeimage(in_handle_Y(ifile));
end
if ~isempty(input_files_sd)
   for ifile=1:n
      closeimage(in_handle_sd(ifile));
   end
end

% End.