[Yt-svn] commit/yt: 10 new changesets
Bitbucket
commits-noreply at bitbucket.org
Thu Jun 9 10:55:24 PDT 2011
10 new changesets in yt:
http://bitbucket.org/yt_analysis/yt/changeset/30f62efd8734/
changeset: 30f62efd8734
branch: yt
user: MatthewTurk
date: 2011-05-07 00:18:15
summary: Adding a very simple initial wrapper of voro++
affected #: 2 files (2.4 KB)
--- a/yt/utilities/setup.py Thu May 05 10:02:16 2011 -0700
+++ b/yt/utilities/setup.py Fri May 06 15:18:15 2011 -0700
@@ -172,6 +172,10 @@
glob.glob("yt/utilities/_amr_utils/*.h") +
glob.glob("yt/utilities/_amr_utils/*.c"),
)
+ #config.add_extension("voropp",
+ # ["yt/utilities/voropp.pyx"],
+ # language="c++",
+ # include_dirs=["yt/utilities/voro++"])
config.add_extension("libconfig_wrapper",
["yt/utilities/libconfig_wrapper.pyx"] +
glob.glob("yt/utilities/_libconfig/*.c"),
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/voropp.pyx Fri May 06 15:18:15 2011 -0700
@@ -0,0 +1,62 @@
+"""
+Wrapping code for voro++
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2011 Matthew Turk. 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/>.
+"""
+
+from cython.operator cimport dereference as deref, preincrement as inc
+from libc.stdlib cimport malloc, free, abs, calloc, labs
+cimport libcpp
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+cdef extern from "voro++.cc":
+ cdef cppclass container:
+ container(double xmin, double xmax, double ymin, double ymax,
+ double zmin, double zmax, int nx, int ny, int nz,
+ libcpp.bool xper, libcpp.bool yper, libcpp.bool zper, int alloc)
+ void put(int n, double x, double y, double z)
+ void store_cell_volumes(double *vols)
+
+cdef class VoronoiVolume:
+ cdef container *my_con
+ cdef int npart
+ def __init__(self, xi, yi, zi):
+ self.my_con = new container(0.0, 1.0, 0.0, 1.0, 0.0, 1.0,
+ xi, yi, zi, False, False, False, 8)
+ self.npart = 0
+
+ def add_array(self, np.ndarray[np.float64_t, ndim=1] xpos,
+ np.ndarray[np.float64_t, ndim=1] ypos,
+ np.ndarray[np.float64_t, ndim=1] zpos):
+ cdef int i
+ for i in range(xpos.shape[0]):
+ self.my_con.put(self.npart, xpos[i], ypos[i], zpos[i])
+ self.npart += 1
+
+ def get_volumes(self):
+ cdef np.ndarray vol = np.zeros(self.npart, 'double')
+ cdef double *vdouble = <double *> vol.data
+ self.my_con.store_cell_volumes(vdouble)
+ return vol
http://bitbucket.org/yt_analysis/yt/changeset/d1c226913621/
changeset: d1c226913621
branch: yt
user: MatthewTurk
date: 2011-05-07 00:35:54
summary: Adding cython decorators and adding a deallocator
affected #: 1 file (174 bytes)
--- a/yt/utilities/voropp.pyx Fri May 06 15:18:15 2011 -0700
+++ b/yt/utilities/voropp.pyx Fri May 06 15:35:54 2011 -0700
@@ -47,6 +47,11 @@
xi, yi, zi, False, False, False, 8)
self.npart = 0
+ def __dealloc__(self):
+ del self.my_con
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
def add_array(self, np.ndarray[np.float64_t, ndim=1] xpos,
np.ndarray[np.float64_t, ndim=1] ypos,
np.ndarray[np.float64_t, ndim=1] zpos):
@@ -55,6 +60,8 @@
self.my_con.put(self.npart, xpos[i], ypos[i], zpos[i])
self.npart += 1
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
def get_volumes(self):
cdef np.ndarray vol = np.zeros(self.npart, 'double')
cdef double *vdouble = <double *> vol.data
http://bitbucket.org/yt_analysis/yt/changeset/5f7577488be7/
changeset: 5f7577488be7
branch: yt
user: MatthewTurk
date: 2011-05-12 17:39:30
summary: Adding a disjoint set structure with a union/find implementation.
affected #: 4 files (4.1 KB)
--- a/yt/utilities/_amr_utils/ContourFinding.pyx Wed May 11 11:30:17 2011 -0400
+++ b/yt/utilities/_amr_utils/ContourFinding.pyx Thu May 12 11:39:30 2011 -0400
@@ -38,6 +38,16 @@
if i0 < i1: return i0
return i1
+cdef extern from "union_find.h":
+ ctypedef struct forest_node:
+ void *value
+ forest_node *parent
+ int rank
+
+ forest_node* MakeSet(void* value)
+ void Union(forest_node* node1, forest_node* node2)
+ forest_node* Find(forest_node* node)
+
#@cython.boundscheck(False)
#@cython.wraparound(False)
def construct_boundary_relationships(
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/union_find.c Thu May 12 11:39:30 2011 -0400
@@ -0,0 +1,66 @@
+/* Copyright (c) 2011 the authors listed at the following URL, and/or
+the authors of referenced articles or incorporated external code:
+http://en.literateprograms.org/Disjoint_set_data_structure_(C)?action=history&offset=20080516180553
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+Retrieved from: http://en.literateprograms.org/Disjoint_set_data_structure_(C)?oldid=13366
+*/
+
+#include <stdlib.h>
+
+#include "union_find.h"
+
+forest_node* MakeSet(void* value) {
+ forest_node* node = malloc(sizeof(forest_node));
+ node->value = value;
+ node->parent = NULL;
+ node->rank = 0;
+ return node;
+}
+
+void Union(forest_node* node1, forest_node* node2) {
+ if (node1->rank > node2->rank) {
+ node2->parent = node1;
+ } else if (node2->rank > node1->rank) {
+ node1->parent = node2;
+ } else { /* they are equal */
+ node2->parent = node1;
+ node1->rank++;
+ }
+}
+
+forest_node* Find(forest_node* node) {
+ forest_node* temp;
+ /* Find the root */
+ forest_node* root = node;
+ while (root->parent != NULL) {
+ root = root->parent;
+ }
+ /* Update the parent pointers */
+ while (node->parent != NULL) {
+ temp = node->parent;
+ node->parent = root;
+ node = temp;
+ }
+ return root;
+}
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/union_find.h Thu May 12 11:39:30 2011 -0400
@@ -0,0 +1,41 @@
+/* Copyright (c) 2011 the authors listed at the following URL, and/or
+the authors of referenced articles or incorporated external code:
+http://en.literateprograms.org/Disjoint_set_data_structure_(C)?action=history&offset=20080516180553
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+Retrieved from: http://en.literateprograms.org/Disjoint_set_data_structure_(C)?oldid=13366
+*/
+
+#ifndef _UNION_FIND_H_
+#define _UNION_FIND_H_
+
+typedef struct forest_node_t {
+ void* value;
+ struct forest_node_t* parent;
+ int rank;
+} forest_node;
+
+forest_node* MakeSet(void* value);
+void Union(forest_node* node1, forest_node* node2);
+forest_node* Find(forest_node* node);
+
+#endif /* #ifndef _UNION_FIND_H_ */
+
--- a/yt/utilities/setup.py Wed May 11 11:30:17 2011 -0400
+++ b/yt/utilities/setup.py Thu May 12 11:39:30 2011 -0400
@@ -161,7 +161,8 @@
config.add_extension("amr_utils",
["yt/utilities/amr_utils.pyx",
"yt/utilities/_amr_utils/FixedInterpolator.c",
- "yt/utilities/_amr_utils/kdtree.c"] +
+ "yt/utilities/_amr_utils/kdtree.c",
+ "yt/utilities/_amr_utils/union_find.c"] +
glob.glob("yt/utilities/_amr_utils/healpix_*.c"),
define_macros=[("PNG_SETJMP_NOT_SUPPORTED", True)],
include_dirs=["yt/utilities/_amr_utils/", png_inc,
http://bitbucket.org/yt_analysis/yt/changeset/2b61626339e0/
changeset: 2b61626339e0
branch: yt
user: MatthewTurk
date: 2011-05-13 15:30:39
summary: Removing a bunch of old, unused code, and starting to add a join tree object
for grids.
affected #: 2 files (9.0 KB)
--- a/yt/analysis_modules/level_sets/contour_finder.py Thu May 12 11:39:30 2011 -0400
+++ b/yt/analysis_modules/level_sets/contour_finder.py Fri May 13 09:30:39 2011 -0400
@@ -30,202 +30,6 @@
import yt.utilities.data_point_utilities as data_point_utilities
import yt.utilities.amr_utils as amr_utils
-class GridConsiderationQueue:
- def __init__(self, white_list, priority_func=None):
- """
- This class exists to serve the contour finder. It ensures that
- we can create a cascading set of queue dependencies, and if
- a grid is touched again ahead of time we can bump it to the top
- of the queue again. It like has few uses.
- """
- self.to_consider = []
- self.considered = set()
- self.n = 0
- self.white_list = set(white_list)
- self.priority_func = priority_func
-
- def add(self, grids, force=False):
- if not hasattr(grids,'size'):
- grids = ensure_list(grids)
- i = self.n
- to_check = self.white_list.intersection(grids)
- if not force: to_check.difference_update(self.considered)
- for g in sorted(to_check, key=self.priority_func):
- try:
- # We only delete from subsequent checks
- del self.to_consider[self.to_consider.index(g, i)]
- self.to_consider.insert(i,g)
- i += 1
- except ValueError:
- self.to_consider.append(g)
-
- def __iter__(self):
- return self
-
- def next(self):
- if self.n >= len(self.to_consider):
- raise StopIteration
- tr = self.to_consider[self.n]
- self.considered.add(tr)
- self.n += 1
- return tr
-
- def progress(self):
- return self.n, len(self.to_consider)
-
-# We want an algorithm that deals with growing a given contour to *all* the
-# cells in a grid.
-
-def old_identify_contours(data_source, field, min_val, max_val, cached_fields=None):
- """
- Given a *data_source*, we will search for topologically connected sets
- in *field* between *min_val* and *max_val*.
- """
- if cached_fields is None: cached_fields = defaultdict(lambda: dict())
- maxn_cells = na.sum([g.ActiveDimensions.prod() for g in data_source._grids])
- contour_ind = na.where( (data_source[field] > min_val)
- & (data_source[field] < max_val))[0]
- np = contour_ind.size
- if np == 0:
- return {}
- cur_max_id = maxn_cells - np
- contour_ids = na.arange(maxn_cells, cur_max_id, -1) + 1 # Minimum of 1
- data_source["tempContours"] = na.ones(data_source[field].shape, dtype='int32')*-1
- mylog.info("Contouring over %s cells with %s candidates", contour_ids[0],np)
- data_source["tempContours"][contour_ind] = contour_ids[:]
- data_source._flush_data_to_grids("tempContours", -1, dtype='int32')
- my_queue = GridConsiderationQueue(white_list = data_source._grids,
- priority_func = lambda g: -1*g["tempContours"].max())
- my_queue.add(data_source._grids)
- for i,grid in enumerate(my_queue):
- mylog.info("Examining %s of %s", *my_queue.progress())
- max_before = grid["tempContours"].max()
- to_get = ["tempContours"]
- if field in cached_fields[grid.id] and \
- not na.any( (cached_fields[grid.id][field] > min_val)
- & (cached_fields[grid.id][field] < max_val)):
- continue
- for f in [field, "GridIndices"]:
- if f not in cached_fields[grid.id]: to_get.append(f)
- cg = grid.retrieve_ghost_zones(1,to_get)
- for f in [field, "GridIndices"]:
- if f in cached_fields[grid.id]:
- cg.data[f] = cached_fields[grid.id][f]
- else:
- cached_fields[grid.id][f] = cg[f]
- local_ind = na.where( (cg[field] > min_val)
- & (cg[field] < max_val)
- & (cg["tempContours"] == -1) )
- if local_ind[0].size > 0:
- kk = na.arange(cur_max_id, cur_max_id-local_ind[0].size, -1)
- cg["tempContours"][local_ind] = kk[:]
- cur_max_id -= local_ind[0].size
- fd = cg["tempContours"].astype('int64')
- fd_original = fd.copy()
- if na.all(fd > -1):
- fd[:] = fd.max()
- else:
- xi_u,yi_u,zi_u = na.where(fd > -1)
- cor_order = na.argsort(-1*fd[(xi_u,yi_u,zi_u)])
- xi = xi_u[cor_order]
- yi = yi_u[cor_order]
- zi = zi_u[cor_order]
- while data_point_utilities.FindContours(fd, xi, yi, zi) < 0: pass
- cg["tempContours"] = fd.copy().astype('float64')
- cg.flush_data("tempContours")
- my_queue.add(cg._grids)
- force_ind = na.unique(cg["GridIndices"][na.where(
- (cg["tempContours"] > fd_original)
- & (cg["GridIndices"] != grid.id-grid._id_offset) )])
- if len(force_ind) > 0:
- my_queue.add(data_source.hierarchy.grids[force_ind.astype('int32')], force=True)
- for ax in 'xyz':
- if not iterable(grid['d%s'%ax]):
- grid['d%s'%ax] = grid['d%s'%ax]*na.ones(grid.ActiveDimensions)
- del data_source.data["tempContours"] # Force a reload from the grids
- data_source.get_data("tempContours", in_grids=True)
- i = 0
- contour_ind = {}
- for contour_id in na.unique(data_source["tempContours"]):
- if contour_id == -1: continue
- contour_ind[i] = na.where(data_source["tempContours"] == contour_id)
- mylog.debug("Contour id %s has %s cells", i, contour_ind[i][0].size)
- i += 1
- mylog.info("Identified %s contours between %0.5e and %0.5e",
- len(contour_ind.keys()),min_val,max_val)
- for grid in chain(data_source._grids, cg._grids):
- grid.data.pop("tempContours", None)
- del data_source.data["tempContours"]
- return contour_ind
-
-def check_neighbors(data_object, field="Contours"):
- """
- This method is a means of error checking in the contour finder.
- """
- n_bad = na.zeros(1, dtype='int32')
- for cid in na.unique(data_object[field]):
- if cid == -1: continue
- ids = na.where(data_object[field] == cid)
- mx = data_object['x'][ids].copy()
- my = data_object['y'][ids].copy()
- mz = data_object['z'][ids].copy()
- mdx = data_object['dx'][ids].copy()
- mdy = data_object['dy'][ids].copy()
- mdz = data_object['dz'][ids].copy()
- grid_ids_m = data_object['GridIndices'][ids].copy()
- grid_levels_m = data_object['GridLevel'][ids].copy()
- mp = mx.size
- ids = na.where( (data_object[field] != cid)
- & (data_object[field] >= 0 ))
- nx = data_object['x'][ids].copy()
- ny = data_object['y'][ids].copy()
- nz = data_object['z'][ids].copy()
- ndx = data_object['dx'][ids].copy()
- ndy = data_object['dy'][ids].copy()
- ndz = data_object['dz'][ids].copy()
- grid_ids_n = data_object['GridIndices'][ids].copy()
- grid_levels_n = data_object['GridLevel'][ids].copy()
- np = nx.size
- weave.inline(check_cell_distance,
- ['mx','my','mz','mdx','mdy','mdz','mp',
- 'nx','ny','nz','ndx','ndy','ndz','np','n_bad',
- 'grid_ids_m', 'grid_levels_m', 'grid_ids_n', 'grid_levels_n'],
- compiler='gcc', type_converters=converters.blitz,
- auto_downcast=0, verbose=2)
- return n_bad[0]
-
-check_cell_distance = \
-r"""
-using namespace std;
-int i, j, k;
-long double cell_dist, rad_m, rad_n;
-k=0;
-for(i=0;i<mp;i++){
- for(j=0;j<np;j++){
- /*
- cell_dist = sqrtl(pow(mx(i)-nx(j),2) +
- pow(my(i)-ny(j),2) +
- pow(mz(i)-nz(j),2));
- rad_m = sqrtl(pow(mdx(i),2) +
- pow(mdy(i),2) +
- pow(mdz(i),2));
- rad_n = sqrtl(pow(ndx(j),2) +
- pow(ndy(j),2) +
- pow(ndz(j),2));
- */
- //if(cell_dist > 1.01 * (rad_n/2.0+rad_m/2.0)) continue;
- if(fabsl(mx(i)-nx(j))>(mdx(i)+ndx(j))/2.0) continue;
- if(fabsl(my(i)-ny(j))>(mdy(i)+ndy(j))/2.0) continue;
- if(fabsl(mz(i)-nz(j))>(mdz(i)+ndz(j))/2.0) continue;
- k++;
- break;
- cout << cell_dist << "\t" << 1.01*(rad_n/2.0+rad_m/2.0) << "\t";
- cout << grid_ids_m(i) << "\t" << grid_ids_n(j) << endl;
- }
-}
-n_bad(0) += k;
-"""
-
def coalesce_join_tree(jtree1):
joins = defaultdict(set)
nj = jtree1.shape[0]
--- a/yt/utilities/_amr_utils/ContourFinding.pyx Thu May 12 11:39:30 2011 -0400
+++ b/yt/utilities/_amr_utils/ContourFinding.pyx Fri May 13 09:30:39 2011 -0400
@@ -26,6 +26,7 @@
import numpy as np
cimport numpy as np
cimport cython
+from stdlib cimport malloc, free
cdef extern from "math.h":
double fabs(double x)
@@ -48,6 +49,40 @@
void Union(forest_node* node1, forest_node* node2)
forest_node* Find(forest_node* node)
+ctypedef struct CellIdentifier:
+ np.int64_t hindex
+ int level
+
+cdef class GridContourContainer:
+ cdef np.int64_t dims[3]
+ cdef np.int64_t start_indices[3]
+ cdef forest_node **join_tree
+ cdef np.int64_t ncells
+
+ def __init__(self, dimensions, indices):
+ cdef int i
+ self.ncells = 1
+ for i in range(3):
+ self.ncells *= dimensions[i]
+ self.dims[i] = dimensions[i]
+ self.start_indices[i] = indices[i]
+ self.join_tree = <forest_node **> malloc(sizeof(forest_node)
+ * self.ncells)
+ for i in range(self.ncells): self.join_tree[i] = NULL
+
+ def __dealloc__(self):
+ cdef int i
+ for i in range(self.ncells):
+ if self.join_tree[i] != NULL: free(self.join_tree[i])
+ free(self.join_tree)
+
+ #def construct_join_tree(self,
+ # np.ndarray[np.float64_t, ndim=3] field,
+ # np.ndarray[np.bool_t, ndim=3] mask):
+ # # This only looks at the components of the grid that are actually
+ # # inside this grid -- boundary conditions are handled later.
+ # pass
+
#@cython.boundscheck(False)
#@cython.wraparound(False)
def construct_boundary_relationships(
http://bitbucket.org/yt_analysis/yt/changeset/67d76649fdcc/
changeset: 67d76649fdcc
branch: yt
user: MatthewTurk
date: 2011-05-18 17:22:51
summary: Merging the voro++ stuff with the experimentation with the new join tree.
affected #: 2 files (2.5 KB)
--- a/yt/utilities/setup.py Fri May 13 09:30:39 2011 -0400
+++ b/yt/utilities/setup.py Wed May 18 11:22:51 2011 -0400
@@ -173,6 +173,10 @@
glob.glob("yt/utilities/_amr_utils/*.h") +
glob.glob("yt/utilities/_amr_utils/*.c"),
)
+ #config.add_extension("voropp",
+ # ["yt/utilities/voropp.pyx"],
+ # language="c++",
+ # include_dirs=["yt/utilities/voro++"])
config.add_extension("libconfig_wrapper",
["yt/utilities/libconfig_wrapper.pyx"] +
glob.glob("yt/utilities/_libconfig/*.c"),
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/voropp.pyx Wed May 18 11:22:51 2011 -0400
@@ -0,0 +1,69 @@
+"""
+Wrapping code for voro++
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2011 Matthew Turk. 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/>.
+"""
+
+from cython.operator cimport dereference as deref, preincrement as inc
+from libc.stdlib cimport malloc, free, abs, calloc, labs
+cimport libcpp
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+cdef extern from "voro++.cc":
+ cdef cppclass container:
+ container(double xmin, double xmax, double ymin, double ymax,
+ double zmin, double zmax, int nx, int ny, int nz,
+ libcpp.bool xper, libcpp.bool yper, libcpp.bool zper, int alloc)
+ void put(int n, double x, double y, double z)
+ void store_cell_volumes(double *vols)
+
+cdef class VoronoiVolume:
+ cdef container *my_con
+ cdef int npart
+ def __init__(self, xi, yi, zi):
+ self.my_con = new container(0.0, 1.0, 0.0, 1.0, 0.0, 1.0,
+ xi, yi, zi, False, False, False, 8)
+ self.npart = 0
+
+ def __dealloc__(self):
+ del self.my_con
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def add_array(self, np.ndarray[np.float64_t, ndim=1] xpos,
+ np.ndarray[np.float64_t, ndim=1] ypos,
+ np.ndarray[np.float64_t, ndim=1] zpos):
+ cdef int i
+ for i in range(xpos.shape[0]):
+ self.my_con.put(self.npart, xpos[i], ypos[i], zpos[i])
+ self.npart += 1
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def get_volumes(self):
+ cdef np.ndarray vol = np.zeros(self.npart, 'double')
+ cdef double *vdouble = <double *> vol.data
+ self.my_con.store_cell_volumes(vdouble)
+ return vol
http://bitbucket.org/yt_analysis/yt/changeset/7fbc91c3db60/
changeset: 7fbc91c3db60
branch: yt
user: MatthewTurk
date: 2011-06-09 00:31:54
summary: Fixed resolution projections now contain a conditional for non-periodicity.
This results in a 2x speedup.
affected #: 2 files (1.5 KB)
--- a/yt/data_objects/data_containers.py Wed Jun 08 11:17:23 2011 -0700
+++ b/yt/data_objects/data_containers.py Wed Jun 08 15:31:54 2011 -0700
@@ -2139,6 +2139,9 @@
self._dls = {}
self.domain_width = na.rint((self.pf.domain_right_edge -
self.pf.domain_left_edge)/self.dds).astype('int64')
+ if not na.any((self.global_startindex + self.ActiveDimensions) >
+ self.domain_width):
+ self.domain_width[:] = 0
self._refresh_data()
def _get_list_of_grids(self):
@@ -2189,6 +2192,7 @@
for i,grid in enumerate(self._get_grids()):
mylog.debug("Getting fields from %s", i)
self._get_data_from_grid(grid, fields_to_get, dls)
+ mylog.info("IO completed; summing")
for field in fields_to_get:
self[field] = self._mpi_Allsum_double(self[field])
conv = self.pf.units[self.pf.field_info[field].projection_conversion]
--- a/yt/utilities/data_point_utilities.c Wed Jun 08 11:17:23 2011 -0700
+++ b/yt/utilities/data_point_utilities.c Wed Jun 08 15:31:54 2011 -0700
@@ -880,7 +880,7 @@
npy_int64 gxs, gys, gzs, gxe, gye, gze;
npy_int64 cxs, cys, czs, cxe, cye, cze;
npy_int64 ixs, iys, izs, ixe, iye, ize;
- npy_int64 gxi, gyi, gzi, cxi, cyi, czi;
+ int gxi, gyi, gzi, cxi, cyi, czi;
npy_int64 cdx, cdy, cdz;
npy_int64 dw[3];
int i;
@@ -1014,17 +1014,17 @@
ci = (cxi % dw[0]);
ci = (ci < 0) ? ci + dw[0] : ci;
if ( ci < gxs*refratio || ci >= gxe*refratio) continue;
- gxi = floor(ci / refratio) - gxs;
+ gxi = ((int) (ci / refratio)) - gxs;
for(cyi=cys;cyi<=cye;cyi++) {
cj = cyi % dw[1];
cj = (cj < 0) ? cj + dw[1] : cj;
if ( cj < gys*refratio || cj >= gye*refratio) continue;
- gyi = floor(cj / refratio) - gys;
+ gyi = ((int) (cj / refratio)) - gys;
for(czi=czs;czi<=cze;czi++) {
ck = czi % dw[2];
ck = (ck < 0) ? ck + dw[2] : ck;
if ( ck < gzs*refratio || ck >= gze*refratio) continue;
- gzi = floor(ck / refratio) - gzs;
+ gzi = ((int) (ck / refratio)) - gzs;
if ((ll) || (*(npy_int32*)PyArray_GETPTR3(mask, gxi,gyi,gzi) > 0))
{
for(n=0;n<n_fields;n++){
@@ -1214,43 +1214,75 @@
cye = (cys + cdy - 1);
cze = (czs + cdz - 1);
- /* It turns out that C89 doesn't define a mechanism for choosing the sign
- of the remainder.
- */
int x_loc, y_loc; // For access into the buffer
- for(cxi=cxs;cxi<=cxe;cxi++) {
- ci = (cxi % dw[0]);
- ci = (ci < 0) ? ci + dw[0] : ci;
- if ( ci < gxs*refratio || ci >= gxe*refratio) continue;
+
+ /* We check here if the domain is important or not.
+ If it's not, then, well, we get to use the fast version. */
+ if (dw[0] == dw[1] == dw[2] == 0) {
+ for(gxi=gxs,cxi=gxs*refratio;gxi<gxe;gxi++,cxi+=refratio) {
+ for(gyi=gys,cyi=gys*refratio;gyi<gye;gyi++,cyi+=refratio) {
+ for(gzi=gzs,czi=gzs*refratio;gzi<gze;gzi++,czi+=refratio) {
+ if ((refratio!=1) &&
+ (*(npy_int32*)PyArray_GETPTR3(mask, gxi,gyi,gzi)==0)) continue;
+ switch (axis) {
+ case 0: x_loc = cyi-cys; y_loc = czi-czs; break;
+ case 1: x_loc = cxi-cxs; y_loc = czi-czs; break;
+ case 2: x_loc = cxi-cys; y_loc = cyi-cys; break;
+ }
+ //fprintf(stderr, "%d %d %d %d %d\n", x_loc, y_loc, gxi, gyi, gzi);
+ for(ri=0;ri<refratio;ri++){
+ for(rj=0;rj<refratio;rj++){
+ for(n=0;n<n_fields;n++){
+ for(n=0;n<n_fields;n++){
+ *(npy_float64*) PyArray_GETPTR2(c_data[n], x_loc+ri, y_loc+rj)
+ += *(npy_float64*) PyArray_GETPTR3(g_data[n],
+ gxi-gxs, gyi-gys, gzi-gzs) * dls[n] / refratio;
+ }
+ }
+ }
+ }
+ total+=1;
+ }
+ }
+ }
+ } else {
+ /* Gotta go the slow route. */
+ for(cxi=gxs*refratio;cxi<=cxe;cxi++) {
+ /* It turns out that C89 doesn't define a mechanism for choosing the sign
+ of the remainder.
+ */
+ ci=cxi;//ci = (cxi % dw[0]);
+ //ci = (ci < 0) ? ci + dw[0] : ci;
+ if ( ci >= gxe*refratio) break;
gxi = floor(ci / refratio) - gxs;
- for(cyi=cys;cyi<=cye;cyi++) {
- cj = cyi % dw[1];
- cj = (cj < 0) ? cj + dw[1] : cj;
- if ( cj < gys*refratio || cj >= gye*refratio) continue;
- gyi = floor(cj / refratio) - gys;
- for(czi=czs;czi<=cze;czi++) {
- ck = czi % dw[2];
- ck = (ck < 0) ? ck + dw[2] : ck;
- if ( ck < gzs*refratio || ck >= gze*refratio) continue;
- gzi = floor(ck / refratio) - gzs;
- if (refratio == 1 || *(npy_int32*)PyArray_GETPTR3(mask, gxi,gyi,gzi) > 0)
- {
- switch (axis) {
- case 0: x_loc = cyi-cys; y_loc = czi-czs; break;
- case 1: x_loc = cxi-cxs; y_loc = czi-czs; break;
- case 2: x_loc = cxi-cys; y_loc = cyi-cys; break;
- }
- for(n=0;n<n_fields;n++){
- *(npy_float64*) PyArray_GETPTR2(c_data[n], x_loc, y_loc)
- += *(npy_float64*) PyArray_GETPTR3(g_data[n], gxi, gyi, gzi)
- * dls[n] / refratio;
- }
- total += 1;
- }
+ for(cyi=gys*refratio;cyi<=cye;cyi++) {
+ cj=cyi;//cj = cyi % dw[1];
+ //cj = (cj < 0) ? cj + dw[1] : cj;
+ if ( cj >= gye*refratio) break;
+ gyi = floor(cj / refratio) - gys;
+ for(czi=gzs*refratio;czi<=cze;czi++) {
+ ck=czi;//ck = czi % dw[2];
+ //ck = (ck < 0) ? ck + dw[2] : ck;
+ if ( ck >= gze*refratio) break;
+ gzi = floor(ck / refratio) - gzs;
+ if (refratio == 1 || *(npy_int32*)PyArray_GETPTR3(mask, gxi,gyi,gzi) > 0)
+ {
+ switch (axis) {
+ case 0: x_loc = cyi-cys; y_loc = czi-czs; break;
+ case 1: x_loc = cxi-cxs; y_loc = czi-czs; break;
+ case 2: x_loc = cxi-cys; y_loc = cyi-cys; break;
+ }
+ for(n=0;n<n_fields;n++){
+ *(npy_float64*) PyArray_GETPTR2(c_data[n], x_loc, y_loc)
+ += *(npy_float64*) PyArray_GETPTR3(g_data[n], gxi, gyi, gzi)
+ * dls[n] / refratio;
+ }
+ total += 1;
}
+ }
}
+ }
}
-
Py_DECREF(g_start);
Py_DECREF(c_start);
Py_DECREF(g_dims);
http://bitbucket.org/yt_analysis/yt/changeset/9ddb217547d3/
changeset: 9ddb217547d3
branch: yt
user: MatthewTurk
date: 2011-06-09 00:48:10
summary: These were not meant to be commented out
affected #: 1 file (33 bytes)
--- a/yt/utilities/data_point_utilities.c Wed Jun 08 15:31:54 2011 -0700
+++ b/yt/utilities/data_point_utilities.c Wed Jun 08 15:48:10 2011 -0700
@@ -1251,18 +1251,18 @@
/* It turns out that C89 doesn't define a mechanism for choosing the sign
of the remainder.
*/
- ci=cxi;//ci = (cxi % dw[0]);
- //ci = (ci < 0) ? ci + dw[0] : ci;
+ ci = (cxi % dw[0]);
+ ci = (ci < 0) ? ci + dw[0] : ci;
if ( ci >= gxe*refratio) break;
gxi = floor(ci / refratio) - gxs;
for(cyi=gys*refratio;cyi<=cye;cyi++) {
- cj=cyi;//cj = cyi % dw[1];
- //cj = (cj < 0) ? cj + dw[1] : cj;
+ cj = cyi % dw[1];
+ cj = (cj < 0) ? cj + dw[1] : cj;
if ( cj >= gye*refratio) break;
gyi = floor(cj / refratio) - gys;
for(czi=gzs*refratio;czi<=cze;czi++) {
- ck=czi;//ck = czi % dw[2];
- //ck = (ck < 0) ? ck + dw[2] : ck;
+ ck = czi % dw[2];
+ ck = (ck < 0) ? ck + dw[2] : ck;
if ( ck >= gze*refratio) break;
gzi = floor(ck / refratio) - gzs;
if (refratio == 1 || *(npy_int32*)PyArray_GETPTR3(mask, gxi,gyi,gzi) > 0)
http://bitbucket.org/yt_analysis/yt/changeset/38eee46e2efa/
changeset: 38eee46e2efa
branch: yt
user: MatthewTurk
date: 2011-06-09 01:00:26
summary: Merging in voro++ branch
affected #: 6 files (15.6 KB)
--- a/yt/analysis_modules/level_sets/contour_finder.py Wed Jun 08 15:48:10 2011 -0700
+++ b/yt/analysis_modules/level_sets/contour_finder.py Wed Jun 08 16:00:26 2011 -0700
@@ -30,202 +30,6 @@
import yt.utilities.data_point_utilities as data_point_utilities
import yt.utilities.amr_utils as amr_utils
-class GridConsiderationQueue:
- def __init__(self, white_list, priority_func=None):
- """
- This class exists to serve the contour finder. It ensures that
- we can create a cascading set of queue dependencies, and if
- a grid is touched again ahead of time we can bump it to the top
- of the queue again. It like has few uses.
- """
- self.to_consider = []
- self.considered = set()
- self.n = 0
- self.white_list = set(white_list)
- self.priority_func = priority_func
-
- def add(self, grids, force=False):
- if not hasattr(grids,'size'):
- grids = ensure_list(grids)
- i = self.n
- to_check = self.white_list.intersection(grids)
- if not force: to_check.difference_update(self.considered)
- for g in sorted(to_check, key=self.priority_func):
- try:
- # We only delete from subsequent checks
- del self.to_consider[self.to_consider.index(g, i)]
- self.to_consider.insert(i,g)
- i += 1
- except ValueError:
- self.to_consider.append(g)
-
- def __iter__(self):
- return self
-
- def next(self):
- if self.n >= len(self.to_consider):
- raise StopIteration
- tr = self.to_consider[self.n]
- self.considered.add(tr)
- self.n += 1
- return tr
-
- def progress(self):
- return self.n, len(self.to_consider)
-
-# We want an algorithm that deals with growing a given contour to *all* the
-# cells in a grid.
-
-def old_identify_contours(data_source, field, min_val, max_val, cached_fields=None):
- """
- Given a *data_source*, we will search for topologically connected sets
- in *field* between *min_val* and *max_val*.
- """
- if cached_fields is None: cached_fields = defaultdict(lambda: dict())
- maxn_cells = na.sum([g.ActiveDimensions.prod() for g in data_source._grids])
- contour_ind = na.where( (data_source[field] > min_val)
- & (data_source[field] < max_val))[0]
- np = contour_ind.size
- if np == 0:
- return {}
- cur_max_id = maxn_cells - np
- contour_ids = na.arange(maxn_cells, cur_max_id, -1) + 1 # Minimum of 1
- data_source["tempContours"] = na.ones(data_source[field].shape, dtype='int32')*-1
- mylog.info("Contouring over %s cells with %s candidates", contour_ids[0],np)
- data_source["tempContours"][contour_ind] = contour_ids[:]
- data_source._flush_data_to_grids("tempContours", -1, dtype='int32')
- my_queue = GridConsiderationQueue(white_list = data_source._grids,
- priority_func = lambda g: -1*g["tempContours"].max())
- my_queue.add(data_source._grids)
- for i,grid in enumerate(my_queue):
- mylog.info("Examining %s of %s", *my_queue.progress())
- max_before = grid["tempContours"].max()
- to_get = ["tempContours"]
- if field in cached_fields[grid.id] and \
- not na.any( (cached_fields[grid.id][field] > min_val)
- & (cached_fields[grid.id][field] < max_val)):
- continue
- for f in [field, "GridIndices"]:
- if f not in cached_fields[grid.id]: to_get.append(f)
- cg = grid.retrieve_ghost_zones(1,to_get)
- for f in [field, "GridIndices"]:
- if f in cached_fields[grid.id]:
- cg.data[f] = cached_fields[grid.id][f]
- else:
- cached_fields[grid.id][f] = cg[f]
- local_ind = na.where( (cg[field] > min_val)
- & (cg[field] < max_val)
- & (cg["tempContours"] == -1) )
- if local_ind[0].size > 0:
- kk = na.arange(cur_max_id, cur_max_id-local_ind[0].size, -1)
- cg["tempContours"][local_ind] = kk[:]
- cur_max_id -= local_ind[0].size
- fd = cg["tempContours"].astype('int64')
- fd_original = fd.copy()
- if na.all(fd > -1):
- fd[:] = fd.max()
- else:
- xi_u,yi_u,zi_u = na.where(fd > -1)
- cor_order = na.argsort(-1*fd[(xi_u,yi_u,zi_u)])
- xi = xi_u[cor_order]
- yi = yi_u[cor_order]
- zi = zi_u[cor_order]
- while data_point_utilities.FindContours(fd, xi, yi, zi) < 0: pass
- cg["tempContours"] = fd.copy().astype('float64')
- cg.flush_data("tempContours")
- my_queue.add(cg._grids)
- force_ind = na.unique(cg["GridIndices"][na.where(
- (cg["tempContours"] > fd_original)
- & (cg["GridIndices"] != grid.id-grid._id_offset) )])
- if len(force_ind) > 0:
- my_queue.add(data_source.hierarchy.grids[force_ind.astype('int32')], force=True)
- for ax in 'xyz':
- if not iterable(grid['d%s'%ax]):
- grid['d%s'%ax] = grid['d%s'%ax]*na.ones(grid.ActiveDimensions)
- del data_source.data["tempContours"] # Force a reload from the grids
- data_source.get_data("tempContours", in_grids=True)
- i = 0
- contour_ind = {}
- for contour_id in na.unique(data_source["tempContours"]):
- if contour_id == -1: continue
- contour_ind[i] = na.where(data_source["tempContours"] == contour_id)
- mylog.debug("Contour id %s has %s cells", i, contour_ind[i][0].size)
- i += 1
- mylog.info("Identified %s contours between %0.5e and %0.5e",
- len(contour_ind.keys()),min_val,max_val)
- for grid in chain(data_source._grids, cg._grids):
- grid.data.pop("tempContours", None)
- del data_source.data["tempContours"]
- return contour_ind
-
-def check_neighbors(data_object, field="Contours"):
- """
- This method is a means of error checking in the contour finder.
- """
- n_bad = na.zeros(1, dtype='int32')
- for cid in na.unique(data_object[field]):
- if cid == -1: continue
- ids = na.where(data_object[field] == cid)
- mx = data_object['x'][ids].copy()
- my = data_object['y'][ids].copy()
- mz = data_object['z'][ids].copy()
- mdx = data_object['dx'][ids].copy()
- mdy = data_object['dy'][ids].copy()
- mdz = data_object['dz'][ids].copy()
- grid_ids_m = data_object['GridIndices'][ids].copy()
- grid_levels_m = data_object['GridLevel'][ids].copy()
- mp = mx.size
- ids = na.where( (data_object[field] != cid)
- & (data_object[field] >= 0 ))
- nx = data_object['x'][ids].copy()
- ny = data_object['y'][ids].copy()
- nz = data_object['z'][ids].copy()
- ndx = data_object['dx'][ids].copy()
- ndy = data_object['dy'][ids].copy()
- ndz = data_object['dz'][ids].copy()
- grid_ids_n = data_object['GridIndices'][ids].copy()
- grid_levels_n = data_object['GridLevel'][ids].copy()
- np = nx.size
- weave.inline(check_cell_distance,
- ['mx','my','mz','mdx','mdy','mdz','mp',
- 'nx','ny','nz','ndx','ndy','ndz','np','n_bad',
- 'grid_ids_m', 'grid_levels_m', 'grid_ids_n', 'grid_levels_n'],
- compiler='gcc', type_converters=converters.blitz,
- auto_downcast=0, verbose=2)
- return n_bad[0]
-
-check_cell_distance = \
-r"""
-using namespace std;
-int i, j, k;
-long double cell_dist, rad_m, rad_n;
-k=0;
-for(i=0;i<mp;i++){
- for(j=0;j<np;j++){
- /*
- cell_dist = sqrtl(pow(mx(i)-nx(j),2) +
- pow(my(i)-ny(j),2) +
- pow(mz(i)-nz(j),2));
- rad_m = sqrtl(pow(mdx(i),2) +
- pow(mdy(i),2) +
- pow(mdz(i),2));
- rad_n = sqrtl(pow(ndx(j),2) +
- pow(ndy(j),2) +
- pow(ndz(j),2));
- */
- //if(cell_dist > 1.01 * (rad_n/2.0+rad_m/2.0)) continue;
- if(fabsl(mx(i)-nx(j))>(mdx(i)+ndx(j))/2.0) continue;
- if(fabsl(my(i)-ny(j))>(mdy(i)+ndy(j))/2.0) continue;
- if(fabsl(mz(i)-nz(j))>(mdz(i)+ndz(j))/2.0) continue;
- k++;
- break;
- cout << cell_dist << "\t" << 1.01*(rad_n/2.0+rad_m/2.0) << "\t";
- cout << grid_ids_m(i) << "\t" << grid_ids_n(j) << endl;
- }
-}
-n_bad(0) += k;
-"""
-
def coalesce_join_tree(jtree1):
joins = defaultdict(set)
nj = jtree1.shape[0]
--- a/yt/utilities/_amr_utils/ContourFinding.pyx Wed Jun 08 15:48:10 2011 -0700
+++ b/yt/utilities/_amr_utils/ContourFinding.pyx Wed Jun 08 16:00:26 2011 -0700
@@ -26,6 +26,7 @@
import numpy as np
cimport numpy as np
cimport cython
+from stdlib cimport malloc, free
cdef extern from "math.h":
double fabs(double x)
@@ -38,6 +39,50 @@
if i0 < i1: return i0
return i1
+cdef extern from "union_find.h":
+ ctypedef struct forest_node:
+ void *value
+ forest_node *parent
+ int rank
+
+ forest_node* MakeSet(void* value)
+ void Union(forest_node* node1, forest_node* node2)
+ forest_node* Find(forest_node* node)
+
+ctypedef struct CellIdentifier:
+ np.int64_t hindex
+ int level
+
+cdef class GridContourContainer:
+ cdef np.int64_t dims[3]
+ cdef np.int64_t start_indices[3]
+ cdef forest_node **join_tree
+ cdef np.int64_t ncells
+
+ def __init__(self, dimensions, indices):
+ cdef int i
+ self.ncells = 1
+ for i in range(3):
+ self.ncells *= dimensions[i]
+ self.dims[i] = dimensions[i]
+ self.start_indices[i] = indices[i]
+ self.join_tree = <forest_node **> malloc(sizeof(forest_node)
+ * self.ncells)
+ for i in range(self.ncells): self.join_tree[i] = NULL
+
+ def __dealloc__(self):
+ cdef int i
+ for i in range(self.ncells):
+ if self.join_tree[i] != NULL: free(self.join_tree[i])
+ free(self.join_tree)
+
+ #def construct_join_tree(self,
+ # np.ndarray[np.float64_t, ndim=3] field,
+ # np.ndarray[np.bool_t, ndim=3] mask):
+ # # This only looks at the components of the grid that are actually
+ # # inside this grid -- boundary conditions are handled later.
+ # pass
+
#@cython.boundscheck(False)
#@cython.wraparound(False)
def construct_boundary_relationships(
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/union_find.c Wed Jun 08 16:00:26 2011 -0700
@@ -0,0 +1,66 @@
+/* Copyright (c) 2011 the authors listed at the following URL, and/or
+the authors of referenced articles or incorporated external code:
+http://en.literateprograms.org/Disjoint_set_data_structure_(C)?action=history&offset=20080516180553
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+Retrieved from: http://en.literateprograms.org/Disjoint_set_data_structure_(C)?oldid=13366
+*/
+
+#include <stdlib.h>
+
+#include "union_find.h"
+
+forest_node* MakeSet(void* value) {
+ forest_node* node = malloc(sizeof(forest_node));
+ node->value = value;
+ node->parent = NULL;
+ node->rank = 0;
+ return node;
+}
+
+void Union(forest_node* node1, forest_node* node2) {
+ if (node1->rank > node2->rank) {
+ node2->parent = node1;
+ } else if (node2->rank > node1->rank) {
+ node1->parent = node2;
+ } else { /* they are equal */
+ node2->parent = node1;
+ node1->rank++;
+ }
+}
+
+forest_node* Find(forest_node* node) {
+ forest_node* temp;
+ /* Find the root */
+ forest_node* root = node;
+ while (root->parent != NULL) {
+ root = root->parent;
+ }
+ /* Update the parent pointers */
+ while (node->parent != NULL) {
+ temp = node->parent;
+ node->parent = root;
+ node = temp;
+ }
+ return root;
+}
+
+
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/union_find.h Wed Jun 08 16:00:26 2011 -0700
@@ -0,0 +1,41 @@
+/* Copyright (c) 2011 the authors listed at the following URL, and/or
+the authors of referenced articles or incorporated external code:
+http://en.literateprograms.org/Disjoint_set_data_structure_(C)?action=history&offset=20080516180553
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+Retrieved from: http://en.literateprograms.org/Disjoint_set_data_structure_(C)?oldid=13366
+*/
+
+#ifndef _UNION_FIND_H_
+#define _UNION_FIND_H_
+
+typedef struct forest_node_t {
+ void* value;
+ struct forest_node_t* parent;
+ int rank;
+} forest_node;
+
+forest_node* MakeSet(void* value);
+void Union(forest_node* node1, forest_node* node2);
+forest_node* Find(forest_node* node);
+
+#endif /* #ifndef _UNION_FIND_H_ */
+
--- a/yt/utilities/setup.py Wed Jun 08 15:48:10 2011 -0700
+++ b/yt/utilities/setup.py Wed Jun 08 16:00:26 2011 -0700
@@ -161,7 +161,8 @@
config.add_extension("amr_utils",
["yt/utilities/amr_utils.pyx",
"yt/utilities/_amr_utils/FixedInterpolator.c",
- "yt/utilities/_amr_utils/kdtree.c"] +
+ "yt/utilities/_amr_utils/kdtree.c",
+ "yt/utilities/_amr_utils/union_find.c"] +
glob.glob("yt/utilities/_amr_utils/healpix_*.c"),
define_macros=[("PNG_SETJMP_NOT_SUPPORTED", True)],
include_dirs=["yt/utilities/_amr_utils/", png_inc,
@@ -172,6 +173,10 @@
glob.glob("yt/utilities/_amr_utils/*.h") +
glob.glob("yt/utilities/_amr_utils/*.c"),
)
+ #config.add_extension("voropp",
+ # ["yt/utilities/voropp.pyx"],
+ # language="c++",
+ # include_dirs=["yt/utilities/voro++"])
config.add_extension("libconfig_wrapper",
["yt/utilities/libconfig_wrapper.pyx"] +
glob.glob("yt/utilities/_libconfig/*.c"),
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/voropp.pyx Wed Jun 08 16:00:26 2011 -0700
@@ -0,0 +1,69 @@
+"""
+Wrapping code for voro++
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2011 Matthew Turk. 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/>.
+"""
+
+from cython.operator cimport dereference as deref, preincrement as inc
+from libc.stdlib cimport malloc, free, abs, calloc, labs
+cimport libcpp
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+cdef extern from "voro++.cc":
+ cdef cppclass container:
+ container(double xmin, double xmax, double ymin, double ymax,
+ double zmin, double zmax, int nx, int ny, int nz,
+ libcpp.bool xper, libcpp.bool yper, libcpp.bool zper, int alloc)
+ void put(int n, double x, double y, double z)
+ void store_cell_volumes(double *vols)
+
+cdef class VoronoiVolume:
+ cdef container *my_con
+ cdef int npart
+ def __init__(self, xi, yi, zi):
+ self.my_con = new container(0.0, 1.0, 0.0, 1.0, 0.0, 1.0,
+ xi, yi, zi, False, False, False, 8)
+ self.npart = 0
+
+ def __dealloc__(self):
+ del self.my_con
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def add_array(self, np.ndarray[np.float64_t, ndim=1] xpos,
+ np.ndarray[np.float64_t, ndim=1] ypos,
+ np.ndarray[np.float64_t, ndim=1] zpos):
+ cdef int i
+ for i in range(xpos.shape[0]):
+ self.my_con.put(self.npart, xpos[i], ypos[i], zpos[i])
+ self.npart += 1
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def get_volumes(self):
+ cdef np.ndarray vol = np.zeros(self.npart, 'double')
+ cdef double *vdouble = <double *> vol.data
+ self.my_con.store_cell_volumes(vdouble)
+ return vol
http://bitbucket.org/yt_analysis/yt/changeset/1aaba47d7105/
changeset: 1aaba47d7105
branch: yt
user: MatthewTurk
date: 2011-06-09 01:01:12
summary: Removing new method for the moment; needs more thought
affected #: 1 file (153 bytes)
--- a/yt/data_objects/data_containers.py Wed Jun 08 16:00:26 2011 -0700
+++ b/yt/data_objects/data_containers.py Wed Jun 08 16:01:12 2011 -0700
@@ -2139,9 +2139,6 @@
self._dls = {}
self.domain_width = na.rint((self.pf.domain_right_edge -
self.pf.domain_left_edge)/self.dds).astype('int64')
- if not na.any((self.global_startindex + self.ActiveDimensions) >
- self.domain_width):
- self.domain_width[:] = 0
self._refresh_data()
def _get_list_of_grids(self):
http://bitbucket.org/yt_analysis/yt/changeset/38c841d2054c/
changeset: 38c841d2054c
branch: yt
user: MatthewTurk
date: 2011-06-09 19:55:14
summary: Making the width in the mapserver set to the entire visible window.
affected #: 3 files (140 bytes)
--- a/yt/gui/reason/html/map_index.html Wed Jun 08 16:01:12 2011 -0700
+++ b/yt/gui/reason/html/map_index.html Thu Jun 09 10:55:14 2011 -0700
@@ -8,6 +8,8 @@
<script type="text/javascript">
$(document).ready(function() {
// initialize the map on the "map" div with a given center and zoom
+ $("#map").width($(window).width());
+ $("#map").height($(window).height());
var map = new L.Map('map', {
center: new L.LatLng(0.0, 0.0),
zoom: 0,
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Wed Jun 08 16:01:12 2011 -0700
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Jun 09 10:55:14 2011 -0700
@@ -58,6 +58,7 @@
cdef extern from "math.h":
double exp(double x)
float expf(float x)
+ long double expl(long double x)
double floor(double x)
double ceil(double x)
double fmod(double x, double y)
@@ -742,7 +743,7 @@
cdef void add_stars(self, kdtree_utils.kdres *ballq,
np.float64_t dt, np.float64_t pos[3], np.float64_t *rgba):
cdef int i, n, ns
- cdef double px, py, pz
+ cdef np.float64_t px, py, pz
cdef np.float64_t gexp, gaussian
cdef np.float64_t* colors = NULL
ns = kdtree_utils.kd_res_size(ballq)
@@ -754,7 +755,7 @@
gexp = (px - pos[0])*(px - pos[0]) \
+ (py - pos[1])*(py - pos[1]) \
+ (pz - pos[2])*(pz - pos[2])
- gaussian = self.star_coeff * exp(-gexp/self.star_sigma_num)
+ gaussian = self.star_coeff * expl(-gexp/self.star_sigma_num)
for i in range(3): rgba[i] += gaussian*dt*colors[i]
kdtree_utils.kd_res_rewind(ballq)
--- a/yt/utilities/data_point_utilities.c Wed Jun 08 16:01:12 2011 -0700
+++ b/yt/utilities/data_point_utilities.c Thu Jun 09 10:55:14 2011 -0700
@@ -1236,7 +1236,7 @@
for(n=0;n<n_fields;n++){
*(npy_float64*) PyArray_GETPTR2(c_data[n], x_loc+ri, y_loc+rj)
+= *(npy_float64*) PyArray_GETPTR3(g_data[n],
- gxi-gxs, gyi-gys, gzi-gzs) * dls[n] / refratio;
+ gxi-gxs, gyi-gys, gzi-gzs) * dls[n];
}
}
}
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