# HG changeset patch
# User Nicholas R. Jankowski <jankowski.nicholas@gmail.com>
# Date 1708715651 18000
#      Fri Feb 23 14:14:11 2024 -0500
# Node ID 655e4efd19f6b37d0844940c73f944760a49c20b
# Parent  e073e676c4162aeb7317beb3ebd9a24009219a09
unqiue.m: Implement third output when 'stable' sort is used. (bug #65176)

* scripts/set/unique.m: Implement third output mapping Y back to X when
'stable' sort option is used. Refactor code for second and third outputs to
remove redundancy.  Replace repeated calls to sort in second and third output
code with index expressions. Remove unimplemented note from docstring. Update
multi-output BISTs for stable with new correct output.  Add additional BISTs to
test 'stable' output with other options.
* etc/NEWS.10.md: Note implementation under Matlab compatibility section.

diff -r e073e676c416 -r 655e4efd19f6 etc/NEWS.10.md
--- a/etc/NEWS.10.md	Thu Feb 22 15:20:16 2024 -0500
+++ b/etc/NEWS.10.md	Fri Feb 23 14:14:11 2024 -0500
@@ -39,6 +39,8 @@
 vector and the other is a column vector, the interpolating points are processed
 through meshgrid and the output is a matrix the same size as the meshgrid.
 
+- Enable the third output for `unique` with the `stable` sort option.
+
 ### Alphabetical list of new functions added in Octave 10
 
 * `rticklabels`
diff -r e073e676c416 -r 655e4efd19f6 scripts/set/unique.m
--- a/scripts/set/unique.m	Thu Feb 22 15:20:16 2024 -0500
+++ b/scripts/set/unique.m	Fri Feb 23 14:14:11 2024 -0500
@@ -84,9 +84,6 @@
 ## outputs @var{i}, @var{j} will follow the shape of the input @var{x} rather
 ## than always being column vectors.
 ##
-## The third output, @var{j}, has not been implemented yet when the sort
-## order is @qcode{"stable"}.
-##
 ## @seealso{uniquetol, union, intersect, setdiff, setxor, ismember}
 ## @end deftypefn
 
@@ -230,36 +227,86 @@
     endif
   endif
 
-  ## Calculate j output (3rd output)
-  if (nargout > 2)
-    j = i;  # cheap way to copy dimensions
-    j(i) = cumsum ([1; ! match(:)]);
-    if (! optsorted)
-      warning ("unique: third output J is not yet implemented");
-      j = [];
-    endif
-
-    if (optlegacy && isrowvec)
-      j = j.';
-    endif
-  endif
-
   ## Calculate i output (2nd output)
   if (nargout > 1)
+
     if (optsorted)
+
       idx = find (match);
+
       if (! optlegacy && optfirst)
         idx += 1;   # in-place is faster than other forms of increment
       endif
+
+      if (nargout > 2)
+        j = i;  # cheap way to copy dimensions
+        j(i) = cumsum (! [false; match(:)]);
+      endif
+
       i(idx) = [];
+
     else
-      i([false; match(:)]) = [];
-      ## FIXME: Is there a way to avoid a call to sort?
-      i = sort (i);
+
+      ## Get inverse of sort index i so that sort(x)(k) = x.
+      k = i;  # cheap way to copy dimensions
+      k(i) = 1:n;
+
+      if (nargout > 2)
+
+        ## Generate logical index of sorted unique value locations.
+        nomatch = ! [false; match(:)];
+
+        ## Calculate i output as those locations remapped to unsorted x.
+        i_out = find (nomatch(k));
+
+        ## Find the linear indexes of the unique elements of sort(x).
+        u = find (nomatch);
+
+        ## Find unique indexes for all element locations .
+        l = u(cumsum (nomatch));
+
+        ## l(k) gives u element locations mapped back to unsorted x. E.g.,
+        ## for x = [40,20,40,20,20,30,10]' # data
+        ## i =   [7,2,4,5,6,1,3]'          # sort(x) index, x(i) = sort(x)
+        ## nomatch = [1,1,0,0,1,1,0]'      # logical sorted index of unique x
+        ## i_out = [1,2,6,7]'              # unique output index, y = x(i_out)
+        ## k = [6,2,7,3,4,5,1]'            # inverse idx of i, sort(x)(k) = x
+        ## l = [1,2,2,2,5,6,6]'            # unique elem. to reproduce sort(x)
+        ## l(k) = [6,2,6,2,2,5,1]'         # unique elements to reproduce (x)
+        ## i(l(k)) =  [1,2,1,2,2,6,7]'     # unique elem. mapped to sort(x)idx
+        ##
+        ## i_out == i(l(k))' broadcasts to:
+        ##  [ 1  1  0  0  0  0  0
+        ##    0  0  1  1  1  0  0
+        ##    0  0  0  0  0  1  0
+        ##    0  0  0  0  0  0  1 ]
+        ## Row value of each column maps i(l(k)) to i_out, gives j for y(j)=x
+
+        ## FIXME: 2-D projection to find j increases largest stored element
+        ##        from n to m x n (m = numel (y)) and uses slower find
+        ##        codepath.  Ideally would be replaced by direct linear or
+        ##        logical indexing.
+
+        [j,~] = find (i_out == i(l(k))');
+
+        ## Replace full x-> sort(x) output with duplicates removed.
+        i = i_out;
+      else
+
+        ## Find i output as unique value locations remapped to unsorted x.
+        i = find (! [false; match(:)](k));
+
+      endif
+
     endif
 
     if (optlegacy && isrowvec)
       i = i.';
+
+      if (nargout > 2)
+        j = j.';
+      endif
+
     endif
   endif
 
@@ -302,11 +349,10 @@
 %! assert (j, [1;1;2;3;3;3;4]);
 
 %!test
-%! [y,i,~] = unique ([4,4,2,2,2,3,1], "stable");
+%! [y,i,j] = unique ([4,4,2,2,2,3,1], "stable");
 %! assert (y, [4,2,3,1]);
 %! assert (i, [1;3;6;7]);
-%! ## FIXME: 'j' input not calculated with stable
-%! ##assert (j, []);
+%! assert (j, [1;1;2;2;2;3;4]);
 
 %!test
 %! [y,i,j] = unique ([1,1,2,3,3,3,4]', "last");
@@ -335,11 +381,10 @@
 
 %!test
 %! A = [4,5,6; 1,2,3; 4,5,6];
-%! [y,i,~] = unique (A, "rows", "stable");
+%! [y,i,j] = unique (A, "rows", "stable");
 %! assert (y, [4,5,6; 1,2,3]);
 %! assert (A(i,:), y);
-%! ## FIXME: 'j' output not calculated correctly with "stable"
-%! ##assert (y(j,:), A);
+%! assert (y(j,:), A);
 
 ## Test "legacy" option
 %!test
@@ -355,6 +400,52 @@
 %! assert (i, [2; 5; 4; 3]);
 %! assert (j, [4; 1; 4; 3; 2]);
 
+%!test <*65176>
+%! a = [3 2 1 2; 2 1 2 1];
+%! [o1, o2, o3] = unique (a);
+%! assert ({o1, o2, o3}, {[1;2;3], [4;2;1], [3;2;2;1;1;2;2;1]});
+%! [o1, o2, o3] = unique (a, "stable");
+%! assert ({o1, o2, o3}, {[3;2;1], [1;2;4], [1;2;2;3;3;2;2;3]})
+
+%!test <*65176>
+%! a = [3 2 1 2; 2 1 2 1];
+%! [o1, o2, o3] = unique (a(1,:), "rows");
+%! assert ({o1, o2, o3}, {a(1,:), 1, 1});
+%! [o1, o2, o3] = unique (a(1,:), "rows", "stable");
+%! assert ({o1, o2, o3}, {a(1,:), 1, 1});
+%! [o1, o2, o3] = unique (a, "rows");
+%! assert ({o1, o2, o3}, {[a(2,:); a(1,:)], [2;1], [2;1]});
+%! [o1, o2, o3] = unique (a, "rows", "stable");
+%! assert ({o1, o2, o3}, {a, [1;2], [1;2]});
+%! [o1, o2, o3] = unique ([a;a], "rows");
+%! assert ({o1, o2, o3}, {[a(2,:); a(1,:)], [2;1], [2;1;2;1]});
+%! [o1, o2, o3] = unique ([a;a], "rows", "stable");
+%! assert ({o1, o2, o3}, {a, [1;2], [1;2;1;2]});
+
+%!test <*65176>
+%! a = gallery ("integerdata", [-100, 100], 6, 6);
+%! a = [a(2,:); a(1:5,:); a(2:6,:)];
+%! [o1, o2, o3] = unique (a);
+%! assert ({o1, o1(o3), o2, o3}, {a(:)(o2), a(:), ...
+%! [26;22;34;45;57; 6;11;17;33;28;35;15;56; 2;59; 4;66; ...
+%!  16;50;49;27;24;37;44;48;39;38;13;23; 5;12;46;55; 1], ...
+%! [34;14;34;16;30; 6;34;16;30; 6; 7;31;28;31;12;18; 8;31;12;18; 8; 2;29; ...
+%!  22;29; 1;21;10;29; 1;21;10; 9; 3;11; 3;23;27;26; 3;23;27;26;24; 4;32; ...
+%!  4; 25;20;19; 4;25;20;19;33;13; 5;13;15; 2;24;13;15; 2;24;17]});
+%! [o1, o2, o3] = unique (a, "stable");
+%! assert ({o1, o1(o3), o2, o3}, {a(:)(o2), a(:), ...
+%! [ 1; 2; 4; 5; 6;11;12;13;15;16;17;22;23;24;26;27;28; ...
+%!  33;34;35;37;38;39;44;45;46;48;49;50;55;56;57;59;66], ...
+%! [ 1; 2; 1; 3; 4; 5; 1; 3; 4; 5; 6; 7; 8; 7; 9;10;11; 7; 9;10;11;12;13; ...
+%!  14;13;15;16;17;13;15;16;17;18;19;20;19;21;22;23;19;21;22;23;24;25;26;...
+%!  25;27;28;29;25;27;28;29;30;31;32;31;33;12;24;31;33;12;24;34]});
+%! [o1, o2, o3] = unique (a, "rows");
+%! assert ({o1, o1(o3,:), o2, o3}, {a(o2,:), a, ...
+%! [6;11;2;4;5;1], [6;3;6;4;5;1;6;4;5;1;2]});
+%! [o1, o2, o3] = unique (a, "rows", "stable");
+%! assert ({o1, o1(o3,:), o2, o3}, {a(o2,:), a, ...
+%! [1;2;4;5;6;11], [1;2;1;3;4;5;1;3;4;5;6]});
+
 ## Test input validation
 %!error <Invalid call> unique ()
 %!error <X must be an array or cell array of strings> unique ({1})
@@ -376,6 +467,4 @@
 %!error <invalid option> unique ({"a", "b", "c"}, "rows", "UnknownOption2")
 %!error <invalid option> unique ({"a", "b", "c"}, "UnknownOption1", "last")
 %!warning <"rows" is ignored for cell arrays> unique ({"1"}, "rows");
-%!warning <third output J is not yet implemented>
-%! [y,i,j] = unique ([2,1], "stable");
-%! assert (j, []);
+
