correspondence between discrete samples and field

Hey @steven-murray, I'm using the create_discrete_samples function in the main PowerBox class, and noticed that if I naively plot the density of tracers, there is a mismatch with the input field. This is just a tranpose issue related to whether one uses indexing='xy' or indexing='ij' in calls to numpy's meshgrid function. For me the latter is more intuitive and provides a one line fix, though if people generally prefer to stick with numpy defaults then maybe just some clarifying language in the docstring is the better solution. What do you think?

Here's a MWE that shows what I mean:

import numpy as np
import powerbox as pbox
import matplotlib.pyplot as plt

box = pbox.LogNormalPowerBox(N=128, dim=3, pk = lambda k: 0.1*k**-2.,
    boxlength=100)

rho = box.delta_x()

fig, axes = plt.subplots(1, 3, figsize=(12,6))

axes[0].imshow(rho.mean(axis=0), origin='lower')
axes[0].set_title('raw box')

for i, indexing in enumerate(['xy', 'ij']):

    ##
    # Do what "create_discrete_sample" is doing but change 'indexing'
    # kwarg in meshgrid call.

    # Get coordinates of each voxel
    args = [box.x + 0.5 * box.boxlength] * box.dim
    xx, yy, zz = np.meshgrid(*args, indexing=indexing)

    # Make some bin edges to prep for histogram
    xe = ye = ze = np.arange(0, box.boxlength+box.dx, box.dx)

    # Histgram the density "catalog", i.e., a thing that looks like the output
    # of "create_discrete_sample", to recover a field.
    data = np.array([xx.ravel(), yy.ravel(), zz.ravel()]).T
    hist, edges = np.histogramdd(data, bins=[xe, ye, ze],
        weights=rho.ravel(), density=False)

    # Compare to raw density field
    axes[i+1].imshow(hist.mean(axis=0), origin='lower')
    axes[i+1].set_title(f'indexing={indexing}')

steven-murray / powerbox

correspondence between discrete samples and field #15