vectorfind() upgrade: short needle, wildcard, hypermats, %nan
[scilab.git] / scilab / modules / elementary_functions / macros / vectorfind.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
3 // Copyright (C) 2010 - DIGITEO - Vincent COUVERT
4 // Copyright (C) 2017 - Samuel GOUGEON
5 //
6 // Copyright (C) 2012 - 2016 - Scilab Enterprises
7 //
8 // This file is hereby licensed under the terms of the GNU GPL v2.0,
9 // pursuant to article 5.3.4 of the CeCILL v.2.1.
10 // This file was originally licensed under the terms of the CeCILL v2.1,
11 // and continues to be available under such terms.
12 // For more information, see the COPYING file which you should have received
13 // along with this program.
14
15 function [ind, matching] = vectorfind(m, v, way, joker, indType)
16     // way: "r" or 1, "c" or 2, 2 < i < ndims(m)
17     // indType: "" (default), "headIJK", "headN"
18
19     rhs = argn(2)
20
21     // Check number of inputs, at least 2 needed
22     if rhs < 2 | rhs > 5 then
23         msg = _("%s: Wrong number of input arguments: %d to %d expected.\n")
24         error(msprintf(msg, "vectorfind", 2, 5))
25     end
26
27     if ~isdef("joker", "l") | type(joker) == 0 then
28         joker = []
29     end
30
31     // Types of m and v must be compatible:
32     ok = (type(m) == type(v)) | ..
33     (or(type(m) == [1 5 8]) & or(type(v) == [1 5 8])) | ..
34     (type(m) == 4 & joker ~=[] & or(type(v) == [1 5 8]));
35
36     if ~ok then
37         msg = _("%s: Incompatible input arguments #%d and #%d: Same type expected.\n");
38         error(msprintf(msg, "vectorfind", 1, 2))
39     end
40
41     // Checking the needle:
42     if v == [] then
43         ind = []
44         return
45     end
46     if ~isvector(v) & ~isscalar(v) then
47         msg = _("%s: Wrong size for input argument #%d: Vector expected.\n")
48         error(msprintf(msg, "vectorfind", 2))
49     else
50         v = v(:);
51     end
52
53     Ndim = ndims(m);
54     // Checking the direction:
55     if ~isdef("way", "l") | type(way) == 0 then
56         way = 1         // Default direction
57     elseif size(way, "*") > 1 then
58         msg = _("%s: Argument #%d: Scalar (1 element) expected.\n");
59         error(msprintf(msg, "vectorfind", 3));
60     elseif and(type(way) ~= [1 10]) then
61         msg = _("%s: Argument #%d: Text or integer decimal expected.\n");
62         error(msprintf(msg, "vectorfind", 3));
63     elseif and(way ~= ["r" "c"]) & and(way ~= (1:Ndim)) then
64         msg = _("%s: Argument #%d: Must be in the set {%s}.\n");
65         tmp = msprintf("""r"",""c"",1,..,%d", Ndim);
66         error(msprintf(msg, "vectorfind",3, tmp))
67     else
68         if way == "r" then
69             way = 1
70         elseif way == "c" then
71             way = 2
72         end
73     end
74
75     // Checking the joker
76     if joker ~= [] & size(joker,"*") ~= 1 then
77         msg = _("%s: Wrong size for input argument #%d: Scalar expected.\n")
78         error(msprintf(msg, "vectorfind",4))
79     end
80
81     if joker ~= [] & ((or(type(joker) == [1 5 8]) & ~or(type(v) == [1 5 8])) | ..
82         (or(type(joker) == [2 10]) & type(joker) ~= type(v)))
83         msg = _("%s: Incompatible input arguments #%d and #%d: Same type expected.\n")
84         error(msprintf(msg, "vectorfind", 2, 4))
85     end
86     if type(m) == 4 & joker ~= [] & joker == 0 then
87         msg = _("%s: Argument #%d: non-zero number expected.\n")
88         error(msprintf(msg, "vectorfind", 4))
89     end
90     // Checkinh indType
91     if ~isdef("indType", "l") | type(indType) == 0 then
92         indType = ""
93     elseif type(indType) ~= 10 then
94         msg = _("%s: Argument #%d: Text expected.\n");
95         error(msprintf(msg, "vectorfind", 5));
96     elseif size(indType, "*") > 1 then
97         msg = _("%s: Argument #%d: Scalar (1 element) expected.\n");
98         error(msprintf(msg, "vectorfind", 5));
99     else
100         indType = convstr(indType);
101         if ~or(indType == ["" "headijk" "headn"]) then
102             msg = _("%s: Argument #%d: Must be in the set {%s}.\n");
103             error(msprintf(msg, "vectorfind", 5,""""", ""headIJK"", ""headN"""));
104         end
105     end
106
107     sm = size(m)
108     sv = size(v,"*")
109     matching = [];
110
111     // Normalizing m and v shapes: setting the direction along COLUMNS
112     if Ndim == 2 then
113         if way == 1 then
114             m = m.';
115         end
116     else    // Ndim>2
117         i = 1:Ndim;
118         if way == 1 then
119             i([1 2]) = [2 1];
120             m = permute(m,i);
121             m = matrix(m, sm(2), -1);
122         elseif way == 2 then
123             m = matrix(m, sm(1), -1)
124         elseif way==Ndim
125             m = matrix(m, -1, sm(way)).';
126         else
127             m = matrix(m, prod(sm(1:way-1)), sm(way), -1);
128             m = permute(m, [2 1 3]);
129             m = matrix(m, sm(way), -1);
130         end
131     end
132     // flag for partial v
133     partial = sv < size(m, 1)
134
135     //
136     if sv > size(m,1)
137         ind = []
138         return
139
140     elseif partial // v shorter than m size
141         cs = size(m,1) - sv
142         k = ndgrid(1:sv, 1:(cs+1)) + ones(sv,1)*(0:cs)
143         M = matrix(m(k(:),:), sv, -1)
144         ind = vectorfind(M, v, "c", joker)
145         if joker~=[]
146             matching = M(:,ind).';
147         end
148         clear M
149         // We compute k = index of start of matching inside the full range
150         if ind ~= [] then
151             k = modulo(ind, cs + 1);
152             k(k == 0) = cs + 1;
153         else
154             k = []
155         end
156         ind = ceil(ind / (cs + 1));  // linearized indices of full vectors
157
158     else
159         // Performing the detection:
160         // ------------------------
161         // Selecting only needle's components not being the joker
162         if joker ~= [] then
163             if type(v) == 1 & isnan(joker) then
164                 c = find(~isnan(v))
165             else
166                 c = find(v ~= joker)
167             end
168         else
169             c = 1:size(v, "*")
170         end
171         if type(m) == 4 & or(type(v) == [1 5 8]) then
172             v = (v ~= 0);
173         end
174
175         ind = 1:size(m, 2);
176         for k = c   // Loop over needle components not being jokers
177             if isnan(v(k)) then
178                 i = find(isnan(m(k, ind)));
179             else
180                 i = find(m(k, ind) == v(k));
181             end
182             ind = ind(i);
183             if ind == [] then
184                 break
185             end
186         end
187         if joker ~= [] then
188             matching = m(:,ind).';
189         end
190     end
191
192     // Post-processing output indices
193     // ------------------------------
194     if indType ~= "" | partial then
195         if way == 1 then
196             way = 2
197         elseif way == 2 then
198             way = 1
199         end
200         i = 1:Ndim;
201         i(way) = [];
202         I = ind2sub(sm(i), ind);
203         I2 = [];
204         if way>1
205             I2 = I(:, 1:way-1);
206         end
207         if ~partial then
208             I2 = [I2 ones(I(:,1))];
209         else
210             I2 = [I2 k(:)];
211         end
212         if way < Ndim then
213             I2 = [I2 I(:,way:$)];
214         end
215         if indType == "headijk" then
216             ind = I2;
217             if partial then
218                 [tmp, k] = gsort(ind(:,$:-1:1), "lr", "i");
219                 ind = ind(k,:);
220                 // otherwise: it should already be sorted
221             end
222         elseif indType=="headn" | partial then
223             ind = gsort(sub2ind(sm,I2)', "g", "i");
224         end
225     end
226 endfunction