numcmp(A,FUN,B,TOL)

function [IND,DT] = numcmp(A,FUN,B,TOL)
% NUMCMP  Compares numbers up to a certain precision.
%
%   Syntax:
%               IND = numcmp(A);           % Compares: A==0.0, by default.
%               IND = numcmp(A,FUN);
%               IND = numcmp(A,FUN,B);
%               IND = numcmp(A,FUN,B,TOL);
%          [IND,BA] = numcmp(...);
%
%   Input:
%     A   - First number(s). It can be a single value or an array.
%     FUN - Function to be use on the comparision. One of: 
%            @eq , @ne , @lt , @gt , @le , @ge , ...   or 
%            'eq', 'ne', 'lt', 'gt', 'le', 'ge', ...   or
%            '==', '~=',  '<',  '>', '<=', '>='  . 
%           Default: '==' (equal to).
%     B   - Second number(s). It can be a single value or an array. If A is
%           also an array, both must be of the same size, or one a vector
%           with the same number of columns or rows.
%           Default: 0 (zero).
%     TOL - Integer defining the precision: 10^(-TOL).
%           Default: 1 (compares up to a 0.1 precision) 
%
%   Output:
%     IND  - Logical array with ones where the relation is true and zeros
%            where it is false.
%     BA   - Difference: (B-A)*10^TOL (optional).
%
%   Description:
%     This programs compares two arrays up to a defined precision by
%     rounding the numbers once multiplied by 10^TOL, and eliminated the
%     minimum value from all of them, in order to avoid floating points
%     errors as possible. The user should keep in mind that 
%                max([A(:); B(:)])*10^TOL 
%     should be less than REALMAX.
%
%   Example:
%     x = pi;
%     y = 3.14;
%     x==y               % up to eps (~2.22e-016)
%     numcmp(x,'==',y,2) % up to 0.01
%     % Results:
%     %      0
%     %      1
%
%   See also: STRCMP, EQ, NE, LT, GT, LE, GE and REALMAX. 

%   Copyright 2008 Carlos Adrian Vargas Aguilera
%   $Release: 1.1 $  $Date: 2008/09/11 11:30:00 $

%   Written by
%   M.S. Carlos Adrian Vargas Aguilera
%   Physical Oceanography PhD candidate
%   CICESE 
%   Mexico, 2008
%   [e-mail]
%
%   Download from:
%   http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objec
%   tType=author&objectId=1093874

%   1.1  2008/08/22:    Release
%   1.2  2008/09/05:    Fixed bug when not finite inputs.

% Checks arguments:
if nargin<1
 error('CVARGAS:numcmp:NotEnoughInputs', ...
 'At least one input is required.')
end
[mA,nA] = size(A);
if nargin<2 || isempty(FUN)
 FUN = 'eq';
end
if nargin<3 || isempty(B)
 B = zeros(mA,nA);
end
if nargin<4 || isempty(TOL)
 TOL = 1;
end
if nargin>4
 error('CVARGAS:numcmp:TooManyInputs', ...
 'Maximum 4 inputs are required.')
end

% Forces integer TOL:
TOL = round(TOL);

% Checks comparision function:
if ischar(FUN)
 switch FUN
  case '==', FUN = 'eq';
  case '~=', FUN = 'ne';
  case  '<', FUN = 'lt';
  case  '>', FUN = 'gt';
  case '<=', FUN = 'le';
  case '>=', FUN = 'ge';
  otherwise, % do nothing, hopping is one of the right column above.
 end
end

% Sets the minumum value as zero:
ABmin = min([A(isfinite(A)); B(isfinite(B))]); % Fixed bug SEP 09, 2008
A     = A-ABmin;
B     = B-ABmin;

% Sets B size equal to A size (I avoided the use of BSXFUN):
[mB,nB] = size(B);
if (mB*nB)~=(mA*nA)
 if mB*nB==1
  B = repmat(B,mA,nA);
 elseif nB == nA
  B = repmat(B,mA,1);
 elseif mB == mA
  B = repmat(B,1,nA);
 else
  error('CVARGAS:numcmp:IncorrectInputSize', ...
   'A and B must be arrays of equal sizes, or equal columns or rows.')
 end
end

% Converts to 10^(TOL) positive integers:
A = round(A*10^TOL);
B = round(B*10^TOL);

% Makes comparisions:
IND = feval(FUN,A,B);

% Diference output in 10^(-TOL) values:
if nargout==2
 BA = B-A;
end