[yt-svn] commit/yt: 9 new changesets
Bitbucket
commits-noreply at bitbucket.org
Mon Sep 10 08:24:05 PDT 2012
9 new commits in yt:
https://bitbucket.org/yt_analysis/yt/changeset/971c4d171969/
changeset: 971c4d171969
branch: yt
user: xarthisius
date: 2012-09-05 12:48:46
summary: Add module for Cartesian domain decomposition
affected #: 1 file
diff -r 0c844cb7df71949206ff497c57ef437d6429d58e -r 971c4d171969e388fdf4f10eb288ae9b67023e3e yt/utilities/decompose.py
--- /dev/null
+++ b/yt/utilities/decompose.py
@@ -0,0 +1,157 @@
+"""
+Automagical cartesian domain decomposition.
+
+Author: Kacper Kowalik <xarthisius.kk at gmail.com>
+Affiliation: CA UMK
+Author: Artur Gawryszczak <gawrysz at gmail.com>
+Affiliation: PCSS
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2012 Kacper Kowalik. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+SIEVE_PRIMES = \
+ lambda l: l and l[:1] + SIEVE_PRIMES([n for n in l if n % l[0]])
+
+
+def decompose_to_primes(max_prime):
+ """ Decompose number into the primes """
+ for prime in SIEVE_PRIMES(range(2, max_prime)):
+ if prime * prime > max_prime:
+ break
+ while max_prime % prime == 0:
+ yield prime
+ max_prime /= prime
+ if max_prime > 1:
+ yield max_prime
+
+
+def decompose_array(arr, psize, bbox):
+ """ Calculate list of product(psize) subarrays of arr, along with their
+ left and right edges
+ """
+ grid_left_edges = np.empty([np.product(psize), 3], dtype=np.float64)
+ grid_right_edges = np.empty([np.product(psize), 3], dtype=np.float64)
+ n_d = arr.shape
+ d_s = (bbox[:, 1] - bbox[:, 0]) / n_d
+ dist = np.mgrid[bbox[0, 0]:bbox[0, 1]:d_s[0],
+ bbox[1, 0]:bbox[1, 1]:d_s[1],
+ bbox[2, 0]:bbox[2, 1]:d_s[2]]
+ for i in range(3):
+ xyz = split_array(dist[i], psize)
+ for j in range(np.product(psize)):
+ grid_left_edges[j, i] = xyz[j][0, 0, 0]
+ grid_right_edges[j, i] = xyz[j][-1, -1, -1] + d_s[i]
+ del xyz
+ del dist
+ patches = split_array(arr, psize)
+ return grid_left_edges, grid_right_edges, patches
+
+
+def evaluate_domain_decomposition(n_d, pieces, ldom):
+ """ Evaluate longest to shortest edge ratio
+ BEWARE: lot's of magic here """
+ ideal_bsize = 3.0 * (pieces * np.product(n_d) ** 2) ** (1.0 / 3.0)
+ bsize = int(np.sum(
+ ldom / np.array(n_d, dtype=np.float64) * np.product(n_d)))
+ load_balance = float(np.product(n_d)) / \
+ (float(pieces) * np.product((n_d - 1) / ldom + 1))
+
+ # 0.25 is magic number
+ quality = load_balance / (1 + 0.25 * (bsize / ideal_bsize - 1.0))
+ # \todo add a factor that estimates lower cost when x-direction is
+ # not chopped too much
+ # \deprecated estimate these magic numbers
+ quality *= (1. - (0.001 * ldom[0] + 0.0001 * ldom[1]) / pieces)
+ if np.any(ldom > n_d):
+ quality = 0
+
+ return quality
+
+
+def factorize_number(pieces):
+ """ Return array consiting of prime, its power and number of different
+ decompositions in three dimensions for this prime
+ """
+ factors = [factor for factor in decompose_to_primes(pieces)]
+ temp = np.bincount(factors)
+ return np.array(
+ [(prime, temp[prime], (temp[prime] + 1) * (temp[prime] + 2) / 2)
+ for prime in np.unique(factors)]
+ )
+
+
+def get_psize(n_d, pieces):
+ """ Calculate the best division of array into px*py*pz subarrays.
+ The goal is to minimize the ratio of longest to shortest edge
+ to minimize the amount of inter-process communication.
+ """
+ fac = factorize_number(pieces)
+ print fac, pieces
+ nfactors = len(fac[:, 2])
+ best = 0.0
+ while np.all(fac[:, 2] > 0):
+ ldom = np.ones(3, dtype=np.int)
+ for nfac in range(nfactors):
+ i = int(np.sqrt(0.25 + 2 * (fac[nfac, 2] - 1)) - 0.5)
+ k = fac[nfac, 2] - int(1 + i * (i + 1) / 2)
+ i = fac[nfac, 1] - i
+ j = fac[nfac, 1] - (i + k)
+ ldom *= fac[nfac, 0] ** np.array([i, j, k])
+
+ quality = evaluate_domain_decomposition(n_d, pieces, ldom)
+ if quality > best:
+ best = quality
+ p_size = ldom
+ # search for next unique combination
+ for j in range(nfactors):
+ if fac[j, 2] > 1:
+ fac[j, 2] -= 1
+ break
+ else:
+ if (j < nfactors - 1):
+ fac[j, 2] = int((fac[j, 1] + 1) * (fac[j, 1] + 2) / 2)
+ else:
+ fac[:, 2] = 0 # no more combinations to try
+
+ return p_size
+
+
+def split_array(tab, psize):
+ """ Split array into px*py*pz subarrays using internal numpy routine. """
+ temp = [np.array_split(array, psize[1], axis=1)
+ for array in np.array_split(tab, psize[2], axis=2)]
+ temp = [item for sublist in temp for item in sublist]
+ temp = [np.array_split(array, psize[0], axis=0) for array in temp]
+ temp = [item for sublist in temp for item in sublist]
+ return temp
+
+
+if __name__ == "__main__":
+
+ NPROC = 12
+ ARRAY = np.zeros((128, 128, 129))
+ BBOX = np.array([[0., 1.0], [-1.5, 1.5], [1.0, 2.5]])
+
+ PROCS = get_psize(np.array(ARRAY.shape), NPROC)
+ LE, RE, DATA = decompose_array(ARRAY, PROCS, BBOX)
+
+ for idx in range(NPROC):
+ print LE[idx, :], RE[idx, :], DATA[idx].shape
https://bitbucket.org/yt_analysis/yt/changeset/f4e3776f28ab/
changeset: f4e3776f28ab
branch: yt
user: xarthisius
date: 2012-09-05 12:49:41
summary: stream::load_uniform_data can handle arrays of any size and split them into smaller chunks for parallelization
affected #: 1 file
diff -r 971c4d171969e388fdf4f10eb288ae9b67023e3e -r f4e3776f28ab22fa25c426be566555a59b57bc4d yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -40,6 +40,8 @@
FieldInfoContainer, NullFunc
from yt.utilities.lib import \
get_box_grids_level
+from yt.utilities.decompose import \
+ decompose_array, get_psize
from yt.utilities.definitions import \
mpc_conversion, sec_conversion
@@ -152,7 +154,6 @@
self.pf.field_info.add_field(
field, lambda a, b: None,
convert_function=cf, take_log=False)
-
def _parse_hierarchy(self):
self.grid_dimensions = self.stream_handler.dimensions
@@ -296,13 +297,13 @@
@property
def all_fields(self): return self[0].keys()
-def load_uniform_grid(data, domain_dimensions, domain_size_in_cm,
- sim_time=0.0, number_of_particles=0):
+def load_uniform_grid(data, domain_dimensions, bbox, bbox_to_cm=1.0,
+ nprocs=1, sim_time=0.0, number_of_particles=0):
r"""Load a uniform grid of data into yt as a
:class:`~yt.frontends.stream.data_structures.StreamHandler`.
This should allow a uniform grid of data to be loaded directly into yt and
- analyzed as would any others. This comes with several caveats:
+ nalyzed as would any others. This comes with several caveats:
* Units will be incorrect unless the data has already been converted to
cgs.
* Some functions may behave oddly, and parallelism will be
@@ -313,55 +314,74 @@
----------
data : dict
This is a dict of numpy arrays, where the keys are the field names.
- domain_dimensiosn : array_like
+ domain_dimensions : array_like
This is the domain dimensions of the grid
- domain_size_in_cm : float
- The size of the domain, in centimeters
+ bbox : array_like (xdim:zdim, LE:RE)
+ Size of computational domain in centimeters
+ bbox_to_cm : float, optional
+ Conversion factor from bbox units to centimeters
+ nprocs: integer, optional
+ If greater than 1, will create this number of subarrays out of data
sim_time : float, optional
The simulation time in seconds
number_of_particles : int, optional
If particle fields are included, set this to the number of particles
-
+
Examples
--------
- >>> arr = na.random.random((256, 256, 256))
+ >>> arr = na.random.random((128, 128, 129))
>>> data = dict(Density = arr)
- >>> pf = load_uniform_grid(data, [256, 256, 256], 3.08e24)
-
+ >>> bbox = np.array([[0., 1.0], [-1.5, 1.5], [1.0, 2.5]])
+ >>> pf = load_uniform_grid(data, arr.shape, bbox, nprocs=12)
+
"""
+
+ domain_dimensions = na.array(domain_dimensions)
+ domain_left_edge = na.array(bbox[:, 0], 'float64')
+ domain_right_edge = na.array(bbox[:, 1], 'float64')
+ grid_levels = na.zeros(nprocs, dtype='int32').reshape((nprocs,1))
+
sfh = StreamDictFieldHandler()
- sfh.update({0:data})
- domain_dimensions = na.array(domain_dimensions)
- if na.unique(domain_dimensions).size != 1:
- print "We don't support variably sized domains yet."
- raise RuntimeError
- domain_left_edge = na.zeros(3, 'float64')
- domain_right_edge = na.ones(3, 'float64')
- grid_left_edges = na.zeros(3, "int64").reshape((1,3))
- grid_right_edges = na.array(domain_dimensions, "int64").reshape((1,3))
- grid_levels = na.array([0], dtype='int32').reshape((1,1))
+ if nprocs > 1:
+ temp = {}
+ new_data = {}
+ for key in data.keys():
+ psize = get_psize(na.array(data[key].shape), nprocs)
+ grid_left_edges, grid_right_edges, temp[key] = \
+ decompose_array(data[key], psize, bbox)
+ for gid in range(nprocs):
+ new_data[gid] = {}
+ for key in temp.keys():
+ new_data[gid].update({key:temp[key][gid]})
+ sfh.update(new_data)
+ del new_data, temp
+ else:
+ sfh.update({0:data})
+ grid_left_edges = na.zeros(3, "int64").reshape((1,3))
+ grid_right_edges = na.array(domain_dimensions, "int64").reshape((1,3))
+
+ grid_left_edges = grid_left_edges.astype("float64")
+ grid_left_edges /= domain_dimensions*2**grid_levels
+ grid_left_edges *= domain_right_edge - domain_left_edge
+ grid_left_edges += domain_left_edge
+
+ grid_right_edges = grid_right_edges.astype("float64")
+ grid_right_edges /= domain_dimensions*2**grid_levels
+ grid_right_edges *= domain_right_edge - domain_left_edge
+ grid_right_edges += domain_left_edge
+
grid_dimensions = grid_right_edges - grid_left_edges
- grid_left_edges = grid_left_edges.astype("float64")
- grid_left_edges /= domain_dimensions*2**grid_levels
- grid_left_edges *= domain_right_edge - domain_left_edge
- grid_left_edges += domain_left_edge
-
- grid_right_edges = grid_right_edges.astype("float64")
- grid_right_edges /= domain_dimensions*2**grid_levels
- grid_right_edges *= domain_right_edge - domain_left_edge
- grid_right_edges += domain_left_edge
-
handler = StreamHandler(
grid_left_edges,
grid_right_edges,
grid_dimensions,
grid_levels,
na.array([-1], dtype='int64'),
- number_of_particles*na.ones(1, dtype='int64').reshape((1,1)),
- na.zeros(1).reshape((1,1)),
+ number_of_particles*na.ones(nprocs, dtype='int64').reshape(nprocs,1),
+ na.zeros(nprocs).reshape((nprocs,1)),
sfh,
)
@@ -375,10 +395,10 @@
handler.cosmology_simulation = 0
spf = StreamStaticOutput(handler)
- spf.units["cm"] = domain_size_in_cm
+ spf.units["cm"] = bbox_to_cm
spf.units['1'] = 1.0
spf.units["unitary"] = 1.0
- box_in_mpc = domain_size_in_cm / mpc_conversion['cm']
+ box_in_mpc = bbox_to_cm / mpc_conversion['cm']
for unit in mpc_conversion.keys():
spf.units[unit] = mpc_conversion[unit] * box_in_mpc
return spf
https://bitbucket.org/yt_analysis/yt/changeset/d00b3ecbda4f/
changeset: d00b3ecbda4f
branch: yt
user: xarthisius
date: 2012-09-05 13:20:42
summary: Fix grid_dimensions and parent_ids
affected #: 1 file
diff -r f4e3776f28ab22fa25c426be566555a59b57bc4d -r d00b3ecbda4f13448a7a729e816451771c391429 yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -332,7 +332,7 @@
>>> arr = na.random.random((128, 128, 129))
>>> data = dict(Density = arr)
- >>> bbox = np.array([[0., 1.0], [-1.5, 1.5], [1.0, 2.5]])
+ >>> bbox = na.array([[0., 1.0], [-1.5, 1.5], [1.0, 2.5]])
>>> pf = load_uniform_grid(data, arr.shape, bbox, nprocs=12)
"""
@@ -351,6 +351,7 @@
psize = get_psize(na.array(data[key].shape), nprocs)
grid_left_edges, grid_right_edges, temp[key] = \
decompose_array(data[key], psize, bbox)
+ grid_dimensions = na.array([grid.shape for grid in temp[key]])
for gid in range(nprocs):
new_data[gid] = {}
for key in temp.keys():
@@ -371,15 +372,15 @@
grid_right_edges /= domain_dimensions*2**grid_levels
grid_right_edges *= domain_right_edge - domain_left_edge
grid_right_edges += domain_left_edge
+ grid_dimensions = grid_right_edges - grid_left_edges
- grid_dimensions = grid_right_edges - grid_left_edges
handler = StreamHandler(
grid_left_edges,
grid_right_edges,
grid_dimensions,
grid_levels,
- na.array([-1], dtype='int64'),
+ -na.ones(nprocs, dtype='int64'),
number_of_particles*na.ones(nprocs, dtype='int64').reshape(nprocs,1),
na.zeros(nprocs).reshape((nprocs,1)),
sfh,
https://bitbucket.org/yt_analysis/yt/changeset/a8600ea68f49/
changeset: a8600ea68f49
branch: yt
user: xarthisius
date: 2012-09-05 14:28:01
summary: [gdf:writer] grid_parent_id should be array, even if filled with dummy number for now
affected #: 1 file
diff -r d00b3ecbda4f13448a7a729e816451771c391429 -r a8600ea68f497baf16074f05caaea4bff11f8969 yt/utilities/grid_data_format/writer.py
--- a/yt/utilities/grid_data_format/writer.py
+++ b/yt/utilities/grid_data_format/writer.py
@@ -139,7 +139,7 @@
f["grid_left_index"] = pf.h.grid_left_edge
f["grid_level"] = pf.h.grid_levels
# @todo: Do we need to loop over the grids for this?
- f["grid_parent_id"] = -1
+ f["grid_parent_id"] = -np.ones(pf.h.grid_dimensions.shape[0])
f["grid_particle_count"] = pf.h.grid_particle_count
###
https://bitbucket.org/yt_analysis/yt/changeset/2cb1493b5aa4/
changeset: 2cb1493b5aa4
branch: yt
user: xarthisius
date: 2012-09-05 15:03:46
summary: remove debug print
affected #: 1 file
diff -r a8600ea68f497baf16074f05caaea4bff11f8969 -r 2cb1493b5aa44adf5e69f9cdb48b607183558f89 yt/utilities/decompose.py
--- a/yt/utilities/decompose.py
+++ b/yt/utilities/decompose.py
@@ -104,7 +104,6 @@
to minimize the amount of inter-process communication.
"""
fac = factorize_number(pieces)
- print fac, pieces
nfactors = len(fac[:, 2])
best = 0.0
while np.all(fac[:, 2] > 0):
https://bitbucket.org/yt_analysis/yt/changeset/40139dff7b92/
changeset: 40139dff7b92
branch: yt
user: xarthisius
date: 2012-09-05 15:05:50
summary: [gdf:writer] add dummy variables to field required by yt reader, fix grid_left_index that was erroneously assigned with grid_left_edge
affected #: 1 file
diff -r 2cb1493b5aa44adf5e69f9cdb48b607183558f89 -r 40139dff7b92fa81cb7fce008618810395f8c0cf yt/utilities/grid_data_format/writer.py
--- a/yt/utilities/grid_data_format/writer.py
+++ b/yt/utilities/grid_data_format/writer.py
@@ -83,11 +83,11 @@
g.attrs["unique_identifier"] = pf.unique_identifier
g.attrs["cosmological_simulation"] = pf.cosmological_simulation
# @todo: Where is this in the yt API?
- #g.attrs["num_ghost_zones"] = pf...
+ g.attrs["num_ghost_zones"] = 0
# @todo: Where is this in the yt API?
- #g.attrs["field_ordering"] = pf...
+ g.attrs["field_ordering"] = 0
# @todo: not yet supported by yt.
- #g.attrs["boundary_conditions"] = pf...
+ g.attrs["boundary_conditions"] = np.array([0, 0, 0, 0, 0, 0], 'int32')
if pf.cosmological_simulation:
g.attrs["current_redshift"] = pf.current_redshift
@@ -136,9 +136,11 @@
# root datasets -- info about the grids
###
f["grid_dimensions"] = pf.h.grid_dimensions
- f["grid_left_index"] = pf.h.grid_left_edge
+ f["grid_left_index"] = np.array(
+ [g.get_global_startindex() for g in pf.h.grids]
+ ).reshape(pf.h.grid_dimensions.shape[0], 3)
f["grid_level"] = pf.h.grid_levels
- # @todo: Do we need to loop over the grids for this?
+ # @todo: Fill with proper values
f["grid_parent_id"] = -np.ones(pf.h.grid_dimensions.shape[0])
f["grid_particle_count"] = pf.h.grid_particle_count
https://bitbucket.org/yt_analysis/yt/changeset/cdc68c51044c/
changeset: cdc68c51044c
branch: yt
user: xarthisius
date: 2012-09-05 15:07:54
summary: [stream:load_uniform_grid] simplify and retain original API with small change: use unit conversion factor rather than demand domain size in cm
affected #: 1 file
diff -r 40139dff7b92fa81cb7fce008618810395f8c0cf -r cdc68c51044c73a9f681a2d2b4c90fba21536d1a yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -297,13 +297,13 @@
@property
def all_fields(self): return self[0].keys()
-def load_uniform_grid(data, domain_dimensions, bbox, bbox_to_cm=1.0,
+def load_uniform_grid(data, domain_dimensions, sim_unit_to_cm, bbox=None,
nprocs=1, sim_time=0.0, number_of_particles=0):
r"""Load a uniform grid of data into yt as a
:class:`~yt.frontends.stream.data_structures.StreamHandler`.
This should allow a uniform grid of data to be loaded directly into yt and
- nalyzed as would any others. This comes with several caveats:
+ analyzed as would any others. This comes with several caveats:
* Units will be incorrect unless the data has already been converted to
cgs.
* Some functions may behave oddly, and parallelism will be
@@ -316,10 +316,11 @@
This is a dict of numpy arrays, where the keys are the field names.
domain_dimensions : array_like
This is the domain dimensions of the grid
- bbox : array_like (xdim:zdim, LE:RE)
- Size of computational domain in centimeters
- bbox_to_cm : float, optional
- Conversion factor from bbox units to centimeters
+ sim_unit_to_cm : float
+ Conversion factor from simulation units to centimeters
+ bbox : array_like (xdim:zdim, LE:RE), optional
+ Size of computational domain in units sim_unit_to_cm, if the latter
+ is not provided centimeters are assumed
nprocs: integer, optional
If greater than 1, will create this number of subarrays out of data
sim_time : float, optional
@@ -333,11 +334,13 @@
>>> arr = na.random.random((128, 128, 129))
>>> data = dict(Density = arr)
>>> bbox = na.array([[0., 1.0], [-1.5, 1.5], [1.0, 2.5]])
- >>> pf = load_uniform_grid(data, arr.shape, bbox, nprocs=12)
+ >>> pf = load_uniform_grid(data, arr.shape, 3.08e24, bbox=bbox, nprocs=12)
"""
domain_dimensions = na.array(domain_dimensions)
+ if bbox is None:
+ bbox = na.array([[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]], 'float64')
domain_left_edge = na.array(bbox[:, 0], 'float64')
domain_right_edge = na.array(bbox[:, 1], 'float64')
grid_levels = na.zeros(nprocs, dtype='int32').reshape((nprocs,1))
@@ -360,20 +363,9 @@
del new_data, temp
else:
sfh.update({0:data})
- grid_left_edges = na.zeros(3, "int64").reshape((1,3))
- grid_right_edges = na.array(domain_dimensions, "int64").reshape((1,3))
-
- grid_left_edges = grid_left_edges.astype("float64")
- grid_left_edges /= domain_dimensions*2**grid_levels
- grid_left_edges *= domain_right_edge - domain_left_edge
- grid_left_edges += domain_left_edge
-
- grid_right_edges = grid_right_edges.astype("float64")
- grid_right_edges /= domain_dimensions*2**grid_levels
- grid_right_edges *= domain_right_edge - domain_left_edge
- grid_right_edges += domain_left_edge
- grid_dimensions = grid_right_edges - grid_left_edges
-
+ grid_left_edges = domain_left_edge
+ grid_right_edges = domain_right_edge
+ grid_dimensions = domain_dimensions.reshape(nprocs,3)
handler = StreamHandler(
grid_left_edges,
@@ -396,10 +388,10 @@
handler.cosmology_simulation = 0
spf = StreamStaticOutput(handler)
- spf.units["cm"] = bbox_to_cm
+ spf.units["cm"] = sim_unit_to_cm
spf.units['1'] = 1.0
spf.units["unitary"] = 1.0
- box_in_mpc = bbox_to_cm / mpc_conversion['cm']
+ box_in_mpc = sim_unit_to_cm / mpc_conversion['cm']
for unit in mpc_conversion.keys():
spf.units[unit] = mpc_conversion[unit] * box_in_mpc
return spf
https://bitbucket.org/yt_analysis/yt/changeset/387113039a96/
changeset: 387113039a96
branch: yt
user: xarthisius
date: 2012-09-05 15:17:14
summary: update docstring
affected #: 1 file
diff -r cdc68c51044c73a9f681a2d2b4c90fba21536d1a -r 387113039a96cb27f8081240801d2de161be04a0 yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -319,8 +319,7 @@
sim_unit_to_cm : float
Conversion factor from simulation units to centimeters
bbox : array_like (xdim:zdim, LE:RE), optional
- Size of computational domain in units sim_unit_to_cm, if the latter
- is not provided centimeters are assumed
+ Size of computational domain in units sim_unit_to_cm
nprocs: integer, optional
If greater than 1, will create this number of subarrays out of data
sim_time : float, optional
https://bitbucket.org/yt_analysis/yt/changeset/64c090a6d57f/
changeset: 64c090a6d57f
branch: yt
user: MatthewTurk
date: 2012-09-10 17:24:04
summary: Merged in xarthisius/yt (pull request #263)
affected #: 3 files
diff -r 5b0682498d834a77ec3779928f7d2d2611554c7f -r 64c090a6d57f30cd7dfa87fe242bb7128f4f791a yt/frontends/stream/data_structures.py
--- a/yt/frontends/stream/data_structures.py
+++ b/yt/frontends/stream/data_structures.py
@@ -40,6 +40,8 @@
FieldInfoContainer, NullFunc
from yt.utilities.lib import \
get_box_grids_level
+from yt.utilities.decompose import \
+ decompose_array, get_psize
from yt.utilities.definitions import \
mpc_conversion, sec_conversion
@@ -152,7 +154,6 @@
self.pf.field_info.add_field(
field, lambda a, b: None,
convert_function=cf, take_log=False)
-
def _parse_hierarchy(self):
self.grid_dimensions = self.stream_handler.dimensions
@@ -296,8 +297,8 @@
@property
def all_fields(self): return self[0].keys()
-def load_uniform_grid(data, domain_dimensions, domain_size_in_cm,
- sim_time=0.0, number_of_particles=0):
+def load_uniform_grid(data, domain_dimensions, sim_unit_to_cm, bbox=None,
+ nprocs=1, sim_time=0.0, number_of_particles=0):
r"""Load a uniform grid of data into yt as a
:class:`~yt.frontends.stream.data_structures.StreamHandler`.
@@ -313,55 +314,66 @@
----------
data : dict
This is a dict of numpy arrays, where the keys are the field names.
- domain_dimensiosn : array_like
+ domain_dimensions : array_like
This is the domain dimensions of the grid
- domain_size_in_cm : float
- The size of the domain, in centimeters
+ sim_unit_to_cm : float
+ Conversion factor from simulation units to centimeters
+ bbox : array_like (xdim:zdim, LE:RE), optional
+ Size of computational domain in units sim_unit_to_cm
+ nprocs: integer, optional
+ If greater than 1, will create this number of subarrays out of data
sim_time : float, optional
The simulation time in seconds
number_of_particles : int, optional
If particle fields are included, set this to the number of particles
-
+
Examples
--------
- >>> arr = na.random.random((256, 256, 256))
+ >>> arr = na.random.random((128, 128, 129))
>>> data = dict(Density = arr)
- >>> pf = load_uniform_grid(data, [256, 256, 256], 3.08e24)
-
+ >>> bbox = na.array([[0., 1.0], [-1.5, 1.5], [1.0, 2.5]])
+ >>> pf = load_uniform_grid(data, arr.shape, 3.08e24, bbox=bbox, nprocs=12)
+
"""
+
+ domain_dimensions = na.array(domain_dimensions)
+ if bbox is None:
+ bbox = na.array([[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]], 'float64')
+ domain_left_edge = na.array(bbox[:, 0], 'float64')
+ domain_right_edge = na.array(bbox[:, 1], 'float64')
+ grid_levels = na.zeros(nprocs, dtype='int32').reshape((nprocs,1))
+
sfh = StreamDictFieldHandler()
- sfh.update({0:data})
- domain_dimensions = na.array(domain_dimensions)
- if na.unique(domain_dimensions).size != 1:
- print "We don't support variably sized domains yet."
- raise RuntimeError
- domain_left_edge = na.zeros(3, 'float64')
- domain_right_edge = na.ones(3, 'float64')
- grid_left_edges = na.zeros(3, "int64").reshape((1,3))
- grid_right_edges = na.array(domain_dimensions, "int64").reshape((1,3))
- grid_levels = na.array([0], dtype='int32').reshape((1,1))
- grid_dimensions = grid_right_edges - grid_left_edges
-
- grid_left_edges = grid_left_edges.astype("float64")
- grid_left_edges /= domain_dimensions*2**grid_levels
- grid_left_edges *= domain_right_edge - domain_left_edge
- grid_left_edges += domain_left_edge
-
- grid_right_edges = grid_right_edges.astype("float64")
- grid_right_edges /= domain_dimensions*2**grid_levels
- grid_right_edges *= domain_right_edge - domain_left_edge
- grid_right_edges += domain_left_edge
+ if nprocs > 1:
+ temp = {}
+ new_data = {}
+ for key in data.keys():
+ psize = get_psize(na.array(data[key].shape), nprocs)
+ grid_left_edges, grid_right_edges, temp[key] = \
+ decompose_array(data[key], psize, bbox)
+ grid_dimensions = na.array([grid.shape for grid in temp[key]])
+ for gid in range(nprocs):
+ new_data[gid] = {}
+ for key in temp.keys():
+ new_data[gid].update({key:temp[key][gid]})
+ sfh.update(new_data)
+ del new_data, temp
+ else:
+ sfh.update({0:data})
+ grid_left_edges = domain_left_edge
+ grid_right_edges = domain_right_edge
+ grid_dimensions = domain_dimensions.reshape(nprocs,3)
handler = StreamHandler(
grid_left_edges,
grid_right_edges,
grid_dimensions,
grid_levels,
- na.array([-1], dtype='int64'),
- number_of_particles*na.ones(1, dtype='int64').reshape((1,1)),
- na.zeros(1).reshape((1,1)),
+ -na.ones(nprocs, dtype='int64'),
+ number_of_particles*na.ones(nprocs, dtype='int64').reshape(nprocs,1),
+ na.zeros(nprocs).reshape((nprocs,1)),
sfh,
)
@@ -375,10 +387,10 @@
handler.cosmology_simulation = 0
spf = StreamStaticOutput(handler)
- spf.units["cm"] = domain_size_in_cm
+ spf.units["cm"] = sim_unit_to_cm
spf.units['1'] = 1.0
spf.units["unitary"] = 1.0
- box_in_mpc = domain_size_in_cm / mpc_conversion['cm']
+ box_in_mpc = sim_unit_to_cm / mpc_conversion['cm']
for unit in mpc_conversion.keys():
spf.units[unit] = mpc_conversion[unit] * box_in_mpc
return spf
diff -r 5b0682498d834a77ec3779928f7d2d2611554c7f -r 64c090a6d57f30cd7dfa87fe242bb7128f4f791a yt/utilities/decompose.py
--- /dev/null
+++ b/yt/utilities/decompose.py
@@ -0,0 +1,156 @@
+"""
+Automagical cartesian domain decomposition.
+
+Author: Kacper Kowalik <xarthisius.kk at gmail.com>
+Affiliation: CA UMK
+Author: Artur Gawryszczak <gawrysz at gmail.com>
+Affiliation: PCSS
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2012 Kacper Kowalik. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+
+SIEVE_PRIMES = \
+ lambda l: l and l[:1] + SIEVE_PRIMES([n for n in l if n % l[0]])
+
+
+def decompose_to_primes(max_prime):
+ """ Decompose number into the primes """
+ for prime in SIEVE_PRIMES(range(2, max_prime)):
+ if prime * prime > max_prime:
+ break
+ while max_prime % prime == 0:
+ yield prime
+ max_prime /= prime
+ if max_prime > 1:
+ yield max_prime
+
+
+def decompose_array(arr, psize, bbox):
+ """ Calculate list of product(psize) subarrays of arr, along with their
+ left and right edges
+ """
+ grid_left_edges = np.empty([np.product(psize), 3], dtype=np.float64)
+ grid_right_edges = np.empty([np.product(psize), 3], dtype=np.float64)
+ n_d = arr.shape
+ d_s = (bbox[:, 1] - bbox[:, 0]) / n_d
+ dist = np.mgrid[bbox[0, 0]:bbox[0, 1]:d_s[0],
+ bbox[1, 0]:bbox[1, 1]:d_s[1],
+ bbox[2, 0]:bbox[2, 1]:d_s[2]]
+ for i in range(3):
+ xyz = split_array(dist[i], psize)
+ for j in range(np.product(psize)):
+ grid_left_edges[j, i] = xyz[j][0, 0, 0]
+ grid_right_edges[j, i] = xyz[j][-1, -1, -1] + d_s[i]
+ del xyz
+ del dist
+ patches = split_array(arr, psize)
+ return grid_left_edges, grid_right_edges, patches
+
+
+def evaluate_domain_decomposition(n_d, pieces, ldom):
+ """ Evaluate longest to shortest edge ratio
+ BEWARE: lot's of magic here """
+ ideal_bsize = 3.0 * (pieces * np.product(n_d) ** 2) ** (1.0 / 3.0)
+ bsize = int(np.sum(
+ ldom / np.array(n_d, dtype=np.float64) * np.product(n_d)))
+ load_balance = float(np.product(n_d)) / \
+ (float(pieces) * np.product((n_d - 1) / ldom + 1))
+
+ # 0.25 is magic number
+ quality = load_balance / (1 + 0.25 * (bsize / ideal_bsize - 1.0))
+ # \todo add a factor that estimates lower cost when x-direction is
+ # not chopped too much
+ # \deprecated estimate these magic numbers
+ quality *= (1. - (0.001 * ldom[0] + 0.0001 * ldom[1]) / pieces)
+ if np.any(ldom > n_d):
+ quality = 0
+
+ return quality
+
+
+def factorize_number(pieces):
+ """ Return array consiting of prime, its power and number of different
+ decompositions in three dimensions for this prime
+ """
+ factors = [factor for factor in decompose_to_primes(pieces)]
+ temp = np.bincount(factors)
+ return np.array(
+ [(prime, temp[prime], (temp[prime] + 1) * (temp[prime] + 2) / 2)
+ for prime in np.unique(factors)]
+ )
+
+
+def get_psize(n_d, pieces):
+ """ Calculate the best division of array into px*py*pz subarrays.
+ The goal is to minimize the ratio of longest to shortest edge
+ to minimize the amount of inter-process communication.
+ """
+ fac = factorize_number(pieces)
+ nfactors = len(fac[:, 2])
+ best = 0.0
+ while np.all(fac[:, 2] > 0):
+ ldom = np.ones(3, dtype=np.int)
+ for nfac in range(nfactors):
+ i = int(np.sqrt(0.25 + 2 * (fac[nfac, 2] - 1)) - 0.5)
+ k = fac[nfac, 2] - int(1 + i * (i + 1) / 2)
+ i = fac[nfac, 1] - i
+ j = fac[nfac, 1] - (i + k)
+ ldom *= fac[nfac, 0] ** np.array([i, j, k])
+
+ quality = evaluate_domain_decomposition(n_d, pieces, ldom)
+ if quality > best:
+ best = quality
+ p_size = ldom
+ # search for next unique combination
+ for j in range(nfactors):
+ if fac[j, 2] > 1:
+ fac[j, 2] -= 1
+ break
+ else:
+ if (j < nfactors - 1):
+ fac[j, 2] = int((fac[j, 1] + 1) * (fac[j, 1] + 2) / 2)
+ else:
+ fac[:, 2] = 0 # no more combinations to try
+
+ return p_size
+
+
+def split_array(tab, psize):
+ """ Split array into px*py*pz subarrays using internal numpy routine. """
+ temp = [np.array_split(array, psize[1], axis=1)
+ for array in np.array_split(tab, psize[2], axis=2)]
+ temp = [item for sublist in temp for item in sublist]
+ temp = [np.array_split(array, psize[0], axis=0) for array in temp]
+ temp = [item for sublist in temp for item in sublist]
+ return temp
+
+
+if __name__ == "__main__":
+
+ NPROC = 12
+ ARRAY = np.zeros((128, 128, 129))
+ BBOX = np.array([[0., 1.0], [-1.5, 1.5], [1.0, 2.5]])
+
+ PROCS = get_psize(np.array(ARRAY.shape), NPROC)
+ LE, RE, DATA = decompose_array(ARRAY, PROCS, BBOX)
+
+ for idx in range(NPROC):
+ print LE[idx, :], RE[idx, :], DATA[idx].shape
diff -r 5b0682498d834a77ec3779928f7d2d2611554c7f -r 64c090a6d57f30cd7dfa87fe242bb7128f4f791a yt/utilities/grid_data_format/writer.py
--- a/yt/utilities/grid_data_format/writer.py
+++ b/yt/utilities/grid_data_format/writer.py
@@ -83,11 +83,11 @@
g.attrs["unique_identifier"] = pf.unique_identifier
g.attrs["cosmological_simulation"] = pf.cosmological_simulation
# @todo: Where is this in the yt API?
- #g.attrs["num_ghost_zones"] = pf...
+ g.attrs["num_ghost_zones"] = 0
# @todo: Where is this in the yt API?
- #g.attrs["field_ordering"] = pf...
+ g.attrs["field_ordering"] = 0
# @todo: not yet supported by yt.
- #g.attrs["boundary_conditions"] = pf...
+ g.attrs["boundary_conditions"] = np.array([0, 0, 0, 0, 0, 0], 'int32')
if pf.cosmological_simulation:
g.attrs["current_redshift"] = pf.current_redshift
@@ -136,10 +136,12 @@
# root datasets -- info about the grids
###
f["grid_dimensions"] = pf.h.grid_dimensions
- f["grid_left_index"] = pf.h.grid_left_edge
+ f["grid_left_index"] = np.array(
+ [g.get_global_startindex() for g in pf.h.grids]
+ ).reshape(pf.h.grid_dimensions.shape[0], 3)
f["grid_level"] = pf.h.grid_levels
- # @todo: Do we need to loop over the grids for this?
- f["grid_parent_id"] = -1
+ # @todo: Fill with proper values
+ f["grid_parent_id"] = -np.ones(pf.h.grid_dimensions.shape[0])
f["grid_particle_count"] = pf.h.grid_particle_count
###
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