[Yt-svn] yt: 3 new changesets
hg at spacepope.org
hg at spacepope.org
Fri Apr 9 15:57:11 PDT 2010
hg Repository: yt
details: yt/rev/b3360e1df793
changeset: 1542:b3360e1df793
user: Stephen Skory <stephenskory at yahoo.com>
date:
Thu Apr 08 12:51:05 2010 -0700
description:
The structure function now passes the unit vector between
the points.
hg Repository: yt
details: yt/rev/ce36cb611817
changeset: 1543:ce36cb611817
user: Stephen Skory <stephenskory at yahoo.com>
date:
Fri Apr 09 15:56:23 2010 -0700
description:
Major update to the SFG. Vectorized many things for a ~23x
speedup. Requires rebuilding the fortran kd tree object.
hg Repository: yt
details: yt/rev/ed60ff36a8d3
changeset: 1544:ed60ff36a8d3
user: Stephen Skory <stephenskory at yahoo.com>
date:
Fri Apr 09 15:56:46 2010 -0700
description:
Merging.
diffstat:
scripts/yt_lodgeit.py | 1 +
yt/extensions/kdtree/fKD.f90 | 27 +++
yt/extensions/kdtree/fKD.v | 2 +
yt/lagos/BaseDataTypes.py | 2 +-
yt/lagos/DerivedQuantities.py | 18 +-
yt/lagos/EnzoFields.py | 2 +
yt/lagos/HaloFinding.py | 3 +-
yt/lagos/StructureFunctionGenerator.py | 318 +++++++++++++++++------------------
yt/raven/PlotCollection.py | 32 +++-
yt/raven/PlotTypes.py | 27 ++-
10 files changed, 255 insertions(+), 177 deletions(-)
diffs (truncated from 733 to 300 lines):
diff -r 6c5ab74488f2 -r ed60ff36a8d3 scripts/yt_lodgeit.py
--- a/scripts/yt_lodgeit.py Wed Apr 07 10:58:04 2010 -0700
+++ b/scripts/yt_lodgeit.py Fri Apr 09 15:56:46 2010 -0700
@@ -222,6 +222,7 @@
data = read_file(sys.stdin)
if not langopt:
mime = get_mimetype(data, '') or ''
+ fname = ""
elif len(filenames) == 1:
fname = filenames[0]
data = read_file(open(filenames[0], 'rb'))
diff -r 6c5ab74488f2 -r ed60ff36a8d3 yt/extensions/kdtree/fKD.f90
--- a/yt/extensions/kdtree/fKD.f90 Wed Apr 07 10:58:04 2010 -0700
+++ b/yt/extensions/kdtree/fKD.f90 Fri Apr 09 15:56:46 2010 -0700
@@ -43,6 +43,33 @@
end subroutine find_nn_nearest_neighbors
+subroutine find_many_nn_nearest_neighbors()
+ ! Given an input array of query vectors (qv_many), find their
+ ! nearest neighbors.
+ use kdtree2_module
+ use fKD_module
+ use kdtree2module
+ use tree_nodemodule
+ use intervalmodule
+
+ integer :: k, number
+ type(kdtree2_result),allocatable :: results(:)
+
+ allocate(results(nn))
+
+ number = size(qv_many,2)
+
+ do k=1, number
+ qv(:) = qv_many(:,k)
+ call kdtree2_n_nearest(tp=tree2,qv=qv,nn=nn,results=results)
+ nn_tags(:, k) = results%idx
+ end do
+
+ deallocate(results)
+ return
+
+end subroutine find_many_nn_nearest_neighbors
+
subroutine find_all_nn_nearest_neighbors()
! for all particles in pos, find their nearest neighbors and return the
! indexes and distances as big arrays
diff -r 6c5ab74488f2 -r ed60ff36a8d3 yt/extensions/kdtree/fKD.v
--- a/yt/extensions/kdtree/fKD.v Wed Apr 07 10:58:04 2010 -0700
+++ b/yt/extensions/kdtree/fKD.v Fri Apr 09 15:56:46 2010 -0700
@@ -12,6 +12,7 @@
dens(:) _real
mass(:) _real
qv(3) real
+qv_many(3,:) _real
nparts integer
nn integer
nMerge integer # number of nearest neighbors used in chain merging
@@ -97,5 +98,6 @@
create_tree() subroutine
free_tree() subroutine
find_all_nn_nearest_neighbors subroutine
+find_many_nn_nearest_neighbors subroutine
find_chunk_nearest_neighbors subroutine
chainHOP_tags_dens subroutine
diff -r 6c5ab74488f2 -r ed60ff36a8d3 yt/lagos/BaseDataTypes.py
--- a/yt/lagos/BaseDataTypes.py Wed Apr 07 10:58:04 2010 -0700
+++ b/yt/lagos/BaseDataTypes.py Fri Apr 09 15:56:46 2010 -0700
@@ -2747,7 +2747,7 @@
elif level == 0 and self.level == 0:
DLE = self.pf["DomainLeftEdge"]
self.global_startindex = na.array(na.floor(LL/ dx), dtype='int64')
- idims = na.ceil((self.right_edge-self.left_edge)/dx)
+ idims = na.rint((self.right_edge-self.left_edge)/dx).astype('int64')
self[field] = na.zeros(idims,dtype='float64')-999
def _refine(self, dlevel, field):
diff -r 6c5ab74488f2 -r ed60ff36a8d3 yt/lagos/DerivedQuantities.py
--- a/yt/lagos/DerivedQuantities.py Wed Apr 07 10:58:04 2010 -0700
+++ b/yt/lagos/DerivedQuantities.py Fri Apr 09 15:56:46 2010 -0700
@@ -372,8 +372,8 @@
p1 = p.sum()
if na.any(na.isnan(p)): raise ValueError
return p1 * (length_scale_factor / (mass_scale_factor**2.0))
-
-def _Extrema(data, fields, filter=None):
+
+def _Extrema(data, fields, non_zero = False, filter=None):
"""
This function returns the extrema of a set of fields
@@ -389,12 +389,18 @@
maxs.append(-1e90)
continue
if filter is None:
- mins.append(data[field].min())
- maxs.append(data[field].max())
+ # Note that we're hijacking an argument here
+ if non_zero: filter = data[field]>0.0
+ mins.append(data[field][filter].min())
+ maxs.append(data[field][filter].max())
else:
if this_filter.any():
- mins.append(data[field][this_filter].min())
- maxs.append(data[field][this_filter].max())
+ if non_zero:
+ nz_filter = ((this_filter) &
+ (data[field][this_filter] > 0.0))
+ else: nz_filter = this_filter
+ mins.append(data[field][nz_filter].min())
+ maxs.append(data[field][nz_filter].max())
else:
mins.append(1e90)
maxs.append(-1e90)
diff -r 6c5ab74488f2 -r ed60ff36a8d3 yt/lagos/EnzoFields.py
--- a/yt/lagos/EnzoFields.py Wed Apr 07 10:58:04 2010 -0700
+++ b/yt/lagos/EnzoFields.py Fri Apr 09 15:56:46 2010 -0700
@@ -138,11 +138,13 @@
def _TotalEnergy(field, data):
return data["Total_Energy"] / _convertEnergy(data)
add_field("TotalEnergy", function=_TotalEnergy,
+ display_name = "\mathrm{Total}\/\mathrm{Energy}",
units=r"\rm{ergs}/\rm{g}", convert_function=_convertEnergy)
def _Total_Energy(field, data):
return data["TotalEnergy"] / _convertEnergy(data)
add_field("Total_Energy", function=_Total_Energy,
+ display_name = "\mathrm{Total}\/\mathrm{Energy}",
units=r"\rm{ergs}/\rm{g}", convert_function=_convertEnergy)
def _NumberDensity(field, data):
diff -r 6c5ab74488f2 -r ed60ff36a8d3 yt/lagos/HaloFinding.py
--- a/yt/lagos/HaloFinding.py Wed Apr 07 10:58:04 2010 -0700
+++ b/yt/lagos/HaloFinding.py Fri Apr 09 15:56:46 2010 -0700
@@ -1269,8 +1269,7 @@
ms = -self.Tot_M.copy()
del self.Tot_M
Cx = self.CoM[:,0].copy()
- indexes = na.arange(self.group_count)
- sorted = na.asarray([indexes[i] for i in na.lexsort([indexes, Cx, ms])])
+ sorted = na.lexsort([Cx, ms])
del indexes, Cx, ms
self._groups = self._groups[sorted]
self._max_dens = self._max_dens[sorted]
diff -r 6c5ab74488f2 -r ed60ff36a8d3 yt/lagos/StructureFunctionGenerator.py
--- a/yt/lagos/StructureFunctionGenerator.py Wed Apr 07 10:58:04 2010 -0700
+++ b/yt/lagos/StructureFunctionGenerator.py Fri Apr 09 15:56:46 2010 -0700
@@ -25,16 +25,17 @@
from yt.lagos import *
from yt.math_utils import *
+from yt.performance_counters import yt_counters, time_function
try:
from yt.extensions.kdtree import *
except ImportError:
mylog.info("The Fortran kD-Tree did not import correctly. The structure function generator will not work correctly.")
-import math, sys, itertools, inspect, types, random
+import math, sys, itertools, inspect, types
from collections import defaultdict
class StructFcnGen(ParallelAnalysisInterface):
- def __init__(self, pf, left_edge=None, right_edge=None,
+ def __init__(self, pf, fields, left_edge=None, right_edge=None,
total_values=1000000, comm_size=10000, length_type="lin",
length_number=10, length_range=None, vol_ratio = 1,
salt=0):
@@ -63,7 +64,7 @@
number of point pairs per ruler length). Default: 0.
"""
self._fsets = []
- self.fields = set([])
+ self.fields = fields
# MPI stuff.
self.size = self._mpi_get_size()
self.mine = self._mpi_get_rank()
@@ -125,20 +126,16 @@
right_edge + self.lengths[-1], rank_ratio=self.vol_ratio)
mylog.info("LE %s RE %s %s" % (str(self.LE), str(self.RE), str(self.ds)))
self.width = self.ds.right_edge - self.ds.left_edge
- random.seed(-1 * self.mine + salt)
+ self.mt = na.random.mtrand.RandomState(seed = 1234 * self.mine + salt)
- def add_function(self, function, fields, out_labels, sqrt):
+ def add_function(self, function, out_labels, sqrt):
"""
Add a structure function to the list that will be evaluated at the
generated pairs of points.
"""
fargs = inspect.getargspec(function)
- if len(fargs.args) != 4:
- raise SyntaxError("The structure function %s needs four arguments." %\
- function.__name__)
- fields = list(fields)
- if len(fields) < 1:
- raise SyntaxError("Please specify at least one field for function %s." %\
+ if len(fargs.args) != 5:
+ raise SyntaxError("The structure function %s needs five arguments." %\
function.__name__)
out_labels = list(out_labels)
if len(out_labels) < 1:
@@ -148,9 +145,8 @@
if len(sqrt) != len(out_labels):
raise SyntaxError("Please have the same number of elements in out_labels as in sqrt for function %s." %\
function.__name__)
- self._fsets.append(StructSet(self, function, self.min_edge, fields,
+ self._fsets.append(StructSet(self, function, self.min_edge,
out_labels, sqrt))
- self.fields.update(fields)
return self._fsets[-1]
def __getitem__(self, key):
@@ -226,6 +222,9 @@
if self.mine == 0:
mylog.info("Length (%d of %d) %1.5e took %d communication cycles to complete." % \
(bigloop+1, len(self.lengths), length, self.comm_cycle_count))
+ if self.nlevels > 1:
+ del fKD.pos, fKD.qv_many, fKD.nn_tags
+ free_tree() # Frees the kdtree object.
self._allsum_bin_hits()
def _init_kd_tree(self):
@@ -241,12 +240,9 @@
fKD.pos[0, :] = xp[:]
fKD.pos[1, :] = yp[:]
fKD.pos[2, :] = zp[:]
- fKD.qv = na.empty(3, dtype='float64')
fKD.nn = 1
fKD.sort = False
fKD.rearrange = True
- fKD.tags = na.empty(1, dtype='int64')
- fKD.dist = na.empty(1, dtype='float64')
create_tree()
def _build_sort_array(self):
@@ -259,9 +255,8 @@
xp = self.ds["x"]
yp = self.ds["y"]
zp = self.ds["z"]
- self.sizes = [na.unique(xp).size, na.unique(yp).size, na.unique(zp).size]
- indices = na.arange(xp.size)
- self.sort = na.asarray([indices[i] for i in na.lexsort([indices, zp, yp, xp])])
+ self.sizes = [na.unique(xp).size, na.unique(yp).size, na.unique(zp).size]
+ self.sort = na.lexsort([zp, yp, xp])
del xp, yp, zp
self.ds.clear_data()
@@ -269,7 +264,7 @@
"""
Builds an array to store the field values array.
"""
- self.fields_vals = na.empty((self.comm_size, len(self.fields)), \
+ self.fields_vals = na.empty((self.comm_size, len(self.fields)*2), \
dtype='float64')
# At the same time build a dict to label the columns.
self.fields_columns = {}
@@ -345,7 +340,7 @@
to the left-hand neighbor.
"""
self.recv_points = na.ones((self.comm_size, 6), dtype='float64') * -1.
- self.recv_fields_vals = na.zeros((self.comm_size, len(self.fields)), \
+ self.recv_fields_vals = na.zeros((self.comm_size, len(self.fields)*2), \
dtype='float64')
self.recv_gen_array = na.zeros(self.size, dtype='int64')
self.recv_hooks.append(self._mpi_Irecv_double(self.recv_points, \
@@ -393,140 +388,140 @@
fset.length_bin_hits[length] = \
fset.length_bin_hits[length].reshape(fset.bin_number)
- def _pick_random_points(self, length):
+ def _pick_random_points(self, length, size):
"""
- Picks out two random points separated by *length*.
+ Picks out size random pairs separated by length *length*.
"""
- # A random point inside this subvolume.
- r1 = na.empty(3, dtype='float64')
- for i in xrange(3):
- r1[i] = random.uniform(self.LE[i], self.RE[i]) # 3 randoms.
- # Pick our theta, phi random pair.
- theta = random.uniform(0, 2.*math.pi) # 1 random.
- phi = random.uniform(-math.pi/2., math.pi/2.) # 1 random.
- # Find our second point.
- r2 = na.array([r1[0] + length * math.cos(theta) * math.cos(phi),
- r1[1] + length * math.sin(theta) * math.cos(phi),
- r1[2] + length * math.sin(phi)], dtype='float64')
- # Reflect it so it's inside the (full) volume.
- r2 = r2 % self.period
- return (r1, na.array(r2))
+ # First make random points inside this subvolume.
+ r1 = na.empty((size,3), dtype='float64')
+ for dim in range(3):
+ r1[:,dim] = self.mt.uniform(low=self.ds.left_edge[dim],
+ high=self.ds.right_edge[dim], size=size)
+ # Next we find the second point, determined by a random
+ # theta, phi angle.
+ theta = self.mt.uniform(low=0, high=2.*math.pi, size=size)
+ phi = self.mt.uniform(low=-math.pi/2., high=math.pi/2., size=size)
+ r2 = na.empty((size,3), dtype='float64')
More information about the yt-svn
mailing list