--- /home/jwe/src/octave/scripts/plot/draw/hist.m 2019-02-13 12:23:38.436417027 +0000 +++ /home/jwe/downloads/hist.m 2019-03-13 16:58:17.689686257 +0000 @@ -104,6 +104,9 @@ min_val = double (min_val); endif + ## Equidistant entries allow much more efficient algorithms. + x_is_range = true; + ## Process possible second argument if (nargin == 1 || ischar (varargin{iarg})) n = 10; @@ -141,6 +144,7 @@ endif x = x.'; # Convert to matrix; elseif (isvector (x)) + x_is_range = isequal (typeinfo (x), "range"); x = x(:); if (! issorted (x)) warning ("hist: bin values X not sorted on input"); @@ -168,29 +172,48 @@ n = rows (x); y_nc = columns (y); - if (n < 30 && columns (x) == 1) - ## The following algorithm works fastest for n less than about 30. + + if (n < 26 * (1 + double (! x_is_range))) + ## The following algorithm works fastest for n less than about 10. + nanind=isnan (y); chist = zeros (n+1, y_nc); for i = 1:n-1 chist(i+1,:) = sum (y <= cutoff(i)); endfor - chist(n+1,:) = sum (! isnan (y)); + chist(n+1,:) = sum (! nanind); + + freq = diff (chist); else - ## The following algorithm works fastest for n greater than about 30. - ## Put cutoff elements between boundaries, integrate over all - ## elements, keep totals at boundaries. - m = (nthargout (2, @sort, [y; repmat(cutoff, 1, y_nc)]) <= rows (y)); - chist = cumsum (m); - chist = [(zeros (1, y_nc)); - (reshape (chist(! m), rows (cutoff), y_nc)); - (chist(end,:) - sum (isnan (y)))]; + ## Lookup is more efficient if y is sorted, but sorting costs. + if (! x_is_range && n > sqrt (rows (y) * 1e4)) + y = sort (y); + end + + nanind = isnan (y); + y(nanind) = 0; + freq = zeros (n, y_nc); + if (x_is_range) + if (n < 3) + d = 1; + else + d = cutoff(2) - cutoff(1); + end + for j = 1:y_nc + freq(:,j) = accumarray (1 + max (0, min (length (cutoff), ceil ((y(:,j) - cutoff(1)) / d))), + double (! nanind(:,j)), [n, 1]); + end + else + for j = 1:y_nc + i = lookup (cutoff, y(:,j)); + i = 1 + i - (cutoff(max (i, 1)) == y(:,j)); + freq(:,j) = accumarray (i, double (! nanind(:,j)), [n, 1]); + end + end endif - freq = diff (chist); - if (norm) ## Normalize the histogram - freq = freq .* norm ./ sum (! isnan (y)); + freq = freq .* norm ./ sum (! nanind); endif if (nargout == 0)