b691b496b00a7a9a9f013c85ef57faba1bcf6016
[scilab.git] / scilab / modules / elementary_functions / macros / members.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2013 - Samuel GOUGEON
3 // Copyright (C) 2009 - Université du Maine - Samuel GOUGEON
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING, which
7 // you should have received as part of this distribution.  The terms
8 // are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
10
11 function [nb, loc] = members(A, S, varargin)
12     //
13     // Looks for how many times in S each element of A occurs
14     // Looks for how many rows of S match each row of A
15     // Can return the position in S of each respective first or last match.
16     //
17     // Calling sequence:
18     // -------------------------
19     // [nb [,loc]] = members(A, S)
20     // [nb [,loc]] = members(A, S, "last")
21     // [nb [,loc]] = members(A, S, "rows"|"cols")
22     // [nb [,loc]] = members(A, S, "rows"|"cols", "last")
23     // [nb [,loc]] = members(A, S, "rows"|"cols", "shuffle")
24     // [nb [,loc]] = members(A, S, "rows"|"cols", "shuffle", "last")
25     //
26     // Input / Output arguments:
27     // -------------------------
28     // A   : Matrix or hypermatrix of booleans, integers, reals, complexes, polynomials or strings:
29     //       entities (needles) to search in S (haystack).
30     // S   : Matrix or hypermatrix of same datatype as A
31     // "last": litteral keyword (optional): if it is provided, indices of last
32     //       occurrences are returned instead of for the first ones.
33     // "rows": litteral keyword (optional). if it is provided, each row of A
34     //       is searched as a whole, among rows of S.
35     // "cols": litteral keyword (optional). if it is provided, each column of A
36     //       is searched as a whole, among columns of S.
37     //       "rows" and "cols" options are exclusive. If both are specified,
38     //       only the last specified one is considered.
39     // "shuffle": litteral keyword (optional), modifying the "rows" or "cols"
40     //       option. If it is provided, detection of -- say -- A rows in S is
41     //       performed without respect of the order of elements in each row.
42     //       This option is ignored when processing polynomials.
43     //
44     // NORMAL mode: neither "rows" nor "cols" is set. Then,
45     // nb  : Matrix of reals: same sizes as A
46     //       nb(i, j, ...): number of occurrences of A(i, j, ...) in S.
47     // loc : Matrix of reals: same sizes as A
48     //       loc(i, j, ...): linear index in S of the first occurrence of A(i, j, ...).
49     //       If "last" is set, the index of the last occurrence is returned instead.
50     //       loc(i, j, ...) returns zero if A(i, j, ...) is not found.
51     // ROW-WISE mode:
52     // nb  : Row of reals with size(A, 1) elements
53     //       nb(i): number of occurrences of A(i, :) found as rows of S
54     // loc : Row of reals with size(A, 1) elements
55     //       loc(i): index of the first row in S which matches A(i, :).
56     //       If "last" is set, the index of the last occurrence is returned instead.
57     //       loc(i, j, ...) returns zero if A(i, :) is not found.
58     // COLUMN-WISE mode:
59     //  same as above. nb and loc are row vectors with size(A, 2) elements.
60     //
61     // REMARK: %inf, -%inf values are supported in A as well as in S.
62     //
63     // LIMITATION: in normal element-wise mode, %nan are supported only in A.
64     //
65     // Examples:
66     // ---------
67     // a) with reals:
68     //   N = [ 7  3
69     //       %inf 0
70     //       %nan 1 ];
71     //   H = [ 5   8    0   4
72     //         3   4    7   7
73     //         3 %inf %inf  2
74     //         7   5    5   8 ];
75     //   [nb, loc] = members(N, H)
76     //   [nb, loc] = members(N, H, "last")
77     //
78     // b) with hypermatrices, from previous N and H:
79     //   N = matrix(N, [3 1 2]);
80     //   H = matrix(H, [4 2 2]);
81     //   [nb, loc] = members(N, H, "last")
82     //
83     // c) with integers:
84     //   N = int8(grand(3, 2, "uin", -5, 5));
85     //   H = int8(grand(4, 4, "uin", -5, 5));
86     //   [nb, loc] = members(N, H)
87     //
88     // d) with polynomials (complex coefficients are accepted):
89     //   z = %z;
90     //   N = [z (1-z)^2 ; -4 %i*z ];
91     //   H = [2  %i*z -z  3-z  z  z^3 z];
92     //   [nb, loc] = members(N, H)
93     //
94     // e) with text:
95     //   N = [ "Hi" "Hu" "Allo"];
96     //   H = [ "Hello" "Bonjour" "Allo"
97     //         "Holà"  "Allo"  "Hallo"
98     //         "Hi"    "Hé"    "Salud" ];
99     //   [nb, loc] = members(N, H, "last")
100     //
101     // f) by rows:
102     //   H = [
103     //    3  3  0
104     //    4  1  0
105     //    2  0  3
106     //    0  1  4
107     //    3  4  3
108     //    0  1  4
109     //    3  1  0 ];
110     //   N = [
111     //    1 2 3
112     //    0 1 4
113     //    3 0 3
114     //    4 1 0
115     //    2 0 2 ];
116     //   N, H
117     //   [nb, loc] = members(N, H, "rows")
118     //   [nb, loc] = members(N, H, "rows", "last")
119     //   [nb, loc] = members(N, H, "rows", "shuffle")
120     //
121     // g) by columns, from previous N and H:
122     //   N = N.', H = H.'
123     //   [nb, loc] = members(N, H, "cols", "shuffle")
124
125     [lhs, rhs] = argn();
126     nb = [];
127     if rhs == 0 then
128         head_comments("members");
129         return
130     end
131     if rhs < 2 then
132         error(msprintf(gettext("%s: Wrong number of input argument(s): at least %d expected.\n"), "members", 2));
133     end
134     if A == [] then
135         if lhs > 1 then
136             loc = [];
137         end
138         return
139     end
140     if  S == []  then
141         nb = zeros(A);
142         if lhs > 1 then
143             loc = zeros(A);
144         end
145         return
146     end
147
148     sA = size(A);
149     type_A = type(A(:));
150     type_S = type(S(:));
151     if type_A ~= type_S then
152         error(msprintf(gettext("%s: Wrong type for input argument #%d: expected same type as first argument.\n"), "members", 2));
153     end
154     if and(type_A ~= [1 2 4 8 10]) then
155         error(msprintf(gettext("%s: Wrong type for input argument #%d: Matrix of encoded-Integers, Reals, Complexes, Booleans, Polynomials or Strings expected.\n"), "members", 1));
156     end
157     if and(type_S ~= [1 2 4 8 10]) then
158         error(msprintf(gettext("%s: Wrong type for input argument #%d: Matrix of encoded-Integers, Reals, Complexes, Booleans, Polynomials or Strings expected.\n"), "members", 2));
159     end
160     if or(isnan(S)) then
161         error(msprintf(gettext("%s: Wrong value for argument #%d: Must not contain NaN.\n"), "members", 2));
162     end
163
164     // Checking optional flags (if any):
165     last = %f;
166     direction = "";
167     shuffle = %f;
168     shuffle_temp = %f;
169     if  rhs > 2 then
170         i = 2;
171         for option=varargin
172             i = i + 1;
173             if typeof(option) ~= "string" | ~isscalar(option) then
174                 msg = _("%s: Wrong type for argument #%d: Scalar string expected => option ignored.\n");
175                 warning(msprintf(msg, "members", i))
176                 continue
177             end
178             o = convstr(option, "l");
179             if o == "last" then
180                 last = %t;
181             elseif o == "rows"
182                 direction = "rows";
183             elseif o == "cols"
184                 direction = "cols";
185             elseif o == "shuffle"
186                 shuffle_temp = %t;
187             else
188                 msg = _("%s: Unknown option ""%s"" => ignored.\n");
189                 warning(msprintf(msg, "members", o))
190             end
191         end
192         if shuffle_temp == %t then
193             if direction == "" then
194                 msg = _("%s: Option ""%s"" only relevant for ""%s"" and ""%s"" modes => ignored.\n");
195                 warning(msprintf(msg, "members", "shuffle", "rows", "cols"))
196             else
197                 shuffle = %t;
198             end
199         end
200     end
201
202     // ------------------------------------------------------------------------
203
204     // Usual processing: searching for matching individual components
205     // --------------------------------------------------------------
206     if direction == "" then
207         // Reals: special faster processing with dsearch().
208         //  %inf, -%inf are supported in A and S. %nan are only supported in A.
209         //  Processing integers and strings requires bugs 6305 & 12778 to be fixed.
210
211         if type_A == 8 then   // Convert integers into reals in order to use dsearch
212             A = double(A);
213             S = double(S);
214         elseif type_A == 4 then   // Convert booleans into reals
215             A = bool2s(A);
216             S = bool2s(S);
217         end
218         if type_A == 1 & isreal(A) & isreal(S) then
219             S = S(:);
220             if last then
221                 S = S($:-1:1);
222             end
223             [Su, kS] = unique(S);
224             [i, nbS] = dsearch(S, Su, "d");
225
226             I = dsearch(A(:), Su, "d");
227             k = find(I~=0);
228             nb = I;
229             nb(k) = nbS(I(k));
230             nb = matrix(nb, size(A));
231             if lhs > 1 then
232                 loc = I;
233                 loc(k) = kS(I(k));
234                 if last then
235                     loc(k) = length(S)-loc(k)+1;
236                 end
237                 loc = matrix(loc, size(A));
238             end
239             // ------------------------------------------------------------------------
240         else
241             // Other cases : polynomials, text, complexes
242             // ==========================================
243             LA = size(A, "*");
244             LS = size(S, "*");
245             A = A(:);
246             if ~last then
247                 S = S($:-1:1);
248             end
249
250             // Slices may be needed in case of memory overflow
251             //
252             // Function returning the memory space in [bytes] occupied by each
253             //  querried variable. If a variable exists in several contexts,
254             //  only the last created one is considered
255             function varargout = varspace(varargin)
256                 [Ln, Ls] = who("local");
257                 vs = list();
258                 for vn = varargin
259                     i = find(Ln==vn);
260                     if i ~= [] then
261                         vs($+1) = Ls(i(1))*8;
262                     else
263                         vs($+1) = 0;
264                     end
265                 end
266                 varargout = vs;
267             endfunction
268             // -----------------------------------
269             // Memory space needed to process a i:j slice of A:
270             //   nb  : j doubles (occupied at the end)
271             //   loc : j doubles (occupied at the end)
272             //   tmp : LS * (j-i+1) doubles needed
273             //   A2  : LS * mem(A(i:j)) needed
274             //   S2  : mem(S)*(j-i+1) needed
275
276             // Function setting the thickness of the slice, according to the
277             //    available and needed memory
278             function j = set_j(i, LS, memS, rA)
279                 // rA : remaining unprocessed part of A
280                 [memA2, memS2, memtmp] = varspace("A2", "S2", "tmp");
281                 // Available memory [bytes]:
282                 m = stacksize();
283                 avail = (m(1) - m(2))*8 + memA2 + memS2 + memtmp - 16*(i-1); // [bytes]
284                 avail = 0.15*avail;  // A lot of memory is used as Intermediate memory
285                 // Diadic loop on e = j-i+1 (slice's thickness):
286                 // Init
287                 e = size(rA, "*");
288                 Vtest = rA(1:e);
289                 memrA = varspace("Vtest");   // [bytes]
290                 // Loop fitting the slice's thickness e for the available memory
291                 while e*(memS+8*LS+16)+LS*memrA > avail & e > 1
292                     e = ceil(e/2);   // Remaining thickness => divided by 2
293                     Vtest = rA(1:e);
294                     memrA = varspace("Vtest");
295                 end
296                 if (e*(memS+8*LS+16)+LS*memrA > avail) then
297                     msg = _("%s: not enough memory to proceed\n");
298                     error(sprintf(msg, "members"));
299                 end
300                 j = i+e-1;
301             endfunction
302
303             memS = varspace("S");  // [bytes]
304
305             // Starting and ending indices in A for the current slice
306             i = 1;   // Starting index: initialization
307             j = 0;   // Ending index: initialization
308             while j<LA      // Slicing loop
309                 // Setting next j (=> thickness of the slice = j-i+1)
310                 j = set_j(i, LS, memS, A(i:$))
311
312                 // Progression bar
313                 if j < LA then
314                     if isdef("waitH") then
315                         waitbar(j/LA, waitH);
316                     elseif j/LA < 0.4   // Else: too few steps => not worthwhile
317                         waitH = waitbar(j/LA);
318                     end
319                 end
320
321                 A2 = repmat(A(i:j), 1, LS);
322                 S2 = repmat(S(:).', j-i+1, 1);
323
324                 tmp = bool2s(A2==S2);
325                 nb(i:j) = sum(tmp, "c");
326                 if lhs > 1 then
327                     tmp = tmp.* ( (1:LS) .*. ones(j-i+1, 1) );
328                     loc(i:j) = max(tmp, "c");
329                 end
330                 i = j+1;
331                 // End of slice processing
332             end
333             // Deleting the waitbar (if any) and clearing transient variables
334             clear A2 S2 tmp
335             if isdef("waitH") then
336                 delete(waitH);
337             end
338
339             // Final operations on the overall result(s)
340             nb = matrix(nb, sA);
341             if lhs > 1
342                 if ~last
343                     k = loc~=0;
344                     loc(k) = LS - loc(k) + 1;
345                 end
346                 loc = matrix(loc, sA);
347             end
348         end
349         // ========================================================================
350
351     else
352         // Row-wise processing: searching for matching rows
353         // ------------------------------------------------
354
355         // Additionnal input checking:
356         A = squeeze(A);
357         if ~ismatrix(A) then
358             msg = _("%s: Wrong type for argument #%d: Matrix expected.\n"); // error #209
359             error(msprintf(msg, "members", 1))
360         end
361
362         S = squeeze(S);
363         if ~ismatrix(S) then
364             msg = _("%s: Wrong type for argument #%d: Matrix expected.\n"); // error #209
365             error(msprintf(msg, "members", 2))
366         end
367
368         if direction == "rows" & size(A, 2) ~= size(S, 2) then
369             msg = _("%s: Incompatible input arguments #%d and #%d: Same number of columns expected.\n");
370             error(msprintf(msg, "members", 1, 2))
371         elseif direction == "cols" & size(A, 1) ~= size(S, 1) then
372             msg = _("%s: Incompatible input arguments #%d and #%d: Same number of rows expected.\n");
373             error(msprintf(msg, "members", 1, 2))
374         end
375
376         // Column-wise = Row-wise after transposition
377         if direction == "cols" then
378             A = A.';
379             S = S.';
380         end
381
382         // If "shuffle" is set, we first sort elements along each row of A and S
383         if  type_A ~= 2 & shuffle then    // Polynomials are not sortable
384             A = gsort(A, "c", "i");
385             S = gsort(S, "c", "i");
386         end
387
388         // Pre-Processing: Preparing initial chart of double indices
389         nLA = size(A, 1);  // Number of rows of A
390         nLS = size(S, 1);  // Number of rows of S
391         I1 = (1:nLA)' .*. ones(nLS, 1); // indices of A rows, grouped by S rows
392         I2 = ones(nLA, 1) .*. (1:nLS)';  // set of S rows indices, replicated by the number of A rows
393         IND = [I1 I2];   // Column #1 gives index in A. column #2 gives index in S
394
395         // Processing: loop over columns
396         for k = 1:size(S, 2)
397             cA = A(IND(:, 1), k);
398             cS = S(IND(:, 2), k);
399             IND = IND(find(cA==cS), :);
400             if IND == [] then
401                 break
402             end
403         end
404
405         // Post-processing
406         if last then
407             [nb, loc] = members((1:nLA)', IND(:, 1), "last");
408         else
409             [nb, loc] = members((1:nLA)', IND(:, 1));
410         end
411         k = find(loc~=0);
412         loc(k) = IND(loc(k), 2);
413         nb = nb';
414         loc = loc';
415     end
416
417 endfunction