* Bug 15841 fixed: intersect() now accepts sparse
[scilab.git] / scilab / modules / elementary_functions / macros / intersect.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA - Serge Steer
3 // Copyright (C) DIGITEO - 2011 - Allan CORNET
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 // Copyright (C) 2019 - 2020 - Samuel GOUGEON
6 //
7 // This file is hereby licensed under the terms of the GNU GPL v2.0,
8 // pursuant to article 5.3.4 of the CeCILL v.2.1.
9 // This file was originally licensed under the terms of the CeCILL v2.1,
10 // and continues to be available under such terms.
11 // For more information, see the COPYING file which you should have received
12 // along with this program.
13
14 function [x_out, ka_out, kb_out] = intersect(a_in, b_in, orient)
15     // returns the vector of common values of two vectors
16
17     [lhs, rhs] = argn()
18     [x_out, ka_out, kb_out] = ([], [], [])
19
20     if rhs < 2 then
21         msg = gettext("%s: Wrong number of input arguments: %d or %d expected.\n")
22         error(msprintf(msg, "intersect", 2, 3))
23     end
24
25     // Special empty cases
26     [esp, espb] = (sparse([]), sparse(%f)); espb(1) = []
27     if a_in == [] || b_in == [] | isequal(a_in, esp) || isequal(b_in, esp) | ..
28         isequal(a_in, espb) | isequal(b_in, espb)
29         return
30     end
31     //
32     aIsComplex = or(type(a_in)==[1 5]) && ~isreal(a_in);
33     bIsComplex = or(type(b_in)==[1 5]) && ~isreal(b_in);
34
35     // Without orientation
36     // -------------------
37     if ~isdef("orient","l") then
38         if lhs == 1
39             a = unique(matrix(a_in, 1, -1));
40             b = unique(matrix(b_in, 1, -1));
41             if aIsComplex | bIsComplex
42                 x = gsort([a, b], "g", ["i" "i"], list(abs, atan));
43             else
44                 x = gsort([a, b], "g", "i");
45             end
46             keq = find( x(2:$) == x(1:$-1) ); // find consecutive equal values index
47             if keq <> [] then
48                 x_out = x(keq); //the intersection values in increasing order
49             end
50             return
51         end
52         // lhs > 1
53         // .......
54         [a, ka] = unique(matrix(a_in, 1, -1));
55         [b, kb] = unique(matrix(b_in, 1, -1));
56         kab = [ka, -kb];
57         // find duplicated values in [a_in,b_in],
58         // sort the array
59         if aIsComplex | bIsComplex then
60             [x, ksort] = gsort([a, b], "g", ["i" "i"], list(abs, atan));
61         else
62             [x, ksort] = gsort([a, b], "g", "i");
63         end
64         kab = kab(ksort); // apply [a_in,b_in] sorting permutation to kab
65         keq = find( x(2:$) == x(1:$-1) ); // find consecutive equal values index
66
67         if keq <> [] then
68             x_out = x(keq); // the intersection values in increasing order
69             // each duplicated value appear twice  and only twice and in
70             // consecutive positions keq(i) and keq(i)+1 in the sorted array x
71             kab = kab([keq keq+1])
72
73             // the positive values correspond to a_in index while the negative to b_in index.
74             ka_out = kab(kab > 0) // select index of intersection elements in a_in
75             // insure that a_in(ka_out)==x_out and b_in(kb_out)==x_out.
76             // I was'nt able to find a simple way.
77             if aIsComplex
78                 [s, k] = gsort(a_in(ka_out), "g", ["i" "i"], list(abs, atan));
79             else
80                 [s, k] = gsort(a_in(ka_out), "g", "i");
81             end
82             ka_out = ka_out(k)
83
84             if lhs > 2
85                 kb_out = -kab(kab < 0); //select index of intersection elements in b_in
86                 if bIsComplex
87                     [s, k] = gsort(b_in(kb_out), "g", ["i" "i"], list(abs, atan));
88                 else
89                     [s, k] = gsort(b_in(kb_out), "g", "i");
90                 end
91                 kb_out = kb_out(k);
92             end
93         end  // keq <> []
94         return
95     end  // rhs < 3
96
97     // WITH an ORIENTATION
98     // -------------------
99     if ndims(a_in) > 2 | ndims(b_in) > 2 then
100         msg = gettext("%s: Argument #%d: Orientation not supported for hypermatrix.\n")
101         error(msprintf(msg, "intersect", 3))
102     end
103     columnwise = orient==2 | orient=="c"
104     if columnwise then
105         a_in = a_in.'
106         b_in = b_in.'
107     end
108     if and(orient <> [1 2]) & and(orient <> ["r" "c"]) then
109         msg = gettext("%s: Argument #%d: Must be in the set {%s}.\n")
110         error(msprintf(msg, "intersect", 3, "1,''r'',2,''c''"))
111     end
112     // row-wise processing
113     // -------------------
114     if lhs == 1 then
115         a = unique(a_in, "r");
116         b = unique(b_in, "r");
117         if aIsComplex | bIsComplex
118             x = gsort([a ; b], "lr", ["i" "i"], list(abs, atan))
119         else
120             x = gsort([a ; b], "lr", "i")
121         end
122         keq = find(and(x(2:$,:) == x(1:$-1,:),"c")) // find index of consecutive equal values
123
124         if keq == [] then
125             return
126         else
127             x_out = x(keq,:); //the intersection values in increasing order
128         end
129
130     else
131         // ka requested : lhs > 1
132         // ............
133         [a, ka] = unique(a_in, "r");
134         [b, kb] = unique(b_in, "r");
135         kab = [ka; -kb];
136         if aIsComplex | bIsComplex
137             [x,ksort] = gsort([a; b], "lr", ["i" "i"], list(abs, atan));
138         else
139             [x,ksort] = gsort([a; b], "lr", "i");
140         end
141
142         kab = kab(ksort);//apply [a_in,b_in] sorting permutation to kab
143         keq = find(and(x(2:$,:) == x(1:$-1,:),"c")) // find index of consecutive equal values
144
145         if keq == [] then
146             return
147         else
148             x_out = x(keq,:); //the intersection values in increasing order
149
150             kab = kab([keq keq+1]);
151
152             //the positive values correspond to a_in index while the negative to b_in index.
153             ka_out = kab(kab>0); //select index of intersection elements in a_in
154             //insure that a_in(ka_out,:)==x_out and b_in(kb_out,:)==x_out.
155             //I was'nt able to find a simple way.
156             if aIsComplex
157                 [s, k] = gsort(a_in(ka_out,:), "lr", ["i" "i"], list(abs, atan));
158             else
159                 [s, k] = gsort(a_in(ka_out,:), "lr", "i");
160             end
161             ka_out = matrix(ka_out(k), 1, -1)
162
163             if lhs > 2 then
164                 kb_out = -kab(kab<0); //select index of intersection elements in b_in
165                 if bIsComplex
166                     [s, k] = gsort(b_in(kb_out,:), "lr", ["i" "i"], list(abs, atan));
167                 else
168                     [s, k] = gsort(b_in(kb_out,:), "lr", "i");
169                 end
170                 kb_out = matrix(kb_out(k),1,-1)
171             end
172         end
173     end
174
175     // Columnwise post-processing
176     // ..........................
177     if columnwise then
178         x_out = x_out.'
179     end
180 endfunction