scilab/modules/elementary_functions/macros/unique.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) INRIA
   3 // Copyright (C) 2012 - 2016 - Scilab Enterprises
   4 // Copyright (C) 2018 - 2020 - Samuel GOUGEON
   5 //
   6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
   7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
   8 // This file was originally licensed under the terms of the CeCILL v2.1,
   9 // and continues to be available under such terms.
  10 // For more information, see the COPYING file which you should have received
  11 // along with this program.
  12
  13 function [x, ki, ko, nb] = unique(x, varargin)
  14     // extract unique components of a vector, matrix, or hypermatrix
  15     // varargin : orient=1|2|"r"|"c", "uniqueNan", "keepOrder"
  16     //
  17     // History:
  18     // * 2019 - S. Gougeon :
  19     //   - add uniqueNan option: http://bugzilla.scilab.org/15522
  20     //   - add keepOrder option: http://bugzilla.scilab.org/15795
  21     //   - add nb output option: http://bugzilla.scilab.org/8418
  22     // * 2020 - S. Gougeon :
  23     //   - add ku output indices: http://bugzilla.scilab.org/16337
  24
  25     keepOrder = %f
  26     uniqueNan = %f
  27     orient = "*"
  28     newInf = [] // init Inf substitute in case of "uniqueNan" and or(x==%inf)
  29     ki = []
  30     ko = []
  31     nb = []
  32     sx = size(x)
  33
  34     // CHECKING INPUT ARGUMENTS
  35     // ------------------------
  36     in = varargin
  37     i = 2;  // index of the current input argument
  38     if size(in)>0 then
  39         a = in(1)
  40         if type(a)==0
  41             a = "*"
  42         end
  43         if typeof(a)=="string"
  44             a = convstr(a)
  45         end
  46         select a
  47         case 1
  48             orient = "r"
  49         case 2
  50             orient = "c"
  51         case "uniquenan"
  52             uniqueNan = %t
  53         case "keeporder"
  54             keepOrder = %t
  55         case "r"
  56             orient = "r"
  57         case "c"
  58             orient = "c"
  59         else
  60             msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
  61             error(msprintf(msg, "unique", i, "1,2,""r"",""c"",""keepOrder"",""uniqueNan"""))
  62         end
  63         in(1) = null()
  64         i = 3
  65     end
  66     if or(orient==["r" "c"]) & ndims(x)>2 then
  67         msg = _("%s: Argument #%d: ''%s'' not allowed for an hypermatrix.\n")
  68         error(msprintf(msg, "unique", 2, orient))
  69     end
  70
  71     while size(in)>0 & i<5 then
  72         a = in(1)
  73         if typeof(a)=="string"
  74             a = convstr(a)
  75         end
  76         select a
  77         case "uniquenan"
  78             uniqueNan = %t
  79         case "keeporder"
  80             keepOrder = %t
  81         else
  82             msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
  83             error(msprintf(msg, "unique", i, """keepOrder"",""uniqueNan"""))
  84         end
  85         in(1) = null()
  86         i = i+1
  87     end
  88     uniqueNan = uniqueNan & or(type(x)==[1 5])
  89
  90     sz = size(x);
  91     if size(x, orient)==1 then
  92         ki = 1
  93         ko = 1
  94         nb = 1
  95         return
  96     end
  97     if uniqueNan
  98         [x, newInf] = uniqueProcessNan(x, [], "removeNan")
  99     end
 100     getK = argn(1) > 1 | keepOrder
 101
 102
 103     // [] trivial case
 104     // ---------------
 105     if isempty(x) then
 106         return  // ki, ko, nb are already []. x is [] or sparse([])
 107     end
 108
 109     // PROCESSING
 110     // ----------
 111     areComplex = or(type(x)==[1 5]) && ~isreal(x, 0)
 112     if orient=="*" then
 113         if getK then
 114             // ki (begin)
 115             if areComplex
 116                 [x,ki] = gsort(x,"g",["i" "i"], list(abs, atan));
 117             else
 118                 [x,ki] = gsort(x,"g","i");
 119             end
 120             keq = x(2:$) == x(1:$-1);
 121             if argn(1) > 2
 122                 // ko (begin)
 123                 ko(ki) = cumsum(~[%F ; keq(:)])
 124                 if ko <> [] then
 125                     ko = matrix(ko, sx)
 126                 end
 127                 // nb
 128                 if argn(1) > 3
 129                     nb = [0 find(~keq) size(x,"*")]
 130                     nb = nb(2:$) - nb(1:$-1)
 131                 end
 132             end
 133             // ki (end)
 134             keq = find(keq);
 135             if keq <> [] then keq = keq + 1; end
 136             x(keq) = [];
 137             ki(keq) = [];
 138         else
 139             if areComplex
 140                 x = gsort(x,"g",["d" "d"], list(abs, atan));
 141             else
 142                 x = gsort(x,"g","d");
 143             end
 144             x = x($:-1:1);
 145             x( find(x(2:$) == x(1:$-1)) ) = [];
 146         end
 147         if sx(1) > 1 | length(sx) > 2
 148             x = x(:)
 149             ki = ki(:)
 150             nb = nb(:)
 151         end
 152
 153     else // orient used
 154         if  orient==2 | orient=="c" then
 155             x = x.'
 156         end
 157         if getK then
 158             if areComplex
 159                 [x,ki] = gsort(x,"lr",["i" "i"], list(abs, atan));
 160             else
 161                 [x,ki] = gsort(x,"lr","i");
 162             end
 163             keq = and(x(2:$,:) == x(1:$-1,:),"c")
 164             if argn(1)>2
 165                 ko(ki) = cumsum(~[%F ; keq(:)])
 166                 // nb
 167                 if argn(1) > 3
 168                     nb = [0 find(~keq) size(x,1)]
 169                     nb = nb(2:$) - nb(1:$-1)
 170                     nb = nb(:)
 171                 end
 172             end
 173             keq = find(keq)
 174             if keq <> [] then keq = keq + 1;end
 175             x(keq,:) = [];
 176             ki(keq,:) = [];
 177         else
 178             if areComplex
 179                 x = gsort(x,"lr",["i" "i"], list(abs, atan));
 180             else
 181                 x = gsort(x,"lr","i");
 182             end
 183             x( find(and(x(2:$,:) == x(1:$-1,:),"c")),:) = [];
 184         end
 185         if  orient==2 | orient=="c" then
 186             x = x.'
 187             ki = matrix(ki, 1, -1)
 188             ko = ko'
 189             nb = nb'
 190         end
 191     end
 192
 193     if uniqueNan
 194         x = uniqueProcessNan(x, newInf, "restoreNan")
 195     end
 196
 197     if keepOrder
 198         [ki, kk] = gsort(ki,"g","i")
 199         select orient
 200         case "*"
 201             x = x(kk)
 202         case "r"
 203             x = x(kk,:)
 204         case "c"
 205             x = x(:,kk)
 206         end
 207         if argn(1)>2
 208             nb = nb(kk)
 209             [?, kk2] = gsort(kk,"g","i")
 210             ko = kk2(ko)
 211         end
 212     end
 213 endfunction
 214
 215 // -------------------------------------------------------------------
 216
 217 // To consider Nan mutually equal, we replace all of them with a "regular" substitute.
 218 // Since Nan are sorted as > Inf, we must use anyway Inf as the Nan substitute.
 219 // If the original array have already some Inf, we must priorly replace them with
 220 // a decimal greater than the finite maximum of the array values.
 221 // After processing, we restore Inf => Nan (, and maxNum => Inf).
 222
 223 function [x, newInf] = uniqueProcessNan(x, newInf, way)
 224
 225     if way=="removeNan" & or(isnan(x)) then
 226         // Replacing Nan
 227         // -------------
 228         if isreal(x)
 229             if or(x==%inf)
 230                 b = x==%inf
 231                 m = max([1 max(x(~b))])
 232                 newInf = m*1.5
 233                 x(b) = newInf
 234             end
 235             x(x<>x) = %inf
 236         else
 237             r = real(x)
 238             i = imag(x)
 239             if or([r i]==%inf)
 240                 br = r==%inf
 241                 m = max(r(~br),1)
 242                 bi = i==%inf
 243                 m = max(i(~bi),m,1)
 244                 newInf = m*1.5
 245                 r(br) = newInf
 246                 i(bi) = newInf
 247             end
 248             r(r<>r) = %inf
 249             i(i<>i) = %inf
 250             x = complex(r,i);
 251         end
 252
 253     // Restoring  NaN
 254     // --------------
 255     elseif way=="restoreNan"
 256         if isreal(x)
 257             x(x==%inf) = %nan
 258             if newInf~=[]
 259                 x(x==newInf) = %inf
 260             end
 261         else
 262             r = real(x)
 263             r(r==%inf) = %nan
 264             i = imag(x)
 265             i(i==%inf) = %nan
 266             if newInf~=[]
 267                 r(r==newInf) = %inf
 268                 i(i==newInf) = %inf
 269             end
 270             x = complex(r, i)
 271         end
 272     end
 273 endfunction