Code covered by the BSD License

# evalit.m v3 (Jul 2013)

### Carlos Adrian Vargas Aguilera (view profile)

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)
%
%   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.

%   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
clear varargin

% Checks FUNC input

% Checks VARS input
VARSnum = length(VARS);

% Checks arguments names

% Checks VALS input
if isempty(VALS) % No values?
if OUT
disp(' ')
fprintf(1,'EVALIT: FUNC(%s) = %s\n',...
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
VALSnum = size(VALS,1);

% Checks ERRS input

%%%%%%%%%%%%%%%%%%%%%
% 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.'],...
throw(ME)
end
end, clear n
end

% Checks PDS arguments

% 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 '...
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
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.');
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.', ...
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);
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.');
warning(ME.identifier,ME.message)
end
end

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

% Displays result
if OUT
disp(' ')
fprintf(1,'EVALIT: FUNC(%s) = %s\n',...
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'
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.'])

%%%%%%%%%%%%%%%%%%%%%
% 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
all(ERRSsize==[1 VARSnum]) || all(ERRSsize==[VARSnum 1]) || ...
'EVALIT:incorrectInput',...
'When ERRS is a matrix, it must have the same size as VALS.')

% Forces vectors to be horizontal
% All variables have the same error
ERRS = repmat(ERRS,VALSnum,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);

%%%%%%%%%%%%%%%%%%%%%
% 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);

%%%%%%%%%%%%%%%%%%%%%
% 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

%%%%%%%%%%%%%%%%%%%%%
% 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.'])

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

% 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.')
otherwise % 'function_handle'
end
end

% Checks for repeated VARS
'EVALIT:incorrectInput','Found repeated arguments on VARS.')

% Checks VARS names
assert(all(cellfun(@(x)~isempty(...
'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);
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
'EVALIT:incorrectInput',...
'FUNC should have the same arguments as in VARS.')

% Gets the order
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.')

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

% 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});
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
for k = 1:PDSnum
assert(length(PDSvar{k})<=VARSnum,'EVALIT:incorrectInput',...
'PDS{%d} has more input arguments than the provided in VARS.',k)