[yt-svn] commit/yt: 3 new changesets
Bitbucket
commits-noreply at bitbucket.org
Wed Dec 14 09:40:33 PST 2011
3 new commits in yt:
https://bitbucket.org/yt_analysis/yt/changeset/1b327229d1e6/
changeset: 1b327229d1e6
branch: yt
user: ngoldbaum
date: 2011-12-12 23:49:13
summary: Allows a user to look up the current timestep via pf['timestep'] in a FLASH simulation
affected #: 1 file
diff -r bdd0b1bfc80dd2963f010fcba4cf41384bb392eb -r 1b327229d1e6c547e9d6065fdb38ade25754f098 yt/frontends/flash/data_structures.py
--- a/yt/frontends/flash/data_structures.py
+++ b/yt/frontends/flash/data_structures.py
@@ -316,6 +316,13 @@
self.current_time = \
float(self._find_parameter("real", "time", scalar=True))
+ if self._flash_version == 7:
+ self.parameters['timestep'] = float(
+ self._handle["simulation parameters"]["timestep"])
+ else:
+ self.parameters['timestep'] = \
+ float(self._find_parameter("real", "timestep", scalar=True))
+
try:
use_cosmo = self._find_parameter("logical", "usecosmology")
except:
https://bitbucket.org/yt_analysis/yt/changeset/7b29ef36edf5/
changeset: 7b29ef36edf5
branch: yt
user: MatthewTurk
date: 2011-12-14 18:39:57
summary: Changing "timestep" to "dt" as per discussion:
https://bitbucket.org/yt_analysis/yt/pull-request/46/access-timestep-in-flash-outputs
affected #: 1 file
diff -r 1b327229d1e6c547e9d6065fdb38ade25754f098 -r 7b29ef36edf51be59b711703c15a0a4bfecd2b54 yt/frontends/flash/data_structures.py
--- a/yt/frontends/flash/data_structures.py
+++ b/yt/frontends/flash/data_structures.py
@@ -321,7 +321,7 @@
self._handle["simulation parameters"]["timestep"])
else:
self.parameters['timestep'] = \
- float(self._find_parameter("real", "timestep", scalar=True))
+ float(self._find_parameter("real", "dt", scalar=True))
try:
use_cosmo = self._find_parameter("logical", "usecosmology")
https://bitbucket.org/yt_analysis/yt/changeset/cce024601183/
changeset: cce024601183
branch: yt
user: MatthewTurk
date: 2011-12-14 18:40:17
summary: Merging.
affected #: 6 files
diff -r 7b29ef36edf51be59b711703c15a0a4bfecd2b54 -r cce024601183277b73f32feada172a5a7a05a077 yt/analysis_modules/star_analysis/sfr_spectrum.py
--- a/yt/analysis_modules/star_analysis/sfr_spectrum.py
+++ b/yt/analysis_modules/star_analysis/sfr_spectrum.py
@@ -96,6 +96,8 @@
self._pf.current_redshift) # seconds
# Build the distribution.
self.build_dist()
+ # Attach some convenience arrays.
+ self.attach_arrays()
def build_dist(self):
"""
@@ -127,6 +129,47 @@
# We will want the time taken between bins.
self.time_bins_dt = self.time_bins[1:] - self.time_bins[:-1]
+ def attach_arrays(self):
+ """
+ Attach convenience arrays to the class for easy access.
+ """
+ if self.mode == 'data_source':
+ try:
+ vol = self._data_source.volume('mpc')
+ except AttributeError:
+ # If we're here, this is probably a HOPHalo object, and we
+ # can get the volume this way.
+ ds = self._data_source.get_sphere()
+ vol = ds.volume('mpc')
+ elif self.mode == 'provided':
+ vol = self.volume
+ tc = self._pf["Time"]
+ self.time = []
+ self.lookback_time = []
+ self.redshift = []
+ self.Msol_yr = []
+ self.Msol_yr_vol = []
+ self.Msol = []
+ self.Msol_cumulative = []
+ # Use the center of the time_bin, not the left edge.
+ for i, time in enumerate((self.time_bins[1:] + self.time_bins[:-1])/2.):
+ self.time.append(time * tc / YEAR)
+ self.lookback_time.append((self.time_now - time * tc)/YEAR)
+ self.redshift.append(self.cosm.ComputeRedshiftFromTime(time * tc))
+ self.Msol_yr.append(self.mass_bins[i] / \
+ (self.time_bins_dt[i] * tc / YEAR))
+ self.Msol_yr_vol.append(self.mass_bins[i] / \
+ (self.time_bins_dt[i] * tc / YEAR) / vol)
+ self.Msol.append(self.mass_bins[i])
+ self.Msol_cumulative.append(self.cum_mass_bins[i])
+ self.time = na.array(self.time)
+ self.lookback_time = na.array(self.lookback_time)
+ self.redshift = na.array(self.redshift)
+ self.Msol_yr = na.array(self.Msol_yr)
+ self.Msol_yr_vol = na.array(self.Msol_yr_vol)
+ self.Msol = na.array(self.Msol)
+ self.Msol_cumulative = na.array(self.Msol_cumulative)
+
def write_out(self, name="StarFormationRate.out"):
r"""Write out the star analysis to a text file *name*. The columns are in
order.
@@ -150,31 +193,21 @@
>>> sfr.write_out("stars-SFR.out")
"""
fp = open(name, "w")
- if self.mode == 'data_source':
- try:
- vol = self._data_source.volume('mpc')
- except AttributeError:
- # If we're here, this is probably a HOPHalo object, and we
- # can get the volume this way.
- ds = self._data_source.get_sphere()
- vol = ds.volume('mpc')
- elif self.mode == 'provided':
- vol = self.volume
- tc = self._pf["Time"]
- # Use the center of the time_bin, not the left edge.
fp.write("#time\tlookback\tredshift\tMsol/yr\tMsol/yr/Mpc3\tMsol\tcumMsol\t\n")
- for i, time in enumerate((self.time_bins[1:] + self.time_bins[:-1])/2.):
+ for i, time in enumerate(self.time):
line = "%1.5e %1.5e %1.5e %1.5e %1.5e %1.5e %1.5e\n" % \
- (time * tc / YEAR, # Time
- (self.time_now - time * tc)/YEAR, # Lookback time
- self.cosm.ComputeRedshiftFromTime(time * tc), # Redshift
- self.mass_bins[i] / (self.time_bins_dt[i] * tc / YEAR), # Msol/yr
- self.mass_bins[i] / (self.time_bins_dt[i] * tc / YEAR) / vol, # Msol/yr/vol
- self.mass_bins[i], # Msol in bin
- self.cum_mass_bins[i]) # cumulative
+ (time, # Time
+ self.lookback_time[i], # Lookback time
+ self.redshift[i], # Redshift
+ self.Msol_yr[i], # Msol/yr
+ self.Msol_yr_vol[i], # Msol/yr/vol
+ self.Msol[i], # Msol in bin
+ self.Msol_cumulative[i]) # cumulative
fp.write(line)
fp.close()
+### Begin Synthetic Spectrum Stuff. ####
+
CHABRIER = {
"Z0001" : "bc2003_hr_m22_chab_ssp.ised.h5", #/* 0.5% */
"Z0004" : "bc2003_hr_m32_chab_ssp.ised.h5", #/* 2% */
diff -r 7b29ef36edf51be59b711703c15a0a4bfecd2b54 -r cce024601183277b73f32feada172a5a7a05a077 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -847,6 +847,38 @@
self[field] = temp_data[field]
def to_frb(self, width, resolution, center = None):
+ r"""This function returns a FixedResolutionBuffer generated from this
+ object.
+
+ A FixedResolutionBuffer is an object that accepts a variable-resolution
+ 2D object and transforms it into an NxM bitmap that can be plotted,
+ examined or processed. This is a convenience function to return an FRB
+ directly from an existing 2D data object.
+
+ Parameters
+ ----------
+ width : width specifier
+ This can either be a floating point value, in the native domain
+ units of the simulation, or a tuple of the (value, unit) style.
+ This will be the width of the FRB.
+ resolution : int or tuple of ints
+ The number of pixels on a side of the final FRB.
+ center : array-like of floats, optional
+ The center of the FRB. If not specified, defaults to the center of
+ the current object.
+
+ Returns
+ -------
+ frb : :class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer`
+ A fixed resolution buffer, which can be queried for fields.
+
+ Examples
+ --------
+
+ >>> proj = pf.h.proj(0, "Density")
+ >>> frb = proj.to_frb( (100.0, 'kpc'), 1024)
+ >>> write_image(na.log10(frb["Density"]), 'density_100kpc.png')
+ """
if center is None:
center = self.get_field_parameter("center")
if center is None:
@@ -1272,6 +1304,52 @@
return "%s/c%s_L%s" % \
(self._top_node, cen_name, L_name)
+ def to_frb(self, width, resolution):
+ r"""This function returns an ObliqueFixedResolutionBuffer generated
+ from this object.
+
+ An ObliqueFixedResolutionBuffer is an object that accepts a
+ variable-resolution 2D object and transforms it into an NxM bitmap that
+ can be plotted, examined or processed. This is a convenience function
+ to return an FRB directly from an existing 2D data object. Unlike the
+ corresponding to_frb function for other AMR2DData objects, this does
+ not accept a 'center' parameter as it is assumed to be centered at the
+ center of the cutting plane.
+
+ Parameters
+ ----------
+ width : width specifier
+ This can either be a floating point value, in the native domain
+ units of the simulation, or a tuple of the (value, unit) style.
+ This will be the width of the FRB.
+ resolution : int or tuple of ints
+ The number of pixels on a side of the final FRB.
+
+ Returns
+ -------
+ frb : :class:`~yt.visualization.fixed_resolution.ObliqueFixedResolutionBuffer`
+ A fixed resolution buffer, which can be queried for fields.
+
+ Examples
+ --------
+
+ >>> v, c = pf.h.find_max("Density")
+ >>> sp = pf.h.sphere(c, (100.0, 'au'))
+ >>> L = sp.quantities["AngularMomentumVector"]()
+ >>> cutting = pf.h.cutting(L, c)
+ >>> frb = cutting.to_frb( (1.0, 'pc'), 1024)
+ >>> write_image(na.log10(frb["Density"]), 'density_1pc.png')
+ """
+ if iterable(width):
+ w, u = width
+ width = w/self.pf[u]
+ if not iterable(resolution):
+ resolution = (resolution, resolution)
+ from yt.visualization.fixed_resolution import ObliqueFixedResolutionBuffer
+ bounds = (-width/2.0, width/2.0, -width/2.0, width/2.0)
+ frb = ObliqueFixedResolutionBuffer(self, bounds, resolution)
+ return frb
+
class AMRFixedResCuttingPlaneBase(AMR2DData):
"""
AMRFixedResCuttingPlaneBase is an oblique plane through the data,
diff -r 7b29ef36edf51be59b711703c15a0a4bfecd2b54 -r cce024601183277b73f32feada172a5a7a05a077 yt/mods.py
--- a/yt/mods.py
+++ b/yt/mods.py
@@ -122,6 +122,10 @@
from yt.convenience import all_pfs, max_spheres, load, projload
+# Import some helpful math utilities
+from yt.utilities.math_utils import \
+ ortho_find, quartiles
+
# We load plugins. Keep in mind, this can be fairly dangerous -
# the primary purpose is to allow people to have a set of functions
diff -r 7b29ef36edf51be59b711703c15a0a4bfecd2b54 -r cce024601183277b73f32feada172a5a7a05a077 yt/utilities/_amr_utils/VolumeIntegrator.pyx
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx
@@ -65,6 +65,9 @@
double log2(double x)
long int lrint(double x)
double fabs(double x)
+ double cos(double x)
+ double sin(double x)
+ double asin(double x)
cdef struct Triangle:
Triangle *next
@@ -238,6 +241,33 @@
tr[i] = ipnest
return tr
+def arr_fisheye_vectors(int resolution, np.float64_t fov):
+ # We now follow figures 4-7 of:
+ # http://paulbourke.net/miscellaneous/domefisheye/fisheye/
+ # ...but all in Cython.
+ cdef np.ndarray[np.float64_t, ndim=3] vp
+ cdef int i, j, k
+ cdef np.float64_t r, phi, theta, px, py
+ cdef np.float64_t pi = 3.1415926
+ cdef np.float64_t fov_rad = fov * pi / 180.0
+ vp = np.zeros((resolution, resolution, 3), dtype="float64")
+ for i in range(resolution):
+ px = 2.0 * i / (resolution) - 1.0
+ for j in range(resolution):
+ py = 2.0 * j / (resolution) - 1.0
+ r = (px*px + py*py)**0.5
+ if r == 0.0:
+ phi = 0.0
+ elif px < 0:
+ phi = pi - asin(py / r)
+ else:
+ phi = asin(py / r)
+ theta = r * fov_rad / 2.0
+ vp[i,j,0] = sin(theta) * cos(phi)
+ vp[i,j,1] = sin(theta) * sin(phi)
+ vp[i,j,2] = cos(theta)
+ return vp
+
cdef class star_kdtree_container:
cdef kdtree_utils.kdtree *tree
cdef public np.float64_t sigma
diff -r 7b29ef36edf51be59b711703c15a0a4bfecd2b54 -r cce024601183277b73f32feada172a5a7a05a077 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -515,3 +515,107 @@
vec2 /= norm2
vec3 = na.cross(vec1, vec2)
return vec1, vec2, vec3
+
+def quartiles(a, axis=None, out=None, overwrite_input=False):
+ """
+ Compute the quartile values (25% and 75%) along the specified axis
+ in the same way that the numpy.median calculates the median (50%) value
+ alone a specified axis. Check numpy.median for details, as it is
+ virtually the same algorithm.
+
+ Returns an array of the quartiles of the array elements [lower quartile,
+ upper quartile].
+
+ Parameters
+ ----------
+ a : array_like
+ Input array or object that can be converted to an array.
+ axis : {None, int}, optional
+ Axis along which the quartiles are computed. The default (axis=None)
+ is to compute the quartiles along a flattened version of the array.
+ out : ndarray, optional
+ Alternative output array in which to place the result. It must
+ have the same shape and buffer length as the expected output,
+ but the type (of the output) will be cast if necessary.
+ overwrite_input : {False, True}, optional
+ If True, then allow use of memory of input array (a) for
+ calculations. The input array will be modified by the call to
+ quartiles. This will save memory when you do not need to preserve
+ the contents of the input array. Treat the input as undefined,
+ but it will probably be fully or partially sorted. Default is
+ False. Note that, if `overwrite_input` is True and the input
+ is not already an ndarray, an error will be raised.
+
+ Returns
+ -------
+ quartiles : ndarray
+ A new 2D array holding the result (unless `out` is specified, in
+ which case that array is returned instead). If the input contains
+ integers, or floats of smaller precision than 64, then the output
+ data-type is float64. Otherwise, the output data-type is the same
+ as that of the input.
+
+ See Also
+ --------
+ numpy.median, numpy.mean, numpy.percentile
+
+ Notes
+ -----
+ Given a vector V of length N, the quartiles of V are the 25% and 75% values
+ of a sorted copy of V, ``V_sorted`` - i.e., ``V_sorted[(N-1)/4]`` and
+ ``3*V_sorted[(N-1)/4]``, when N is odd. When N is even, it is the average
+ of the two values bounding these values of ``V_sorted``.
+
+ Examples
+ --------
+ >>> a = na.arange(100).reshape(10,10)
+ >>> a
+ array([[ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
+ [10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
+ [20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
+ [30, 31, 32, 33, 34, 35, 36, 37, 38, 39],
+ [40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
+ [50, 51, 52, 53, 54, 55, 56, 57, 58, 59],
+ [60, 61, 62, 63, 64, 65, 66, 67, 68, 69],
+ [70, 71, 72, 73, 74, 75, 76, 77, 78, 79],
+ [80, 81, 82, 83, 84, 85, 86, 87, 88, 89],
+ [90, 91, 92, 93, 94, 95, 96, 97, 98, 99]])
+ >>> mu.quartiles(a)
+ array([ 24.5, 74.5])
+ >>> mu.quartiles(a,axis=0)
+ array([[ 15., 16., 17., 18., 19., 20., 21., 22., 23., 24.],
+ [ 65., 66., 67., 68., 69., 70., 71., 72., 73., 74.]])
+ >>> mu.quartiles(a,axis=1)
+ array([[ 1.5, 11.5, 21.5, 31.5, 41.5, 51.5, 61.5, 71.5, 81.5,
+ 91.5],
+ [ 6.5, 16.5, 26.5, 36.5, 46.5, 56.5, 66.5, 76.5, 86.5,
+ 96.5]])
+ """
+ if overwrite_input:
+ if axis is None:
+ sorted = a.ravel()
+ sorted.sort()
+ else:
+ a.sort(axis=axis)
+ sorted = a
+ else:
+ sorted = na.sort(a, axis=axis)
+ if axis is None:
+ axis = 0
+ indexer = [slice(None)] * sorted.ndim
+ indices = [int(sorted.shape[axis]/4), int(sorted.shape[axis]*.75)]
+ result = []
+ for index in indices:
+ if sorted.shape[axis] % 2 == 1:
+ # index with slice to allow mean (below) to work
+ indexer[axis] = slice(index, index+1)
+ else:
+ indexer[axis] = slice(index-1, index+1)
+ # special cases for small arrays
+ if sorted.shape[axis] == 2:
+ # index with slice to allow mean (below) to work
+ indexer[axis] = slice(index, index+1)
+ # Use mean in odd and even case to coerce data type
+ # and check, use out array.
+ result.append(na.mean(sorted[indexer], axis=axis, out=out))
+ return na.array(result)
diff -r 7b29ef36edf51be59b711703c15a0a4bfecd2b54 -r cce024601183277b73f32feada172a5a7a05a077 yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -32,7 +32,7 @@
from yt.utilities.amr_utils import TransferFunctionProxy, VectorPlane, \
arr_vec2pix_nest, arr_pix2vec_nest, AdaptiveRaySource, \
- arr_ang2pix_nest
+ arr_ang2pix_nest, arr_fisheye_vectors
from yt.visualization.image_writer import write_bitmap
from yt.data_objects.data_containers import data_object_registry
from yt.utilities.parallel_tools.parallel_analysis_interface import \
@@ -793,6 +793,59 @@
oc.sub_samples, oc.pf)
return (left_camera, right_camera)
+class FisheyeCamera(Camera):
+ def __init__(self, center, radius, fov, resolution,
+ transfer_function = None, fields = None,
+ sub_samples = 5, log_fields = None, volume = None,
+ pf = None, no_ghost=False):
+ ParallelAnalysisInterface.__init__(self)
+ if pf is not None: self.pf = pf
+ self.center = na.array(center, dtype='float64')
+ self.radius = radius
+ self.fov = fov
+ if iterable(resolution):
+ raise RuntimeError("Resolution must be a single int")
+ self.resolution = resolution
+ if transfer_function is None:
+ transfer_function = ProjectionTransferFunction()
+ self.transfer_function = transfer_function
+ if fields is None: fields = ["Density"]
+ self.fields = fields
+ self.sub_samples = sub_samples
+ self.log_fields = log_fields
+ if volume is None:
+ volume = AMRKDTree(self.pf, fields=self.fields, no_ghost=no_ghost,
+ log_fields=log_fields)
+ self.volume = volume
+
+ def snapshot(self):
+ image = na.zeros((self.resolution**2,1,3), dtype='float64', order='C')
+ # We now follow figures 4-7 of:
+ # http://paulbourke.net/miscellaneous/domefisheye/fisheye/
+ # ...but all in Cython.
+ vp = arr_fisheye_vectors(self.resolution, self.fov)
+ vp.shape = (self.resolution**2,1,3)
+ uv = na.ones(3, dtype='float64')
+ positions = na.ones((self.resolution**2, 1, 3), dtype='float64') * self.center
+ vector_plane = VectorPlane(positions, vp, self.center,
+ (0.0, 1.0, 0.0, 1.0), image, uv, uv)
+ tfp = TransferFunctionProxy(self.transfer_function)
+ tfp.ns = self.sub_samples
+ self.volume.initialize_source()
+ mylog.info("Rendering fisheye of %s^2", self.resolution)
+ pbar = get_pbar("Ray casting",
+ (self.volume.brick_dimensions + 1).prod(axis=-1).sum())
+
+ total_cells = 0
+ for brick in self.volume.traverse(None, self.center, image):
+ brick.cast_plane(tfp, vector_plane)
+ total_cells += na.prod(brick.my_data[0].shape)
+ pbar.update(total_cells)
+ pbar.finish()
+ image.shape = (self.resolution, self.resolution, 3)
+ return image
+
+
def off_axis_projection(pf, center, normal_vector, width, resolution,
field, weight = None, volume = None, no_ghost = True):
r"""Project through a parameter file, off-axis, and return the image plane.
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list