[Yt-svn] yt: 2 new changesets
hg at spacepope.org
hg at spacepope.org
Wed May 5 07:58:36 PDT 2010
hg Repository: yt
details: yt/rev/6061669dba57
changeset: 1640:6061669dba57
user: Sam Skillman <samuel.skillman at colorado.edu>
date:
Tue May 04 18:15:32 2010 -0600
description:
Adding amr_kdtree extension to be used in volume rendering decomposition. At the moment it is working, but only supports the Density field and not any of the cool multivariate stuff, which is the next goal. Also modifying the north vector slightly which can be changed back if people don't like the new behavior. Now it always makes sure the north vector is in the up direction in the image plane, though it can point in/out of it. This seems to make rotations easier to handle for me.
hg Repository: yt
details: yt/rev/857c1f4280f5
changeset: 1641:857c1f4280f5
user: Sam Skillman <samuel.skillman at colorado.edu>
date:
Wed May 05 08:57:31 2010 -0600
description:
Multivariate now working in kd-tree ray cast.
diffstat:
.hgignore | 5 +
yt/extensions/amr_kdtree/__init__.py | 5 +
yt/extensions/amr_kdtree/amr_kdtree.py | 116 +++++++++++++++++++++++
yt/extensions/setup.py | 1 +
yt/extensions/volume_rendering/software_sampler.py | 115 ++++++++++++++++++++++-
yt/lagos/ParallelTools.py | 30 ++++++
6 files changed, 270 insertions(+), 2 deletions(-)
diffs (truncated from 354 to 300 lines):
diff -r 0f1c190e5e47 -r 857c1f4280f5 .hgignore
--- a/.hgignore Tue May 04 07:51:43 2010 -0700
+++ b/.hgignore Wed May 05 08:57:31 2010 -0600
@@ -1,3 +1,8 @@
syntax: glob
*.pyc
.*.swp
+*.o
+hdf5.cfg
+build
+dist
+yt.egg-info
\ No newline at end of file
diff -r 0f1c190e5e47 -r 857c1f4280f5 yt/extensions/amr_kdtree/__init__.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/extensions/amr_kdtree/__init__.py Wed May 05 08:57:31 2010 -0600
@@ -0,0 +1,5 @@
+"""
+
+Initialize amr_kdtree
+"""
+from amr_kdtree import *
diff -r 0f1c190e5e47 -r 857c1f4280f5 yt/extensions/amr_kdtree/amr_kdtree.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/extensions/amr_kdtree/amr_kdtree.py Wed May 05 08:57:31 2010 -0600
@@ -0,0 +1,116 @@
+import numpy as np
+from copy import *
+from yt.mods import *
+from time import time
+
+class AMRKDTree(object):
+ def __init__(self, pf, L_MAX=None, extra_dim_cost=0, le=None, re=None):
+ self.pf = pf
+ if L_MAX is None:
+ self.L_MAX = self.pf.hierarchy.max_level+1
+ else:
+ self.L_MAX = np.min([L_MAX,self.pf.hierarchy.max_level+1])
+
+ self.extra_dim_cost=extra_dim_cost
+ if le is None: self.domain_left_edge = np.array([0.0]*3)
+ else: self.domain_left_edge = le
+ if re is None: self.domain_right_edge = np.array([1.0]*3)
+ else: self.domain_right_edge = re
+ print 'Making kd tree from le ', self.domain_left_edge, 'to ', self.domain_right_edge
+
+ self.root_grids = pf.hierarchy.get_levels().next()
+
+ self.leaf_count = 0
+ self.current_parent = -1
+ self.dim = 0
+ tt1 = time()
+ self.tree = self.__build(self.root_grids, None, self.domain_left_edge, self.domain_right_edge)
+ print 'It took %e seconds to build the tree' % (time()-tt1)
+ self.total_cost = self.count_cost(self.tree)
+
+ class node(object):
+ pass
+
+ class leafnode(node):
+ def __init__(self,leaf_id, grid_id, leaf_l_corner, leaf_r_corner):
+ self.leaf_id = leaf_id
+ self.grid = grid_id
+ self.leaf_l_corner = leaf_l_corner
+ self.leaf_r_corner = leaf_r_corner
+ dds = self.grid.dds
+ gle = self.grid.LeftEdge
+ gre = self.grid.RightEdge
+ self.li = ((leaf_l_corner-gle)/dds).astype('int32')
+ self.ri = ((leaf_r_corner-gle)/dds).astype('int32')
+ self.dims = (self.ri - self.li).astype('int32')
+ # Here the cost is actually inversely proportional to 4**Level (empirical)
+ self.cost = (np.prod(self.dims)/4.**self.grid.Level).astype('int64')
+ # Here is the old way
+ # self.cost = np.prod(self.dims).astype('int64')
+ del dds, gle, gre
+
+ class dividing_node(node):
+ def __init__(self, split_ax, split_pos, left_children, right_children):
+ self.split_ax = split_ax
+ self.split_pos = split_pos
+ self.left_children = left_children
+ self.right_children = right_children
+
+ def count_cost(self,node):
+ if isinstance(node,AMRKDTree.leafnode):
+ return node.cost
+ else:
+ return self.count_cost(node.left_children) + self.count_cost(node.right_children)
+
+ def __build(self, grids, parent, l_corner, r_corner):
+ if len(grids) == 0:
+ self.leaf_count += 1
+ return AMRKDTree.leafnode(self.leaf_count-1, parent, l_corner, r_corner)
+
+ if len(grids) == 1:
+ thisgrid = grids[0]
+ if ( (thisgrid.LeftEdge[0] <= l_corner[0]) and (thisgrid.RightEdge[0] >= r_corner[0]) and
+ (thisgrid.LeftEdge[1] <= l_corner[1]) and (thisgrid.RightEdge[1] >= r_corner[1]) and
+ (thisgrid.LeftEdge[2] <= l_corner[2]) and (thisgrid.RightEdge[2] >= r_corner[2]) ):
+ children = []
+ if (len(thisgrid.Children) > 0) and (thisgrid.Level < self.L_MAX):
+ children = np.unique(thisgrid._get_child_index_mask()[thisgrid._get_child_index_mask() >= 0])
+ children = self.pf.hierarchy.grids[children -1]
+ child_l_data = np.array([child.LeftEdge for child in children])
+ child_r_data = np.array([child.RightEdge for child in children])
+ children_we_want = (child_l_data[:,0] < r_corner[0])*(child_r_data [:,0] > l_corner[0])* \
+ (child_l_data[:,1] < r_corner[1])*(child_r_data [:,1] > l_corner[1])* \
+ (child_l_data[:,2] < r_corner[2])*(child_r_data [:,2] > l_corner[2])
+ children = children[children_we_want]
+ return self.__build(children, thisgrid, l_corner, r_corner)
+
+ l_data = np.array([child.LeftEdge for child in grids])
+ r_data = np.array([child.RightEdge for child in grids])
+
+ best_dim = 0
+ best_choices = []
+ for d in range(3):
+ choices = np.unique(np.clip(np.concatenate(([l_corner[d]],l_data[:,d],r_data[:,d],[r_corner[d]])),l_corner[d],r_corner[d]))
+ if len(choices) > len(best_choices):
+ best_choices = choices
+ best_dim = d
+ best_choices.sort()
+ split = best_choices[len(best_choices)/2]
+
+ less_ids = np.nonzero(l_data[:,best_dim] < split)
+ greater_ids = np.nonzero(split < r_data[:,best_dim])
+ return AMRKDTree.dividing_node(d, split,
+ self.__build(grids[less_ids], parent,
+ l_corner, self.__corner_bounds(best_dim, split, current_right=r_corner)),
+ self.__build(grids[greater_ids], parent,
+ self.__corner_bounds(best_dim, split, current_left=l_corner), r_corner))
+
+ def __corner_bounds(self, split_dim, split, current_left = None, current_right = None):
+ if(current_left is not None):
+ new_left = deepcopy(current_left)
+ new_left[split_dim] = split
+ return new_left
+ elif(current_right is not None):
+ new_right = deepcopy(current_right)
+ new_right[split_dim] = split
+ return new_right
diff -r 0f1c190e5e47 -r 857c1f4280f5 yt/extensions/setup.py
--- a/yt/extensions/setup.py Tue May 04 07:51:43 2010 -0700
+++ b/yt/extensions/setup.py Wed May 05 08:57:31 2010 -0600
@@ -9,5 +9,6 @@
config.add_subpackage("lightcone")
config.add_subpackage("volume_rendering")
config.add_subpackage("kdtree")
+ config.add_subpackage("amr_kdtree")
config.add_subpackage("image_panner")
return config
diff -r 0f1c190e5e47 -r 857c1f4280f5 yt/extensions/volume_rendering/software_sampler.py
--- a/yt/extensions/volume_rendering/software_sampler.py Tue May 04 07:51:43 2010 -0700
+++ b/yt/extensions/volume_rendering/software_sampler.py Wed May 05 08:57:31 2010 -0600
@@ -28,6 +28,10 @@
from yt.extensions.volume_rendering import *
from yt.funcs import *
from yt.lagos import data_object_registry, ParallelAnalysisInterface
+from yt.extensions.amr_kdtree import *
+from yt.lagos.ParallelTools import *
+import time
+from copy import deepcopy
# We're going to register this class, but it does not directly inherit from
# AMRData.
@@ -37,7 +41,8 @@
resolution, transfer_function,
fields = None, whole_box = False,
sub_samples = 5, north_vector = None,
- pf = None):
+ pf = None, kd=False, l_max = None, kd_tree=None,
+ memory_factor=0.5, le=None, re=None):
# Now we replicate some of the 'cutting plane' logic
if not iterable(resolution):
resolution = (resolution, resolution)
@@ -50,6 +55,24 @@
if fields is None: fields = ["Density"]
self.fields = fields
self.transfer_function = transfer_function
+ self.pf = pf
+ self.l_max = l_max
+ self.total_cells = 0
+ self.total_cost = 0
+ self.time_in_vertex = 0.0
+ self.current_saved_grids = []
+ self.current_vcds = []
+ self.memory_per_process = 0.0
+ self.my_total = 0
+ self.my_storage = 0.
+ self.memory_factor = memory_factor
+ self.le = le
+ self.re = re
+ self.log_field=[]
+ for field in self.fields:
+ self.log_field.append(field in self.pf.field_info and
+ self.pf.field_info[field].take_log)
+
# Now we set up our various vectors
normal_vector /= na.sqrt( na.dot(normal_vector, normal_vector))
@@ -58,6 +81,8 @@
t = na.cross(normal_vector, vecs).sum(axis=1)
ax = t.argmax()
north_vector = na.cross(vecs[ax,:], normal_vector).ravel()
+ else:
+ north_vector = north_vector - na.dot(north_vector,normal_vector)*normal_vector
north_vector /= na.sqrt(na.dot(north_vector, north_vector))
east_vector = -na.cross(north_vector, normal_vector).ravel()
east_vector /= na.sqrt(na.dot(east_vector, east_vector))
@@ -72,9 +97,24 @@
self.back_center = center - 0.5*width[0]*self.unit_vectors[2]
self.front_center = center + 0.5*width[0]*self.unit_vectors[2]
- self._initialize_source()
+ if kd: self._kd_initialize_source(l_max=self.l_max, kd_tree=kd_tree)
+ else: self._initialize_source()
+
self._construct_vector_array()
+ def _kd_initialize_source(self,l_max=None, kd_tree=None):
+ self._kd_brick_tree = None
+ if self.pf is not None:
+ if kd_tree is None:
+ self._kd_brick_tree = AMRKDTree(self.pf,L_MAX=l_max, le=self.le, re=self.re)
+ else:
+ self._kd_brick_tree = kd_tree
+ self.total_cost = self._kd_brick_tree.total_cost
+ source = self.pf.hierarchy.inclined_box(self.origin, self.box_vectors)
+ self.source = source
+ self._base_source = source
+ self.res_fac = (1.0,1.0)
+
def _initialize_source(self):
check, source, rf = self._partition_hierarchy_2d_inclined(
self.unit_vectors, self.origin, self.width, self.box_vectors)
@@ -91,6 +131,77 @@
# _distributed can't be overridden in that case.
self._brick_collection = HomogenizedBrickCollection(self.source)
+ def cast_tree(self,tfp, node,viewpoint,pbar=None,up_every=10000):
+ if isinstance(node,AMRKDTree.leafnode):
+ my_rank = self._mpi_get_rank()
+ nprocs = self._mpi_get_size()
+ if ((self.total_cells >= 1.0*my_rank*self.total_cost/nprocs) and
+ (self.total_cells < 1.0*(my_rank+1)*self.total_cost/nprocs)):
+ #node_time = time.time()
+ if node.grid in self.current_saved_grids:
+ dds = self.current_vcds[self.current_saved_grids.index(node.grid)]
+ else:
+ #print 'Casting leaf %d, grid %s' % (node.leaf_id, node.grid)
+ t1 = time.time()
+ dds = []
+ for i,field in enumerate(self.fields):
+ vcd = node.grid.get_vertex_centered_data(field).astype('float64')
+ if self.log_field[i]: vcd = na.log10(vcd)
+ dds.append(vcd)
+ self.time_in_vertex += time.time()-t1
+ self.current_saved_grids.append(node.grid)
+ self.current_vcds.append(dds)
+ self.my_storage += na.prod(node.grid.ActiveDimensions)*len(self.fields)*8
+
+ while( self.my_storage*len(self.fields)*8 >= self.memory_per_process*self.memory_factor ):
+ self.my_storage -= len(self.fields)*8*na.prod(self.current_saved_grids[0].ActiveDimensions)
+ del self.current_saved_grids[0]
+ del self.curren_vcds[0]
+
+ data = [d[node.li[0]:node.ri[0]+1,
+ node.li[1]:node.ri[1]+1,
+ node.li[2]:node.ri[2]+1].copy() for d in dds]
+
+ b = PartitionedGrid(node.grid.id, len(self.fields), data,
+ node.leaf_l_corner.copy(),
+ node.leaf_r_corner.copy(),
+ node.dims.astype('int64'))
+ b.cast_plane(tfp, self.vector_plane)
+ #if my_rank == 0:
+ # print '[%04i] node_level: %d node_cost: %e node_time: %e'%(my_rank, node.grid.Level, node.cost, time.time()-node_time)
+ del data, b
+ self.my_total += node.cost
+ self.total_cells += node.cost
+ if pbar is not None:
+ pbar.update(self.total_cells)
+ else:
+ if viewpoint[node.split_ax] <= node.split_pos:
+ self.cast_tree(tfp,node.right_children,viewpoint,pbar=pbar)
+ self.cast_tree(tfp,node.left_children,viewpoint,pbar=pbar)
+ else:
+ self.cast_tree(tfp,node.left_children,viewpoint,pbar=pbar)
+ self.cast_tree(tfp,node.right_children,viewpoint,pbar=pbar)
+
+ def kd_ray_cast(self, finalize=False, pbar=None,up_every=100, memory_per_process=2**28):
+ self.memory_per_process = memory_per_process
+ if self._kd_brick_tree is None:
+ print 'No KD Tree Exists'
+ return
+ tfp = TransferFunctionProxy(self.transfer_function)
+ tfp.ns = self.sub_samples
+ self.total_cells = 0
+ pbar = get_pbar("Ray casting ", self._kd_brick_tree.total_cost)
More information about the yt-svn
mailing list