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

%FMRILM_ARP
%
% Reads in an fMRI minc file and writes out a T statistic minc file. 
% The method is based on linear models with correlated AR(p) errors:
% Y=hrf*Xb+e, e_t=a_1 e_(t-1) + ... + a_p e_(t-p) + white noise_t. 
% 
% [DF,P] = FMRILM_ARP(INPUT_FILE, OUTPUT_FILE_BASE, X_CACHE, CONTRAST [,
% EXCLUDE [,WHICH_STATS [,FWHM_RHO [,N_POLY [,REFERENCE_DELAY 
% [, NUMLAGS]]]]]] )
% 
% 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 non-excluded  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 parameters, OUTPUT_FILE_BASE_rho.mnc. The first 
%      NUMLAGS frames are the autocorrelations at lags 1:NUMLAGS,
%      the next NUMLAGS frames are the AR parameters a_1 ... a_NUMLAGS.
%      All parameters are bias corrected and smoothed (see below).
%   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.
%  
% NUMLAGS is the order (p) of the autoregressive model. Default is 1.
%
% 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
if nargin < 10
   numlags=1
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

% Check for estimability:

tolerance=0.0000001
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];
      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
   
% 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

Coradj_vol=zeros(numpix*numlags,numslices);
if fwhm_rho < Inf 
   
   if which_stats(1,5)
      out_handle_rho=newimage( ...
         [deblank(output_file_base(1,:)) '_rho.mnc'], ...
         [2*numlags 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:
       
   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;
         M=zeros(numlags+1);
         for i=1:(numlags+1)
            for j=1:(numlags+1)
               Di=(diag(ones(1,n-i+1),i-1)+diag(ones(1,n-i+1),-i+1))/(1+(i==1));
               Dj=(diag(ones(1,n-j+1),j-1)+diag(ones(1,n-j+1),-j+1))/(1+(j==1));
               M(i,j)=trace(R*Di*R*Dj)/(1+(i>1));
            end
         end
         invM=inv(M);
      end
      
      % Loop over pixels in slice to find rho: 
      
      Coradj_slice=zeros(numpix,numlags);
      for y=1:numys 
         img = getimages (in_handle, slice, 1:numframes, [], y);
         Y=img(:,keep)';
         betahat_ls=pinvX*Y;
         resid=Y-X*betahat_ls;
         for lag=0:numlags
            Cov(lag+1,:)=sum(resid(1:(n-lag),:).*resid((lag+1):n,:));
         end
         Covadj=invM*Cov;
         ind=find(Covadj(1,:)>0);
         Coradj_slice(ind+(y-1)*numxs,:)=  ...
            (Covadj(2:(numlags+1),ind)./(ones(numlags,1)*Covadj(1,ind)))';
      end
      for lag=1:numlags
         Cor_xy=reshape(Coradj_slice(:,lag),numxs,numys);
         ind=(1:numpix)+(lag-1)*numpix;
         Coradj_vol(ind,slice)=reshape(conv2(Cor_xy,ker_xy,'same'),numpix,1);   
      end
      
   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);
      Coradj_vol=Coradj_vol*K;
   end
   
   if which_stats(1,5)
      for lag=1:numlags
         ind=(1:numpix)+(lag-1)*numpix;
         putimages(out_handle_rho,Coradj_vol(ind,:),1:numslices,lag);
      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

sump=0;

% 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
      p=rank(X);
      df=n-p;
   end
   sump=sump+p;
   A_slice=zeros(numpix,numlags-1);
  
   for y=1:numys 
      img = getimages (in_handle, slice, 1:numframes, [], y);
      Y=img(:,keep)';
      
      for x=1:numxs
         pix=x+(y-1)*numxs;
         Coradj_pix=Coradj_vol(pix+((1:numlags)-1)*numpix,slice);
         [Ainvt posdef]=chol(toeplitz([1 Coradj_pix']));
         nl=size(Ainvt,1);
         A=inv(Ainvt');
         A_slice(pix,1:(nl-1))=-A(nl,(nl-1):-1:1)/A(nl,nl);
         B=ones(n-nl,1)*A(nl,:);
         Vmhalf=spdiags(B,1:nl,n-nl,n);
         
         Ystar=zeros(n,1);
         Ystar(1:nl)=A*Y(1:nl,x);
         Ystar((nl+1):n)=Vmhalf*Y(:,x);
         
         Xstar(1:nl,:)=A*X(1:nl,:);
         Xstar((nl+1):n,:)=Vmhalf*X;
         pinvXstar=pinv(Xstar);
         if isdelay
            X2star(1:nl,:)=A*X2(1:nl,:);
            X2star((nl+1):n,:)=Vmhalf*X2;
            D=sum(pinvXstar(colsX1,:).*X2star',2);
            pinvXstar01=pinvXstar(colsX01,:);
            V=pinvXstar01*pinvXstar01';
         end
         CpinvX=Contrast*pinvXstar;
         vareffect=CpinvX*CpinvX';
         sdeffect_mult=sqrt(diag(vareffect)/df);
         if which_stats(1,4)
            vareffinv=pinv(vareffect)*df;
         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,5)
      putimages(out_handle_rho,A_slice,slice,(1:numlags)+numlags);
   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

p=round(sump/numslices)
df=n-p

% 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.