[yt-svn] commit/yt: 11 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Jun 19 05:53:07 PDT 2014
11 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/87039982afe8/
Changeset: 87039982afe8
Branch: yt-3.0
User: jzuhone
Date: 2014-05-28 20:31:04
Summary: Merging
Affected #: 22 files
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 doc/source/analyzing/objects.rst
--- a/doc/source/analyzing/objects.rst
+++ b/doc/source/analyzing/objects.rst
@@ -236,7 +236,7 @@
-------------------------------
Data objects can be cut by their field values using the ``cut_region``
-method. For example, this could be used to compute the total mass within
+method. For example, this could be used to compute the total gas mass within
a certain temperature range, as in the following example.
.. notebook-cell::
@@ -244,11 +244,11 @@
from yt.mods import *
ds = load("enzo_tiny_cosmology/DD0046/DD0046")
ad = ds.all_data()
- total_mass = ad.quantities.total_mass()
+ total_mass = ad.quantities.total_quantity('cell_mass')
# now select only gas with 1e5 K < T < 1e7 K.
new_region = ad.cut_region(['obj["temperature"] > 1e5',
'obj["temperature"] < 1e7'])
- cut_mass = new_region.quantities.total_mass()
+ cut_mass = new_region.quantities.total_quantity('cell_mass')
print "The fraction of mass in this temperature range is %f." % \
(cut_mass / total_mass)
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/data_objects/profiles.py
--- a/yt/data_objects/profiles.py
+++ b/yt/data_objects/profiles.py
@@ -1154,7 +1154,10 @@
for bin_field in bin_fields:
bf_units = data_source.pf._get_field_info(bin_field[0],
bin_field[1]).units
- field_ex = list(extrema[bin_field[-1]])
+ try:
+ field_ex = list(extrema[bin_field[-1]])
+ except KeyError:
+ field_ex = list(extrema[bin_field])
if iterable(field_ex[0]):
field_ex[0] = data_source.pf.quan(field_ex[0][0], field_ex[0][1])
field_ex[0] = field_ex[0].in_units(bf_units)
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -24,7 +24,10 @@
YTSelectionContainer1D, YTSelectionContainer2D, YTSelectionContainer3D
from yt.data_objects.derived_quantities import \
DerivedQuantityCollection
-from yt.utilities.exceptions import YTSphereTooSmall
+from yt.utilities.exceptions import \
+ YTSphereTooSmall, \
+ YTIllDefinedCutRegion, \
+ YTMixedCutRegion
from yt.utilities.linear_interpolators import TrilinearFieldInterpolator
from yt.utilities.minimal_representation import \
MinimalSliceData
@@ -683,6 +686,9 @@
self.base_object.get_data(fields)
ind = self._cond_ind
for field in fields:
+ f = self.base_object[field]
+ if f.shape != ind.shape:
+ raise YTMixedCutRegion(self.conditionals, field)
self.field_data[field] = self.base_object[field][ind]
@property
@@ -693,6 +699,8 @@
for cond in self.conditionals:
res = eval(cond)
if ind is None: ind = res
+ if ind.shape != res.shape:
+ raise YTIllDefinedCutRegion(self.conditionals)
np.logical_and(res, ind, ind)
return ind
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -519,16 +519,28 @@
def find_max(self, field):
"""
- Returns (value, center) of location of maximum for a given field.
+ Returns (value, location) of the maximum of a given field.
"""
mylog.debug("Searching for maximum value of %s", field)
source = self.all_data()
max_val, maxi, mx, my, mz = \
- source.quantities["MaxLocation"](field)
+ source.quantities.max_location(field)
mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f",
max_val, mx, my, mz)
return max_val, np.array([mx, my, mz], dtype="float64")
+ def find_min(self, field):
+ """
+ Returns (value, location) for the minimum of a given field.
+ """
+ mylog.debug("Searching for minimum value of %s", field)
+ source = self.all_data()
+ min_val, maxi, mx, my, mz = \
+ source.quantities.min_location(field)
+ mylog.info("Min Value is %0.5e at %0.16f %0.16f %0.16f",
+ min_val, mx, my, mz)
+ return min_val, np.array([mx, my, mz], dtype="float64")
+
# Now all the object related stuff
def all_data(self, find_max=False):
if find_max: c = self.find_max("density")[1]
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -44,6 +44,7 @@
rho_crit_g_cm3_h2, cm_per_mpc
from yt.utilities.io_handler import io_registry
from yt.utilities.logger import ytLogger as mylog
+from yt.utilities.pyparselibconfig import libconfig
from .fields import \
EnzoFieldInfo
@@ -277,7 +278,15 @@
si, ei, LE, RE, fn, npart = [], [], [], [], [], []
all = [si, ei, LE, RE, fn]
pbar = get_pbar("Parsing Hierarchy ", self.num_grids)
- if self.parameter_file.parameters["VersionNumber"] > 2.0:
+ version = self.parameter_file.parameters.get("VersionNumber", None)
+ params = self.parameter_file.parameters
+ if version is None and "Internal" in params:
+ version = float(params["Internal"]["Provenance"]["VersionNumber"])
+ if version >= 3.0:
+ active_particles = True
+ nap = dict((ap_type, []) for ap_type in
+ params["Physics"]["ActiveParticles"]["ActiveParticlesEnabled"])
+ elif version > 2.0:
active_particles = True
nap = {}
for type in self.parameters.get("AppendActiveParticleType", []):
@@ -731,10 +740,53 @@
# Let's read the file
self.unique_identifier = \
int(os.stat(self.parameter_filename)[stat.ST_CTIME])
- data_labels = {}
- data_label_factors = {}
- conversion_factors = {}
- for line in (l.strip() for l in open(self.parameter_filename)):
+ with open(self.parameter_filename, "r") as f:
+ line = f.readline().strip()
+ f.seek(0)
+ if line == "Internal:":
+ self._parse_enzo3_parameter_file(f)
+ else:
+ self._parse_enzo2_parameter_file(f)
+
+ def _parse_enzo3_parameter_file(self, f):
+ self.parameters = p = libconfig(f)
+ sim = p["SimulationControl"]
+ internal = p["Internal"]
+ phys = p["Physics"]
+ self.refine_by = sim["AMR"]["RefineBy"]
+ self.periodicity = tuple(a == 3 for a in
+ sim["Domain"]["LeftFaceBoundaryCondition"])
+ self.dimensionality = sim["Domain"]["TopGridRank"]
+ self.domain_dimensions = np.array(sim["Domain"]["TopGridDimensions"],
+ dtype="int64")
+ self.domain_left_edge = np.array(sim["Domain"]["DomainLeftEdge"],
+ dtype="float64")
+ self.domain_right_edge = np.array(sim["Domain"]["DomainRightEdge"],
+ dtype="float64")
+ self.gamma = phys["Hydro"]["Gamma"]
+ self.unique_identifier = internal["Provenance"]["CurrentTimeIdentifier"]
+ self.current_time = internal["InitialTime"]
+ self.cosmological_simulation = phys["Cosmology"]["ComovingCoordinates"]
+ if self.cosmological_simulation == 1:
+ cosmo = phys["Cosmology"]
+ self.current_redshift = internal["CosmologyCurrentRedshift"]
+ self.omega_lambda = cosmo["OmegaLambdaNow"]
+ self.omega_matter = cosmo["OmegaMatterNow"]
+ self.hubble_constant = cosmo["HubbleConstantNow"]
+ else:
+ self.current_redshift = self.omega_lambda = self.omega_matter = \
+ self.hubble_constant = self.cosmological_simulation = 0.0
+ self.particle_types = ["DarkMatter"] + \
+ phys["ActiveParticles"]["ActiveParticlesEnabled"]
+ self.particle_types = tuple(self.particle_types)
+ self.particle_types_raw = self.particle_types
+ if self.dimensionality == 1:
+ self._setup_1d()
+ elif self.dimensionality == 2:
+ self._setup_2d()
+
+ def _parse_enzo2_parameter_file(self, f):
+ for line in (l.strip() for l in f):
if len(line) < 2: continue
param, vals = (i.strip() for i in line.split("=",1))
# First we try to decipher what type of value it is.
@@ -835,8 +887,11 @@
if self.cosmological_simulation:
k = self.cosmology_get_units()
# Now some CGS values
- self.length_unit = \
- self.quan(self.parameters["CosmologyComovingBoxSize"], "Mpccm/h")
+ box_size = self.parameters.get("CosmologyComovingBoxSize", None)
+ if box_size is None:
+ box_size = self.parameters["Physics"]["Cosmology"]\
+ ["CosmologyComovingBoxSize"]
+ self.length_unit = self.quan(box_size, "Mpccm/h")
self.mass_unit = \
self.quan(k['urho'], 'g/cm**3') * (self.length_unit.in_cgs())**3
self.time_unit = self.quan(k['utim'], 's')
@@ -846,6 +901,11 @@
length_unit = self.parameters["LengthUnits"]
mass_unit = self.parameters["DensityUnits"] * length_unit**3
time_unit = self.parameters["TimeUnits"]
+ elif "SimulationControl" in self.parameters:
+ units = self.parameters["SimulationControl"]["Units"]
+ length_unit = units["Length"]
+ mass_unit = units["Density"] * length_unit**3
+ time_unit = units["Time"]
else:
mylog.warning("Setting 1.0 in code units to be 1.0 cm")
mylog.warning("Setting 1.0 in code units to be 1.0 s")
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -96,7 +96,10 @@
)
def __init__(self, pf, field_list):
- if pf.parameters["HydroMethod"] == 2:
+ hydro_method = pf.parameters.get("HydroMethod", None)
+ if hydro_method is None:
+ hydro_method = pf.parameters["Physics"]["Hydro"]["HydroMethod"]
+ if hydro_method == 2:
sl_left = slice(None,-2,None)
sl_right = slice(1,-1,None)
div_fac = 1.0
@@ -142,7 +145,11 @@
def setup_fluid_fields(self):
# Now we conditionally load a few other things.
- if self.pf.parameters["MultiSpecies"] > 0:
+ params = self.pf.parameters
+ multi_species = params.get("MultiSpecies", None)
+ if multi_species is None:
+ multi_species = params["Physics"]["AtomicPhysics"]["MultiSpecies"]
+ if multi_species > 0:
self.setup_species_fields()
self.setup_energy_field()
@@ -150,6 +157,16 @@
# We check which type of field we need, and then we add it.
ge_name = None
te_name = None
+ params = self.pf.parameters
+ multi_species = params.get("MultiSpecies", None)
+ if multi_species is None:
+ multi_species = params["Physics"]["AtomicPhysics"]["MultiSpecies"]
+ hydro_method = params.get("HydroMethod", None)
+ if hydro_method is None:
+ hydro_method = params["Physics"]["Hydro"]["HydroMethod"]
+ dual_energy = params.get("DualEnergyFormalism", None)
+ if dual_energy is None:
+ dual_energy = params["Physics"]["Hydro"]["DualEnergyFormalism"]
if ("enzo", "Gas_Energy") in self.field_list:
ge_name = "Gas_Energy"
elif ("enzo", "GasEnergy") in self.field_list:
@@ -159,12 +176,12 @@
elif ("enzo", "TotalEnergy") in self.field_list:
te_name = "TotalEnergy"
- if self.pf.parameters["HydroMethod"] == 2:
+ if hydro_method == 2:
self.add_output_field(("enzo", te_name),
units="code_velocity**2")
self.alias(("gas", "thermal_energy"), ("enzo", te_name))
- elif self.pf.parameters["DualEnergyFormalism"] == 1:
+ elif dual_energy == 1:
self.add_output_field(
("enzo", ge_name),
units="code_velocity**2")
@@ -172,7 +189,7 @@
("gas", "thermal_energy"),
("enzo", ge_name),
units = "erg/g")
- elif self.pf.parameters["HydroMethod"] in (4, 6):
+ elif hydro_method in (4, 6):
self.add_output_field(
("enzo", te_name),
units="code_velocity**2")
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -62,6 +62,7 @@
("Temperature", ("K", ["temperature"], None)),
("Epsilon", ("code_length", [], None)),
("Metals", ("code_metallicity", ["metallicity"], None)),
+ ("Metallicity", ("code_metallicity", ["metallicity"], None)),
("Phi", ("code_length", [], None)),
("FormationTime", ("code_time", ["creation_time"], None)),
# These are metallicity fields that get discovered for FIRE simulations
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/frontends/stream/io.py
--- a/yt/frontends/stream/io.py
+++ b/yt/frontends/stream/io.py
@@ -134,24 +134,29 @@
def _initialize_index(self, data_file, regions):
# self.fields[g.id][fname] is the pattern here
- pos = np.column_stack(self.fields[data_file.filename][
- ("io", "particle_position_%s" % ax)]
- for ax in 'xyz')
- if np.any(pos.min(axis=0) < data_file.pf.domain_left_edge) or \
- np.any(pos.max(axis=0) > data_file.pf.domain_right_edge):
- raise YTDomainOverflow(pos.min(axis=0), pos.max(axis=0),
- data_file.pf.domain_left_edge,
- data_file.pf.domain_right_edge)
- regions.add_data_file(pos, data_file.file_id)
- morton = compute_morton(
- pos[:,0], pos[:,1], pos[:,2],
- data_file.pf.domain_left_edge,
- data_file.pf.domain_right_edge)
- return morton
+ morton = []
+ for ptype in self.pf.particle_types_raw:
+ pos = np.column_stack(self.fields[data_file.filename][
+ (ptype, "particle_position_%s" % ax)]
+ for ax in 'xyz')
+ if np.any(pos.min(axis=0) < data_file.pf.domain_left_edge) or \
+ np.any(pos.max(axis=0) > data_file.pf.domain_right_edge):
+ raise YTDomainOverflow(pos.min(axis=0), pos.max(axis=0),
+ data_file.pf.domain_left_edge,
+ data_file.pf.domain_right_edge)
+ regions.add_data_file(pos, data_file.file_id)
+ morton.append(compute_morton(
+ pos[:,0], pos[:,1], pos[:,2],
+ data_file.pf.domain_left_edge,
+ data_file.pf.domain_right_edge))
+ return np.concatenate(morton)
def _count_particles(self, data_file):
- npart = self.fields[data_file.filename]["io", "particle_position_x"].size
- return {'io': npart}
+ pcount = {}
+ for ptype in self.pf.particle_types_raw:
+ d = self.fields[data_file.filename]
+ pcount[ptype] = d[ptype, "particle_position_x"].size
+ return pcount
def _identify_fields(self, data_file):
return self.fields[data_file.filename].keys(), {}
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/funcs.py
--- a/yt/funcs.py
+++ b/yt/funcs.py
@@ -622,6 +622,22 @@
@contextlib.contextmanager
def parallel_profile(prefix):
+ r"""A context manager for profiling parallel code execution using cProfile
+
+ This is a simple context manager that automatically profiles the execution
+ of a snippet of code.
+
+ Parameters
+ ----------
+ prefix : string
+ A string name to prefix outputs with.
+
+ Examples
+ --------
+
+ >>> with parallel_profile('my_profile'):
+ ... yt.PhasePlot(ds.all_data(), 'density', 'temperature', 'cell_mass')
+ """
import cProfile
from yt.config import ytcfg
fn = "%s_%04i_%04i.cprof" % (prefix,
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -627,8 +627,9 @@
domain_width = DW[i]
if region_width <= 0:
- print "Error: region right edge < left edge", region_width
- raise RuntimeError
+ raise RuntimeError(
+ "Region right edge < left edge: width = %s" % region_width
+ )
if dobj.pf.periodicity[i]:
# shift so left_edge guaranteed in domain
@@ -641,10 +642,13 @@
else:
if dobj.left_edge[i] < dobj.pf.domain_left_edge[i] or \
dobj.right_edge[i] > dobj.pf.domain_right_edge[i]:
- print "Error: bad Region in non-periodic domain:", dobj.left_edge[i], \
- dobj.pf.domain_left_edge[i], dobj.right_edge[i], dobj.pf.domain_right_edge[i]
- raise RuntimeError
-
+ raise RuntimeError(
+ "Error: bad Region in non-periodic domain along dimension %s. "
+ "Region left edge = %s, Region right edge = %s"
+ "Dataset left edge = %s, Dataset right edge = %s" % \
+ (i, dobj.left_edge[i], dobj.right_edge[i],
+ dobj.pf.domain_left_edge[i], dobj.pf.domain_right_edge[i])
+ )
# Already ensured in code
self.left_edge[i] = dobj.left_edge[i]
self.right_edge[i] = dobj.right_edge[i]
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/units/yt_array.py
--- a/yt/units/yt_array.py
+++ b/yt/units/yt_array.py
@@ -252,7 +252,7 @@
"Perhaps you meant to do something like this instead: \n"
"ds.arr(%s, \"%s\")" % (input_array, input_units)
)
- if _astropy.units is not None:
+ if _astropy._units is not None:
if isinstance(input_array, _astropy.units.quantity.Quantity):
return cls.from_astropy(input_array)
if isinstance(input_array, YTArray):
@@ -1029,9 +1029,9 @@
def get_binary_op_return_class(cls1, cls2):
if cls1 is cls2:
return cls1
- if cls1 is np.ndarray or issubclass(cls1, numeric_type):
+ if cls1 is np.ndarray or issubclass(cls1, (numeric_type, np.number)):
return cls2
- if cls2 is np.ndarray or issubclass(cls2, numeric_type):
+ if cls2 is np.ndarray or issubclass(cls2, (numeric_type, np.number)):
return cls1
if issubclass(cls1, YTQuantity):
return cls2
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/utilities/exceptions.py
--- a/yt/utilities/exceptions.py
+++ b/yt/utilities/exceptions.py
@@ -377,3 +377,26 @@
r = """Position arrays must be length and shape (N,3).
But this one has %s and %s.""" % (self.dimensions, self.shape)
return r
+
+class YTIllDefinedCutRegion(Exception):
+ def __init__(self, conditions):
+ self.conditions = conditions
+
+ def __str__(self):
+ r = """Can't mix particle/discrete and fluid/mesh conditions or
+ quantities. Conditions specified:
+ """
+ r += "\n".join([c for c in self.conditions])
+ return r
+
+class YTMixedCutRegion(Exception):
+ def __init__(self, conditions, field):
+ self.conditions = conditions
+ self.field = field
+
+ def __str__(self):
+ r = """Can't mix particle/discrete and fluid/mesh conditions or
+ quantities. Field: %s and Conditions specified:
+ """ % (self.field,)
+ r += "\n".join([c for c in self.conditions])
+ return r
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/utilities/pyparselibconfig/__init__.py
--- /dev/null
+++ b/yt/utilities/pyparselibconfig/__init__.py
@@ -0,0 +1,2 @@
+from api import libconfig
+
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/utilities/pyparselibconfig/api.py
--- /dev/null
+++ b/yt/utilities/pyparselibconfig/api.py
@@ -0,0 +1,1 @@
+from libconfig import libconfig
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/utilities/pyparselibconfig/libconfig.py
--- /dev/null
+++ b/yt/utilities/pyparselibconfig/libconfig.py
@@ -0,0 +1,102 @@
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, Samuel Skillman
+#
+# Distributed under the terms of the MIT License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+class libconfig(dict):
+ def __init__(self, config=None):
+ if config is not None:
+ self.read(config)
+
+ def read(self, config):
+ if not hasattr(config, "read"):
+ cfile = open(config, 'r')
+ else:
+ cfile = config
+
+ # Strip out spaces and blanks
+ lines = [line.strip() for line in cfile if len(line) > 0]
+
+ # Strip out comments
+ lines = [line for line in lines if not (line.startswith('#') or
+ line.startswith('/'))]
+
+ # Concatenate
+ oneline = ''
+ for line in lines:
+ oneline += line
+
+ statements = oneline.split(';')
+
+ self.parse_statements(self, statements)
+
+ def parse_statements(self, this_dict, statements):
+ while len(statements) > 0:
+ statement = statements.pop(0)
+ if len(statement) == 0:
+ continue
+ if ':' in statement:
+ # DICTIONARY
+ new_level_lines = []
+ k, v = statement.split(':', 1)
+ k = k.strip(' :')
+ v = v.strip(' {')
+ level = 1 + v.count(':')
+ this_dict[k] = {}
+ new_level_lines.append(v)
+
+ while(level > 0 and len(statements) > 0):
+ nextline = statements.pop(0).strip()
+ level += nextline.count('{')
+ if nextline == '}' and level == 1:
+ level = 0
+ break
+ new_level_lines.append(nextline)
+ level -= nextline.count('}')
+ self.parse_statements(this_dict[k], new_level_lines)
+ else:
+ k, v = statement.split('=')
+ k = k.strip()
+ v = v.strip()
+ this_dict[k] = self.correct_type(v)
+
+ def correct_type(self, v):
+ if v == "true":
+ v = "True"
+ elif v == "false":
+ v = "False"
+ # ...are we really evaling this? We should work around that somehow.
+ return eval(v)
+
+ def write(self, filename):
+ f = file(filename, 'w')
+
+ self.write_dict(f, self, 0)
+
+ def write_dict(self, f, this_dict, level):
+ tab = ' '*4
+
+ dict_dict = {}
+ for k, v in this_dict.iteritems():
+ if type(v) == dict:
+ dict_dict[k] = v
+ else:
+ if type(v) == str:
+ v = '"%s"' % v
+ f.writelines(tab*level + '%s = %s;\n' % (k, v))
+
+ for k, v in dict_dict.iteritems():
+ f.writelines('\n')
+ f.writelines(tab*level + '%s :\n' % k)
+ f.writelines(tab*level+'{\n')
+ self.write_dict(f, v, level+1)
+ f.writelines(tab*level+'};\n')
+
+if __name__ == '__main__':
+ cfg = libconfig()
+ cfg.read('test_config.cfg')
+ print cfg
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/utilities/setup.py
--- a/yt/utilities/setup.py
+++ b/yt/utilities/setup.py
@@ -163,6 +163,7 @@
"yt/utilities/data_point_utilities.c",
libraries=["m"])
config.add_subpackage("tests")
+ config.add_subpackage("pyparselibconfig")
config.make_config_py() # installs __config__.py
# config.make_svn_version_py()
return config
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/visualization/_colormap_data.py
--- a/yt/visualization/_colormap_data.py
+++ b/yt/visualization/_colormap_data.py
@@ -1,5 +1,6 @@
### Auto-generated colormap tables, taken from Matplotlib ###
+import numpy as np
from numpy import array
color_map_luts = {}
@@ -7757,6 +7758,64 @@
1.0, 1.0, 1.0, 1.0, 1.0, 1.0]),
)
+color_map_luts["doom"] = (
+array([
+ 0, 31, 23, 75, 255, 27, 19, 11, 7, 47, 35, 23, 15, 79, 71, 63,
+ 255, 247, 243, 235, 231, 223, 219, 211, 203, 199, 191, 187, 179, 175, 167, 163,
+ 155, 151, 143, 139, 131, 127, 119, 115, 107, 103, 95, 91, 83, 79, 71, 67,
+ 255, 255, 255, 255, 255, 255, 255, 255, 255, 247, 239, 231, 223, 215, 207, 203,
+ 191, 179, 171, 163, 155, 143, 135, 127, 119, 107, 95, 83, 75, 63, 51, 43,
+ 239, 231, 223, 219, 211, 203, 199, 191, 183, 179, 171, 167, 159, 151, 147, 139,
+ 131, 127, 119, 111, 107, 99, 91, 87, 79, 71, 67, 59, 55, 47, 39, 35,
+ 119, 111, 103, 95, 91, 83, 75, 67, 63, 55, 47, 39, 31, 23, 19, 11,
+ 191, 183, 175, 167, 159, 155, 147, 139, 131, 123, 119, 111, 103, 95, 87, 83,
+ 159, 143, 131, 119, 103, 91, 79, 67, 123, 111, 103, 91, 83, 71, 63, 55,
+ 255, 235, 215, 195, 175, 155, 135, 115, 255, 255, 255, 255, 255, 255, 255, 255,
+ 255, 239, 227, 215, 203, 191, 179, 167, 155, 139, 127, 115, 103, 91, 79, 67,
+ 231, 199, 171, 143, 115, 83, 55, 27, 0, 0, 0, 0, 0, 0, 0, 0,
+ 255, 255, 255, 255, 255, 255, 255, 255, 243, 235, 223, 215, 203, 195, 183, 175,
+ 255, 255, 255, 255, 255, 255, 255, 255, 167, 159, 147, 135, 79, 67, 55, 47,
+ 0, 0, 0, 0, 0, 0, 0, 0, 255, 255, 255, 255, 207, 159, 111,
+ 167]) / 255.0,
+array([
+ 0, 23, 15, 75, 255, 27, 19, 11, 7, 55, 43, 31, 23, 59, 51, 43,
+ 183, 171, 163, 151, 143, 135, 123, 115, 107, 99, 91, 87, 79, 71, 63, 59,
+ 51, 47, 43, 35, 31, 27, 23, 19, 15, 11, 7, 7, 7, 0, 0, 0,
+ 235, 227, 219, 211, 207, 199, 191, 187, 179, 171, 163, 155, 147, 139, 131, 127,
+ 123, 115, 111, 107, 99, 95, 87, 83, 79, 71, 67, 63, 55, 47, 43, 35,
+ 239, 231, 223, 219, 211, 203, 199, 191, 183, 179, 171, 167, 159, 151, 147, 139,
+ 131, 127, 119, 111, 107, 99, 91, 87, 79, 71, 67, 59, 55, 47, 39, 35,
+ 255, 239, 223, 207, 191, 175, 159, 147, 131, 115, 99, 83, 67, 51, 35, 23,
+ 167, 159, 151, 143, 135, 127, 123, 115, 107, 99, 95, 87, 83, 75, 67, 63,
+ 131, 119, 107, 95, 83, 71, 59, 51, 127, 115, 107, 99, 87, 79, 71, 63,
+ 255, 219, 187, 155, 123, 91, 67, 43, 255, 219, 187, 155, 123, 95, 63, 31,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 231, 199, 171, 143, 115, 83, 55, 27, 0, 0, 0, 0, 0, 0, 0, 0,
+ 255, 235, 215, 199, 179, 163, 143, 127, 115, 111, 103, 95, 87, 79, 71, 67,
+ 255, 255, 255, 255, 255, 255, 255, 255, 63, 55, 47, 35, 59, 47, 35, 27,
+ 0, 0, 0, 0, 0, 0, 0, 0, 159, 231, 123, 0, 0, 0, 0,
+ 107]) / 255.0,
+array([
+ 0, 11, 7, 75, 255, 27, 19, 11, 7, 31, 15, 7, 0, 43, 35, 27,
+ 183, 171, 163, 151, 143, 135, 123, 115, 107, 99, 91, 87, 79, 71, 63, 59,
+ 51, 47, 43, 35, 31, 27, 23, 19, 15, 11, 7, 7, 7, 0, 0, 0,
+ 223, 211, 199, 187, 179, 167, 155, 147, 131, 123, 115, 107, 99, 91, 83, 79,
+ 75, 71, 67, 63, 59, 55, 51, 47, 43, 39, 35, 31, 27, 23, 19, 15,
+ 239, 231, 223, 219, 211, 203, 199, 191, 183, 179, 171, 167, 159, 151, 147, 139,
+ 131, 127, 119, 111, 107, 99, 91, 87, 79, 71, 67, 59, 55, 47, 39, 35,
+ 111, 103, 95, 87, 79, 71, 63, 55, 47, 43, 35, 27, 23, 15, 11, 7,
+ 143, 135, 127, 119, 111, 107, 99, 91, 87, 79, 75, 67, 63, 55, 51, 47,
+ 99, 83, 75, 63, 51, 43, 35, 27, 99, 87, 79, 71, 59, 51, 43, 39,
+ 115, 87, 67, 47, 31, 19, 7, 0, 255, 219, 187, 155, 123, 95, 63, 31,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 255, 255, 255, 255, 255, 255, 255, 255, 255, 227, 203, 179, 155, 131, 107, 83,
+ 255, 219, 187, 155, 123, 91, 59, 27, 23, 15, 15, 11, 7, 0, 0, 0,
+ 255, 215, 179, 143, 107, 71, 35, 0, 0, 0, 0, 0, 39, 27, 19, 11,
+ 83, 71, 59, 47, 35, 23, 11, 0, 67, 75, 255, 255, 207, 155, 107,
+ 107]) / 255.0,
+np.ones(256),
+)
+
color_map_luts['B-W LINEAR'] = color_map_luts['idl00']
color_map_luts['BLUE'] = color_map_luts['idl01']
color_map_luts['GRN-RED-BLU-WHT'] = color_map_luts['idl02']
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/visualization/color_maps.py
--- a/yt/visualization/color_maps.py
+++ b/yt/visualization/color_maps.py
@@ -145,12 +145,12 @@
add_cmap("cubehelix", _cubehelix_data)
# Add colormaps in _colormap_data.py that weren't defined here
-_vs = np.linspace(0,1,255)
+_vs = np.linspace(0,1,256)
for k,v in _cm.color_map_luts.iteritems():
if k not in yt_colormaps and k not in mcm.cmap_d:
- cdict = { 'red': izip(_vs,v[0],v[0]),
- 'green': izip(_vs,v[1],v[1]),
- 'blue': izip(_vs,v[2],v[2]) }
+ cdict = { 'red': zip(_vs,v[0],v[0]),
+ 'green': zip(_vs,v[1],v[1]),
+ 'blue': zip(_vs,v[2],v[2]) }
add_cmap(k, cdict)
def _extract_lookup_table(cmap_name):
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/visualization/plot_container.py
--- a/yt/visualization/plot_container.py
+++ b/yt/visualization/plot_container.py
@@ -30,10 +30,10 @@
from yt.funcs import \
defaultdict, get_image_suffix, \
- get_ipython_api_version, iterable
+ get_ipython_api_version, iterable, \
+ ensure_list
from yt.utilities.exceptions import \
YTNotInsideNotebook
-from ._mpl_imports import FigureCanvasAgg
def invalidate_data(f):
@wraps(f)
@@ -241,7 +241,7 @@
if field is 'all':
fields = self.plots.keys()
else:
- fields = [field]
+ fields = ensure_list(field)
for field in self.data_source._determine_fields(fields):
myzmin = zmin
myzmax = zmax
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -761,7 +761,7 @@
image = self.frb[f]
- if image.max() == image.min():
+ if image.max() == image.min() and zlim == (None, None):
if self._field_transform[f] == log_transform:
mylog.warning("Plot image for field %s has zero dynamic "
"range. Min = Max = %d." % (f, image.max()))
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/visualization/profile_plotter.py
--- a/yt/visualization/profile_plotter.py
+++ b/yt/visualization/profile_plotter.py
@@ -24,7 +24,6 @@
import matplotlib
import numpy as np
import cStringIO
-import __builtin__
from .base_plot_types import ImagePlotMPL
@@ -69,7 +68,7 @@
class FigureContainer(dict):
def __init__(self):
- super(dict, self).__init__()
+ super(FigureContainer, self).__init__()
def __missing__(self, key):
figure = mpl.matplotlib.figure.Figure((10, 8))
@@ -79,13 +78,18 @@
class AxesContainer(dict):
def __init__(self, fig_container):
self.fig_container = fig_container
- super(dict, self).__init__()
+ self.ylim = {}
+ super(AxesContainer, self).__init__()
def __missing__(self, key):
figure = self.fig_container[key]
self[key] = figure.add_subplot(111)
return self[key]
+ def __setitem__(self, key, value):
+ super(AxesContainer, self).__setitem__(key, value)
+ self.ylim[key] = (None, None)
+
def sanitize_label(label, nprofiles):
label = ensure_list(label)
@@ -156,16 +160,17 @@
This creates profiles of a single dataset.
- >>> pf = load("enzo_tiny_cosmology/DD0046/DD0046")
- >>> ad = pf.h.all_data()
+ >>> import yt
+ >>> ds = yt.load("enzo_tiny_cosmology/DD0046/DD0046")
+ >>> ad = ds.all_data()
>>> plot = ProfilePlot(ad, "density", ["temperature", "velocity_x"],
- weight_field="cell_mass",
- plot_spec=dict(color='red', linestyle="--"))
+ ... weight_field="cell_mass",
+ ... plot_spec=dict(color='red', linestyle="--"))
>>> plot.save()
This creates profiles from a time series object.
-
- >>> es = simulation("AMRCosmology.enzo", "Enzo")
+
+ >>> es = yt.simulation("AMRCosmology.enzo", "Enzo")
>>> es.get_time_series()
>>> profiles = []
@@ -218,7 +223,9 @@
self.plot_spec = [dict() for p in self.profiles]
if not isinstance(self.plot_spec, list):
self.plot_spec = [self.plot_spec.copy() for p in self.profiles]
-
+
+ self.figures = FigureContainer()
+ self.axes = AxesContainer(self.figures)
self._setup_plots()
def save(self, name=None):
@@ -235,7 +242,7 @@
self._setup_plots()
unique = set(self.figures.values())
if len(unique) < len(self.figures):
- figiter = izip(xrange(len(unique)), sorted(unique))
+ iters = izip(xrange(len(unique)), sorted(unique))
else:
iters = self.figures.iteritems()
if name is None:
@@ -256,8 +263,7 @@
if isinstance(uid, types.TupleType):
uid = uid[1]
canvas = canvas_cls(fig)
- fn = "%s_1d-Profile_%s_%s%s" % \
- (prefix, xfn, uid, suffix)
+ fn = "%s_1d-Profile_%s_%s%s" % (prefix, xfn, uid, suffix)
mylog.info("Saving %s", fn)
canvas.print_figure(fn)
return self
@@ -276,8 +282,10 @@
Examples
--------
- >>> slc = SlicePlot(pf, "x", ["Density", "VelocityMagnitude"])
- >>> slc.show()
+ >>> import yt
+ >>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+ >>> pp = ProfilePlot(ds.all_data(), 'density', 'temperature')
+ >>> pp.show()
"""
if "__IPYTHON__" in dir(__builtin__):
@@ -290,14 +298,13 @@
else:
raise YTNotInsideNotebook
-
def _repr_html_(self):
"""Return an html representation of the plot object. Will display as a
png for each WindowPlotMPL instance in self.plots"""
ret = ''
unique = set(self.figures.values())
if len(unique) < len(self.figures):
- figiter = izip(xrange(len(unique)), sorted(unique))
+ iters = izip(xrange(len(unique)), sorted(unique))
else:
iters = self.figures.iteritems()
for uid, fig in iters:
@@ -310,14 +317,12 @@
return ret
def _setup_plots(self):
- self.figures = FigureContainer()
- self.axes = AxesContainer(self.figures)
for i, profile in enumerate(self.profiles):
for field, field_data in profile.items():
- self.axes[field].plot(np.array(profile.x),
- np.array(field_data),
- label=self.label[i],
- **self.plot_spec[i])
+ if field in self.axes:
+ self.axes[field].cla()
+ self.axes[field].plot(np.array(profile.x), np.array(field_data),
+ label=self.label[i], **self.plot_spec[i])
# This relies on 'profile' leaking
for fname, axes in self.axes.items():
@@ -327,6 +332,7 @@
axes.set_yscale(yscale)
axes.set_xlabel(xtitle)
axes.set_ylabel(ytitle)
+ axes.set_ylim(*self.axes.ylim[fname])
if any(self.label):
axes.legend(loc="best")
self._plot_valid = True
@@ -443,7 +449,7 @@
def set_unit(self, field, unit):
"""Sets a new unit for the requested field
- parameters
+ Parameters
----------
field : string
The name of the field that is to be changed.
@@ -460,6 +466,97 @@
raise KeyError("Field %s not in profile plot!" % (field))
return self
+ @invalidate_plot
+ def set_xlim(self, xmin=None, xmax=None):
+ """Sets the limits of the bin field
+
+ Parameters
+ ----------
+
+ xmin : float or None
+ The new x minimum. Defaults to None, which leaves the xmin
+ unchanged.
+
+ xmax : float or None
+ The new x maximum. Defaults to None, which leaves the xmax
+ unchanged.
+
+ Examples
+ --------
+
+ >>> import yt
+ >>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+ >>> pp = yt.ProfilePlot(ds.all_data(), 'density', 'temperature')
+ >>> pp.set_xlim(1e-29, 1e-24)
+ >>> pp.save()
+
+ """
+ for i, p in enumerate(self.profiles):
+ if xmin is None:
+ xmi = p.x_bins.min()
+ else:
+ xmi = xmin
+ if xmax is None:
+ xma = p.x_bins.max()
+ else:
+ xma = xmax
+ extrema = {p.x_field: ((xmi, str(p.x.units)), (xma, str(p.x.units)))}
+ units = {p.x_field: str(p.x.units)}
+ for field in p.field_map.values():
+ units[field] = str(p.field_data[field].units)
+ self.profiles[i] = \
+ create_profile(p.data_source, p.x_field,
+ n_bins=len(p.x_bins)-2,
+ fields=p.field_map.values(),
+ weight_field=p.weight_field,
+ accumulation=p.accumulation,
+ fractional=p.fractional,
+ extrema=extrema, units=units)
+ return self
+
+ @invalidate_plot
+ def set_ylim(self, field, ymin=None, ymax=None):
+ """Sets the plot limits for the specified field we are binning.
+
+ Parameters
+ ----------
+
+ field : string or field tuple
+
+ The field that we want to adjust the plot limits for.
+
+ ymin : float or None
+ The new y minimum. Defaults to None, which leaves the ymin
+ unchanged.
+
+ ymax : float or None
+ The new y maximum. Defaults to None, which leaves the ymax
+ unchanged.
+
+ Examples
+ --------
+
+ >>> import yt
+ >>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+ >>> pp = yt.ProfilePlot(ds.all_data(), 'density', ['temperature', 'x-velocity'])
+ >>> pp.set_ylim('temperature', 1e4, 1e6)
+ >>> pp.save()
+
+ """
+ for i, p in enumerate(self.profiles):
+ if field is 'all':
+ fields = self.axes.keys()
+ else:
+ fields = ensure_list(field)
+ for profile in self.profiles:
+ for field in profile.data_source._determine_fields(fields):
+ if field in profile.field_map:
+ field = profile.field_map[field]
+ self.axes.ylim[field] = (ymin, ymax)
+ # Continue on to the next profile.
+ break
+ return self
+
def _get_field_log(self, field_y, profile):
pf = profile.data_source.pf
yf, = profile.data_source._determine_fields([field_y])
@@ -508,7 +605,6 @@
self._get_field_label(field_y, yfi, y_unit, fractional)
return (x_title, y_title)
-
class PhasePlot(ImagePlotContainer):
r"""
@@ -570,10 +666,11 @@
Examples
--------
- >>> pf = load("enzo_tiny_cosmology/DD0046/DD0046")
+ >>> import yt
+ >>> pf = yt.load("enzo_tiny_cosmology/DD0046/DD0046")
>>> ad = pf.h.all_data()
>>> plot = PhasePlot(ad, "density", "temperature", ["cell_mass"],
- weight_field=None)
+ ... weight_field=None)
>>> plot.save()
>>> # Change plot properties.
@@ -706,7 +803,7 @@
self.plots[f] = PhasePlotMPL(self.profile.x, self.profile.y, data,
x_scale, y_scale, z_scale,
- self._colormaps[f], np.array(zlim),
+ self._colormaps[f], zlim,
self.figure_size, fp.get_size(),
fig, axes, cax)
@@ -750,6 +847,8 @@
self._setup_plots()
if mpl_kwargs is None:
mpl_kwargs = {}
+ if name is None:
+ name = str(self.profile.pf)
xfn = self.profile.x_field
yfn = self.profile.y_field
if isinstance(xfn, types.TupleType):
@@ -761,11 +860,10 @@
if isinstance(f, types.TupleType):
_f = _f[1]
middle = "2d-Profile_%s_%s_%s" % (xfn, yfn, _f)
- if name[-1] == os.sep and not os.path.isdir(name):
- os.mkdir(name)
- if name is None:
- prefix = self.profile.pf
- elif os.path.isdir(name) and name != str(self.profile.pf):
+ splitname = os.path.split(name)
+ if splitname[0] != '' and not os.path.isdir(splitname[0]):
+ os.makedirs(splitname[0])
+ if os.path.isdir(name) and name != str(self.profile.pf):
prefix = name + (os.sep if name[-1] != os.sep else '')
prefix += str(self.profile.pf)
else:
@@ -775,7 +873,6 @@
for k, v in self.plots.iteritems():
names.append(v.save(name, mpl_kwargs))
return names
-
fn = "%s_%s%s" % (prefix, middle, suffix)
names.append(fn)
self.plots[f].save(fn, mpl_kwargs)
@@ -837,7 +934,7 @@
def set_unit(self, field, unit):
"""Sets a new unit for the requested field
- parameters
+ Parameters
----------
field : string
The name of the field that is to be changed.
@@ -852,10 +949,111 @@
self.profile.set_y_unit(unit)
elif field in fields:
self.profile.set_field_unit(field, unit)
+ self.plots[field].zmin, self.plots[field].zmax = (None, None)
else:
raise KeyError("Field %s not in phase plot!" % (field))
return self
+ @invalidate_plot
+ def set_xlim(self, xmin=None, xmax=None):
+ """Sets the limits of the x bin field
+
+ Parameters
+ ----------
+
+ xmin : float or None
+ The new x minimum. Defaults to None, which leaves the xmin
+ unchanged.
+
+ xmax : float or None
+ The new x maximum. Defaults to None, which leaves the xmax
+ unchanged.
+
+ Examples
+ --------
+
+ >>> import yt
+ >>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+ >>> pp = yt.PhasePlot(ds.all_data(), 'density', 'temperature', 'cell_mass')
+ >>> pp.set_xlim(1e-29, 1e-24)
+ >>> pp.save()
+
+ """
+ p = self.profile
+ if xmin is None:
+ xmin = p.x_bins.min()
+ if xmax is None:
+ xmax = p.x_bins.max()
+ units = {p.x_field: str(p.x.units),
+ p.y_field: str(p.y.units)}
+ zunits = dict((field, str(p.field_units[field])) for field in p.field_units)
+ extrema = {p.x_field: ((xmin, str(p.x.units)), (xmax, str(p.x.units))),
+ p.y_field: ((p.y_bins.min(), str(p.y.units)),
+ (p.y_bins.max(), str(p.y.units)))}
+ self.profile = create_profile(
+ p.data_source,
+ [p.x_field, p.y_field],
+ p.field_map.values(),
+ n_bins=[len(p.x_bins)-2, len(p.y_bins)-2],
+ weight_field=p.weight_field,
+ accumulation=p.accumulation,
+ fractional=p.fractional,
+ units=units,
+ extrema=extrema)
+ for field in zunits:
+ self.profile.set_field_unit(field, zunits[field])
+ return self
+
+ @invalidate_plot
+ def set_ylim(self, ymin=None, ymax=None):
+ """Sets the plot limits for the y bin field.
+
+ Parameters
+ ----------
+
+ ymin : float or None
+ The new y minimum. Defaults to None, which leaves the ymin
+ unchanged.
+
+ ymax : float or None
+ The new y maximum. Defaults to None, which leaves the ymax
+ unchanged.
+
+ Examples
+ --------
+
+ >>> import yt
+ >>> ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
+ >>> pp = yt.PhasePlot(ds.all_data(), 'density', 'temperature', 'cell_mass')
+ >>> pp.set_ylim(1e4, 1e6)
+ >>> pp.save()
+
+ """
+ p = self.profile
+ if ymin is None:
+ ymin = p.y_bins.min()
+ if ymax is None:
+ ymax = p.y_bins.max()
+ units = {p.x_field: str(p.x.units),
+ p.y_field: str(p.y.units)}
+ zunits = dict((field, str(p.field_units[field])) for field in p.field_units)
+ extrema = {p.x_field: ((p.x_bins.min(), str(p.x.units)),
+ (p.x_bins.max(), str(p.x.units))),
+ p.y_field: ((ymin, str(p.y.units)), (ymax, str(p.y.units)))}
+ self.profile = create_profile(
+ p.data_source,
+ [p.x_field, p.y_field],
+ p.field_map.values(),
+ n_bins=[len(p.x_bins), len(p.y_bins)],
+ weight_field=p.weight_field,
+ accumulation=p.accumulation,
+ fractional=p.fractional,
+ units=units,
+ extrema=extrema)
+ for field in zunits:
+ self.profile.set_field_unit(field, zunits[field])
+ return self
+
def run_callbacks(self, *args):
raise NotImplementedError
def setup_callbacks(self, *args):
@@ -900,10 +1098,10 @@
norm = matplotlib.colors.Normalize(zlim[0], zlim[1])
self.image = None
self.cb = None
- self.image = self.axes.pcolormesh(np.array(x_data),
- np.array(y_data),
+ self.image = self.axes.pcolormesh(np.array(x_data),
+ np.array(y_data),
np.array(image_data.T),
- norm=norm,
+ norm=norm,
cmap=cmap)
self.axes.set_xscale(x_scale)
self.axes.set_yscale(y_scale)
diff -r a76ff96ec74ff9bd30919f8ef2b50707a7ce16c2 -r 87039982afe8b3bb4c43d177471ecbbb6b527416 yt/visualization/volume_rendering/transfer_functions.py
--- a/yt/visualization/volume_rendering/transfer_functions.py
+++ b/yt/visualization/volume_rendering/transfer_functions.py
@@ -466,10 +466,10 @@
ax.fill_between(np.arange(self.alpha.y.size), self.alpha.x.size * self.alpha.y, y2=self.alpha.x.size, color='white')
ax.set_xlim(0, self.alpha.x.size)
xticks = np.arange(np.ceil(self.alpha.x[0]), np.floor(self.alpha.x[-1]) + 1, 1) - self.alpha.x[0]
- xticks *= self.alpha.x.size / (self.alpha.x[-1] - self.alpha.x[0])
+ xticks *= (self.alpha.x.size-1) / (self.alpha.x[-1] - self.alpha.x[0])
ax.xaxis.set_ticks(xticks)
def x_format(x, pos):
- return "%.1f" % (x * (self.alpha.x[-1] - self.alpha.x[0]) / (self.alpha.x.size) + self.alpha.x[0])
+ return "%.1f" % (x * (self.alpha.x[-1] - self.alpha.x[0]) / (self.alpha.x.size-1) + self.alpha.x[0])
ax.xaxis.set_major_formatter(FuncFormatter(x_format))
yticks = np.linspace(0,1,5) * self.alpha.y.size
ax.yaxis.set_ticks(yticks)
@@ -507,12 +507,12 @@
ax.fill_between(np.arange(self.alpha.y.size), self.alpha.x.size * self.alpha.y, y2=self.alpha.x.size, color='white')
ax.set_xlim(0, self.alpha.x.size)
xticks = np.arange(np.ceil(self.alpha.x[0]), np.floor(self.alpha.x[-1]) + 1, 1) - self.alpha.x[0]
- xticks *= self.alpha.x.size / (self.alpha.x[-1] - self.alpha.x[0])
+ xticks *= (self.alpha.x.size-1) / (self.alpha.x[-1] - self.alpha.x[0])
if len(xticks) > 5:
xticks = xticks[::len(xticks)/5]
ax.xaxis.set_ticks(xticks)
def x_format(x, pos):
- return "%.1f" % (x * (self.alpha.x[-1] - self.alpha.x[0]) / (self.alpha.x.size) + self.alpha.x[0])
+ return "%.1f" % (x * (self.alpha.x[-1] - self.alpha.x[0]) / (self.alpha.x.size-1) + self.alpha.x[0])
ax.xaxis.set_major_formatter(FuncFormatter(x_format))
yticks = np.linspace(0,1,5) * self.alpha.y.size
ax.yaxis.set_ticks(yticks)
@@ -558,18 +558,26 @@
# Set TF limits based on what is visible
visible = np.argwhere(self.alpha.y > 1.0e-3*self.alpha.y.max())
+
# Display colobar values
xticks = np.arange(np.ceil(self.alpha.x[0]), np.floor(self.alpha.x[-1]) + 1, 1) - self.alpha.x[0]
- xticks *= self.alpha.x.size / (self.alpha.x[-1] - self.alpha.x[0])
+ xticks *= (self.alpha.x.size-1) / (self.alpha.x[-1] - self.alpha.x[0])
if len(xticks) > 5:
xticks = xticks[::len(xticks)/5]
# Add colorbar limits to the ticks (May not give ideal results)
xticks = np.append(visible[0], xticks)
xticks = np.append(visible[-1], xticks)
+ # remove dupes
+ xticks = list(set(xticks))
ax.yaxis.set_ticks(xticks)
def x_format(x, pos):
- return "%.1f" % (x * (self.alpha.x[-1] - self.alpha.x[0]) / (self.alpha.x.size) + self.alpha.x[0])
+ val = x * (self.alpha.x[-1] - self.alpha.x[0]) / (self.alpha.x.size-1) + self.alpha.x[0]
+ if abs(val) < 1.e-3 or abs(val) > 1.e4:
+ e = np.floor(np.log10(abs(val)))
+ return r"${:.2f}\times 10^{:d}$".format(val/10.0**e, int(e))
+ else:
+ return "%.1g" % (val)
ax.yaxis.set_major_formatter(FuncFormatter(x_format))
yticks = np.linspace(0,1,2,endpoint=True) * max_alpha
https://bitbucket.org/yt_analysis/yt/commits/84ad27924c60/
Changeset: 84ad27924c60
Branch: yt-3.0
User: jzuhone
Date: 2014-06-12 06:25:27
Summary: Merge
Affected #: 0 files
https://bitbucket.org/yt_analysis/yt/commits/ad397535f0ca/
Changeset: ad397535f0ca
Branch: yt-3.0
User: jzuhone
Date: 2014-06-12 07:52:07
Summary: Fixing Windows test failures
Affected #: 2 files
diff -r 84ad27924c602deb9649d55de29b0c63586f8ec1 -r ad397535f0caf0af2a89ab010fe7c21c6538c189 yt/units/yt_array.py
--- a/yt/units/yt_array.py
+++ b/yt/units/yt_array.py
@@ -563,7 +563,7 @@
unit_lut = pickle.loads(dataset.attrs['unit_registry'])
else:
unit_lut = None
-
+ f.close()
registry = UnitRegistry(lut=unit_lut, add_default_symbols=False)
return cls(data, units, registry=registry)
diff -r 84ad27924c602deb9649d55de29b0c63586f8ec1 -r ad397535f0caf0af2a89ab010fe7c21c6538c189 yt/utilities/lib/tests/test_ragged_arrays.py
--- a/yt/utilities/lib/tests/test_ragged_arrays.py
+++ b/yt/utilities/lib/tests/test_ragged_arrays.py
@@ -13,7 +13,7 @@
def test_index_unop():
np.random.seed(0x4d3d3d3)
- indices = np.arange(1000)
+ indices = np.arange(1000, dtype="int64")
np.random.shuffle(indices)
sizes = np.array([
200, 50, 50, 100, 32, 32, 32, 32, 32, 64, 376], dtype="int64")
https://bitbucket.org/yt_analysis/yt/commits/b6f0110bc38e/
Changeset: b6f0110bc38e
Branch: yt-3.0
User: jzuhone
Date: 2014-06-18 15:25:24
Summary: This should resolve the FITS answer testing issues.
Affected #: 1 file
diff -r ad397535f0caf0af2a89ab010fe7c21c6538c189 -r b6f0110bc38ec3b9400146d0609e71383f1ee922 yt/frontends/fits/tests/test_outputs.py
--- a/yt/frontends/fits/tests/test_outputs.py
+++ b/yt/frontends/fits/tests/test_outputs.py
@@ -19,13 +19,14 @@
small_patch_amr, \
big_patch_amr, \
data_dir_load
+from ..data_structures import FITSDataset
_fields = ("intensity")
m33 = "radio_fits/m33_hi.fits"
@requires_pf(m33, big_data=True)
def test_m33():
- pf = data_dir_load(m33, nan_mask=0.0)
+ pf = data_dir_load(m33, cls=FITSDataset, nan_mask=0.0)
yield assert_equal, str(pf), "m33_hi.fits"
for test in small_patch_amr(m33, _fields):
test_m33.__name__ = test.description
@@ -36,7 +37,7 @@
grs = "radio_fits/grs-50-cube.fits"
@requires_pf(grs)
def test_grs():
- pf = data_dir_load(grs, nan_mask=0.0)
+ pf = data_dir_load(grs, cls=FITSDataset, nan_mask=0.0)
yield assert_equal, str(pf), "grs-50-cube.fits"
for test in small_patch_amr(grs, _fields):
test_grs.__name__ = test.description
@@ -47,7 +48,7 @@
vf = "UniformGrid/velocity_field_20.fits"
@requires_pf(vf)
def test_velocity_field():
- pf = data_dir_load(bf)
+ pf = data_dir_load(bf, cls=FITSDataset)
yield assert_equal, str(pf), "velocity_field_20.fits"
for test in small_patch_amr(vf, _fields):
test_velocity_field.__name__ = test.description
https://bitbucket.org/yt_analysis/yt/commits/8f8bb4c481d1/
Changeset: 8f8bb4c481d1
Branch: yt-3.0
User: jzuhone
Date: 2014-06-18 15:25:46
Summary: Merge
Affected #: 31 files
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 doc/README
--- a/doc/README
+++ b/doc/README
@@ -5,6 +5,6 @@
http://sphinx.pocoo.org/
Because the documentation requires a number of dependencies, we provide
-pre-build versions online, accessible here:
+pre-built versions online, accessible here:
-http://yt-project.org/docs/
+http://yt-project.org/docs/dev-3.0/
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 doc/source/cookbook/amrkdtree_to_uniformgrid.py
--- /dev/null
+++ b/doc/source/cookbook/amrkdtree_to_uniformgrid.py
@@ -0,0 +1,33 @@
+import numpy as np
+import yt
+
+#This is an example of how to map an amr data set
+#to a uniform grid. In this case the highest
+#level of refinement is mapped into a 1024x1024x1024 cube
+
+#first the amr data is loaded
+ds = yt.load("~/pfs/galaxy/new_tests/feedback_8bz/DD0021/DD0021")
+
+#next we get the maxium refinement level
+lmax = ds.parameters['MaximumRefinementLevel']
+
+#calculate the center of the domain
+domain_center = (ds.domain_right_edge - ds.domain_left_edge)/2
+
+#determine the cellsize in the highest refinement level
+cell_size = pf.domain_width/(pf.domain_dimensions*2**lmax)
+
+#calculate the left edge of the new grid
+left_edge = domain_center - 512*cell_size
+
+#the number of cells per side of the new grid
+ncells = 1024
+
+#ask yt for the specified covering grid
+cgrid = pf.h.covering_grid(lmax, left_edge, np.array([ncells,]*3))
+
+#get a map of the density into the new grid
+density_map = cgrid["density"].astype(dtype="float32")
+
+#save the file as a numpy array for convenient future processing
+np.save("/pfs/goldbaum/galaxy/new_tests/feedback_8bz/gas_density_DD0021_log_densities.npy", density_map)
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 doc/source/cookbook/ffmpeg_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/ffmpeg_volume_rendering.py
@@ -0,0 +1,99 @@
+#This is an example of how to make videos of
+#uniform grid data using Theia and ffmpeg
+
+#The Scene object to hold the ray caster and view camera
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+
+#GPU based raycasting algorithm to use
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+
+#These will be used to define how to color the data
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+#This will be used to launch ffmpeg
+import subprocess as sp
+
+#Of course we need numpy for math magic
+import numpy as np
+
+#Opacity scaling function
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, (v-mi)/(ma-mi) + 0.0)
+
+#load the uniform grid from a numpy array file
+bolshoi = "/home/bogert/log_densities_1024.npy"
+density_grid = np.load(bolshoi)
+
+#Set the TheiaScene to use the density_grid and
+#setup the raycaster for a resulting 1080p image
+ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (1920,1080) ))
+
+#the min and max values in the data to color
+mi, ma = 0.0, 3.6
+
+#setup colortransferfunction
+bins = 5000
+tf = ColorTransferFunction( (mi, ma), bins)
+tf.map_to_colormap(0.5, ma, colormap="spring", scale_func = scale_func)
+
+#pass the transfer function to the ray caster
+ts.source.raycaster.set_transfer(tf)
+
+#Initial configuration for start of video
+#set initial opacity and brightness values
+#then zoom into the center of the data 30%
+ts.source.raycaster.set_opacity(0.03)
+ts.source.raycaster.set_brightness(2.3)
+ts.camera.zoom(30.0)
+
+#path to ffmpeg executable
+FFMPEG_BIN = "/usr/local/bin/ffmpeg"
+
+pipe = sp.Popen([ FFMPEG_BIN,
+ '-y', # (optional) overwrite the output file if it already exists
+ #This must be set to rawvideo because the image is an array
+ '-f', 'rawvideo',
+ #This must be set to rawvideo because the image is an array
+ '-vcodec','rawvideo',
+ #The size of the image array and resulting video
+ '-s', '1920x1080',
+ #This must be rgba to match array format (uint32)
+ '-pix_fmt', 'rgba',
+ #frame rate of video
+ '-r', '29.97',
+ #Indicate that the input to ffmpeg comes from a pipe
+ '-i', '-',
+ # Tells FFMPEG not to expect any audio
+ '-an',
+ #Setup video encoder
+ #Use any encoder you life available from ffmpeg
+ '-vcodec', 'libx264', '-preset', 'ultrafast', '-qp', '0',
+ '-pix_fmt', 'yuv420p',
+ #Name of the output
+ 'bolshoiplanck2.mkv' ],
+ stdin=sp.PIPE,stdout=sp.PIPE)
+
+
+#Now we loop and produce 500 frames
+for k in range (0,500) :
+ #update the scene resulting in a new image
+ ts.update()
+
+ #get the image array from the ray caster
+ array = ts.source.get_results()
+
+ #send the image array to ffmpeg
+ array.tofile(pipe.stdin)
+
+ #rotate the scene by 0.01 rads in x,y & z
+ ts.camera.rotateX(0.01)
+ ts.camera.rotateZ(0.01)
+ ts.camera.rotateY(0.01)
+
+ #zoom in 0.01% for a total of a 5% zoom
+ ts.camera.zoom(0.01)
+
+
+#Close the pipe to ffmpeg
+pipe.terminate()
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 doc/source/cookbook/opengl_stereo_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/opengl_stereo_volume_rendering.py
@@ -0,0 +1,370 @@
+from OpenGL.GL import *
+from OpenGL.GLUT import *
+from OpenGL.GLU import *
+from OpenGL.GL.ARB.vertex_buffer_object import *
+
+import sys, time
+import numpy as np
+import pycuda.driver as cuda_driver
+import pycuda.gl as cuda_gl
+
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+import numexpr as ne
+
+window = None # Number of the glut window.
+rot_enabled = True
+
+#Theia Scene
+ts = None
+
+#RAY CASTING values
+c_tbrightness = 1.0
+c_tdensity = 0.05
+
+output_texture = None # pointer to offscreen render target
+
+leftButton = False
+middleButton = False
+rightButton = False
+
+#Screen width and height
+width = 1920
+height = 1080
+
+eyesep = 0.1
+
+(pbo, pycuda_pbo) = [None]*2
+(rpbo, rpycuda_pbo) = [None]*2
+
+#create 2 PBO for stereo scopic rendering
+def create_PBO(w, h):
+ global pbo, pycuda_pbo, rpbo, rpycuda_pbo
+ num_texels = w*h
+ array = np.zeros((num_texels, 3),np.float32)
+
+ pbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, pbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pycuda_pbo = cuda_gl.RegisteredBuffer(long(pbo))
+
+ rpbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, rpbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ rpycuda_pbo = cuda_gl.RegisteredBuffer(long(rpbo))
+
+def destroy_PBO(self):
+ global pbo, pycuda_pbo, rpbo, rpycuda_pbo
+ glBindBuffer(GL_ARRAY_BUFFER, long(pbo))
+ glDeleteBuffers(1, long(pbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pbo,pycuda_pbo = [None]*2
+
+ glBindBuffer(GL_ARRAY_BUFFER, long(rpbo))
+ glDeleteBuffers(1, long(rpbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ rpbo,rpycuda_pbo = [None]*2
+
+#consistent with C initPixelBuffer()
+def create_texture(w,h):
+ global output_texture
+ output_texture = glGenTextures(1)
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+ # set basic parameters
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
+ # buffer data
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ w, h, 0, GL_RGB, GL_FLOAT, None)
+
+#consistent with C initPixelBuffer()
+def destroy_texture():
+ global output_texture
+ glDeleteTextures(output_texture);
+ output_texture = None
+
+def init_gl(w = 512 , h = 512):
+ Width, Height = (w, h)
+
+ glClearColor(0.1, 0.1, 0.5, 1.0)
+ glDisable(GL_DEPTH_TEST)
+
+ #matrix functions
+ glViewport(0, 0, Width, Height)
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+
+ #matrix functions
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+ glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)
+
+def resize(Width, Height):
+ global width, height
+ (width, height) = Width, Height
+ glViewport(0, 0, Width, Height) # Reset The Current Viewport And Perspective Transformation
+ glMatrixMode(GL_PROJECTION)
+ glLoadIdentity()
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+
+
+def do_tick():
+ global time_of_last_titleupdate, frame_counter, frames_per_second
+ if ((time.clock () * 1000.0) - time_of_last_titleupdate >= 1000.):
+ frames_per_second = frame_counter # Save The FPS
+ frame_counter = 0 # Reset The FPS Counter
+ szTitle = "%d FPS" % (frames_per_second )
+ glutSetWindowTitle ( szTitle )
+ time_of_last_titleupdate = time.clock () * 1000.0
+ frame_counter += 1
+
+oldMousePos = [ 0, 0 ]
+def mouseButton( button, mode, x, y ):
+ """Callback function (mouse button pressed or released).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively"""
+
+ global leftButton, middleButton, rightButton, oldMousePos
+
+ if button == GLUT_LEFT_BUTTON:
+ if mode == GLUT_DOWN:
+ leftButton = True
+ else:
+ leftButton = False
+
+ if button == GLUT_MIDDLE_BUTTON:
+ if mode == GLUT_DOWN:
+ middleButton = True
+ else:
+ middleButton = False
+
+ if button == GLUT_RIGHT_BUTTON:
+ if mode == GLUT_DOWN:
+ rightButton = True
+ else:
+ rightButton = False
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def mouseMotion( x, y ):
+ """Callback function (mouse moved while button is pressed).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively.
+ The global translation vector is updated according to
+ the movement of the mouse pointer."""
+
+ global ts, leftButton, middleButton, rightButton, oldMousePos
+ deltaX = x - oldMousePos[ 0 ]
+ deltaY = y - oldMousePos[ 1 ]
+
+ factor = 0.001
+
+ if leftButton == True:
+ ts.camera.rotateX( - deltaY * factor)
+ ts.camera.rotateY( - deltaX * factor)
+ if middleButton == True:
+ ts.camera.translateX( deltaX* 2.0 * factor)
+ ts.camera.translateY( - deltaY* 2.0 * factor)
+ if rightButton == True:
+ ts.camera.scale += deltaY * factor
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def keyPressed(*args):
+ global c_tbrightness, c_tdensity, eyesep
+ # If escape is pressed, kill everything.
+ if args[0] == '\033':
+ print 'Closing..'
+ destroy_PBOs()
+ destroy_texture()
+ exit()
+
+ #change the brightness of the scene
+ elif args[0] == ']':
+ c_tbrightness += 0.025
+ elif args[0] == '[':
+ c_tbrightness -= 0.025
+
+ #change the density scale
+ elif args[0] == ';':
+ c_tdensity -= 0.001
+ elif args[0] == '\'':
+ c_tdensity += 0.001
+
+ #change the transfer scale
+ elif args[0] == '-':
+ eyesep -= 0.01
+ elif args[0] == '=':
+ eyesep += 0.01
+
+def idle():
+ glutPostRedisplay()
+
+def display():
+ try:
+ #process left eye
+ process_image()
+ display_image()
+
+ #process right eye
+ process_image(eye = False)
+ display_image(eye = False)
+
+
+ glutSwapBuffers()
+
+ except:
+ from traceback import print_exc
+ print_exc()
+ from os import _exit
+ _exit(0)
+
+def process(eye = True):
+ global ts, pycuda_pbo, rpycuda_pbo, eyesep, c_tbrightness, c_tdensity
+ """ Use PyCuda """
+
+ ts.get_raycaster().set_opacity(c_tdensity)
+ ts.get_raycaster().set_brightness(c_tbrightness)
+
+ if (eye) :
+ ts.camera.translateX(-eyesep)
+ dest_mapping = pycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ dest_mapping.unmap()
+ ts.camera.translateX(eyesep)
+ else :
+ ts.camera.translateX(eyesep)
+ dest_mapping = rpycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ dest_mapping.unmap()
+ ts.camera.translateX(-eyesep)
+
+
+def process_image(eye = True):
+ global output_texture, pbo, rpbo, width, height
+ """ copy image and process using CUDA """
+ # run the Cuda kernel
+ process(eye)
+ # download texture from PBO
+ if (eye) :
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(pbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ width, height, 0,
+ GL_RGB, GL_FLOAT, None)
+ else :
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(rpbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ width, height, 0,
+ GL_RGB, GL_FLOAT, None)
+
+def display_image(eye = True):
+ global width, height
+ """ render a screen sized quad """
+ glDisable(GL_DEPTH_TEST)
+ glDisable(GL_LIGHTING)
+ glEnable(GL_TEXTURE_2D)
+ glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE)
+
+ #matix functions should be moved
+ glMatrixMode(GL_PROJECTION)
+ glPushMatrix()
+ glLoadIdentity()
+ glOrtho(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)
+ glMatrixMode( GL_MODELVIEW)
+ glLoadIdentity()
+ glViewport(0, 0, width, height)
+
+ if (eye) :
+ glDrawBuffer(GL_BACK_LEFT)
+ else :
+ glDrawBuffer(GL_BACK_RIGHT)
+
+ glBegin(GL_QUADS)
+ glTexCoord2f(0.0, 0.0)
+ glVertex3f(-1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 0.0)
+ glVertex3f(1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 1.0)
+ glVertex3f(1.0, 1.0, 0.5)
+ glTexCoord2f(0.0, 1.0)
+ glVertex3f(-1.0, 1.0, 0.5)
+ glEnd()
+
+ glMatrixMode(GL_PROJECTION)
+ glPopMatrix()
+
+ glDisable(GL_TEXTURE_2D)
+ glBindTexture(GL_TEXTURE_2D, 0)
+ glBindBuffer(GL_PIXEL_PACK_BUFFER, 0)
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0)
+
+
+#note we may need to init cuda_gl here and pass it to camera
+def main():
+ global window, ts, width, height
+ (width, height) = (1920, 1080)
+
+ glutInit(sys.argv)
+ glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH | GLUT_STEREO)
+ glutInitWindowSize(*initial_size)
+ glutInitWindowPosition(0, 0)
+ window = glutCreateWindow("Stereo Volume Rendering")
+
+
+ glutDisplayFunc(display)
+ glutIdleFunc(idle)
+ glutReshapeFunc(resize)
+ glutMouseFunc( mouseButton )
+ glutMotionFunc( mouseMotion )
+ glutKeyboardFunc(keyPressed)
+ init_gl(width, height)
+
+ # create texture for blitting to screen
+ create_texture(width, height)
+
+ import pycuda.gl.autoinit
+ import pycuda.gl
+ cuda_gl = pycuda.gl
+
+ create_PBO(width, height)
+ # ----- Load and Set Volume Data -----
+
+ density_grid = np.load("/home/bogert/dd150_log_densities.npy")
+
+ mi, ma= 21.5, 24.5
+ bins = 5000
+ tf = ColorTransferFunction( (mi, ma), bins)
+ tf.map_to_colormap(mi, ma, colormap="algae", scale_func = scale_func)
+
+ ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (width, height), tf = tf))
+
+ ts.get_raycaster().set_sample_size(0.01)
+ ts.get_raycaster().set_max_samples(5000)
+
+ glutMainLoop()
+
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, np.abs((v)-ma)/np.abs(mi-ma) + 0.0)
+
+# Print message to console, and kick off the main to get it rolling.
+if __name__ == "__main__":
+ print "Hit ESC key to quit, 'a' to toggle animation, and 'e' to toggle cuda"
+ main()
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 doc/source/cookbook/opengl_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/opengl_volume_rendering.py
@@ -0,0 +1,322 @@
+from OpenGL.GL import *
+from OpenGL.GLUT import *
+from OpenGL.GLU import *
+from OpenGL.GL.ARB.vertex_buffer_object import *
+
+import sys, time
+import numpy as np
+import pycuda.driver as cuda_driver
+import pycuda.gl as cuda_gl
+
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+import numexpr as ne
+
+window = None # Number of the glut window.
+rot_enabled = True
+
+#Theia Scene
+ts = None
+
+#RAY CASTING values
+c_tbrightness = 1.0
+c_tdensity = 0.05
+
+output_texture = None # pointer to offscreen render target
+
+leftButton = False
+middleButton = False
+rightButton = False
+
+#Screen width and height
+width = 1024
+height = 1024
+
+eyesep = 0.1
+
+(pbo, pycuda_pbo) = [None]*2
+
+def create_PBO(w, h):
+ global pbo, pycuda_pbo
+ num_texels = w*h
+ array = np.zeros((w,h,3),np.uint32)
+
+ pbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, pbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pycuda_pbo = cuda_gl.RegisteredBuffer(long(pbo))
+
+def destroy_PBO(self):
+ global pbo, pycuda_pbo
+ glBindBuffer(GL_ARRAY_BUFFER, long(pbo))
+ glDeleteBuffers(1, long(pbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pbo,pycuda_pbo = [None]*2
+
+#consistent with C initPixelBuffer()
+def create_texture(w,h):
+ global output_texture
+ output_texture = glGenTextures(1)
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+ # set basic parameters
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
+ # buffer data
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA,
+ w, h, 0, GL_RGBA, GL_UNSIGNED_INT_8_8_8_8, None)
+
+#consistent with C initPixelBuffer()
+def destroy_texture():
+ global output_texture
+ glDeleteTextures(output_texture);
+ output_texture = None
+
+def init_gl(w = 512 , h = 512):
+ Width, Height = (w, h)
+
+ glClearColor(0.1, 0.1, 0.5, 1.0)
+ glDisable(GL_DEPTH_TEST)
+
+ #matrix functions
+ glViewport(0, 0, Width, Height)
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+
+ #matrix functions
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+ glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)
+
+def resize(Width, Height):
+ global width, height
+ (width, height) = Width, Height
+ glViewport(0, 0, Width, Height) # Reset The Current Viewport And Perspective Transformation
+ glMatrixMode(GL_PROJECTION)
+ glLoadIdentity()
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+
+
+def do_tick():
+ global time_of_last_titleupdate, frame_counter, frames_per_second
+ if ((time.clock () * 1000.0) - time_of_last_titleupdate >= 1000.):
+ frames_per_second = frame_counter # Save The FPS
+ frame_counter = 0 # Reset The FPS Counter
+ szTitle = "%d FPS" % (frames_per_second )
+ glutSetWindowTitle ( szTitle )
+ time_of_last_titleupdate = time.clock () * 1000.0
+ frame_counter += 1
+
+oldMousePos = [ 0, 0 ]
+def mouseButton( button, mode, x, y ):
+ """Callback function (mouse button pressed or released).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively"""
+
+ global leftButton, middleButton, rightButton, oldMousePos
+
+ if button == GLUT_LEFT_BUTTON:
+ if mode == GLUT_DOWN:
+ leftButton = True
+ else:
+ leftButton = False
+
+ if button == GLUT_MIDDLE_BUTTON:
+ if mode == GLUT_DOWN:
+ middleButton = True
+ else:
+ middleButton = False
+
+ if button == GLUT_RIGHT_BUTTON:
+ if mode == GLUT_DOWN:
+ rightButton = True
+ else:
+ rightButton = False
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def mouseMotion( x, y ):
+ """Callback function (mouse moved while button is pressed).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively.
+ The global translation vector is updated according to
+ the movement of the mouse pointer."""
+
+ global ts, leftButton, middleButton, rightButton, oldMousePos
+ deltaX = x - oldMousePos[ 0 ]
+ deltaY = y - oldMousePos[ 1 ]
+
+ factor = 0.001
+
+ if leftButton == True:
+ ts.camera.rotateX( - deltaY * factor)
+ ts.camera.rotateY( - deltaX * factor)
+ if middleButton == True:
+ ts.camera.translateX( deltaX* 2.0 * factor)
+ ts.camera.translateY( - deltaY* 2.0 * factor)
+ if rightButton == True:
+ ts.camera.scale += deltaY * factor
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def keyPressed(*args):
+ global c_tbrightness, c_tdensity
+ # If escape is pressed, kill everything.
+ if args[0] == '\033':
+ print 'Closing..'
+ destroy_PBOs()
+ destroy_texture()
+ exit()
+
+ #change the brightness of the scene
+ elif args[0] == ']':
+ c_tbrightness += 0.025
+ elif args[0] == '[':
+ c_tbrightness -= 0.025
+
+ #change the density scale
+ elif args[0] == ';':
+ c_tdensity -= 0.001
+ elif args[0] == '\'':
+ c_tdensity += 0.001
+
+def idle():
+ glutPostRedisplay()
+
+def display():
+ try:
+ #process left eye
+ process_image()
+ display_image()
+
+ glutSwapBuffers()
+
+ except:
+ from traceback import print_exc
+ print_exc()
+ from os import _exit
+ _exit(0)
+
+def process(eye = True):
+ global ts, pycuda_pbo, eyesep, c_tbrightness, c_tdensity
+
+ ts.get_raycaster().set_opacity(c_tdensity)
+ ts.get_raycaster().set_brightness(c_tbrightness)
+
+ dest_mapping = pycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ # ts.get_raycaster().cast()
+ dest_mapping.unmap()
+
+
+def process_image():
+ global output_texture, pbo, width, height
+ """ copy image and process using CUDA """
+ # run the Cuda kernel
+ process()
+ # download texture from PBO
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(pbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA,
+ width, height, 0, GL_RGBA, GL_UNSIGNED_INT_8_8_8_8_REV, None)
+
+def display_image(eye = True):
+ global width, height
+ """ render a screen sized quad """
+ glDisable(GL_DEPTH_TEST)
+ glDisable(GL_LIGHTING)
+ glEnable(GL_TEXTURE_2D)
+ glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE)
+
+ #matix functions should be moved
+ glMatrixMode(GL_PROJECTION)
+ glPushMatrix()
+ glLoadIdentity()
+ glOrtho(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)
+ glMatrixMode( GL_MODELVIEW)
+ glLoadIdentity()
+ glViewport(0, 0, width, height)
+
+ glBegin(GL_QUADS)
+ glTexCoord2f(0.0, 0.0)
+ glVertex3f(-1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 0.0)
+ glVertex3f(1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 1.0)
+ glVertex3f(1.0, 1.0, 0.5)
+ glTexCoord2f(0.0, 1.0)
+ glVertex3f(-1.0, 1.0, 0.5)
+ glEnd()
+
+ glMatrixMode(GL_PROJECTION)
+ glPopMatrix()
+
+ glDisable(GL_TEXTURE_2D)
+ glBindTexture(GL_TEXTURE_2D, 0)
+ glBindBuffer(GL_PIXEL_PACK_BUFFER, 0)
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0)
+
+
+#note we may need to init cuda_gl here and pass it to camera
+def main():
+ global window, ts, width, height
+ (width, height) = (1024, 1024)
+
+ glutInit(sys.argv)
+ glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH )
+ glutInitWindowSize(width, height)
+ glutInitWindowPosition(0, 0)
+ window = glutCreateWindow("Stereo Volume Rendering")
+
+
+ glutDisplayFunc(display)
+ glutIdleFunc(idle)
+ glutReshapeFunc(resize)
+ glutMouseFunc( mouseButton )
+ glutMotionFunc( mouseMotion )
+ glutKeyboardFunc(keyPressed)
+ init_gl(width, height)
+
+ # create texture for blitting to screen
+ create_texture(width, height)
+
+ import pycuda.gl.autoinit
+ import pycuda.gl
+ cuda_gl = pycuda.gl
+
+ create_PBO(width, height)
+ # ----- Load and Set Volume Data -----
+
+ density_grid = np.load("/home/bogert/dd150_log_densities.npy")
+
+ mi, ma= 21.5, 24.5
+ bins = 5000
+ tf = ColorTransferFunction( (mi, ma), bins)
+ tf.map_to_colormap(mi, ma, colormap="algae", scale_func = scale_func)
+
+ ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (width, height), tf = tf))
+
+ ts.get_raycaster().set_sample_size(0.01)
+ ts.get_raycaster().set_max_samples(5000)
+ ts.update()
+
+ glutMainLoop()
+
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, np.abs((v)-ma)/np.abs(mi-ma) + 0.0)
+
+# Print message to console, and kick off the main to get it rolling.
+if __name__ == "__main__":
+ print "Hit ESC key to quit, 'a' to toggle animation, and 'e' to toggle cuda"
+ main()
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 doc/source/developing/developing.rst
--- a/doc/source/developing/developing.rst
+++ b/doc/source/developing/developing.rst
@@ -74,6 +74,8 @@
this manner, but still want to contribute, please consider creating an external
package, which we'll happily link to.
+.. _requirements-for-code-submission:
+
Requirements for Code Submission
++++++++++++++++++++++++++++++++
@@ -88,22 +90,22 @@
* New Features
* New unit tests (possibly new answer tests) (See :ref:`testing`)
- * Docstrings for public API
- * Addition of new feature to the narrative documentation
- * Addition of cookbook recipe
+ * Docstrings in the source code for the public API
+ * Addition of new feature to the narrative documentation (See :ref:`writing_documentation`)
+ * Addition of cookbook recipe (See :ref:`writing_documentation`)
* Issue created on issue tracker, to ensure this is added to the changelog
* Extension or Breakage of API in Existing Features
- * Update existing narrative docs and docstrings
- * Update existing cookbook recipes
+ * Update existing narrative docs and docstrings (See :ref:`writing_documentation`)
+ * Update existing cookbook recipes (See :ref:`writing_documentation`)
* Modify of create new unit tests (See :ref:`testing`)
* Issue created on issue tracker, to ensure this is added to the changelog
* Bug fixes
* Unit test is encouraged, to ensure breakage does not happen again in the
- future.
+ future. (See :ref:`testing`)
* Issue created on issue tracker, to ensure this is added to the changelog
When submitting, you will be asked to make sure that your changes meet all of
@@ -178,14 +180,14 @@
$ python2.7 setup.py build --compiler=mingw32 install
+.. _sharing-changes:
+
Making and Sharing Changes
++++++++++++++++++++++++++
The simplest way to submit changes to yt is to commit changes in your
``$YT_DEST/src/yt-hg`` directory, fork the repository on BitBucket, push the
-changesets to your fork, and then issue a pull request. If you will be
-developing much more in-depth features for yt, you will also
-likely want to edit the paths in your
+changesets to your fork, and then issue a pull request.
Here's a more detailed flowchart of how to submit changes.
@@ -224,25 +226,72 @@
hg push https://bitbucket.org/YourUsername/yt/
- #. Update your pull request by visiting
- https://bitbucket.org/YourUsername/yt/pull-request/new
+ #. Your pull request will be automatically updated.
.. _writing_documentation:
How to Write Documentation
-++++++++++++++++++++++++++
+--------------------------
-The process for writing documentation is identical to the above, except that
-you're modifying source files in the doc directory (i.e. ``$YT_DEST/src/yt-hg/doc``)
-instead of the src directory (i.e. ``$YT_DEST/src/yt-hg/yt``) of the yt repository.
+Writing documentation is one of the most important but often overlooked tasks
+for increasing yt's impact in the community. It is the way in which the
+world will understand how to use our code, so it needs to be done concisely
+and understandably. Typically, when a developer submits some piece of code
+with new functionality, she should also include documentation on how to use
+that functionality (as per :ref:`requirements-for-code-submission`).
+Depending on the nature of the code addition, this could be a new narrative
+docs section describing how the new code works and how to use it, it could
+include a recipe in the cookbook section, or it could simply be adding a note
+in the relevant docs text somewhere.
+
+The documentation exists in the main mercurial code repository for yt in the
+``doc`` directory (i.e. ``$YT_DEST/src/yt-hg/doc/source`` on systems installed
+using the installer script). It is organized hierarchically into the main
+categories of:
+
+ * Visualizing
+ * Analyzing
+ * Examining
+ * Cookbook
+ * Bootcamp
+ * Developing
+ * Reference
+ * Help
+
+You will have to figure out where your new/modified doc fits into this, but
+browsing through the pre-built documentation is a good way to sort that out.
+
All the source for the documentation is written in
-`Sphinx <http://sphinx-doc.org/>`_, which uses ReST for markup.
+`Sphinx <http://sphinx-doc.org/>`_, which uses ReST for markup. ReST is very
+straightforward to markup in a text editor, and if you are new to it, we
+recommend just using other .rst files in the existing yt documentation as
+templates or checking out the
+`ReST reference documentation <http://sphinx-doc.org/rest.html>`_.
-Cookbook recipes go in ``source/cookbook/`` and must be added to one of the
-``.rst`` files in that directory.
+New cookbook recipes (see :ref:`cookbook`) are very helpful for the community
+as they provide simple annotated recipes on how to use specific functionality.
+To add one, create a concise python script which demonstrates some
+functionality and pare it down to its minimum. Add some comment lines to
+describe what it is that you're doing along the way. Place this ``.py`` file
+in the ``source/cookbook/`` directory, and then link to it explicitly in one
+of the relevant ``.rst`` files in that directory (e.g. ``complex_plots.rst``,
+``cosmological_analysis.rst``, etc.), and add some description of what the script
+actually does. We recommend that you use one of the
+`sample data sets <http://yt-project.org/data>`_ in your recipe. When the full
+docs are built, each of the cookbook recipes are executed dynamically on
+a system which has access to all of the sample datasets. Any output images
+generated by your script will then be attached inline in the built documentation
+directly following your script.
-For more information on how to build the documentation to make sure it looks
-the way you expect it to after modifying it, see :ref:`docs_build`.
+After you have made your modifications to the docs, you will want to make sure
+that they render the way you expect them to render. For more information on
+this, see the section on :ref:`docs_build`. Unless you're contributing cookbook
+recipes or notebooks which require a dynamical build, you can probably get
+away with just doing a 'quick' docs build.
+
+When you have completed your documentation additions, commit your changes
+to your repository and make a pull request in the same way you would contribute
+a change to the codebase, as described in the section on :ref:`sharing-changes`.
How To Get The Source Code For Editing
--------------------------------------
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 doc/source/visualizing/volume_rendering.rst
--- a/doc/source/visualizing/volume_rendering.rst
+++ b/doc/source/visualizing/volume_rendering.rst
@@ -466,3 +466,90 @@
your homogenized volume to then be passed in to the camera. A sample usage is shown
in :ref:`cookbook-amrkdtree_downsampling`.
+Hardware Volume Rendering on NVidia Graphics cards
+--------------------------------------------------
+.. versionadded:: 3.0
+
+Theia is a hardware volume renderer that takes advantage of NVidias CUDA language
+to peform ray casting with GPUs instead of the CPU.
+
+Only unigrid rendering is supported, but yt provides a grid mapping function
+ to get unigrid data from amr or sph formats :
+ :ref:`cookbook-amrkdtree_to_uniformgrid`.
+
+System Requirements
+-------------------
+.. versionadded:: 3.0
+
+Nvidia graphics card - The memory limit of the graphics card sets the limit
+ on the size of the data source.
+
+CUDA 5 or later and
+
+The environment variable CUDA_SAMPLES must be set pointing to
+the common/inc samples shipped with CUDA. The following shows an example
+in bash with CUDA 5.5 installed in /usr/local :
+
+export CUDA_SAMPLES=/usr/local/cuda-5.5/samples/common/inc
+
+PyCUDA must also be installed to use Theia.
+
+PyCUDA can be installed following these instructions :
+
+ git clone --recursive http://git.tiker.net/trees/pycuda.git
+
+ python configure.py
+ python setup.py install
+
+
+Tutorial
+--------
+.. versionadded:: 3.0
+
+Currently rendering only works on uniform grids. Here is an example
+on a 1024 cube of float32 scalars.
+
+.. code-block:: python
+ from yt.visualization.volume_rendering.theia.scene import TheiaScene
+ from yt.visualization.volume_rendering.algorithms.front_to_back import FrontToBackRaycaster
+ import numpy as np
+
+ #load 3D numpy array of float32
+ volume = np.load("/home/bogert/log_densities_1024.npy")
+
+ scene = TheiaScene( volume = volume, raycaster = FrontToBackRaycaster() )
+
+ scene.camera.rotateX(1.0)
+ scene.update()
+
+ surface = scene.get_results()
+ #surface now contains an image array 2x2 int32 rbga values
+
+.. _the-theiascene-interface:
+
+The TheiaScene Interface
+--------------------
+.. versionadded:: 3.0
+
+A TheiaScene object has been created to provide a high level entry point for
+controlling the raycaster's view onto the data. The class
+:class:`~yt.visualization.volume_rendering.theia.TheiaScene` encapsulates
+ a Camera object and a TheiaSource that intern encapsulates
+a volume. The :class:`~yt.visualization.volume_rendering.theia.Camera`
+provides controls for rotating, translating, and zooming into the volume.
+Using the :class:`~yt.visualization.volume_rendering.theia.TheiaSource`
+automatically transfers the volume to the graphic's card texture memory.
+
+Example Cookbooks
+---------------
+
+OpenGL Example for interactive volume rendering:
+:ref:`cookbook-opengl_volume_rendering`.
+
+OpenGL Stereoscopic Example :
+.. warning:: Frame rate will suffer significantly from stereoscopic rendering.
+ ~2x slower since the volume must be rendered twice.
+:ref:`cookbook-opengl_stereo_volume_rendering`.
+
+Pseudo-Realtime video rendering with ffmpeg :
+:ref:`cookbook-ffmpeg_volume_rendering`.
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 yt/data_objects/octree_subset.py
--- a/yt/data_objects/octree_subset.py
+++ b/yt/data_objects/octree_subset.py
@@ -102,10 +102,11 @@
new_shape = (nz, nz, nz, n_oct, 3)
else:
raise RuntimeError
- # This will retain units now.
- arr.shape = new_shape
- if not arr.flags["F_CONTIGUOUS"]:
- arr = arr.reshape(new_shape, order="F")
+ # Note that if arr is already F-contiguous, this *shouldn't* copy the
+ # data. But, it might. However, I can't seem to figure out how to
+ # make the assignment to .shape, which *won't* copy the data, make the
+ # resultant array viewed in Fortran order.
+ arr = arr.reshape(new_shape, order="F")
return arr
_domain_ind = None
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
@@ -0,0 +1,99 @@
+//-----------------------------------------------------------------------------
+// Copyright (c) 2014, yt Development Team.
+//
+// Distributed under the terms of the Modified BSD License.
+//
+// The full license is in the file COPYING.txt, distributed with this software.
+//-----------------------------------------------------------------------------
+
+#include <helpers/cuda.cu>
+
+//A 3d texture holding scalar values to be volume rendered
+texture<float , cudaTextureType3D, cudaReadModeElementType> volume;
+//A 1d texture holding color transfer function values to act on raycaster results
+texture<float4, cudaTextureType1D, cudaReadModeElementType> transfer;
+
+
+extern "C"
+__global__ void front_to_back(uint *buffer, float *modelviewmatrix, int buffer_w,
+ int buffer_h, int max_steps, float density,
+ float offset, float scale,
+ float brightness, float step_size
+ ) {
+ //Rays will terminate when opacity_threshold is reached
+ const float opacity_threshold = 0.95f;
+
+ //We don't use the far and near clipping information from the modelviewmatrix
+ float3x4 matrix = mat_array_to_3x4(modelviewmatrix);
+
+ //The X,Y coordinate of the surface (image) array
+ int x = BLOCK_X;
+ int y = BLOCK_Y;
+
+ if ((x >= buffer_w) || (y >= buffer_h)) return;
+
+ float u = COORD_STD(x, buffer_w);
+ float v = COORD_STD(y, buffer_h);
+
+ // calculate eye ray in world space
+ Ray eye_ray = make_ray_from_eye(matrix, u, v);
+
+ // find intersection with box
+ float tnear, tfar;
+ int hit = intersect_std_box(eye_ray, &tnear, &tfar);
+
+ // If the ray misses the box set the coloro to black
+ if (!hit) {
+ buffer[y*buffer_w + x] = 0;
+ return;
+ }
+
+ if (tnear < 0.0f) tnear = 0.0f; // clamp to near plane
+
+ // march along ray from front to back, accumulating color
+ float4 sum = make_float4(0.0f);
+ float t = tnear;
+ float3 pos = eye_ray.o + eye_ray.d*tnear;
+ float3 step = eye_ray.d*step_size;
+
+ // Cast Rays
+ float4 col = make_float4(0.0,0.0,0.0,0.0);
+ float top = (1.0/scale) + offset;
+ for (int i=0; i<=max_steps; i++) {
+ float sample = tex3D(volume, pos.x*0.5f+0.5f, pos.y*0.5f+0.5f, pos.z*0.5f+0.5f);
+
+
+ if (sample > offset) {
+ col = tex1D(transfer, (sample-offset)*scale);
+ }
+ else {
+ col = make_float4(0.0,0.0,0.0,0.0);
+ }
+
+ col.w *= density;
+
+ // pre-multiply alpha
+ col.x *= col.w;
+ col.y *= col.w;
+ col.z *= col.w;
+
+ // "over" operator for front-to-back blending
+ sum = sum + col*(1.0f - sum.w);
+
+ // exit early if opaque
+ if (sum.w > opacity_threshold)
+ break;
+
+ // Increment step size and position
+ t += step_size;
+ pos += step;
+
+ // If ray has cast too far, end
+ if (t > tfar) break;
+ }
+
+ sum *= brightness;
+
+ buffer[SCREEN_XY(x,y,buffer_w)] = rgbaFloatToInt(sum);
+}
+
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
@@ -0,0 +1,283 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+
+import pycuda.driver as drv
+import yt.visualization.volume_rendering.theia.helpers.cuda as cdh
+import numpy as np
+
+from yt.visualization.volume_rendering.theia.transfer.linear_transfer import LinearTransferFunction
+from yt.visualization.volume_rendering.theia.transfer.helper import *
+
+
+from yt.visualization.volume_rendering.theia.volumes.unigrid import Unigrid
+from yt.visualization.volume_rendering.theia.surfaces.array_surface import ArraySurface
+import os
+
+class FrontToBackRaycaster:
+ r"""
+ This is the python wrapper for the CUDA raycaster. This
+ raycaster can act on unigrid data only. Use yt to map
+ AMR data to unigrid prior to using this algorithm.
+
+
+ Parameters
+ ----------
+ matrix : 4x4 numpy.matriix
+ ModelViewMatrix defines view onto volume
+
+ surface: yt.visualization.volume_rendering.theia.surface.array_surface.ArraySurface
+ The surface to contain the image information from
+ the results of the raycasting
+
+ volume: 3D numpy array Float32
+ Scalar values to be volume rendered
+
+ tf : yt.visualization.volume_rendering.transfer_functions.ColorTransferFunction
+ Transfer function defining how the raycasting results are
+ to be colored.
+
+ size : Tuple of 2 Floats
+ If no surface is specified creates an ArraySurface of the specified size.
+
+ """
+
+ def __init__(self, matrix = None, surface = None,
+ volume = None, tf = None,
+ size = (864, 480)) :
+
+ # kernel module
+ self.kernel_module = cdh.source_file(
+ self.base_directory(__file__, "front_to_back.cu"))
+
+ #kernel function definitions
+ self.cast_rays_identifier = \
+ self.kernel_module.get_function("front_to_back")
+ self.transfer_identifier = \
+ self.kernel_module.get_texref("transfer")
+ self.volume_identifier = \
+ self.kernel_module.get_texref("volume")
+
+ self.set_matrix(matrix)
+
+ if (surface == None) :
+ self.set_surface(ArraySurface(size = size))
+ else :
+ self.set_surface(surface)
+
+
+ self.volume = None
+ self.set_transfer(tf)
+
+ self.set_sample_size()
+ self.set_max_samples()
+ self.set_opacity()
+ self.set_brightness()
+
+ """
+ Envoke the ray caster to cast rays
+ """
+ def __call__(self):
+ self.cast()
+ """
+ Returns
+ ----------
+ A 2d numpy array containing the image of the
+ volumetric rendering
+ """
+ def get_surface(self) :
+ return self.surface.get_array()
+
+ """
+ Returns
+ ----------
+ The sample size the ray caster is
+ configured to take.
+ """
+ def get_sample_size(self) :
+ return self.sample_size
+
+ """
+ Returns
+ ----------
+ The Max number of samples per ray the ray caster is
+ configured to take.
+ """
+ def get_max_samples(self) :
+ return self.max_samples
+
+ """
+ Returns
+ ----------
+ The Global density scalar.
+ """
+ def get_opacity(self) :
+ return self.opacity
+
+ """
+ Returns
+ ----------
+ The Global brightness scalar.
+
+ """
+ def get_brightness(self):
+ return self.brightness
+
+ """
+ Parameters
+ ----------
+ size : scalra Float
+ The distance between each sample, a smaller size will result in
+ more samples and can cause performance loss.
+ """
+ def set_sample_size(self, size = 0.01) :
+ self.sample_size = size
+
+ """
+ Parameters
+ ----------
+ max : scalar Int
+ The limit on the number of samples the ray caster will
+ take per ray.
+ """
+ def set_max_samples(self, max = 5000) :
+ self.max_samples = max
+
+ """
+ Parameters
+ ----------
+ scale : scalar Float
+ Global multiplier on volume data
+ """
+ def set_opacity(self, scale = 0.05) :
+ self.opacity = scale
+
+ """
+ Parameters
+ ----------
+ brightness : scalar Float
+ Global multiplier on returned alpha values
+ """
+ def set_brightness(self, brightness = 1.0):
+ self.brightness = brightness
+
+ """
+ Causes the ray caster to act on the volume data
+ and puts the results in the surface array.
+ """
+ def cast(self):
+ w, h = self.surface.bounds
+
+ self.cast_rays_identifier(np.uint64(self.surface.device_ptr),
+ drv.In(self.matrix),
+ np.int32(w), np.int32(h),
+ np.int32(self.max_samples),
+ np.float32(self.opacity),
+ np.float32(self.transfer.offset),
+ np.float32(self.transfer.scale),
+ np.float32(self.brightness),
+ np.float32(self.sample_size),
+ block=self.block,
+ grid=self.grid
+ )
+
+
+ """
+ Set the ModelViewMatrix, this does not send data to the GPU
+ Parameters
+ ----------
+ matrix : 4x4 numpy.matrix
+ ModelViewMatrix
+
+ """
+ def set_matrix(self, matrix):
+ self.matrix = matrix
+
+ """
+ Setup the image array for ray casting results
+ Parameters
+ ----------
+ surface : yt.visualization..volume_rendering.theia.surfaces.ArraySurface
+ Surface to contain results of the ray casting
+ """
+ def set_surface(self, surface = None, block_size = 32):
+ if surface == None:
+ self.surface = None
+ return
+
+ self.surface = surface
+ self.grid = cdh.block_estimate(block_size, self.surface.bounds)
+ self.block = (block_size, block_size, 1)
+
+ """
+ This function will convert a 3d numpy array into a unigrid volume
+ and move the data to the gpu texture memory
+ Parameters
+ ----------
+ volume : 3D numpy array float32
+ Contains scalar volumes to be acted on by the ray caster
+
+ """
+ def send_volume_to_gpu(self, volume = None) :
+ if (volume != None) :
+ self.set_volume(Unigrid(array = volume, allocate = True))
+
+ """
+ Parameters
+ ----------
+ volume : yt.visualization..volume_rendering.theia.volumes.Unigrid
+ Contains scalar volumes to be acted on by the ray caster
+
+ """
+ def set_volume(self, volume = None):
+ if volume == None:
+ self.volume = None
+ return
+
+ self.volume = volume
+ self.volume_identifier.set_flags(self.volume.flags)
+ self.volume_identifier.set_filter_mode(self.volume.filter_mode)
+ self.volume_identifier.set_address_mode(0, self.volume.address_mode)
+ self.volume_identifier.set_address_mode(1, self.volume.address_mode)
+ self.volume_identifier.set_array(self.volume.cuda_volume_array)
+
+ """
+ Parameters
+ ----------
+ tf : yt.visualization.volume_rendering.transfer_functions.ColorTransferFunction
+ Used to color the results of the raycasting
+
+ """
+ def set_transfer(self, tf = None):
+ if tf == None:
+ self.transfer = None
+ return
+
+ self.transfer = LinearTransferFunction(arrays = yt_to_rgba_transfer(tf), range = tf.x_bounds)
+ self.transfer_identifier.set_flags(self.transfer.flags)
+ self.transfer_identifier.set_filter_mode(self.transfer.filter_mode)
+ self.transfer_identifier.set_address_mode(0, self.transfer.address_mode)
+ self.transfer_identifier.set_array(self.transfer.cuda_transfer_array)
+
+ """
+ Attach the base directory path to the desired source file.
+ Parameters
+ ----------
+ dir : string
+ Directory where source file is located
+
+ file : string
+ Source file name
+
+
+ """
+ def base_directory(self, dir, file):
+ base = os.path.dirname(dir)
+ src = os.path.join(base, file)
+ return src
+
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 yt/visualization/volume_rendering/theia/camera.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/camera.py
@@ -0,0 +1,139 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+#Matix helpers for camera control
+import yt.visualization.volume_rendering.theia.helpers.matrix_transformation as mth
+
+import numpy as np
+
+class Camera:
+ r"""
+ This is a basic camera intended to be a base class.
+ The camera provies utilities for controling the view point onto the data.
+
+
+ Parameters
+ ----------
+ pos : numpy array
+ The 3d position vector of the camera
+ rot : numpy array
+ The 3d rotation vector of the camera
+ scale : float
+ The scalar acting on camera position to zoom
+
+ Example:
+
+ from yt.visualization.volume_rendering.cameras.camera import Camera
+
+ cam = Camera()
+
+ cam.zoom(10.0)
+ cam.rotateX(np.pi)
+ cam.rotateY(np.pi)
+ cam.rotateZ(np.pi)
+ cam.translateX(-0.5)
+ cam.translateY(-0.5)
+ cam.translateZ(-0.5)
+
+ view = cam.get_matrix()
+
+ """
+ def __init__(self, rot = np.array([0.0, 0.0, 0.0]), scale = 1.0,
+ pos = np.array([0.0, 0.0, 0.0]), ):
+ #Values to control camera
+ self.rot = rot
+ self.scale = scale
+ self.pos = pos
+
+ #TODO : users should have control over perspective frustrum
+ self.stdmatrix = mth.perspective( [0.0, 0.0, 0.25] )
+
+
+ def zoom(self, percentage = 10.0):
+ r"""This will zoom the camera by percentage specified.
+
+ Parameters
+ ----------
+ percentage : float
+ If percentage is postive the camera will zoom in, if negative
+ then the camera will zoom out. Cannot exceed 100%.
+ """
+ if (percentage > 100.0) :
+ percentage = 100.0
+ self.scale = ((100.0 - percentage)/100.0)*self.scale
+
+ def rotateX(self, degrees = 0.1):
+ r"""This will rotate the camera around its X axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[0] += degrees
+
+ def rotateY(self, degrees = 0.1):
+ r"""This will rotate the camera around its Y axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[1] += degrees
+
+ def rotateZ(self, degrees = 0.1):
+ r"""This will rotate the camera around its Z axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[2] += degrees
+
+ def translateX(self, dist = 0.1):
+ r"""This will move the camera's position along its X axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[0] += dist
+
+ def translateY(self, dist = 0.1):
+ r"""This will move the camera's position along its Y axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[1] += dist
+
+ def translateZ(self, dist = 0.1):
+ r"""This will move the camera's position along its Z axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[2] += dist
+
+ def get_matrix(self):
+ r"""This will calculate the curent modelview matrix
+ from the camera's position, rotation, and zoom scale.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ return mth.rotate_x(self.rot[0]) * mth.rotate_y(self.rot[1]) * mth.rotate_z(self.rot[2]) * mth.scale((self.scale, self.scale, self.scale)) * mth.perspective([0.0, 0.0, 0.25]) *mth.translate((self.pos[0], self.pos[1], self.pos[2]))
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 yt/visualization/volume_rendering/theia/helpers/cuda.cu
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda.cu
@@ -0,0 +1,189 @@
+//-----------------------------------------------------------------------------
+// Copyright (c) 2014, yt Development Team.
+//
+// Distributed under the terms of the Modified BSD License.
+//
+// The full license is in the file COPYING.txt, distributed with this software.
+//-----------------------------------------------------------------------------
+
+
+#include <helper_cuda.h>
+#include <helper_math.h>
+
+#define BLOCK_X (blockIdx.x*blockDim.x + threadIdx.x)
+#define BLOCK_Y (blockIdx.y*blockDim.y + threadIdx.y)
+#define SCREEN_XY(a, b, w) (b*w+a)
+#define COORD_STD(x, w) ((x / (float) w)*2.0f-1.0f)
+
+#define COLOR_BLACK make_float3(0.0f, 0.0f, 0.0f)
+
+typedef struct {
+ float4 m[3];
+} float3x4;
+
+typedef struct {
+ float4 m[4];
+} float4x4;
+
+
+struct Ray {
+ float3 o; // origin
+ float3 d; // direction
+};
+
+
+
+__device__ float4 mul(const float3x4 &M, const float4 &v);
+__device__ float3 mul(const float3x4 &M, const float3 &v);
+__device__ float3x4 mat_array_to_3x4(float *m);
+__device__ Ray make_ray_from_eye(float3x4 m, float u, float v);
+__device__ int intersect_box(Ray r, float3 boxmin, float3 boxmax,
+ float *tnear, float *tfar);
+__device__ int intersect_std_box(Ray r, float *tnear, float *tfar);
+
+__device__ float mag(const float3 &v);
+__device__ float avg(const float3 &v);
+__device__ uint rgbaFloatToInt(float4 rgba);
+__device__ uint rgbFloatToInt(float3 rgb);
+__device__ float4 intToRGBAFloat(uint irgba);
+__device__ float3 intToRGBFloat(uint irgba);
+
+
+
+__device__ uint rgbFloatToInt(float3 rgba)
+{
+ rgba.x = __saturatef(rgba.x); // clamp to [0.0, 1.0]
+ rgba.y = __saturatef(rgba.y);
+ rgba.z = __saturatef(rgba.z);
+ return (uint(rgba.z*255)<<16) | (uint(rgba.y*255)<<8) | uint(rgba.x*255);
+}
+
+__device__ uint rgbaFloatToInt(float4 rgba)
+{
+ rgba.x = __saturatef(rgba.x); // clamp to [0.0, 1.0]
+ rgba.y = __saturatef(rgba.y);
+ rgba.z = __saturatef(rgba.z);
+ rgba.w = __saturatef(rgba.w);
+ return (uint(rgba.w*255)<<24) | (uint(rgba.z*255)<<16) | (uint(rgba.y*255)<<8) | uint(rgba.x*255);
+}
+
+__device__ float4 intToRGBAFloat(uint irgba)
+{
+ float4 rgba;
+ rgba.w = (irgba >> 24 & 0xFF)/255;
+ rgba.z = (irgba >> 16 & 0xFF)/255;
+ rgba.y = (irgba >> 8 & 0xFF)/255;
+ rgba.x = (irgba & 0xFF)/255;
+ return rgba;
+}
+
+__device__ float3 intToRGBFloat(uint irgba)
+{
+ float3 rgb;
+ rgb.z = (irgba >> 16 & 0xFF)/255;
+ rgb.y = (irgba >> 8 & 0xFF)/255;
+ rgb.x = (irgba & 0xFF)/255;
+ return rgb;
+}
+
+
+__device__ float3x4 mat_array_to_3x4(float *m) {
+ float3x4 matrix;
+ matrix.m[0] = make_float4(m[0], m[1], m[2], m[3]);
+ matrix.m[1] = make_float4(m[4], m[5], m[6], m[7]);
+ matrix.m[2] = make_float4(m[8], m[9], m[10], m[11]);
+ return matrix;
+}
+
+__device__ Ray make_ray_from_eye(float3x4 m, float u, float v) {
+ Ray eyeRay;
+ eyeRay.o = make_float3(mul(m, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
+ eyeRay.d = normalize(make_float3(u, v, -2.0f));
+ eyeRay.d = mul(m, eyeRay.d);
+ return eyeRay;
+}
+
+__device__
+int intersect_box(Ray r, float3 boxmin, float3 boxmax, float *tnear, float *tfar)
+{
+ // compute intersection of ray with all six bbox planes
+ float3 invR = make_float3(1.0f) / r.d;
+ float3 tbot = invR * (boxmin - r.o);
+ float3 ttop = invR * (boxmax - r.o);
+
+ // re-order intersections to find smallest and largest on each axis
+ float3 tmin = fminf(ttop, tbot);
+ float3 tmax = fmaxf(ttop, tbot);
+
+ // find the largest tmin and the smallest tmax
+ float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
+ float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));
+
+ *tnear = largest_tmin;
+ *tfar = smallest_tmax;
+
+ return smallest_tmax > largest_tmin;
+}
+
+__device__
+int intersect_std_box(Ray r, float *tnear, float *tfar)
+{
+ const float3 boxmin = make_float3(-1.0f, -1.0f, -1.0f);
+ const float3 boxmax = make_float3( 1.0f, 1.0f, 1.0f);
+
+ // compute intersection of ray with all six bbox planes
+ float3 invR = make_float3(1.0f) / r.d;
+ float3 tbot = invR * (boxmin - r.o);
+ float3 ttop = invR * (boxmax - r.o);
+
+ // re-order intersections to find smallest and largest on each axis
+ float3 tmin = fminf(ttop, tbot);
+ float3 tmax = fmaxf(ttop, tbot);
+
+ // find the largest tmin and the smallest tmax
+ float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
+ float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));
+
+ *tnear = largest_tmin;
+ *tfar = smallest_tmax;
+
+ return smallest_tmax > largest_tmin;
+}
+
+__device__
+float4 mul(const float3x4 &M, const float4 &v)
+{
+ float4 r;
+ r.x = dot(v, M.m[0]);
+ r.y = dot(v, M.m[1]);
+ r.z = dot(v, M.m[2]);
+ r.w = 1.0f;
+ return r;
+}
+
+// transform vector by matrix with translation
+
+__device__
+float3 mul(const float3x4 &M, const float3 &v)
+{
+ float3 r;
+ r.x = dot(v, make_float3(M.m[0]));
+ r.y = dot(v, make_float3(M.m[1]));
+ r.z = dot(v, make_float3(M.m[2]));
+ return r;
+}
+
+
+__device__
+float mag(const float3 &v) {
+ return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
+}
+
+
+__device__
+float avg(const float3 &v) {
+ return (v.x + v.y + v.z) / 3.0;
+}
+
+
+
diff -r b6f0110bc38ec3b9400146d0609e71383f1ee922 -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 yt/visualization/volume_rendering/theia/helpers/cuda.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda.py
@@ -0,0 +1,58 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+
+import pycuda.driver as drv
+from pycuda.compiler import SourceModule
+import pycuda.gpuarray as gpuarray
+import pycuda.autoinit
+import math
+import os
+
+def sample_path():
+ #return '/pfs/sw/cuda-5.5/samples/common/inc'
+ return os.environ.get('CUDA_SAMPLES')
+
+def cuda_helpers_path():
+ return os.path.dirname(__file__) + "/../"
+
+def source_file(path):
+ return SourceModule(open(path).read(),
+ include_dirs=[sample_path(),cuda_helpers_path()],
+ no_extern_c=True, keep=True)
+
+def block_estimate(thread_size, shape):
+ (x, y) = shape
+ return (int(math.ceil(float(x) / thread_size)), int(math.ceil(float(y) / thread_size)), 1)
+
+def np3d_to_device_array(np_array, allow_surface_bind=True):
+ d, h, w = np_array.shape
+
+ descr = drv.ArrayDescriptor3D()
+ descr.width = w
+ descr.height = h
+ descr.depth = d
+ descr.format = drv.dtype_to_array_format(np_array.dtype)
+ descr.num_channels = 1
+ descr.flags = 0
+
+ if allow_surface_bind:
+ descr.flags = drv.array3d_flags.SURFACE_LDST
+
+ device_array = drv.Array(descr)
+
+ copy = drv.Memcpy3D()
+ copy.set_src_host(np_array)
+ copy.set_dst_array(device_array)
+ copy.width_in_bytes = copy.src_pitch = np_array.strides[1]
+ copy.src_height = copy.height = h
+ copy.depth = d
+
+ copy()
+
+ return device_array
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/yt_analysis/yt/commits/87b8587c79c8/
Changeset: 87b8587c79c8
Branch: yt-3.0
User: jzuhone
Date: 2014-06-18 18:45:23
Summary: Keyword arguments should go in this way.
Affected #: 1 file
diff -r 8f8bb4c481d1efefacb549ecd68ea3f30b7e5343 -r 87b8587c79c8f2999c883bb7e551b32142212ba9 yt/frontends/fits/tests/test_outputs.py
--- a/yt/frontends/fits/tests/test_outputs.py
+++ b/yt/frontends/fits/tests/test_outputs.py
@@ -26,7 +26,7 @@
m33 = "radio_fits/m33_hi.fits"
@requires_pf(m33, big_data=True)
def test_m33():
- pf = data_dir_load(m33, cls=FITSDataset, nan_mask=0.0)
+ pf = data_dir_load(m33, cls=FITSDataset, kwargs={"nan_mask":0.0})
yield assert_equal, str(pf), "m33_hi.fits"
for test in small_patch_amr(m33, _fields):
test_m33.__name__ = test.description
@@ -37,7 +37,7 @@
grs = "radio_fits/grs-50-cube.fits"
@requires_pf(grs)
def test_grs():
- pf = data_dir_load(grs, cls=FITSDataset, nan_mask=0.0)
+ pf = data_dir_load(grs, cls=FITSDataset, kwargs={"nan_mask":0.0})
yield assert_equal, str(pf), "grs-50-cube.fits"
for test in small_patch_amr(grs, _fields):
test_grs.__name__ = test.description
https://bitbucket.org/yt_analysis/yt/commits/36f776e37ed9/
Changeset: 36f776e37ed9
Branch: yt-3.0
User: jzuhone
Date: 2014-06-18 19:04:24
Summary: More flexible options for the center of the testing sphere and the weight field for the projection tests.
Affected #: 1 file
diff -r 87b8587c79c8f2999c883bb7e551b32142212ba9 -r 36f776e37ed9c7e18ce3a2bdfa51068c114b1c59 yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -691,32 +691,32 @@
else:
return ftrue
-def small_patch_amr(pf_fn, fields):
+def small_patch_amr(pf_fn, fields, input_center="max", input_weight="density"):
if not can_run_pf(pf_fn): return
- dso = [ None, ("sphere", ("max", (0.1, 'unitary')))]
+ dso = [ None, ("sphere", (input_center, (0.1, 'unitary')))]
yield GridHierarchyTest(pf_fn)
yield ParentageRelationshipsTest(pf_fn)
for field in fields:
yield GridValuesTest(pf_fn, field)
for axis in [0, 1, 2]:
for ds in dso:
- for weight_field in [None, "density"]:
+ for weight_field in [None, input_weight]:
yield ProjectionValuesTest(
pf_fn, axis, field, weight_field,
ds)
yield FieldValuesTest(
pf_fn, field, ds)
-def big_patch_amr(pf_fn, fields):
+def big_patch_amr(pf_fn, fields, input_center="max", input_weight="density"):
if not can_run_pf(pf_fn): return
- dso = [ None, ("sphere", ("max", (0.1, 'unitary')))]
+ dso = [ None, ("sphere", (input_center, (0.1, 'unitary')))]
yield GridHierarchyTest(pf_fn)
yield ParentageRelationshipsTest(pf_fn)
for field in fields:
yield GridValuesTest(pf_fn, field)
for axis in [0, 1, 2]:
for ds in dso:
- for weight_field in [None, "density"]:
+ for weight_field in [None, input_weight]:
yield PixelizedProjectionValuesTest(
pf_fn, axis, field, weight_field,
ds)
https://bitbucket.org/yt_analysis/yt/commits/874265672c37/
Changeset: 874265672c37
Branch: yt-3.0
User: jzuhone
Date: 2014-06-18 19:04:46
Summary: These lists were not set up correctly
Affected #: 1 file
diff -r 36f776e37ed9c7e18ce3a2bdfa51068c114b1c59 -r 874265672c372b25b870e98b40098a84cb2923c8 yt/frontends/fits/tests/test_outputs.py
--- a/yt/frontends/fits/tests/test_outputs.py
+++ b/yt/frontends/fits/tests/test_outputs.py
@@ -21,35 +21,35 @@
data_dir_load
from ..data_structures import FITSDataset
-_fields = ("intensity")
+_fields_m33 = ("intensity",)
m33 = "radio_fits/m33_hi.fits"
@requires_pf(m33, big_data=True)
def test_m33():
pf = data_dir_load(m33, cls=FITSDataset, kwargs={"nan_mask":0.0})
yield assert_equal, str(pf), "m33_hi.fits"
- for test in small_patch_amr(m33, _fields):
+ for test in small_patch_amr(m33, _fields_m33, input_center="c", input_weight="ones"):
test_m33.__name__ = test.description
yield test
-_fields = ("temperature")
+_fields_grs = ("temperature",)
grs = "radio_fits/grs-50-cube.fits"
@requires_pf(grs)
def test_grs():
pf = data_dir_load(grs, cls=FITSDataset, kwargs={"nan_mask":0.0})
yield assert_equal, str(pf), "grs-50-cube.fits"
- for test in small_patch_amr(grs, _fields):
+ for test in small_patch_amr(grs, _fields_grs, input_center="c", input_weight="ones"):
test_grs.__name__ = test.description
yield test
-_fields = ("x-velocity","y-velocity","z-velocity")
+_fields_vels = ("x-velocity","y-velocity","z-velocity")
vf = "UniformGrid/velocity_field_20.fits"
@requires_pf(vf)
def test_velocity_field():
pf = data_dir_load(bf, cls=FITSDataset)
yield assert_equal, str(pf), "velocity_field_20.fits"
- for test in small_patch_amr(vf, _fields):
+ for test in small_patch_amr(vf, _fields_vels):
test_velocity_field.__name__ = test.description
yield test
https://bitbucket.org/yt_analysis/yt/commits/65714677a970/
Changeset: 65714677a970
Branch: yt-3.0
User: jzuhone
Date: 2014-06-19 12:24:38
Summary: Replacing pf with ds, fixing a typo
Affected #: 1 file
diff -r 874265672c372b25b870e98b40098a84cb2923c8 -r 65714677a970e73c5667dde59a2cff42d6b10b4f yt/frontends/fits/tests/test_outputs.py
--- a/yt/frontends/fits/tests/test_outputs.py
+++ b/yt/frontends/fits/tests/test_outputs.py
@@ -17,7 +17,6 @@
from yt.utilities.answer_testing.framework import \
requires_pf, \
small_patch_amr, \
- big_patch_amr, \
data_dir_load
from ..data_structures import FITSDataset
@@ -26,8 +25,8 @@
m33 = "radio_fits/m33_hi.fits"
@requires_pf(m33, big_data=True)
def test_m33():
- pf = data_dir_load(m33, cls=FITSDataset, kwargs={"nan_mask":0.0})
- yield assert_equal, str(pf), "m33_hi.fits"
+ ds = data_dir_load(m33, cls=FITSDataset, kwargs={"nan_mask":0.0})
+ yield assert_equal, str(ds), "m33_hi.fits"
for test in small_patch_amr(m33, _fields_m33, input_center="c", input_weight="ones"):
test_m33.__name__ = test.description
yield test
@@ -37,19 +36,19 @@
grs = "radio_fits/grs-50-cube.fits"
@requires_pf(grs)
def test_grs():
- pf = data_dir_load(grs, cls=FITSDataset, kwargs={"nan_mask":0.0})
- yield assert_equal, str(pf), "grs-50-cube.fits"
+ ds = data_dir_load(grs, cls=FITSDataset, kwargs={"nan_mask":0.0})
+ yield assert_equal, str(ds), "grs-50-cube.fits"
for test in small_patch_amr(grs, _fields_grs, input_center="c", input_weight="ones"):
test_grs.__name__ = test.description
yield test
-_fields_vels = ("x-velocity","y-velocity","z-velocity")
+_fields_vels = ("velocity_x","velocity_y","velocity_z")
-vf = "UniformGrid/velocity_field_20.fits"
+vf = "UnigridData/velocity_field_20.fits"
@requires_pf(vf)
def test_velocity_field():
- pf = data_dir_load(bf, cls=FITSDataset)
- yield assert_equal, str(pf), "velocity_field_20.fits"
- for test in small_patch_amr(vf, _fields_vels):
+ ds = data_dir_load(vf, cls=FITSDataset)
+ yield assert_equal, str(ds), "velocity_field_20.fits"
+ for test in small_patch_amr(vf, _fields_vels, input_center="c", input_weight="ones"):
test_velocity_field.__name__ = test.description
yield test
https://bitbucket.org/yt_analysis/yt/commits/e17556a9d270/
Changeset: e17556a9d270
Branch: yt-3.0
User: jzuhone
Date: 2014-06-19 12:34:59
Summary: M33 test is taking too long...
Affected #: 1 file
diff -r 65714677a970e73c5667dde59a2cff42d6b10b4f -r e17556a9d270d158341c7fb5ca07507eedc9c624 yt/frontends/fits/tests/test_outputs.py
--- a/yt/frontends/fits/tests/test_outputs.py
+++ b/yt/frontends/fits/tests/test_outputs.py
@@ -20,17 +20,6 @@
data_dir_load
from ..data_structures import FITSDataset
-_fields_m33 = ("intensity",)
-
-m33 = "radio_fits/m33_hi.fits"
- at requires_pf(m33, big_data=True)
-def test_m33():
- ds = data_dir_load(m33, cls=FITSDataset, kwargs={"nan_mask":0.0})
- yield assert_equal, str(ds), "m33_hi.fits"
- for test in small_patch_amr(m33, _fields_m33, input_center="c", input_weight="ones"):
- test_m33.__name__ = test.description
- yield test
-
_fields_grs = ("temperature",)
grs = "radio_fits/grs-50-cube.fits"
https://bitbucket.org/yt_analysis/yt/commits/5088ddc41bd7/
Changeset: 5088ddc41bd7
Branch: yt-3.0
User: MatthewTurk
Date: 2014-06-19 14:52:56
Summary: Merged in jzuhone/yt-3.x/yt-3.0 (pull request #961)
FITS answer tests, fixing Windows test failures
Affected #: 4 files
diff -r 534e2b948eda6907fa640eabbf4f307e8a0360ef -r 5088ddc41bd7c54f01f22aa5a0c217728268dbe0 yt/frontends/fits/tests/test_outputs.py
--- a/yt/frontends/fits/tests/test_outputs.py
+++ b/yt/frontends/fits/tests/test_outputs.py
@@ -17,38 +17,27 @@
from yt.utilities.answer_testing.framework import \
requires_pf, \
small_patch_amr, \
- big_patch_amr, \
data_dir_load
+from ..data_structures import FITSDataset
-_fields = ("intensity")
-
-m33 = "radio_fits/m33_hi.fits"
- at requires_pf(m33, big_data=True)
-def test_m33():
- pf = data_dir_load(m33, nan_mask=0.0)
- yield assert_equal, str(pf), "m33_hi.fits"
- for test in small_patch_amr(m33, _fields):
- test_m33.__name__ = test.description
- yield test
-
-_fields = ("temperature")
+_fields_grs = ("temperature",)
grs = "radio_fits/grs-50-cube.fits"
@requires_pf(grs)
def test_grs():
- pf = data_dir_load(grs, nan_mask=0.0)
- yield assert_equal, str(pf), "grs-50-cube.fits"
- for test in small_patch_amr(grs, _fields):
+ ds = data_dir_load(grs, cls=FITSDataset, kwargs={"nan_mask":0.0})
+ yield assert_equal, str(ds), "grs-50-cube.fits"
+ for test in small_patch_amr(grs, _fields_grs, input_center="c", input_weight="ones"):
test_grs.__name__ = test.description
yield test
-_fields = ("x-velocity","y-velocity","z-velocity")
+_fields_vels = ("velocity_x","velocity_y","velocity_z")
-vf = "UniformGrid/velocity_field_20.fits"
+vf = "UnigridData/velocity_field_20.fits"
@requires_pf(vf)
def test_velocity_field():
- pf = data_dir_load(bf)
- yield assert_equal, str(pf), "velocity_field_20.fits"
- for test in small_patch_amr(vf, _fields):
+ ds = data_dir_load(vf, cls=FITSDataset)
+ yield assert_equal, str(ds), "velocity_field_20.fits"
+ for test in small_patch_amr(vf, _fields_vels, input_center="c", input_weight="ones"):
test_velocity_field.__name__ = test.description
yield test
diff -r 534e2b948eda6907fa640eabbf4f307e8a0360ef -r 5088ddc41bd7c54f01f22aa5a0c217728268dbe0 yt/units/yt_array.py
--- a/yt/units/yt_array.py
+++ b/yt/units/yt_array.py
@@ -563,7 +563,7 @@
unit_lut = pickle.loads(dataset.attrs['unit_registry'])
else:
unit_lut = None
-
+ f.close()
registry = UnitRegistry(lut=unit_lut, add_default_symbols=False)
return cls(data, units, registry=registry)
diff -r 534e2b948eda6907fa640eabbf4f307e8a0360ef -r 5088ddc41bd7c54f01f22aa5a0c217728268dbe0 yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -691,32 +691,32 @@
else:
return ftrue
-def small_patch_amr(pf_fn, fields):
+def small_patch_amr(pf_fn, fields, input_center="max", input_weight="density"):
if not can_run_pf(pf_fn): return
- dso = [ None, ("sphere", ("max", (0.1, 'unitary')))]
+ dso = [ None, ("sphere", (input_center, (0.1, 'unitary')))]
yield GridHierarchyTest(pf_fn)
yield ParentageRelationshipsTest(pf_fn)
for field in fields:
yield GridValuesTest(pf_fn, field)
for axis in [0, 1, 2]:
for ds in dso:
- for weight_field in [None, "density"]:
+ for weight_field in [None, input_weight]:
yield ProjectionValuesTest(
pf_fn, axis, field, weight_field,
ds)
yield FieldValuesTest(
pf_fn, field, ds)
-def big_patch_amr(pf_fn, fields):
+def big_patch_amr(pf_fn, fields, input_center="max", input_weight="density"):
if not can_run_pf(pf_fn): return
- dso = [ None, ("sphere", ("max", (0.1, 'unitary')))]
+ dso = [ None, ("sphere", (input_center, (0.1, 'unitary')))]
yield GridHierarchyTest(pf_fn)
yield ParentageRelationshipsTest(pf_fn)
for field in fields:
yield GridValuesTest(pf_fn, field)
for axis in [0, 1, 2]:
for ds in dso:
- for weight_field in [None, "density"]:
+ for weight_field in [None, input_weight]:
yield PixelizedProjectionValuesTest(
pf_fn, axis, field, weight_field,
ds)
diff -r 534e2b948eda6907fa640eabbf4f307e8a0360ef -r 5088ddc41bd7c54f01f22aa5a0c217728268dbe0 yt/utilities/lib/tests/test_ragged_arrays.py
--- a/yt/utilities/lib/tests/test_ragged_arrays.py
+++ b/yt/utilities/lib/tests/test_ragged_arrays.py
@@ -13,7 +13,7 @@
def test_index_unop():
np.random.seed(0x4d3d3d3)
- indices = np.arange(1000)
+ indices = np.arange(1000, dtype="int64")
np.random.shuffle(indices)
sizes = np.array([
200, 50, 50, 100, 32, 32, 32, 32, 32, 64, 376], dtype="int64")
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