qsubcellfun.m

function varargout = qsubcellfun(fname, varargin)

% QSUBCELLFUN applies a function to each element of a cell-array. The
% function execution is done in parallel using the Torque, SGE, PBS or
% SLURM batch queue system.
%
% Use as
%   argout = qsubcellfun(fname, x1, x2, ...)
%
% This function has a number of optional arguments that have to passed
% as key-value pairs at the end of the list of input arguments. All other
% input arguments (including other key-value pairs) will be passed to the
% function to be evaluated.
%   UniformOutput  = boolean (default = false)
%   StopOnError    = boolean (default = true)
%   diary          = string, can be 'always', 'never', 'warning', 'error' (default = 'error')
%   timreq         = number, the time in seconds required to run a single job
%   timoverhead    = number in seconds, how much time to allow MATLAB to start (default = 180 seconds)
%   memreq         = number, the memory in bytes required to run a single job
%   memoverhead    = number in bytes, how much memory to account for MATLAB itself (default depends on the MATLAB version)
%   stack          = number, stack multiple jobs in a single qsub job (default = 'auto')
%   backend        = string, can be 'torque', 'sge', 'slurm', 'lsf', 'system', 'local' (default is automatic)
%   batchid        = string, to identify the jobs in the queue (default is user_host_pid_batch)
%   compile        = string, can be 'auto', 'yes', 'no' (default = 'no')
%   queue          = string, which queue to submit the job in (default is empty)
%   options        = string, additional options that will be passed to qsub/srun (default is empty)
%   matlabcmd      = string, the Linux command line to start MATLAB on the compute nodes (default is automatic
%   display        = 'yes' or 'no', whether the nodisplay option should be passed to MATLAB (default = 'no', meaning nodisplay)
%   jvm            = 'yes' or 'no', whether the nojvm option should be passed to MATLAB (default = 'yes', meaning with jvm)
%   rerunable      = 'yes' or 'no', whether the job can be restarted on a torque/maui/moab cluster (default = 'no')
%   sleep          = number, time in seconds to wait between checks for job completion (default = 0.5 s)
%
% It is required to give an estimate of the time and memory requirements of
% the individual jobs. The memory requirement of the MATLAB executable
% itself will automatically be added, just as the time required to start
% up a new MATLAB process. If you don't know what the memory and time
% requirements of your job are, you can get an estimate for them using
% TIC/TOC and MEMTIC/MEMTOC around a single execution of one of the jobs in
% your interactive MATLAB session. You can also start with very large
% estimates, e.g. 4*1024^3 bytes for the memory (which is 4GB) and 28800
% seconds for the time (which is 8 hours) and then run a single job through
% qsubcellfun. When the job returns, it will print the memory and time it
% required.
%
% Example
%   fname = 'power';
%   x1    = {1, 2, 3, 4, 5};
%   x2    = {2, 2, 2, 2, 2};
%   y     = qsubcellfun(fname, x1, x2, 'memreq', 1024^3, 'timreq', 300);
%
% Using the compile=yes or compile=auto option, you can compile your
% function into a stand-alone executable that can be executed on the cluster
% without requiring additional MATLAB licenses. You can also call the
% QSUBCOMPILE function prior to calling QSUBCELLFUN. If you plan multiple
% batches of the same function, compiling it prior to QSUBCELLFUN is more
% efficient. In that case you will have to delete the compiled executable
% yourself once you are done.
%
% In case you abort your call to qsubcellfun by pressing ctrl-c,
% the already submitted jobs will be canceled. Some small temporary
% files might remain in your working directory.
%
% To check the the status and healthy execution of the jobs on the Torque
% batch queuing system, you can use
%   qstat
%   qstat -an1
%   qstat -Q
% comands on the linux command line. To delete jobs from the Torque batch
% queue and to abort already running jobs, you can use
%   qdel <jobnumber>
%   qdel all
%
% See also QSUBCOMPILE, QSUBFEVAL, CELLFUN, PEERCELLFUN, FEVAL, BATCH

% -----------------------------------------------------------------------
% Copyright (C) 2011-2016, Robert Oostenveld
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/
%
% $Id$
% -----------------------------------------------------------------------

if ft_platform_supports('onCleanup')
  % switch to zombie when finished or when Ctrl-C gets pressed
  % the onCleanup function does not exist for older versions
  onCleanup(@cleanupfun);
end

% remove the persistent lists with job and pbs identifiers
clear qsublist

stopwatch = tic;

% locate the begin of the optional key-value arguments
optbeg = find(cellfun(@ischar, varargin));
optarg = varargin(optbeg:end);

% get the optional input arguments
UniformOutput = ft_getopt(optarg, 'UniformOutput', false   );
StopOnError   = ft_getopt(optarg, 'StopOnError',   true    );
diary         = ft_getopt(optarg, 'diary',         'error' ); % 'always', 'never', 'warning', 'error'
timreq        = ft_getopt(optarg, 'timreq');
memreq        = ft_getopt(optarg, 'memreq');
timoverhead   = ft_getopt(optarg, 'timoverhead', []); % default is set in QSUBFEVAL
memoverhead   = ft_getopt(optarg, 'memoverhead', []); % default is set in QSUBFEVAL
stack         = ft_getopt(optarg, 'stack',   'auto'); % 'auto' or a number
compile       = ft_getopt(optarg, 'compile', 'no');   % can be 'auto', 'yes' or 'no'
backend       = ft_getopt(optarg, 'backend', []);     % the default will be determined by qsubfeval
queue         = ft_getopt(optarg, 'queue',   []);
submitoptions = ft_getopt(optarg, 'options', []);
batch         = ft_getopt(optarg, 'batch',   getbatch());               % this is a number that is automatically incremented
batchid       = ft_getopt(optarg, 'batchid', generatebatchid(batch));   % this is a string like user_host_pid_batch
display       = ft_getopt(optarg, 'display', 'no');
matlabcmd     = ft_getopt(optarg, 'matlabcmd', []);
jvm           = ft_getopt(optarg, 'jvm', 'yes');
whichfunction = ft_getopt(optarg, 'whichfunction');   % the complete filename to the function, including path
rerunable     = ft_getopt(optarg, 'rerunable');       % the default is determined in qsubfeval
sleep         = ft_getopt(optarg, 'sleep', 0.5);      % in seconds

% skip the optional key-value arguments
if ~isempty(optbeg)
  varargin = varargin(1:(optbeg-1));
end

if isstruct(fname)
  % the function has been compiled by qsubcompile
  fcomp = fname;
  % continue with the original function name
  fname = fcomp.fname;
else
  fcomp = [];
end

% determine which function it is
if isempty(whichfunction)
  if ischar(fname)
    whichfunction = which(fname);
  elseif isa(fname, 'function_handle')
    whichfunction = which(func2str(fname));
  end
end

% if the first attempt failed, it might be due a function that is private to the calling function
if isempty(whichfunction)
  s = dbstack('-completenames');
  s = s(2); % qsubcellfun is the first, the calling function is the second
  if ischar(fname)
    whichfunction = which(fullfile(fileparts(s.file), 'private', fname));
  elseif isa(fname, 'function_handle')
    whichfunction = which(fullfile(fileparts(s.file), 'private', func2str(fname)));
  end
  if ~isempty(whichfunction)
    warning('assuming %s as full function name', whichfunction);
  end
  clear s
end

% there are potentially errors to catch from the which() function
if isempty(whichfunction) && ischar(fname)
  error('not a valid M-file "%s"', fname);
end

% determine the number of input arguments and the number of jobs
numargin    = numel(varargin);
numjob      = numel(varargin{1});

% determine the number of MATLAB jobs to "stack" together into separate qsub jobs
if isequal(stack, 'auto')
  if ~isempty(timreq)
    stack = floor(180/timreq);
  else
    stack = 1;
  end
end

% ensure that the stacking is not higher than the number of jobs
stack = min(stack, numjob);

% give some feedback about the stacking
if stack>1
  fprintf('stacking %d MATLAB jobs in each qsub job\n', stack);
end

% prepare some arrays that are used for bookkeeping
jobid       = cell(1, numjob);
puttime     = nan(1, numjob);
timused     = nan(1, numjob);
memused     = nan(1, numjob);
submitted   = false(1, numjob);
collected   = false(1, numjob);
submittime  = inf(1, numjob);
collecttime = inf(1, numjob);

% it can be difficult to determine the number of output arguments
try
  if isequal(fname, 'cellfun') || isequal(fname, @cellfun)
    if isa(varargin{1}{1}, 'char') || isa(varargin{1}{1}, 'function_handle')
      numargout = nargout(varargin{1}{1});
    elseif isa(varargin{1}{1}, 'struct')
      % the function to be executed has been compiled
      fcomp = varargin{1}{1};
      numargout = nargout(fcomp.fname);
    end
  else
    numargout = nargout(fname);
  end
catch
  % the "catch me" syntax is broken on MATLAB74, this fixes it
  nargout_err = lasterror;
  if strcmp(nargout_err.identifier, 'MATLAB:narginout:doesNotApply')
    % e.g. in case of nargin('plus')
    numargout = 1;
  else
    rethrow(nargout_err);
  end
end

if numargout<0
  % the nargout function returns -1 in case of a variable number of output arguments
  numargout = 1;
elseif numargout>nargout
  % the number of output arguments is constrained by the users' call to this function
  numargout = nargout;
elseif nargout>numargout
  error('Too many output arguments.');
end

% running a compiled version in parallel takes no MATLAB licenses
% auto compilation will be attempted if the total batch takes more than 30 minutes
if (strcmp(compile, 'auto') && (numjob*timreq/3600)>0.5) || istrue(compile)
  try
    % try to compile into a stand-allone application
    fcomp = qsubcompile(fname, 'batch', batch, 'batchid', batchid);
  catch
    if istrue(compile)
      % the error that was caught is critical
      rethrow(lasterror);
    elseif strcmp(compile, 'auto')
      % compilation was only optional, the caught error is not critical
      warning(lasterr);
    end
  end % try-catch
end % if compile

if stack>1
  % combine multiple jobs in one, the idea is to use recursion like this
  % a = {{@plus, @plus}, {{1}, {2}}, {{3}, {4}}}
  % b = cellfun(@cellfun, a{:})

  % these options will be passed to the recursive call after being modified further down
  if ~any(strcmpi(optarg, 'timreq'))
    optarg{end+1} = 'timreq';
    optarg{end+1} = timreq;
  end
  if ~any(strcmpi(optarg, 'stack'))
    optarg{end+1} = 'stack';
    optarg{end+1} = stack;
  end
  if ~any(strcmpi(optarg, 'UniformOutput'))
    optarg{end+1} = 'UniformOutput';
    optarg{end+1} = UniformOutput;
  end
  if ~any(strcmpi(optarg, 'whichfunction'))
    optarg{end+1} = 'whichfunction';
    optarg{end+1} = whichfunction;
  end
  if ~any(strcmpi(optarg, 'compile'))
    optarg{end+1} = 'compile';
    optarg{end+1} = compile;
  end

  % update these settings for the recursive call
  optarg{find(strcmpi(optarg, 'timreq'))+1}        = timreq*stack;
  optarg{find(strcmpi(optarg, 'stack'))+1}         = 1;
  optarg{find(strcmpi(optarg, 'UniformOutput'))+1} = false;
  optarg{find(strcmpi(optarg, 'compile'))+1}       = false;

  % FIXME the partitioning can be further perfected
  partition     = floor((0:numjob-1)/stack)+1;
  numpartition  = partition(end);

  stackargin = cell(1,numargin+3); % include the fname, uniformoutput, false
  if istrue(compile)
    if ischar(fcomp.fname)
      % it should contain function handles, not strings
      stackargin{1} = repmat({str2func(fcomp.fname)}, 1, numpartition);
    else
      stackargin{1} = repmat({fcomp.fname}, 1, numpartition);
    end
  else
    if ischar(fname)
      % it should contain function handles, not strings
      stackargin{1} = repmat({str2func(fname)}, 1, numpartition);
    else
      stackargin{1} = repmat({fname}, 1, numpartition);
    end
  end
  stackargin{end-1} = repmat({'uniformoutput'},1,numpartition);  % uniformoutput
  stackargin{end}   = repmat({false},1,numpartition);            % false

  % reorganize the original input into the stacked format
  for i=1:numargin
    tmp = cell(1,numpartition);
    for j=1:numpartition
      tmp{j} = {varargin{i}{partition==j}};
    end
    stackargin{i+1} = tmp; % note that the first element is the fname
    clear tmp
  end

  stackargout = cell(1,numargout);
  [stackargout{:}] = qsubcellfun(@cellfun, stackargin{:}, optarg{:});

  % reorganise the stacked output into the original format
  for i=1:numargout
    tmp = cell(size(varargin{1}));
    for j=1:numpartition
      tmp(partition==j) = stackargout{i}{j};
    end
    varargout{i} = tmp;
    clear tmp
  end

  if numargout>0 && UniformOutput
    [varargout{:}] = makeuniform(varargout{:});
  end

  return;
end

% check the input arguments
for i=1:numargin
  if ~isa(varargin{i}, 'cell')
    error('input argument #%d should be a cell-array', i+1);
  end
  if numel(varargin{i})~=numjob
    error('inconsistent number of elements in input #%d', i+1);
  end
end

for submit=1:numjob
  % redistribute the input arguments
  argin = cell(1, numargin);
  for j=1:numargin
    argin{j} = varargin{j}{submit};
  end

  % submit the job
  if ~isempty(fcomp)
    % use the compiled version
    [curjobid, curputtime] = qsubfeval(fcomp, argin{:}, 'memreq', memreq, 'timreq', timreq, 'memoverhead', memoverhead, 'timoverhead', timoverhead, 'diary', diary, 'batch', batch, 'batchid', batchid, 'backend', backend, 'options', submitoptions, 'queue', queue, 'matlabcmd', matlabcmd, 'display', display, 'jvm', jvm, 'nargout', numargout, 'whichfunction', whichfunction, 'rerunable', rerunable);
  else
    % use the non-compiled version
    [curjobid, curputtime] = qsubfeval(fname, argin{:}, 'memreq', memreq, 'timreq', timreq, 'memoverhead', memoverhead, 'timoverhead', timoverhead, 'diary', diary, 'batch', batch, 'batchid', batchid, 'backend', backend, 'options', submitoptions, 'queue', queue, 'matlabcmd', matlabcmd, 'display', display, 'jvm', jvm, 'nargout', numargout, 'whichfunction', whichfunction, 'rerunable', rerunable);
  end

  % fprintf('submitted job %d\n', submit);
  jobid{submit}      = curjobid;
  puttime(submit)    = curputtime;
  submitted(submit)  = true;
  submittime(submit) = toc(stopwatch);
  clear curjobid curputtime
end % for

while (~all(collected))
  % try to collect the jobs that have finished

  for collect=find(~collected)
    % this will return empty arguments if the job has not finished
    ws = warning('off', 'FieldTrip:qsub:jobNotAvailable');
    [argout, options] = qsubget(jobid{collect}, 'output', 'cell', 'diary', diary, 'StopOnError', StopOnError);
    warning(ws);

    if ~isempty(argout) || ~isempty(options)
      % fprintf('collected job %d\n', collect);
      collected(collect)   = true;
      collecttime(collect) = toc(stopwatch);

      if isempty(argout) && StopOnError==false
        % this happens if an error was detected in qsubget and StopOnError is false
        % replace the output of the failed jobs with []
        argout = repmat({[]}, 1, numargout);
      end

      % redistribute the output arguments
      for j=1:numargout
        varargout{j}{collect} = argout{j};
      end

      % gather the job statistics
      % these are empty in case an error happened during remote evaluation, therefore the default value of NaN is specified
      timused(collect) = ft_getopt(options, 'timused', nan);
      memused(collect) = ft_getopt(options, 'memused', nan);

    end  % if
  end % for

  pausejava(sleep);
end % while

% ensure the output to have the same size/dimensions as the input
for i=1:numargout
  varargout{i} = reshape(varargout{i}, size(varargin{1}));
end

if numargout>0 && UniformOutput
  [varargout{:}] = makeuniform(varargout{:});
end

% clean up the remains of the compilation
if (strcmp(compile, 'yes') || strcmp(compile, 'auto')) && ~isempty(fcomp)
  % the extension might be .app or .exe or none
  system(sprintf('rm -rf %s',         fcomp.batchid)); % on Linux
  system(sprintf('rm -rf %s.app',     fcomp.batchid)); % on Apple OS X
  system(sprintf('rm -rf %s.exe',     fcomp.batchid)); % on Windows
  system(sprintf('rm -rf run_%s*.sh', fcomp.batchid));
end

% compare the time used inside this function with the total execution time
fprintf('computational time = %.1f sec, elapsed = %.1f sec, speedup %.1f x\n', nansum(timused), toc(stopwatch), nansum(timused)/toc(stopwatch));

if all(puttime>timused)
  warning('the job submission took more time than the actual execution');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function varargout = makeuniform(varargin)
varargout = varargin;
numargout = numel(varargin);
% check whether the output can be converted to a uniform one
for i=1:numel(varargout)
  for j=1:numel(varargout{i})
    if numel(varargout{i}{j})~=1
      % this error message is consistent with the one from cellfun
      error('Non-scalar in Uniform output, at index %d, output %d. Set ''UniformOutput'' to false.', j, i);
    end
  end
end
% convert the output to a uniform one
for i=1:numargout
  varargout{i} = [varargout{i}{:}];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = nanmax(x)
y = max(x(~isnan(x(:))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = nanmin(x)
y = min(x(~isnan(x(:))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = nanmean(x)
x = x(~isnan(x(:)));
y = mean(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = nanstd(x)
x = x(~isnan(x(:)));
y = std(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = nansum(x)
x = x(~isnan(x(:)));
y = sum(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cleanupfun
% the qsublist function maintains a persistent list with all jobs
% request it to kill all the jobs and to cleanup all the files
qsublist('killall');