scilab/modules/statistics/macros/median.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) 1999 - INRIA - Carlos Klimann
   3 // Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS
   4 // Copyright (C) 2016 - 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
  16 function y = median(x,orient)
  17     //
  18     // NOTES
  19     //    - modified by farid.belahcene:the case when x is an hypermatrix
  20     //    - new syntaxes: median(x,'m') and median(x,dim)
  21     //    S. Gougeon 2016:
  22     //     - Fixes overflow issues with encoded integers when interpolating
  23     //       http://bugzilla.scilab.org/14640
  24
  25     // CHECKING ARGUMENTS
  26     // ==================
  27     [lhs, rhs] = argn(0);
  28
  29     if rhs == 0 then
  30         msg = _("%s: Wrong number of input argument(s): %d or %d expected.\n")
  31         error(msprintf(msg, "median", 1, 2))
  32     end
  33
  34     if rhs<2 then
  35         orient = 0
  36     else
  37         if orient=="r" then
  38             orient = 1
  39         elseif orient=="c" then
  40             orient = 2
  41         elseif orient=="*" then
  42             orient = 0
  43         elseif orient=="m" then
  44             orient = find(size(x)>1,1)
  45             if orient==[] then
  46                 orient = 1
  47             end
  48         else
  49             if type(orient)<>1|size(orient,"*")<>1|~isreal(orient)|orient<=0 then
  50                 msg = _("%s: Wrong value for input argument #%d: ''%s'', ''%s'',''%s'' or a positive number expected.\n")
  51                 error(msprintf(msg, "median",2,"r","c","m"))
  52             end
  53         end
  54     end
  55
  56     // PROCESSING
  57     // ==========
  58     // median on all components
  59     // ------------------------
  60     if orient==0 then
  61         if x==[] then
  62             y = %nan
  63             return
  64         end
  65         n = size(x,"*");
  66         x = gsort(x(:),"g","i")
  67         if 2*int(n/2)==n then
  68             // avoid overflow: http://bugzilla.scilab.org/14640
  69             a = x(n/2)
  70             b = x(n/2+1)
  71             y = a/2 + b/2 + ((a-(a/2)*2) + (b-(b/2)*2))/2
  72         else
  73             y = x((n+1)/2);
  74         end
  75
  76     // Projection along a given direction / dimension
  77     // ----------------------------------------------
  78     else
  79         if x==[] then
  80             y = []
  81             return
  82         end
  83         if orient>ndims(x) then
  84             y = x
  85             return
  86         end
  87         xsize = size(x);
  88         if xsize(orient)==1 then
  89             y = x
  90             return
  91         end
  92         orient_above_size = xsize(orient+1:$);
  93         N = prod(orient_above_size)
  94         orient_below_size = xsize(1:orient-1);
  95         M = prod(orient_below_size)
  96         orient_size = xsize(1:orient);
  97         P = prod(orient_size)
  98         y = [];
  99         n = xsize(orient)
 100         for k = 1:N
 101             for i = 1:M
 102                 ytemp = gsort(x(i+(0:n-1)*M+(k-1)*P),"r","i")
 103                 if 2*int(n/2)==n then
 104                     // avoid overflow: http://bugzilla.scilab.org/14640
 105                     a = ytemp(n/2)
 106                     b = ytemp(n/2+1)
 107                     y = [ y ; a/2 + b/2 + ((a-(a/2)*2) + (b-(b/2)*2))/2]
 108                 else
 109                     y = [ y ; ytemp((n+1)/2)]
 110                 end
 111             end
 112         end
 113         xsize(orient) = 1
 114         y = matrix(y, xsize)
 115     end
 116 endfunction