######################################################################## ## ## Copyright (C) 2020 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or . ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## . ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{c} =} uniquetol (@var{A}) ## @deftypefnx {} {@var{c} =} uniquetol (@var{A}, @var{tol}) ## @deftypefnx {} {@var{c} =} uniquetol (@dots{}, @var{property}, @var{value}) ## @deftypefnx {} {[@var{c}, @var{ia}, @var{ic}] =} uniquetol (@dots{}) ## Return the unique elements of @var{A} within tolerance @var{tol}. ## ## Two values, @var{x} and @var{y}, are within relative tolerance if ## @code{abs (@var{x} - @var{y}) <= @var{tol} * max (abs (@var{A}(:)))}. ## ## The input @var{A} must be a floating point type (double or single). ## ## If @var{tol} is unspecified, the default tolerance is 1e-12 for double ## precision input or 1e-6 for single precision input. ## ## The function may also be called with the following optional property/value ## pairs. Property/value pairs must be passed after other input arguments: ## ## @table @asis ## @item @qcode{"ByRows"} (default: @code{false}) ## When true, return the unique rows of @var{A}. @var{A} must be a 2-D array ## to use this option. For rows, the criteria for uniqueness is changed to ## @code{all (abs (@var{x} - @var{y}) <= @var{tol}*max (abs (@var{A}),[],1))} ## which compares each column component of a row against a column-specific ## tolerance. ## ## @item @qcode{"DataScale"} ## The tolerance test is changed to ## @code{abs (@var{x} - @var{y}) <= @var{tol}*@var{DS}} where @var{DS} is a ## scalar unless the property @qcode{"ByRows"} is true. In that case, @var{DS} ## can either be a scalar or a vector with a length equal to the number of ## columns in @var{A}. Using a value of @code{1.0} for @var{DS} will change ## the tolerance from a relative one to an absolute tolerance. Using a value ## of @code{Inf} will disable testing. ## ## @item @qcode{"OutputAllIndices"} (default: @code{false}) ## When true, @var{ia} is a cell array (not a vector) that contains the indices ## for @emph{all} elements in @var{A} that are within tolerance of a value in ## @var{C}. That is, each cell in @var{ia} corresponds to a single unique ## value in @var{C}, and the values in each cell correspond to locations in ## @var{A}. ## @end table ## ## The output @var{c} is a row vector if the input @var{A} is a row vector. ## For all other cases, a column vector is returned. ## ## The optional output @var{ia} is a column index vector such that ## @code{@var{c} = @var{A}(@var{ia})}. If the @qcode{"ByRows"} property is ## true, the condition is @code{@var{c} = @var{A}(@var{ia}, :)}. If the ## @qcode{"OutputAllIndices"} property is true, then the values ## @code{@var{A}(@var{ia}@{@var{i}@})} are all within tolerance of the unique ## value @code{@var{c}(@var{i})}. ## ## The optional output @var{ic} is a column index vector such that ## @code{@var{A} = @var{c}(@var{ic})} when @var{A} is a vector. When @var{A} ## is a matrix, @code{@var{A}(:) = @var{c}(@var{ic})}. If the @qcode{"ByRows"} ## property is true then @code{@var{A} = @var{c}(@var{ic},:)}. ## ## Example: small round-off errors require @code{uniquetol}, not @code{unique} ## ## @example ## @group ## x = [1:5]; ## ## Inverse_Function (Function (x)) should return exactly x ## y = exp (log (x)); ## D = unique ([x, y]) ## @result{} [1.0000 2.0000 3.0000 3.0000 4.0000 5.0000 5.0000] ## C = uniquetol ([x, y]) ## @result{} [1 2 3 4 5] ## @end group ## @end example ## ## @seealso{unique, union, intersect, setdiff, setxor, ismember} ## @end deftypefn function [c, ia, ic] = uniquetol (A, varargin) if (nargin < 1) print_usage (); endif if (isempty (A)) c = A; ia = []; ic = []; return; endif if (! isfloat (A)) error ("Octave:uniquetol:unsupported-type", "uniquetol: A must be a double or single precision array"); endif if (nargin == 1 || ischar (varargin{1})) tol = ifelse (isa (A, "double"), 1e-12, 1e-6); elseif (! (isfloat (varargin{1}) && isscalar (varargin{1}))) error ("Octave:uniquetol:unsupported-type", "uniquetol: TOL must be a double or single precision scalar"); else tol = varargin{1}; varargin(1) = []; endif if (mod (numel (varargin), 2)) error ("uniquetol: PROPERTY/VALUE arguments must be passed in pairs"); endif by_rows = false; output_all_indices = false; data_scale = []; for k = 1:2:numel (varargin) if (! ischar (varargin{k})) error ("uniquetol: PROPERTY must be a string"); endif if (strcmpi (varargin{k}, "ByRows")) by_rows = logical (varargin{k+1}); if (by_rows && ndims (A) > 2) error ('uniquetol: A must be a 2-D array when "ByRows" is true'); endif elseif (strcmpi (varargin{k}, "OutputAllIndices")) output_all_indices = logical (varargin{k+1}); elseif (strcmpi (varargin{k}, "DataScale")) data_scale = varargin{k+1}(:).'; if (! isfloat (data_scale) || any (data_scale(:) < 0) || any (isnan (data_scale(:)))) error ("uniquetol: DataScale must be a non-NaN, positive floating point scalar or vector"); endif cols_data_scale = columns (data_scale); if (cols_data_scale != 1 && cols_data_scale != columns (A)) error ("uniquetol: invalid DataScale size"); endif else error ("uniquetol: unknown property '%s'", varargin{k}); endif endfor if (isempty (data_scale)) data_scale = max (abs (A(! isinf (A))(:))); endif tol_data_scale = tol * data_scale; if (by_rows) points = rows (A); d = columns (A); Iall = zeros (points, 1); I = NaN (d, 1); ia = {}; J = NaN (d, 1); j = 1; ii = 0; for i = 1:points if (any (Iall == i)) continue; else equ = all (abs (A - A(i,:)) <= tol_data_scale, 2); equ(i,1) = equ(i,1) || any (! isfinite (A(i,:)), 2); sumeq = sum (equ); ia_tmp = find (equ); if (output_all_indices) ia{end+1} = ia_tmp; endif Iall(ii+(1:sumeq)) = ia_tmp; I(j) = ia_tmp(1); J(equ) = j; ii += sumeq; j += 1; endif endfor I(isnan (I)) = []; J(isnan (J)) = []; c = A(I,:); if (! output_all_indices) ia = I(1:j-1); endif ic = J; else isrowvec = isrow (A); A = A(:); lengthA = length(A); isnanA = isnan (A); anyisnanA = any (isnanA); [sortA, sAi] = sort (A); diffsortA = diff (sortA); isinfsortA = isinf(sortA); isnansortA = isnan(sortA); numnan = sum (isnansortA); if any (isinfsortA) sAnin=sortA(!(isinfsortA|isnansortA)); diffsortA(isinf (diffsortA)) = abs (sAnin(end)-sAnin(1))+10; end csdx = cumsum (diffsortA); ue = [true;diff([0;csdx-mod(csdx,tol_data_scale)])>eps(max(csdx))]; ueold = nan; while any (ueold != ue) ueold = ue; belowtol = [false;diff(sortA(ue))= lengthA-numnan) = []; ue(needstomove+1) = true; endif endwhile c = sortA(ue); [~, sortsAi] = sort (sAi); cumsumue = cumsum(ue); ic = cumsumue(sortsAi); if anyisnanA findisnanA = find (isnanA); else findisnanA = []; endif if (output_all_indices) nu=cumsumue(end); ia=cell (1,nu); for k = 1:nu ia{k} = setdiff (sAi(cumsumue==k),findisnanA); endfor else ia = sAi(ue); endif if anyisnanA rowsc1=rows (c)+1; c(rowsc1) = nan; if (output_all_indices) ia(rowsc1)=findisnanA; else ia(rowsc1)=findisnanA(1); endif ic(isnanA)=rowsc1; endif ## Matlab-compatible orientation of output if (isrowvec) c = c.'; endif endif endfunction %!assert (uniquetol ([1 1 2; 1 2 1; 1 1 2+10*eps]), [1;2]) %!assert (uniquetol ([1 1 2; 1 0 1; 1 1 2+10*eps], "byrows", true), %! [1 1 2; 1 0 1]) %!assert (uniquetol ([]), []) %!assert (uniquetol ([1]), [1]) %!xtest <59850> %! ## FIXME: Matlab returns values sorted %! assert (uniquetol ([2, 1]), [1, 2]); %!assert (uniquetol ([1; 2]), [1; 2]) %!xtest <59850> %! ## FIXME: Matlab returns only one unique value for Inf. %! assert (uniquetol ([-Inf, 1, NaN, Inf, NaN, Inf]), [-Inf, 1, Inf, NaN]); %!xtest <59850> %! ## FIXME: Matlab returns empty column vectors. %! ## Do we want to bother with that? %! assert (uniquetol (zeros (1,0)), zeros (0,1)); %!assert (uniquetol (zeros (1,0), "byrows", true), zeros (1,0)) %!assert (uniquetol ([1,2,2,3,2,4], "byrows", true), [1,2,2,3,2,4]) %!assert (uniquetol ([1,2,2,3,2,4]), [1,2,3,4]) %!assert (uniquetol ([1,2,2,3,2,4].', "byrows", true), [1;2;3;4]) %!assert (uniquetol (sparse ([2,0;2,0])), sparse ([0;2])) %!assert (uniquetol (sparse ([1,2;2,3])), sparse ([1;2;3])) %!assert (uniquetol (single ([1,2,2,3,2,4]), "byrows", true), %! single ([1,2,2,3,2,4])) %!assert (uniquetol (single ([1,2,2,3,2,4])), single ([1,2,3,4])) %!assert (uniquetol (single ([1,2,2,3,2,4].'), "byrows", true), %! single ([1;2;3;4])) ## Test index vector return arguments %!test %! [c, ia, ic] = uniquetol ([1,1,2,3,3,3,4]); %! assert (c, [1,2,3,4]); %! assert (ia, [1;3;4;7]); %! assert (ic, [1;1;2;3;3;3;4]); ## Test index vector return arguments with "ByRows" %!test %! A = [2, 3, 4; 2, 3, 4]; %! [c, ia, ic] = uniquetol (A, "byrows", true); %! assert (c, [2, 3, 4]); %! assert (A(ia,:), c); %! assert (c(ic,:), A); %!test %! x = (2:7)'*pi; %! y = exp (log (x)); %! C = uniquetol ([x; y]); %! assert (C, x,1e-12); ## Test "ByRows" Property %!test %! A = [0.06, 0.21, 0.38; 0.38, 0.21, 0.39; 0.54, 0.56, 0.41; 0.46, 0.52, 0.95]; %! B = log (exp (A)); %! C = uniquetol ([A; B], "ByRows", true); %! assert (C, A); ## Test "DataScale" Property %!test %! x = 10^11; %! C = uniquetol ([x, exp(log(x))], 1e-6, "DataScale", 1); %! assert (C, [x, exp(log(x))]); ## Test "OutputAllIndices" Property %!test %! A = [.1 .2 .3 10]; %! [C, ia, ic] = uniquetol (A, .1, "OutputAllIndices", true); %! assert (C, [.1, 10]); %! assert (ia, {(1:3)', 4}); %! assert (ic, [1; 1; 1; 2]); ## Test input validation %!error uniquetol () %!error uniquetol (int8 (1)) %!error uniquetol (1, int8 (1)) %!error uniquetol (1, [1, 2]) %!error uniquetol (1, 2, "byrows") %!error uniquetol (1, 2, 3, "bar") %!error uniquetol (ones(2,2,2), "byrows", true) %!error uniquetol (1, "DataScale", '1') %!error uniquetol (1, "DataScale", -1) %!error uniquetol (1, "DataScale", NaN) %!error uniquetol (1, "DataScale", [1 2]) %!error uniquetol (1, "foo", "bar") %!error uniquetol (1, 2, "foo", "bar")