Return automatically chosen histogram edges as UniformEdges and fix silently dropped observations

[andreasnoack]

Return automatically chosen histogram edges as UniformEdges and fix silently dropped observations

Fixes #1009. Includes #1010, which can be closed if this is merged.

Problem

fit(Histogram, v; nbins) could silently drop observations at the extremes of the data (#1009). Two independent causes:

• For Float32 and Float16 data, histrange did its endpoint checks in the element type while returning Float64 edges. Float32(0.7) is 0.69999998 as a Float64, so an edge at 0.7 passed the check in Float32 and then excluded the observation in binindex.
• For Float64 data with closed=:right, the last edge was checked as (start + (len-1)*step)/divisor but the returned TwicePrecision range evaluated it one ulp lower, e.g. 0.19999999999999998 for data [0.0, 0.2] with 9 bins.

Both come from the same source: the edges were represented as a float range whose elements are computed by arithmetic that does not reproduce the decimal numbers the edges are meant to be, and the checks were done against something other than what binindex compares against.

Approach

The automatically chosen edges are now the decimal numbers k * 10^e for consecutive multiples k of a "nice" width (1, 2 or 5 times a power of ten, as before), each rounded to the nearest value of the data's float type F, and stored in a Vector{F}. There are only two representations involved, integers on the decimal side and F on the binary side, with one conversion per edge and no arithmetic on the converted values.

The edges are returned as UniformEdges{F} <: AbstractVector{F}, which wraps the stored vector together with the bin width. The type records what histrange knows and a plain vector would lose, namely that the bins have equal width and what that width is: step(edges) returns it, binvolume and normalize use it so bin volumes are exact rather than differences of rounded edges, show prints the edges as first edge, width, last edge like a range, and binindex estimates the bin from the width and corrects against the stored edges, which is faster than both the previous range arithmetic and a binary search. The elements are the stored F values, so no arithmetic progression is represented and the numerics do not depend on the type. User-supplied edges are untouched, so a non-range vector of edges continues to mean bins of possibly unequal width.

• Rounding each edge individually means no arithmetic progression has to be represented in F, so the same code works for Float16, Float32, Float64, BigFloat and any other AbstractFloat.
• An observation that is the rounding of the same decimal as an edge compares equal to that edge, so decimal-rounded data such as round.(x, digits=2) lands in the bin one expects: 0.7f0 is in the bin starting at 0.7f0.
• The endpoints are adjusted by comparing lo and hi against the F edges themselves, i.e. against exactly what binindex sees, so containment holds by construction for every type.
• The width is clamped to three times the floating-point spacing of the data, which guarantees strictly increasing edges (consecutive decimals then round to distinct floats, also where the last edge lands in the next binade). Data spanning only a few ulps get correspondingly few bins.
• Identical values get a single bin of unit width with decimal edges, e.g. [1.0, 2.0] for data equal to 1.05, instead of [1.05, 2.05].

The conversion of k * 10^e to F is a single IEEE multiplication or division when k and 10^|e| are exactly representable, which is correctly rounded and covers everything up to about 1e±22 for Float64 and 1e±10 for Float32. Beyond that a correctly rounded power of ten is taken from a small table and at most three roundings occur, so edges are within two ulps of the decimal. Exact conversion for all exponents is the job of a decimal parser and is not attempted here. Note that F(k // 10^d) cannot be used: Base performs that division in F with rounded operands (JuliaLang/julia#49749).

binindex now compares with < instead of isless (#1010), which treats -0.0 and 0.0 as equal without the _normalize_zero workaround and lets NaN fall outside the edges as before.

Behaviour changes

• histrange, and therefore edges of histograms fitted with nbins, returns UniformEdges{F} instead of a StepRangeLen. It is an AbstractVector{F} with step, but not an AbstractRange. For Float32 and Float16 data the element type is now F rather than Float64.
• All observations are inside the automatically chosen edges, for every float type.
• Identical values produce decimal edges around the value rather than edges at the value.
• -0.0 in user-supplied ranges is accepted instead of throwing.
• show(::Histogram) prints user-supplied edge vectors with :limit => true; UniformEdges print as first edge, width, last edge.
• midpoints of UniformEdges are averages of neighbouring edges and may differ from the previous range-based midpoints by an ulp.
• fit with nbins is faster: on 1e6 Float64 observations, 13.3 → 10.0 ms with 10 bins and 13.2 → 10.0 ms with 100 bins. histrange itself is about 20x faster (0.5 µs for 100 bins).

The docstring of fit now states that observations outside supplied edges are not counted, which has been the behaviour since 2014 but was undocumented, and that edges should have the element type of the data for decimal-rounded observations and edges to compare as intended.

Downstream packages

Checked by reading their sources. Consumers that call step(h.edges[1]) on nbins histograms keep working because UniformEdges supports step: PairPlots, HistTools, EvoId, Octofitter's legacy Makie extension. Would break:

• RadiationSpectra (@assert isa(h.edges[1], AbstractRange) on histograms passed to its plot recipes)

Ulp-level output changes only: AlgebraOfGraphics (its midpoints(::AbstractVector) method is now taken; its tests use ≈), Makie's datashader, UnicodePlots labels for Float32 data. Transparent: Plots, StatsPlots, Makie hist/stephist, BAT, FHist, LegendSpecFits. Gadfly and the Plotly packages do not use StatsBase.Histogram.

Before merging, RadiationSpectra should be fixed to not require an AbstractRange, and this should be released as a minor version.

Tests

• Regression tests for Data loss with floating-point Histograms #1009 over Float16, Float32, Float64 and BigFloat, both closed values and 1 to 12 bins: edges have element type F, contain the data, and sum(weights) equals the number of observations.
• Element type, step and exact edges for Float32, Float16 and BigFloat input; strictly increasing edges for data spanning 0 to 40 ulps at several magnitudes and up to 1000 requested bins; the two-ulp bound at extreme magnitudes; non-finite data throws; binindex on UniformEdges agrees with the generic search for all inputs including signed zeros and infinities, and NaN is outside the bins on both paths; binvolume equals step; show of UniformEdges and of long user-supplied edge vectors.
• The Fix OOM/hang in histrange when bin width is below floating-point res #1004 test asserting first(r) == x for identical values is replaced by containment and strictness checks, since the edge coinciding with the value was an artifact of the old implementation.
• Existing tests comparing histrange to ranges with == are unchanged and pass, as nearest-rounded edges coincide with the previous values.

🤖 Generated with Claude Code

[nalimilan]

I'm not really able to check the core logic, but here are a few comments.

[nalimilan]

Any idea why Int64 is used here? Same below.

[nalimilan]

Suggested change

	@test StatsBase.histrange(parse.(BigFloat, ["0.7", "0.8"]), 12, :left) ==
	parse.(BigFloat, ["0.7", "0.71", "0.72", "0.73", "0.74", "0.75", "0.76", "0.77", "0.78", "0.79", "0.8", "0.81"])
	@test StatsBase.histrange(BigFloat.(["0.7", "0.8"]), 12, :left) ==
	BigFloat.(["0.7", "0.71", "0.72", "0.73", "0.74", "0.75", "0.76", "0.77", "0.78", "0.79", "0.8", "0.81"])

[nalimilan]

Also test that first(r) < x && nextfloat(x, k) <= last(r)` as above?

[nalimilan]

Test length to check the "fewer" part of the comment?

[nalimilan]

Would there be a way to check whether Int would be enough for custom float types instead of using BigInt? Maybe by comparing maxintfloat to typemax(Int)?

[nalimilan]

Suggested change

	const _POW10_FLOAT64 = Float64[Float64(big(10)^d) for d in 0:308]
	const _POW10_FLOAT32 = Float32[Float32(big(10)^d) for d in 0:38]
	const _POW10_FLOAT16 = Float16[Float16(big(10)^d) for d in 0:4]
	const _POW10_FLOAT64 = [Float64(big(10)^d) for d in 0:308]
	const _POW10_FLOAT32 = [Float32(big(10)^d) for d in 0:38]
	const _POW10_FLOAT16 = [Float16(big(10)^d) for d in 0:4]

[nalimilan]

Shouldn't this be

Suggested change

	_pow10(::Type{F}, d::Integer) where F<:AbstractFloat = F(10)^d
	_pow10(::Type{F}, d::Integer) where F<:AbstractFloat = F(big(10)^d)

[nalimilan]

This reference makes the docs fail. You probably need to include this somewhere in the manual.