[Yt-svn] yt-commit r1696 - in trunk/yt: . _amr_utils extensions/volume_rendering lagos raven
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Sun Apr 25 10:00:34 PDT 2010
Author: mturk
Date: Sun Apr 25 10:00:23 2010
New Revision: 1696
URL: http://yt.enzotools.org/changeset/1696
Log:
Backport from hg:
* New multivariate volume rendering. This allows for having multiple
different ways of controlling (independently) the emission, absorption, and
colors of fields in a volume rendering -- for instance, Planck emission,
density weighted, with approximate rayleigh scattering.
* New off-axis projections in the volume rendering module
* New PNG writer to directly output PNG files, rather than using matplotlib
* Support for the Tiger code
* Colored logging output (optional)
* Some fixes for the fixed res projections
* John's fix for setting specific marker types in the particle plots
The off-axis projections are currently not "fixed" with units, but that should
be coming soon. The direct PNG imaging does not (and probably will not, at
least for a while) support colorbars or any other compositing of information.
Added:
trunk/yt/_amr_utils/fortran_reader.pyx
trunk/yt/_amr_utils/png_writer.pyx
Modified:
trunk/yt/_amr_utils/VolumeIntegrator.pyx
trunk/yt/amr_utils.pyx
trunk/yt/config.py
trunk/yt/extensions/volume_rendering/TransferFunction.py
trunk/yt/extensions/volume_rendering/__init__.py
trunk/yt/extensions/volume_rendering/grid_partitioner.py
trunk/yt/extensions/volume_rendering/image_handling.py
trunk/yt/extensions/volume_rendering/software_sampler.py
trunk/yt/funcs.py
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/BaseGridType.py
trunk/yt/lagos/DataReadingFuncs.py
trunk/yt/lagos/FieldInfoContainer.py
trunk/yt/lagos/HierarchyType.py
trunk/yt/lagos/OutputTypes.py
trunk/yt/lagos/ParallelTools.py
trunk/yt/logger.py
trunk/yt/mods.py
trunk/yt/physical_constants.py
trunk/yt/raven/Callbacks.py
trunk/yt/raven/ColorMaps.py
trunk/yt/setup.py
Modified: trunk/yt/_amr_utils/VolumeIntegrator.pyx
==============================================================================
--- trunk/yt/_amr_utils/VolumeIntegrator.pyx (original)
+++ trunk/yt/_amr_utils/VolumeIntegrator.pyx Sun Apr 25 10:00:23 2010
@@ -73,87 +73,160 @@
cdef class VectorPlane
+cdef struct FieldInterpolationTable:
+ # Note that we make an assumption about retaining a reference to values
+ # externally.
+ np.float64_t *values
+ np.float64_t bounds[2]
+ np.float64_t dbin
+ np.float64_t idbin
+ int field_id
+ int weight_field_id
+ int weight_table_id
+ int nbins
+ int pass_through
+
+cdef void FIT_initialize_table(FieldInterpolationTable *fit, int nbins,
+ np.float64_t *values, np.float64_t bounds1, np.float64_t bounds2,
+ int field_id, int weight_field_id = -1, int weight_table_id = -1,
+ int pass_through = 0):
+ fit.bounds[0] = bounds1; fit.bounds[1] = bounds2
+ fit.nbins = nbins
+ fit.dbin = (fit.bounds[1] - fit.bounds[0])/fit.nbins
+ fit.idbin = 1.0/fit.dbin
+ # Better not pull this out from under us, yo
+ fit.values = values
+ fit.field_id = field_id
+ fit.weight_field_id = weight_field_id
+ fit.weight_table_id = weight_table_id
+ fit.pass_through = pass_through
+
+cdef np.float64_t FIT_get_value(FieldInterpolationTable *fit,
+ np.float64_t *dvs):
+ cdef np.float64_t bv, dy, dd, tf
+ cdef int bin_id
+ if fit.pass_through == 1: return dvs[fit.field_id]
+ bin_id = <int> ((dvs[fit.field_id] - fit.bounds[0]) * fit.idbin)
+ dd = dvs[fit.field_id] - (fit.bounds[0] + bin_id * fit.dbin) # x - x0
+ if bin_id > fit.nbins - 2 or bin_id < 0: return 0.0
+ bv = fit.values[bin_id]
+ dy = fit.values[bin_id + 1] - bv
+ if fit.weight_field_id != -1:
+ return dvs[fit.weight_field_id] * (bv + dd*dy*fit.idbin)
+ return (bv + dd*dy*fit.idbin)
+
cdef class TransferFunctionProxy:
- cdef np.float64_t x_bounds[2]
- cdef np.float64_t *vs[4]
- cdef int nbins
+ cdef int n_fields
+ cdef int n_field_tables
cdef public int ns
- cdef np.float64_t dbin, idbin
- cdef np.float64_t light_color[3]
- cdef np.float64_t light_dir[3]
- cdef int use_light
+
+ # These are the field tables and their affiliated storage.
+ # We have one field_id for every table. Note that a single field can
+ # correspond to multiple tables, and each field table will only have
+ # interpolate called once.
+ cdef FieldInterpolationTable field_tables[6]
+ cdef np.float64_t istorage[6]
+
+ # Here are the field tables that correspond to each of the six channels.
+ # We have three emission channels, three absorption channels.
+ cdef int field_table_ids[6]
+
+ # We store a reference to the transfer function object and to the field
+ # interpolation tables
cdef public object tf_obj
+ cdef public object my_field_tables
+
def __cinit__(self, tf_obj):
- self.tf_obj = tf_obj
+ # We have N fields. We have 6 channels. We have M field tables.
+ # The idea is that we can have multiple channels corresponding to the
+ # same field table. So, we create storage for the outputs from all the
+ # field tables. We need to know which field value to pass in to the
+ # field table, and we need to know which table to use for each of the
+ # six channels.
+ cdef int i
cdef np.ndarray[np.float64_t, ndim=1] temp
- temp = tf_obj.red.y
- self.vs[0] = <np.float64_t *> temp.data
- temp = tf_obj.green.y
- self.vs[1] = <np.float64_t *> temp.data
- temp = tf_obj.blue.y
- self.vs[2] = <np.float64_t *> temp.data
- temp = tf_obj.alpha.y
- self.vs[3] = <np.float64_t *> temp.data
- self.x_bounds[0] = tf_obj.x_bounds[0]
- self.x_bounds[1] = tf_obj.x_bounds[1]
- self.nbins = tf_obj.nbins
- self.dbin = (self.x_bounds[1] - self.x_bounds[0])/self.nbins
- self.idbin = 1.0/self.dbin
- self.light_color[0] = tf_obj.light_color[0]
- self.light_color[1] = tf_obj.light_color[1]
- self.light_color[2] = tf_obj.light_color[2]
- self.light_dir[0] = tf_obj.light_dir[0]
- self.light_dir[1] = tf_obj.light_dir[1]
- self.light_dir[2] = tf_obj.light_dir[2]
- cdef np.float64_t normval = 0.0
- for i in range(3): normval += self.light_dir[i]**2
- normval = normval**0.5
- for i in range(3): self.light_dir[i] /= normval
- self.use_light = tf_obj.use_light
+ cdef FieldInterpolationTable fit
- @cython.boundscheck(False)
- @cython.wraparound(False)
- cdef void interpolate(self, np.float64_t dv, np.float64_t *trgba):
- cdef int bin_id, channel
- cdef np.float64_t bv, dy, dd, tf
- bin_id = <int> ((dv - self.x_bounds[0]) * self.idbin)
- # Recall that linear interpolation is y0 + (x-x0) * dx/dy
- dd = dv-(self.x_bounds[0] + bin_id * self.dbin) # x - x0
- if bin_id > self.nbins - 2 or bin_id < 0:
- for channel in range(4): trgba[channel] = 0.0
- return
- for channel in range(4):
- bv = self.vs[channel][bin_id] # This is x0
- dy = self.vs[channel][bin_id+1]-bv # dy
- # This is our final value for transfer function on the entering face
- trgba[channel] = bv+dd*dy*self.idbin
+ self.tf_obj = tf_obj
- @cython.boundscheck(False)
- @cython.wraparound(False)
- cdef void eval_transfer(self, np.float64_t dt, np.float64_t dv,
- np.float64_t *rgba, np.float64_t *grad):
- cdef int i
- cdef np.float64_t ta, tf, trgba[4], dot_prod
- self.interpolate(dv, trgba)
- # get source alpha first
- # First locate our points
- dot_prod = 0.0
- if self.use_light:
- for i in range(3):
- dot_prod += self.light_dir[i] * grad[i]
- dot_prod = fmax(0.0, dot_prod)
- for i in range(3):
- trgba[i] += dot_prod*self.light_color[i]
- # alpha blending
- ta = (1.0 - rgba[3])*dt*trgba[3]
- for i in range(4):
- rgba[i] += ta*trgba[i]
+ self.n_field_tables = tf_obj.n_field_tables
+ for i in range(6): self.istorage[i] = 0.0
+
+ self.my_field_tables = []
+ for i in range(self.n_field_tables):
+ temp = tf_obj.tables[i].y
+ FIT_initialize_table(&self.field_tables[i],
+ temp.shape[0],
+ <np.float64_t *> temp.data,
+ tf_obj.tables[i].x_bounds[0],
+ tf_obj.tables[i].x_bounds[1],
+ tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
+ tf_obj.weight_table_ids[i],
+ tf_obj.tables[i].pass_through)
+ self.my_field_tables.append((tf_obj.tables[i],
+ tf_obj.tables[i].y))
+ self.field_tables[i].field_id = tf_obj.field_ids[i]
+ self.field_tables[i].weight_field_id = tf_obj.weight_field_ids[i]
+ print "Field table", i, "corresponds to",
+ print self.field_tables[i].field_id,
+ print "(Weighted with ", self.field_tables[i].weight_field_id,
+ print ")"
+
+ for i in range(6):
+ self.field_table_ids[i] = tf_obj.field_table_ids[i]
+ print "Channel", i, "corresponds to", self.field_table_ids[i]
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ cdef void eval_transfer(self, np.float64_t dt, np.float64_t *dvs,
+ np.float64_t *rgba, np.float64_t *grad):
+ cdef int i, fid, use
+ cdef np.float64_t ta, tf, trgba[6], dot_prod
+ # This very quick pass doesn't hurt us too badly, and it helps for
+ # early-cutoff. We check all the field tables, because we want to be
+ # able to attenuate even in the presence of no emissivity.
+ use = 0
+ for i in range(self.n_field_tables):
+ fid = self.field_tables[i].field_id
+ if (dvs[fid] >= self.field_tables[i].bounds[0]) and \
+ (dvs[fid] <= self.field_tables[i].bounds[1]):
+ use = 1
+ break
+ for i in range(self.n_field_tables):
+ self.istorage[i] = FIT_get_value(&self.field_tables[i], dvs)
+ # We have to do this after the interpolation
+ for i in range(self.n_field_tables):
+ fid = self.field_tables[i].weight_table_id
+ if fid != -1: self.istorage[i] *= self.istorage[fid]
+ for i in range(6):
+ trgba[i] = self.istorage[self.field_table_ids[i]]
+ #print i, trgba[i],
+ #print
+ # A few words on opacity. We're going to be integrating equation 1.23
+ # from Rybicki & Lightman. dI_\nu / ds = -\alpha_\nu I_\nu + j_\nu
+ # \alpha_nu = \kappa \rho , but we leave that up to the input
+ # transfer function.
+ # SOoooooOOOooo, the upshot is that we are doing a rectangular
+ # integration here:
+ # I_{i+1} = ds * C_i + (1.0 - ds*alpha_i) * I_i
+ for i in range(3):
+ # This is the new way: alpha corresponds to opacity of a given
+ # slice. Previously it was ill-defined, but represented some
+ # measure of emissivity.
+ ta = fmax((1.0 - dt*trgba[i+3]), 0.0)
+ rgba[i ] = dt*trgba[i ] + ta * rgba[i ]
+ #rgba[i+3] = dt*trgba[i+3] + ta * rgba[i+3]
+ # This is the old way:
+ #rgba[i ] += trgba[i] * (1.0 - rgba[i+3])*dt*trgba[i+3]
+ #rgba[i+3] += trgba[i] * (1.0 - rgba[i+3])*dt*trgba[i+3]
cdef class VectorPlane:
cdef public object avp_pos, avp_dir, acenter, aimage
cdef np.float64_t *vp_pos, *vp_dir, *center, *image,
cdef np.float64_t pdx, pdy, bounds[4]
cdef int nv[2]
+ cdef int vp_strides[3]
+ cdef int im_strides[3]
cdef public object ax_vec, ay_vec
cdef np.float64_t *x_vec, *y_vec
@@ -183,6 +256,9 @@
for i in range(4): self.bounds[i] = bounds[i]
self.pdx = (self.bounds[1] - self.bounds[0])/self.nv[0]
self.pdy = (self.bounds[3] - self.bounds[2])/self.nv[1]
+ for i in range(3):
+ self.vp_strides[i] = vp_pos.strides[i] / 8
+ self.im_strides[i] = image.strides[i] / 8
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -200,41 +276,45 @@
rv[3] = rv[2] + lrint((ex[3] - ex[2])/self.pdy)
cdef inline void copy_into(self, np.float64_t *fv, np.float64_t *tv,
- int i, int j, int nk):
+ int i, int j, int nk, int strides[3]):
# We know the first two dimensions of our from-vector, and our
# to-vector is flat and 'ni' long
cdef int k
+ cdef int offset = strides[0] * i + strides[1] * j
for k in range(nk):
- tv[k] = fv[(((k*self.nv[1])+j)*self.nv[0]+i)]
+ tv[k] = fv[offset + k]
cdef inline void copy_back(self, np.float64_t *fv, np.float64_t *tv,
- int i, int j, int nk):
+ int i, int j, int nk, int strides[3]):
cdef int k
+ cdef int offset = strides[0] * i + strides[1] * j
for k in range(nk):
- tv[(((k*self.nv[1])+j)*self.nv[0]+i)] = fv[k]
+ tv[offset + k] = fv[k]
cdef class PartitionedGrid:
cdef public object my_data
cdef public object LeftEdge
cdef public object RightEdge
- cdef np.float64_t *data
+ cdef np.float64_t *data[6]
+ cdef np.float64_t dvs[6]
cdef np.float64_t left_edge[3]
cdef np.float64_t right_edge[3]
cdef np.float64_t dds[3]
cdef np.float64_t idds[3]
cdef int dims[3]
cdef public int parent_grid_id
+ cdef public int n_fields
@cython.boundscheck(False)
@cython.wraparound(False)
def __cinit__(self,
- int parent_grid_id,
- np.ndarray[np.float64_t, ndim=3] data,
+ int parent_grid_id, int n_fields, data,
np.ndarray[np.float64_t, ndim=1] left_edge,
np.ndarray[np.float64_t, ndim=1] right_edge,
np.ndarray[np.int64_t, ndim=1] dims):
# The data is likely brought in via a slice, so we copy it
cdef int i, j, k, size
+ cdef np.ndarray[np.float64_t, ndim=3] tdata
self.parent_grid_id = parent_grid_id
self.LeftEdge = left_edge
self.RightEdge = right_edge
@@ -245,7 +325,10 @@
self.dds[i] = (self.right_edge[i] - self.left_edge[i])/dims[i]
self.idds[i] = 1.0/self.dds[i]
self.my_data = data
- self.data = <np.float64_t*> data.data
+ self.n_fields = n_fields
+ for i in range(n_fields):
+ tdata = data[i]
+ self.data[i] = <np.float64_t *> tdata.data
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -255,7 +338,7 @@
# like http://courses.csusm.edu/cs697exz/ray_box.htm
cdef int vi, vj, hit, i, ni, nj, nn
cdef int iter[4]
- cdef np.float64_t v_pos[3], v_dir[3], rgba[4], extrema[4]
+ cdef np.float64_t v_pos[3], v_dir[3], rgba[6], extrema[4]
self.calculate_extent(vp, extrema)
vp.get_start_stop(extrema, iter)
iter[0] = iclip(iter[0], 0, vp.nv[0])
@@ -265,10 +348,10 @@
hit = 0
for vi in range(iter[0], iter[1]):
for vj in range(iter[2], iter[3]):
- vp.copy_into(vp.vp_pos, v_pos, vi, vj, 3)
- vp.copy_into(vp.image, rgba, vi, vj, 4)
+ vp.copy_into(vp.vp_pos, v_pos, vi, vj, 3, vp.vp_strides)
+ vp.copy_into(vp.image, rgba, vi, vj, 3, vp.im_strides)
self.integrate_ray(v_pos, vp.vp_dir, rgba, tf)
- vp.copy_back(rgba, vp.image, vi, vj, 4)
+ vp.copy_back(rgba, vp.image, vi, vj, 3, vp.im_strides)
return hit
@cython.boundscheck(False)
@@ -423,12 +506,11 @@
for dti in range(tf.ns):
for i in range(3):
dp[i] += ds[i]
- dv = trilinear_interpolate(self.dims, ci, dp, self.data)
- if (dv < tf.x_bounds[0]) or (dv > tf.x_bounds[1]):
- continue
- if tf.use_light == 1:
- eval_gradient(self.dims, ci, dp, self.data, grad)
- tf.eval_transfer(dt, dv, rgba, grad)
+ for i in range(self.n_fields):
+ self.dvs[i] = trilinear_interpolate(self.dims, ci, dp, self.data[i])
+ #if (dv < tf.x_bounds[0]) or (dv > tf.x_bounds[1]):
+ # continue
+ tf.eval_transfer(dt, self.dvs, rgba, grad)
cdef class GridFace:
cdef int direction
@@ -526,7 +608,6 @@
@cython.wraparound(False)
def get_brick(self, np.ndarray[np.float64_t, ndim=1] grid_left_edge,
np.ndarray[np.float64_t, ndim=1] grid_dds,
- np.ndarray[np.float64_t, ndim=3] data,
child_mask):
# We get passed in the left edge, the dds (which gives dimensions) and
# the data, which is already vertex-centered.
@@ -540,8 +621,8 @@
cdef np.ndarray[np.int64_t, ndim=1] dims = np.empty(3, dtype='int64')
for i in range(3):
dims[i] = idims[i]
- cdef np.ndarray[np.float64_t, ndim=3] new_data
- new_data = data[li[0]:ri[0]+1,li[1]:ri[1]+1,li[2]:ri[2]+1].copy()
- PG = PartitionedGrid(self.parent_grid_id, new_data,
- self.LeftEdge, self.RightEdge, dims)
- return [PG]
+ #cdef np.ndarray[np.float64_t, ndim=3] new_data
+ #new_data = data[li[0]:ri[0]+1,li[1]:ri[1]+1,li[2]:ri[2]+1].copy()
+ #PG = PartitionedGrid(self.parent_grid_id, new_data,
+ # self.LeftEdge, self.RightEdge, dims)
+ return ((li[0], ri[0]), (li[1], ri[1]), (li[2], ri[2]), dims)
Added: trunk/yt/_amr_utils/fortran_reader.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/fortran_reader.pyx Sun Apr 25 10:00:23 2010
@@ -0,0 +1,70 @@
+"""
+Simple readers for fortran unformatted data, specifically for the Tiger code.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2010 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+from stdio cimport fopen, fclose, FILE
+
+cdef extern from "stdio.h":
+ cdef int SEEK_SET
+ cdef int SEEK_CUR
+ cdef int SEEK_END
+ int fseek(FILE *stream, long offset, int whence)
+ size_t fread(void *ptr, size_t size, size_t nmemb, FILE *stream)
+ long ftell(FILE *stream)
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+def read_tiger_section(
+ char *fn,
+ np.ndarray[np.int64_t, ndim=1] slab_start,
+ np.ndarray[np.int64_t, ndim=1] slab_size,
+ np.ndarray[np.int64_t, ndim=1] root_size,
+ int offset = 36):
+ cdef int strides[3]
+ strides[0] = 1
+ strides[1] = root_size[0] * strides[0]
+ strides[2] = strides[1] * root_size[1] + 2
+ cdef np.int64_t i, j, k
+ cdef np.ndarray buffer = np.zeros(slab_size, dtype='float32', order='F')
+ cdef FILE *f = fopen(fn, "rb")
+ #for i in range(3): offset += strides[i] * slab_start[i]
+ cdef np.int64_t pos = 0
+ cdef np.int64_t moff = 0
+ cdef float *data = <float *> buffer.data
+ fseek(f, offset, 0)
+ # If anybody wants to convert this loop to a SEEK_CUR, that'd be great.
+ for i in range(slab_size[2]):
+ for j in range(slab_size[1]):
+ moff = (slab_start[0] ) * strides[0] \
+ + (slab_start[1] + j) * strides[1] \
+ + (slab_start[2] + i) * strides[2]
+ #print offset + 4 * moff, pos
+ fseek(f, offset + 4 * moff, SEEK_SET)
+ fread(<void *> (data + pos), 4, slab_size[0], f)
+ pos += slab_size[0]
+ return buffer
Added: trunk/yt/_amr_utils/png_writer.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/png_writer.pyx Sun Apr 25 10:00:23 2010
@@ -0,0 +1,171 @@
+"""
+A light interface to libpng
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2010 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+from stdio cimport fopen, fclose, FILE
+
+cdef extern from "stdlib.h":
+ # NOTE that size_t might not be int
+ void *alloca(int)
+
+cdef extern from "png.h":
+ ctypedef unsigned long png_uint_32
+ ctypedef long png_int_32
+ ctypedef unsigned short png_uint_16
+ ctypedef short png_int_16
+ ctypedef unsigned char png_byte
+ ctypedef void *png_voidp
+ ctypedef png_byte *png_bytep
+ ctypedef png_uint_32 *png_uint_32p
+ ctypedef png_int_32 *png_int_32p
+ ctypedef png_uint_16 *png_uint_16p
+ ctypedef png_int_16 *png_int_16p
+ ctypedef char *png_charp
+ ctypedef char *png_const_charp
+ ctypedef FILE *png_FILE_p
+
+ ctypedef struct png_struct:
+ pass
+ ctypedef png_struct *png_structp
+
+ ctypedef struct png_info:
+ pass
+ ctypedef png_info *png_infop
+
+ ctypedef struct png_color_8:
+ png_byte red
+ png_byte green
+ png_byte blue
+ png_byte gray
+ png_byte alpha
+ ctypedef png_color_8 *png_color_8p
+
+ cdef png_const_charp PNG_LIBPNG_VER_STRING
+
+ # Note that we don't support error or warning functions
+ png_structp png_create_write_struct(
+ png_const_charp user_png_ver, png_voidp error_ptr,
+ void *error_fn, void *warn_fn)
+
+ png_infop png_create_info_struct(png_structp png_ptr)
+
+ void png_init_io(png_structp png_ptr, png_FILE_p fp)
+
+ void png_set_IHDR(png_structp png_ptr, png_infop info_ptr,
+ png_uint_32 width, png_uint_32 height, int bit_depth,
+ int color_type, int interlace_method, int compression_method,
+ int filter_method)
+
+ cdef int PNG_COLOR_TYPE_RGB_ALPHA, PNG_INTERLACE_NONE
+ cdef int PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE
+
+ void png_set_pHYs(png_structp png_ptr, png_infop info_ptr,
+ png_uint_32 res_x, png_uint_32 res_y, int unit_type)
+
+ cdef int PNG_RESOLUTION_METER
+
+ void png_set_sBIT(png_structp png_ptr, png_infop info_ptr,
+ png_color_8p sig_bit)
+
+ void png_write_info(png_structp png_ptr, png_infop info_ptr)
+ void png_write_image(png_structp png_ptr, png_bytep *image)
+ void png_write_end(png_structp png_ptr, png_infop info_ptr)
+
+ void png_destroy_write_struct(
+ png_structp *png_ptr_ptr, png_infop *info_ptr_ptr)
+
+def write_png(np.ndarray[np.uint8_t, ndim=3] buffer,
+ char *filename, int dpi=100):
+
+ # This is something of a translation of the matplotlib _png module
+ cdef png_byte *pix_buffer = <png_byte *> buffer.data
+ cdef int width = buffer.shape[0]
+ cdef int height = buffer.shape[1]
+
+ cdef FILE* fileobj = fopen(filename, "wb")
+ cdef png_bytep *row_pointers
+ cdef png_structp png_ptr
+ cdef png_infop info_ptr
+
+ cdef png_color_8 sig_bit
+ cdef png_uint_32 row
+
+ row_pointers = <png_bytep *> alloca(sizeof(png_bytep) * height)
+
+ for row in range(height):
+ row_pointers[row] = pix_buffer + row * width * 4
+ png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL)
+ info_ptr = png_create_info_struct(png_ptr)
+
+ # Um we are ignoring setjmp sorry guys
+
+ png_init_io(png_ptr, fileobj)
+
+ png_set_IHDR(png_ptr, info_ptr, width, height, 8,
+ PNG_COLOR_TYPE_RGB_ALPHA, PNG_INTERLACE_NONE,
+ PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE)
+
+ cdef size_t dots_per_meter = <size_t> (dpi / (2.54 / 100.0))
+ png_set_pHYs(png_ptr, info_ptr, dots_per_meter, dots_per_meter,
+ PNG_RESOLUTION_METER)
+
+ sig_bit.gray = 0
+ sig_bit.red = sig_bit.green = sig_bit.blue = sig_bit.alpha = 8
+
+ png_set_sBIT(png_ptr, info_ptr, &sig_bit)
+
+ png_write_info(png_ptr, info_ptr)
+ png_write_image(png_ptr, row_pointers)
+ png_write_end(png_ptr, info_ptr)
+
+ fclose(fileobj)
+ png_destroy_write_struct(&png_ptr, &info_ptr)
+
+def add_points_to_image(
+ np.ndarray[np.uint8_t, ndim=3] buffer,
+ np.ndarray[np.float64_t, ndim=1] px,
+ np.ndarray[np.float64_t, ndim=1] py,
+ np.float64_t pv):
+ cdef int i, j, k, pi
+ cdef int np = px.shape[0]
+ cdef int xs = buffer.shape[0]
+ cdef int ys = buffer.shape[1]
+ cdef int v
+ v = iclip(<int>(pv * 255), 0, 255)
+ print "VALUE CONTRIBUTION", v
+ for pi in range(np):
+ j = <int> (xs * px[pi])
+ i = <int> (ys * py[pi])
+ for k in range(3):
+ buffer[i, j, k] = 0
+ return
+ for i in range(xs):
+ for j in range(ys):
+ for k in range(3):
+ v = buffer[i, j, k]
+ buffer[i, j, k] = iclip(v, 0, 255)
Modified: trunk/yt/amr_utils.pyx
==============================================================================
--- trunk/yt/amr_utils.pyx (original)
+++ trunk/yt/amr_utils.pyx Sun Apr 25 10:00:23 2010
@@ -42,3 +42,5 @@
include "_amr_utils/VolumeIntegrator.pyx"
include "_amr_utils/CICDeposit.pyx"
include "_amr_utils/ContourFinding.pyx"
+include "_amr_utils/png_writer.pyx"
+include "_amr_utils/fortran_reader.pyx"
Modified: trunk/yt/config.py
==============================================================================
--- trunk/yt/config.py (original)
+++ trunk/yt/config.py Sun Apr 25 10:00:23 2010
@@ -54,6 +54,7 @@
'time_functions': 'False',
'LogFile': 'False',
'LogFileName': 'yt.log',
+ 'coloredlogs': 'False',
'suppressStreamLogging': 'False',
'LogLevel': '20',
'timefunctions':'False',
Modified: trunk/yt/extensions/volume_rendering/TransferFunction.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/TransferFunction.py (original)
+++ trunk/yt/extensions/volume_rendering/TransferFunction.py Sun Apr 25 10:00:23 2010
@@ -25,9 +25,12 @@
import numpy as na
from matplotlib.cm import get_cmap
+from yt.funcs import *
+from yt.physical_constants import *
class TransferFunction(object):
def __init__(self, x_bounds, nbins=256):
+ self.pass_through = 0
self.nbins = nbins
self.x_bounds = x_bounds
self.x = na.linspace(x_bounds[0], x_bounds[1], nbins).astype('float64')
@@ -51,6 +54,25 @@
vals[(self.x >= start) & (self.x <= stop)] = value
self.y = na.clip(na.maximum(vals, self.y), 0.0, 1.0)
+ def add_filtered_planck(self, wavelength, trans):
+ vals = na.zeros(self.x.shape, 'float64')
+ nu = clight/(wavelength*1e-8)
+ nu = nu[::-1]
+
+ for i,logT in enumerate(self.x):
+ T = 10**logT
+ # Black body at this nu, T
+ Bnu = ((2.0 * hcgs * nu**3) / clight**2.0) / \
+ (na.exp(hcgs * nu / (kboltz * T)) - 1.0)
+ # transmission
+ f = Bnu * trans[::-1]
+ # integrate transmission over nu
+ vals[i] = na.trapz(f,nu)
+
+ # normalize by total transmission over filter
+ self.y = vals/trans.sum() #/na.trapz(trans[::-1],nu)
+ #self.y = na.clip(na.maximum(vals, self.y), 0.0, 1.0)
+
def plot(self, filename):
import matplotlib;matplotlib.use("Agg");import pylab
pylab.clf()
@@ -59,18 +81,48 @@
pylab.ylim(0.0, 1.0)
pylab.savefig(filename)
-class ColorTransferFunction(object):
+class MultiVariateTransferFunction(object):
+ def __init__(self):
+ self.n_field_tables = 0
+ self.tables = [] # Tables are interpolation tables
+ self.field_ids = [0] * 6 # This correlates fields with tables
+ self.weight_field_ids = [-1] * 6 # This correlates
+ self.field_table_ids = [0] * 6
+ self.weight_table_ids = [-1] * 6
+
+ def add_field_table(self, table, field_id, weight_field_id = -1,
+ weight_table_id = -1):
+ self.tables.append(table)
+ self.field_ids[self.n_field_tables] = field_id
+ self.weight_field_ids[self.n_field_tables] = weight_field_id
+ self.weight_table_ids[self.n_field_tables] = weight_table_id
+ self.n_field_tables += 1
+
+ def link_channels(self, table_id, channels = 0):
+ channels = ensure_list(channels)
+ for c in channels:
+ self.field_table_ids[c] = table_id
+
+class ColorTransferFunction(MultiVariateTransferFunction):
def __init__(self, x_bounds, nbins=256):
+ MultiVariateTransferFunction.__init__(self)
self.x_bounds = x_bounds
self.nbins = nbins
+ # This is all compatibility and convenience.
self.red = TransferFunction(x_bounds, nbins)
self.green = TransferFunction(x_bounds, nbins)
self.blue = TransferFunction(x_bounds, nbins)
self.alpha = TransferFunction(x_bounds, nbins)
self.funcs = (self.red, self.green, self.blue, self.alpha)
- self.light_dir = (0.3,-0.2,0.5)
- self.light_color = (0.10, 0.10, 0.10)
- self.use_light = 0
+
+ # Now we do the multivariate stuff
+ # We assign to Density, but do not weight
+ for i,tf in enumerate(self.funcs[:3]):
+ self.add_field_table(tf, 0, weight_table_id = 3)
+ self.link_channels(i, i)
+ self.add_field_table(self.funcs[3], 0)
+ # We don't have a fifth table, so the value will *always* be zero.
+ self.link_channels(4, [3,4,5])
def add_gaussian(self, location, width, height):
for tf, v in zip(self.funcs, height):
@@ -126,6 +178,65 @@
for v, a in zip(na.mgrid[mi:ma:N*1j], alpha):
self.sample_colormap(v, w, a, colormap=colormap)
+class ProjectionTransferFunction(MultiVariateTransferFunction):
+ def __init__(self, x_bounds = (-1e30, 1e30)):
+ MultiVariateTransferFunction.__init__(self)
+ self.x_bounds = x_bounds
+ self.nbins = 2
+ self.linear_mapping = TransferFunction(x_bounds, 2)
+ self.linear_mapping.pass_through = 1
+ self.add_field_table(self.linear_mapping, 0)
+ self.alpha = TransferFunction(x_bounds, 2)
+ self.alpha.y *= 0.0
+ self.alpha.y += 1.0
+ self.add_field_table(self.alpha, 0)
+ self.link_channels(0, [0,1,2]) # same emission for all rgb
+ self.link_channels(2, [3,4,5]) # this will remove absorption
+
+class PlanckTransferFunction(MultiVariateTransferFunction):
+ def __init__(self, T_bounds, rho_bounds, nbins=256,
+ red='R', green='V', blue='B',
+ anorm = 1e6):
+ """
+ This sets up a planck function for multivariate emission and
+ absorption. We assume that the emission is black body, which is then
+ convolved with appropriate Johnson filters for *red*, *green* and
+ *blue*. *T_bounds* and *rho_bounds* define the limits of tabulated
+ emission and absorption functions. *anorm* is a "fudge factor" that
+ defines the somewhat arbitrary normalization to the scattering
+ approximation: because everything is done largely unit-free, and is
+ really not terribly accurate anyway, feel free to adjust this to change
+ the relative amount of reddenning. Maybe in some future version this
+ will be unitful.
+ """
+ MultiVariateTransferFunction.__init__(self)
+ mscat = -1
+ from UBVRI import johnson_filters
+ for i, f in enumerate([red, green, blue]):
+ jf = johnson_filters[f]
+ tf = TransferFunction(T_bounds)
+ tf.add_filtered_planck(jf['wavelen'], jf['trans'])
+ self.add_field_table(tf, 0, 1)
+ self.link_channels(i, i) # 0 => 0, 1 => 1, 2 => 2
+ mscat = max(mscat, jf["Lchar"]**-4)
+
+ for i, f in enumerate([red, green, blue]):
+ # Now we set up the scattering
+ scat = (johnson_filters[f]["Lchar"]**-4 / mscat)*anorm
+ tf = TransferFunction(rho_bounds)
+ print "Adding: %s with relative scattering %s" % (f, scat)
+ tf.y *= 0.0; tf.y += scat
+ self.add_field_table(tf, 1, weight_field_id = 1)
+ self.link_channels(i+3, i+3)
+
+ self._normalize()
+
+ def _normalize(self):
+ fmax = na.array([f.y for f in self.tables[:3]])
+ normal = fmax.max(axis=0)
+ for f in self.tables[:3]:
+ f.y = f.y/normal
+
if __name__ == "__main__":
tf = ColorTransferFunction((-20, -5))
tf.add_gaussian(-16.0, 0.4, [0.2, 0.3, 0.1])
Modified: trunk/yt/extensions/volume_rendering/__init__.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/__init__.py (original)
+++ trunk/yt/extensions/volume_rendering/__init__.py Sun Apr 25 10:00:23 2010
@@ -25,7 +25,10 @@
import numpy as na
-from TransferFunction import TransferFunction, ColorTransferFunction
+from TransferFunction import TransferFunction, ColorTransferFunction, \
+ PlanckTransferFunction, \
+ MultiVariateTransferFunction, \
+ ProjectionTransferFunction
from yt.amr_utils import PartitionedGrid, VectorPlane, \
TransferFunctionProxy
from grid_partitioner import HomogenizedBrickCollection, \
Modified: trunk/yt/extensions/volume_rendering/grid_partitioner.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/grid_partitioner.py (original)
+++ trunk/yt/extensions/volume_rendering/grid_partitioner.py Sun Apr 25 10:00:23 2010
@@ -54,19 +54,33 @@
def write_hierarchy(self, base_filename):
pass
- def _partition_grid(self, grid, field, log_field = True):
- vcd = grid.get_vertex_centered_data(field).astype('float64')
- if log_field: vcd = na.log10(vcd)
+ def _partition_grid(self, grid, fields, log_field = None):
+ fields = ensure_list(fields)
+ if log_field is None: log_field = [True] * len(fields)
+
+ # This is not super efficient, as it re-fills the regions once for each
+ # field.
+ vcds = []
+ for i,field in enumerate(fields):
+ vcd = grid.get_vertex_centered_data(field).astype('float64')
+ if log_field[i]: vcd = na.log10(vcd)
+ vcds.append(vcd)
GF = GridFaces(grid.Children + [grid])
PP = ProtoPrism(grid.id, grid.LeftEdge, grid.RightEdge, GF)
pgs = []
for P in PP.sweep(0):
- pgs += P.get_brick(grid.LeftEdge, grid.dds, vcd, grid.child_mask)
+ sl = P.get_brick(grid.LeftEdge, grid.dds, grid.child_mask)
+ if len(sl) == 0: continue
+ dd = [d[sl[0][0]:sl[0][1]+1,
+ sl[1][0]:sl[1][1]+1,
+ sl[2][0]:sl[2][1]+1].copy() for d in vcds]
+ pgs.append(PartitionedGrid(grid.id, len(fields), dd,
+ P.LeftEdge, P.RightEdge, sl[-1]))
return pgs
- def _partition_local_grids(self, fields = "Density", log_field = True):
+ def _partition_local_grids(self, fields = "Density", log_field = None):
fields = ensure_list(fields)
bricks = []
# We preload.
@@ -80,7 +94,7 @@
print "THIS MANY GRIDS!", len(grid_list)
for i, g in enumerate(self.source._grids):
pbar.update(i)
- bricks += self._partition_grid(g, fields[0], log_field)
+ bricks += self._partition_grid(g, fields, log_field)
pbar.finish()
bricks = na.array(bricks, dtype='object')
NB = len(bricks)
@@ -98,7 +112,7 @@
self.brick_left_edges[i,:] = b.LeftEdge
self.brick_right_edges[i,:] = b.RightEdge
self.brick_parents[i] = b.parent_grid_id
- self.brick_dimensions[i,:] = b.my_data.shape
+ self.brick_dimensions[i,:] = b.my_data[0].shape
# Vertex-centered means we subtract one from the shape
self.brick_dimensions -= 1
self.bricks = na.array(bricks, dtype='object')
Modified: trunk/yt/extensions/volume_rendering/image_handling.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/image_handling.py (original)
+++ trunk/yt/extensions/volume_rendering/image_handling.py Sun Apr 25 10:00:23 2010
@@ -152,7 +152,7 @@
def plot_rgb(image, name, label=None, label_color='w', label_size='large'):
Nvec = image.shape[0]
image[na.isnan(image)] = 0.0
- if image.shape[2] == 4:
+ if image.shape[2] >= 4:
image = image[:,:,:3]
pylab.clf()
pylab.gcf().set_dpi(100)
Modified: trunk/yt/extensions/volume_rendering/software_sampler.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/software_sampler.py (original)
+++ trunk/yt/extensions/volume_rendering/software_sampler.py Sun Apr 25 10:00:23 2010
@@ -98,7 +98,7 @@
for b in self.bricks:
LE.append(b.LeftEdge)
RE.append(b.RightEdge)
- total_cells += na.prod(b.my_data.shape)
+ total_cells += na.prod(b.my_data[0].shape)
LE = na.array(LE) - self.back_center
RE = na.array(RE) - self.back_center
LE = na.sum(LE * self.unit_vectors[2], axis=1)
@@ -111,7 +111,7 @@
tfp.ns = self.sub_samples
for i, b in enumerate(self.bricks[ind]):
pos = b.cast_plane(tfp, self.vector_plane)
- total_cells += na.prod(b.my_data.shape)
+ total_cells += na.prod(b.my_data[0].shape)
pbar.update(total_cells)
pbar.finish()
if finalize: self._finalize()
@@ -145,8 +145,10 @@
plot_rgb(image, prefix)
def partition_grids(self):
- log_field = (self.fields[0] in self.pf.field_info and
- self.pf.field_info[self.fields[0]].take_log)
+ log_field = []
+ for field in self.fields:
+ log_field.append(field in self.pf.field_info and
+ self.pf.field_info[field].take_log)
self._brick_collection._partition_local_grids(self.fields, log_field)
# UNCOMMENT FOR PARALLELISM
#self._brick_collection._collect_bricks(self.source)
@@ -157,7 +159,7 @@
ry = self.resolution[1] * self.res_fac[1]
# We should move away from pre-generation of vectors like this and into
# the usage of on-the-fly generation in the VolumeIntegrator module
- self.image = na.zeros((rx,ry,4), dtype='float64', order='F')
+ self.image = na.zeros((rx,ry,6), dtype='float64', order='C')
# We might have a different width and back_center
bl = self.source.box_lengths
px = na.linspace(-bl[0]/2.0, bl[0]/2.0, rx)[:,None]
@@ -166,7 +168,7 @@
bc = self.source.origin + 0.5*self.source.box_vectors[0] \
+ 0.5*self.source.box_vectors[1]
vectors = na.zeros((rx, ry, 3),
- dtype='float64', order='F')
+ dtype='float64', order='C')
vectors[:,:,0] = inv_mat[0,0]*px + inv_mat[0,1]*py + bc[0]
vectors[:,:,1] = inv_mat[1,0]*px + inv_mat[1,1]*py + bc[1]
vectors[:,:,2] = inv_mat[2,0]*px + inv_mat[2,1]*py + bc[2]
Modified: trunk/yt/funcs.py
==============================================================================
--- trunk/yt/funcs.py (original)
+++ trunk/yt/funcs.py Sun Apr 25 10:00:23 2010
@@ -24,7 +24,7 @@
"""
import time, types, signal, inspect, traceback, sys, pdb, rpdb, os
-import warnings
+import warnings, struct
import progressbar as pb
from math import floor, ceil
from yt.logger import ytLogger as mylog
@@ -82,6 +82,13 @@
return [obj]
return obj
+def read_struct(f, fmt):
+ """
+ This reads a struct, and only that struct, from an open file.
+ """
+ s = f.read(struct.calcsize(fmt))
+ return struct.unpack(fmt, s)
+
def just_one(obj):
# If we have an iterable, sometimes we only want one item
if hasattr(obj,'flat'):
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Sun Apr 25 10:00:23 2010
@@ -1508,15 +1508,16 @@
self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)][::-1]
def _generate_coords(self):
- xi, yi, zi = self.left_edge + self.dds*0.5
- xf, yf, zf = self.left_edge + self.dds*(self.ActiveDimensions-0.5)
- coords = na.mgrid[xi:xf:self.ActiveDimensions[0]*1j,
- yi:yf:self.ActiveDimensions[1]*1j,
- zi:zf:self.ActiveDimensions[2]*1j]
xax = x_dict[self.axis]
yax = y_dict[self.axis]
- self['px'] = coords[xax]
- self['py'] = coords[yax]
+ ci = self.left_edge + self.dds*0.5
+ cf = self.left_edge + self.dds*(self.ActiveDimensions-0.5)
+ cx = na.mgrid[ci[xax]:cf[xax]:self.ActiveDimensions[xax]*1j]
+ cy = na.mgrid[ci[yax]:cf[yax]:self.ActiveDimensions[yax]*1j]
+ blank = na.ones( (self.ActiveDimensions[xax],
+ self.ActiveDimensions[yax]), dtype='float64')
+ self['px'] = cx[None,:] * blank
+ self['py'] = cx[:,None] * blank
self['pdx'] = self.dds[xax]
self['pdy'] = self.dds[yax]
@@ -1536,10 +1537,11 @@
for field in fields_to_get:
self[field] = na.zeros(self.dims, dtype='float64')
dls = self.__setup_dls(fields_to_get)
- for grid in self._get_grids():
+ for i,grid in enumerate(self._get_grids()):
+ mylog.debug("Getting fields from %s", i)
self._get_data_from_grid(grid, fields_to_get, dls)
for field in fields_to_get:
- self[field] = self._mpi_allsum(self[field])
+ self[field] = self._mpi_Allsum_double(self[field])
conv = self.pf.units[self.pf.field_info[field].projection_conversion]
self[field] *= conv
Modified: trunk/yt/lagos/BaseGridType.py
==============================================================================
--- trunk/yt/lagos/BaseGridType.py (original)
+++ trunk/yt/lagos/BaseGridType.py Sun Apr 25 10:00:23 2010
@@ -669,3 +669,27 @@
self.Parent = []
self.Children = []
self.Level = level
+
+class TigerGrid(AMRGridPatch):
+ _id_offset = 0
+
+ def __init__(self, id, hierarchy, left_edge, right_edge, left_dims, right_dims):
+ AMRGridPatch.__init__(self, id, hierarchy = hierarchy)
+ self.LeftEdge = left_edge
+ self.RightEdge = right_edge
+ self.Level = 0
+ self.NumberOfParticles = 0
+ self.left_dims = na.array(left_dims, dtype='int32')
+ self.right_dims = na.array(right_dims, dtype='int32')
+ self.ActiveDimensions = self.right_dims - self.left_dims
+
+ self.Parent = None
+ self.Children = []
+
+ @property
+ def child_mask(self):
+ return na.ones(self.ActiveDimensions, dtype='int32')
+
+ def __repr__(self):
+ return "TigerGrid_%04i (%s)" % (self.id, self.ActiveDimensions)
+
Modified: trunk/yt/lagos/DataReadingFuncs.py
==============================================================================
--- trunk/yt/lagos/DataReadingFuncs.py (original)
+++ trunk/yt/lagos/DataReadingFuncs.py Sun Apr 25 10:00:23 2010
@@ -23,6 +23,7 @@
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
from yt.lagos import *
+import yt.amr_utils as au
import exceptions
_axis_ids = {0:2,1:1,2:0}
@@ -473,3 +474,28 @@
sl[axis] = slice(coord, coord + 1)
return self._read_data_set(grid,field)[sl]
+class IOHandlerTiger(BaseIOHandler):
+ _data_style = "tiger"
+ _offset = 36
+
+ def __init__(self, *args, **kwargs):
+ BaseIOHandler.__init__(self, *args, **kwargs)
+ self._memmaps = {}
+
+ def _read_data_set(self, grid, field):
+ fn = grid.pf.basename + grid.hierarchy.file_mapping[field]
+ LD = na.array(grid.left_dims, dtype='int64')
+ SS = na.array(grid.ActiveDimensions, dtype='int64')
+ RS = na.array(grid.pf.root_size, dtype='int64')
+ data = au.read_tiger_section(fn, LD, SS, RS).astype("float64")
+ return data
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ fn = grid.pf.basename + grid.hierarchy.file_mapping[field]
+ LD = na.array(grid.left_dims, dtype='int64').copy()
+ SS = na.array(grid.ActiveDimensions, dtype='int64').copy()
+ RS = na.array(grid.pf.root_size, dtype='int64').copy()
+ LD[axis] += coord
+ SS[axis] = 1
+ data = au.read_tiger_section(fn, LD, SS, RS).astype("float64")
+ return data
Modified: trunk/yt/lagos/FieldInfoContainer.py
==============================================================================
--- trunk/yt/lagos/FieldInfoContainer.py (original)
+++ trunk/yt/lagos/FieldInfoContainer.py Sun Apr 25 10:00:23 2010
@@ -123,6 +123,15 @@
GadgetFieldInfo = GadgetFieldContainer()
add_gadget_field = GadgetFieldInfo.add_field
+class TigerFieldContainer(CodeFieldInfoContainer):
+ """
+ This is a container for Tiger-specific fields.
+ """
+ _shared_state = {}
+ _field_list = {}
+TigerFieldInfo = TigerFieldContainer()
+add_tiger_field = TigerFieldInfo.add_field
+
class ValidationException(Exception):
pass
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Sun Apr 25 10:00:23 2010
@@ -1325,3 +1325,79 @@
mask[grid_ind] = True
return [g for g in self.grids[mask] if g.Level == grid.Level + 1]
+class TigerHierarchy(AMRHierarchy):
+
+ grid = TigerGrid
+
+ def __init__(self, pf, data_style):
+ self.directory = pf.fullpath
+ self.data_style = data_style
+ AMRHierarchy.__init__(self, pf, data_style)
+
+ def _count_grids(self):
+ # Tiger is unigrid
+ self.ngdims = [i/j for i,j in
+ izip(self.pf.root_size, self.pf.max_grid_size)]
+ self.num_grids = na.prod(self.ngdims)
+ self.max_level = 0
+
+ def _setup_classes(self):
+ dd = self._get_data_reader_dict()
+ AMRHierarchy._setup_classes(self, dd)
+ self.object_types.sort()
+
+ def _parse_hierarchy(self):
+ grids = []
+ # We need to fill in dims, LE, RE, level, count
+ dims, LE, RE, levels, counts = [], [], [], [], []
+ DLE = self.pf["DomainLeftEdge"]
+ DRE = self.pf["DomainRightEdge"]
+ DW = DRE - DLE
+ gds = DW / self.ngdims
+ rd = [self.pf.root_size[i]-self.pf.max_grid_size[i] for i in range(3)]
+ glx, gly, glz = na.mgrid[DLE[0]:DRE[0]-gds[0]:self.ngdims[0]*1j,
+ DLE[1]:DRE[1]-gds[1]:self.ngdims[1]*1j,
+ DLE[2]:DRE[2]-gds[2]:self.ngdims[2]*1j]
+ gdx, gdy, gdz = na.mgrid[0:rd[0]:self.ngdims[0]*1j,
+ 0:rd[1]:self.ngdims[1]*1j,
+ 0:rd[2]:self.ngdims[2]*1j]
+ LE, RE, levels, counts = [], [], [], []
+ i = 0
+ for glei, gldi in izip(izip(glx.flat, gly.flat, glz.flat),
+ izip(gdx.flat, gdy.flat, gdz.flat)):
+ gld = na.array(gldi)
+ gle = na.array(glei)
+ gre = gle + gds
+ g = self.grid(i, self, gle, gre, gld, gld+self.pf.max_grid_size)
+ grids.append(g)
+ dims.append(self.pf.max_grid_size)
+ LE.append(g.LeftEdge)
+ RE.append(g.RightEdge)
+ levels.append(g.Level)
+ counts.append(g.NumberOfParticles)
+ i += 1
+ self.grids = na.array(grids, dtype='object')
+ self.grid_dimensions[:] = na.array(dims, dtype='int64')
+ self.grid_left_edge[:] = na.array(LE, dtype='float64')
+ self.grid_right_edge[:] = na.array(RE, dtype='float64')
+ self.grid_levels.flat[:] = na.array(levels, dtype='int32')
+ self.grid_particle_count.flat[:] = na.array(counts, dtype='int32')
+
+ def _populate_grid_objects(self):
+ # We don't need to do anything here
+ for g in self.grids: g._setup_dx()
+
+ def _detect_fields(self):
+ self.file_mapping = {"Density" : "rhob",
+ "Temperature" : "temp"}
+
+ @property
+ def field_list(self):
+ return self.file_mapping.keys()
+
+ def _setup_unknown_fields(self):
+ for field in self.field_list:
+ add_tiger_field(field, lambda a, b: None)
+
+ def _setup_derived_fields(self):
+ self.derived_field_list = []
Modified: trunk/yt/lagos/OutputTypes.py
==============================================================================
--- trunk/yt/lagos/OutputTypes.py (original)
+++ trunk/yt/lagos/OutputTypes.py Sun Apr 25 10:00:23 2010
@@ -612,8 +612,10 @@
class GadgetStaticOutput(StaticOutput):
_hierarchy_class = GadgetHierarchy
_fieldinfo_class = GadgetFieldContainer
- def __init__(self, filename, particles, dimensions = None):
+ def __init__(self, filename, particles, dimensions = None,
+ storage_filename = None):
StaticOutput.__init__(self, filename, 'gadget_binary')
+ self.storage_filename = storage_filename
self.particles = particles
if "InitialTime" not in self.parameters:
self.parameters["InitialTime"] = 0.0
@@ -640,8 +642,10 @@
_hierarchy_class = ChomboHierarchy
_fieldinfo_class = ChomboFieldContainer
- def __init__(self,filename,data_style='chombo_hdf5'):
+ def __init__(self, filename, data_style='chombo_hdf5',
+ storage_filename = None):
StaticOutput.__init__(self,filename,data_style)
+ self.storage_filename = storage_filename
self.field_info = self._fieldinfo_class()
# hardcoded for now
@@ -700,4 +704,71 @@
pass
return False
-
+class TigerStaticOutput(StaticOutput):
+ _hierarchy_class = TigerHierarchy
+ _fieldinfo_class = TigerFieldContainer
+
+ def __init__(self, rhobname, root_size, max_grid_size=128,
+ data_style='tiger', storage_filename = None):
+ StaticOutput.__init__(self, rhobname, data_style)
+ self.storage_filename = storage_filename
+ self.basename = rhobname[:-4]
+ if not os.path.exists(self.basename + "rhob"):
+ print "%s doesn't exist, don't know how to handle this!" % (
+ self.basename + "rhob")
+ raise IOError
+ if not iterable(root_size): root_size = (root_size,) * 3
+ self.root_size = root_size
+ if not iterable(max_grid_size): max_grid_size = (max_grid_size,) * 3
+ self.max_grid_size = max_grid_size
+
+ self.field_info = self._fieldinfo_class()
+
+ # We assume that we have basename + "rhob" and basename + "temp"
+ # to get at our various parameters.
+
+ # First we get our our header:
+
+ header = [
+ ('i', 'dummy0'),
+ ('f', 'ZR'),
+ ('f', 'OMEGA0'),
+ ('f', 'FLAM0'),
+ ('f', 'OMEGAB'),
+ ('f', 'H0'),
+ ('f', 'BOXL0'),
+ ('i', 'dummy1'),
+ ]
+
+ h_fmt, h_key = zip(*header)
+ header_string = "".join(h_fmt)
+
+ fs = open(self.basename + "rhob")
+ header_raw = read_struct(fs, header_string)
+ self.parameters.update(dict(zip(h_key, header_raw)))
+
+ if "InitialTime" not in self.parameters:
+ self.parameters["InitialTime"] = 0.0
+ self.parameters["CurrentTimeIdentifier"] = \
+ int(os.stat(self.parameter_filename)[ST_CTIME])
+ self.parameters['TopGridDimensions'] = root_size
+ self.parameters['TopGridRank'] = 3
+ self.units["Density"] = 1.0
+ self.parameters['RefineBy'] = 2
+
+ def _set_units(self):
+ self.parameters["DomainLeftEdge"] = na.zeros(3, dtype='float64')
+ self.parameters["DomainRightEdge"] = na.ones(3, dtype='float64')
+ self.units = {}
+ self.time_units = {}
+ self.time_units['1'] = 1
+ self.units['1'] = 1.0
+ self.units['cm'] = 1.0 # This is just plain false
+ self.units['unitary'] = 1.0 / (self["DomainRightEdge"] - self["DomainLeftEdge"]).max()
+
+ def _parse_parameter_file(self):
+ pass
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ return os.path.exists(args[0] + "rhob")
Modified: trunk/yt/lagos/ParallelTools.py
==============================================================================
--- trunk/yt/lagos/ParallelTools.py (original)
+++ trunk/yt/lagos/ParallelTools.py Sun Apr 25 10:00:23 2010
@@ -848,7 +848,7 @@
if MPI.COMM_WORLD.rank == 0:
cc = MPI.Compute_dims(MPI.COMM_WORLD.size, 2)
nsize = final[0]/cc[0], final[1]/cc[1]
- new_image = na.zeros((final[0], final[1], 4), dtype='float64')
+ new_image = na.zeros((final[0], final[1], 6), dtype='float64')
new_image[0:nsize[0],0:nsize[1],:] = data[:]
for i in range(1,MPI.COMM_WORLD.size):
cy, cx = na.unravel_index(i, cc)
@@ -856,7 +856,7 @@
i, nsize[0]*cx,nsize[0]*(cx+1),
nsize[1]*cy,nsize[1]*(cy+1))
buf = _recv_array(source=i, tag=0).reshape(
- (nsize[0],nsize[1],4))
+ (nsize[0],nsize[1],6))
new_image[nsize[0]*cy:nsize[0]*(cy+1),
nsize[1]*cx:nsize[1]*(cx+1),:] = buf[:]
data = new_image
@@ -1031,6 +1031,51 @@
return data
@parallel_passthrough
+ def _mpi_concatenate_array_on_root_double(self, data):
+ self._barrier()
+ size = 0
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1, MPI.COMM_WORLD.size):
+ size = MPI.COMM_WORLD.recv(source=i, tag=0)
+ new_data = na.empty(size, dtype='float64')
+ MPI.COMM_WORLD.Recv([new_data, MPI.DOUBLE], i, 0)
+ data = na.concatenate((data, new_data))
+ else:
+ MPI.COMM_WORLD.send(data.size, 0, 0)
+ MPI.COMM_WORLD.Send([data, MPI.DOUBLE], 0, 0)
+ return data
+
+ @parallel_passthrough
+ def _mpi_concatenate_array_on_root_int(self, data):
+ self._barrier()
+ size = 0
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1, MPI.COMM_WORLD.size):
+ size = MPI.COMM_WORLD.recv(source=i, tag=0)
+ new_data = na.empty(size, dtype='int32')
+ MPI.COMM_WORLD.Recv([new_data, MPI.INT], i, 0)
+ data = na.concatenate((data, new_data))
+ else:
+ MPI.COMM_WORLD.send(data.size, 0, 0)
+ MPI.COMM_WORLD.Send([data, MPI.INT], 0, 0)
+ return data
+
+ @parallel_passthrough
+ def _mpi_concatenate_array_on_root_long(self, data):
+ self._barrier()
+ size = 0
+ if MPI.COMM_WORLD.rank == 0:
+ for i in range(1, MPI.COMM_WORLD.size):
+ size = MPI.COMM_WORLD.recv(source=i, tag=0)
+ new_data = na.empty(size, dtype='int64')
+ MPI.COMM_WORLD.Recv([new_data, MPI.LONG], i, 0)
+ data = na.concatenate((data, new_data))
+ else:
+ MPI.COMM_WORLD.send(data.size, 0, 0)
+ MPI.COMM_WORLD.Send([data, MPI.LONG], 0, 0)
+ return data
+
+ @parallel_passthrough
def _mpi_minimum_array_long(self, data):
"""
Specifically for parallelHOP. For the identical array on each task,
Modified: trunk/yt/logger.py
==============================================================================
--- trunk/yt/logger.py (original)
+++ trunk/yt/logger.py Sun Apr 25 10:00:23 2010
@@ -28,6 +28,32 @@
import logging.handlers as handlers
from yt.config import ytcfg
+# This next bit is grabbed from:
+# http://stackoverflow.com/questions/384076/how-can-i-make-the-python-logging-output-to-be-colored
+def add_coloring_to_emit_ansi(fn):
+ # add methods we need to the class
+ def new(*args):
+ levelno = args[1].levelno
+ if(levelno>=50):
+ color = '\x1b[31m' # red
+ elif(levelno>=40):
+ color = '\x1b[31m' # red
+ elif(levelno>=30):
+ color = '\x1b[33m' # yellow
+ elif(levelno>=20):
+ color = '\x1b[32m' # green
+ elif(levelno>=10):
+ color = '\x1b[35m' # pink
+ else:
+ color = '\x1b[0m' # normal
+ args[1].msg = color + args[1].msg + '\x1b[0m' # normal
+ #print "after"
+ return fn(*args)
+ return new
+
+if ytcfg.getboolean("yt","coloredlogs"):
+ logging.StreamHandler.emit = add_coloring_to_emit_ansi(logging.StreamHandler.emit)
+
level = min(max(ytcfg.getint("yt", "loglevel"), 0), 50)
fstring = "%(name)-10s %(levelname)-10s %(asctime)s %(message)s"
logging.basicConfig(
Modified: trunk/yt/mods.py
==============================================================================
--- trunk/yt/mods.py (original)
+++ trunk/yt/mods.py Sun Apr 25 10:00:23 2010
@@ -43,6 +43,7 @@
BinnedProfile1D, BinnedProfile2D, BinnedProfile3D, \
derived_field, \
add_field, FieldInfo, EnzoFieldInfo, Enzo2DFieldInfo, OrionFieldInfo, \
+ GadgetFieldInfo, TigerFieldInfo, \
Clump, write_clump_hierarchy, find_clumps, write_clumps, \
OrionStaticOutput, HaloFinder, HOPHaloFinder, FOFHaloFinder, parallelHF, \
axis_names, x_dict, y_dict
Modified: trunk/yt/physical_constants.py
==============================================================================
--- trunk/yt/physical_constants.py (original)
+++ trunk/yt/physical_constants.py Sun Apr 25 10:00:23 2010
@@ -14,7 +14,7 @@
# Physical Constants
boltzmann_constant_cgs = 1.3806504e-16 # erg K^-1
gravitational_constant_cgs = 6.67428e-8 # cm^3 g^-1 s^-2
-
+planck_constant_cgs = 6.62606896e-27 # erg s
rho_crit_now = 1.8788e-29 # g times h^2 (critical mass for closure, Cosmology)
# Misc. Approximations
@@ -49,4 +49,5 @@
mp = mass_hydrogen_cgs
clight = speed_of_light_cgs
kboltz = boltzmann_constant_cgs
+hcgs = planck_constant_cgs
sigma_thompson = cross_section_thompson_cgs
Modified: trunk/yt/raven/Callbacks.py
==============================================================================
--- trunk/yt/raven/Callbacks.py (original)
+++ trunk/yt/raven/Callbacks.py Sun Apr 25 10:00:23 2010
@@ -767,7 +767,7 @@
_type_name = "particles"
region = None
_descriptor = None
- def __init__(self, width, p_size=1.0, col='k', stride=1.0, ptype=None, stars_only=False, dm_only=False):
+ def __init__(self, width, p_size=1.0, col='k', marker='o', stride=1.0, ptype=None, stars_only=False, dm_only=False):
"""
Adds particle positions, based on a thick slab along *axis* with a
*width* along the line of sight. *p_size* controls the number of
@@ -778,6 +778,7 @@
self.width = width
self.p_size = p_size
self.color = col
+ self.marker = marker
self.stride = stride
self.ptype = ptype
self.stars_only = stars_only
@@ -808,7 +809,7 @@
px, py = self.convert_to_pixels(plot,
[reg[field_x][gg][::self.stride],
reg[field_y][gg][::self.stride]])
- plot._axes.scatter(px, py, edgecolors='None',
+ plot._axes.scatter(px, py, edgecolors='None', marker=self.marker,
s=self.p_size, c=self.color)
plot._axes.set_xlim(xx0,xx1)
plot._axes.set_ylim(yy0,yy1)
Modified: trunk/yt/raven/ColorMaps.py
==============================================================================
--- trunk/yt/raven/ColorMaps.py (original)
+++ trunk/yt/raven/ColorMaps.py Sun Apr 25 10:00:23 2010
@@ -101,3 +101,12 @@
(1.0, 0.0, 0.0))}
add_cmap('black_green', cdict)
+
+def _extract_lookup_table(cmap_name):
+ cmap = mcm.get_cmap(cmap_name)
+ if not cmap._isinit: cmap._init()
+ r = cmap._lut[:-3, 0]
+ g = cmap._lut[:-3, 1]
+ b = cmap._lut[:-3, 2]
+ a = na.ones(b.shape)
+ return [r, g, b, a]
Modified: trunk/yt/setup.py
==============================================================================
--- trunk/yt/setup.py (original)
+++ trunk/yt/setup.py Sun Apr 25 10:00:23 2010
@@ -1,5 +1,54 @@
#!/usr/bin/env python
import setuptools
+import os, sys
+
+def check_for_png():
+ # First up: HDF5_DIR in environment
+ if "PNG_DIR" in os.environ:
+ png_dir = os.environ["PNG_DIR"]
+ png_inc = os.path.join(png_dir, "include")
+ png_lib = os.path.join(png_dir, "lib")
+ print "PNG_LOCATION: PNG_DIR: %s, %s" % (png_inc, png_lib)
+ return (png_inc, png_lib)
+ # Next up, we try png.cfg
+ elif os.path.exists("png.cfg"):
+ png_dir = open("png.cfg").read().strip()
+ png_inc = os.path.join(png_dir, "include")
+ png_lib = os.path.join(png_dir, "lib")
+ print "PNG_LOCATION: png.cfg: %s, %s" % (png_inc, png_lib)
+ return (png_inc, png_lib)
+ # Now we see if ctypes can help us:
+ try:
+ import ctypes.util
+ png_libfile = ctypes.util.find_library("png")
+ if png_libfile is not None and os.path.isfile(png_libfile):
+ # Now we've gotten a library, but we'll need to figure out the
+ # includes if this is going to work. It feels like there is a
+ # better way to pull off two directory names.
+ png_dir = os.path.dirname(os.path.dirname(png_libfile))
+ if os.path.isdir(os.path.join(png_dir, "include")) and \
+ os.path.isfile(os.path.join(png_dir, "include", "png.h")):
+ png_inc = os.path.join(png_dir, "include")
+ png_lib = os.path.join(png_dir, "lib")
+ print "PNG_LOCATION: png found in: %s, %s" % (png_inc, png_lib)
+ return png_inc, png_lib
+ except ImportError:
+ pass
+ # X11 is where it's located by default on OSX, although I am slightly
+ # reluctant to link against that one.
+ for png_dir in ["/usr/", "/usr/local/", "/usr/X11/"]:
+ if os.path.isfile(os.path.join(png_dir, "include", "png.h")):
+ if os.path.isdir(os.path.join(png_dir, "include")) and \
+ os.path.isfile(os.path.join(png_dir, "include", "png.h")):
+ png_inc = os.path.join(png_dir, "include")
+ png_lib = os.path.join(png_dir, "lib")
+ print "PNG_LOCATION: png found in: %s, %s" % (png_inc, png_lib)
+ return png_inc, png_lib
+ print "Reading png location from png.cfg failed."
+ print "Please place the base directory of your png install in png.cfg and restart."
+ print "(ex: \"echo '/usr/local/' > png.cfg\" )"
+ sys.exit(1)
+
def configuration(parent_package='', top_path=None):
from numpy.distutils.misc_util import Configuration
@@ -11,10 +60,19 @@
config.add_subpackage('extensions')
config.add_subpackage('parallel_tools')
config.add_subpackage('data_objects')
+ png_inc, png_lib = check_for_png()
+ include_dirs=[png_inc]
+ library_dirs=[png_lib]
+ # Because setjmp.h is included by lots of things, and because libpng hasn't
+ # always properly checked its header files (see
+ # https://bugzilla.redhat.com/show_bug.cgi?id=494579 ) we simply disable
+ # support for setjmp.
config.add_extension("amr_utils",
["yt/amr_utils.c", "yt/_amr_utils/FixedInterpolator.c"],
- include_dirs=["yt/_amr_utils/"],
- libraries=["m"])
+ define_macros=[("PNG_SETJMP_NOT_SUPPORTED", True)],
+ include_dirs=["yt/_amr_utils/", png_inc],
+ library_dirs=[png_lib],
+ libraries=["m", "png"])
config.make_config_py()
config.make_svn_version_py()
return config
More information about the yt-svn
mailing list