function [df, p]=fmrilm(input_file,output_file_base,X_cache, ...
   contrast,exclude,which_stats,fwhm_rho,n_poly,reference_delay)

%FMRILM
%
% Reads in an fMRI minc file and writes out a T statistic minc file. 
% The method is based on linear models with correlated errors.
% 
% [DF,P] = FMRILM(INPUT_FILE, OUTPUT_FILE_BASE, X_CACHE, CONTRAST [,
% EXCLUDE [,WHICH_STATS [,FWHM_RHO [,N_POLY [,REFERENCE_DELAY ]]]]] )
% 
% INPUT_FILE is the input fmri file, the dependent variable.
% Scans must be equally spaced in time, unless FWHM_RHO=Inf (see below).
% 
% OUTPUT_FILE_BASE: matrix whose rows are the base for output statistics,
% one base for each row of CONTRAST, padded with extra blanks if 
% necessary; they will be removed before use.
%
% X_CACHE: The design matrices for all slices usually from fmridesign.m. 
% Rows are the frames, columns are all the regressor variables 
% for slice 1, then slice 2, i.e. with slices running slowest. If X_CACHE  
% only contains variables for just one slice, it is applied to all slices.
%
% CONTRAST is a matrix whose rows are contrasts for the statistic images,
% with row length equal to the number regressor variables in X_CACHE
% plus N_POLY+1; the extra values allow you to set contrasts
% in the drift terms as well as the regressor variables.
% 
% EXCLUDE is a list of frames that should be excluded from the
% analysis. This must be used with Siemens EPI scans to remove the
% first few frames, which do not represent steady-state images.
% Default is [1].
% 
% WHICH_STATS: logical matrix indicating output statistics by 1.
% Rows correspond to rows of CONTRAST, 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: F-statistic image, OUTPUT_FILE_BASE_Fstat.mnc, for testing
%      all rows of CONTRAST simultaneously. The degrees 
%      of freedom are [DF, P].
%   5: the AR parameter, OUTPUT_FILE_BASE_rho.mnc.
%   6: the residuals from the model, OUTPUT_FILE_BASE_resid.mnc,
%      only for the non-excluded frames (warning: uses lots of memory). 
%   7: the whitened residuals from the model, OUTPUT_FILE_BASE_wresid.mnc,
%      only for the non-excluded frames (warning: uses lots of memory). 
% If WHICH_STATS to a row vector of length 7, then it is used for all
% contrasts; if the number of columns is less than 7 it is padded
% with zeros. Default is 1. 
%   
% FWHM_RHO is the fwhm in mm of a 3D Gaussian kernel used to smooth the
% autocorrelation. This is NOT the fwhm of the image data as used
% by tstat_threshold, but should be fairly large. Setting it to Inf
% smooths the autocorrelation to 0, i.e. it assumes the scans
% are uncorrelated (useful for TR>10 seconds). Default is 15. 
% 
% N_POLY is the order of polynomial trend terms to remove
% (N_POLY>=0); =0 will remove just the constant term and no trends.
% -1 will not remove anything, in which case the design matrix is
% completely determined by X. Default is 3.
%
% REFERENCE_DELAY: If the first element >0 then delays are estimated
% as well as effects. We start with a reference hrf h(t) close
% to the true hrf, e.g. the Glover (1999) hrf used by FMRIDESIGN.
% The hrf for each variable is then modeled as  h( t / c ) / c where 
% c = DELAY / REFERENCE_DELAY, i.e. the time is re-scaled by c.
% The program then estimates contrasts in the c's (specified by 
% CONTRAST; contrasts for the polynomial drift terms are ignored),
% their standard deviations, and T statistics for testing c=1 
% (no change). Output is DELAY = c * REFERENCE_DELAY in 
% OUTPUT_FILE_BASE_delay.mnc, sd is in OUTPUT_FILE_BASE_sddelay.mnc,
% and t statistic in OUTPUT_FILE_BASE_tdelay.mnc. Note that
% REFERENCE_DELAY is only used to multiply c and its sd for
% output. If REFERENCE_DELAY is set equal to the first 
% element of HRF_PARAMETERS (5.4 seconds by default) then DELAY is 
% the estimated time to the first peak of the hrf (in seconds).
% To do the estimation, the stimuli must be convolved with first 
% and second derivatives of h with respect to log(c) at c=1
% and added to X_CACHE. FMRIDESIGN with the last element
% of HRF_PARAMETERS set to 2 does this for the Glover (1999) hrf.
% If REFERENCE_DELAY is a scalar, it is used for all variables.
% Note that F statistics are not yet available with this option.
% Default is 0, i.e. only effects are estimated.
%  
% DF is the (residual) degrees of freedom of the statistics.
%
% P is the degrees of freedom of the set of CONTRASTs.

%############################################################################
% 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 < 5
   exclude=[1]
end
if nargin < 6
   which_stats=1
end
if nargin < 7
   fwhm_rho= 15
end
if nargin < 8
   n_poly=3
end
if nargin < 9
   reference_delay=0.0
end

isdelay=(reference_delay>0)

% Open images:

in_handle=openimage(input_file);
numframes=getimageinfo(in_handle, 'NumFrames')
numslices=getimageinfo(in_handle, 'NumSlices')
numys=getimageinfo(in_handle, 'ImageHeight')
numxs=getimageinfo(in_handle, 'ImageWidth')
numpix=numxs*numys;
Steps=getimageinfo(in_handle, 'Steps')

% Keep time points that are not excluded:

allpts = 1:numframes;
allpts(exclude) = zeros(1,length(exclude));
keep = allpts( find( allpts ) );
n=length(keep);

% Create polynomial trends:

if n_poly>=0 
   trend=((2*keep-(max(keep)+min(keep)))./(max(keep)-min(keep)))';
   Trend=(trend*ones(1,n_poly+1)).^(ones(n,1)*(0:n_poly));
else
   Trend=[];
end 
% Sizes of contrasts:

numcontrasts=size(contrast,1)
numcolX_cache0=size(contrast,2)-n_poly-1
if isdelay
   numcolX_cache1=numcolX_cache0*2
   numcolX_cache2=numcolX_cache0*3
   colsX0=1:numcolX_cache0
   colsX1=numcolX_cache0+1:numcolX_cache1
   colsX01=1:numcolX_cache1
   Contrast=[contrast(:,colsX0) ... 
         zeros(numcontrasts,numcolX_cache0) ...
         contrast(:,numcolX_cache0+(1:n_poly+1))]
   contr=contrast(:,colsX0)
   contr2=[contr contr]
   is_single_X=sum(contr==0,2)==(numcolX_cache0-1)
   findX=zeros(numcontrasts,1);
   for i=1:numcontrasts
      findX(i)=min(find(contr(i,:)~=0));
   end
   findX=findX+numcolX_cache0
   if numslices*numcolX_cache2~=size(X_cache,2) & ...
         numcolX_cache2~=size(X_cache,2)
      fprintf(['Error: ( number of columns of contrast - n_poly - 1 ) * 3 = ', ...
            num2str(numcolX_cache2) '\n']);
      fprintf(['       multiplied by the number of slices = ', ...
            num2str(numslices) '\n']);
      fprintf(['       must be equal to the number of columns of X_cache = ', ...
            num2str(size(X_cache,2)) '\n']);
      fprintf('       (or the ( number of columns of contrast - n_poly - 1 ) * 3 \n');
      fprintf('       must be equal to the number of columns of X_cache).\n');
      return
   end
else
   numcolX_cache1=numcolX_cache0
   numcolX_cache2=numcolX_cache0
   Contrast=contrast
   if numslices*numcolX_cache2~=size(X_cache,2) & ...
         numcolX_cache2~=size(X_cache,2)
      fprintf(['Error: number of columns of contrast - n_poly - 1 = ', ...
            num2str(numcolX_cache2) '\n']);
      fprintf(['       multiplied by the number of slices = ', ...
            num2str(numslices) '\n']);
      fprintf(['       must be equal to the number of columns of X_cache = ', ...
            num2str(size(X_cache,2)) '\n']);
      fprintf('       (or the number of columns of contrast - n_poly - 1 \n');
      fprintf('       must be equal to the number of columns of X_cache).\n');
      return
   end
end
numcolX=numcolX_cache1+n_poly+1
p=rank(Contrast)

% Check for estimability:

tolerance=0.0000001;
sumdf=0;
for slice = 1:numslices
   if slice*numcolX_cache2 <= size(X_cache,2)
      whichcolsX=(1:numcolX_cache1)+(slice-1)*numcolX_cache2;
      X=[X_cache(keep,whichcolsX) Trend];
      df=n-rank(X);
      sumdf=sumdf+df;
      NullSpaceX=null(X);
      Cmhalf=diag(1./sqrt(diag(Contrast*Contrast')));
      ContrastNullSpaceX=Cmhalf*Contrast*NullSpaceX;
      nonest=sum(abs(ContrastNullSpaceX)>tolerance,2);
      if sum(nonest)>0
         fprintf(['Error: the following contrasts are nonestimable in slice ' ...
               num2str(slice) ':']);
         RowsNonEstContrast=find(nonest>0)
         NonEstContrast=Contrast(RowsNonEstContrast,:)
         NullSpaceX
         ContrastNullSpaceX
         return
      end
   end
end
df=round(sumdf/numslices)
   
% Open files for output:

if size(which_stats,1)
   which_stats=kron(which_stats,ones(numcontrasts,1));
end
which_stats=[which_stats zeros(numcontrasts,7-size(which_stats,2))]

if isdelay
   for i=1:numcontrasts
      if which_stats(i,1) 
         out_handle_tdelay(i)=newimage( ... 
            [deblank(output_file_base(i,:)) '_tdelay.mnc'], ...
            [0 numslices],input_file,'float');
      end
      if which_stats(i,2)
         out_handle_delay(i)=newimage( ... 
            [deblank(output_file_base(i,:)) '_delay.mnc'], ...
            [0 numslices],input_file,'float');
      end
      if which_stats(i,3)
         out_handle_sddelay(i)=newimage( ...
            [deblank(output_file_base(i,:)) '_sddelay.mnc'], ...
            [0 numslices],input_file,'float');
      end
   end
end

for i=1:numcontrasts
   if which_stats(i,1) 
      out_handle_tstat(i)=newimage( ... 
         [deblank(output_file_base(i,:)) '_tstat.mnc'], ...
         [0 numslices],input_file,'float');
   end
   if which_stats(i,2)
      out_handle_effect(i)=newimage( ... 
         [deblank(output_file_base(i,:)) '_effect.mnc'], ...
         [0 numslices],input_file,'float');
   end
   if which_stats(i,3)
      out_handle_sdeffect(i)=newimage( ...
         [deblank(output_file_base(i,:)) '_sdeffect.mnc'], ...
         [0 numslices],input_file,'float');
   end
end

if which_stats(1,4)
   out_handle_Fstat=newimage( ...
      [deblank(output_file_base(1,:)) '_Fstat.mnc'], ...
      [0 numslices],input_file,'float');
end
if which_stats(1,6)
   out_handle_resid=newimage( ...
      [deblank(output_file_base(1,:)) '_resid.mnc'], ...
      [n numslices],input_file);
end
if which_stats(1,7)
   out_handle_wresid=newimage( ...
      [deblank(output_file_base(1,:)) '_wresid.mnc'], ...
      [n numslices],input_file);
end

rho_vol=zeros(numpix, numslices);

indk1=((keep(2:n)-keep(1:n-1))==1);
k1=find(indk1)+1;

if fwhm_rho < Inf 
   
   if which_stats(1,5)
      out_handle_rho=newimage( ...
         [deblank(output_file_base(1,:)) '_rho.mnc'], ...
         [0 numslices],input_file,'float');
   end
   
   % Smoothing rho in slice is done using conv2 with a kernel ker_xy.
   
   if fwhm_rho>0
      fwhm_x=fwhm_rho/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_rho/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;
   else
      ker_xy=1;
   end
   
   % First loop over slices, then pixels, to get the AR parameter:
   
   D1=diag(indk1,1)+diag(indk1,-1);
    
   for slice = 1:numslices
      slice
      if slice*numcolX_cache2 <= size(X_cache,2)
         whichcolsX=(1:numcolX_cache1)+(slice-1)*numcolX_cache2;
         X=[X_cache(keep,whichcolsX) Trend];
         pinvX=pinv(X);
         
         % Preliminary calculations for unbiased estimates of autocorrelation:
         
         R=eye(n)-X*pinvX;
         M11=trace(R);
         M12=trace(R*D1);
         M21=M12/2;
         M22=trace(R*D1*R*D1)/2;
         M=[M11 M12; M21 M22];
         invM=inv(M);
         %invM=eye(2); this will undo the correction.
      end
      
      % Loop over pixels in slice to find rho: 
      
      rho_slice=zeros(numxs,numys);
      for y=1:numys 
         img=getimages(in_handle,slice,1:numframes,[], y);
         Y=img(:,keep)';
         betahat_ls=pinvX*Y;
         resid=Y-X*betahat_ls;
         Cov0=sum(resid.*resid,1);
         Cov1=sum(resid(k1,:).*resid(k1-1,:),1);
         Covadj=invM*[Cov0; Cov1];
         ind=find(Covadj(1,:)>0);
         rho_slice(ind,y)=(Covadj(2,ind)./Covadj(1,ind))';
      end
      
      rho_vol(:,slice)=reshape(conv2(rho_slice,ker_xy,'same'),numpix,1);   
      
   end
   
   % Smoothing rho betwen slices is done by straight matrix multiplication
   % by a toeplitz matrix K normalized so that the column sums are 1.
   
   if fwhm_rho>0
      fwhm_z=fwhm_rho/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);
      rho_vol=rho_vol*K;
   end
   
   if which_stats(1,5)
      for slice=1:numslices
         putimages(out_handle_rho,rho_vol(:,slice),slice);
      end
   end
   
end

tstat_slice=zeros(numpix,numcontrasts);
effect_slice=zeros(numpix,numcontrasts);
sdeffect_slice=zeros(numpix,numcontrasts);
if isdelay
   tdelay_slice=zeros(numpix,numcontrasts);
   delay_slice=zeros(numpix,numcontrasts);
   sddelay_slice=zeros(numpix,numcontrasts);
end
if which_stats(1,4)
   Fstat_slice=zeros(numpix,1);
end
if which_stats(1,6)
   resid_slice=zeros(numpix,n);
end
if which_stats(1,7)
   wresid_slice=zeros(numpix,n);
end

% Set up caches indexed by rho:

ntable=100;
ntable2=2*ntable-1;
Xstar_cache=zeros(n,numcolX*ntable2);
pinvXstar_cache=zeros(n,numcolX*ntable2);
if isdelay
   D_cache=zeros(numcolX_cache0,ntable2);
   V_cache=zeros(numcolX_cache1,numcolX_cache1*ntable2);
end
sdeffect_mult_cache=zeros(numcontrasts,ntable2);
if which_stats(1,4)
   vareffinv_cache=zeros(numcontrasts,numcontrasts*ntable2);
end

% Second loop over voxels to get statistics:

for slice=1:numslices
   slice
   if slice*numcolX_cache2 <= size(X_cache,2)
      whichcolsX=(1:numcolX_cache1)+(slice-1)*numcolX_cache2;
      X=[X_cache(keep,whichcolsX) Trend];
      Xstar=X;
      if isdelay
         whichcolsX2=(numcolX_cache1+1:numcolX_cache2)+ ...
            (slice-1)*numcolX_cache2;
         X2=X_cache(keep,whichcolsX2);
         X2star=X2;
      end
      rho_is_in_cache=zeros(1,ntable2);
      df=n-rank(X);
   end
  
   for y=1:numys 
      img = getimages (in_handle, slice, 1:numframes, [], y);
      Y=img(:,keep)';
      
      for x=1:numxs
         pix=x+(y-1)*numxs;
         rho_pix=rho_vol(pix,slice);
         factor=1./sqrt(1-rho_pix.*rho_pix);
         Ystar=Y(:,x);
         Ystar(k1,1)=(Y(k1,x)-rho_pix.*Y(k1-1,x)).*factor;
         
         % Look in cache for rho; if it's not there, calculate pinv, 
         % sd effect multiplier and Xstar and put them in the caches:
         
         irho=round((rho_pix+1)*ntable);
         if irho<1 
            irho=1;
         end
         if irho>ntable2
            irho=ntable2;
         end
         colsX=(1:numcolX)+(irho-1)*numcolX;
         colsv=(1:numcontrasts)+(irho-1)*numcontrasts;
         colsV=(1:numcolX_cache1)+(irho-1)*numcolX_cache1;
         if ~rho_is_in_cache(irho)
            
            % either calculate new values for the caches:
            
            rho_pix=(irho-ntable)/ntable;
            factor=1./sqrt(1-rho_pix.*rho_pix);
            Xstar(k1,:)=(X(k1,:)-rho_pix.*X(k1-1,:)).*factor;
            Xstar_cache(:,colsX)=Xstar;
            pinvXstar=pinv(Xstar);
            pinvXstar_cache(:,colsX)=pinvXstar';
            if isdelay
               X2star(k1,:)=(X2(k1,:)-rho_pix.*X2(k1-1,:)).*factor;
               D=sum(pinvXstar(colsX1,:).*X2star',2);
               D_cache(:,irho)=D;
               pinvXstar01=pinvXstar(colsX01,:);
               V=pinvXstar01*pinvXstar01';
               V_cache(:,colsV)=V;
            end
            CpinvX=Contrast*pinvXstar;
            vareffect=CpinvX*CpinvX';
            sdeffect_mult=sqrt(diag(vareffect)/df);
            sdeffect_mult_cache(:,irho)=sdeffect_mult;
            if which_stats(1,4)
               vareffinv=pinv(vareffect)*df;
               vareffinv_cache(:,colsv)=vareffinv;
            end
            rho_is_in_cache(irho)=1;
            
         else
            
            % or use old values already in caches:
            
            Xstar=Xstar_cache(:,colsX);
            pinvXstar=pinvXstar_cache(:,colsX)';
            if isdelay
               D=D_cache(:,irho);
               V=V_cache(:,colsV);
            end
            sdeffect_mult=sdeffect_mult_cache(:,irho);
            if which_stats(1,4)
               vareffinv=vareffinv_cache(:,colsv);
            end
            
         end
         
         % Calculate statistics:
         
         betahat=pinvXstar*Ystar;
         resid=Ystar-Xstar*betahat;
         if which_stats(1,6)
            resid_slice(pix,:)=(Y(:,x)-X*betahat)';
         end
         if which_stats(1,7)
            wresid_slice(pix,:)=resid';
         end
         SSE=sum(resid.^2);
         if SSE>0
            
            if isdelay
               bhat=betahat(colsX0);
               dbhat=betahat(colsX1);
               diagV=diag(V);
               sd=sqrt(SSE/df);
               T2=(bhat.^2)./(SSE/df*diagV(colsX0));
               invbhat=(bhat~=0)./(bhat+(bhat==0));
               dhat_L=dbhat.*invbhat;
               dhat_LU=dhat_L-0.5*dhat_L.^2.*D;
               dhat=dhat_LU./(T2+1).*T2;
               delta=contr*dhat;
               gdot=[(1-(1-D.*dhat_L).*T2).*T2./(T2+1).^2 ...
                  .*dhat_L.*invbhat; ... 
                  (1-D.*dhat_L)./(T2+1).*T2.*invbhat]; 
               gVg=(gdot*gdot').*V; 
               cVc=contr2*gVg*contr2';
               sddelta=sqrt(diag(cVc))*sd; 
               expdhat=exp(dhat);
               delay=contr*expdhat*reference_delay;
               expdhat2=[expdhat; expdhat];
               Vdelay=contr2*((expdhat2*expdhat2').*gVg)*contr2';
               sddelay=sqrt(diag(Vdelay))*sd*reference_delay; 
               delay_slice(pix,:)=delay';
               sddelay_slice(pix,:)=sddelay';
               tdelay_slice(pix,:)= ( is_single_X.* ...                  
                  betahat(findX)./sqrt(diagV(findX))/sd ...
                  +(~is_single_X).*delta./sddelta )';
            end
            effect=Contrast*betahat;
            sdeffect=sqrt(SSE)*sdeffect_mult;
            effect_slice(pix,:)=effect';
            sdeffect_slice(pix,:)=sdeffect';
            tstat_slice(pix,:)=(effect./sdeffect)';
            if which_stats(1,4)
               Fstat_slice(pix,1)=effect'*vareffinv*effect/p/SSE;
            end
            
         else
            if isdelay
               tdelay_slice(pix,:)=zeros(1,numcontrasts);
               delay_slice(pix,:)=zeros(1,numcontrasts);
               sddelay_slice(pix,:)=zeros(1,numcontrasts);
            end   
            tstat_slice(pix,:)=zeros(1,numcontrasts);
            effect_slice(pix,:)=zeros(1,numcontrasts);
            sdeffect_slice(pix,:)=zeros(1,numcontrasts);
            if which_stats(1,4)
               Fstat_slice(pix,1)=0.0;
            end
         end
      end
   end
   
   % Write out results:
      
   for i=1:numcontrasts
      if isdelay
         if which_stats(i,1)
            putimages(out_handle_tdelay(i),tdelay_slice(:,i),slice);
         end
         if which_stats(i,2)
            putimages(out_handle_delay(i),delay_slice(:,i),slice);
         end
         if which_stats(i,3)
            putimages(out_handle_sddelay(i),sddelay_slice(:,i),slice);
         end
      end
      if which_stats(i,1)
         putimages(out_handle_tstat(i),tstat_slice(:,i),slice);
      end
      if which_stats(i,2)
         putimages(out_handle_effect(i),effect_slice(:,i),slice);
      end
      if which_stats(i,3)
         putimages(out_handle_sdeffect(i),sdeffect_slice(:,i),slice);
      end
   end
   if which_stats(1,4)
      putimages(out_handle_Fstat,Fstat_slice,slice);
   end  
   if which_stats(1,6)
      putimages(out_handle_resid,resid_slice,slice,1:n);
   end  
   if which_stats(1,7)
      putimages(out_handle_wresid,wresid_slice,slice,1:n);
   end  

end

% Tidy up:

closeimage(in_handle);
for i=1:numcontrasts
   if isdelay
      if which_stats(i,1)
         closeimage(out_handle_tdelay(i));
      end
      if which_stats(i,2)
         closeimage(out_handle_delay(i));
      end
      if which_stats(i,3)
         closeimage(out_handle_sddelay(i));
      end
   end
   if which_stats(i,1)
      closeimage(out_handle_tstat(i));
   end
   if which_stats(i,2)
      closeimage(out_handle_effect(i));
   end
   if which_stats(i,3)
      closeimage(out_handle_sdeffect(i));
   end
end
if which_stats(1,4)
   closeimage(out_handle_Fstat);
end
if which_stats(1,5)
   closeimage(out_handle_rho);
end
if which_stats(1,6)
   closeimage(out_handle_resid);
end
if which_stats(1,7)
   closeimage(out_handle_wresid);
end

% End.