feat(analyses): compute analyses in transformed scale space

CHANGELOG.md

@@ -2,6 +2,9 @@
 
 ## Unreleased
 
+- Analyses (`linear`, `smooth`, `density`, `histogram`, `expectation`, `filled_contours`) now fit in transformed space and back-transform their output when the relevant aesthetic carries a `scale` function set via `scales` (e.g. `scales(Y = (; scale = log10))`), so a fit on log-scaled data is computed in log space. This is distinct from `axis = (; yscale = ...)`, which only transforms the display [#773](https://github.com/MakieOrg/AlgebraOfGraphics.jl/pull/773).
+- `histogram` now accepts a `direction` (`:x` or `:y`) keyword for 1D histograms, which sets the bar orientation and the aesthetic (and thus scale space) the bins are computed in [#773](https://github.com/MakieOrg/AlgebraOfGraphics.jl/pull/773).
+
 ## v0.12.13 - 2026-06-16
 
 - `draw!` now honors `facet = (; size = FacetSize(...))` instead of erroring [#770](https://github.com/MakieOrg/AlgebraOfGraphics.jl/pull/770).

docs/src/reference/draw.md

@@ -598,6 +598,27 @@ draw(
 )
 ```
 
+Unlike `ticks` and `tickformat`, the `scale` function also affects [analyses](@ref "Analyses"). When a `scale` is set on an aesthetic, analyses like `linear`, `smooth`, `density` and `histogram` fit in the transformed space and back-transform their output, so the fit matches the scaled display. This differs from setting `axis = (; yscale = log10)`, which is a display-only Makie axis attribute applied after the analysis has already run in data space. In the left panel below the line is fit in linear space and curves away from the log-linear data, in the right panel it is fit in log space and tracks it:
+
+```@example
+using AlgebraOfGraphics
+using CairoMakie
+
+t = repeat(0.0:1:8, inner = 2)
+conc = @. 100 * 10^(-0.09 * t)
+base = data((; t, conc)) * mapping(:t => "time (h)", :conc => "concentration (mg/L)")
+spec = base * visual(Scatter) + base * linear(interval = nothing) * visual(color = :firebrick)
+
+fig = Figure(size = (800, 350))
+draw!(fig[1, 1], spec; axis = (; yscale = log10, title = "axis = (; yscale = log10)"))
+draw!(fig[1, 2], spec, scales(Y = (; scale = log10)); axis = (; title = "scales(Y = (; scale = log10))"))
+fig
+```
+
+The scale function must have an inverse registered with `Makie.inverse_transform` (e.g. `log10`, `log2`, `log`, `sqrt`) so analyses can back-transform their fit.
+
+Only statistics computed by AlgebraOfGraphics analyses participate in this. A statistic computed inside a Makie recipe (e.g. `visual(Density)` or `visual(Hist)`, which run their own kernel density estimate or binning on the data AlgebraOfGraphics hands them) does not see the scale and is computed in data space, only its display is transformed. To compute in the scaled space, use the AlgebraOfGraphics analysis (`AlgebraOfGraphics.density()`, `histogram()`, ...) instead of the recipe.
+
 
 ## Legend options
 

src/algebra/layer.jl

@@ -135,6 +135,8 @@ Base.@kwdef struct ProcessedLayer <: AbstractDrawable
     attributes::NamedArguments = NamedArguments()
     scale_mapping::Dictionary{KeyType, Symbol} = Dictionary{KeyType, Symbol}() # maps mapping entries to scale ids for use of additional scales
     label::Union{Nothing, Symbol} = nothing # for selective styling via `visual(target(label))`
+    axis_transforms::Dictionary{Type{<:Aesthetic}, Base.Callable} = Dictionary{Type{<:Aesthetic}, Base.Callable}() # aesthetic type => forward scale function so analyses can fit in scale space (inverse via `Makie.inverse_transform`)
+    scale_assumed_aes::Dictionary{Int, Type{<:Aesthetic}} = Dictionary{Int, Type{<:Aesthetic}}() # positional => aesthetic the analysis assumed when fitting in scale space, checked against the resolved layer
 end
 
 function ProcessedLayer(processedlayer::ProcessedLayer; kwargs...)
@@ -147,6 +149,8 @@ function ProcessedLayer(processedlayer::ProcessedLayer; kwargs...)
         processedlayer.attributes,
         processedlayer.scale_mapping,
         processedlayer.label,
+        processedlayer.axis_transforms,
+        processedlayer.scale_assumed_aes,
     )
     return ProcessedLayer(; merge(nt, values(kwargs))...)
 end

src/algebra/layers.jl

@@ -49,9 +49,9 @@ struct ProcessedLayers <: AbstractDrawable
     layers::Vector{ProcessedLayer}
 end
 
-function ProcessedLayers(a::AbstractAlgebraic)
+function ProcessedLayers(a::AbstractAlgebraic, axis_transforms = Dictionary{Type{<:Aesthetic}, Base.Callable}())
     layers::Layers = a
-    processedlayers_array = map(process, layers)
+    processedlayers_array = map(layer -> process(layer, axis_transforms), layers)
     return ProcessedLayers(reduce(vcat, [processedlayer.layers for processedlayer in processedlayers_array]))
 end
 
@@ -70,6 +70,9 @@ ProcessedLayers(layer::Layer) = process(layer)
 ProcessedLayers(p::ProcessedLayer) = ProcessedLayers([p])
 ProcessedLayers(p::ProcessedLayers) = p
 
+ProcessedLayers(p::ProcessedLayer, ::Any) = ProcessedLayers([p])
+ProcessedLayers(p::ProcessedLayers, ::Any) = p
+
 function compute_processedlayers_grid(processedlayers, categoricalscales)
     gridpositions = compute_grid_positions(categoricalscales)
     indices = CartesianIndices((maximum(first, gridpositions), maximum(last, gridpositions)))
@@ -89,6 +92,7 @@ function compute_entries_continuousscales(pls_grid, categoricalscales, scale_pro
     continuousscales_grid = map(_ -> MultiAesScaleDict{ContinuousScale}(), pls_grid)
 
     for idx in eachindex(pls_grid), pl in pls_grid[idx]
+        check_scale_aesthetics(pl)
         # Apply continuous transformations
         positional = map(contextfree_rescale, pl.positional)
         named = map(contextfree_rescale, pl.named)
@@ -214,6 +218,68 @@ function compute_entries_continuousscales(pls_grid, categoricalscales, scale_pro
     return entries_grid, continuousscales_grid, merged_continuousscales
 end
 
+function axis_transforms_from_scales(scales::Scales)
+    transforms = Dictionary{Type{<:Aesthetic}, Base.Callable}()
+    for (sym, aes) in ((:X, AesX), (:Y, AesY), (:Z, AesZ))
+        haskey(scales.dict, sym) || continue
+        sub = scales.dict[sym]
+        haskey(sub, :scale) || continue
+        forward = sub[:scale]
+        forward === identity && continue
+        Makie.inverse_transform(forward) === nothing && error("Scale function `$forward` set for aesthetic `$sym` has no inverse registered with `Makie.inverse_transform`, so analyses cannot back-transform their fit into data space. Use a scale with a known inverse (e.g. `log10`, `log2`, `log`, `sqrt`).")
+        insert!(transforms, aes, forward)
+    end
+    return transforms
+end
+
+positional_transform(transforms, ::Nothing) = nothing
+positional_transform(transforms, aes::DataType) = haskey(transforms, aes) ? (aes => transforms[aes]) : nothing
+
+# The positional-index => aesthetic mapping a layer resolves to, via the same `aesthetic_mapping` the final
+# scale resolution uses. Computed with continuous scitypes of the given arity, which is what determines the
+# position aesthetics; empty if the plot type has no mapping for that arity. This is the single source for
+# both picking a scale space and checking that the resolved layer still agrees with it.
+function position_aesthetics(plottype, attributes, arity::Int)
+    scitypes = ScientificType[Continuous() for _ in 1:arity]
+    am = aesthetic_mapping(plottype, merge(mandatory_attributes(plottype), attributes), scitypes)
+    out = Dictionary{Int, Type{<:Aesthetic}}()
+    for (k, v) in pairs(am)
+        k isa Int && v isa Type{<:Aesthetic} && insert!(out, k, v)
+    end
+    return out
+end
+
+output_aesthetic(plottype, attributes, arity::Int, idx::Int) =
+    get(position_aesthetics(plottype, attributes, arity), idx, nothing)
+
+position_transform(transforms, plottype, attributes, arity::Int, idx::Int) =
+    isempty(transforms) ? nothing : positional_transform(transforms, output_aesthetic(plottype, attributes, arity, idx))
+
+tag_scale_aesthetics(pl::ProcessedLayer, fit::Bool) =
+    fit ? ProcessedLayer(pl; scale_assumed_aes = position_aesthetics(pl.plottype, pl.attributes, length(pl.positional))) : pl
+tag_scale_aesthetics(pls::ProcessedLayers, fit::Bool) =
+    fit ? ProcessedLayers(map(pl -> tag_scale_aesthetics(pl, true), pls.layers)) : pls
+
+# A downstream transformation (a `visual` switching plot type or orientation) can resolve a positional to a
+# different aesthetic than the analysis assumed when it chose its scale space. Compare the assumed mapping to
+# the resolved one and error if a positional moved on or off a scaled aesthetic, since the fit and the display
+# would then disagree. Aesthetic-preserving changes (color, marker, a plot type with the same position mapping)
+# don't trip it.
+function check_scale_aesthetics(pl::ProcessedLayer)
+    isempty(pl.scale_assumed_aes) && return
+    resolved = position_aesthetics(pl.plottype, pl.attributes, length(pl.positional))
+    scaled = keys(pl.axis_transforms)
+    for (idx, assumed) in pairs(pl.scale_assumed_aes)
+        actual = get(resolved, idx, nothing)
+        if actual !== assumed && (assumed in scaled || (actual !== nothing && actual in scaled))
+            error(
+                "A downstream transformation (e.g. a `visual` changing plot type or orientation) altered the aesthetic mapping of an analysis layer that was fit in a transformed scale space: positional $idx was fit assuming aesthetic `$(nameof(assumed))` but the layer resolves it to `$(actual === nothing ? "none" : nameof(actual))`, so the fit no longer matches the display. Set the plot type and orientation on the analysis itself instead of via a downstream `visual`."
+            )
+        end
+    end
+    return
+end
+
 function aesthetic_for_symbol(s::Symbol)
     aessym = Symbol("Aes", s)
     t = isdefined(AlgebraOfGraphics, aessym) ? getproperty(AlgebraOfGraphics, aessym) : nothing
@@ -377,7 +443,7 @@ function compute_axes_grid(d::AbstractDrawable, scales::Scales = scales(); axis
     axis_dict::NamedArguments = axis isa NamedArguments ? copy(axis) : dictionary(pairs(axis))
     facet_size = get(facet, :size, nothing)
 
-    processedlayers = ProcessedLayers(d).layers
+    processedlayers = ProcessedLayers(d, axis_transforms_from_scales(scales)).layers
 
     scale_props = compute_scale_properties(processedlayers, scales)
 

src/algebra/processing.jl

@@ -166,7 +166,10 @@ function process_mappings(layer::Layer)
     return ProcessedLayer(; primary, positional, named, labels, scale_mapping)
 end
 
-function process(layer::Layer)
+analysis_typename(f::ComposedFunction) = analysis_typename(f.inner)
+analysis_typename(f) = nameof(typeof(f))
+
+function process(layer::Layer, axis_transforms = Dictionary{Type{<:Aesthetic}, Base.Callable}())
     processedlayer = process_mappings(layer)
     grouped_entry = if layer.data === Pregrouped()
         # For pregrouped data, apply shiftdims to match the structure of grouped data
@@ -181,7 +184,15 @@ function process(layer::Layer)
         group(processedlayer)
     end
     primary = map(vs -> map(getuniquevalue, vs), grouped_entry.primary)
-    transformed_processlayers = ProcessedLayers(layer.transformation(ProcessedLayer(grouped_entry; primary)))
+    transformed = try
+        layer.transformation(ProcessedLayer(grouped_entry; primary, axis_transforms))
+    catch e
+        (e isa ScaleDomainError && e.analysis === nothing) ? e : rethrow()
+    end
+    if transformed isa ScaleDomainError
+        throw(ScaleDomainError(transformed.forward, transformed.value, analysis_typename(layer.transformation), transformed.aes))
+    end
+    transformed_processlayers = ProcessedLayers(transformed)
     return ProcessedLayers(
         map(transformed_processlayers.layers) do transformed_processlayer
             attributes = merge(mandatory_attributes(transformed_processlayer.plottype), transformed_processlayer.attributes)

src/scales.jl

@@ -616,6 +616,58 @@ function from_unitless_numerical(x̂::AbstractVector{<:Real}, ::AbstractVector{<
 end
 from_unitless_numerical(x̂, ::AbstractVector) = x̂
 
+apply_scale_forward(::Nothing, xn) = xn
+function apply_scale_forward(t::Pair, xn)
+    forward = last(t)
+    x̂n = try
+        forward.(xn)
+    catch e
+        e isa DomainError ? _scale_domain_error(t, e.val) : rethrow()
+    end
+    bad = findfirst(i -> isfinite(xn[i]) && !isfinite(x̂n[i]), eachindex(xn))
+    bad === nothing || _scale_domain_error(t, xn[bad])
+    return x̂n
+end
+apply_scale_inverse(::Nothing, x̂n) = x̂n
+apply_scale_inverse(t::Pair, x̂n) = Makie.inverse_transform(last(t)).(x̂n)
+
+struct ScaleDomainError <: Exception
+    forward::Any
+    value::Any
+    analysis::Union{Nothing, Symbol}
+    aes::Union{Nothing, Type}
+end
+ScaleDomainError(forward, value) = ScaleDomainError(forward, value, nothing, nothing)
+
+function Base.showerror(io::IO, e::ScaleDomainError)
+    on_aes = e.aes === nothing ? "" : " set for aesthetic `$(aesname(e.aes))`"
+    fit_by = e.analysis === nothing ? "Analyses fit" : "`$(e.analysis)` fits"
+    hint = e.aes === nothing ? "a display-only axis scale attribute (e.g. `axis = (; yscale = log10)`)" :
+        "`axis = (; $(lowercase(aesname(e.aes)))scale = $(e.forward))`"
+    return print(
+        io,
+        "The scale function `$(e.forward)`$on_aes is not finite at the data value `$(e.value)`. ",
+        "$fit_by in transformed scale space, so every value mapped onto a scale-transformed aesthetic must lie in the function's domain (e.g. strictly positive for `log`). ",
+        "Filter the offending rows, or use $hint instead, which transforms only the display and leaves the analysis in data space.",
+    )
+end
+
+_scale_domain_error(t::Pair, value) = throw(ScaleDomainError(last(t), value, nothing, first(t)))
+
+to_transformed_numerical(v, t) = apply_scale_forward(t, to_unitless_numerical(v))
+from_transformed_numerical(v̂, ref, t) = from_unitless_numerical(apply_scale_inverse(t, v̂), ref)
+
+to_transformed_nested(::Nothing, v) = map(to_unitless_numerical, v)
+to_transformed_nested(t::Pair, v) = map(g -> apply_scale_forward(t, to_unitless_numerical(g)), v)
+
+forward_datalimits(dl, ts) = dl
+forward_datalimits(dl::Tuple{Real, Real}, ts) =
+    (only(apply_scale_forward(first(ts), [dl[1]])), only(apply_scale_forward(first(ts), [dl[2]])))
+forward_datalimits(dl::Tuple, ts) = ntuple(length(dl)) do i
+    lo, hi = dl[i]
+    (only(apply_scale_forward(ts[i], [lo])), only(apply_scale_forward(ts[i], [hi])))
+end
+
 # `datalimits` accepts a single `(lo, hi)` (broadcast to every dim) or a tuple-of-pairs (one per dim);
 # strip units from the scalar limit values without disturbing that shape.
 _strip_datalimits_units(dl) = dl

src/transformations/density.jl

@@ -16,34 +16,43 @@ applydatalimits(f::Function, d) = map(f, d)
 applydatalimits(limits::Tuple{Real, Real}, d) = map(_ -> limits, d)
 applydatalimits(limits::Tuple, _) = limits
 
-function _density(vs::Tuple, vs_orig::Tuple; datalimits, npoints, kwargs...)
+function _density(vs::Tuple, vs_orig::Tuple, transforms::Tuple; datalimits, npoints, kwargs...)
     k = _kde(vs; kwargs...)
     intervals = applydatalimits(datalimits, vs)
     rgs = map(intervals) do (min, max)
         return range(min, max; length = npoints)
     end
     res = pdf(k, rgs...)
     rgs_reapplied = ntuple(length(vs_orig)) do i
-        from_unitless_numerical(collect(rgs[i]), vs_orig[i])
+        from_transformed_numerical(collect(rgs[i]), vs_orig[i], transforms[i])
     end
     return (rgs_reapplied..., res)
 end
 
 function (d::DensityAnalysis)(input::ProcessedLayer)
+    N = length(input.positional)
+    direction = N == 1 ? (d.direction === automatic ? :x : d.direction) : nothing
+    nd_plottype = N == 1 ? nothing : Makie.plottype(input.plottype, [Heatmap, Volume][N - 1])
+    dimtransforms = ntuple(N) do i
+        if N == 1
+            positional_transform(input.axis_transforms, direction === :y ? AesY : AesX)
+        else
+            position_transform(input.axis_transforms, nd_plottype, input.attributes, N + 1, i)
+        end
+    end
+    scales_active = !isempty(input.axis_transforms)
     datalimits = if d.datalimits === automatic
-        defaultdatalimits(map(v -> map(to_unitless_numerical, v), input.positional))
+        defaultdatalimits(ntuple(i -> to_transformed_nested(dimtransforms[i], input.positional[i]), N))
     else
-        _strip_datalimits_units(d.datalimits)
+        forward_datalimits(_strip_datalimits_units(d.datalimits), dimtransforms)
     end
     options = valid_options(; datalimits, d.npoints, d.kernel, d.bandwidth)
     output = map(input) do p, n
         p, n = _drop_missing_nan_rows(p, n)
-        pn = map(to_unitless_numerical, p)
-        return _density(Tuple(pn), Tuple(p); pairs(n)..., pairs(options)...), (;)
+        pn = ntuple(i -> to_transformed_numerical(p[i], dimtransforms[i]), length(p))
+        return _density(Tuple(pn), Tuple(p), dimtransforms; pairs(n)..., pairs(options)...), (;)
     end
-    N = length(input.positional)
     if N == 1
-        direction = d.direction === automatic ? :x : d.direction
         labels_base = set(input.labels, N + 1 => "pdf")
         # When direction is :y, Lines reverses the positional arguments, so we need to swap labels 1 and 2
         labels_lines = if direction === :y
@@ -73,13 +82,12 @@ function (d::DensityAnalysis)(input::ProcessedLayer)
                 (p[1], zero(p[2]), p[2]), n
             end; plottype = Band, label = :area, attributes = dictionary([:alpha => 0.15, :direction => direction]), labels = labels_base
         )
-        return ProcessedLayers([bandlayer, linelayer])
+        return tag_scale_aesthetics(ProcessedLayers([bandlayer, linelayer]), scales_active)
     else
         d.direction === automatic || error("The direction = $(repr(d.direction)) keyword in a density analysis may only be set for the 1-dimensional case")
         labels = set(input.labels, N + 1 => "pdf")
-        default_plottype = [Heatmap, Volume][N - 1]
-        plottype = Makie.plottype(input.plottype, default_plottype)
-        return ProcessedLayer(output; plottype, labels)
+        plottype = nd_plottype
+        return tag_scale_aesthetics(ProcessedLayer(output; plottype, labels), scales_active)
     end
 end
 

src/transformations/expectation.jl

@@ -11,20 +11,23 @@ function (e::ExpectationAnalysis)(input::ProcessedLayer)
     # Strip units from the value column before reducing; the Mean aggregator's
     # `(0, 0.0)` init isn't dimensionally compatible with unit-bearing values.
     # Units are reapplied to the per-group means afterwards.
+    N = length(input.positional)
+    plottype = Makie.plottype(input.plottype, categoricalplottypes[N - 1])
+    t = position_transform(input.axis_transforms, plottype, input.attributes, N, N)
     y_originals = last(input.positional)
     input_filtered = map(input) do p, n
         return _drop_missing_nan_rows(p, n)
     end
     positional_stripped = copy(input_filtered.positional)
-    positional_stripped[end] = map(to_unitless_numerical, positional_stripped[end])
+    positional_stripped[end] = to_transformed_nested(t, positional_stripped[end])
     input_stripped = ProcessedLayer(input_filtered; positional = positional_stripped)
 
     reduced = groupreduce(Mean, input_stripped)
     positional = copy(reduced.positional)
     positional[end] = map(positional[end], y_originals) do m, y_orig
-        from_unitless_numerical(m, y_orig)
+        from_transformed_numerical(m, y_orig, t)
     end
-    return ProcessedLayer(reduced; positional)
+    return tag_scale_aesthetics(ProcessedLayer(reduced; positional), !isempty(input.axis_transforms))
 end
 
 """