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