* Bug #12705 fixed - Elementary_functions & m2sci: members() function
[scilab.git] / scilab / modules / elementary_functions / macros / members.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2009,2013 - Université du Maine - Samuel Gougeon
3 //
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10 function [nb, loc] = members(A, S, last)
11     //
12     // Input / Output arguments:
13     // -------------------------
14     // A   : Matrix or hypermatrix of booleans, integers, reals, complexes, polynomials or strings:
15     //       entities (needles) to search in S (haystack).
16     // S   : Matrix or hypermatrix of same datatype as A
17     // last: Scalar boolean
18     // nb  : Matrix of reals: same sizes as A
19     //       nb(i, j, ...): number of occurrences of A(i, j, ...) in S.
20     // loc : Matrix of reals: same sizes as A
21     //       loc(i, j, ...): linear index in S of the first occurrence of A(i, j, ...).
22     //       If last==%t, the index of the last occurrence is returned instead.
23     //       loc(i, j, ...) returns zero if A(i, j, ...) is not found.
24     //
25     // %inf, -%inf values are supported in A as well as in S.
26     // %nan are supported only in A.
27     //
28     // Examples:
29     // ---------
30     // a) with reals:
31     // N = [ 7  3
32     //     %inf 0
33     //     %nan 1 ];
34     // H = [ 5   8    0   4
35     //       3   4    7   7
36     //       3 %inf %inf  2
37     //       7   5    5   8 ];
38     // [nb, loc] = members(N, H)
39     // [nb, loc] = members(N, H, %T)
40     //
41     // b) with hypermatrices, from previous N and H:
42     // N = matrix(N, [3 1 2]);
43     // H = matrix(H, [4 2 2]);
44     // [nb, loc] = members(N, H, %T)
45     //
46     // c) with integers:
47     // N = int8(grand(3, 2, "uin", -5, 5));
48     // H = int8(grand(4, 4, "uin", -5, 5));
49     // [nb, loc] = members(N, H)
50     //
51     // d) with polynomials (complex coefficients are accepted):
52     // z = %z;
53     // N = [z (1-z)^2 ; -4 %i*z ];
54     // H = [2  %i*z -z  3-z  z  z^3 z];
55     // [nb, loc] = members(N, H)
56     //
57     // e) with text:
58     // N = [ "Hi" "Hu" "Allo"];
59     // H = [ "Hello" "Bonjour" "Allo"
60     //       "Holà"  "Allo"  "Hallo"
61     //       "Hi"    "Hé"    "Salud" ];
62     // [nb, loc] = members(N, H, %t)
63
64     [lhs, rhs] = argn();
65     if rhs < 2 then
66         error(msprintf(gettext("%s: Wrong number of input argument(s): at least %d expected.\n"), "members", 2));
67     end
68     if rhs == 2 then
69         last = %f;
70     end
71     if A == [] then
72         nb = [];
73         if lhs > 1 then
74             loc = [];
75         end
76         return
77     end
78     if  S == []  then
79         nb = zeros(A);
80         if lhs > 1 then
81             loc = zeros(A);
82         end
83         return
84     end
85
86     type_A = type(A(:));
87     type_S = type(S(:));
88     if type_A <> type_S then
89         error(msprintf(gettext("%s: Wrong type for input argument #%d: expected same type as first argument.\n"), "members", 2));
90     end
91     if and(type_A <> [1 2 4 8 10]) then
92         error(msprintf(gettext("%s: Wrong type for input argument #%d: Matrix of Integers, Reals, Complexes, Booleans, Polynomials or Strings expected.\n"), "members", 1));
93     end
94     if and(type_S <> [1 2 4 8 10]) then
95         error(msprintf(gettext("%s: Wrong type for input argument #%d: Matrix of Integers, Reals, Complexes, Booleans, Polynomials or Strings expected.\n"), "members", 2));
96     end
97     if or(isnan(S)) then
98         error(msprintf(gettext("%s: Wrong value for argument #%d: Must not contain NaN.\n"), "members", 2));
99     end
100     if type(last) <> 4 then
101         error(msprintf(gettext("%s: Wrong type for input argument #%d: Boolean matrix expected.\n"), "members", 3));
102     end
103     if ~isscalar(last) then
104         error(msprintf(gettext("%s: Wrong size for input argument #%d: Scalar expected.\n"), "members", 3));
105     end
106
107     // ------------------------------------------------------------------------
108
109     // Reals: special faster processing with dsearch().
110     //  %inf, -%inf are supported in A and S. %nan are only supported in A.
111     //  Processing integers and strings requires bugs 6305 & 12778 to be fixed.
112
113     if type_A == 8 then   // Convert integers into reals in order to use dsearch
114         A = double(A);
115         S = double(S);
116     end
117     if type_A == 1 & isreal(A) & isreal(S) then
118         S = S(:);
119         if last then
120             S = S($:-1:1);
121         end
122         [Su, kS] = unique(S);
123         [i, nbS] = dsearch(S, Su, "d");
124
125         I = dsearch(A(:), Su, "d");
126         k = find(I~=0);
127         nb = I;
128         nb(k) = nbS(I(k));
129         nb = matrix(nb, size(A));
130         if lhs > 1 then
131             loc = I;
132             loc(k) = kS(I(k));
133             if last then
134                 loc(k) = length(S)-loc(k)+1;
135             end
136             loc = matrix(loc, size(A));
137         end
138         // ------------------------------------------------------------------------
139     else
140         // Other cases : polynomials, text, complexes
141         // ==========================================
142         sA = size(A);
143         LA = size(A, "*");
144         LS = size(S, "*");
145         A = A(:);
146         if ~last then
147             S = S($:-1:1);
148         end
149
150         // Slices may be needed in case of memory overflow
151         //
152         // Function returning the memory space in [bytes] occupied by each
153         //  querried variable. If a variable exists in several contexts,
154         //  only the last created one is considered
155         function varargout = varspace(varargin)
156             [Ln, Ls] = who("local");
157             vs = list();
158             for vn = varargin
159                 i = find(Ln==vn);
160                 if i ~= [] then
161                     vs($+1) = Ls(i(1))*8;
162                 else
163                     vs($+1) = 0;
164                 end
165             end
166             varargout = vs;
167         endfunction
168         // -----------------------------------
169         // Memory space needed to process a i:j slice of A:
170         // nb  : j doubles (occupied at the end)
171         // loc : j doubles (occupied at the end)
172         // tmp : LS * (j-i+1) doubles needed
173         // A2  : LS * mem(A(i:j)) needed
174         // S2  : mem(S)*(j-i+1) needed
175
176         // Function setting the thickness of the slice, according to the
177         //    available and needed memory
178         function j = set_j(i, LS, memS, rA)
179             // rA : remaining unprocessed part of A
180             [memA2, memS2, memtmp] = varspace("A2", "S2", "tmp");
181             // Available memory [bytes]:
182             m = stacksize();
183             avail = (m(1)-m(2))*8 + memA2 + memS2 + memtmp -16*(i-1); // [bytes]
184             avail = 0.15*avail;  // A lot of memory is used as Intermediate memory
185             // Diadic loop on e = j-i+1 (slice's thickness):
186             // Init
187             e = size(rA, "*");
188             Vtest = rA(1:e);
189             memrA = varspace("Vtest");   // [bytes]
190             // Loop fitting the slice's thickness e for the available memory
191             while (e*(memS+8*LS+16)+LS*memrA > avail & e > 1)
192                 e = ceil(e/2);   // Remaining thickness => divided by 2
193                 Vtest = rA(1:e);
194                 memrA = varspace("Vtest");
195             end
196             if (e*(memS+8*LS+16)+LS*memrA > avail) then
197                 msg = _("%s: not enough memory to proceed");
198                 error(sprintf(msg, "members"));
199             end
200             j = i+e-1;
201         endfunction
202
203         memS = varspace("S");  // [bytes]
204
205         // Starting and ending indices in A for the current slice
206         i = 1;   // Starting index: initialization
207         j = 0;   // Ending index: initialization
208         while j < LA      // Slicing loop
209             // Setting next j (=> thickness of the slice = j-i+1)
210             j = set_j(i, LS, memS, A(i:$));
211
212             // Progression bar
213             if j < LA then
214                 if isdef("waitH") then
215                     waitbar(j/LA, waitH);
216                 elseif j/LA < 0.4   // Else: too few steps => not worthwhile
217                     waitH = waitbar(j/LA);
218                 end
219             end
220
221             A2 = repmat(A(i:j), 1, LS);
222             S2 = repmat(S(:).', j-i+1, 1);
223
224             tmp = bool2s(A2==S2);
225             nb(i:j) = sum(tmp, "c");
226             if lhs > 1 then
227                 tmp = tmp.* ( (1:LS) .*. ones(j-i+1, 1) );
228                 loc(i:j) = max(tmp, "c");
229             end
230             i = j+1;
231             // End of slice processing
232         end
233         // Deleting the waitbar (if any) and clearing transient variables
234         clear A2 S2 tmp
235         if isdef("waitH") then
236             delete(waitH);
237         end
238
239         // Final operations on the overall result(s)
240         nb = matrix(nb, sA);
241         if lhs > 1
242             if ~last
243                 k = loc~=0;
244                 loc(k) = LS - loc(k) + 1;
245             end
246             loc = matrix(loc, sA);
247         end
248     end
249
250 endfunction