* Bugs 16337 16455 fixed: [..,..,ku] = unique(..) implemented
[scilab.git] / scilab / modules / elementary_functions / macros / unique.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
3 // Copyright (C) 2012 - 2016 - Scilab Enterprises
4 // Copyright (C) 2018 - 2020 - Samuel GOUGEON
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function [x, ki, ko, nb] = unique(x, varargin)
14     // extract unique components of a vector, matrix, or hypermatrix
15     // varargin : orient=1|2|"r"|"c", "uniqueNan", "keepOrder"
16     //
17     // History:
18     // * 2019 - S. Gougeon :
19     //   - add uniqueNan option: http://bugzilla.scilab.org/15522
20     //   - add keepOrder option: http://bugzilla.scilab.org/15795
21     //   - add nb output option: http://bugzilla.scilab.org/8418
22     // * 2020 - S. Gougeon :
23     //   - add ku output indices: http://bugzilla.scilab.org/16337
24
25     keepOrder = %f
26     uniqueNan = %f
27     orient = "*"
28     newInf = [] // init Inf substitute in case of "uniqueNan" and or(x==%inf)
29     ki = []
30     ko = []
31     nb = []
32     sx = size(x)
33
34     // CHECKING INPUT ARGUMENTS
35     // ------------------------
36     in = varargin
37     i = 2;  // index of the current input argument
38     if size(in)>0 then
39         a = in(1)
40         if type(a)==0
41             a = "*"
42         end
43         if typeof(a)=="string"
44             a = convstr(a)
45         end
46         select a
47         case 1
48             orient = "r"
49         case 2
50             orient = "c"
51         case "uniquenan"
52             uniqueNan = %t
53         case "keeporder"
54             keepOrder = %t
55         case "r"
56             orient = "r"
57         case "c"
58             orient = "c"
59         else
60             msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
61             error(msprintf(msg, "unique", i, "1,2,""r"",""c"",""keepOrder"",""uniqueNan"""))
62         end
63         in(1) = null()
64         i = 3
65     end
66     if or(orient==["r" "c"]) & ndims(x)>2 then
67         msg = _("%s: Argument #%d: ''%s'' not allowed for an hypermatrix.\n")
68         error(msprintf(msg, "unique", 2, orient))
69     end
70
71     while size(in)>0 & i<5 then
72         a = in(1)
73         if typeof(a)=="string"
74             a = convstr(a)
75         end
76         select a
77         case "uniquenan"
78             uniqueNan = %t
79         case "keeporder"
80             keepOrder = %t
81         else
82             msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
83             error(msprintf(msg, "unique", i, """keepOrder"",""uniqueNan"""))
84         end
85         in(1) = null()
86         i = i+1
87     end
88     uniqueNan = uniqueNan & or(type(x)==[1 5])
89
90     sz = size(x);
91     if size(x, orient)==1 then
92         ki = 1
93         ko = 1
94         nb = 1
95         return
96     end
97     if uniqueNan
98         [x, newInf] = uniqueProcessNan(x, [], "removeNan")
99     end
100     getK = argn(1) > 1 | keepOrder
101
102
103     // [] trivial case
104     // ---------------
105     if isempty(x) then
106         return  // ki, ko, nb are already []. x is [] or sparse([])
107     end
108
109     // PROCESSING
110     // ----------
111     areComplex = or(type(x)==[1 5]) && ~isreal(x, 0)
112     if orient=="*" then
113         if getK then
114             // ki (begin)
115             if areComplex
116                 [x,ki] = gsort(x,"g",["i" "i"], list(abs, atan));
117             else
118                 [x,ki] = gsort(x,"g","i");
119             end
120             keq = x(2:$) == x(1:$-1);
121             if argn(1) > 2
122                 // ko (begin)
123                 ko(ki) = cumsum(~[%F ; keq(:)])
124                 if ko <> [] then
125                     ko = matrix(ko, sx)
126                 end
127                 // nb
128                 if argn(1) > 3
129                     nb = [0 find(~keq) size(x,"*")]
130                     nb = nb(2:$) - nb(1:$-1)
131                 end
132             end
133             // ki (end)
134             keq = find(keq);
135             if keq <> [] then keq = keq + 1; end
136             x(keq) = [];
137             ki(keq) = [];
138         else
139             if areComplex
140                 x = gsort(x,"g",["d" "d"], list(abs, atan));
141             else
142                 x = gsort(x,"g","d");
143             end
144             x = x($:-1:1);
145             x( find(x(2:$) == x(1:$-1)) ) = [];
146         end
147         if sx(1) > 1 | length(sx) > 2
148             x = x(:)
149             ki = ki(:)
150             nb = nb(:)
151         end
152
153     else // orient used
154         if  orient==2 | orient=="c" then
155             x = x.'
156         end
157         if getK then
158             if areComplex
159                 [x,ki] = gsort(x,"lr",["i" "i"], list(abs, atan));
160             else
161                 [x,ki] = gsort(x,"lr","i");
162             end
163             keq = and(x(2:$,:) == x(1:$-1,:),"c")
164             if argn(1)>2
165                 ko(ki) = cumsum(~[%F ; keq(:)])
166                 // nb
167                 if argn(1) > 3
168                     nb = [0 find(~keq) size(x,1)]
169                     nb = nb(2:$) - nb(1:$-1)
170                     nb = nb(:)
171                 end
172             end
173             keq = find(keq)
174             if keq <> [] then keq = keq + 1;end
175             x(keq,:) = [];
176             ki(keq,:) = [];
177         else
178             if areComplex
179                 x = gsort(x,"lr",["i" "i"], list(abs, atan));
180             else
181                 x = gsort(x,"lr","i");
182             end
183             x( find(and(x(2:$,:) == x(1:$-1,:),"c")),:) = [];
184         end
185         if  orient==2 | orient=="c" then
186             x = x.'
187             ki = matrix(ki, 1, -1)
188             ko = ko'
189             nb = nb'
190         end
191     end
192
193     if uniqueNan
194         x = uniqueProcessNan(x, newInf, "restoreNan")
195     end
196
197     if keepOrder
198         [ki, kk] = gsort(ki,"g","i")
199         select orient
200         case "*"
201             x = x(kk)
202         case "r"
203             x = x(kk,:)
204         case "c"
205             x = x(:,kk)
206         end
207         if argn(1)>2
208             nb = nb(kk)
209             [?, kk2] = gsort(kk,"g","i")
210             ko = kk2(ko)
211         end
212     end
213 endfunction
214
215 // -------------------------------------------------------------------
216
217 // To consider Nan mutually equal, we replace all of them with a "regular" substitute.
218 // Since Nan are sorted as > Inf, we must use anyway Inf as the Nan substitute.
219 // If the original array have already some Inf, we must priorly replace them with
220 // a decimal greater than the finite maximum of the array values.
221 // After processing, we restore Inf => Nan (, and maxNum => Inf).
222
223 function [x, newInf] = uniqueProcessNan(x, newInf, way)
224
225     if way=="removeNan" & or(isnan(x)) then
226         // Replacing Nan
227         // -------------
228         if isreal(x)
229             if or(x==%inf)
230                 b = x==%inf
231                 m = max([1 max(x(~b))])
232                 newInf = m*1.5
233                 x(b) = newInf
234             end
235             x(x<>x) = %inf
236         else
237             r = real(x)
238             i = imag(x)
239             if or([r i]==%inf)
240                 br = r==%inf
241                 m = max(r(~br),1)
242                 bi = i==%inf
243                 m = max(i(~bi),m,1)
244                 newInf = m*1.5
245                 r(br) = newInf
246                 i(bi) = newInf
247             end
248             r(r<>r) = %inf
249             i(i<>i) = %inf
250             x = complex(r,i);
251         end
252
253     // Restoring  NaN
254     // --------------
255     elseif way=="restoreNan"
256         if isreal(x)
257             x(x==%inf) = %nan
258             if newInf~=[]
259                 x(x==newInf) = %inf
260             end
261         else
262             r = real(x)
263             r(r==%inf) = %nan
264             i = imag(x)
265             i(i==%inf) = %nan
266             if newInf~=[]
267                 r(r==newInf) = %inf
268                 i(i==newInf) = %inf
269             end
270             x = complex(r, i)
271         end
272     end
273 endfunction