[Yt-svn] yt: 3 new changesets
hg at spacepope.org
hg at spacepope.org
Mon Aug 9 12:22:52 PDT 2010
hg Repository: yt
details: yt/rev/7d4b1935a88f
changeset: 1914:7d4b1935a88f
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Aug 03 23:51:38 2010 -0700
description:
Adding some initial passes at the signature calculation in Cython
hg Repository: yt
details: yt/rev/3127297bff7d
changeset: 1915:3127297bff7d
user: Matthew Turk <matthewturk at gmail.com>
date:
Mon Aug 09 07:54:01 2010 -0700
description:
A clearer signature and second derivative calculation, pulled from Enzo
directly
hg Repository: yt
details: yt/rev/cc5baf021d41
changeset: 1916:cc5baf021d41
user: Matthew Turk <matthewturk at gmail.com>
date:
Mon Aug 09 12:21:54 2010 -0700
description:
Another set of iterations on the ramses reader protosubgrid calculation. I can
now split and coalesce, but it requires some code that currently only exists
external to this wrapper.
diffstat:
yt/ramses_reader.pyx | 147 ++++++++++++++++++++++++++++++++++++++++++++++++-
1 files changed, 145 insertions(+), 2 deletions(-)
diffs (173 lines):
diff -r 007544606172 -r cc5baf021d41 yt/ramses_reader.pyx
--- a/yt/ramses_reader.pyx Sun Aug 01 22:53:37 2010 -0700
+++ b/yt/ramses_reader.pyx Mon Aug 09 12:21:54 2010 -0700
@@ -27,12 +27,20 @@
# Cython wrapping code for Oliver Hahn's RAMSES reader
from cython.operator cimport dereference as deref, preincrement as inc
-from libc.stdlib cimport malloc, free, abs
+from libc.stdlib cimport malloc, free, abs, calloc, labs
import numpy as np
cimport numpy as np
cimport cython
+cdef inline np.int64_t i64max(np.int64_t i0, np.int64_t i1):
+ if i0 > i1: return i0
+ return i1
+
+cdef inline np.int64_t i64min(np.int64_t i0, np.int64_t i1):
+ if i0 < i1: return i0
+ return i1
+
cdef extern from "<vector>" namespace "std":
cdef cppclass vector[T]:
cppclass iterator:
@@ -569,7 +577,7 @@
cdef int i
cdef np.ndarray[np.float64_t, ndim=3] tr = np.empty((2,2,2), dtype='float64',
- order='F")
+ order='F')
cdef tree_iterator grid_it, grid_end
cdef double* data = <double*> tr.data
@@ -587,3 +595,138 @@
for i in range(8):
data[i] = local_hydro_data.m_var_array[level][8*grid_id+i]
return tr
+
+def identify_new_subgrids_by_signature(
+ np.ndarray[np.int64_t, ndim=2] left_edges, # In integer indices
+ np.ndarray[np.int64_t, ndim=2] right_edges, # In integer indices
+ np.ndarray[np.int64_t, ndim=2] grid_dimensions):
+ # We operate slightly differently than Enzo does. We can't afford to store
+ # all the flagging fields in memory; so, we operate only on 1D signatures.
+ # So, we have a list of grids that we want to cluster, and then we pass
+ # that around and determine which grids are appropriate.
+ # We start with a proto subgrid that contains the entire domain.
+ #cdef ProtoSubgrid *psg = ProtoSubgrid()
+ pass
+
+cdef class ProtoSubgrid:
+ cdef np.int64_t *signature[3]
+ cdef np.int64_t left_edge[3]
+ cdef np.int64_t right_edge[3]
+ cdef np.int64_t dimensions[3]
+ cdef np.float64_t efficiency
+ cdef public object sigs
+
+ #@cython.boundscheck(False)
+ #@cython.wraparound(False)
+ def __cinit__(self,
+ np.ndarray[np.int64_t, ndim=1] left_index,
+ np.ndarray[np.int64_t, ndim=1] dimensions,
+ np.ndarray[np.int64_t, ndim=2] left_edges,
+ np.ndarray[np.int64_t, ndim=2] right_edges,
+ np.ndarray[np.int64_t, ndim=2] grid_dimensions):
+ # This also includes the shrinking step.
+ cdef int i, ci, ng = left_edges.shape[0]
+ cdef np.ndarray temp_arr
+ cdef int l0, r0, l1, r1, l2, r2, i0, i1, i2
+ cdef np.int64_t temp_l[3], temp_r[3], ncells
+ cdef np.float64_t efficiency
+ self.sigs = []
+ for i in range(3):
+ temp_l[i] = left_index[i] + dimensions[i]
+ temp_r[i] = left_index[i]
+ self.signature[i] = NULL
+ for gi in range(ng):
+ if left_edges[gi,0] > left_index[0]+dimensions[0] or \
+ right_edges[gi,0] < left_index[0] or \
+ left_edges[gi,1] > left_index[1]+dimensions[1] or \
+ right_edges[gi,1] < left_index[1] or \
+ left_edges[gi,2] > left_index[2]+dimensions[2] or \
+ right_edges[gi,2] < left_index[2]:
+ #print "Skipping grid", gi, "which lies outside out box"
+ continue
+ for i in range(3):
+ temp_l[i] = i64min(left_edges[gi,i], temp_l[i])
+ temp_r[i] = i64max(right_edges[gi,i], temp_r[i])
+ for i in range(3):
+ self.left_edge[i] = i64max(temp_l[i], left_index[i])
+ self.right_edge[i] = i64min(temp_r[i], left_index[i] + dimensions[i])
+ self.dimensions[i] = self.right_edge[i] - self.left_edge[i]
+ if self.dimensions[i] <= 0:
+ self.efficiency = -1.0
+ return
+ self.sigs.append(np.zeros(self.dimensions[i], 'int64'))
+ #print self.sigs[0].size, self.sigs[1].size, self.sigs[2].size
+
+ # My guess is that this whole loop could be done more efficiently.
+ # However, this is clear and straightforward, so it is a good first
+ # pass.
+ cdef np.ndarray[np.int64_t, ndim=1] sig0, sig1, sig2
+ sig0 = self.sigs[0]
+ sig1 = self.sigs[1]
+ sig2 = self.sigs[2]
+ efficiency = 0.0
+ for gi in range(ng):
+ nnn = 0
+ for l0 in range(grid_dimensions[gi, 0]):
+ i0 = left_edges[gi, 0] + l0
+ if i0 < self.left_edge[0]: continue
+ if i0 >= self.right_edge[0]: break
+ for l1 in range(grid_dimensions[gi, 1]):
+ i1 = left_edges[gi, 1] + l1
+ if i1 < self.left_edge[1]: continue
+ if i1 >= self.right_edge[1]: break
+ for l2 in range(grid_dimensions[gi, 2]):
+ i2 = left_edges[gi, 2] + l2
+ if i2 < self.left_edge[2]: continue
+ if i2 >= self.right_edge[2]: break
+ i = i0 - self.left_edge[0]
+ sig0[i] += 1
+ i = i1 - self.left_edge[1]
+ sig1[i] += 1
+ i = i2 - self.left_edge[2]
+ sig2[i] += 1
+ efficiency += 1
+
+ for i in range(3): efficiency /= self.dimensions[i]
+ #print "Efficiency is %0.3e" % (efficiency)
+ self.efficiency = efficiency
+
+ def get_efficiency(self):
+ return self.efficiency
+
+ def find_split(self, int ax):
+ # First look for zeros
+ cdef int i, center
+ center = self.dimensions[i] / 2
+ cdef np.int64_t strength, zcstrength, zcp
+ cdef np.ndarray[np.int64_t] sig = self.sigs[ax]
+ for i in range(self.dimensions[ax]):
+ if sig[i] == 0:
+ #print "zero: %s (%s)" % (i, self.dimensions[ax])
+ return i
+ cdef np.int64_t *sig2d = <np.int64_t *> calloc(
+ sizeof(np.int64_t), self.dimensions[ax])
+ sig2d[0] = sig2d[self.dimensions[ax]-1] = 0
+ for i in range(1, self.dimensions[ax] - 1):
+ sig2d[i] = (sig[i-1] - 2*sig[i] + sig[i+1])
+ zcstrength = 0
+ zcp = 0
+ for i in range(1, self.dimensions[ax] - 1):
+ if sig2d[i] * sig2d[i+1] <= 0:
+ strength = labs(sig2d[i] - sig2d[i+1])
+ if (strength > zcstrength) or \
+ (strength == zcstrength and (abs(center - i) <
+ abs(center - zcp))):
+ zcstrength = strength
+ zcp = i
+ #print "zcp: %s (%s)" % (zcp, self.dimensions[ax])
+ return zcp
+
+ def get_properties(self):
+ cdef np.ndarray[np.int64_t, ndim=2] tr = np.empty((3,3), dtype='int64')
+ cdef int i
+ for i in range(3):
+ tr[0,i] = self.left_edge[i]
+ tr[1,i] = self.right_edge[i]
+ tr[2,i] = self.dimensions[i]
+ return tr
More information about the yt-svn
mailing list