yt 4.0 adding arbitrary_grid to interpolate SPH data #1814

AshKelly · 2018-06-04T11:02:29Z

PR Summary

PR to add a new feature to interpolate SPH fields onto an arbitrary grid.

Add the cython pixelize function which takes data in 3D and a returns a 3D buffer
wire into the arbitrary grid container
Add test to compare the arbitrary_grid with the slice plot

An example use case of this new wiring is as follows,

import yt
from yt.testing import \
    assert_equal, \
    fake_sph_orientation_ds, fake_random_ds

ds = fake_sph_orientation_ds()

agrid = ds.arbitrary_grid([0.0, 0.0, 0.0], [2.99, 2.99, 2.99],dims=[5, 2, 2])
dens = agrid['gas', 'density']

print(dens)

PR Checklist

Code passes flake8 checker
New features are documented, with docstrings and narrative docs
Adds a test for any bugs fixed. Adds tests for new features.

…identifying any neighbouring pixels? will attempt to debug later

…se) always debug cython code with bounds checking on

… function

…kernel_arbitrary_grid

…ph arbitrary_grid

…d sph field

…ata container

…rbitary grid class. Appears to be almost working -- other than returning an array of zeros. will debug

…pears to be working

…-arbitrary-grid pulling remote changes

ngoldbaum · 2018-06-06T21:04:28Z

yt/utilities/lib/pixelization_routines.pyx

+                            buff[xi, yi, zi] += coeff * kernel_func(qxy)
+
+    if use_norm:
+        # now we can calculate the normalized image buffer we want to


i'd just say buffer since this is a 3D grid, not an image.

ngoldbaum · 2018-06-06T21:05:08Z

yt/data_objects/construction_data_containers.py

+            if(is_sph_field):
+                self._fill_sph_particles(fields)
+            else:
+                self._fill_particles(part)


Let's work through this together during our call tomorrow, I'm a little confused about your concern in the PR description and I think it might be easier to chat about it over a telecon.

ngoldbaum · 2018-06-06T21:06:24Z

yt/data_objects/construction_data_containers.py

+                pz = chunk[(ptype,'particle_position_z')].in_units('cm')
+                hsml = chunk[(ptype,'smoothing_length')].in_units('cm')*10.0
+                mass = chunk[(ptype,'particle_mass')].in_units('g')
+                dens = chunk[(ptype,'density')].in_units('g/cm**3')


we could probably be doing this in code units. both here and in the image pixelizers. it doesn't really matter which unit system we choose as long as we stick to a single system, but code units is probably the most natural choice.

ngoldbaum · 2018-06-06T21:06:51Z

yt/data_objects/construction_data_containers.py

+                bounds[4] = self.left_edge[2]
+                bounds[1] = self.right_edge[0]
+                bounds[3] = self.right_edge[1]
+                bounds[5] = self.right_edge[2]


Are these always in 'cm'?

Does it even matter? it might not if the cython function ignores units.

ngoldbaum · 2018-06-06T21:07:39Z

yt/data_objects/construction_data_containers.py

+        for field in fields:
+            dest = np.zeros(self.ActiveDimensions, dtype="float64")
+
+            for chunk in self._data_source.chunks(fields, "io"):


i think you should pass in [field] here?

ngoldbaum · 2018-06-06T21:08:39Z

yt/data_objects/construction_data_containers.py

+                bounds[5] = self.right_edge[2]
+
+                pixelize_sph_kernel_arbitrary_grid(dest,px,py,pz,hsml,mass,
+                                                   dens,mass,bounds)


hmm, did you forget to pass in the field values for the field we're smoothing onto the grid?

ngoldbaum · 2018-06-06T21:09:13Z

yt/utilities/lib/pixelization_routines.pyx

@@ -1154,6 +1154,130 @@ def pixelize_sph_kernel_slice(
                    continue
                buff[xi, yi] += buff_num[xi, yi] / buff_denom[xi, yi]

+@cython.initializedcheck(True)
+@cython.boundscheck(True)


these should be turned off when we're sure everything is correct

ngoldbaum · 2018-06-06T21:09:35Z

yt/utilities/lib/pixelization_routines.pyx

+        for j in prange(0, posx.shape[0]):
+            if j % 1000 == 0:
+                with gil:
+                    PyErr_CheckSignals()


i'd like to see a progress bar here as well

ngoldbaum · 2018-06-06T21:09:50Z

yt/utilities/lib/pixelization_routines.pyx

+            w_j = pmass[j] / pdens[j] / hsml[j]**3
+            coeff = w_j * quantity_to_smooth[j]
+
+            # Now we know which pixels to deposit onto for this particle,


voxels, not pixels

ngoldbaum · 2018-06-12T22:29:29Z

Not at all :)

In general you don’t need to ask permission to use code I share publicly, just assume the answer is yes.

… work with derived fields as advised by Nathan

…utput for the field being interpolated

qobilidop · 2018-06-13T00:43:29Z

@ngoldbaum Just tried your script on FIRE_M12i_ref11/snapshot_600.hdf5, and I see the VR of this snapshot for the first time! I can't wait to do some VR on the FIRE data I'm currently analyzing!

@AshKelly Thanks for your fantastic work! This is super useful.

qobilidop · 2018-06-13T01:21:38Z

I'm trying to compare a slice buffer from the arbitrary grid and a slice buffer made from the slice plot directly using the following script:

import yt
import numpy as np
from yt.units import kpc

ds = yt.load('GadgetDiskGalaxy/snapshot_200.hdf5')
_, c = ds.find_max(('gas', 'density'))
width = 50 * kpc
buff_size = 256
field = ('gas', 'density')

# buffer from arbitrary grid
ag = ds.arbitrary_grid(c - width / 2,
                       c + width / 2,
                       [buff_size]*2 + [1])
buff_ag = ag[field][:, :, 0].to('g*cm**-3').d

# buffer from slice
p = yt.SlicePlot(ds, 'z', field, center=c, width=width)
p.set_buff_size(buff_size)
buff_slc = p.frb.data[field].to('g*cm**-3').d

# comparison plot
from matplotlib.colors import LogNorm
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.imshow(buff_ag, norm=LogNorm())
plt.colorbar()
plt.title('arbitrary grid')

plt.subplot(1, 2, 2)
plt.imshow(buff_slc, norm=LogNorm())
plt.colorbar()
plt.title('slice')

plt.tight_layout()
plt.savefig('ag_slc_cmp.png')

And here's the result:

It looks like a lot of particles are missed in the arbitrary grid.

chummels · 2018-06-13T02:22:00Z

Holy cow--that VR on Gadget data is awesome!

AshKelly · 2018-06-13T09:06:57Z

I disagree @qobilidop. I have a suspicion the difference is due to normalization. In the slice plot we don't normalize because the results look better without it, but since the arbitrary grid may be used for analysis we went with normalization by default. I can just check this on my machine though.

update:
This doesn't seem to be normalization causing this.

I suspect the confusion is coming from how the arbitrary_grid works. You can select a huge region, but, just because you have selected every particle in that region does not mean their smoothing lengths will overlap with the center of the voxels, remember the smoothing length uses separation in x,y and z for the arbitrary_grid. In the slice plot, the smoothing length does not use contributions from the z-direction and so particles far away in the z direction will still be included, I believe. (Please correct if I'm wrong here @ngoldbaum)

There was certainly a bug, I'd missed a sqrt which could have been missing particles. I have now fixed this. There also seems to be a difference in the orientation of the return arrays from slice and arbitrary grid. Maybe slice has a transpose I missed or something. If I modify your script to this

import yt
import numpy as np
from yt.units import kpc

ds = yt.load('GadgetDiskGalaxy/snapshot_200.hdf5')
_, c = ds.find_max(('gas', 'density'))
width = 50 * kpc
buff_size = 256
field = ('gas', 'density')

# buffer from arbitrary grid
ag = ds.arbitrary_grid(c - width / 2,
                       c + width / 2,
                       [buff_size]*2 + [1])
buff_ag = ag[field][:, :, 0].to('g*cm**-3').d
buff_ag = buff_ag.T

# buffer from slice
p = yt.SlicePlot(ds, 'z', field, center=c, width=width)
p.set_buff_size(buff_size)
buff_slc = p.frb.data[field].to('g*cm**-3').d

# comparison plot
from matplotlib.colors import LogNorm
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.imshow(buff_ag, norm=LogNorm())
plt.colorbar()
plt.title('arbitrary grid')

plt.subplot(1, 2, 2)
plt.imshow(buff_slc, norm=LogNorm())
plt.colorbar()
plt.title('slice')

plt.tight_layout()
plt.savefig('ag_slc_cmp.png')

then I get the following:

but I don't know why values are off.

I think the following line is causing the strange results in the grid,
h_j3 = fmax(hsml[j]*hsml[j]*hsml[j], dx*dy*dz)
as dz is huge, this will make the smoothing length huge so the contributions from most particles will be very small! Thoughts on the best way to deal with this? If we change this then particles far away won't be included in the contribution even if we are in the voxel - but in its current state the slice images look strange.

I also think we can avoid the sqrt I had to add by saying
2 * x < 1
is equivalent to
4 * x**2 < 1
just wanted to check that this logic is sound?

dnarayanan · 2018-06-13T12:03:51Z

this is amazing stuff! well done!

ngoldbaum · 2018-06-13T14:46:35Z

@qobilidop I suspect if you plot that on the same colorbar scale the differences will appear smaller.

AshKelly · 2018-06-13T15:09:39Z

They definitely seem to be getting the same results, but the grid looks massively over smoothed due to the
h_j3 = fmax(hsml[j]*hsml[j]*hsml[j], dx*dy*dz)

ngoldbaum · 2018-06-13T15:39:23Z

Do you understand why there are zones with no data in the arbitrary grid but not in the slice?

AshKelly · 2018-06-13T17:39:29Z

no I don't, I'm going to look into that now

AshKelly · 2018-06-13T17:58:04Z

It looks like we aren't passing in as many particles in the arbitrary_grid call as we are in the slice call. Any reason why this might be?

Every chunk passes more particles to pixelize_slice and than pixelize_arbitrary ...

ngoldbaum · 2018-06-13T18:10:11Z

Hmm can you share the test script you used to check? Just the one Bili shared upthread?

AshKelly · 2018-06-13T18:32:05Z

I'm using this script (I think its the same as upthread)

import yt
import numpy as np
from yt.units import kpc

ds = yt.load('GadgetDiskGalaxy/snapshot_200.hdf5')
_, c = ds.find_max(('gas', 'density'))
width = 50 * kpc
buff_size = 256
field = ('gas', 'density')

# buffer from arbitrary grid
ag = ds.arbitrary_grid(c - width / 2,
                       c + width / 2,
                       [buff_size]*2 + [1])
buff_ag = ag[field][:, :, 0].to('g*cm**-3').d
buff_ag = buff_ag.T

# buffer from slice
p = yt.SlicePlot(ds, 'z', field, center=c, width=width)
p.set_buff_size(buff_size)
buff_slc = p.frb.data[field].to('g*cm**-3').d

# comparison plot
from matplotlib.colors import LogNorm
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
plt.imshow(buff_ag, norm=LogNorm(), vmin=10**(-29), vmax=10**(-22))
plt.colorbar()
plt.title('arbitrary grid')

plt.subplot(1, 2, 2)
plt.imshow(buff_slc, norm=LogNorm(), vmin=10**(-29), vmax=10**(-22))
plt.colorbar()
plt.title('slice')

plt.tight_layout()
plt.savefig('ag_slc_cmp.png')

with the extra prints I've just committed. You can see that a differing number of particles is passed. I presume something is going on with the region/slice code?

ngoldbaum · 2018-06-13T20:54:51Z

yt/utilities/lib/pixelization_routines.pyx

+                        # yt's kernel functions use a different convention than
+                        # the SPLASH paper (following Gadget-2), and qxy is
+                        # twice the value of q used in the SPLASH paper
+                        qxy = 2.0 * math.sqrt(posx_diff + posy_diff + posz_diff)


I think we should just call this q here and in the other SPH pixelizers. It's called q in the paper, so it makes it a little easier to compare with the paper. In addition, here z is included as well, so qxy is a bit of a misnomer.

I'd also add in the comment above to see equation 5 in the SPLASH paper.

It also should be one half the value in the SPLASH paper, and the factor of 2.0 should be deleted. See the patch I sent you over slack.

…tor of 2 from two of the pixelizers

AshKelly · 2018-06-13T21:31:31Z

Just need to check this still passes my test... will wait for the results

…id results with those of the slice

ngoldbaum · 2018-06-18T14:59:55Z

Great job on seeing this one through Ash! 🎆 💯 🥇 🎉

AshKelly added 21 commits June 1, 2018 17:29

beginning to add a cython function to pixelize sph data onto a 3D grid

998b623

Identified bug - do not seem to be traversing the for loops. i.e not …

3f66202

…identifying any neighbouring pixels? will attempt to debug later

Fixed the bug in the sph arbritrary grid! note to self (and anyone el…

48c4f5c

…se) always debug cython code with bounds checking on

updated the comments on the pixelize_sph_kernel_arbitrary_grid cython…

7677613

… function

re structured to reduce the number of sqrts required in pixelize_sph_…

5ad475f

…kernel_arbitrary_grid

removed an unused variable, ran flake8 and tightened up comments in s…

ab56000

…ph arbitrary_grid

attempting to add an sph interface to the arbitrary grid

fcade25

make progress with arbitrary grid, but the interface is very hacky

bab24d4

_generate_sph_field is returning a non-zero array for the interpolate…

c17277e

…d sph field

fixed an array initialiser in _generate_sph_field

cd1c946

neatened up the arbitrary grid api for sph. removed the hardcode field

634ad5f

moved the generate_sph_field method to the arbitrary grid class.

df9961e

modified to only call the method for arbitrary grids and not any yt d…

5821d69

…ata container

refactored the front end such that the sph particle field is in the a…

961cbb1

…rbitary grid class. Appears to be almost working -- other than returning an array of zeros. will debug

finally wired up the _fill_sph_particles for YTArbitraryGrid. this ap…

55dad77

…pears to be working

cleaning up some accidental changes

8d869ce

frontend for 3d pixelisation of sph data working

2798633

Merge branch 'yt-4.0' of https://github.com/yt-project/yt into yt-4.0…

a02426f

…-arbitrary-grid pulling remote changes

pulled remote changes and merged

3d92947

hopefully fixed the global check on sph attribute

eefb104

ensuring is sph init, however, believe I could move the sph check

df7fa2f

ngoldbaum reviewed Jun 6, 2018

View reviewed changes

AshKelly added 2 commits June 12, 2018 23:39

added sph field to tqdm progress bar and updated some flow control to…

48c5b5a

… work with derived fields as advised by Nathan

changed variable to improve readability and used a more descriptive o…

cee185a

…utput for the field being interpolated

fixed sqrt bug in arb

5cf82c6

re-enabled normalizatio

746c80c

added prints to show the particles being passed to the pixelizers

54448f9

ngoldbaum reviewed Jun 13, 2018

View reviewed changes

AshKelly added 3 commits June 13, 2018 22:03

replaced qxy with the more familiar q. also removed the incorrect fac…

67948d0

…tor of 2 from two of the pixelizers

removed the kwarg over-riding the default normalization option

10c38a5

Merge branch 'yt-4.0-arbitrary-grid' into HEAD

6a2993c

AshKelly added 3 commits June 18, 2018 10:47

removed prints and added normalization

e59cea6

set normalization by default. also added a test to compare the arb gr…

508d494

…id results with those of the slice

neatened up some formatting and fixed an import in the sph tests

06e1d31

ngoldbaum merged commit 53cac3f into yt-project:yt-4.0 Jun 18, 2018

AshKelly deleted the yt-4.0-arbitrary-grid branch June 19, 2018 12:48

AshKelly mentioned this pull request Aug 13, 2018

Add the ability to smooth SPH fields onto a covering grid or arbitrary grid #1436

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

yt 4.0 adding arbitrary_grid to interpolate SPH data #1814

yt 4.0 adding arbitrary_grid to interpolate SPH data #1814

AshKelly commented Jun 4, 2018 •

edited

Loading

ngoldbaum Jun 6, 2018

ngoldbaum Jun 6, 2018

ngoldbaum Jun 6, 2018

ngoldbaum Jun 6, 2018

ngoldbaum Jun 6, 2018

ngoldbaum Jun 6, 2018 •

edited

Loading

ngoldbaum Jun 6, 2018

ngoldbaum Jun 6, 2018

ngoldbaum Jun 6, 2018 •

edited

Loading

ngoldbaum Jun 6, 2018

ngoldbaum commented Jun 12, 2018

qobilidop commented Jun 13, 2018

qobilidop commented Jun 13, 2018 •

edited

Loading

chummels commented Jun 13, 2018

AshKelly commented Jun 13, 2018 •

edited

Loading

dnarayanan commented Jun 13, 2018

ngoldbaum commented Jun 13, 2018

AshKelly commented Jun 13, 2018

ngoldbaum commented Jun 13, 2018

AshKelly commented Jun 13, 2018

AshKelly commented Jun 13, 2018

ngoldbaum commented Jun 13, 2018

AshKelly commented Jun 13, 2018 •

edited

Loading

ngoldbaum Jun 13, 2018 •

edited

Loading

ngoldbaum Jun 13, 2018

AshKelly commented Jun 13, 2018

ngoldbaum commented Jun 18, 2018

yt 4.0 adding arbitrary_grid to interpolate SPH data #1814

yt 4.0 adding arbitrary_grid to interpolate SPH data #1814

Conversation

AshKelly commented Jun 4, 2018 • edited Loading

PR Summary

PR Checklist

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ngoldbaum Jun 6, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ngoldbaum Jun 6, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ngoldbaum commented Jun 12, 2018

qobilidop commented Jun 13, 2018

qobilidop commented Jun 13, 2018 • edited Loading

chummels commented Jun 13, 2018

AshKelly commented Jun 13, 2018 • edited Loading

dnarayanan commented Jun 13, 2018

ngoldbaum commented Jun 13, 2018

AshKelly commented Jun 13, 2018

ngoldbaum commented Jun 13, 2018

AshKelly commented Jun 13, 2018

AshKelly commented Jun 13, 2018

ngoldbaum commented Jun 13, 2018

AshKelly commented Jun 13, 2018 • edited Loading

ngoldbaum Jun 13, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

AshKelly commented Jun 13, 2018

ngoldbaum commented Jun 18, 2018

AshKelly commented Jun 4, 2018 •

edited

Loading

ngoldbaum Jun 6, 2018 •

edited

Loading

ngoldbaum Jun 6, 2018 •

edited

Loading

qobilidop commented Jun 13, 2018 •

edited

Loading

AshKelly commented Jun 13, 2018 •

edited

Loading

AshKelly commented Jun 13, 2018 •

edited

Loading

ngoldbaum Jun 13, 2018 •

edited

Loading