bdEval.m

%bdEval  Evaluate the solution of a differential equation.
%   [y,yp] = bdEval(sol,xint,idx) evaluates the solution of the
%   differential equation (sol) at the time points (xint). It is
%   equivalent to the matlab DEVAL function except that it also
%   works for SOL structures generated by third-party solvers.
%   It returns both the values of the solution (y) and its first
%   derivative (yp). The optional idx parameter specifies which
%   rows of the solution (sol.y) to return.
%
%   [y,yp] = bdEval(sol,xint,sys,varname) is an alternative syntax
%   that extracts a system variable by its name and reshapes it
%   accordingly.
%
%EXAMPLE 1
%   sys = LinearODE;                         % Linear system of ODEs
%   tspan = [0 10];                          % Integration time domain
%   sol = bdSolve(sys,tspan,@ode45);         % Apply the ode45 solver
%   tplot = 0:0.1:10;                        % Interpolation time points
%   [y,yp] = bdEval(sol,tplot);              % Interpolate the solution
%   figure; plot(tplot,y); legend('y1','y2');
%   figure; plot(tplot,yp); legend('dy1/dt','dy2/dt');
%
%EXAMPLE 2
%   sys = LinearODE;                         % Linear system of ODEs
%   tspan = [0 10];                          % Integration time domain
%   sol = bdSolve(sys,tspan,@ode45);         % Apply the ode45 solver
%   tplot = 0:0.1:10;                        % Interpolation time points
%   [x,xp] = bdEval(sol,tplot,sys,'x');      % Interpolate the solution in 'x'
%   [y,yp] = bdEval(sol,tplot,sys,'y');      % Interpolate the solution in 'y'
%   figure; plot(tplot,x, tplot,y); legend('x','y')
%   figure; plot(tplot,xp, tplot,yp); legend('dx/dt','dy/dt');
% 
%AUTHORS
%  Stewart Heitmann (2016a,2019a)

% Copyright (C) 2016-2019 QIMR Berghofer Medical Research Institute
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions
% are met:
%
% 1. Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
% 
% 2. Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in
%    the documentation and/or other materials provided with the
%    distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
% FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
% COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
% LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
% ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
function [y,yp] = bdEval(sol,xint,varargin)
    if ~isstruct(sol) || ~isfield(sol,'x') || ~isfield(sol,'y')
        throwAsCaller(MException('bdEval:InvalidSol','bdEval(sol,xint): sol must be a structure returned by a solver'));
    end
    
    % number of variables in sol
    n = size(sol.y,1);
    
    % post-processing flag
    reshapeflag = false;
    
    switch nargin
        case 1
            throwAsCaller(MException('bdEval:Syntax','bdEval(sol): Not enough input arguments'));
            
        case 2
            % bdEval(sol,xint) where indx is an optional argument
            idx = 1:n;           

        case 3
            % bdEval(sol,xint,idx)
            idx = varargin{1};
            if ~isnumeric(idx)
                throwAsCaller(MException('bdEval:InvalidIdx','bdEval(sol,xint,idx): idx must be numeric'));
            end
            % verify the indexes are within bounds
            if any(idx<1) || any(idx>n)
                throwAsCaller(MException('bdEval:InvalidIdx','bdEval(sol,xint,idx): idx is out of bounds'));
            end 
            
        case 4
            % bdEval(sol,xint,sys,varname)
            sys = varargin{1};
            varname = varargin{2};
            if ~isstruct(sys) || ~isfield(sys,'vardef')
                throwAsCaller(MException('bdEval:InvalidSys','bdEval(sol,xint,sys,varname): sys must be a bdtoolbox system structure'));
            end
            [~,varindx,idx] = bdGetVar(sys,varname);
            if isempty(idx)
                throwAsCaller(MException('bdEval:InvalidName','bdEval(sol,xint,sys,varname): ''%s'' is not a valid variable name for sys',varname));
            end
            varsize = size(sys.vardef(varindx).value);
            reshapeflag = true;
            
        otherwise
            throwAsCaller(MException('bdEval:Syntax','bdEval: Too many input arguments'));
    end
    
    % verify the xint values are within bounds
    if any(xint<sol.x(1)) || any(xint>sol.x(end))
        error('bdEval(sol,xint): xint is out of bounds');
    end
    
    switch sol.solver
        case {'ode45','ode23','ode113','ode15s','ide23s','ode23t','ode23tb','dde23'}
            % Use MATLAB deval for MATLAB solvers
            [y,yp] = deval(sol,xint,idx);       
            
        otherwise
            % Use interpolation for third-party solvers.
            % Annoyingly, interp1() transposes the output when the input is
            % a matrix but not a vector. So we treat both cases separately.
            if size(idx,2)==1
                % Here the input is a VECTOR so we DON'T transpose the output.
                y  = interp1(sol.x, sol.y(idx,:)', xint); 
            else
                % Here the input is a MATRIX so we DO transpose the output.
                y  = interp1(sol.x, sol.y(idx,:)', xint)'; 
            end
            
            % Return the gradient vector (if requested)
            if nargout>1
                if isfield(sol,'yp')
                    % The solver has already computed the gradient, so we
                    % only need to interpolate those values.
                    if size(idx,2)==1
                        % Here the input is a VECTOR so we DON'T transpose the output.
                        yp = interp1(sol.x, sol.yp(idx,:)', xint);
                    else
                        % Here the input is a MATRIX so we DO transpose the output.
                        yp = interp1(sol.x, sol.yp(idx,:)', xint)'; 
                    end
                else
                    % Compute the gradient from scratch
                    dt = sol.x(2)-sol.x(1);
                    yp = gradient(sol.y(idx,:),dt);
                end
            end
    end
    
    % post-processing 
    if reshapeflag
        % reshape the output to match that in sys.vardef
        nr = varsize(1);
        nc = varsize(2);
        nt = numel(xint);
        
        % compute the new shape of y
        if nr==1
            % var has one row
            if nc==1
                % var has one row and one column
                % Return y as (1 x nt).
                newshape = [1 nt];
            else
                % var has one row and multiple columns
                % Return y as (nc x nt)
                newshape = [nc nt];
            end 
        else
            % var has multiple rows
            if nc==1
                % var has multiple rows and one column
                % Return y as (nr x nt)
                newshape = [nr nt];
            else
                % var has multiple rows and multiple columns
                % Return y as (nr x nc x nt)
                newshape = [nr nc nt];
            end 
        end
        
        % apply the new shape
        y = reshape(y,newshape);
        if nargout>1
            yp = reshape(yp,newshape);
        end
    end
end