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