* Bug 15734 fixed: intersect() dit not support complex numbers
[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 = []
19     ka_out = []
20     kb_out = []
21
22     if rhs < 2 then
23         msg = gettext("%s: Wrong number of input arguments: %d or %d expected.\n")
24         error(msprintf(msg, "intersect", 2, 3))
25     end
26
27     if a_in == [] | b_in == [] then
28         return
29     end
30     aIsComplex = type(a_in)==1 && ~isreal(a_in);
31     bIsComplex = type(b_in)==1 && ~isreal(b_in);
32
33     // Without orientation
34     // -------------------
35     if ~isdef("orient","l") then
36         if lhs == 1
37             a = unique(matrix(a_in, 1, -1));
38             b = unique(matrix(b_in, 1, -1));
39             if aIsComplex | bIsComplex
40                 x = gsort([a, b], "g", ["i" "i"], list(abs, atan));
41             else
42                 x = gsort([a, b], "g", "i");
43             end
44             keq = find( x(2:$) == x(1:$-1) ); // find consecutive equal values index
45             if keq <> [] then
46                 x_out = x(keq); //the intersection values in increasing order
47             end
48             return
49         end
50         // lhs > 1
51         // .......
52         [a, ka] = unique(matrix(a_in, 1, -1));
53         [b, kb] = unique(matrix(b_in, 1, -1));
54         kab = [ka, -kb];
55         // find duplicated values in [a_in,b_in],
56         // sort the array
57         if aIsComplex | bIsComplex then
58             [x, ksort] = gsort([a, b], "g", ["i" "i"], list(abs, atan));
59         else
60             [x, ksort] = gsort([a, b], "g", "i");
61         end
62         kab = kab(ksort); // apply [a_in,b_in] sorting permutation to kab
63         keq = find( x(2:$) == x(1:$-1) ); // find consecutive equal values index
64
65         if keq <> [] then
66             x_out = x(keq); // the intersection values in increasing order
67             // each duplicated value appear twice  and only twice and in
68             // consecutive positions keq(i) and keq(i)+1 in the sorted array x
69             kab = kab([keq keq+1])
70
71             // the positive values correspond to a_in index while the negative to b_in index.
72             ka_out = kab(kab > 0) // select index of intersection elements in a_in
73             // insure that a_in(ka_out)==x_out and b_in(kb_out)==x_out.
74             // I was'nt able to find a simple way.
75             if aIsComplex
76                 [s, k] = gsort(a_in(ka_out), "g", ["i" "i"], list(abs, atan));
77             else
78                 [s, k] = gsort(a_in(ka_out), "g", "i");
79             end
80             ka_out = ka_out(k)
81
82             if lhs > 2
83                 kb_out = -kab(kab < 0); //select index of intersection elements in b_in
84                 if bIsComplex
85                     [s, k] = gsort(b_in(kb_out), "g", ["i" "i"], list(abs, atan));
86                 else
87                     [s, k] = gsort(b_in(kb_out), "g", "i");
88                 end
89                 kb_out = kb_out(k);
90             end
91         end  // keq <> []
92         return
93     end  // rhs < 3
94
95     // WITH an ORIENTATION
96     // -------------------
97     if ndims(a_in) > 2 | ndims(b_in) > 2 then
98         msg = gettext("%s: Argument #%d: Orientation not supported for hypermatrix.\n")
99         error(msprintf(msg, "intersect", 3))
100     end
101     columnwise = orient==2 | orient=="c"
102     if columnwise then
103         a_in = a_in.'
104         b_in = b_in.'
105     end
106     if and(orient <> [1 2]) & and(orient <> ["r" "c"]) then
107         msg = gettext("%s: Argument #%d: Must be in the set {%s}.\n")
108         error(msprintf(msg, "intersect", 3, "1,''r'',2,''c''"))
109     end
110     // row-wise processing
111     // -------------------
112     if lhs == 1 then
113         a = unique(a_in, "r");
114         b = unique(b_in, "r");
115         if aIsComplex | bIsComplex
116             x = gsort([a ; b], "lr", ["i" "i"], list(abs, atan))
117         else
118             x = gsort([a ; b], "lr", "i")
119         end
120         keq = find(and(x(2:$,:) == x(1:$-1,:),"c")) // find index of consecutive equal values
121
122         if keq == [] then
123             return
124         else
125             x_out = x(keq,:); //the intersection values in increasing order
126         end
127
128     else
129         // ka requested : lhs > 1
130         // ............
131         [a, ka] = unique(a_in, "r");
132         [b, kb] = unique(b_in, "r");
133         kab = [ka; -kb];
134         if aIsComplex | bIsComplex
135             [x,ksort] = gsort([a; b], "lr", ["i" "i"], list(abs, atan));
136         else
137             [x,ksort] = gsort([a; b], "lr", "i");
138         end
139
140         kab = kab(ksort);//apply [a_in,b_in] sorting permutation to kab
141         keq = find(and(x(2:$,:) == x(1:$-1,:),"c")) // find index of consecutive equal values
142
143         if keq == [] then
144             return
145         else
146             x_out = x(keq,:); //the intersection values in increasing order
147
148             kab = kab([keq keq+1]);
149
150             //the positive values correspond to a_in index while the negative to b_in index.
151             ka_out = kab(kab>0); //select index of intersection elements in a_in
152             //insure that a_in(ka_out,:)==x_out and b_in(kb_out,:)==x_out.
153             //I was'nt able to find a simple way.
154             if aIsComplex
155                 [s, k] = gsort(a_in(ka_out,:), "lr", ["i" "i"], list(abs, atan));
156             else
157                 [s, k] = gsort(a_in(ka_out,:), "lr", "i");
158             end
159             ka_out = matrix(ka_out(k), 1, -1)
160
161             if lhs > 2 then
162                 kb_out = -kab(kab<0); //select index of intersection elements in b_in
163                 if bIsComplex
164                     [s, k] = gsort(b_in(kb_out,:), "lr", ["i" "i"], list(abs, atan));
165                 else
166                     [s, k] = gsort(b_in(kb_out,:), "lr", "i");
167                 end
168                 kb_out = matrix(kb_out(k),1,-1)
169             end
170         end
171     end
172
173     // Columnwise post-processing
174     // ..........................
175     if columnwise then
176         x_out = x_out.'
177     end
178 endfunction