From 4c0ca07e095fa74f9ddb4146c3a645aa0a187396 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Sun, 22 Sep 2024 09:28:12 -0600 Subject: [PATCH] Formatting and tweaks --- docs/src/versions.md | 1 + ext/PlottingExt.jl | 128 ++++++++--- src/Binning/Binning.jl | 421 +++++++++++++++++----------------- src/Binning/ExperimentData.jl | 161 +++++++------ src/Sunny.jl | 4 + test/test_binning.jl | 71 ++++++ 6 files changed, 461 insertions(+), 325 deletions(-) create mode 100644 test/test_binning.jl diff --git a/docs/src/versions.md b/docs/src/versions.md index b5071f453..75add4c10 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -9,6 +9,7 @@ polynomial order according to an error tolerance. * Rename mode `:dipole_large_S` to `:dipole_uncorrected` to emphasize that corrections are missing. +* Binning features re-introduced internally for testing. ## v0.7.2 (Sep 11, 2024) diff --git a/ext/PlottingExt.jl b/ext/PlottingExt.jl index dba3e62b3..559594a37 100644 --- a/ext/PlottingExt.jl +++ b/ext/PlottingExt.jl @@ -1010,8 +1010,8 @@ function get_unit_energy(units, into) end function colorrange_from_data(; data, saturation, sensitivity, allpositive) - cmax = Statistics.quantile(vec(maximum(data; dims=1)), saturation) - cmin = Statistics.quantile(vec(minimum(data; dims=1)), 1 - saturation) + cmax = Statistics.quantile(filter(!isnan, vec(maximum(data; dims=1))), saturation) + cmin = Statistics.quantile(filter(!isnan, vec(minimum(data; dims=1))), 1 - saturation) # The returned pair (lo, hi) should be strictly ordered, lo < hi, for use in # Makie.Colorbar @@ -1192,56 +1192,51 @@ function Sunny.plot_intensities!(panel, res::Sunny.PowderIntensities{Float64}; c return ax end -function Sunny.plot_intensities!(panel, res::Sunny.BinnedIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, title="", axisopts=NamedTuple(),divide_counts = true) - +# Axes will currently be labeled as a linear combination of crystal lattice +# vectors. See https://github.com/SunnySuite/Sunny.jl/pull/310 for details. +function Sunny.plot_intensities!(panel, res::Sunny.BinnedIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, title="", axisopts=NamedTuple(), divide_counts=true) axisopts = (; title, axisopts...) unit_energy, elabel = get_unit_energy(units, into) - data = if divide_counts - res.data ./ res.counts - else - res.data - end - - data_nonan = copy(data) - good_dats = data_nonan[isfinite.(data_nonan)] - data_nonan[.!isfinite.(data_nonan)] .= sum(good_dats) / length(good_dats) + data = divide_counts ? (res.data ./ res.counts) : res.data - colorrange_suggest = colorrange_from_data(; data=reshape(data_nonan, 1, size(data)...), saturation, sensitivity=0, allpositive) + colorrange_suggest = colorrange_from_data(; data, saturation, sensitivity=0, allpositive) colormap = @something colormap (allpositive ? :gnuplot2 : :bwr) colorrange = @something colorrange colorrange_suggest # Low-dimension cases n_dims_resolved = count(res.params.numbins .!= 1) - ax = if n_dims_resolved == 0 - # No resolved data: Just display the one value! - ax = Makie.Axis(panel[1,1];axisopts...) - Makie.text!(ax,0,0;text = "$(Sunny.number_to_simple_string(data[1];digits = 4))") - ax + if n_dims_resolved == 0 + # No resolved data: Just display the one value! + ax = Makie.Axis(panel[1,1]; axisopts...) + text = Sunny.number_to_simple_string(data[1]; digits=4) + Makie.text!(ax, 0, 0; text) + return ax elseif n_dims_resolved == 1 - # Only resolved on one axis! - x_axis = findfirst(res.params.numbins .!= 1) - ax = Makie.Axis(panel[1,1];xlabel = x_axis == 4 ? elabel : Sunny.covector_name(res.params.covectors[x_axis,:]),ylabel = "Integrated Intensity",axisopts...) - bcs = Sunny.axes_bincenters(res.params) - bcs[4] /= unit_energy - Makie.barplot!(ax,bcs[x_axis],data[:];colormap, colorrange) - ax + # Only resolved on one axis! + x_axis = findfirst(res.params.numbins .!= 1) + xlabel = (x_axis == 4) ? elabel : Sunny.covector_name(res.params.covectors[x_axis, :]) + ax = Makie.Axis(panel[1,1]; xlabel, ylabel="Integrated Intensity", axisopts...) + bcs = Sunny.axes_bincenters(res.params) + bcs[4] /= unit_energy + Makie.barplot!(ax, bcs[x_axis], data[:]; colormap, colorrange) + return ax elseif n_dims_resolved == 2 - x_axis = findfirst(res.params.numbins .!= 1) - y_axis = x_axis + findfirst(res.params.numbins[x_axis+1:end] .!= 1) - ax = Makie.Axis(panel[1,1]; - xlabel = Sunny.covector_name(res.params.covectors[x_axis,:]), - ylabel = y_axis == 4 ? elabel : Sunny.covector_name(res.params.covectors[y_axis,:]),axisopts...) - bcs = Sunny.axes_bincenters(res.params) - bcs[4] /= unit_energy - hm = Makie.heatmap!(ax,bcs[x_axis],bcs[y_axis],reshape(data,size(data,x_axis),size(data,y_axis));colormap, colorrange) - Makie.Colorbar(panel[1, 2], hm) - ax + x_axis = findfirst(res.params.numbins .!= 1) + y_axis = x_axis + findfirst(res.params.numbins[x_axis+1:end] .!= 1) + xlabel = Sunny.covector_name(res.params.covectors[x_axis, :]) + ylabel = (y_axis == 4) ? elabel : Sunny.covector_name(res.params.covectors[y_axis, :]) + ax = Makie.Axis(panel[1,1]; xlabel, ylabel, axisopts...) + bcs = Sunny.axes_bincenters(res.params) + bcs[4] /= unit_energy + data = reshape(data, size(data, x_axis), size(data, y_axis)) + hm = Makie.heatmap!(ax, bcs[x_axis], bcs[y_axis], data; colormap, colorrange) + Makie.Colorbar(panel[1, 2], hm) + return ax elseif n_dims_resolved > 2 - error("Higher-dimensional binned data visualization not yet implemented!") + error("Higher-dimensional binned data visualization not yet implemented!") end - ax end #= @@ -1292,4 +1287,61 @@ function Sunny.plot_intensities(res::Sunny.AbstractIntensities; opts...) return fig end + +function Sunny.viz_qqq_path(params::Sunny.BinningParameters; kwargs...) + f = Makie.Figure() + ax = Makie.LScene(f[1,1]; show_axis=false) + Makie.cam3d!(ax.scene; projectiontype=Makie.Orthographic) + viz_qqq_path!(ax, params; kwargs...) + + aabb_lo, aabb_hi = Sunny.binning_parameters_aabb(params) + lo = min.(0, floor.(Int64, aabb_lo)) + hi = max.(0, ceil.(Int64, aabb_hi)) + Makie.scatter!(ax, map(x -> Makie.Point3f(lo .+ x.I .- 1), CartesianIndices(ntuple(i -> 1 + hi[i] - lo[i], 3)))[:], color=:black) + global_axes = [(Makie.Point3f(-1,0,0), Makie.Point3f(1,0,0)), + (Makie.Point3f(0,-1,0), Makie.Point3f(0,1,0)), + (Makie.Point3f(0,0,-1), Makie.Point3f(0,0,1))] + Makie.linesegments!(ax, global_axes, color=:black) + Makie.text!(1.1, 0, 0; text="𝐛₁ [R.L.U.]") + Makie.text!(0, 1.1, 0; text="𝐛₂ [R.L.U.]") + Makie.text!(0, 0, 1.1; text="𝐛₃ [R.L.U.]") + display(f) + ax +end + +function viz_qqq_path!(ax, params::Sunny.BinningParameters; line_alpha=0.3, color=nothing, colorrange=nothing, bin_colors=[:red,:blue,:green], bin_line_width=0.5) + @assert iszero(params.covectors[1:3, 4]) && iszero(params.covectors[4, 1:3]) + bes = Sunny.axes_binedges(params) + M = inv(params.covectors[1:3, 1:3]) + for dir in 1:3 + ix = [2, 3, 1][dir] + iy = [3, 1, 2][dir] + iz = dir + + # The grid of q points making up the lowest side of the histogram + # along the iz direction + grid = Vector{Float64}[] + grid_sparse = Vector{Float64}[] + for i in 1:length(bes[ix]), j in 1:length(bes[iy]) + is_i_edge = (i == 1 || i == length(bes[ix])) + is_j_edge = (j == 1 || j == length(bes[iy])) + grid_point = [bes[ix][i], bes[iy][j], bes[iz][1]][invperm([ix, iy, iz])] + if is_i_edge && is_j_edge # Corner case; render special for outline + push!(grid_sparse, grid_point) + continue + end + push!(grid,grid_point) + end + offset = [0, 0, bes[iz][end] - bes[iz][1]][invperm([ix, iy, iz])] + + if !isempty(grid) + segs = map(x -> (Makie.Point3f(M * x), Makie.Point3f(M * (x .+ offset))), grid[:]) + Makie.linesegments!(ax, segs, color=bin_colors[dir], linewidth=bin_line_width, alpha=line_alpha) + end + + segs = map(x -> (Makie.Point3f(M * x), Makie.Point3f(M * (x .+ offset))), grid_sparse[:]) + Makie.linesegments!(ax, segs; color=isnothing(color) ? :black : color, linewidth=2.5, colorrange) + end +end + end diff --git a/src/Binning/Binning.jl b/src/Binning/Binning.jl index 53ccfeff3..8c672f312 100644 --- a/src/Binning/Binning.jl +++ b/src/Binning/Binning.jl @@ -1,7 +1,7 @@ """ - BinningParameters(binstart,binend,binwidth;covectors = I(4)) - BinningParameters(binstart,binend;numbins,covectors = I(4)) + BinningParameters(binstart, binend, binwidth; covectors=I(4)) + BinningParameters(binstart, binend; numbins, covectors=I(4)) Describes a 4D parallelepided histogram in a format compatible with experimental Inelasitic Neutron Scattering data. See @@ -12,8 +12,8 @@ software](https://www.mantidproject.org/), or [`load_nxs`](@ref) to load The coordinates of the histogram axes are specified by multiplication of `(q,ω)` with each row of the `covectors` matrix, with `q` given in [R.L.U.]. Since the -default `covectors` matrix is the identity matrix, the default axes are -`(qx,qy,qz,ω)` in absolute units. +default `covectors` matrix is the identity matrix, the default axes are `(qx, +qy, qz, ω)` in absolute units. The convention for the binning scheme is that: - The left edge of the first bin starts at `binstart` @@ -24,17 +24,20 @@ The convention for the binning scheme is that: A `value` can be binned by computing its bin index: +```jl coords = covectors * value - bin_ix = 1 .+ floor.(Int64,(coords .- binstart) ./ binwidth) + bin_ix = 1 .+ floor.(Int64, (coords .- binstart) ./ binwidth) +``` """ mutable struct BinningParameters - binstart::MVector{4,Float64} - binend::MVector{4,Float64} - binwidth::MVector{4,Float64} - covectors::MMatrix{4,4,Float64} + binstart :: MVector{4, Float64} + binend :: MVector{4, Float64} + binwidth :: MVector{4, Float64} + covectors :: MMatrix{4, 4, Float64} end -# TODO: Use the more efficient three-argument `div(a,b,RoundDown)` instead of `floor(a/b)` -# to implement binning. Both performance and correctness need to be checked. +# TODO: Use the more efficient three-argument `fld(a, b)` instead of +# `floor(a/b)` to implement binning. Both performance and correctness need to be +# checked. function Base.show(io::IO, ::MIME"text/plain", params::BinningParameters) printstyled(io, "Binning Parameters\n"; bold=true, color=:underline) @@ -46,140 +49,140 @@ function Base.show(io::IO, ::MIME"text/plain", params::BinningParameters) printstyled(io, @sprintf("⊡ %5d bins",nbin[k]); bold=true) end bin_edges = axes_binedges(params) - first_edges = map(x -> x[1],bin_edges) - last_edges = map(x -> x[end],bin_edges) - @printf(io," from %+.3f to %+.3f along [", first_edges[k], last_edges[k]) - print(io,covector_name(params.covectors[k,:])) + first_edges = map(x -> x[1], bin_edges) + last_edges = map(x -> x[end], bin_edges) + @printf(io, " from %+.3f to %+.3f along [", first_edges[k], last_edges[k]) + print(io, covector_name(params.covectors[k,:])) @printf(io,"] (Δ = %.3f)", params.binwidth[k]/norm(params.covectors[k,:])) println(io) end end function covector_name(cov) - str = "" - axes_names = ["x","y","z","E"] - inMiddle = false - for j = 1:4 - if cov[j] != 0. - if(inMiddle) - str *= " " - end - str *= @sprintf("%+.2f d%s",cov[j],axes_names[j]) - inMiddle = true - end - end - str + str = "" + axes_names = ["𝐚₁/2π", "𝐚₂/2π", "𝐚₃/2π", "ω"] + inMiddle = false + for j in 1:4 + if cov[j] != 0. + if(inMiddle) + str *= " " + end + str *= @sprintf("%+.2f %s", cov[j], axes_names[j]) + inMiddle = true + end + end + str end # Creates a binning scheme centered on the q_path, with the specified transverse # binning directions and bin widths. function specify_transverse_binning(q_path::QPath, first_transverse_axis, second_transverse_axis, first_transverse_binwidth, second_transverse_binwidth) - # Ensure path is non-empty and single-segment - if length(q_path.qs) < 2 - error("Q path must have at least two points to infer bin width") - end + # Ensure path is non-empty and single-segment + if length(q_path.qs) < 2 + error("Q path must have at least two points to infer bin width") + end - Δq = q_path.qs[2] - q_path.qs[1] + Δq = q_path.qs[2] - q_path.qs[1] - if !all([Δq ≈ dq for dq = diff(q_path.qs)]) - error("Q path is multi-segment or irregular!") - end + if !all(isapprox(Δq), diff(q_path.qs)) + error("Q path is multi-segment or irregular!") + end - binstart = zero(MVector{4,Float64}) - binstart[4] = -Inf # Default to integrate over all energies + binstart = zero(MVector{4, Float64}) + binstart[4] = -Inf # Default to integrate over all energies - binend = zero(MVector{4,Float64}) - binend[4] = 0 + binend = zero(MVector{4, Float64}) + binend[4] = 0 - covectors = zero(MMatrix{4,4,Float64}) - recip_directions = zeros(Float64,3,3) - recip_directions[:,1] .= Δq ./ norm(Δq) - recip_directions[:,2] .= first_transverse_axis - recip_directions[:,3] .= second_transverse_axis + covectors = zero(MMatrix{4, 4, Float64}) + recip_directions = zeros(Float64, 3, 3) + recip_directions[:, 1] .= Δq ./ norm(Δq) + recip_directions[:, 2] .= first_transverse_axis + recip_directions[:, 3] .= second_transverse_axis - if minimum(svd(recip_directions).S) < 1e-8 - error("Axes are collinear!") - end - covectors[1:3,1:3] .= inv(recip_directions) + if minimum(svd(recip_directions).S) < 1e-8 + error("Axes are collinear!") + end + covectors[1:3, 1:3] .= inv(recip_directions) - coords_start = covectors[1:3,1:3] * q_path.qs[1] - coords_end = covectors[1:3,1:3] * q_path.qs[end] + coords_start = covectors[1:3, 1:3] * q_path.qs[1] + coords_end = covectors[1:3, 1:3] * q_path.qs[end] - binwidth = zero(MVector{4,Float64}) - binwidth[1] = (covectors[1:3,1:3] * Δq)[1] - binwidth[2] = first_transverse_binwidth - binwidth[3] = second_transverse_binwidth - binwidth[4] = Inf + binwidth = zero(MVector{4, Float64}) + binwidth[1] = (covectors[1:3, 1:3] * Δq)[1] + binwidth[2] = first_transverse_binwidth + binwidth[3] = second_transverse_binwidth + binwidth[4] = Inf - binstart[1:3] .= coords_start[1:3] .- binwidth[1:3]/2 - binend[1:3] .= coords_end[1:3] + binstart[1:3] .= coords_start[1:3] .- binwidth[1:3]/2 + binend[1:3] .= coords_end[1:3] - # Check the original Q points end up in bin centers - in_bin(q) = (covectors[1:3,1:3] * q .- binstart[1:3]) ./ binwidth[1:3] - centering_error(q) = (in_bin(q) .- 0.5) .- (round.(Int64,in_bin(q) .- 0.5)) - @assert all([norm(centering_error(q_path.qs[i])) < 1e-12 for i = 1:length(q_path.qs)]) + # Check the original Q points end up in bin centers + in_bin(q) = (covectors[1:3, 1:3] * q .- binstart[1:3]) ./ binwidth[1:3] + centering_error(q) = (in_bin(q) .- 0.5) .- (round.(Int64,in_bin(q) .- 0.5)) + @assert all(norm(centering_error(q)) < 1e-12 for q in q_path.qs) - # Energy axis - covectors[4,:] .= [0,0,0,1] + # Energy axis + covectors[4,:] .= [0, 0, 0, 1] - BinningParameters(binstart,binend,binwidth,covectors) + BinningParameters(binstart, binend, binwidth, covectors) end # Creates a binning scheme centered on the q_grid, with the specified transverse # binning direction and bin width. function specify_transverse_binning(q_grid::QGrid{2}, transverse_axis, transverse_binwidth) - # Ensure grid is non-empty and single-segment - if size(q_grid.qs,1) < 2 || size(q_grid.qs,2) < 2 - error("2D Q grid must have at least two points in each direction") - end - - Δq1 = q_grid.qs[2,1] - q_grid.qs[1,1] - Δq2 = q_grid.qs[1,2] - q_grid.qs[1,1] - - if !all([Δq1 ≈ dq for dq = diff(q_grid.qs,dims=1)]) || !all([Δq2 ≈ dq for dq = diff(q_grid.qs,dims = 2)]) - error("2D Q grid is irregular!") - end - - binstart = zero(MVector{4,Float64}) - binstart[4] = -Inf # Default to integrate over all energies - - binend = zero(MVector{4,Float64}) - binend[4] = 0 - - covectors = zero(MMatrix{4,4,Float64}) - recip_directions = zeros(Float64,3,3) - recip_directions[:,1] .= Δq1 ./ norm(Δq1) - recip_directions[:,2] .= Δq2 ./ norm(Δq2) - recip_directions[:,3] .= transverse_axis - - if minimum(svd(recip_directions).S) < 1e-8 - error("Transverse axis is in-plane!") - end - covectors[1:3,1:3] .= inv(recip_directions) - - coords_start = covectors[1:3,1:3] * q_grid.qs[1] - coords_end = covectors[1:3,1:3] * q_grid.qs[end] - - first_binwidth = (covectors[1:3,1:3] * Δq1)[1] - second_binwidth = (covectors[1:3,1:3] * Δq2)[2] - - binwidth = zero(MVector{4,Float64}) - binwidth[1] = first_binwidth - binwidth[2] = second_binwidth - binwidth[3] = transverse_binwidth - binwidth[4] = Inf - - binstart[1:3] .= coords_start[1:3] .- binwidth[1:3]/2 - binend[1:3] .= coords_end[1:3] - - # Check the original Q points end up in bin centers - in_bin(q) = (covectors[1:3,1:3] * q .- binstart[1:3]) ./ binwidth[1:3] - centering_error(q) = (in_bin(q) .- 0.5) .- (round.(Int64,in_bin(q) .- 0.5)) - @assert all([norm(centering_error(q_grid.qs[i])) < 1e-8 for i = 1:length(q_grid.qs)]) - - # Energy axis - covectors[4,:] .= [0,0,0,1] - BinningParameters(binstart,binend,binwidth,covectors) + # Ensure grid is non-empty and single-segment + if size(q_grid.qs, 1) < 2 || size(q_grid.qs, 2) < 2 + error("2D Q grid must have at least two points in each direction") + end + + Δq1 = q_grid.qs[2,1] - q_grid.qs[1,1] + Δq2 = q_grid.qs[1,2] - q_grid.qs[1,1] + + if !all(isapprox(Δq1), diff(q_grid.qs, dims=1)) || !all(isapprox(Δq2), diff(q_grid.qs, dims=2)) + error("2D Q grid is irregular!") + end + + binstart = zero(MVector{4, Float64}) + binstart[4] = -Inf # Default to integrate over all energies + + binend = zero(MVector{4, Float64}) + binend[4] = 0 + + covectors = zero(MMatrix{4, 4, Float64}) + recip_directions = zeros(Float64, 3, 3) + recip_directions[:, 1] .= Δq1 ./ norm(Δq1) + recip_directions[:, 2] .= Δq2 ./ norm(Δq2) + recip_directions[:, 3] .= transverse_axis + + if minimum(svd(recip_directions).S) < 1e-8 + error("Transverse axis is in-plane!") + end + covectors[1:3, 1:3] .= inv(recip_directions) + + coords_start = covectors[1:3, 1:3] * q_grid.qs[1] + coords_end = covectors[1:3, 1:3] * q_grid.qs[end] + + first_binwidth = (covectors[1:3, 1:3] * Δq1)[1] + second_binwidth = (covectors[1:3, 1:3] * Δq2)[2] + + binwidth = zero(MVector{4,Float64}) + binwidth[1] = first_binwidth + binwidth[2] = second_binwidth + binwidth[3] = transverse_binwidth + binwidth[4] = Inf + + binstart[1:3] .= coords_start[1:3] .- binwidth[1:3]/2 + binend[1:3] .= coords_end[1:3] + + # Check the original Q points end up in bin centers + in_bin(q) = (covectors[1:3, 1:3] * q .- binstart[1:3]) ./ binwidth[1:3] + centering_error(q) = (in_bin(q) .- 0.5) .- (round.(Int64, in_bin(q) .- 0.5)) + @assert all(norm(centering_error(q)) < 1e-12 for q in q_grid.qs) + + # Energy axis + covectors[4,:] .= [0,0,0,1] + BinningParameters(binstart, binend, binwidth, covectors) end """ @@ -188,23 +191,23 @@ end Create [`BinningParameters`](@ref) which place one histogram bin centered at each possible `(q,ω)` scattering vector of the crystal. This is the finest possible binning without creating bins with zero scattering vectors in them. """ -function unit_resolution_binning_parameters(sc;negative_energies = true) - ωvals = available_energies_including_zero(sc;negative_energies) +function unit_resolution_binning_parameters(sc; negative_energies=true) + ωvals = available_energies_including_zero(sc; negative_energies) good_qs = available_wave_vectors(sc) - numbins = (size(good_qs)...,length(ωvals)) + numbins = (size(good_qs)..., length(ωvals)) # Bin centers should be at Sunny scattering vectors maxQ = 1 .- (1 ./ numbins) - min_val = (0.,0.,0.,minimum(ωvals)) - max_val = (maxQ[1],maxQ[2],maxQ[3],maximum(ωvals)) + min_val = (0., 0., 0., minimum(ωvals)) + max_val = (maxQ[1], maxQ[2], maxQ[3], maximum(ωvals)) total_size = max_val .- min_val binwidth = total_size ./ (numbins .- 1) - binstart = (0.,0.,0.,minimum(ωvals)) .- (binwidth ./ 2) - binend = (maxQ[1],maxQ[2],maxQ[3],maximum(ωvals)) # bin end is well inside of last bin + binstart = (0., 0., 0., minimum(ωvals)) .- (binwidth ./ 2) + binend = (maxQ[1], maxQ[2], maxQ[3], maximum(ωvals)) # bin end is well inside of last bin - params = BinningParameters(binstart,binend,binwidth,I(4)) + params = BinningParameters(binstart, binend, binwidth, I(4)) # Special case for when there is only one bin in a direction for i = 1:4 @@ -217,10 +220,18 @@ function unit_resolution_binning_parameters(sc;negative_energies = true) params end -Base.copy(p::BinningParameters) = BinningParameters(copy(p.binstart),copy(p.binend),copy(p.binwidth),copy(p.covectors)) +function Base.copy(p::BinningParameters) + BinningParameters(copy(p.binstart), copy(p.binend), copy(p.binwidth), copy(p.covectors)) +end # Support numbins as a (virtual) property, even though only the binwidth is stored -Base.getproperty(params::BinningParameters, sym::Symbol) = sym == :numbins ? [count_bins(params.binstart[i],params.binend[i],params.binwidth[i]) for i = 1:4] : getfield(params,sym) +function Base.getproperty(params::BinningParameters, sym::Symbol) + if sym == :numbins + [count_bins(params.binstart[i], params.binend[i], params.binwidth[i]) for i in 1:4] + else + getfield(params, sym) + end +end function Base.setproperty!(params::BinningParameters, sym::Symbol, numbins) if sym == :numbins @@ -241,69 +252,70 @@ This function defines how partial bins are handled, so it should be used prefere computing the number of bins manually. """ function count_bins(bin_start,bin_end,bin_width) - if !isfinite(bin_width) - 1 - else - ceil(Int64,(bin_end - bin_start) / bin_width) - end + if !isfinite(bin_width) + 1 + else + ceil(Int64,(bin_end - bin_start) / bin_width) + end end """ axes_bincenters(params::BinningParameters) -Returns tick marks which label the bins of the histogram described by [`BinningParameters`](@ref) by their bin centers. +Returns tick marks which label the bins of the histogram described by +[`BinningParameters`](@ref) by their bin centers. -The following alternative syntax can be used to compute bin centers for a single axis: +The following alternative syntax can be used to compute bin centers for a single +axis: - axes_bincenters(binstart,binend,binwidth) + axes_bincenters(binstart, binend, binwidth) """ -function axes_bincenters(binstart,binend,binwidth) - bincenters = Vector{AbstractVector{Float64}}(undef,0) - for k = eachindex(binstart) +function axes_bincenters(binstart, binend, binwidth) + bincenters = AbstractVector{Float64}[] + for k in eachindex(binstart) if isfinite(binwidth[k]) - first_center = binstart[k] .+ binwidth[k] ./ 2 - nbin = count_bins(binstart[k],binend[k],binwidth[k]) - push!(bincenters,range(first_center,step = binwidth[k],length = nbin)) + first_center = binstart[k] .+ binwidth[k] ./ 2 + nbin = count_bins(binstart[k], binend[k], binwidth[k]) + push!(bincenters, range(first_center, step=binwidth[k], length=nbin)) else - push!(bincenters,[binstart[k]]) + push!(bincenters, [binstart[k]]) end end bincenters end -axes_bincenters(params::BinningParameters) = axes_bincenters(params.binstart,params.binend,params.binwidth) +axes_bincenters(params::BinningParameters) = axes_bincenters(params.binstart, params.binend, params.binwidth) function axes_binedges(binstart,binend,binwidth) - binedges = Vector{AbstractVector{Float64}}(undef,0) + binedges = AbstractVector{Float64}[] for k = eachindex(binstart) if isfinite(binwidth[k]) - nbin = count_bins(binstart[k],binend[k],binwidth[k]) - push!(binedges,range(binstart[k],step = binwidth[k],length = nbin + 1)) + nbin = count_bins(binstart[k], binend[k], binwidth[k]) + push!(binedges,range(binstart[k], step=binwidth[k], length=nbin+1)) else - push!(binedges,[-Inf,Inf]) + push!(binedges, [-Inf, Inf]) end end binedges end -axes_binedges(params::BinningParameters) = axes_binedges(params.binstart,params.binend,params.binwidth) +axes_binedges(params::BinningParameters) = axes_binedges(params.binstart, params.binend, params.binwidth) # Find an axis-aligned bounding box containing the histogram function binning_parameters_aabb(params) - (; binstart, binend, covectors) = params bin_edges = axes_binedges(params) - first_edges = map(x -> x[1],bin_edges) - last_edges = map(x -> x[end],bin_edges) - bin_edges = [first_edges last_edges] - this_corner = MVector{4,Float64}(undef) - q_corners = MMatrix{4,16,Float64}(undef) - for j = 1:16 # The sixteen corners of a 4-cube - for k = 1:4 # The four axes - this_corner[k] = bin_edges[k,1 + (j >> (k-1) & 1)] + first_edges = map(x -> x[1], bin_edges) + last_edges = map(x -> x[end], bin_edges) + bin_edges = hcat(first_edges, last_edges) + this_corner = zero(MVector{4, Float64}) + q_corners = zero(MMatrix{4, 16, Float64}) + for j in 1:16 # The sixteen corners of a 4-cube + for k in 1:4 # The four axes + this_corner[k] = bin_edges[k, 1 + (j >> (k-1) & 1)] end this_corner[.!isfinite.(this_corner)] .= 0 - q_corners[:,j] = covectors \ this_corner + q_corners[:, j] = params.covectors \ this_corner end - lower_aabb_q = minimum(q_corners,dims=2)[1:3] - upper_aabb_q = maximum(q_corners,dims=2)[1:3] + lower_aabb_q = minimum(q_corners, dims=2)[1:3] + upper_aabb_q = maximum(q_corners, dims=2)[1:3] return lower_aabb_q, upper_aabb_q end @@ -323,10 +335,10 @@ function Base.show(io::IO, res::BinnedIntensities) print(io, string(typeof(res)) * " ($sz bins)") end -function binned_intensities(sc,params::BinningParameters;kT = nothing,integrated_kernel = nothing) +function binned_intensities(sc, params::BinningParameters; kT=nothing, integrated_kernel=nothing) static_mode = sc isa SampledCorrelationsStatic if !isnothing(integrated_kernel) && static_mode - error("Can't broaden if data is not energy-resolved") + error("Can't broaden if data is not energy-resolved") end # Decide on which Q points can possibly contribute (depends on geometry of @@ -335,63 +347,62 @@ function binned_intensities(sc,params::BinningParameters;kT = nothing,integrated # Round the axis-aligned bounding box *outwards* to lattice sites # SQTODO: are these bounds optimal? Ls = size((static_mode ? sc.parent : sc).data)[4:6] - lower_aabb_cell = floor.(Int64,lower_aabb_q .* Ls .+ 1) - upper_aabb_cell = ceil.(Int64,upper_aabb_q .* Ls .+ 1) - cells = CartesianIndices(Tuple(((:).(lower_aabb_cell,upper_aabb_cell))))[:] + lower_aabb_cell = floor.(Int64, lower_aabb_q .* Ls .+ 1) + upper_aabb_cell = ceil.(Int64, upper_aabb_q .* Ls .+ 1) + cells = CartesianIndices(Tuple(((:).(lower_aabb_cell, upper_aabb_cell))))[:] qpts = QPoints([Vec3((cell.I .- 1) ./ Ls) for cell = cells]) # Calculate intensity at prepared qpts. No broadening yet because # that depends on the bin layout and uses an integrated_kernel! - energies = if static_mode - [0.0] - else - sort(available_energies_including_zero(sc;negative_energies=true)) - end - res = if static_mode - intensities_static(sc, qpts) + if static_mode + energies = [0.0] else - intensities(sc, qpts; energies, kT) + energies = sort(available_energies_including_zero(sc; negative_energies=true)) + end + if static_mode + res = intensities_static(sc, qpts) + else + res = intensities(sc, qpts; energies, kT) end # Bin (and broaden) those intensities according to BinningParameters - k = MVector{3,Float64}(undef) - v = MVector{4,Float64}(undef) + k = zero(MVector{3, Float64}) + v = zero(MVector{4, Float64}) q = view(v,1:3) - coords = MVector{4,Float64}(undef) - xyztBin = MVector{4,Int64}(undef) + coords = zero(MVector{4, Float64}) + xyztBin = zero(MVector{4, Int64}) xyzBin = view(xyztBin,1:3) - (; binwidth, binstart, binend, covectors, numbins) = params + (; binwidth, binstart, covectors, numbins) = params return_type = typeof(res).parameters[1] - output_intensities = zeros(return_type,numbins...) - output_counts = zeros(Float64,numbins...) + output_intensities = zeros(return_type, numbins...) + output_counts = zeros(Float64, numbins...) # Pre-compute discrete broadening kernel from continuous one provided if !isnothing(integrated_kernel) # Upgrade to 2-argument kernel if needed integrated_kernel_edep = try - integrated_kernel(0.,0.) + integrated_kernel(0., 0.) integrated_kernel catch MethodError - (ω,Δω) -> integrated_kernel(Δω) + (ω, Δω) -> integrated_kernel(Δω) end fraction_in_bin = Vector{Vector{Float64}}(undef,length(energies)) for (iω,ω) in enumerate(energies) fraction_in_bin[iω] = Vector{Float64}(undef,numbins[4]) - for iωother = 1:numbins[4] - ci_other = CartesianIndex(xyzBin[1],xyzBin[2],xyzBin[3],iωother) + for iωother in 1:numbins[4] # Start and end points of the target bin a = binstart[4] + (iωother - 1) * binwidth[4] b = binstart[4] + iωother * binwidth[4] # P(ω picked up in bin [a,b]) = ∫ₐᵇ Kernel(ω' - ω) dω' - fraction_in_bin[iω][iωother] = integrated_kernel_edep(ω,b - ω) - integrated_kernel_edep(ω,a - ω) + fraction_in_bin[iω][iωother] = integrated_kernel_edep(ω, b - ω) - integrated_kernel_edep(ω, a - ω) end end end - for cell_ix = 1:length(cells), (iω,ω) in enumerate(energies) + for cell_ix in eachindex(cells), (iω, ω) in enumerate(energies) cell = cells[cell_ix] q .= ((cell.I .- 1) ./ Ls) # q is in R.L.U. k .= (static_mode ? sc.parent : sc).crystal.recipvecs * q @@ -401,7 +412,7 @@ function binned_intensities(sc,params::BinningParameters;kT = nothing,integrated mul!(coords,covectors,v) coords .= (coords .- binstart) ./ binwidth coords[.!isfinite.(binwidth)] .= 0 - xyztBin .= 1 .+ floor.(Int64,coords) + xyztBin .= 1 .+ floor.(Int64, coords) # Check this bin is within the 4D histogram bounds if all(xyztBin .<= numbins) && all(xyztBin .>= 1) @@ -413,23 +424,23 @@ function binned_intensities(sc,params::BinningParameters;kT = nothing,integrated end else # `Energy broadening into bins' logic # For now, only support broadening for `simple' energy axes - if covectors[4,:] == [0,0,0,1] && norm(covectors[1:3,:] * [0,0,0,1]) == 0 + if covectors[4, :] == [0, 0, 0, 1] && norm(covectors[1:3, :] * [0, 0, 0, 1]) == 0 # Check this bin is within the *spatial* 3D histogram bounds # If we are energy-broadening, then scattering vectors outside the histogram # in the energy direction need to be considered - mul!(view(coords,1:3),view(covectors,1:3,1:3), view(v,1:3)) - xyzBin .= 1 .+ floor.(Int64,(view(coords,1:3) .- view(binstart,1:3)) ./ view(binwidth,1:3)) - if all(xyzBin .<= view(numbins,1:3)) && all(xyzBin .>= 1) + mul!(view(coords, 1:3),view(covectors, 1:3, 1:3), view(v, 1:3)) + xyzBin .= 1 .+ floor.(Int64, (view(coords, 1:3) .- view(binstart, 1:3)) ./ view(binwidth, 1:3)) + if all(xyzBin .<= view(numbins, 1:3)) && all(xyzBin .>= 1) # Calculate source scattering vector intensity only once - intensity = res.data[iω,cell_ix] + intensity = res.data[iω, cell_ix] # Broaden from the source scattering vector (k,ω) to # each target bin ci_other - ci_other = CartesianIndex(xyzBin[1],xyzBin[2],xyzBin[3]) - view(output_intensities,ci_other,:) .+= fraction_in_bin[iω] .* Ref(intensity) - view(output_counts,ci_other,:) .+= fraction_in_bin[iω] + ci_other = CartesianIndex(xyzBin[1], xyzBin[2], xyzBin[3]) + view(output_intensities, ci_other, :) .+= fraction_in_bin[iω] .* Ref(intensity) + view(output_counts, ci_other, :) .+= fraction_in_bin[iω] end else error("Energy broadening not yet implemented for histograms with complicated energy axes") @@ -438,7 +449,7 @@ function binned_intensities(sc,params::BinningParameters;kT = nothing,integrated end N_bins_in_BZ = abs(det(covectors[1:3,1:3])) / prod(binwidth[1:3]) output_data = output_intensities ./ N_bins_in_BZ ./ length(energies) - BinnedIntensities(params,output_data,output_counts) + BinnedIntensities(params, output_data, output_counts) end """ @@ -446,8 +457,8 @@ end Integrate over one or more axes of the histogram by setting the number of bins in that axis to 1. Examples: - integrate_axes!(params; axes = [2,3]) - integrate_axes!(params; axes = 2) + integrate_axes!(params; axes=[2,3]) + integrate_axes!(params; axes=2) """ function integrate_axes!(params::BinningParameters;axes) for k in axes @@ -459,17 +470,15 @@ function integrate_axes!(params::BinningParameters;axes) end function energy_resolve!(params::BinningParameters,energies) - energies = sort(energies) - params.binend[4] = maximum(energies) - params.binwidth[4] = energies[2] - energies[1] - params.binstart[4] = energies[1] - params.binwidth[4]/2 - params + energies = sort(energies) + params.binend[4] = maximum(energies) + params.binwidth[4] = energies[2] - energies[1] + params.binstart[4] = energies[1] - params.binwidth[4]/2 + params end function available_energies_including_zero(x; kwargs...) - ωs = available_energies(x;kwargs...) + ωs = available_energies(x; kwargs...) # Special case due to NaN definition of instant_correlations (length(ωs) == 1 && isnan(ωs[1])) ? [0.] : ωs end - - diff --git a/src/Binning/ExperimentData.jl b/src/Binning/ExperimentData.jl index c83cc6f53..1c0c5296f 100644 --- a/src/Binning/ExperimentData.jl +++ b/src/Binning/ExperimentData.jl @@ -6,19 +6,19 @@ Generate a Mantid script which bins data according to the given !!! warning "Units" - Take care to ensure the units are correct (R.L.U. or absolute). - You may want to call `Sunny.bin_rlu_as_absolute_units!` or - `Sunny.bin_absolute_units_as_rlu!` first. + Take care to ensure the units are correct (R.L.U. or absolute). You may want + to call `Sunny.bin_rlu_as_absolute_units!` or + `Sunny.bin_absolute_units_as_rlu!` first. """ function generate_mantid_script_from_binning_parameters(params) covectorsK = params.covectors # Please call rlu_to_absolute_units! first if needed - #function bin_string(k) - #if params.numbins[k] == 1 - #return "$(params.binsstart[k]),$(params.binend[k])" - #else - #return "$(params.binsstart[k]),$(params.binend[k])" - #end - #end + # function bin_string(k) + # if params.numbins[k] == 1 + # return "$(params.binsstart[k]),$(params.binend[k])" + # else + # return "$(params.binsstart[k]),$(params.binend[k])" + # end + # end return """MakeSlice(InputWorkspace="merged_mde_INPUT", QDimension0="$(covectorsK[1,1]),$(covectorsK[1,2]),$(covectorsK[1,3])", QDimension1="$(covectorsK[2,1]),$(covectorsK[2,2]),$(covectorsK[2,3])", @@ -45,78 +45,77 @@ To load another field instead of the signal, specify e.g. `num_events`, and `signal`. """ function load_nxs(filename; field="signal") - JLD2.jldopen(filename,"r") do file + JLD2.jldopen(filename, "r") do file read_covectors_from_axes_labels = false # This variable is basically the "Mantid W-Matrix". See discussion on # Github: https://github.com/SunnySuite/Sunny.jl/pull/256. - spatial_covectors = Matrix{Float64}(undef,3,3) + spatial_covectors = zeros(3, 3) try - try - w_matrix = file["MDHistoWorkspace"]["experiment0"]["logs"]["W_MATRIX"]["value"] - - # Transpose to arrange axes labels as columns - spatial_covectors .= transpose(reshape(w_matrix,3,3)) - catch e - printstyled("Warning",color=:yellow) - print(": failed to load W_MATRIX from Mantid file $filename due to:\n") - println(e) - println("Falling back to reading transform_from_orig") - - coordinate_system = file["MDHistoWorkspace"]["coordinate_system"][1] - - # Permalink to where this enum is defined: - # https://github.com/mantidproject/mantid/blob/057df5b2de1c15b819c7dd06e50bdbf5461b09c6/Framework/Kernel/inc/MantidKernel/SpecialCoordinateSystem.h#L14 - system_name = ["None", "QLab", "QSample", "HKL"][coordinate_system + 1] - - # The matrix stored in the file is transpose of the actual - # transform_from_orig matrix - transform_from_orig = file["MDHistoWorkspace"]["transform_from_orig"] - transform_from_orig = reshape(transform_from_orig,5,5) - transform_from_orig = transform_from_orig' - - # U: Orthogonal rotation matrix - # B: inv(lattice_vectors(...)) matrix, as in Sunny - # The matrix stored in the file is transpose of U * B - ub_matrix = file["MDHistoWorkspace"]["experiment0"]["sample"]["oriented_lattice"]["orientation_matrix"] - ub_matrix = ub_matrix' - - # Following: https://docs.mantidproject.org/nightly/concepts/Lattice.html - # It can be verified that the matrix G^* = (ub_matrix' * ub_matrix) - # is equal to B' * B, where B = inv(lattice_vectors(...)), and the diagonal - # entries of the inverse of this matrix are the lattice parameters squared - # - #abcMagnitude = sqrt.(diag(inv(ub_matrix' * ub_matrix))) - #println("[a,b,c] = ",abcMagnitude) - - # This is how you extract the covectors from `transform_from_orig` and `ub_matrix` - # TODO: Parse this from the `long_name` of the data_dims instead - spatial_covectors .= 2π .* transform_from_orig[1:3,1:3] * ub_matrix - end + try + w_matrix = file["MDHistoWorkspace"]["experiment0"]["logs"]["W_MATRIX"]["value"] + # Transpose to arrange axes labels as columns + spatial_covectors .= transpose(reshape(w_matrix, 3, 3)) + catch e + printstyled("Warning", color=:yellow) + print(": failed to load W_MATRIX from Mantid file $filename due to:\n") + println(e) + println("Falling back to reading transform_from_orig") + + coordinate_system = file["MDHistoWorkspace"]["coordinate_system"][1] + + # Permalink to where this enum is defined: + # https://github.com/mantidproject/mantid/blob/057df5b2de1c15b819c7dd06e50bdbf5461b09c6/Framework/Kernel/inc/MantidKernel/SpecialCoordinateSystem.h#L14 + system_name = ["None", "QLab", "QSample", "HKL"][coordinate_system + 1] + + # The matrix stored in the file is transpose of the actual + # transform_from_orig matrix + transform_from_orig = file["MDHistoWorkspace"]["transform_from_orig"] + transform_from_orig = reshape(transform_from_orig, 5, 5) + transform_from_orig = transform_from_orig' + + # U: Orthogonal rotation matrix + # B: inv(lattice_vectors(...)) matrix, as in Sunny + # The matrix stored in the file is transpose of U * B + ub_matrix = file["MDHistoWorkspace"]["experiment0"]["sample"]["oriented_lattice"]["orientation_matrix"] + ub_matrix = ub_matrix' + + # Following: https://docs.mantidproject.org/nightly/concepts/Lattice.html + # It can be verified that the matrix G^* = (ub_matrix' * ub_matrix) + # is equal to B' * B, where B = inv(lattice_vectors(...)), and the diagonal + # entries of the inverse of this matrix are the lattice parameters squared + # + #abcMagnitude = sqrt.(diag(inv(ub_matrix' * ub_matrix))) + #println("[a,b,c] = ",abcMagnitude) + + # This is how you extract the covectors from `transform_from_orig` and `ub_matrix` + # TODO: Parse this from the `long_name` of the data_dims instead + spatial_covectors .= 2π .* transform_from_orig[1:3,1:3] * ub_matrix + end catch e - printstyled("Warning",color=:yellow) - print(": failed to read histogram axes from Mantid file $filename due to:\n") - println(e) - println("Defaulting to low-accuracy reading of axes labels!") - read_covectors_from_axes_labels = true + printstyled("Warning", color=:yellow) + print(": failed to read histogram axes from Mantid file $filename due to:\n") + println(e) + println("Defaulting to low-accuracy reading of axes labels!") + read_covectors_from_axes_labels = true end signal = file["MDHistoWorkspace"]["data"][field] - axes = Dict(JLD2.load_attributes(file,"MDHistoWorkspace/data/signal"))[:axes] + axes = Dict(JLD2.load_attributes(file, "MDHistoWorkspace/data/signal"))[:axes] # Axes are just stored backwards in Mantid .nxs for some reason axes_names = reverse(split(axes,":")) - data_dims = Vector{Vector{Float64}}(undef,length(axes_names)) - binwidth = Vector{Float64}(undef,4) - binstart = Vector{Float64}(undef,4) - binend = Vector{Float64}(undef,4) + data_dims = Vector{Vector{Float64}}(undef, length(axes_names)) + binwidth = zeros(4) + binstart = zeros(4) + binend = zeros(4) covectors = zeros(4, 4) spatial_covector_ixs = [0,0,0] std = x -> sqrt(sum((x .- sum(x) ./ length(x)).^2)) - for (i,name) in enumerate(axes_names) - long_name = Dict(JLD2.load_attributes(file,"MDHistoWorkspace/data/$name"))[:long_name] + for (i, name) in enumerate(axes_names) + long_name = Dict(JLD2.load_attributes(file, "MDHistoWorkspace/data/$name"))[:long_name] if long_name == "DeltaE" covectors[i,:] .= [0,0,0,1] # energy covector @@ -143,32 +142,32 @@ function load_nxs(filename; field="signal") covectors[spatial_covector_ixs,1:3] .= inv(spatial_covectors) - return BinningParameters(binstart,binend,binwidth,covectors), signal + return BinningParameters(binstart, binend, binwidth, covectors), signal end end -function Base.permutedims(params::BinningParameters,perm) - out = copy(params) - out.covectors .= params.covectors[perm,:] - out.binwidth .= params.binwidth[perm] - out.binstart .= params.binstart[perm] - out.binend .= params.binend[perm] - out +function Base.permutedims(params::BinningParameters, perm) + out = copy(params) + out.covectors .= params.covectors[perm,:] + out.binwidth .= params.binwidth[perm] + out.binstart .= params.binstart[perm] + out.binend .= params.binend[perm] + out end # Parse the "[0.5H,0.3H,0.1H]" type of Mantid string describing # a histogram axis function parse_long_name(long_name) - # You're welcome - found = match(r"\[(|0|(?:-?[0-9]?(?:\.[0-9]+)?))[HKL]?,(|0|(?:-?[0-9]?(?:\.[0-9]+)?))[HKL]?,(|0|(?:-?[0-9]?(?:\.[0-9]+)?))[HKL]?\]",long_name) - if isnothing(found) - return nothing - end - return map(x -> isempty(x) ? 1. : x == "-" ? -1. : parse(Float64,x),found) + # You're welcome + found = match(r"\[(|0|(?:-?[0-9]?(?:\.[0-9]+)?))[HKL]?,(|0|(?:-?[0-9]?(?:\.[0-9]+)?))[HKL]?,(|0|(?:-?[0-9]?(?:\.[0-9]+)?))[HKL]?\]", long_name) + if isnothing(found) + return nothing + end + return map(x -> isempty(x) ? 1. : x == "-" ? -1. : parse(Float64, x), found) end -function quick_view_nxs(filename,keep_ax) - integration_axes = setdiff(1:4,keep_ax) +function quick_view_nxs(filename, keep_ax) + integration_axes = setdiff(1:4, keep_ax) params, signal = load_nxs(filename) integrate_axes!(params, axes=integration_axes) int_signal = dropdims(sum(signal, dims=integration_axes); dims=Tuple(integration_axes)) diff --git a/src/Sunny.jl b/src/Sunny.jl index 1f5c13d2e..5d0b8d9a9 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -156,6 +156,10 @@ function plot_intensities(args...; opts...) end export view_crystal, plot_spins, plot_spins!, plot_intensities, plot_intensities! +function viz_qqq_path(args...; opts...) + error(isloaded("Makie") ? "Invalid method call" : "Import GLMakie to enable plotting") +end + ### ext/ExportVTKExt.jl, dependent on WriteVTK function export_vtk(args...) error(isloaded("WriteVTK") ? "Invalid method call" : "Import WriteVTK to enable exporting") diff --git a/test/test_binning.jl b/test/test_binning.jl new file mode 100644 index 000000000..151159c2d --- /dev/null +++ b/test/test_binning.jl @@ -0,0 +1,71 @@ +@testitem "Binning" begin + using LinearAlgebra + + # Setup an arbitrary qgrid + latvecs = lattice_vectors(1, 1, 10, 90, 90, 90) + cryst = Crystal(latvecs, [[0, 0, 0]]) + qgrid = q_space_grid(cryst, [0, 1.2, 0.5], range(0, 1, 100), [0, -8, 0.1], (-3,-1)) + + dqx = diff(qgrid.qs, dims=1)[1] + dqy = diff(qgrid.qs, dims=2)[1] + + @test_throws "in-plane" Sunny.specify_transverse_binning(qgrid, dqx, 0.1) + + params = Sunny.specify_transverse_binning(qgrid, dqx × dqy, 0.1) + + # There should be exactly one transverse bin + @test params.numbins[3] == 1 + + # Energy axis should be fully integrated by default + @test isinf(params.binwidth[4]) + + # Test that moving by one gridpoint in the QGrid moves by one binwidth in + # binning coordinate space + @test (params.covectors[1:3,1:3] * dqx)[1] ≈ params.binwidth[1] + @test (params.covectors[1:3,1:3] * dqy)[2] ≈ params.binwidth[2] + + # Make a system with only a few momenta available + sys = System(cryst, [1 => Moment(s=1/2, g=2)], :dipole; dims=(5,5,1)) + randomize_spins!(sys) + + # This should be added to Sunny (oops) + function Sunny.unit_resolution_binning_parameters(isc::SampledCorrelationsStatic; kwargs...) + params = Sunny.unit_resolution_binning_parameters(isc.parent; kwargs...) + # Integrate over all energies + params.binstart[4] = -Inf + params.binwidth[4] = Inf + params + end + + static_corrs = SampledCorrelationsStatic(sys; measure=ssf_trace(sys)) + add_sample!(static_corrs,sys) + + res_interp = intensities_static(static_corrs, qgrid) + res_bin = Sunny.binned_intensities(static_corrs, params) + + # The binning and interpolating results should be very different because + # only a few momenta are available. The binned data will be very sparse, + # while interpolation makes it dense + @test count(iszero.(res_bin.data)) > 10 * count(iszero.(res_interp.data)) + + # Returning to a better behaved q-grid, we should get agreement + unit_params = Sunny.unit_resolution_binning_parameters(static_corrs) + res_bin_unit = Sunny.binned_intensities(static_corrs, unit_params) + + # Construct the equivalent grid of q points + bcs = Sunny.axes_bincenters(unit_params) + unit_qgrid_qpts = Sunny.QPoints([[qx,qy,qz] for qx = bcs[1], qy = bcs[2], qz = bcs[3]][:]) + res_interp_unit = intensities_static(static_corrs,unit_qgrid_qpts) + + ratio = res_bin_unit.data[:,:,1,1] ./ reshape(res_interp_unit.data,5,5) + + # The results should be proportional + @test sum((ratio .- ratio[1]).^2) < 1e-12 + + # The constant of proportionality should be the bin volume. + @test ratio[1] ≈ prod(unit_params.binwidth[1:3]) + + # This is the expected sum rule for this setup (TODO: expression in terms of + # N, S, etc) + @test sum(res_bin_unit.data) ≈ 1 +end