On Monday I rewrote the (adaptive) projection backend to use Quad
Trees for point placement instead of the current system. As it
stands, the current system for projecting works pretty well, except
for a couple caveats:
- It has high peak memory usage
- It can only be decomposed in the image plane, so inline projections
can't be made, and load balancing AMR was poor
- It's complicated, but that complication comes from wanting things to be fast
Essentially, it's fast by wall-clock standards, but you need to run in
parallel for big datasets. Plus, if you want to run with big
datasets, you probably need to fully occupy nodes, because of the peak
Anyway, I was getting very frustrated with points #1 and #3 (but not
#2, oddly) and so I wrote a new backend using Quad Trees. We did
everything using int math anyway, but because joining grids was
O(NxM), things got a bit slow in some domains. The conversion wasn't
I ran some simple benchmarks. These were all conducted on the Triton
resource, in serial, on a single node. They do not reflect time for
disk IO, but they should roughly reflect real memory load.
* For my 31 level PopIII binary star run (875 million cells, 8000
grids), it took 9 seconds with peak memory usage of 176 megs.
* For a large, 768^3 with (peak) 12 levels of refinement (874 million
cells in 125,000 grids), it takes 81 seconds with peak memory load of
* For the 512^3 L7 lightcone run, at a time when the data had 300,000
grids, it took 170 seconds with peak memory usage of 2.5 gigs.
In all cases, the results are correct according to the standard check
of projecting "Ones". These are crazy improvements, especially
considering that these are times for serial projections.
And important note here is that this *will* parallelize independent of
the image plane, although I have yet to see the need to parallelize
just yet. The load balancing may be better optimized based on image
plane decomposition, but the algorithm should work with independent
joins. Additionally in parallel the final join at the end of the
calculation may have higher peak memory usage if the decomposition is
not via image plane decomposition, but it should still work.
Unfortunately, right now my time is very constrained; so I am going to
put out there that if someone would be willing to help me break apart
the existing code and insert this code, then we can swap out the
engines. But in the absence of that, I'm not sure that it'll get
inserted into the main line before the Summer, or at least the late
Spring. As is, I can definitely provide a very simple and manual
interface to use it (and I intend to do so) but integrating it into
the mainline is not going to be my priority right now. But, again, I
will provide *an* interface, even if it's not a *final* interface and
likely not a *nice* interface.
Anyway, I'm pretty excited about this -- definitely a huge improvement
over the previous algorithm, and we still retain the *adaptive* nature
of the projections: project once, plot many.
I wrote a remote image viewer, so that you can interactively explore
things like projections without having to copy *any* of the data
locally. You can use it inside the Chaco-based explorer; all
pixelization is done remotely and only pixelized buffers are sent back
over the wire. I've stuck an example here:
but it'll be documented (along with the rest of the image panner
stuff) in the 1.7 docs, specifically how to get the FURL info
(~/.ipython/security/ipcontroller-mec.furl). Anyway, I'm finding it
incredibly useful for exploring a huge amount of remote-living data.
Hopefully somebody else will find it useful, too!
I pulled the new hg version of yt today, when installing (on both
ranger and my local mac) I get an error in the compile for _MPL.c in
error: Command "gcc -arch ppc -arch i386 -isysroot /Developer/SDKs/
MacOSX10.4u.sdk -fno-strict-aliasing -Wno-long-double -no-cpp-precomp -
mno-fused-madd -fno-common -dynamic -DNDEBUG -g -O3 -I/Library/
include/python2.5 -c yt/raven/_MPL.c -o build/temp.macosx-10.3-
i386-2.5/yt/raven/_MPL.o" failed with exit status 1
Dr. Eric J. Hallman
NSF Astronomy and Astrophysics Postdoctoral Fellow
Center for Astrophysics and Space Astronomy
University of Colorado at Boulder
hallman (at) casa.colorado.edu
Phone: (312) 725-4626
Just a note: with mercurial 1.5, the meaning of "hg heads" has changed
slightly -- it is no longer neutral to branch mechanics. This means
that if a named branch has been merged, but not "closed", it will
still show up in the default output of hg heads. For the most part I
don't think this is really what we want, because nobody has been
closing branches. So, if you use hg 1.5 right now on any repo with
multiple named branches, you're going to see way more than you
expected. To get around this,
hg heads -t
will only show revisions with no children *at all*, rather than the
new default behavior, which was to show heads with no children *in
that named branch*. You can alias this to "hg heads" by putting this
in the [alias] section of your ~/.hgrc :
heads = heads -t
(Another fun alias is "nudge = push -r ." which will only push the
necessary changesets to fully describe the current revision -- any
other trees of development will not be included in the push. Useful
for when doing local development with multiple heads!!)
The vr-multivariate branch has been merged into the main yt-hg
repository. I worked a little bit on speed today and for
ColorTransferFunction it's now competitive. Additionally, I've added
a Projection mechanism that should enable off-axis projections (this
will be documented very soon.) The off-axis projections are pretty
nice, I think -- the only hold up on that front is that I want to
ensure that I understand the units. I've attached an example image,
which was created by the volume renderer using this script:
note that we output directly to PNG here -- no matplotlib intervention at all!
This is from data supplied by Cameron Hummels, which he was previously
using a stacked cutting plane method to look at off-axis. I'm excited
about the off-axis projections, because my protostars are rarely
aligned with the grids -- so now off-axis projections of baryon
quantities are going to become a part of my standard analysis. Again,
the units need to be examined before this is finished, and I'll handle
Let me know if there are any issues. I'll probably port this back to
trunk this weekend. By merging this branch we now have in the main hg
* Off-axis projections
* Multivariate volume rendering
* Planck transfer functions
* Approximate Rayleigh scattering
Over the last little while, the volume renderer has really taken off
in popularity! I'm a bit conflicted about this, since I'd rather
people were super duper into things like phase plots, radial profiles
and clump finding, but it's also good to know people are using yt for
volume visualization as well as other types of analysis. I used it
myself to what I think was great effect at a conference last week, and
I hope that other people are getting similar use out of it.
Anyway, a couple months ago John wrote to yt-users and suggested that
we try doing multivariate volume rendering, similar to what's in
Kaehler, Wise, & Abel. For a while mostly nothing happened, but over
the Friday->Sunday, Jeff and I hammered it out and it's now working.
Jeff took care of the emission spectra (with some assistance from
John's RGB code), I took a whack at adding multiple variables to the
partitioned grids and the volume rendering itself, and we can now do
multivariate rendering -- although absorption still needs some work,
which I'll discuss below.
The new mechanism for grid partitioning only returns the indices into
the vertex-centered data. The actual processing of the
vertex-centered data. This means the actual locating of sweeps is
much faster -- right now we don't take advantage of that, but I can
see it being useful in the future. A few weeks ago I wrote a
distributed object collection so that we could do 3D domain decomp to
get the bricks, then pass them around between processes for the actual
volume rendering, which was parallelized via 2D image plane decomp
into decomposed inclined boxes. Now the 3D decomp will only pass
around single arrays that describe the indices into the bricks, and
then each image-decomp processor will generate and slice up the data
as necessary -- this isn't yet working, but it's a clear process and
I'll try to take care of it soon.
Transfer functions now have several more parameters, and I'm working
to convert the old-style transfer functions to this. The idea is that
now you add field interpolation tables which correspond to a specific
field, link them to channels (R, Ralpha, G, Galpha, B, Balpha), and
optionally have their values weighted by other fields. So for an
example, here's how the Planck function sets up its various channels:
(there's a bit of cruft where I tested adding lines for density absorption.)
It creates a TransferFunction, links it to a field (0) and a weighting
field (2), then links it to channels. There's a lot of complexity in
this, because it's a complex thing and we want to enable very complex
behavior, but for the most part it will be hidden in the
ColorTransferFunction, which will abstract this all away.
Absorption works slightly differently now. Because we do
back-to-front integration, I believe we do not have to store the
accumulated opacity -- so while right now the volume renderer has
6xNxM arrays, I think it should just be 3. Additionally, we're
presented with the radiative transfer equation:
dI/ds = -alpha * I_0 + j
Here j would be our emission at a given sample, I_0 is the current
value of the image plane, and ds is the distance between samples.
Since we're doing back-to-front, I believe we are directly solving
this, and so our alpha fed in directly corresponds to the absorption.
However, because ds is so low, careful calibration between alpha and
emissivity is necessary. Alternately, converting everything back to
cgs is probably not difficult... (approximate) Scattering should be
easy to put in, it just hasn't been yet.
However, that's neither here nor there. I believe that some exciting
things can be done with this (Britton and Jeff and I have chatted
about line transfer as well as LyA maps) and I think it's rapidly
converging on a workable implementation.
I made up some images with it. The first is a sample dataset of
John's with an HII region:
The second is one of my protostars from my thesis (I had to
aggressively clip the image here) just at around 1.5e4K. This image
is 3AU on a side:
And, just for fun, I got the transfer function widget to steal the
camera position from the VTK interface, so I made up an image of the
two playing nicely together:
Anyway, that's my status update. All of the MultiVariate stuff went
in the "vr-multivariate" branch, so feel free to take a look. (You'll
probably have to re-cython.) An example script is here:
I'll shoot another, much shorter, email to the list when it's working
and merged back into the main hg branch. We'll then have a parallel
software volume renderer that is fully multi-channel, multi-variate,
and comes with black body transfer functions. All of this is only
achievable thanks to the development community we have, and our team
efforts. So, thank you all for your hard work and your enthusiasm.
I wrote a new module (in branch yt of hg) that will do simple image
writing; it only depends on libpng (no matplotlib dependencies) and
does not do anything fancy: no text, no callbacks, no colormaps. But,
it's really quite fast and it comes with a couple colormaps -- algae,
jet, gist_stern and hot. The API is pretty easy. You can either
import yt.extensions.image_writer as iw
iw.write_image(image, filename, color_bounds = None, cmap = "algae")
where image is a numpy array of values to be colormapped, or you can
directly call the image writer:
where some_image_data is uint8 and is (N,M,4).
Unfortunately, generating FixedResolutionBuffers still requires that
matplotlib be imported; I'll work on getting around this sometime,
because if this ends up working we should be able to run on Kraken in
the compute nodes without importing matplotlib. Although, Cray has
now announced they're going to support dynamically linked libraries on
their compute nodes, so maybe we won't need this for long.
That being said, it should also be pretty handy for things like
writing out volume renderings and anywhere else we previously used
pylab to save images rather than plots. Let me know if you have any
ideas or problems!
I know several of you are or will be involved in test problem coverage
for the new release of Enzo. I would like to underscore that
problems, cosmetic or real, that show up in yt as a result of
increasing test coverage of Enzo are *critically important* to fix.
This goes even to things like the display of field names and so on.
Please feel free to fix any problems that you find along the way.