From d83b1177d6047937901a5c89a6ee638725d3d62c Mon Sep 17 00:00:00 2001 From: Sam Quinn Date: Wed, 25 Sep 2024 08:23:34 -0700 Subject: [PATCH] 0.7-style binning (#310) Co-authored-by: Kipton Barros --- ext/PlottingExt.jl | 108 +++++- src/Binning/Binning.jl | 607 ++++++++++++++++++++-------------- src/Binning/ExperimentData.jl | 164 +++++---- src/Sunny.jl | 4 + test/test_binning.jl | 69 ++++ 5 files changed, 625 insertions(+), 327 deletions(-) create mode 100644 test/test_binning.jl diff --git a/ext/PlottingExt.jl b/ext/PlottingExt.jl index 69213b864..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,6 +1192,53 @@ function Sunny.plot_intensities!(panel, res::Sunny.PowderIntensities{Float64}; c return ax end +# 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 = divide_counts ? (res.data ./ res.counts) : res.data + + 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) + + 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) + 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) + 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!") + end +end + #= * `axisopts`: An additional collection of named arguments that will be passed to the `Makie.Axis` constructor. This allows to override the axis `xlabel`, @@ -1240,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 6f76b377c..5c395a4ed 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,163 +49,167 @@ 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]) - axes_names = ["x","y","z","E"] - inMiddle = false - for j = 1:4 - if params.covectors[k,j] != 0. - if(inMiddle) - print(io," ") - end - @printf(io,"%+.2f d%s",params.covectors[k,j],axes_names[j]) - inMiddle = true - end - end + 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 -#= - -Base.copy(p::BinningParameters) = BinningParameters(copy(p.binstart),copy(p.binend),copy(p.binwidth),copy(p.covectors)) - -# 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,params.binend,params.binwidth) : getfield(params,sym) - -function Base.setproperty!(params::BinningParameters, sym::Symbol, numbins) - if sym == :numbins - # *Ensure* that the last bin contains params.binend - params.binwidth .= (params.binend .- params.binstart) ./ (numbins .- 0.5) - else - setfield!(params,sym,numbins) +function covector_name(cov) + 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 -""" - count_bins(binstart,binend,binwidth) +# 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 -Returns the number of bins in the binning scheme implied by `binstart`, `binend`, and `binwidth`. -To count the bins in a [`BinningParameters`](@ref), use `params.numbins`. + Δq = q_path.qs[2] - q_path.qs[1] -This function defines how partial bins are handled, so it should be used preferentially over -computing the number of bins manually. -""" -count_bins(bin_start,bin_end,bin_width) = ceil.(Int64,(bin_end .- bin_start) ./ bin_width) + if !all(isapprox(Δq), diff(q_path.qs)) + error("q-path is multi-segment or irregular!") + end -function BinningParameters(binstart,binend,binwidth;covectors = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]) - return BinningParameters(binstart,binend,binwidth,covectors) -end + binstart = zero(MVector{4, Float64}) + binstart[4] = -Inf # Default to integrate over all energies -function BinningParameters(binstart,binend;numbins,kwargs...) - params = BinningParameters(binstart,binend,[0.,0,0,0];kwargs...) - params.numbins = numbins # Use the setproperty to do it correctly - params -end + binend = zero(MVector{4, Float64}) + binend[4] = 0 -""" - integrate_axes!(params::BinningParameters; axes) -Integrate over one or more axes of the histogram by setting the number of bins -in that axis to 1. Examples: + 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 - integrate_axes!(params; axes = [2,3]) - integrate_axes!(params; axes = 2) -""" -function integrate_axes!(params::BinningParameters;axes) - for k in axes - nbins = [params.numbins.data...] - nbins[k] = 1 - params.numbins = SVector{4}(nbins) + if minimum(svd(recip_directions).S) < 1e-8 + error("Axes are collinear!") end - return params -end + covectors[1:3, 1:3] .= inv(recip_directions) -# 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)] - end - q_corners[:,j] = covectors \ this_corner - end - 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 + coords_start = covectors[1:3, 1:3] * q_path.qs[1] + coords_end = covectors[1:3, 1:3] * q_path.qs[end] -""" -If `params` expects to bin values `(k,ω)` in absolute units, then calling + 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 - bin_rlu_as_absolute_units!(params::BinningParameters,[reciprocal lattice vectors]) + binstart[1:3] .= coords_start[1:3] .- binwidth[1:3]/2 + binend[1:3] .= coords_end[1:3] -will modifiy the `covectors` in `params` so that they will accept `(q,ω)` in Reciprocal Lattice Units (R.L.U.) instead. -Conversly, if `params` expects `(q,ω)` R.L.U., calling + # 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) - bin_absolute_units_as_rlu!(params::BinningParameters,[reciprocal lattice vectors]) + # Energy axis + covectors[4,:] .= [0, 0, 0, 1] -will adjust `params` to instead accept `(k,ω)` absolute units. + BinningParameters(binstart, binend, binwidth, covectors) +end -The second argument may be a 3x3 matrix specifying the reciprocal lattice vectors, or a [`Crystal`](@ref). -""" -bin_absolute_units_as_rlu!, bin_rlu_as_absolute_units! +# 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 -function bin_rlu_as_absolute_units!(params::BinningParameters,recip_vecs::AbstractMatrix) - covectorsK = params.covectors + Δq1 = q_grid.qs[2,1] - q_grid.qs[1,1] + Δq2 = q_grid.qs[1,2] - q_grid.qs[1,1] - # covectorsQ * q = covectorsK * recip_vecs * q = covectorsK * k - # covectorsQ = covectorsK * recip_vecs - covectorsQ = covectorsK * [recip_vecs [0;0;0]; [0 0 0] 1] - params.covectors = MMatrix{4,4}(covectorsQ) - params -end + 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 -# covectorsK * k = covectorsQ * inv(recip_vecs) * k = covectorsQ * q -# covectorsK = covectorsQ * inv(recip_vecs) -bin_absolute_units_as_rlu!(params::BinningParameters,recip_vecs::AbstractMatrix) = bin_rlu_as_absolute_units!(params,inv(recip_vecs)) + binstart = zero(MVector{4, Float64}) + binstart[4] = -Inf # Default to integrate over all energies -bin_absolute_units_as_rlu!(params::BinningParameters,crystal::Crystal) = bin_absolute_units_as_rlu!(params,2π*inv(crystal.latvecs)') + binend = zero(MVector{4, Float64}) + binend[4] = 0 -bin_rlu_as_absolute_units!(params::BinningParameters,crystal::Crystal) = bin_absolute_units_as_rlu!(params,crystal.latvecs'/2π) + 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 -""" - unit_resolution_binning_parameters(sc::SampledCorrelations) + if minimum(svd(recip_directions).S) < 1e-8 + error("Transverse axis is in-plane!") + end + covectors[1:3, 1:3] .= inv(recip_directions) -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. + coords_start = covectors[1:3, 1:3] * q_grid.qs[1] + coords_end = covectors[1:3, 1:3] * q_grid.qs[end] -This function can be used without reference to a [`SampledCorrelations`](@ref) using an alternate syntax to manually specify the bin centers for the energy axis and the lattice size: + first_binwidth = (covectors[1:3, 1:3] * Δq1)[1] + second_binwidth = (covectors[1:3, 1:3] * Δq2)[2] - unit_resolution_binning_parameters(ω_bincenters,dims,[reciprocal lattice vectors]) + binwidth = zero(MVector{4,Float64}) + binwidth[1] = first_binwidth + binwidth[2] = second_binwidth + binwidth[3] = transverse_binwidth + binwidth[4] = Inf -The last argument may be a 3x3 matrix specifying the reciprocal lattice vectors, or a [`Crystal`](@ref). + binstart[1:3] .= coords_start[1:3] .- binwidth[1:3]/2 + binend[1:3] .= coords_end[1:3] -Lastly, binning parameters for a single axis may be specifed by their bin centers: + # 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)) + # TODO: Understand loss of precision here (1e-13 will fail tests) + @assert all(norm(centering_error(q)) < 1e-10 for q in q_grid.qs) + + # Energy axis + covectors[4,:] .= [0,0,0,1] + BinningParameters(binstart, binend, binwidth, covectors) +end - (binstart,binend,binwidth) = unit_resolution_binning_parameters(bincenters::Vector{Float64}) """ -function unit_resolution_binning_parameters(ωvals,dims,args...) - numbins = (dims...,length(ωvals)) + unit_resolution_binning_parameters(sc::SampledCorrelations) + +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::SampledCorrelations; negative_energies=true) + ωvals = available_energies_including_zero(sc; negative_energies) + + good_qs = available_wave_vectors(sc) + numbins = (size(good_qs)..., length(ωvals)) # Bin centers should be at Sunny scattering vectors - maxQ = 1 .- (1 ./ numbins) + 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) + params = BinningParameters(binstart, binend, binwidth, I(4)) # Special case for when there is only one bin in a direction for i = 1:4 @@ -215,157 +222,273 @@ function unit_resolution_binning_parameters(ωvals,dims,args...) params end -unit_resolution_binning_parameters(sc::SampledCorrelations; negative_energies=false,kwargs...) = unit_resolution_binning_parameters(available_energies_including_zero(sc;negative_energies),sc.dims,sc;kwargs...) +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 -function unit_resolution_binning_parameters(ωvals::AbstractVector{Float64}) - if !all(abs.(diff(diff(ωvals))) .< 1e-12) - @warn "Non-uniform bins will be re-spaced into uniform bins" - end - if length(ωvals) == 1 - error("Can not infer bin width given only one bin center") +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 +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 - ωbinwidth = (maximum(ωvals) - minimum(ωvals)) / (length(ωvals) - 1) - ωstart = minimum(ωvals) - ωbinwidth / 2 - ωend = maximum(ωvals) +end - return ωstart, ωend, ωbinwidth +function Base.setproperty!(params::BinningParameters, sym::Symbol, numbins) + if sym == :numbins + # *Ensure* that the last bin contains params.binend + params.binwidth .= (params.binend .- params.binstart) ./ (numbins .- 0.5) + else + setfield!(params,sym,numbins) + end end """ - slice_2D_binning_parameter(sc::SampledCorrelations, cut_from_q, cut_to_q, cut_bins::Int64, cut_width::Float64; plane_normal = [0,0,1],cut_height = cutwidth) - -Creates [`BinningParameters`](@ref) which make a cut along one dimension of Q-space. - -The x-axis of the resulting histogram consists of `cut_bins`-many bins ranging -from `cut_from_q` to `cut_to_q`. -The width of the bins in the transverse direciton is controlled by `cut_width` and `cut_height`. - -The binning in the transverse directions is defined in the following way, which sets their normalization and orthogonality properties: - - cut_covector = normalize(cut_to_q - cut_from_q) - transverse_covector = normalize(plane_normal × cut_covector) - cotransverse_covector = normalize(transverse_covector × cut_covector) - -In other words, the axes are orthonormal with respect to the Euclidean metric. - -If the cut is too narrow, there will be very few scattering vectors per bin, or -the number per bin will vary substantially along the cut. -If the output appears under-resolved, try increasing `cut_width`. - -The four axes of the resulting histogram are: - 1. Along the cut - 2. Fist transverse Q direction - 3. Second transverse Q direction - 4. Energy - -This function can be used without reference to a [`SampledCorrelations`](@ref) using this alternate syntax to manually specify the bin centers for the energy axis: + count_bins(binstart,binend,binwidth) - slice_2D_binning_parameter(ω_bincenters, cut_from, cut_to,...) +Returns the number of bins in the binning scheme implied by `binstart`, `binend`, and `binwidth`. +To count the bins in a [`BinningParameters`](@ref), use `params.numbins`. -where `ω_bincenters` specifies the energy axis, and both `cut_from` and `cut_to` are arbitrary covectors, in any units. +This function defines how partial bins are handled, so it should be used preferentially over +computing the number of bins manually. """ -function slice_2D_binning_parameters(ωvals::Vector{Float64},cut_from_q,cut_to_q,cut_bins::Int64,cut_width;plane_normal = [0,0,1],cut_height = cut_width) - # This covector should measure progress along the cut in r.l.u. - cut_covector = normalize(cut_to_q - cut_from_q) - # These two covectors should be perpendicular to the cut, and to each other - transverse_covector = normalize(plane_normal × cut_covector) - cotransverse_covector = normalize(transverse_covector × cut_covector) - - start_x = cut_covector ⋅ cut_from_q - end_x = cut_covector ⋅ cut_to_q - - transverse_center = transverse_covector ⋅ cut_from_q # Equal to using cut_to_q - cotransverse_center = cotransverse_covector ⋅ cut_from_q - - ωstart, ωend, ωbinwidth = unit_resolution_binning_parameters(ωvals) - xstart, xend, xbinwidth = unit_resolution_binning_parameters(range(start_x,end_x,length = cut_bins)) - - binstart = [xstart,transverse_center - cut_width/2,cotransverse_center - cut_height/2,ωstart] - binend = [xend,transverse_center,cotransverse_center,ωend] - numbins = [cut_bins,1,1,length(ωvals)] - covectors = [cut_covector... 0; transverse_covector... 0; cotransverse_covector... 0; 0 0 0 1] - - BinningParameters(binstart,binend;numbins = numbins, covectors = covectors) -end - -function slice_2D_binning_parameters(sc::SampledCorrelations,cut_from_q,cut_to_q,args...;kwargs...) - slice_2D_binning_parameters(available_energies_including_zero(sc),cut_from_q,cut_to_q,args...;kwargs...) +function count_bins(bin_start,bin_end,bin_width) + 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{AbstractRange{Float64}}(undef,0) - for k = eachindex(binstart) - 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)) +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)) + else + 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{AbstractRange{Float64}}(undef,0) + binedges = AbstractVector{Float64}[] for k = eachindex(binstart) - nbin = count_bins(binstart[k],binend[k],binwidth[k]) - push!(binedges,range(binstart[k],step = binwidth[k],length = nbin + 1)) + if isfinite(binwidth[k]) + 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]) + 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) + bin_edges = axes_binedges(params) + 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] = 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] + return lower_aabb_q, upper_aabb_q +end +struct BinnedIntensities{T} <: AbstractIntensities + # Original chemical cell + #crystal :: Crystal + # BinningParameters in RLU + params :: BinningParameters + # Intensity data as bin-integrated values + data :: Array{T, 4} # collect(size(data)) == params.numbins + # Number of individually binned contributions (useful for some normalizations) + counts :: Array{Float64, 4} +end -""" - reciprocal_space_path_bins(sc,qs,density,args...;kwargs...) +function Base.show(io::IO, res::BinnedIntensities) + sz = join(res.params.numbins, "×") + print(io, string(typeof(res)) * " ($sz bins)") +end + +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") + end + + # Decide on which q points can possibly contribute (depends on geometry of + # supercell and params) + lower_aabb_q, upper_aabb_q = binning_parameters_aabb(params) + # 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))))[:] + 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! + if static_mode + energies = [0.0] + else + 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 -Takes a list of wave vectors, `qs` in R.L.U., and builds a series of histogram [`BinningParameters`](@ref) -whose first axis traces a path through the provided points. -The second and third axes are integrated over according to the `args` and `kwargs`, -which are passed through to [`slice_2D_binning_parameters`](@ref). + # Bin (and broaden) those intensities according to BinningParameters + k = zero(MVector{3, Float64}) + v = zero(MVector{4, Float64}) + q = view(v,1:3) + coords = zero(MVector{4, Float64}) + xyztBin = zero(MVector{4, Int64}) + xyzBin = view(xyztBin,1:3) + + (; binwidth, binstart, covectors, numbins) = params + return_type = typeof(res).parameters[1] + 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 + catch MethodError + (ω, Δω) -> 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 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 - ω) + end + end + end + + 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 + if isnothing(integrated_kernel) # `Delta-function energy' logic + # Figure out which bin this goes in + v[4] = ω + mul!(coords,covectors,v) + coords .= (coords .- binstart) ./ binwidth + coords[.!isfinite.(binwidth)] .= 0 + xyztBin .= 1 .+ floor.(Int64, coords) + + # Check this bin is within the 4D histogram bounds + if all(xyztBin .<= numbins) && all(xyztBin .>= 1) + intensity = static_mode ? res.data[cell_ix] : res.data[iω,cell_ix] + + ci = CartesianIndex(xyztBin.data) + output_intensities[ci] += intensity + output_counts[ci] += 1 + 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 + + # 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) + + # Calculate source scattering vector intensity only once + 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ω] + end + else + error("Energy broadening not yet implemented for histograms with complicated energy axes") + end + end + 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) +end -Also returned is a list of marker indices corresponding to the input points, and -a list of ranges giving the indices of each histogram `x`-axis within a concatenated histogram. -The `density` parameter is given in samples per reciprocal lattice unit (R.L.U.). """ -function q_space_path_bins(ωvals,qs,density,args...;kwargs...) - nPts = length(qs) - params = [] - markers = [] - ranges = [] - total_bins_so_far = 0 - push!(markers, total_bins_so_far+1) - for k = 1:(nPts-1) - startPt = qs[k] - endPt = qs[k+1] - # Density is taken in R.L.U. since that's where the - # scattering vectors are equally spaced! - nBins = round(Int64,density * norm(endPt - startPt)) - # SQTODO: Automatic density that adjusts itself lower - # if there are not enough (e.g. zero) counts in some bins - - param = slice_2D_binning_parameters(ωvals,startPt,endPt,nBins,args...;kwargs...) - push!(params,param) - push!(ranges, total_bins_so_far .+ (1:nBins)) - total_bins_so_far = total_bins_so_far + nBins - push!(markers, total_bins_so_far+1) + integrate_axes!(params::BinningParameters; axes) +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) +""" +function integrate_axes!(params::BinningParameters;axes) + for k in axes + nbins = [params.numbins.data...] + nbins[k] = 1 + params.numbins = SVector{4}(nbins) end - return params, markers, ranges + return params +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 end -q_space_path_bins(sc::SampledCorrelations, qs::Vector, density,args...;kwargs...) = q_space_path_bins(available_energies_including_zero(sc), qs, density,args...;kwargs...) 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..b158d9dcf 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,31 @@ 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 +# Parse "[0.5H,0.3H,0.1H]" into the Julia vector [0.5, 0.3, 0.1]. These strings +# typically arises as the Mantid histogram axis label. 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) + 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..0180c7275 --- /dev/null +++ b/test/test_binning.jl @@ -0,0 +1,69 @@ + +@testitem "Parse nxs" begin + @test Sunny.parse_long_name("[0.5H,0.3K,0.1H]") == [0.5, 0.3, 0.1] + # TODO: Add stripped down .nxs files of various flavors +end + + +@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) + + 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