0.7-style binning

[Lazersmoke]

In Sunny v0.7, there are QGrid{2}s and QPaths. These each contain slightly less information than is needed to fully specify the regions of Q-space which the data comes from. In particular, you can have multiple transverse binning configurations (blue box, orange box) which are centered on exactly the same QPath (red Xs), see image:

This PR adds specify_transverse_binning which "upgrades" a QGrid{2} or QPath into the corresponding BinningParameters, taking the extra information about the transverse binning as input. This function should be used in conjunction with integrate_axes! and energy_resolve! to further add energy-binning information or to integrate out unneeded axes.

It also adds binned_intensities, which does the following:

1. It decides which of the finitely-many available Q scattering vectors in a SampledCorrelations could possibly contribute to a given BinningParameters histogram
2. It constructs a QPoints and uses the existing intensities(...) function to compute the intensities
3. It takes those intensities and sorts them into the bins, returning res with type BinnedIntensities, which is this:

struct BinnedIntensities{T} <: AbstractIntensities
    # 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

Then, the user can call plot_intensities to get a plot with automatically determined axes and axes labels (supports the cases where all but up to two dimensions are integrated out currently). By default, this divides data by counts to best match the idea of a continuous S(q,w) as used by the other intensities types. As an example, here's how you would modify the Example#2 to pick a transverse binning direction and plot:

# q-space orthogonal grid of q vectors
grid = q_space_grid(cryst, [1, 0, 0], range(-10, 10, 200), [0, 1, 0], (-10, 10))

# No binning
res = intensities_static(sc, grid)
plot_intensities(res; saturation=1.0, title="Static Intensities at T = 16 K")

# Yes binning
params = Sunny.specify_transverse_binning(grid,[0,0,1],0.1)
bres = Sunny.binned_intensities(sc,params)
plot_intensities(bres; saturation=1.0, title="Static Intensities at T = 16 K")

The two plots (not shown) are identical up to an overall scale because transverse binning didn't matter in this case.

Note: The momentum axis labels of these plots are different from the non-binned case; they describe the real-space direction of dispersion represented by each axis. This gives a slightly different set of information than the [H,K,L] labels when fewer than three axes are displayed, but this is often the more relevant information in situations where you care about transverse binning effects anyway. This becomes important in the case where things are not orthogonal:

# Non-orthogonal grid of q vectors
grid = q_space_grid(cryst, [1, 0, 0], range(-10, 10, 75), [1, 1, 0], (-10, 10))

# No binning
res = intensities_static(sc, grid)
# Not supported! throws an error
# plot_intensities(res; saturation=1.0, title="Static Intensities at T = 16 K")

# Yes binning
params = Sunny.specify_transverse_binning(grid,[0,0,1],0.1) # Same line of code as before
bres = Sunny.binned_intensities(sc,params)
plot_intensities(bres; saturation=1.0, title="Static Intensities at T = 16 K")

Notice the axes labels are NOT [H,0,0] and [K,K,0], even though that's how the q-grid was specified. If I told you that the transverse binning direction was + 1.0 dz, then you could infer the HKL labels from the given real-space information, but not knowing that and just looking at the plot you couldn't infer that. Instead, you get to know that the x-axis has to do with dispersion in the -b+a real-space direction, which is not something you can know from the HKL labels directly (unless you also knew that [0,0,L] was the third direction).

The PR also un-comments and cleans up a few internal binning methods used by binned_intensities to better support the static case (where energy is integrated out). Static and dynamic correlations should be fully supported!

This PR should be thought of as a reference implementation/problem statement for what's missing in terms of binning support from QGrid, rather than a finished product, so happy for any edits or suggestions :)

[kbarros]

Thanks Sam, this looks like a really nice start. I will need to spend some time digesting it. We should also collect feedback from people who used binning in 0.6.

[Lazersmoke]

I don't want to conflict with the changes in progress, but here are some unit tests! This tests the codepath involving constructing a QGrid followed by specify_transverse_binning, and also using unit_resolution_binning_parameters, but there is no load_nxs test.

It also resulted in a nice graphic :)

Click to show code

# 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,cross(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