d00379029ecab4d4b21d9becd3a4f8df0c85598e
[scilab.git] / 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