db2d19418045269cc5628df11b0d761eb0be0c50
[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, 2019 - 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
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
23     keepOrder = %f
24     uniqueNan = %f
25     orient = "*"
26     newInf = [] // init Inf substitute in case of "uniqueNan" and or(x==%inf)
27     ki = []
28     ko = []
29     nb = []
30
31     // CHECKING INPUT ARGUMENTS
32     // ------------------------
33     in = varargin
34     i = 2;  // index of the current input argument
35     if size(in)>0 then
36         a = in(1)
37         if typeof(a)=="string"
38             a = convstr(a)
39         end
40         select a
41         case 1
42             orient = "r"
43         case 2
44             orient = "c"
45         case "uniquenan"
46             uniqueNan = %t
47         case "keeporder"
48             keepOrder = %t
49         case "r"
50             orient = "r"
51         case "c"
52             orient = "c"
53         else
54             msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
55             error(msprintf(msg, "unique", i, "1,2,""r"",""c"",""keepOrder"",""uniqueNan"""))
56         end
57         in(1) = null()
58         i = 3
59     end
60     while size(in)>0 & i<5 then
61         a = in(1)
62         if typeof(a)=="string"
63             a = convstr(a)
64         end
65         select a
66         case "uniquenan"
67             uniqueNan = %t
68         case "keeporder"
69             keepOrder = %t
70         else
71             msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
72             error(msprintf(msg, "unique", i, """keepOrder"",""uniqueNan"""))
73         end
74         in(1) = null()
75         i = i+1
76     end
77     uniqueNan = uniqueNan & or(type(x)==[1 5])
78
79     sz = size(x);
80     if size(x, orient)==1 then
81         ki = 1
82         return
83     end
84     if uniqueNan
85         [x, newInf] = uniqueProcessNan(x, [], "removeNan")
86     end
87     getK = argn(1)>1 | keepOrder
88
89
90     // [] trivial case
91     // ---------------
92     if isempty(x) then
93         return  // ki, nb are already []. x is [] or sparse([])
94     end
95
96     // PROCESSING complex numbers
97     // --------------------------
98     if or(type(x)==[1 5]) then
99         if ~isreal(x)
100             if isreal(x,0)
101                 x = real(x);
102             else
103                 if orient=="*"
104                     x = [real(x(:)) imag(x(:))]
105                     if ~getK
106                         x = unique(x,"r")
107                     else
108                         [x, ki, ko, nb] = unique(x,"r")
109                     end
110                     x = complex(x(:,1),x(:,2));
111                     if sz(1)==1 // => put results in row
112                         x = x.'
113                         if getK
114                             ki = ki'
115                             nb = nb'
116                         end
117                     end
118                 elseif orient=="r" | orient==1
119                     x = [real(x) imag(x)]
120                     if ~getK
121                         x = unique(x,"r")
122                     else
123                         [x, ki, ko, nb] = unique(x,"r")
124                     end
125                     x = complex(x(:,1:sz(2)), x(:,sz(2)+1:$));
126                 elseif orient=="c" | orient==2
127                     x = [real(x) ; imag(x)]
128                     if ~getK
129                         x = unique(x,"c")
130                     else
131                         [x, ki, ko, nb] = unique(x,"c")
132                     end
133                     x = complex(x(1:sz(1),:), x(sz(1)+1:$,:));
134                 end
135                 if uniqueNan
136                     x = uniqueProcessNan(x, newInf, "restoreNan")
137                 end
138                 if keepOrder
139                     [ki, kk] = gsort(ki,"g","i")
140                     select orient
141                     case "*"
142                         x = x(kk)
143                     case "r"
144                         x = x(kk,:)
145                     case "c"
146                         x = x(:,kk)
147                     end
148                     nb = nb(kk)
149                 end
150                 return
151             end
152         end
153     end
154
155     // PROCESSING text and other numerical types
156     // -----------------------------------------
157     if orient=="*" then
158         if getK then
159             [x,ki] = gsort(x,"g","i");
160             keq = x(2:$) == x(1:$-1);
161             if argn(1)>2
162                 nb = [0 find(~keq) size(x,"*")]
163                 nb = nb(2:$) - nb(1:$-1)
164             end
165             keq = find(keq);
166             if keq<>[] then keq = keq+1;end
167             x(keq) = [];
168             ki(keq) = [];
169             if size(x,1)>1 | ndims(x)>2
170                 x = x(:)
171                 ki = ki(:)
172                 nb = nb(:)
173             end
174         else
175             x = gsort(x,"g","d");
176             x = x($:-1:1);
177             x( find(x(2:$) == x(1:$-1)) ) = [];
178         end
179     elseif  orient==1|orient=="r" then
180         if getK then
181             [x,ki] = gsort(x,"lr","i");
182             keq = and(x(2:$,:) == x(1:$-1,:),"c")
183             if argn(1)>2
184                 nb = [0 find(~keq) size(x,1)]
185                 nb = nb(2:$) - nb(1:$-1)
186                 nb = nb(:)
187             end
188             keq = find(keq)
189             if keq<>[] then keq = keq+1;end
190             x(keq,:) = [];
191             ki(keq,:) = [];
192         else
193             x = gsort(x,"lr","i");
194             x( find(and(x(2:$,:) == x(1:$-1,:),"c")),:) = [];
195         end
196     elseif  orient==2|orient=="c" then
197         if getK then
198             [x,ki] = gsort(x,"lc","i");
199             keq = and(x(:,2:$) == x(:,1:$-1),"r")
200             if argn(1)>2
201                 nb = [0 find(~keq) size(x,2)]
202                 nb = nb(2:$) - nb(1:$-1)
203             end
204             keq = find(keq)
205             if keq<>[] then keq = keq+1;end
206             x(:,keq) = [];
207             ki(:,keq) = [];
208         else
209             x = gsort(x,"lc","i");
210             x(:, find(and(x(:,2:$) == x(:,1:$-1),"r")) ) = [];
211         end
212     end
213     if uniqueNan
214         x = uniqueProcessNan(x, newInf, "restoreNan")
215     end
216     if keepOrder
217         [ki, kk] = gsort(ki,"g","i")
218         select orient
219         case "*"
220             x = x(kk)
221         case "r"
222             x = x(kk,:)
223         case "c"
224             x = x(:,kk)
225         end
226         if argn(1)>2
227             nb = nb(kk)
228         end
229     end
230 endfunction
231 // -------------------------------------------------------------------
232
233 // To consider Nan mutually equal, we replace all of them with a "regular" substitute.
234 // Since Nan are sorted as > Inf, we must use anyway Inf as the Nan substitute.
235 // If the original array have already some Inf, we must priorly replace them with
236 // a decimal greater than the finite maximum of the array values.
237 // After processing, we restore Inf => Nan (, and maxNum => Inf).
238
239 function [x, newInf] = uniqueProcessNan(x, newInf, way)
240
241     if way=="removeNan" & or(isnan(x)) then 
242         // Replacing Nan
243         // -------------
244         if isreal(x)
245             if or(x==%inf)
246                 b = x==%inf
247                 m = max([1 max(x(~b))])
248                 newInf = m*1.5
249                 x(b) = newInf
250             end
251             x(x<>x) = %inf
252         else
253             r = real(x)
254             i = imag(x)
255             if or([r i]==%inf)
256                 br = r==%inf
257                 m = max(r(~br),1)
258                 bi = i==%inf
259                 m = max(i(~bi),m,1)
260                 newInf = m*1.5
261                 r(br) = newInf
262                 i(bi) = newInf
263             end
264             r(r<>r) = %inf
265             i(i<>i) = %inf
266             x = complex(r,i);
267         end
268
269     // Restoring  NaN
270     // --------------
271     elseif way=="restoreNan"
272         if isreal(x)
273             x(x==%inf) = %nan
274             if newInf~=[]
275                 x(x==newInf) = %inf
276             end
277         else
278             r = real(x)
279             r(r==%inf) = %nan
280             i = imag(x)
281             i(i==%inf) = %nan
282             if newInf~=[]
283                 r(r==newInf) = %inf
284                 i(i==newInf) = %inf
285             end
286             x = complex(r, i)
287         end
288     end
289 endfunction