image thumbnail

evalit.m v3 (Jul 2013)

by

 

20 Jun 2013 (Updated )

Evaluates a mathematical expression and rounds it accordingly to the error propagation theory.

evalit.m
function [VAL,ERR,PRE,SIG,ORD,ECS,PDS] = evalit(FUNC,VARS,VALS,ERRS,varargin)
%EVALIT   Evaluates and rounds a function accordingly to error propagation.
%
%SYNTAX:
%                                     evalit(FUNC,VARS,VALS,ERRS)
%                               VAL = evalit(FUNC,VARS,VALS,ERRS)
%                         [VAL,ERR] = evalit(FUNC,VARS,VALS,ERRS)
%                     [VAL,ERR,PRE] = evalit(FUNC,VARS,VALS,ERRS)
%                 [VAL,ERR,PRE,SIG] = evalit(FUNC,VARS,VALS,ERRS)
%             [VAL,ERR,PRE,SIG,ORD] = evalit(FUNC,VARS,VALS,ERRS)
%         [VAL,ERR,PRE,SIG,ORD,ECS] = evalit(FUNC,VARS,VALS,ERRS)
%     [VAL,ERR,PRE,SIG,ORD,ECS,PDS] = evalit(FUNC,VARS,VALS,ERRS)
%                                 ... evalit(FUNC,VARS,VALS,ERRS,...,true)
%                                 ... evalit(FUNC,VARS,VALS,ERRS,...,1)
%                                 ... evalit(FUNC,VARS,VALS,ERRS,...,'sum')
%                                 ... evalit(FUNC,VARS,VALS,ERRS,...,PDS)
%
%INPUTS:
%   FUNC  - Function to be evaluated. It can be a symbolic expression or an
%                                     anonymous function_handle.
%   VARS  - Arguments' names. A symbolic or strings array. Like {'x','y'}.
%   VALS  - Arguments' values. A numerical matrix with arguments' values
%                              arranged by columns, that is, its number of
%                              columns must be equal to the number of
%                              arguments. Like [x(:) y(:)].
%   ERRS  - Arguments' errors. A numerical matrix with corresponding errors
%                              or significant figures (negative integers),
%                              of the same size as VALS.
%   true  - Specifies if the result should be printed or not.
%           DEFAULT: true (when no output arguments)
%   1     - ERRS' number of significant figures or digits.
%           DEFAULT: 2.
%   'sum' - Calculates errors by simple SUMmation of error contributions.
%           DEFAULT: 'rss' (via square Root of the Sum of Square errors).
%   PDS   - FUNC's partial derivatives. Array of symbolic or function
%                                       handles, with the same length as
%                                       VARS.
%           DEFAULT: [] (calculated via DIFF from the Symbolic Toolbox).
%
%OUTPUTS:
%   VAL - Values of the evaluated expression rounded to the same precision
%         as its error.
%   ERR - VAL's error.
%   PRE - VAL's precision.
%   SIG - VAL's significant figures.
%   ORD - VAL's orders of magnitude.
%   ECS - Rounded errors contributions to ERR from each argument VARS,
%         ordered by columns, that is, its number of columns is equal to
%         the number of arguments.
%   PDS - Corresponding symbolic partial derivatives. Cell array of
%         symbolic expressions or anonymous function_handles.
%
%DESCRIPTION:
%   This program evaluates a function and rounds it accordingly to the
%   error propagation formula (if 'rss'):
%      dF = ( (dX*(dF/dX))^2 + (dY*(dF/dY))^2 + ... )^(1/2),
%   where
%              X,Y,... - are the (independent or uncorrelated) arguments,
%            dX,dY,... - their errors,
%                    F - the function,
%                   dF - its error, and
%      dF/dX,dF/dY,... - its evaluated partial derivatives.
%
%   Or via
%       dF = dX*abs(dF/dX) + dY*abs(dF/dY) + ...
%   when option 'sum' is used.
%
%EXAMPLE 1 (Use of significant figures with ideal gas equation state):
%   syms p R T      % Variables
%   d = p/(R*T);    % Density from the ideal gas state equation, kg/m3
%   vp = 8.40e4;    % Pressure with 3 significant figures, Pa.
%   vR = 287.0;     % Dry air's constant with 4 significant figures, J/Kkg.
%   vT = 292.13;    % Absolute air temperature with an error of +/-0.32 K.
%   evalit(d,[p R T],[vp vR vT],[-3 -4 0.32])
%
%   DISPLAYS:
%
%      EVALIT: FUNC(p,R,T) = p/(R*T)
%      --------------------------------------------------
%       value +/- error (arguments' error contributions)
%      --------------------------------------------------
%        1.0019 +/- 0.0013 ( 0.00060 + 0.00017 + 0.0011 )
%      --------------------------------------------------
%
%   CONCLUSION: In this example air's density is 1.0019 kg/m3, with an
%               error of +/-0.0013 kg/m3. The biggest contribution to this
%               error is from the temperature T (0.0011 kg/m3), and the
%               smallest one from the constant R (0.00017 kg/m3). Without
%               using this function, you'll get
%                  >> fprintf('%30.29f\n',vp/(vR*vT))
%                  1.00189274237246530000000000000...,
%               that is, like having a super precision of 16 decimal places
%               (DOUBLE class floating point precision). When actually you
%               have a maximum certainty of only 3 decimal places (with the
%               error rounded to 1 significant figure):
%                  (1.002 +/- 0.001) kg/m3.
%
%EXAMPLE 2 (Matrix input, with error rounded to 1 significant figure):
%   syms p R T
%   d  = p/(R*T);
%   vp = 8.40e4; ep =  0;     % Double-float precision error for pressure.
%   vR =  287.0; eR = -4;     % 4 significant figures for constant.
%   vT = [273.15:5:298.15];   % Temperature array.
%   eT = 0.01;                % All temperature's error. Can be an array.
%   evalit(d,[p R T],{vp vR vT},{ep eR eT},1)
%
%   DISPLAYS:
%
%      EVALIT: FUNC(p,R,T) = p/(R*T)
%      --------------------------------------------------
%       value +/- error (arguments' error contributions)
%      --------------------------------------------------
%        1.0715 +/- 0.0002 ( 0.0000000000000001  + 0.0002 + 0.00004 )
%        1.0522 +/- 0.0002 ( 0.0000000000000001  + 0.0002 + 0.00004 )
%        1.0337 +/- 0.0002 ( 0.0000000000000001  + 0.0002 + 0.00004 )
%        1.0157 +/- 0.0002 ( 0.0000000000000001  + 0.0002 + 0.00004 )
%        0.9984 +/- 0.0002 ( 0.0000000000000001  + 0.0002 + 0.00003 )
%        0.9817 +/- 0.0002 ( 0.00000000000000009 + 0.0002 + 0.00003 )
%      --------------------------------------------------
%
%EXAMPLE 3 (Same as before but without using the Symbolic Math Toolbox)
%   d  = @(p,R,T)  p./(R   .*T);
%   dp = @(R,T)    1./(R   .*T);
%   dR = @(p,R,T) -p./(R.^2.*T);
%   dT = @(p,R,T) -p./(R   .*T.^2);
%   vp = 8.40e4;             ep =    0;
%   vR =  287.0;             eR =   -4;
%   vT = [273.15:5:298.15];  eT = 0.01;
%   evalit(d,{'p','R','T'},{vp vR vT},{ep eR eT},{dp dR dT},1)
%
%SEE ALSO:
%   SYM/SYMS, SYM/VPA, SYM/DIFF, SYM/SUBS, ROUND, DOUBLE.
%
%   -----
%   MFILE:   evalit.m
%   VERSION: 3.0 (JUL 21, 2013) (<a href="matlab:web('http://www.mathworks.com/matlabcentral/fileexchange/authors/11258')">download</a>)
%   MATLAB:  7.7.0.471 (R2008b)
%   AUTHOR:  Carlos Adrian Vargas Aguilera (MEXICO)
%   CONTACT: nubeobscura@hotmail.com

%   REVISIONS:
%   1.0      Released. (JUN 20, 2013)
%   2.0      Now handles DOUBLE's precision error (EVAL). Changed 'max' and
%            'rms' inputs for 'sum' and 'rss' respectively. Forces DOUBLE
%            class for numerical inputs. Other minor changes. (JUL 03,
%            2013)
%   3.0      Now handles COMPLEX VALS. Fixed bug when VARS dos not coincide
%            with FUNC's arguments. VARS and PDS input arrays can be cells
%            or matrix. ERRS input can be empty, a single value, a single
%            row or a cell array. Empty, zero or not finite errors are set
%            to double floating point precision error. New display format
%            with ordered dot. Uses  higher order approx. to compare and to
%            change zero errors. (JUL 21, 2013)

%   DISCLAIMER:
%   evalit.m is provided "as is" without warranty of any kind, under the
%   revised BSD license.
%
%   Copyright (c) 2013 Carlos Adrian Vargas Aguilera.
%   All rights reserved.

%   NOTES:
%      * This program have so much lines to make it optionally work
%        without the Symbolic Math Toolbox and to find out any input
%        mistake.
%      * Input arrays can be given as cell o matrix. When they are cells
%        they should have as many elements as the number of elements in
%        VARS. When matrix, then its number of columns should be equal to
%        the number of elements in VARS.
%      * Except for the error, which can be even an empty input. When this
%        is the case, then the floating point error for every argument is
%        considered.
%      * Same happens when the error input is zero.
%      * Arguments names provided in VARS do not need to have the same
%        order as in the function FUNC, which generally is abc ordered.

%%%%%%%%%%%%%%%%%%%%%
% PARAMETERS %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

% DEFAULTS
OUT = (nargout==0); % true or false   : Prints result?
ESF = 2;            % 1 or 2          : ERRS significant figures.
MET = 'rss';        % 'rss'  or 'sum' : Root of squares or addition.
PDS = {};           % Empty cell      : Partial derivatives.

% CONSTANTS
% Note: When the error of a rounded value is unknown, it may be equal to
%   one unit of value's precision (position of its last digit), that is,
%   one half below and one half above:
%                           rounded value: 3.1416
%      value range before rounding (unit): 3.14155 < 3.1416 < 3.14165
%                               so, error: 3.1416  0.00005
%   That is so because any number outside that range is rounded to another
%   number. For example, 3.14153 is rounded to 3.1415.
%
%   Anyway, if you thing the error is twice bigger, then set MIN as 1, and
%   remember this is used only when instead of the error the significant
%   figures of the value are provided by a negative value on ERRS.
MIN         = 0.5;  % 0.5 or 1.0.
VPA         = 18;   % >15. Precision for the VPA.
checkSecond = true; % Checks if the first order error aprox. is smallest.
higherOrder = 5;    % Checks zero error up to this order.

%%%%%%%%%%%%%%%%%%%%%
% CHECKS INPUTS %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

% Checks number of in/out-puts
assert(nargin>=4,'EVALIT:incorrectInput',...
    'Too few input arguments. At least 4.')
assert(nargout<=7,'EVALIT:incorrectOutput',...
    'Too many output arguments. At most 7.')

% Reads optional inputs
[OUT,ESF,MET,PDS] = readOptions(OUT,ESF,MET,PDS,varargin{:});
clear varargin

% Checks FUNC input
FUNC = readFUNC(FUNC);

% Checks VARS input
VARS    = readVARS(VARS);
VARSnum = length(VARS);

% Checks arguments names
[VARSstr,FUNCind,FUNCstr] = readVARSnames(FUNC,VARS,length(PDS));

% Checks VALS input
if isempty(VALS) % No values?
    if OUT
        disp(' ')
        fprintf(1,'EVALIT: FUNC(%s) = %s\n',...
            regexprep(sprintf(',%s',VARSstr{:}),'^.',''),FUNCstr)
        disp('--------------------------------------------------')
        disp(' value +/- error (arguments'' error contributions) ')
        disp('--------------------------------------------------')
        disp('--------------------------------------------------')
        disp(' ')
    end
    VAL=[]; ERR=[]; PRE=[]; SIG=[]; ORD=[]; ECS=[];
    switch nargout
        case 0, clear VAL ERR PRE SIG ORD ECS PDS
        case 1, clear ERR PRE SIG ORD ECS PDS
        case 2, clear PRE SIG ORD ECS PDS
        case 3, clear SIG ORD ECS PDS
        case 4, clear ORD ECS PDS
        case 5, clear ECS PDS
        case 6, clear PDS
    end
    return
end
VALS    = readVALS(VALS,VARSnum);
VALSnum = size(VALS,1);

% Checks ERRS input
[ERRS,ERRSneg] = readERRS(ERRS,VARSnum,VALSnum);

%%%%%%%%%%%%%%%%%%%%%
% MAIN %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

% Gets absolute data by the maximum value of its real or imaginary part,
% instead of from its magnitude (VALSmax = max(abs(VALS)))
VALSmax = max(abs(real(VALS)),abs(imag(VALS)));

% Changes significant figures on ERRS (negative values) to the 'real' error
% in VALS' last digit
VALSord       = floor(log10(VALSmax));
if (VALSnum>1)
    % Checks for -Inf values
    [rflag,cflag] = find(VALSmax==0);
    for col = unique(cflag(:)')
        ind = find(cflag==col);
        ind(ERRS(rflag(ind),col)>0) = [];
        VALSordMIN = min(VALSord(isfinite(VALSord(:,col)),col));
        if isempty(VALSordMIN) || isempty(ind), continue, end
        VALSord(rflag(ind),col) = VALSordMIN;
    end
end
ERRS(ERRSneg) = MIN * 10.^(VALSord(ERRSneg)+ERRS(ERRSneg)+1);
clear ERRSneg VALSord

% Gets VALS' double-precision floating point error (single digit)
ERRSfp = floatPointError(VALSmax,MIN);
clear VALSmax

% Sets very smalls errors to floating-point precision
ERRSmin = (ERRS<=ERRSfp);
if any(ERRS(ERRSmin)>0)
    warning('EVALIT:doublePrecision',...
        ['Changed some ERRS smaller than eps(VALS), that is, smaller '...
        'than the double precision error.'])
end
ERRS(ERRSmin) = ERRSfp(ERRSmin);
clear ERRSfp

% Rounds input ERRS to specified ESF (number of digits)
ERRSord          = floor(log10(ERRS));
ERRSpre          = ERRSord-ESF+1;
ERRSpre(ERRSmin) = ERRSord(ERRSmin);
ERRSten          = 10.^ERRSpre;
flag             = (ERRSten>realmin);
ERRS(flag)       = ceil(ERRS(flag)./ERRSten(flag)).*ERRSten(flag);
clear ERRSord ERRSpre

% Rounds also the input VALS to the same precision as its error ERRS
flag       = flag & isfinite(VALS) & ~ERRSmin;
VALS(flag) = round(VALS(flag)./ERRSten(flag)).*ERRSten(flag);
clear ERRSmin ERRSten flag

% Derivates partials (PDSs) if not given
if isempty(PDS)
    assert(exist('syms','file')==2,'EVALIT:symbolicToolboxNotFound',...
        'Couldn''t be found the Symbolic Math Toolbox.')
    PDS = cell(1,VARSnum);
    for n = 1:VARSnum
        try
            PDS{n} = diff(FUNC,VARS{n},1);
        catch ME
            NE = MException('EVALIT:errorWhileRunning',...
                ['MATLAB could not differentiate FUNC(%s)=%s, with '...
                'respect to %s.'],...
                regexprep(sprintf(',%s',VARSstr{:}),'^.',''),...
                FUNCstr,VARSstr{n});
            ME = addCause(NE,ME);
            throw(ME)
        end
    end, clear n
end

% Checks PDS arguments
[PDSind,PDSstr] = readPDS(PDS,VARSstr);

% Evaluates partials (PDSs)
VALS    = num2cell(VALS,1); % Cell of VALS' columns
partial = zeros(VALSnum,VARSnum);
try
    if isa(PDS{1},'sym')
        partial          = vpa(partial,VPA);
        for n = 1:VARSnum
            partial(:,n) = subs(PDS{n},VARS(PDSind{n}),VALS(PDSind{n}));
        end
        partial          = double(partial);
    else % function_handle
        for n = 1:VARSnum
            partial(:,n) = PDS{n}(VALS{PDSind{n}});
        end
    end
    clear n
catch ME
    if isa(PDS{1},'sym')
        NE = MException('EVALIT:errorWhileRunning',...
            ['An error occurred while evaluating PDS{%d} = %s. That is '...
            'd(FUNC)/d(%s).'],n,PDSstr{n},VARSstr{n});
    else
        switch n
            case 1,    pref = 'st';
            case 2,    pref = 'nd';
            case 3,    pref = 'rd';
            otherwise, pref = 'th';
        end
        NE = MException('EVALIT:errorWhileRunning',...
            ['An error occurred while evaluating PDS{%d} = %s. '...
            '(FUNC''s partial derivative with respect to '...
            'its %d%s argument).'],n,PDSstr{n},n,pref);
    end
    ME = addCause(NE,ME);
    throw(ME)
end

% If is zero, tries higher order approximation
if (higherOrder>2) && isa(PDS{1},'sym')
    try
        partial2 = partial;
        for k = 2:higherOrder
            [r,c] = find(partial2==0);
            if ~isempty(r)
                itHelps = false;
                for col = unique(c(:)')
                    if PDS{col}==0, continue, end;
                    ind         = find(c==col);
                    PDS2        = diff(PDS{col},VARS{col},1);
                    partial2    = vpa(zeros(length(ind),1),VPA);
                    partial2(:) = subs(PDS2,VARS(PDSind{col}),...
                        cellfun(@(x)x(r(ind)),VALS(PDSind{col}),...
                        'UniformOutput',false));
                    partial(r(ind),col) = ...
                        abs(double(partial2)).*ERRS(r(ind),col).^(k-1)/...
                        factorial(k);
                    itHelps = itHelps || any(partial(r(ind),col));
                end, clear col
                if itHelps
                    warning('EVALIT:errorWhileRunning',...
                        ['Higher order error approximations were used '...
                        'to avoid zeros.'])
                end
            end
        end, clear k
        clear r c itHelps ind PDS2 partial2
    catch ME
        NE = MException('EVALIT:errorWhileRunning',...
            'Could not use higher order approximations.');
        ME = addCause(NE,ME);
        warning(ME.identifier,ME.message)
    end
end

% Evaluates FUNC
VAL = zeros(VALSnum,1);
try
    if isa(FUNC,'sym') % Symbolic expression
        VAL    = vpa(VAL,VPA);
        VAL(:) = subs(FUNC,VARS(FUNCind),VALS(FUNCind));
        VAL    = double(VAL);
    else               % Function_handle
        VAL    = FUNC(VALS{FUNCind});
    end
catch ME
    NE = MException('EVALIT:errorWhileRunning',...
        'An error occurred while evaluating FUNC(%s) = %s.', ...
        regexprep(sprintf(',%s',VARSstr{:}),'^.',''),FUNCstr);
    ME = addCause(NE,ME);
    throw(ME)
end

% Gets errors contributions
ECS         = abs(partial).*ERRS;
ECSmin      = (ECS==0);
ECS(ECSmin) = eps(0);
clear partial

% Gets propagated error
if strcmp(MET,'sum')
    ERR = sum(ECS,2);
else %        'rss'
    ERR = sqrt(sum(ECS.^2,2));
end
ERR(ERR==0) = eps(0);

% Checks errors propagated to INFinity
if (strcmp(MET,'rss') && all(isfinite(ECS(:))) && any(~isfinite(ERR(:))))
    warning('EVALIT:infiniteError',...
        'Error propagated beyond REALMAX. Try ''sum'' instead.')
end

% Rounds error contributions (ESF digits)
ECSord         = floor(log10(ECS));
ECSpre         = ECSord-ESF+1;
ECSpre(ECSmin) = ECSord(ECSmin);
ECSten         = 10.^ECSpre;
flag           = (ECSten>realmin);
ECS(flag)      = round(ECS(flag)./ECSten(flag)).*ECSten(flag);
ECSord         = floor(log10(ECS));
ECSpre         = ECSord-ESF+1;
ECSpre(ECSmin) = ECSord(ECSmin);
clear flag ECSten ECSmin

% Gets VAL's biggest magnitude
noComplex = isreal(VAL);
if noComplex
    VALmax = abs(VAL);
else
    VALmax = max(abs(real(VAL)),abs(imag(VAL)));
end

% Gets VAL's double-precision floating point error (single digit)
ERRfp = floatPointError(VALmax,MIN);

% Sets very small propagated error to floating-point precision
ERRmin      = (ERR<=ERRfp);
ERR(ERRmin) = ERRfp(ERRmin);
clear ERRfp

% Rounds output ERR (ESF digits)
ERRord         = floor(log10(ERR));
ERRpre         = ERRord-ESF+1;
ERRpre(ERRmin) = ERRord(ERRmin);
TEN            = 10.^ERRpre;
flag           = (TEN>realmin);
ERR(flag)      = ceil(ERR(flag)./TEN(flag)).*TEN(flag);
ERRord         = floor(log10(ERR));
ERRpre         = ERRord-ESF+1;
ERRpre(ERRmin) = ERRord(ERRmin);
TEN            = 10.^ERRpre;
flag           = (TEN>realmin);
clear ERRmin

% Rounds output VAL to the same precision as its error ERR
PRE       = ERRpre;
flag      = flag & isfinite(VAL);
VAL(flag) = round(VAL(flag)./TEN(flag)).*TEN(flag);
clear TEN flag

% If the error is unknown, neither the result
if noComplex
    VAL(~isfinite(ERR)) = NaN;
else
    VAL(~isfinite(ERR)|isnan(VAL)) = NaN+NaN*1i;
    VAR = real(VAL);
    VAI = imag(VAL);
end
PRE(~isfinite(VAL)) = NaN;

% All unknowns errors and results have unknowns errors (even when Infs)
ERR(~isfinite(ERR)|isnan(VAL)) = NaN;
ECS(~isfinite(ECS))            = NaN;

% Gets VAL's orders of magnitude
ORD       = floor(log10(VALmax));
flag      = (VALmax==0);
ORD(flag) = 0;
PRE(flag) = 0;
if ~noComplex
    ORR = floor(log10(abs(real(VAL))));
    ORI = floor(log10(abs(imag(VAL))));
    ORR(real(VAL)==0) = 0;
    ORI(imag(VAL)==0) = 0;
end
clear VALmax flag

% Gets VAL's significant figures
SIG = ORD-PRE+1;

% Checks second order error
if checkSecond && isa(PDS{1},'sym')
    try
        % Gets second derivative
        PDS2 =  PDS;
        for n = 1:VARSnum
            PDS2{n} = diff(PDS{n},VARS{n},1);
        end, clear n
        % Evalues second order derivatives
        partial2          = zeros(VALSnum,VARSnum);
        partial2          = vpa(partial2,VPA);
        for n = 1:VARSnum
            partial2(:,n) = subs(PDS2{n},VARS(PDSind{n}),VALS(PDSind{n}));
        end, clear n
        partial2          = double(partial2);
        % Gets second order error approximation
        ECS2 = (abs(partial2).*ERRS.^2)/2;
        % Compares 2nd with 1st order approximation
        if any((ECS2>=ECS) & (ECS2>0) & isfinite(ECS2) & isfinite(ECS))
            warning('EVALIT:incorrectApproximation',...
                ['Some second order error approximations are greater than '...
                'the first one.']);
        end
        clear PDS2 partial2 ECS2
    catch ME
        NE = MException('EVALIT:errorWhileRunning',...
            'Could not check the second order error approximation.');
        ME = addCause(NE,ME);
        warning(ME.identifier,ME.message)
    end
end

%%%%%%%%%%%%%%%%%%%%%
% CHECKS OUTPUTS %%%%%
%%%%%%%%%%%%%%%%%%%%%%%

% Displays result
if OUT
    disp(' ')
    fprintf(1,'EVALIT: FUNC(%s) = %s\n',...
        regexprep(sprintf(',%s',VARSstr{:}),'^.',''),FUNCstr)
    disp('--------------------------------------------------')
    disp(' value +/- error (arguments'' error contributions) ')
    disp('--------------------------------------------------')
    % Results and errors format
    if noComplex
        VALcharWidth = ORD.*(ORD>0) - PRE.*(PRE<0) + (PRE<0) + 2;
        VALcharWidth(isnan(VAL)) = 4;
    else
        VARcharWidth = ORR.*(ORR>0) - PRE.*(PRE<0) + (PRE<0) + 2;
        VAIcharWidth = ORI.*(ORI>0) - PRE.*(PRE<0) + (PRE<0) + 1;
        VARcharWidth(isnan(VAL)) = 4;
        VAIcharWidth(isnan(VAL)) = 4;
    end
    VALdecPlaces             = -PRE.*(PRE<0);
    VALdecPlaces(isnan(VAL)) = 0;
    ERRcharWidth             = ERRord.*(ERRord>0)-PRE.*(PRE<0)+(PRE<0)+1;
    ERRcharWidth(isnan(ERR)) = 3;
    ERRdecPlaces             = -ERRpre.*(ERRpre<0);
    ERRdecPlaces(isnan(ERR)) = 0;
    % Extra spaces to align quantities
    if noComplex
        VALbef = VALcharWidth-VALdecPlaces-(VALdecPlaces>0)-2;
        VALbef = max(VALbef) - VALbef + 1;
    else
        VARbef = VARcharWidth-VALdecPlaces-(VALdecPlaces>0)-2;
        VARbef = max(VARbef) - VARbef + 1;
        VAIbef = VAIcharWidth-VALdecPlaces-(VALdecPlaces>0)-2;
        VAIbef = max(VAIbef) - VAIbef + 1;
    end
    VALaft = max(VALdecPlaces) - VALdecPlaces + ...
        (VALdecPlaces==0).*any(VALdecPlaces>0)+1;
    ERRbef = ERRcharWidth-ERRdecPlaces-(ERRdecPlaces>0)-2;
    ERRbef = max(ERRbef) - ERRbef + 1;
    ERRaft = max(ERRdecPlaces) - ERRdecPlaces + ...
        (ERRdecPlaces==0).*any(ERRdecPlaces>0)+1;
    % Errors contribution format
    ECScharWidth =  ECSord.*(ECSord>0)-ECSpre.*(ECSpre<0)+(ECSpre<0)+1;
    ECScharWidth(isnan(ECS)) = 3;
    ECSdecPlaces             = -ECSpre.*(ECSpre<0);
    ECSdecPlaces(isnan(ECS)) = 0;
    % Extra space
    ECSbef = ECScharWidth-ECSdecPlaces-(ECSdecPlaces>0);
    ECSbef = repmat(max(ECSbef,[],1),VALSnum,1) - ECSbef + 1;
    ECSaft = repmat(max(ECSdecPlaces,[],1),VALSnum,1) - ECSdecPlaces + ...
        (ECSdecPlaces==0).*repmat(any(ECSdecPlaces,1),VALSnum,1)+1;
    for n = 1:VALSnum
        if noComplex
            theFormat = sprintf('%*s%%%d.%df%*s+/-%*s%%%d.%df%*s(',...
                VALbef(n),' ',...
                VALcharWidth(n),VALdecPlaces(n),...
                VALaft(n),' ',ERRbef(n),' ',...
                ERRcharWidth(n),ERRdecPlaces(n),...
                ERRaft(n),' ');
        else
            theFormat = sprintf(...
                '(%*s%%%d.%df%*s%%s%*s%%%d.%dfi%*s) +/-%*s%%%d.%df%*s(',...
                VARbef(n),' ',...
                VARcharWidth(n),VALdecPlaces(n),...
                VALaft(n),' ',VAIbef(n),' ',...
                VAIcharWidth(n),VALdecPlaces(n),...
                VALaft(n),' ',ERRbef(n),' ',...
                ERRcharWidth(n),ERRdecPlaces(n),...
                ERRaft(n),' ');
        end
        for m = 1:VARSnum
            theFormat = sprintf('%s%*s%%%d.%df%*s+',...
                theFormat,ECSbef(n,m),' ',...
                ECScharWidth(n,m),ECSdecPlaces(n,m),...
                ECSaft(n,m),' ');
        end, clear m
        theFormat(end) = [];
        theFormat = sprintf('%s)\n',theFormat);
        ecs = num2cell(ECS(n,:));
        % Prints
        if noComplex
            fprintf(1,theFormat,VAL(n),ERR(n),ecs{:})
        else
            if (VAI(n)<0)
                isign = '-';
            else
                isign = '+';
            end
            fprintf(1,theFormat,VAR(n),isign,abs(VAI(n)),ERR(n),ecs{:})
        end
    end, clear n
    disp('--------------------------------------------------')
    disp(' ')
    clear VALcharWidth VARcharWidth VAIcharWidth VALdecPlaces
    clear ERRcharWidth ERRdecPlaces
    clear ECScharWidth ECSdecPlaces ecs noComplex
    clear VAR VAI ORR ORI VARbef VARaft VAIbef VAIaft isign theFormat
    clear VALbef VALaft ERRbef ERRaft ECSbef ECSaft
end

% How many outputs?
switch nargout
    case 0, clear VAL ERR PRE SIG ORD ECS PDS
    case 1, clear ERR PRE SIG ORD ECS PDS
    case 2, clear PRE SIG ORD ECS PDS
    case 3, clear SIG ORD ECS PDS
    case 4, clear ORD ECS PDS
    case 5, clear ECS PDS
    case 6, clear PDS
end

end % evalit (MAIN FUNCTION)


%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTION  %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function E = floatPointError(X,MIN)
% Gets floating point errors from number X

E       = MIN*eps(X);     % EPS(X) is the separation between X and the next
flag    = (E==0);         % floating point number. So X's true value is
E(flag) = eps(X(flag));   % inside the range +/-0.5*eps(X).
Eord    = floor(log10(E));
Eten    = 10.^Eord;
Etmp    = E;
E       = ceil(E./Eten).*Eten;
flag    = ~isfinite(E);
E(flag) = Etmp(flag);

end % floatPointError

%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTION  %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function [OUT,ESF,MET,PDS] = readOptions(OUT,ESF,MET,PDS,varargin)
% Read optional inputs. Each one are distinguished just by its type.

% Main loop to read optional inputs
for k = 1:length(varargin)
    if isempty(varargin{k}), continue, end
    switch class(varargin{k})
        case 'logical'                             % Print it?
            OUT = varargin{k}(:).';
        case {'double','single','int8','uint8','int16','uint16',...
                'int32','uint32','int64','uint64'} % Sig. Figures?
            ESF = double(varargin{k}(:).');
        case 'char'                                % Method?
            MET = lower(varargin{k}(:).');
        case 'sym'                                 % Partial derivative?
            PDS = num2cell(varargin{k}(:).');
        case 'function_handle'                     % Partial derivative?
            PDS = {varargin{k}};
        case 'cell'                                % Partial derivative?
            PDS = varargin{k}(:).';
        otherwise
            warning('EVALIT:incorrectInput',...
                'Unrecognized optional input number %d.',k)
    end
end

% Checks OUT to be a single logical
assert(length(OUT)==1,'EVALIT:incorrectInput',...
    'Logical optional input must be a single value.')

% Checks ESF to be an integer number
assert((length(ESF)==1) && isfinite(ESF) && isreal(ESF) && (ESF>0) && ...
    (ESF==round(ESF)),'EVALIT:incorrectInput',...
    'Numerical optional input must be a single natural number.')

% Checks MET to be 'rss' or 'sum'
assert(strncmp(MET,'rss',min(length(MET),3)) || ...
    strncmp(MET,'sum',min(length(MET),3)),'EVALIT:incorrectInput',...
    'String optional input must be one of ''rss'' or ''sum''.')

% Checks PDS to be cell of only unique symbolic or only anonymous
% function_handles
if isempty(PDS), return, end
assert((all(cellfun(@(x)isa(x,'sym'),PDS)) && ...
    all(cellfun(@(x)numel(x)==1,PDS))) || ...
    (all(cellfun(@(x)isa(x,'function_handle'),PDS)) && ...
    all(cellfun(@(x)strcmp(getfield(functions(x),'type'),'anonymous'),...
    PDS))),...
    'EVALIT:incorrectInput',...
    ['PDS optional input must be a cell array of only single symbolic '...
    'expressions or only anonymous function_handles.'])

end % readOptions

%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTION  %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function [ERRS,ERRSneg] = readERRS(ERRS,VARSnum,VALSnum)
% Checks ERRS to be a real numerical matrix with integer negative values,
% same size as VALS, and forces it to be double precision.  Or it can be a
% cell array with the arguments', VARS, errors.

% Empty input? Double precision error
if isempty(ERRS)
    ERRS    = zeros(VALSnum,VARSnum);
    ERRSneg = false(VALSnum,VARSnum);
    return
end

% Cell or matrix input?
assert(iscell(ERRS) || isnumeric(ERRS),...
    'EVALIT:incorrectInput',...
    'ERRS must be empty, a cell array or a numerical one.')

% Checks size and changes cell to matrix
if isnumeric(ERRS)
    % Matrix of numbers input
    ERRSsize = size(ERRS);
    assert((length(ERRSsize)==2) && (all(ERRSsize==[1 1]) || ...
        all(ERRSsize==[1 VARSnum]) || all(ERRSsize==[VARSnum 1]) || ...
        all(ERRSsize==[VALSnum VARSnum])),...
        'EVALIT:incorrectInput',...
        'When ERRS is a matrix, it must have the same size as VALS.')
    
    % Forces vectors to be horizontal
    if all(ERRSsize==1)
        % All variables have the same error
        ERRS = repmat(ERRS,VALSnum,VARSnum);
    elseif (prod(ERRSsize)==VARSnum)
        % Each variable have the same error
        ERRS = repmat(ERRS(:).',VALSnum,1);
    end
    
else
    % Cell of numerical inputs
    ERRS = ERRS(:).';
    assert(any(length(ERRS)==[1 VARSnum]) && ...
        all(cellfun(@(x)isnumeric(x),ERRS)),...
        'EVALIT:incorrectInput',...
        ['When ERRS is a cell, it must have as many elements as '...
        'variables in VARS, and all of them should be numerical.'])
    
    % Check the elements' lengths
    lengths = cellfun(@(x)numel(x),ERRS);
    assert(all(ismember(lengths,[0 1 VALSnum])),...
        'EVALIT:incorrectInput',...
        ['When ERRS is a cell, its elements must have a single value '...
        'or at least have the same rows as VALS when not.'])
    
    % Sets empty elements to zero
    ERRS(lengths==0) = repmat({0},1,sum(lengths==0));
    
    % Changes cell to matrix, expanding the singletons elements, arranging
    % them elements by columns
    ERRS =  cell2mat(reshape(cellfun(@(x)x(:).*ones(VALSnum,1),ERRS,...
        'UniformOutput',false),1,VARSnum));
    
end

% Are they real numbers?
assert(isreal(ERRS),'EVALIT:incorrectInput','ERRS data must be real.')

% Negative values are integers?
ERRSneg = (ERRS<0);
assert(all(ERRS(ERRSneg)==round(ERRS(ERRSneg))),...
    'EVALIT:incorrectInput',...
    ['When VALS''s significant figures are provided in ERRS via '...
    'negative values, they must be integer negative numbers.'])

% Not finite errors (INFs and NANs) are set to zero (floating point error)
ERRS(~isfinite(ERRS)) = 0;

% Force it to be double
ERRS = double(ERRS);

end % readERRS

%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTION  %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function VALS = readVALS(VALS,VARSnum)
% Checks VALS to be a numeric matrix with equal columns as number of
% arguments, and forces it to be double precision. Or it can be a cell
% array with the arguments', VARS, values.

% Cell or matrix input?
assert(iscell(VALS)||isnumeric(VALS),'EVALIT:incorrectInput',...
    'VALS must be a cell or a numerical array.')

% Checks size and change from cell to matrix
if isnumeric(VALS)
    % Matrix of numbers input
    assert((ndims(VALS)==2) && ...
        ((size(VALS,2)==VARSnum) || all(size(VALS)==[VARSnum 1])),...
        'EVALIT:incorrectInput',...
        ['When VALS is a matrix, it must have as many columns as '...
        'variables in VARS.'])
    % Forces vectors to be horizontal
    if (size(VALS,2) == 1)
        VALS = VALS.';
    end
else
    % Cell of numerical inputs
    VALS = VALS(:).';
    assert((length(VALS)==VARSnum) && ...
        all(cellfun(@(x)~isempty(x) && isnumeric(x),VALS)),...
        'EVALIT:incorrectInput',...
        ['When VALS is a cell, it must have as many elements as '...
        'variables in VARS, and all of them should be numerical.'])
    
    % Check the elements' lengths
    lengths = cellfun(@(x)numel(x),VALS);
    VALSnum = max(lengths);
    assert(all(ismember(lengths,[1 VALSnum])),...
        'EVALIT:incorrectInput',...
        ['When VALS is a cell, its elements must have a single value '...
        'or at least have the same elements when not.'])
    
    % Changes cell to matrix, expanding the singletons elements, arranging
    % them elements by columns
    VALS =  cell2mat(reshape(cellfun(@(x)x(:).*ones(VALSnum,1),VALS,...
        'UniformOutput',false),1,VARSnum));
    
end

% Forces it to be double
VALS = double(VALS);

end % readVALS

%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTION  %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function VARS = readVARS(VARS)
% Checks VARS to be a vector cell of only symbolic or only strings names.

if ~isempty(VARS)
    switch class(VARS)
        case 'sym'                        % Symbolic array?
            VARS = num2cell(VARS(:).');
            assert(all(cellfun(@(x)numel(x)==1,VARS)),...
                'EVALIT:incorrectInput',...
                'Symbolic VARS elements must be unique.')
        case 'char'                       % Single argument name?
            VARS = {VARS};
            assert(numel(VARS{1})==size(VARS{1},2),...
                'EVALIT:incorrectInput',...
                'Char VARS must be a single line.')
        case 'cell'                       % Cell of syms or chars?
            VARS = VARS(:).';
            assert((iscellstr(VARS) && ...
                all(cellfun(@(x)numel(x)==size(x,2),VARS))) || ...
                (all(cellfun(@(x)isa(x,'sym'),VARS)) && ...
                all(cellfun(@(x)numel(x)==1,VARS))),...
                'EVALIT:incorrectInput',...
                ['VARS must be a cell of symbolic arguments only or of '...
                'strings lines only.'])
        otherwise                         % Other?
            error('EVALIT:incorrectInput','Unrecognized VARS input type.')
    end
end

end % readVARS

%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTION  %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function FUNC = readFUNC(FUNC)
% Checks FUNC to be a single symbolic expression or an anonymous
% function_handle.

assert((numel(FUNC)==1) && ((isa(FUNC,'sym')&&(numel(FUNC)==1)) || ...
    (isa(FUNC,'function_handle') && ...
    strcmp(getfield(functions(FUNC),'type'),'anonymous'))),...
    'EVALIT:incorrectInput',...
    ['FUNC must be a unique symbolic expression or an anonymous '...
    'function_handle.'])

end % readFUNC

%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTION  %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function [VARSstr,FUNCind,FUNCstr] = readVARSnames(FUNC,VARS,PDSnum)
% Checks arguments names

% Check VARS type and get the strings
if ~isempty(VARS)
    switch class(VARS{1})
        case 'sym'
            assert(exist('syms','file')==2,...
                'EVALIT:symbolicToolboxNotFound',...
                'Couldn''t be found the Symbolic Math Toolbox.')
            VARSstr = cellfun(@(x)char(x),VARS,'UniformOutput',false);
        otherwise % 'function_handle'
            VARSstr = VARS;
    end
end
VARSnum = length(VARSstr);

% Checks for repeated VARS
assert(length(unique(VARSstr))==VARSnum,...
    'EVALIT:incorrectInput','Found repeated arguments on VARS.')

% Checks VARS names
assert(all(cellfun(@(x)~isempty(...
    regexp(x,'^[a-zA-Z][a-zA-Z0-9_]*?$','once')),VARSstr)),...
    'EVALIT:incorrectInput','Found incorrect VARS name(s).')

% Gets FUNC arguments and string
switch class(FUNC)
    case 'sym'
        assert(exist('syms','file')==2,'EVALIT:symbolicToolboxNotFound',...
            'Couldn''t be found the Symbolic Math Toolbox.')
        FUNCstr = char(FUNC);
        try
            FUNCvar = symvar(FUNC);
        catch ME
            NE = MException('EVALIT:errorWhileRunning',...
                'Couldn''t get arguments from FUNC = %s.',FUNCstr);
            ME = addCause(NE,ME);
            throw(ME)
        end
        FUNCvar = cellfun(@(x)char(x),num2cell(FUNCvar),...
            'UniformOutput',false);
    otherwise % 'function_handle'
        tmp = functions(FUNC);
        % Get everything after this char
        aftThis = '\)';
        FUNCstr = regexp(tmp.function,...
            ['(?<=' aftThis ').+$'],'once','match');
        % Get everything between these chars
        aftThis = '\(';
        befThis = '\)';
        FUNCvar = regexp(tmp.function,...
            ['(?<=' aftThis ')(.*?)(?=' befThis ')'],'match','once');
        % Split delimited text by this flag
        flag = ',';
        FUNCvar = regexp(FUNCvar,flag,'split');
end
FUNCnum = length(FUNCvar);

% Compares FUNCvar with VARSstr
assert((FUNCnum==VARSnum) && all(ismember(FUNCvar,VARSstr)),...
    'EVALIT:incorrectInput',...
    'FUNC should have the same arguments as in VARS.')

% Gets the order
VARSmat = strvcat(VARSstr);
FUNCind = cellfun(@(x)strmatch(x,VARSmat,'exact'),FUNCvar);

% Checks PDS and FUNC
assert((PDSnum>0) || isa(FUNC,'sym'),'EVALIT:incorrectInput',...
    'When PDSs are not provided FUNC must be a Symbolic expression.')

% Checks PDS and VARS
assert((PDSnum==0) || (PDSnum==VARSnum),'EVALIT:incorrectInput',...
    'When provided, the number of PDSs must be equal to VARS'' number.')

end % readVARSnames

%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTION  %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%

function [PDSind,PDSstr] = readPDS(PDS,VARSstr)
% Check PDS arguments

% Number of partials
PDSnum = length(PDS);

% Initializations
PDSind = cell(1,PDSnum);
PDSstr = PDSind;
PDSvar = PDSind;

% Reads PDS arguments
switch class(PDS{1})
    case 'sym'
        assert(exist('syms','file')==2,'EVALIT:symbolicToolboxNotFound',...
            'Couldn''t be found the Symbolic Math Toolbox.')
        for k = 1:PDSnum
            PDSstr{k} = char(PDS{k});
            try
                PDSvar{k} = symvar(PDS{k});
            catch ME
                NE = MException('EVALIT:errorWhileRunning',...
                    'Couldn''t get arguments from PDS{%d} = %s.',k,...
                    PDSstr{k});
                ME = addCause(NE,ME);
                throw(ME)
            end
            PDSvar{k} = cellfun(@(x)char(x),num2cell(PDSvar{k}),...
                'UniformOutput',false);
        end
    otherwise % 'function_handle'
        for k = 1:PDSnum
            tmp = functions(PDS{k});
            % Get everything after this char
            aftThis   = '\)';
            PDSstr{k} = regexp(tmp.function,...
                ['(?<=' aftThis ').+$'],'once','match');
            % Get everything between these chars
            aftThis   = '\(';
            befThis   = '\)';
            PDSvar{k} = regexp(tmp.function,...
                ['(?<=' aftThis ')(.*?)(?=' befThis ')'],'match','once');
            % Split delimited text by this char
            delThis   = ',';
            PDSvar{k} = regexp(PDSvar{k},delThis,'split');
        end
end

% Compares PDSvar with VARSstr
VARSnum =  length(VARSstr);
VARSmat = strvcat(VARSstr);
for k = 1:PDSnum
    assert(length(PDSvar{k})<=VARSnum,'EVALIT:incorrectInput',...
        'PDS{%d} has more input arguments than the provided in VARS.',k)
    assert(all(ismember(PDSvar{k},VARSstr)),'EVALIT:incorrectInput',...
        'PDS have different arguments than the provided in VARS.')
    
    % Gets the order
    PDSind{k} = cellfun(@(x)strmatch(x,VARSmat,'exact'),PDSvar{k});
end

end % readVARSnames

% [EOF] EVALIT.M by CARLOS VARGAS

Contact us