[yt-dev] RFC: (yt-3.0) Geometry handling redesign

[Matthew Turk]

Hi all,

This email is going to be long, but the changes I am proposing in it
will require feedback.  Next Monday at 1PM Eastern I'll be holding a
hangout on Google+, if you would like to join in and discuss it with
lower latency.  (Of course, if nobody shows up, that's okay too.)  We
don't really have a detailed review process for design decisions like
this, and I hope that this discussion will guide developing that kind
of review process in the future.

Everything described here is in the branch geometry_handling at
https://bitbucket.org/MatthewTurk/yt-refactor

= Background =

Geometry in yt has always been designed around two things:

1) Grid patches
2) Spatial location based on bounding boxes

This has presented a number of excellent opportunities for optimizing
patch-based AMR data.  However, it has posed challenges for supporting
particle-based datasets.  (Particles inside patch-based data are
somewhat easier, although if you want to see one of the worst pieces
of yt, look no further than FakeGridForParticles, which mocks up
particles as grids with varying x,y,z, so that existing
fluid-selection routines can be used.)  Furthermore, the state of
Octree codes is not good.  For relatively small (up to about 64)
domains in RAMSES, we can do analysis with relative ease -- volume
rendering, projections, fluid analysis, etc, all work very nicely.
However, what is required is a costly step of regridding octs at the
outset, to form grid patches with some efficiency.  This regridding is
not stored between runs, and it's not *efficient*.

This has worked so far, but it is showing its strain.  As it stands,
yt cannot directly support SPH codes without somehow making SPH look
like grids.  Octree codes can't be supported without making them look
like grid codes.  Both of these are very slow operations.  And on top
of that, the current state of how data is read in, even for grid
codes, is a bit gross.  There are a number of conditionals, mechanisms
for striding over data (i.e., "_get_grids()" which is called in
parallel, and iterates over grids, and sometimes you have to call
_get_point_indices and sometimes you don't), and probably worst of all
a lot of unnecessary memory duplication because arrays in
regions/spheres/etc are constructed by iterating over grids, filling
in a list of arrays, and then concatenating that list into a second
array.

Furthermore, several times recently the question has come up of
handling and supporting more generic geometries: spherical, polar,
cylindrical.  Unfortunately, because we're so heavily based around the
idea of grids as is, these are difficult to support.  All of our
collision detection has been designed with rectilinear surfaces in
mind.  That being said, I'm really quite happy with how yt *functions*
for patch-based grids.  I'm just not sure that as it stands it is
maintainable, and I think that addressing the geometry issues could
ultimately also help out a great deal with scalability.  (More on that
below.)

= Geometry Overhaul =

For the last couple weeks, I have been working on an overhaul of the
geometry system in yt, attempting to handle parallel striding in a
format-agnostic way, deal with generating fields in as simple a way as
possible, and also attempting to speed things up.  I'm going to
summarize below the changes that I've made, the new design, and the
current status.

== Design: Selectors ==

As it stands, yt data objects (such as Region, Sphere, etc) are mostly
thought of as conduits for data.  They apply a criteria to the domain
and select all of the fluid points or particles within a region or
that fit a specific non-spatial criteria.  I have extended this
concept to construct very lightweight objects called selectors,
similar to how ParticleIO is now conducted for Enzo.  These are
objects in Cython that accept either a set of bounding-box arrays or a
set of spatial points and decide if those points are included in a
given region.  This is essentially a combination of the
get_point_indices and get_cut_mask methods that currently exist.
These have all been designed to be GIL-free, so that in the future
they can be threaded using OpenMP (although it turns out not very much
time is spent on them.)  To date I have implemented spheres, regions,
disks, cutting planes, slices, ortho rays, and grid collections.  They
all give bitwise identical results to the old method and are almost
all substantially faster.  (Cutting planes are particularly better,
regions are slightly slower, and the others are all about 25-40%
faster.)

One primary advantage is that these can (and do) trivially provide
counting methods, masking methods, and point-evaluation methods.  So
we can use the same routines to select particles and individual points
(for collision detection) that we use for selecting grid cells.  This
will allow reuse of data object code with particle outputs, and is a
primary advantage of the refactor.

This also led to separating selection-based data containers from
constructed data containers.

== Design: Ditching the Hierarchy ==

(Hopefully that heading caught some eyes.)  In line with attempting to
reduce Enzo-isms in the code, I have been moving away from
specifically calling something a "hierarchy," and instead moving
toward referring to geometry handlers.  grids, grid_left_edges,
grid_right_edges, and so on will no longer be guaranteed to appear
hanging off these objects.  The new class structure:

GeometryHandler
GridGeometryHandler(GeometryHandler)
OctreeGeometryHandler(GeometryHandler)

and eventually, other methods.  The concept behind this is to isolate
those items that *must* require domain knowledge, mesh knowledge, or
particle/fluid distribution knowledge, and place those in the
GeometryHandler.  All items specific to patch based or block based AMR
will go in GridGeometryHandler and so on.  This will allow for
fundamentally different method of IO for octrees from grids, for
instance.  The geometry handlers are all required to implement the
following methods:

_read_selection(self, fields, data_object, chunk)
_chunk(self, data_object, chunking_style, ngz = 0)

I will describe chunks later.  Note that the first routine,
_read_selection, accepts the names of the fields, and the data object
that describes how to select which portions of a dataset participate
in a given selection.  For the specific case of Enzo (as I am still in
design phase, I have only been working with Enzo data, although I will
move on to FLASH and Boxlib next) read_selection calls a routine that
counts and masks grids, batches operations to files, and then returns
only the necessary data for a selection.  The data object does no
discarding of out-of-bounds or child-masked data.

== Design: Chunking ==

Currently, yt uses de facto chunking of data in the form of grids.
When Derived Quantities want to handle individual grids, one at a
time, they "preload" the data from whatever grids the
ParallelAnalysisInterface thinks they deserve.  These grids are
iterated over, and handled individually, then the result is combined
at the end.  Profiles do something similar.  However, both of these
are de facto, and not really designed.  They rely on calls to
semi-private functions on data objects, manually masking data, on and
on.

To get around this ugly spot, I have started working on redesigning IO
and data handling to function in "chunks."  In this method, data
objects will request that the geometry handler supply a set of chunks
(the function _chunk, described above, conducts this task.)  Chunks
are of the form (IO_unit, size), where IO_unit is only ever managed or
handled by _read_selection, above.

To start the Chunks shuffling over the output, the code calls
data_source.chunks(fields, chunking_style).  Right now only "spatial",
"io" and "all" are supported for chunking styles.  This corresponds to
spatially-oriented division, IO-conserving, and all-at-once (not
usually relevant.)  What happens next is kind of fun.  The chunks
function looks like this:

    def chunks(self, fields, chunking_style, **kwargs):
        for chunk in self.hierarchy._chunk(self, chunking_style, **kwargs):
            with self._chunked_read(chunk):
                self.get_data(fields)
                yield self

Note what it does here -- it actually yields *itself*.  However,
inside the chunked_read function, what happens is that the attributes
corresponding to the size, the current data source, and so on, are set
by the geometry handler (still called a hierarchy here.)  So, for
instance, execution might look like this:

for ds in my_obj.chunks(["Density"], "spatial"):
  print ds is my_obj
  print ds["Density"].size

The first line will actually print True, but the results from the
second one will be the size of (for instance) the grid it's currently
iterating over.  In this way, it becomes much easier to stride over
subsets of data.  Derived quantities now look like this:

        chunks = self._data_source.chunks([], chunking_style="io")
        for ds in parallel_objects(chunks, -1):
            rv = self.func(ds, *args, **kwargs)

It chunks things, evaluates and then stores intermediate results.

This is not meant to replace spatial decomposition in parallel jobs,
but it *is* designed to enable much easier and *mesh-neutral* division
of labor for parallelism and for IO.  If we were to call chunk on an
octree, it no longer has to make things look like grids; it just makes
them look like flattened arrays (unless you chunk over spatial, which
I haven't gotten into yet.)

Essentially, by making the method of subsetting and striding over
subsetted data more compartmentalized, the code becomes more clear and
more maintainable.

= Status =

As it stands, I am able to select data very nicely -- this means
slices, profiles, derived quantities all work.  However (see below)
particles still need some work and handling, and constructed data
objects (smoothed covering grids, projections) have yet to be ported.

= Open Questions =

However, the main reason I am writing -- when the code is not yet
fully functional! -- is to ask for advice and consideration on a few
points.

== Particle Types ==

In the past, we've talked about handling particle types differently
depending on the particle.  Because this is destined for yt-3.0, I
would like to propose that we break compatibility with old scripts in
favor of clarity in our particle types.  Specifically, rather than
doing:

obj["particle_position_x"][obj["particle_type"] == 1]

We do,

obj[("particle_positions_x", 1)]

(or obj[("particle_positions_x", "PopIICluster")] if you are using a
code that supports such expression.)

And if a field is found to be a particle type field, but without a
tuple, we simply raise a KeyError.  We can emulate old behavior with:

obj[("particle_position_x",None)]

But I am inclined to say that I would rather *mandate* particle type
selection, rather than merely allow it.  This will lead to much
greater clarity, I think.

[+-][01] on

== Particle Access ==

Currently, we access particles off a data_object:

sphere["particle_position_x"]

Do we want to keep this, or move to being much more explicit, that we
require sphere.particles["particle_position_x"] ?  I can see both
sides to this.  Explicit is good, and would allow us to be both more
clear with users and about fields (and with our IO methods) but
allowing old behavior encourages not having to necessarily worry about
if something is a fluid or a particle, which is a good thing I think.

== User Survey ==

Should we conduct a brief (only a handful of questions) polling users
for their feelings on the changes or potential changes to
nomenclature, particles, and so on?  We could use it also as a time to
identify other outstanding issues with yt.

= Conclusion =

Does this all make a sort of sense?

[+-][01] on the design presented so far

If you're interested in helping out, please get in touch.  I'm still
laying groundwork in some areas, but in other areas I would appreciate
help and advice moving forward.

Thanks to everyone who made it this far,

Matt