5c8b7583e57abaa3f3a2af22d11ef65f16b89b09
[scilab.git] / scilab / modules / statistics / macros / histc.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2017 - 2020 - Samuel GOUGEON
3 //
4 // This file is hereby licensed under the terms of the GNU GPL v2.0,
5 // pursuant to article 5.3.4 of the CeCILL v.2.1.
6 // This file was originally licensed under the terms of the CeCILL v2.1,
7 // and continues to be available under such terms.
8 // For more information, see the COPYING file which you should have received
9 // along with this program.
10
11
12 function [heights, jokers, binsOut, ind] = histc(data, bins, options)
13 // Scilab 5.5:
14 //  1) histc(bins, data, normalization)
15 //  2) By défaut, histc() worked in "densityNorm,normWith:all" mode
16 // Scilab 6.0:
17 //  1) histc(data, bins, options)
18 //  2) By default, histc() works in "counts" mode
19
20 //    histc(data)
21 //    histc(data, nbins)
22 //    histc(data, binsWidth<0)
23 //    histc(data, binsEdges)
24 //    histc(data, binsValues [, "discrete"])
25 //    histc(data, binsAlgo)
26 //    histc(data, ..., options)
27 //    histc(nbins, data)                        // deprecated
28 //    histc(binsEdges, data)                    // deprecated
29 //    histc(binsEdges, data, normalization=%t)  // deprecated
30 //     heights                        = histc(data ..)
31 //    [heights, jokers]               = histc(data ..)
32 //    [heights, jokers, binsOut]      = histc(data ..)
33 //    [heights, jokers, binsOut, ind] = histc(data ..)
34 //
35 //    Binning algorithm:
36 //      "sqrt" | "sturges" | "freediac"
37 //    Binning options:
38 //      "discrete"  // , "smartbins"
39 //    Scale options:
40 //      "counts" (default) | "countsNorm" | "density" | "densityNorm"
41 //    Normalization "normWith: " options:  Priorities: "all" > "out" > ..
42 //     decimal, complex numbers: "leftout rightout out inf nan zeros all"
43 //                     integers: "leftout rightout (out=) all"
44 //                         Text: "empty out all"
45 //                  polynomials: "(out=) all"
46 //       vector or single string with space-separated flags (case-unsensitive)
47 //
48 //    Complex numbers are sorted by real parts, then by imaginary parts
49 //    jokers =
50 //      * integers: [leftout, rightout] (continuous) or [out] (discrete)
51 //      * decimal or complex numbers, full or sparse:
52 //            continuous: [leftout, rightout, zeros, nan, inf]
53 //              discrete: [out, 0, zeros, nan, inf]
54 //      * polynomials: [out, 0, 0, nan, inf]
55 //      * texts: continuous: [leftout, rightout, ""]; discrete: [out, 0, ""]
56 //
57
58     [lhs, rhs] = argn();
59     fname = "histc"
60     [heights, jokers, binsOut, ind] = ([],[],[],[])
61     [Nout, Nleftout, Nzeros, Nrightout, Ninf, Nnan] = (0,0,0,0,0,0)
62
63     // PARSING INPUT ARGUMENTS
64     // =======================
65     if rhs<1 | rhs>3 then
66         msg = _("%s: Wrong number of input arguments: %d to %d expected.\n")
67         error(msprintf(msg, fname, 1, 3))
68     end
69     // Old 5.5.x mode management
70     // -------------------------
71     mode55 = %f
72     // Check for Scilab 5.5 back-compatibility
73     if (isdef("data","l") & isdef("bins", "l") & type(bins)>0) & ..
74         size(data,"*")< size(bins,"*") & or(type(data)==[1 8]) & type(bins)~=10
75         msg = _("%s: data are now expected in argument #1 => arguments #1 and #2 switched\n")
76         warning(msprintf(msg, fname))
77         mode55 = %t
78         [data, bins] = (bins, data)
79         tmp = isdef("options","l") && type(options)==4
80         normarg = %f
81         norma = %f
82         if isdef("normalization","l")
83             normarg = %t
84             if type(normalization)==4 then
85                 norma = normalization(1)
86             end
87         else
88             if rhs<3 then
89                 norma = %t
90             else
91                 if isdef("options","l") & type(options)==4
92                     normarg = %t
93                     norma = options(1)
94                 end
95             end
96         end
97         if norma then
98             options = ["densityNorm" "normWith:all"]
99         else
100             options = "counts"
101         end
102         if normarg
103             msg = _("%s: ""normalization"" boolean argument #3 is obsolete. ""%s"" option used.\n")
104             warning(msprintf(msg, fname, options(1)))
105         end
106     end
107     // -----------------------
108     if ~isdef("data","l") | type(data)==0 then
109         msg = _("%s: Argument #%d: %s\n")
110         error(msprintf(msg, fname, 1, "data array missing"))
111     end
112     if ~isdef("bins","l") | type(bins)==0 then
113         bins = []
114     end
115     if ~isdef("options","l") then
116         options = "counts"
117     end
118     // Now, bins, data and options are defined. Let's check their types and sizes
119
120     // CHECKING INPUT ARGUMENTS AND SETTING DEFAULTS
121     // =============================================
122     // Options
123     // -------
124     if options~=[]
125         if type(options)~=10 then
126             msg = _("%s: Argument #%d: Text expected.\n")
127             error(msprintf(msg, fname, 3))
128         end
129         options = tokens(strcat(convstr(options),","),",")
130     end
131     discrete = grep(options, "discrete")~=[] | type(data)==2 // polynomials
132
133     // Data
134     // ----
135     if isempty(data) then
136         return
137     end
138     if ~or(type(data)==[1 2 5 8 10]) then
139         msg = _("%s: Data in argument #%d: Numbers or polynomials or texts expected.\n")
140         error(msprintf(msg, fname, 1))
141     end
142     data0 = data    // Saving
143     // If encoded integers (including u-int64): => decimal (for dsearch(), etc)
144     if type(data) == 8 then
145         data = double(data)
146     end
147     // If sparse numbers: we remove zeros
148     if type(data)==5
149         [IJ, data] = spget(data)
150     end
151     // If there are only %inf or %nan => return ==!!!!!!!!!!!!!!!!!!!!!!!!!!!
152     if ~discrete & or(type(data0)==[1 5]) then
153         L = data
154         tmp = isnan(data)
155         Nnan = sum(tmp)
156         L(tmp) = []
157         tmp = isinf(data)
158         Ninf = sum(tmp)
159         L(tmp) = []
160         if data==[]
161             msg = _("%s: data have only %inf or/and %nan values.\n")
162             warning(msprintf(msg, fname))
163             if type(data0)==1
164                 Nzeros = sum(data==0)
165             else
166                 Nzeros = size(data0,"*") - Nnan - Ninf
167             end
168             Nleftout = 0
169             Nrightout = 0
170             if length(bins)>1
171                 b = bins(~isinf(bins))
172                 if 0<min(b)
173                     Nleftout = Nzeros
174                 end
175                 if 0 > max(b)
176                     Nrightout = Nzeros
177                 end
178             end
179             jokers = [Ninf Nnan Nzeros Nleftout Nrightout]
180             return
181         end
182     end
183     // Polynomials: setting Nan and Ninf
184     if type(data)==2
185         L = data
186         tmp = isnan(L)
187         Nnan = sum(tmp)
188         L(tmp) = []     // to ensure excluding any intersection
189         tmp = isinf(L)
190         Ninf = sum(tmp)
191     end
192
193     // Bins
194     // ----
195     // Setting a default binning mode
196     if bins==[] then
197         // histc(data): Default bins:
198         //  * strings: discrete, from unique()
199         //  * polynomials: discrete, from defBins set hereabove
200         //  * encoded integers: "sqrt". if strange() < "sqrt" => binsWidth = 1
201         //  * decimal, complex: binMethod option or default "sqrt"
202         if type(data)==10       // texts
203             discrete = %t
204             bins = unique(data)(:)'
205             bins = bins(:)'
206
207         elseif type(data)==2    // Polynomials
208             L = data
209             tmp = isnan(L)
210             L(tmp) = []
211             bins = []
212             while L~=[]
213                 bins = [bins L(1)]
214                 L(L==L(1)) = []
215             end
216             clear L
217
218         else    // decimals, complexes, encoded integers
219             if ~discrete
220                 bins = "sqrt"
221             else
222                 bins = unique(data)(:).'
223             end
224         end
225     end
226     //
227     if ~discrete
228         if type(bins)==8 | (type(bins)==1 & ~isreal(bins))
229             bins = real(double(bins))
230         end
231         if size(bins,"*")==1 & type(bins)==10 then
232             // binsAlgo
233             // --------
234             if ~or(convstr(bins)==["sqrt" "sturges" "freediac"])
235                 msg = _("%s: Argument #%d: wrong value for binning algo")
236                 error(msprintf(msg, fname, 2))
237             end
238             // We set the number of bins. Their edges are set later.
239             Nvalid = size(data,"*") - Nnan - Ninf
240             select bins
241             case "sturges"
242                 bins = ceil(1 + log2(Nvalid))
243             case "freediac"
244                 tmp = real(data)
245                 binWidth = 2*iqr(tmp)* Nvalid^(-1/3)
246                 if binWidth~=0
247                     bins = round((max(tmp) - min(tmp)) / binWidth)
248                 else
249                     bins = 2
250                 end
251             else  // "sqrt"
252                 bins = max(1, round(sqrt(Nvalid)))
253             end
254         end
255
256         if length(bins)==1 & type(bins)==1  then
257             // The number of bins or their width is provided
258             if type(data)==10   // texts
259                 msg = _("%s: Argument #2: Please provide bins edges or values or leave choosing default bins.\n")
260                 error(msprintf(msg, fname))
261             end
262             // Case polynomials: not possible here <= discrete previously set to %t
263             // Cases: Encoded integers, decimal or complex numbers
264             if bins>=0 & (bins~=floor(bins) | bins==0) | isnan(bins) | bins==%inf
265                 msg = _("%s: Argument #%d: non-zero decimal integer expected.\n")
266                 error(msprintf(msg, fname, 2))
267             elseif bins==-%inf
268                 msg = _("%s: Argument #%d: decimal number > -Inf expected.\n")
269                 error(msprintf(msg, fname, 2))
270             end
271
272             // If data are complex-encoded, they will be sorted by increasing
273             //  real parts, then by increasing imaginary parts.
274             b = data(~isnan(data) & ~isinf(data))
275             b = real(b)
276             mind = min(b)
277             maxd = max(b)
278             if type(data0)==1
279                 if bins > 0
280                     // Required number of bins
281                     if (mind == maxd) then
282                         mind = mind*(1 - 0.5/bins);
283                         maxd = maxd*(1 + 0.5/bins);
284                     end
285                     bins = linspace(mind, maxd, bins+1);          // Bins edges
286                 else
287                     // Required bins width
288                     if (mind == maxd) then
289                         bins = mind + [1 -1]*bins/2
290                     else
291                         // Edges set à 0.0, width = |bins|, from mind- to maxd+
292                         tmp = -bins
293                         bins = (floor(mind/tmp):ceil(maxd/tmp))*tmp
294                     end
295                 end
296             else
297                 // Encoded integers
298                 if bins > 0
299                     // Required number of bins
300                     tmp = maxd - mind + 1
301                     if tmp <= bins & tmp <= 2^16  // integer-centered bins of width==1
302                         bins = (mind:(maxd+1)) - 0.5              // Bins edges
303                     else
304                         ds = (maxd - mind)/bins/2
305                         bins = linspace(mind-ds, maxd+ds, bins+1) // Bins edges
306                     end
307                 else
308                     // Required bins width
309                     tmp = max(round(-bins),1) // Bins edges @ -0.5 modulo |int(bins)|
310                     bins = (floor(mind/tmp):ceil(maxd/tmp))*tmp - 0.5 // Bins edges
311                 end
312             end
313
314         else
315             // Bins edges are provided: we check that edges are sorted
316             // Texts
317             // -----
318             if type(data0)==10
319                 if type(bins)~=10
320                     msg = _("%s: Arguments #%d and #%d: Same types expected.\n")
321                     error(msprintf(msg, fname, 1, 2))
322                 end
323                 if ~and(bins==unique(bins))
324                     msg = _("%s: Argument #%d: Elements must be sorted in increasing order.\n")
325                     error(msprintf(msg, fname, 1))
326                 end
327
328             // Numbers
329             // -------
330             else
331                 if ~or(type(bins)==[1 8])
332                     msg = _("%s: Arguments #%d and #%d: Same types expected.\n")
333                     error(msprintf(msg, fname, 1, 2))
334                 end
335                 bins = bins(~isnan(bins))
336                 if min(diff(bins)) <= 0 then
337                     msg = _("%s: Argument #%d: Elements must be in increasing order.\n")
338                     error(msprintf(msg, fname, 1))
339                 end
340                 // We manage -%inf open left bin
341                 if bins(1)==-%inf
342                     tmp = min(data(~isinf(data)))
343                     if bins(2)<=tmp
344                         bins(1) = []
345                     else
346                         bins(1) = tmp
347                     end
348                 end
349                 // We manage %inf open right bin
350                 if bins($)==%inf
351                     tmp = max(data(~isinf(data)))
352                     if bins($-1)>=tmp
353                         bins($) = []
354                     else
355                         bins($) = tmp
356                     end
357                 end
358             end
359         end
360          // if ~discrete
361
362     else    // if discrete
363         if type(data0)==10
364             if type(bins)~=10
365                 msg = _("%s: Arguments #%d and #%d: Same types expected.\n")
366                 error(msprintf(msg, fname, 1, 2))
367             end
368             [bins, dataOrder] = unique(bins)
369
370         elseif or(type(data)==[1 8])
371             if ~or(type(bins)==[1 8])
372                 msg = _("%s: Arguments #%d and #%d: Same types expected.\n")
373                 error(msprintf(msg, fname, 1, 2))
374             end
375             bins = real(double(bins))
376             bins(isnan(bins)) = []
377             [bins, dataOrder] = unique(bins)
378         end
379     end
380     binsOut = bins
381
382     // PROCESSING
383     // ==========
384     // jokers =
385     //  * integers: [leftout, rightout] (continuous) or [out] (discrete)
386     //  * decimal or complex numbers, full or sparse:
387     //       continuous: [inf, nan, zeros, leftout, rightout]
388     //         discrete: [inf, nan, zeros, out]
389     //  * polynomials: [inf, nan, out]
390     //  * texts: continuous: [empty, leftout, rightout]; discrete: [empty, out]
391     if type(data0)~=2 then
392         hmode = "c"            // continuous
393         if discrete then
394             hmode = "d"        // discrete
395         end
396         if or(type(data0)==[1 5])
397             if type(data0)==1
398                 tmp = data0
399             else
400                 tmp = data     // sparse
401             end
402             tmp1 = isinf(real(tmp))
403             tmp2 = isinf(imag(tmp))
404             tmp3 = isnan(tmp) | (tmp1 & tmp2)
405             if ~discrete & ~isreal(tmp,0)
406                 data(tmp1 | tmp2) = %inf
407                 data(tmp3) = %nan
408             end
409             clear tmp tmp1 tmp2 tmp3
410             data = real(data)
411         end
412         [ind, heights] = dsearch(data, bins, hmode)
413                                      // ind: bin's number of each data component
414
415         if discrete
416             // We restore the original order of binsValues:
417             [tmp, k]  = gsort(dataOrder, "g","i")
418             [tmp, k2] = gsort(k,"g","i")
419             bins = bins(k)
420             binsOut = bins
421             heights = heights(k)
422             tmp2 = ind~=0
423             ind(tmp2) = k2(ind(tmp2))
424         end
425
426     else // Polynomials
427         // Counts (and location)
428          ind = zeros(data)
429          heights = []
430          for i = 1:size(bins,"*")
431              tmp = data==bins(i)
432              heights(1,i) = sum(tmp)
433              ind(tmp) = i
434          end
435          tmp = Ninf - sum(heights(isinf(bins)))     // inf out of the selection
436          Nout = size(data,"*") - sum(heights) - Nnan - tmp// excludes all nan and inf not in the selection
437          jokers = [Nout, 0, 0, Nnan, Ninf]
438     end
439
440     // POST-PROCESSING
441     // ===============
442     // Counting jokers
443     // ---------------
444     if or(type(data0)==[1 5])        // Real and complex numbers, sparse or not
445         inf = isinf(data)
446         Ninf = sum(inf)
447         Nnan = sum(isnan(data))
448         if type(data0)==1
449             Nzeros = sum(data==0)
450         else
451             Nzeros = size(data0,"*") - nnz(data0)
452         end
453         if ~discrete                 // Sparses: 0 is not taken into account
454             Nleftout  = sum(data(~inf) < bins(1))
455             Nrightout = sum(data(~inf) > bins($))
456             jokers = [Nleftout Nrightout Nzeros Nnan Ninf]
457         else
458             Nout = size(data0,"*") - sum(heights) - Nnan - Ninf
459             binsinf = isinf(bins)
460             if or(binsinf)
461                 Nout = Nout + sum(heights(binsinf))
462             end
463             if type(data0)==5
464                 Nout = Nout - Nzeros
465             end
466             jokers = [Nout 0 Nzeros Nnan Ninf]
467         end
468
469     elseif type(data0)==8            // Encoded integers
470         if ~discrete
471             Nleftout  = sum(data < bins(1))
472             Nrightout = sum(data > bins($))
473             jokers = [Nleftout Nrightout]
474             Nout = Nleftout + Nrightout
475         else
476             Nout = size(data,"*") - sum(heights)
477             jokers = Nout
478         end
479
480     elseif type(data0)==10           // Texts
481         Nempty = sum(data0=="")
482         Nout = size(data0,"*") - sum(heights)
483         if ~discrete
484             tmp = gsort(data,"g","i")
485             Nleftout = sum(find(strcmp(tmp,bins(1))<0)($)) - Nempty
486             tmp = sum(find(strcmp(tmp,bins($))>0)(1))
487             if ~tmp
488                 Nrightout = 0
489             else
490                 Nrightout = size(data,"*") - tmp + 1
491             end
492             clear tmp
493             jokers = [Nleftout, Nrightout, Nempty]
494         else
495             if ~or(bins=="")
496                 Nout = Nout - Nempty
497             end
498             jokers = [Nout, 0, Nempty]
499         end
500
501     // Polynomials: jokers already set when setting bins
502     end
503
504     // Tuning the histogram scale
505     // --------------------------
506     histmode = "counts"
507     if discrete then
508         tmp = grep(options, ["counts" "countsnorm"])
509     else
510         tmp = grep(options, ["counts" "countsnorm" "density" "densitynorm"])
511     end
512     if tmp~=[] then
513         histmode = options(tmp($))
514     end
515     // jokers =
516     //  * integers: [leftout, rightout] (continuous) or [out] (discrete)
517     //  * decimal or complex numbers, full or sparse:
518     //       continuous: [leftout, rightout, zeros, nan, inf]
519     //         discrete: [out, 0, zeros, nan, inf]
520     //  * polynomials: [out, 0, 0, nan, inf].
521     //  * texts: continuous: [leftout, rightout, empty]; discrete: [out, empty]
522
523     // Number of components to take into account:
524     // Parsing the option "normWith:
525     //  Continuous: " leftout rightout inf nan zeros"
526     //    Discrete: " out inf nan zeros"
527     if grep(histmode, "norm") then
528         N = sum(heights)
529         tmp = grep(options, "normwith")
530         normMode = options(tmp)
531         if normMode~=[]
532             gh = "/[\:|\s+]"    // grep() header
533             if ~discrete
534                 if or(type(data0)==[1 2 5])
535                     if grep(normMode, gh+["inf" "all"]+"/", "r")~=[]
536                         N = N + Ninf
537                     end
538                     if grep(normMode, gh+["nan" "all"]+"/", "r")~=[]
539                         N = N + Nnan
540                     end
541                     if type(data0)==5 & ..
542                        grep(normMode,gh+["zeros" "all"]+"/", "r")~=[]
543                         N = N + Nzeros
544                     end
545                 end
546                 if grep(normMode, gh+["leftout" "out" "all"]+"/", "r")~=[]
547                     N = N + Nleftout
548                 end
549                 if grep(normMode, gh+["rightout" "out" "all"]+"/", "r")~=[]
550                     N = N + Nrightout
551                 end
552                 if type(data0)==10 & grep(normMode, gh+["empty" "all"]+"/", "r")~=[]
553                     N = N + Nempty
554                 end
555             // Normalized discrete mode:
556             else
557                 if grep(normMode, gh+"all/", "r")
558                     N = size(data0,"*")
559                 else
560                     if grep(normMode, gh+"out/", "r")
561                         N = N + Nout
562                     end
563                     if or(type(data0)==[1 2 5])
564                         if grep(normMode, gh+"nan/", "r")
565                             N = N + Nnan
566                         end
567                         if grep(normMode, gh+"inf/", "r")
568                             N = N + Ninf
569                             // Inf already in bins must not be counted twice:
570                             binsinf = isinf(bins)
571                             if or(binsinf)
572                                 N = N - sum(heights(binsinf))
573                             end
574                         end
575                         if type(data0)==5 & and(bins~=0) & ..
576                            grep(normMode, gh+"zeros/", "r")~=[]
577                             N = N + Nzeros
578                         end
579                     end
580                     if type(data0)==10 & ~or(bins=="") & ..
581                         grep(normMode, gh+"empty/", "r")~=[]
582                         N = N + Nempty
583                     end
584                 end
585             end
586         //else    // all (but zeros for sparses)
587         //    N = size(data,"*")
588         end
589         heights = heights / N
590         jokers = jokers / N
591     end
592     if grep(histmode, "density") then
593         heights = heights ./ (bins(2:$) - bins(1:($-1)))
594     end
595
596     // Formatting outputs for sparse input data
597     // ----------------------------------------
598     if type(data0)==5 then
599         ind = sparse(IJ, ind, size(data0))
600     end
601
602     // Shifting argouts in Scilab 5.5 compatibility mode
603     // -------------------------------------------------
604     if mode55 then
605         jokers = ind
606     end
607 endfunction