function [ t, vectorY, schemeData ] = ...
        odeCFLvector(intFunc, schemeFunc, tspan, vectorY0, options, schemeData)
% odeCFLvector: integrate a set of CFL constrained ODEs in lockstep.
%
% NOTE: This routine is an old attempt to implement vector level sets.  With
% the completion of version 1.1beta, the regular odeCFLn integrators were
% designed to accept vector level sets, so this routine is no longer
% necessary.
%
%   [ t, vectorY, schemeData, stats ] = ...
%       odeCFLvector(intFunc, schemeFunc, tspan, vectorY0, options, schemeData)
%
% Integrates a set of ODEs forward in time by CFL constrained timesteps
%   such that all ODEs take the same timesteps.
%
% Note that this integration routine uses an odeCFLn routine for the
%   underlying timestepping, but does not itself satisfy the prototype
%   of an odeCFLn function.
%
% The primary purpose of this routine is to allow integration of vector
%   level sets.  Most of the code is devoted to keeping the separate
%   elements of the vector in temporal lockstep, and (if necessary)
%   allowing each element to see the others.
%
% parameters:
%   intFunc      Function handle to an odeCFLn integration routine.
%                  This routine wil perform the actual timesteps.
%                  Usually one of odeCFL1, odeCFL2, odeCFL3.
%   schemeFunc	 Function handle to a CFL constrained ODE system.
%                  Typically an approximation to an HJ term.
%                  May be a cell vector of function handles.  See below.
%   tspan        Range of time over which to integrate.  Must contain
%                  only two entries tspan = [ t0 tf ].  Passed to intFunc.
%   vectorY0     A cell vector, each element of which is 
%                  an initial condition vector to pass to intFunc
%                  (typically the data array in vector form).
%   options      An option structure generated by odeCFLset; passed to intFunc
%                  (use [] as a placeholder if necessary).
%   schemeData   Structure passed through to schemeFunc.
%                  May be a cell vector of structures (see below).
%
%
%   t            Output time(s).
%   vectorY      Cell vector, each element of which is 
%                  an output state (see below).
%   schemeData   Output version of schemeData.
%
% The input parameters tspan and options are passed directly to the 
%   integration routine specified by intFunc.  See the help entry for the
%   integration routine to find out how these parameters are interpreted.
%
% Apart from the distinctions mentioned below, schemeFunc and schemeData
%   are also passed to intFunc.  See the help entry for the integration 
%   routine to find out how these parameters are interpreted.
%
% The parameter schemeFunc may be given as
%
%   1) A function handle, in which case it is passed to intFunc
%      for all the elements of vectorY0.  Use this option if every
%      element of the vector level set is updated using the same scheme.
%
%   2) A cell vector the same length as vectorY0.  In that case,
%      corresponding elements of schemeFunc and vectorY0 are passed
%      to intFunc together.  Use this option if every element of the 
%      vector level set is updated using a different scheme.
%
% The parameter schemeData may be given as
%
%   1) A structure, in which case it is passed to intFunc for all
%      elements of vectorY0.  Use this option if every element of the
%      vector level set is updated using (almost) the same scheme data.
%      See below for explanation of "almost".
%
%   2) A cell vector the same length as vectorY0.  In that case,
%      corresponding elements of schemeData and vectorY0 are passed
%      to intFunc together.  Use this option if every element of the
%      vector level set is updated using a different set of scheme data.
%
% The parameter vectorY0 must be given as a cell vector.  Each element 
%   should be a vector which can be interpreted as an initial condition
%   for that element of the vector level set function by the temporal
%   integration routine.
%
% The output parameters are taken from the intFunc integration routine,
%   except that some may be cell vectors to store multiple elements.
%   See the help entry for the integration routine for more information
%   on what these parameters contain.
%
% In order for the members of the vector level set to see one another,
%   odeCFLvector detects the presence of the following field names in
%   the structure(s) passed in by schemeData.  If they are present, 
%   they are updated at each step of the integration in the following way:
%
%   .vectorY     A cell vector of length vectorY0, each element of which
%                  is the state of the corresponding element of the vector
%                  level set as of the end of the last full timestep.
%   .element     A scalar indicating which element of the vector level set
%                  is currently being integrated.
%
% If these fields are not present, they are not added.  In general, 
%   if the user wants to use this feature, then these fields can be set
%   to the empty vector [] so that they are present in the initial
%   schemeData structure(s).
%
% Note that inside schemeFunc the following equality should hold:
%
%               y = schemeData.vectorY{schemeData.element}
%
% Users should note that odeCFLvector may modify the CFLfactor
%   in the options structure in order to achieve the lockstep behavior.
%   Every attempt is made to choose a timestep small enough to maintain
%   the most restrictive CFL timestep, but occasionally the CFL factor
%   may be violated to avoid redoing a timestep.

% Copyright 2004 Ian M. Mitchell (mitchell@cs.ubc.ca).
% This software is used, copied and distributed under the licensing 
%   agreement contained in the file LICENSE in the top directory of 
%   the distribution.
%
% Ian Mitchell, 8/18/04.

  %---------------------------------------------------------------------------
  % How close is "temporal lockstep"? (relative)
  small = 100 * eps;

  %---------------------------------------------------------------------------
  % By how much do we allow the CFL factor to be violated before we
  %   repeat a step?
  % Note that values greater than 1 could conceivably lead to instability if
  %   every timestep is too big, but the occasional overly large timestep
  %   may be compensated for if most timesteps use a conservative factor.
  maximumFactorCFL = 4.0;

  %---------------------------------------------------------------------------
  % How many elements in the vector level set?
  numElements = length(vectorY0);

  %---------------------------------------------------------------------------
  % Check input parameters.
  if(iscell(schemeFunc) && (length(schemeFunc) ~= numElements))
    error('Cell vector schemeFunc must be the same length as vectorY0');
  end

  if(length(tspan) ~= 2)
    error('tspan must contain two elements: tspan = [ t0 tf ]');
  end

  % Make sure we have the default options settings
  if((nargin < 4) | isempty(options))
    options = odeCFLset;
  end

  if(iscell(schemeData) && (length(schemeData) ~= numElements))
    error('Cell vector schemeData must be the same length as vectorY0');
  end

  %---------------------------------------------------------------------------
  % Keeping track of time progression.
  t0 = tspan(1);
  tf = tspan(2);

  %---------------------------------------------------------------------------
  % Modify the options to keep the elements in temporal lockstep.
  %   Also turn off statistics to avoid messy output.
  vectorOptions = options;
  vectorOptions = odeCFLset(vectorOptions, 'SingleStep', 'on', 'Stats', 'off');

  % Use the user supplied CFL factor for estimating timestep size.
  userFactorCFL = odeCFLget(options, 'FactorCFL');

  % User may also have placed a bound on the maximum step size.
  userMaxStep = odeCFLget(options, 'MaxStep');

  % If user wants screen stats, we'll summarize.
  userStats = strcmp(odeCFLget(options, 'Stats'), 'on');

  % If user wants only a single timestep.
  userSingleStep = strcmp(odeCFLget(options, 'SingleStep'), 'on');

  %---------------------------------------------------------------------------
  % Estimate first timestep size by taking a step with the first element.
  startTime = cputime;
  stepCount = 1;
  vectorY = cell(numElements, 1);
  t = zeros(numElements, 1);

  [ thisSchemeData, thisSchemeFunc ] = ...
                                addFields(schemeData, schemeFunc, vectorY0, 1);
  [ t(1), vectorY{1}, thisSchemeData, stats ] = ...
                     feval(intFunc, thisSchemeFunc, [ t0 tf ], vectorY0{1}, ...
                           vectorOptions, thisSchemeData);

  % Store the (possibly modified) schemeData structure.
  if(iscell(schemeData))
    schemeData{1} = thisSchemeData;
  else
    schemeData = thisSchemeData;
  end

  % Record the timestep bound (for estimating the next timestep's size).
  stepBound = zeros(numElements, 1);
  stepBound(1) = stats.stepBounds;

  %---------------------------------------------------------------------------
  % Complete a step of the same size for the other elements.
  vectorStepSize = min(t(1) - t0, userMaxStep);
  vectorOptions = odeCFLset(vectorOptions, 'MaxStep', vectorStepSize, ...
                                           'FactorCFL', maximumFactorCFL);

  for i = 2 : numElements;

    [ thisSchemeData, thisSchemeFunc ] = ...
                                addFields(schemeData, schemeFunc, vectorY0, i);
    [ newT, vectorY{i}, thisSchemeData, stats ] = ...
                     feval(intFunc, thisSchemeFunc, [ t0 tf ], vectorY0{i}, ...
                           vectorOptions, thisSchemeData);

    % THIS SHOULD BE FIXED TO PERFORM ANOTHER TIMESTEP!
    if((tf - newT > small * tf) & ...
       (abs(newT - (t0 + vectorStepSize)) > small * vectorStepSize))
      warning('maximumFactorCFL was too small');
    end

    t(i) = newT;

    % Store the (possibly modified) schemeData structure.
    if(iscell(schemeData))
      schemeData{i} = thisSchemeData;
    else
      schemeData = thisSchemeData;
    end

    % Record the timestep bound (for estimating the next timestep's size).
    stepBound(i) = stats.stepBounds;
  end

  %---------------------------------------------------------------------------
  % Now loop until tf is achieved, using the timestep bound from
  %   the previous timestep to estimate the size of the next timestep.
  %   No looping if the user requested only a single timestep.
  while(~userSingleStep & (tf - min(t) > small * tf))

    stepCount = stepCount + 1;

    % Keep a copy of vectorY in case schemeData needs it.
    vectorYlast = vectorY;

    % Estimate the timestep for this step
    vectorStepsize = min(userFactorCFL * min(stepBound), userMaxStep);

    % Set the options to respect this bound.
    vectorOptions = odeCFLset(vectorOptions, 'MaxStep', vectorStepsize);

    for i = 1 : numElements;
      if(tf - t(i) > small * tf)
        [ thisSchemeData, thisSchemeFunc ] = ...
                             addFields(schemeData, schemeFunc, vectorYlast, i);
        [ newT, vectorY{i}, thisSchemeData, stats ] = ...
                   feval(intFunc, thisSchemeFunc, [ t(i) tf ], vectorY{i}, ...
                         vectorOptions, thisSchemeData);

        % THIS SHOULD BE FIXED TO PERFORM ANOTHER TIMESTEP!
        if((tf - newT > small * tf) & ...
           (abs(newT - (t(i) + vectorStepSize)) > small * vectorStepSize))
          warning('maximumFactorCFL was too small');
        end

        t(i) = newT;

        % Store the (possibly modified) schemeData structure.
        if(iscell(schemeData))
          schemeData{i} = thisSchemeData;
        else
          schemeData = thisSchemeData;
        end

        % Record the timestep bound (for estimating the next timestep's size).
        stepBound(i) = stats.stepBounds;
      end
    end

  end

  %---------------------------------------------------------------------------
  % Finally, check that all elements are at the same final time.
  if(max(t) - min(t) > small * tf)
    warning('Gap of %f between vector elements'' final time', max(t) - min(t));
  else
    t = median(t);
  end

  endTime = cputime;

  if(userStats)
    fprintf('\t%d steps in %g seconds from %g to %g for %d %s\n', ...
            stepCount, endTime - startTime, t0, t, ...
            numElements, 'element vector');
  end



%---------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------------------------------------------------
function [ thisSchemeData, thisSchemeFunc ] = ...
                           addFields(schemeData, schemeFunc, vectorY, element);
% addFields: Add fields vectorY and element to appropriate schemeData structure
%
%   [ thisSchemeData, thisSchemeFunc ] = ...
%                          addFields(schemeData, schemeFunc, vectorY, element);
%
% Finds the appropriate schemeData and schemeFunc elements for member element
%   of the vector level set.  If necessary, updates fields vectorY and
%   element of the schemeData structure.
%
% parameters:
%   schemeData	   The schemeData structure or cell vector passed 
%                    to odeCFLvector.
%   schemeFunc     The schemeFunc function handle or cell vector passed
%                    to odeCFLvector.
%   vectorY	   The state of the vector level set as of the last
%                    completed timestep.
%   element        The index of the element of the vector level set currently
%                    under consideration.
%
%   thisSchemeData The schemeData structure to use for this element of the
%                    vector level set integration.
%   thisSchemeFunc The schemeFunc function handle to use for this element
%                    of the vector level set integration.

  %---------------------------------------------------------------------------
  % Extract appropriate element of schemeData cell vector if necessary.
  if(iscell(schemeData))
    thisSchemeData = schemeData{element};
  else
    thisSchemeData = schemeData;
  end

  %---------------------------------------------------------------------------
  % Extract appropriate element of schemeFunc cell vector if necessary.
  if(iscell(schemeFunc))
    thisSchemeFunc = schemeFunc{element};
  else
    thisSchemeFunc = schemeFunc;
  end

  %---------------------------------------------------------------------------
  % Update fields of schemeData structure if necessary.
  if(isfield(thisSchemeData, 'vectorY'))
    thisSchemeData.vectorY = vectorY;
  end
  if(isfield(thisSchemeData, 'element'))
    thisSchemeData.element = element;
  end

  % We also need to check if this is a schemeData structure for a compound
  %   term (ie termSum, termRestrictUpdate).
  % Ideally this check should be performed recursively.
  if(isfield(thisSchemeData, 'innerData'))
    for i = 1 : length(thisSchemeData.innerData)
      if(isfield(thisSchemeData.innerData{i}, 'vectorY'))
        thisSchemeData.innerData{i}.vectorY = vectorY;
      end
      if(isfield(thisSchemeData.innerData{i}, 'element'))
        thisSchemeData.innerData{i}.element = element;
      end
    end
  end

