#! /bin/csh -f

#
# glmest
#
# Original Author: Doug Greve
# CVS Revision Info:
#    $Author: nicks $
#    $Date: 2007/01/09 22:41:17 $
#    $Revision: 1.2 $
#
# Copyright (C) 2002-2007,
# The General Hospital Corporation (Boston, MA). 
# All rights reserved.
#
# Distribution, usage and copying of this software is covered under the
# terms found in the License Agreement file named 'COPYING' found in the
# FreeSurfer source code root directory, and duplicated here:
# https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferOpenSourceLicense
#
# General inquiries: freesurfer@nmr.mgh.harvard.edu
# Bug reports: analysis-bugs@nmr.mgh.harvard.edu
#


set inputargs = ($argv);
set VERSION = '$Id: glmest,v 1.2 2007/01/09 22:41:17 nicks Exp $';

set PrintHelp = 0;
set instem = ();
set maskstem = ();
set xmatfile = ();
set acfstem = ();
set acfnmax = ();
set betastem = ();
set resstem = ();
set yhatstem = ();
set synth = 0;
set monly = 0;
set MLF = ();
set LF = ();

## If there are no arguments, just print useage and exit ##
if ( $#argv == 0  )  goto usage_exit;
set n = `echo $argv | grep help | wc -l` 
if($n != 0) then
  set PrintHelp = 1;
  goto usage_exit;
endif

set n = `echo $argv | grep version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif

goto parse_args;
parse_args_return:

goto check_params;
check_params_return:

set StartDate = `date`;
set CurDir = `pwd`;

if($#MLF == 0) set MLF = $betadir/glmest.m
if(-e $MLF) mv $MLF $MLF.bak
echo "Matlab file is $MLF"

if($#LF == 0) set LF = $betastem.log
if(-e $LF) mv $LF $LF.bak
echo "Log file is $LF"

echo "Log file for glmest" >> $LF
echo $StartDate >> $LF
echo $CurDir    >> $LF
which glmest  >> $LF
echo $inputargs >> $LF
df $betadir      >> $LF
echo "matlab file $MLF" >> $LF

#---------------------------------------------#
      #--------------------------------#
tee $MLF <<EOF

tic;
ystem    = '$instem';
maskstem = '$maskstem';
acfstem  = '$acfstem';
xmatfile = '$xmatfile';
acfnmax  = [$acfnmax];
betastem = '$betastem';
resstem  = '$resstem';
yhatstem = '$yhatstem';
betaxmatfile = '$betaxmatfile';

rvarstem = sprintf('%s-var',betastem);

acf = [];
mask = [];
fprintf('\n\n');

%--------- Load design matrix ---------------------------%
XX = load(xmatfile);
if(isempty(XX))
  fprintf('ERROR: could not load %s\n',xmatfile);
  return;
end
if(~isfield(XX,'Xfinal'))
  fprintf('ERROR: %s does not have Xfinal field\n',xmatfile);
  return;quit;
end
if(~isfield(XX,'Nnnc'))
  fprintf('ERROR: %s does not have Nnnc field\n',xmatfile);
  return;quit;
end
if(~isfield(XX,'Navgs_per_cond'))
  fprintf('ERROR: %s does not have Navgs_per_cond field\n',xmatfile);
  return;quit;
end
X = XX.Xfinal;
nbeta = size(X,2);
% Should check condition here %

%--------- Load input --------------------------------%
fprintf('Loading %s (%g)\n',ystem,toc);
y = fmri_ldbvolume(ystem);
if(isempty(y))
  fprintf('ERROR: could not load %s\n',ystem);
  return;quit;
end
[ns nr nc nf] = size(y);
nv = ns*nr*nc;
y = reshape(y,[nv nf])'; %'

if(nf ~= size(X,1))
  fprintf('ERROR: frame mismatch (y=%d,X=%d)\n',nf,size(X,1));
  return; quit;
end
save(betaxmatfile,'ystem','-v4','-append');

% ------------------ Read in acf -----------------------%
if(~isempty(acfstem))
  fprintf('Loading %s (%g)\n',acfstem,toc);
  acf = fmri_ldbvolume(acfstem);
  if(isempty(acf))
    fprintf('ERROR: could not load %s\n',acfstem);
    return; quit;
  end
  sz = size(acf);
  if(prod(sz(1:3)) ~= nv)
    fprintf('ERROR: input/acf dimension mismatch\n');
    return; quit;
  end
  nacf = sz(4);
  acf = reshape(acf,[nv nacf])'; %';

  if(isempty(acfnmax)) acfnmax = nacf; end
  if(acfnmax > nacf) 
    fprintf('ERROR: acfnmax=%d, but cannot exceed %d\n',nmax,nacf);
    return;quit;
  end
  save(betaxmatfile,'acfstem','acfnmax','-v4','-append');
end

% ------------------ Read in mask -----------------------%
if(~isempty(maskstem))
  fprintf('Loading %s (%g)\n',maskstem,toc);
  mask = fmri_ldbvolume(maskstem);
  if(isempty(mask))
    fprintf('ERROR: could not load %s\n',maskstem);
    return;
  end
  imask = find(mask);
  if(isempty(imask))
    fprintf('ERROR: no voxels found in mask\n');
    return;
  end
  nmask = length(imask);
  fprintf('Found %d/%d (%4.1f%%) voxels in mask\n',nmask,nv,100*nmask/nv);
  y   = y(:,imask);
  acf = acf(:,imask);
  save(betaxmatfile,'maskstem','-v4','-append');
  save(betaxmatfile,'imask','-v4','-append');
else
  nmask = nv;
end

DOF = size(X,1)-size(X,2);

if(isempty(acf))
  fprintf('Begining estimation without whitening (%g)\n',toc);
  B = inv(X'*X)*X';
  fprintf('Computing beta (%g)\n',toc);
  beta = B*y;
  fprintf('Signal Estimate (%g)\n',toc);
  yhat = X*beta;
  fprintf('Residual  (%g)\n',toc);
  res = y-yhat;
  fprintf('Residual Variance (%g)\n',toc);
  rvar = sum(res.^2,1)/DOF;  
else
  fprintf('Begining estimation with per-vox whitening (%g)\n',toc);
  beta   = zeros(nbeta,nmask);
  res    = zeros(nf,nmask);
  resvar = zeros(1,nmask);
  fprintf('nvoxels = %d\n',nmask);
  for v = 1:nmask
    if(mod(v,1000)==0) fprintf('%d %g\n',v,toc); end
    acfv = acf(:,v);
    yv = fast_conv_invirf(y(:,v),acfv);
    Xv = fast_conv_invirf(X,acfv);
    betav = (inv(Xv'*Xv)*Xv')*yv; 
    Rv = eye(nf) - Xv*inv(Xv'*Xv)*Xv';    
    resv = Rv*yv;
    rvarv = sum(resv.^2)/DOF;
    res(:,v)    = resv;
    beta(:,v) = betav;
    rvar(v)   = rvarv;
  end % end loop over voxel
end
fprintf('Estimation finished (%g)\n',toc);

if(~isempty(yhatstem))
  fprintf('Computing PMF of Signal Estimate (%g)\n',toc);
  % Only compute the estimate from the task components %
  Xpmf = X;
  n0 = XX.Nnnc*XX.Navgs_per_cond + 1;
  Xpmf(:,n0:nbeta) = 0;
  yhat = Xpmf*beta;   
end

%------ Save Beta and Variance of the Residual ----------%
if(~isempty(mask))
  fprintf('Unmasking Beta (%g)\n',toc);
  betasave = zeros(nbeta,nv);
  betasave(:,imask) = beta;
  rvarsave = zeros(1,nv);
  rvarsave(imask) = rvar;
else
  betasave = beta;
  rvarsave = rvar;
end
clear beta rvar;
fprintf('Saving Beta (%g)\n',toc);
betasave = reshape(betasave',[ns nr nc nbeta]); %'
fmri_svbvolume(betasave,betastem);
clear betasave;
rvarsave = reshape(rvarsave',[ns nr nc 1]); %'
fmri_svbvolume(rvarsave,rvarstem);

%------ Save Residual ----------%
if(~isempty(resstem))
  if(~isempty(mask))
    fprintf('Unmasking Residual (%g)\n',toc);
    ressave = zeros(nf,nv);
    ressave(:,imask) = res;
  else
    ressave = res;
  end
  clear res;
  fprintf('Saving Residual (%g)\n',toc);
  ressave = reshape(ressave',[ns nr nc nf]); %'
  fmri_svbvolume(ressave,resstem);
  clear ressave;
end

%------ Save Signal Estimate ----------%
if(~isempty(yhatstem))
  if(~isempty(mask))
    fprintf('Unmasking Signal Estimate (%g)\n',toc);
    yhatsave = zeros(nf,nv);
    yhatsave(:,imask) = yhat;
  else
    yhatsave = yhat;
  end
  clear yhat;
  fprintf('Saving Residual (%g)\n',toc);
  yhatsave = reshape(yhatsave',[ns nr nc nf]); %'
  fmri_svbvolume(yhatsave,yhatstem);
  clear yhatsave;
end


return;
quit;


EOF
      #--------------------------------#
#---------------------------------------------#

cp $xmatfile $betaxmatfile;
if(! $monly) then
  cat $MLF | matlab -display iconic | tee -a $LF
  rm $MLF
  if(-e $instem.bhdr) then
    cp $instem.bhdr $betastem.bhdr
    cp $instem.bhdr $betastem-var.bhdr
    if($#resstem)  cp $instem.bhdr $resstem.bhdr
    if($#yhatstem) cp $instem.bhdr $yhatstem.bhdr
  endif
  echo "" | tee -a $LF
  echo "" | tee -a $LF
endif


set EndDate = `date`;
echo "Started at $StartDate" | tee -a $LF
echo "Ended   at $EndDate"   | tee -a $LF
echo " " | tee -a $LF
echo " " | tee -a $LF
echo "glmest: finished" | tee -a $LF
echo " " | tee -a $LF
 #--------------------------------#

exit 0;
############################################################

############--------------##################
parse_args:

set cmdline = ($argv);
while( $#argv != 0 )

  set flag = $argv[1]; shift;
  
  switch($flag)

    case "-i":
      if ( $#argv == 0) goto arg1err;
      set instem = $argv[1]; shift;
      set instem = `getfullpath $instem`;
      if($status) exit 1;
      breaksw

    case "-mask":
      if ( $#argv == 0) goto arg1err;
      set maskstem = $argv[1]; shift;
      set maskstem = `getfullpath $maskstem`;
      if($status) exit 1;
      breaksw

    case "-x":
      if ( $#argv == 0) goto arg1err;
      set xmatfile = $argv[1]; shift;
      set xmatfile = `getfullpath $xmatfile`;
      if($status) exit 1;
      breaksw

    case "-acf":
      if ( $#argv == 0) goto arg1err;
      set acfstem = $argv[1]; shift;
      set acfstem = `getfullpath $acfstem`;
      if($status) exit 1;
      breaksw

    case "-beta":
      if ( $#argv == 0) goto arg1err;
      set betastem = $argv[1]; shift;
      breaksw

    case "-res":
      if ( $#argv == 0) goto arg1err;
      set resstem = $argv[1]; shift;
      breaksw

    case "-yhat":
      if ( $#argv == 0) goto arg1err;
      set yhatstem = $argv[1]; shift;
      breaksw

    case "-lf":
      if ( $#argv == 0) goto arg1err;
      set LF = $argv[1]; shift;
      breaksw

    case "-acfnmax":
      if ( $#argv == 0) goto arg1err;
      set acfnmax = $argv[1]; shift;
      breaksw

    case "-monly":
      if ( $#argv == 0) goto arg1err;
      set MLF = $argv[1]; shift;
      set monly = 1; 
      breaksw

    case "-synth":
      if ( $#argv == 0) goto arg1err;
      set seed  = $argv[1]; shift;
      set synth = 1;
      breaksw

    case "-umask":
      if ( $#argv == 0) goto arg1err;
      umask $1; shift;
      breaksw

    case "-debug":
      set verbose = 1;
      set echo = 1;
      breaksw

    default:
      echo ERROR: Flag $flag unrecognized.
      echo $cmdline
      exit 1
      breaksw
  endsw

end
goto parse_args_return;
############--------------##################

############--------------##################
check_params:

  if($#instem == 0) then
    echo "ERROR: must specify an input";
    exit 1
  endif

  if($#xmatfile == 0) then
    echo "ERROR: must specify an xmatfile";
    exit 1
  endif

  set betadir = `dirname $betastem`;
  mkdir -p $betadir
  if($status) then
    echo "ERROR: could not create $betadir"
    exit 1;
  endif
  set betaxmatfile = $betastem.X.mat;

  if($#resstem) then
    set resdir = `dirname $resstem`;
    mkdir -p $resdir
    if($status) then
      echo "ERROR: could not create $resdir"
      exit 1;
    endif
  endif

  if($#yhatstem) then
    set yhatdir = `dirname $yhatstem`;
    mkdir -p $yhatdir
    if($status) then
      echo "ERROR: could not create $yhatdir"
      exit 1;
    endif
  endif

  if($#LF != 0) then
    set d = `dirname $LF`;
    mkdir -p $d
    if($status) then
      echo "ERROR: could not create $d"
      exit 1;
    endif
  endif

  if($#MLF != 0) then
    set d = `dirname $MLF`;
    mkdir -p $d
    if($status) then
      echo "ERROR: could not create $d"
      exit 1;
    endif
  endif

goto check_params_return;
############--------------##################


############--------------##################
arg1err:
  echo "ERROR: flag $flag requires one argument"
  exit 1
############--------------##################

#--------------------------------------------------------------------#
usage_exit:
  echo "USAGE: glmest"
  echo ""
  echo " -i instem"
  echo " -x matfile     : matlab file with design matrix "
  echo " -acf acfstem   : per-vox pre-filtering by inverting ACF"
  echo " -acfnmax n     : truncate ACF at nmax lags"
  echo " -mask maskstem : only perform computations on masked voxels"
  echo ""
  echo " -beta betastem : GLM parameters (var saved in betastem-var)"
  echo " -res  resstem  : save residuals"
  echo " -yhat yhatstem : save signal estimate"
  echo ""
  echo " -umask umask          : set unix file permission mask"
  echo " -version              : print version and exit"
  echo " -lf logfile           : default is betastem.log "
  echo " -help "
  echo ""


  if($PrintHelp) \
  cat $0 | awk 'BEGIN{prt=0}{if(prt) print $0; if($1 == "BEGINHELP") prt = 1 }'

exit 1;
#--------------------------------------------------------------------#

#---- Everything below here is printed out as part of help -----#
BEGINHELP