[h5py] 102/455: Rewrite indexing code; boolean-array indexing allowed. Change h5s list functions to use ndarrays.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Thu Jul 2 18:19:22 UTC 2015
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to annotated tag 1.3.0
in repository h5py.
commit 83ec9d5bfd4e137498d0785802bf9577327dffd3
Author: andrewcollette <andrew.collette at gmail.com>
Date: Tue Aug 19 01:54:56 2008 +0000
Rewrite indexing code; boolean-array indexing allowed. Change h5s list functions to use ndarrays.
---
h5py/h5s.pxd | 1 +
h5py/h5s.pyx | 146 ++++++++++++++++---------------------------
h5py/highlevel.py | 85 ++++++++++++-------------
h5py/numpy.pxd | 6 +-
h5py/tests/test_h5s.py | 7 ++-
h5py/tests/test_highlevel.py | 7 ++-
h5py/utils.pxd | 2 +
h5py/utils.pyx | 51 +++++++++++++++
h5py/utils_hl.py | 71 ++++++++++-----------
9 files changed, 195 insertions(+), 181 deletions(-)
diff --git a/h5py/h5s.pxd b/h5py/h5s.pxd
index 6476555..d4f4bab 100644
--- a/h5py/h5s.pxd
+++ b/h5py/h5s.pxd
@@ -20,6 +20,7 @@
include "std_defs.pxi"
from h5 cimport class ObjectID
+from numpy cimport ndarray
cdef class SpaceID(ObjectID):
pass
diff --git a/h5py/h5s.pyx b/h5py/h5s.pyx
index d9bd34c..f440f9e 100644
--- a/h5py/h5s.pyx
+++ b/h5py/h5s.pyx
@@ -16,7 +16,7 @@
# Pyrex compile-time imports
from utils cimport require_tuple, require_list, convert_dims, convert_tuple, \
- emalloc, efree, pybool
+ emalloc, efree, pybool, create_numpy_hsize, create_hsize_array
# Runtime imports
import h5
@@ -355,80 +355,56 @@ cdef class SpaceID(ObjectID):
return H5Sget_select_elem_npoints(self.id)
def get_select_elem_pointlist(self):
- """ () => LIST elements_list
+ """ () => NDARRAY elements_list
- Get a list of all selected elements, in point-selection mode.
- List entries are <rank>-length tuples containing point coordinates.
+ Get a list of all selected elements. Return is a Numpy array of
+ unsigned ints, with shape (<npoints>, <space rank).
"""
- cdef int rank
- cdef hssize_t npoints
- cdef hsize_t *buf
- cdef int i_point
- cdef int i_entry
+ cdef hsize_t dims[2]
+ cdef ndarray buf
- npoints = H5Sget_select_elem_npoints(self.id)
- if npoints == 0:
- return []
+ dims[0] = H5Sget_select_elem_npoints(self.id)
+ dims[1] = H5Sget_simple_extent_ndims(self.id)
- rank = H5Sget_simple_extent_ndims(self.id)
-
- buf = <hsize_t*>emalloc(sizeof(hsize_t)*rank*npoints)
+ buf = create_numpy_hsize(2, dims)
- try:
- H5Sget_select_elem_pointlist(self.id, 0, <hsize_t>npoints, buf)
+ H5Sget_select_elem_pointlist(self.id, 0, dims[0], <hsize_t*>buf.data)
- retlist = []
- for i_point from 0<=i_point<npoints:
- tmp_tpl = []
- for i_entry from 0<=i_entry<rank:
- tmp_tpl.append( long( buf[i_point*rank + i_entry] ) )
- retlist.append(tuple(tmp_tpl))
+ return buf
- finally:
- efree(buf)
+ def select_elements(self, object coords, int op=H5S_SELECT_SET):
+ """ (SEQUENCE coords, INT op=SELECT_SET)
- return retlist
+ Select elements by specifying coordinates points. The argument
+ "coords" may be an ndarray or any nested sequence which can be
+ converted to an array of uints with the shape:
- def select_elements(self, object coord_list, int op=H5S_SELECT_SET):
- """ (LIST coord_list, INT op=SELECT_SET)
+ (<npoints>, <space rank>)
+
+ Examples:
+ >>> obj.shape
+ (10, 10)
+ >>> obj.select_elements([(1,2), (3,4), (5,9)])
- Select elements using a list of points. List entries should be
- tuples containing point coordinates. A zero-length list is
- apparently not allowed by the HDF5 library.
+ A zero-length selection (i.e. shape (0, <rank>)) is not allowed by
+ the HDF5 library.
"""
- cdef size_t nelements # Number of point coordinates
- cdef hsize_t *coords # Contiguous 2D array nelements x rank x sizeof(hsize_t)
-
- cdef int rank
- cdef int i_point
- cdef int i_entry
- coords = NULL
-
- require_list(coord_list, 0, -1, "coord_list")
- nelements = len(coord_list)
-
- rank = H5Sget_simple_extent_ndims(self.id)
-
- # The docs say this should be an hsize_t**, but it seems that
- # HDF5 expects the coordinates to be a static, contiguous
- # array. We'll simulate that by malloc'ing a contiguous chunk
- # and using pointer arithmetic to initialize it.
- coords = <hsize_t*>emalloc(sizeof(hsize_t)*rank*nelements)
+ cdef ndarray hcoords
+ cdef size_t nelements
- try:
- for i_point from 0<=i_point<nelements:
+ # The docs say the selection list should be an hsize_t**, but it seems
+ # that HDF5 expects the coordinates to be a static, contiguous
+ # array. We simulate that by creating a contiguous NumPy array of
+ # a compatible type and initializing it to the input.
- tpl = coord_list[i_point]
- lmsg = "List element %d" % i_point
- require_tuple(tpl, 0, rank, lmsg)
+ hcoords = create_hsize_array(coords)
- for i_entry from 0<=i_entry<rank:
- coords[(i_point*rank) + i_entry] = tpl[i_entry]
+ if hcoords.nd != 2 or hcoords.dimensions[1] != H5Sget_simple_extent_ndims(self.id):
+ raise ValueError("Coordinate array must have shape (<npoints>, %d)" % self.get_simple_extent_ndims())
- H5Sselect_elements(self.id, <H5S_seloper_t>op, nelements, <hsize_t**>coords)
+ nelements = hcoords.dimensions[0]
- finally:
- efree(coords)
+ H5Sselect_elements(self.id, <H5S_seloper_t>op, nelements, <hsize_t**>hcoords.data)
# === Hyperslab selection functions =======================================
@@ -440,46 +416,30 @@ cdef class SpaceID(ObjectID):
return H5Sget_select_hyper_nblocks(self.id)
def get_select_hyper_blocklist(self):
- """ () => LIST hyperslab_blocks
+ """ () => NDARRAY hyperslab_blocks
+
+ Get the current hyperslab selection. The returned array has shape
+
+ (<npoints>, 2, <rank>)
+
+ and can be interpreted as a nested sequence:
+
+ [ (corner_coordinate_1, opposite_coordinate_1), ... ]
- Get a Python list containing selected hyperslab blocks.
- List entries are 2-tuples in the form:
- ( corner_coordinate, opposite_coordinate )
- where corner_coordinate and opposite_coordinate are <rank>-length
- tuples.
+ with length equal to the total number of blocks.
"""
- cdef hssize_t nblocks
- cdef herr_t retval
- cdef hsize_t *buf
+ cdef hsize_t dims[3] # 0=nblocks 1=(#2), 2=rank
+ cdef ndarray buf
- cdef int rank
- cdef int i_block
- cdef int i_entry
+ dims[0] = H5Sget_select_hyper_nblocks(self.id)
+ dims[1] = 2
+ dims[2] = H5Sget_simple_extent_ndims(self.id)
- rank = H5Sget_simple_extent_ndims(self.id)
- nblocks = H5Sget_select_hyper_nblocks(self.id)
+ buf = create_numpy_hsize(3, dims)
- buf = <hsize_t*>emalloc(sizeof(hsize_t)*2*rank*nblocks)
-
- try:
- H5Sget_select_hyper_blocklist(self.id, 0, nblocks, buf)
-
- outlist = []
- for i_block from 0<=i_block<nblocks:
- corner_list = []
- opposite_list = []
- for i_entry from 0<=i_entry<(2*rank):
- entry = long(buf[ i_block*(2*rank) + i_entry])
- if i_entry < rank:
- corner_list.append(entry)
- else:
- opposite_list.append(entry)
- outlist.append( (tuple(corner_list), tuple(opposite_list)) )
- finally:
- efree(buf)
+ H5Sget_select_hyper_blocklist(self.id, 0, dims[0], <hsize_t*>buf.data)
- return outlist
-
+ return buf
def select_hyperslab(self, object start, object count, object stride=None,
object block=None, int op=H5S_SELECT_SET):
diff --git a/h5py/highlevel.py b/h5py/highlevel.py
index 506c94f..c8b773d 100644
--- a/h5py/highlevel.py
+++ b/h5py/highlevel.py
@@ -50,11 +50,10 @@ import os
import numpy
import inspect
import threading
-from weakref import WeakValueDictionary
from h5py import h5, h5f, h5g, h5s, h5t, h5d, h5a, h5p, h5z, h5i, config
from h5py.h5 import H5Error
-from utils_hl import slicer, hbasename, strhdr, strlist
+from utils_hl import slice_select, hbasename, strhdr, strlist
from browse import _H5Browser
@@ -528,46 +527,35 @@ class Dataset(HLObject):
ds[0:5:2, ..., 0:2, "a", "b"]
"""
with self._lock:
- start, count, stride, names = slicer(self.shape, args)
- if not (len(start) == len(count) == len(stride) == self.id.rank):
- raise ValueError("Indices do not match dataset rank (%d)" % self.id.rank)
+ args = args if isinstance(args, tuple) else (args,)
- htype = self.id.get_type()
- if len(names) > 0:
- if htype.get_class() == h5t.COMPOUND:
+ # Sort field indices from the slicing
+ names = tuple(x for x in args if isinstance(x, str))
+ slices = tuple(x for x in args if not isinstance(x, str))
- subtypes = {}
- for idx in range(htype.get_nmembers()):
- subtypes[htype.get_member_name(idx)] = htype.get_member_type(idx)
-
- for name in names:
- if name not in subtypes:
- raise ValueError("Field %s does not appear in this type." % name)
-
- insertlist = [(name, subtypes[name].get_size()) for name in names]
- totalsize = sum([x[1] for x in insertlist])
+ fspace = self.id.get_space()
- mtype = h5t.create(h5t.COMPOUND, totalsize)
+ # Perform selection on the dataset and retrieve the
+ # dataspace for NumPy to use
+ mspace = slice_select(fspace, slices)
- offset = 0
- for name, size in insertlist:
- mtype.insert(name, offset, subtypes[name])
- offset += size
- else:
- raise ValueError("This dataset has no named fields.")
+ # Create NumPy datatype for read, using the named type restrictions
+ basetype = self.id.dtype
+
+ if len(names) == 0:
+ new_dtype = basetype
else:
- mtype = htype
+ for name in names:
+ if not name in basetype.names:
+ raise ValueError("Field %s does not appear in this type." % name)
- fspace = self.id.get_space()
- if fspace.get_simple_extent_type() == h5s.SCALAR:
- fspace.select_all()
- else:
- fspace.select_hyperslab(start, count, stride)
- mspace = h5s.create_simple(count)
+ new_dtype = numpy.dtype([(name, basetype.fields[name]) for name in names])
- arr = numpy.ndarray(count, mtype.dtype)
+ # Create the holder array
+ arr = numpy.ndarray(mspace.shape, new_dtype)
+ # Perform the actual read
self.id.read(mspace, fspace, arr)
if len(names) == 1:
@@ -575,30 +563,35 @@ class Dataset(HLObject):
return arr[names[0]].squeeze()
return arr.squeeze()
+
def __setitem__(self, args, val):
""" Write to the HDF5 dataset from an Numpy array. The shape of the
Numpy array must match the shape of the selection, and the Numpy
array's datatype must be convertible to the HDF5 datatype.
"""
with self._lock:
- start, count, stride, names = slicer(self.shape, args)
- if len(names) != 0:
- raise ValueError("Field name selections are not allowed for write.")
- val = numpy.array(val, dtype=self.dtype)
+ args = args if isinstance(args, tuple) else (args,)
- if count != val.shape:
- # Allow assignments (1,10) => (10,)
- if numpy.product(count) != numpy.product(val.shape):
- raise ValueError("Selection (%s) must be compatible with target (%s)" % (str(count), str(val.shape)))
- else:
- val = val.reshape(count)
+ # Sort field indices from the slicing
+ names = tuple(x for x in args if isinstance(x, str))
+ slices = tuple(x for x in args if not isinstance(x, str))
+
+ if len(names) != 0:
+ raise NotImplementedError("Field name selections are not yet allowed for write.")
+
+ val = numpy.array(val) # So that you can assign scalars, sequences
fspace = self.id.get_space()
- fspace.select_hyperslab(start, count, stride)
- mspace = h5s.create_simple(val.shape)
- self.id.write(mspace, fspace, numpy.array(val))
+ if val.shape == ():
+ mspace = h5s.create(h5s.SCALAR)
+ else:
+ mspace = h5s.create_simple(val.shape)
+
+ slice_select(fspace, args)
+
+ self.id.write(mspace, fspace, val)
def __str__(self):
with self._lock:
diff --git a/h5py/numpy.pxd b/h5py/numpy.pxd
index e8e1784..b8681fa 100644
--- a/h5py/numpy.pxd
+++ b/h5py/numpy.pxd
@@ -62,7 +62,8 @@ cdef extern from "numpy/arrayobject.h":
NPY_FLOAT32, NPY_FLOAT64, NPY_COMPLEX64, NPY_COMPLEX128
cdef enum:
- NPY_WRITEABLE, NPY_ALIGNED, NPY_C_CONTIGUOUS
+ NPY_WRITEABLE, NPY_ALIGNED, NPY_C_CONTIGUOUS, NPY_CONTIGUOUS,
+ NPY_FORCECAST, NPY_NOTSWAPPED
# Classes
ctypedef extern class numpy.dtype [object PyArray_Descr]:
@@ -98,6 +99,9 @@ cdef extern from "numpy/arrayobject.h":
int PyArray_CheckScalar(object sclr)
void PyArray_ScalarAsCtype(object sclr, void* ptr)
+ object PyArray_SimpleNew(int nd, npy_intp* dims, int typenum)
+ object PyArray_ContiguousFromAny(object arr, int typenum, int min_depth, int max_depth)
+ object PyArray_FROM_OTF(object arr, int typenum, int requirements)
# The NumPy initialization function
void import_array()
diff --git a/h5py/tests/test_h5s.py b/h5py/tests/test_h5s.py
index 72c8d78..cebca56 100644
--- a/h5py/tests/test_h5s.py
+++ b/h5py/tests/test_h5s.py
@@ -11,6 +11,7 @@
#-
import unittest
+import numpy
from h5py import *
from h5py.h5 import H5Error
@@ -129,12 +130,12 @@ class TestH5S(unittest.TestCase):
def test_elements(self):
- pointlist= [(0,0), (15,98), (4,17), (67,32)]
+ pointlist= numpy.array([(0,0), (15,98), (4,17), (67,32)])
sid = h5s.create_simple((100,100))
sid.select_elements(pointlist)
self.assertEqual(sid.get_select_elem_npoints(), len(pointlist))
- self.assertEqual(sid.get_select_elem_pointlist(), pointlist)
+ self.assert_(numpy.all(sid.get_select_elem_pointlist() == pointlist))
def test_get_blocks(self):
@@ -147,7 +148,7 @@ class TestH5S(unittest.TestCase):
self.assertEqual(sid.get_select_hyper_nblocks(), 2)
blocklist = sid.get_select_hyper_blocklist()
- self.assertEqual(blocklist, [( (0,0), (4,4) ), ( (50,60), (62,76) )])
+ self.assert_(numpy.all(blocklist == numpy.array([( (0,0), (4,4) ), ( (50,60), (62,76) )])))
def test_py(self):
diff --git a/h5py/tests/test_highlevel.py b/h5py/tests/test_highlevel.py
index 35a90cf..320cf62 100644
--- a/h5py/tests/test_highlevel.py
+++ b/h5py/tests/test_highlevel.py
@@ -206,7 +206,8 @@ class TestDataset(unittest.TestCase):
slices += [ s[0, ..., 49], s[...], s[..., 49], s[9,...] ]
slices += [ s[0:7:2,0:9:3,15:43:5], s[2:8:2,...] ]
slices += [ s[0], s[1], s[9], s[:] ] # Numpy convention
-
+ slices += [ numpy.random.random((10,10,50)) > 0.5 ] # Truth array
+
for dt in TYPES1:
srcarr = numpy.arange(10*10*50, dtype=dt).reshape(10,10,50)
@@ -221,7 +222,7 @@ class TestDataset(unittest.TestCase):
self.assertEqual(d.dtype, srcarr.dtype)
for argtpl in slices:
# Test read
- print " Checking read %.20s %s" % (dt, argtpl,)
+ print " Checking read %.20s %s" % (dt, argtpl if not isinstance(argtpl, numpy.ndarray) else 'ARRAY')
hresult = d[argtpl]
nresult = srcarr[argtpl]
self.assertEqual(hresult.shape, nresult.shape)
@@ -232,7 +233,7 @@ class TestDataset(unittest.TestCase):
d = Dataset(f, "NewDataset", data=srcarr)
for argtpl in slices:
# Test assignment
- print " Checking write %.20s %s" % (dt, argtpl,)
+ print " Checking write %.20s %s" % (dt, argtpl if not isinstance(argtpl, numpy.ndarray) else 'ARRAY')
srcarr[argtpl] = numpy.cos(srcarr[argtpl])
d[argtpl] = srcarr[argtpl]
self.assert_(numpy.all(d.value == srcarr))
diff --git a/h5py/utils.pxd b/h5py/utils.pxd
index 96bd74e..65f8660 100644
--- a/h5py/utils.pxd
+++ b/h5py/utils.pxd
@@ -37,4 +37,6 @@ cdef extern from "utils_low.h":
cdef int require_tuple(object tpl, int none_allowed, int size, char* name) except -1
cdef int require_list(object lst, int none_allowed, int size, char* name) except -1
cdef object pybool(long long val)
+cdef object create_numpy_hsize(int rank, hsize_t* dims)
+cdef object create_hsize_array(object arr)
diff --git a/h5py/utils.pyx b/h5py/utils.pyx
index 9c964ac..6d3901b 100644
--- a/h5py/utils.pyx
+++ b/h5py/utils.pyx
@@ -11,7 +11,58 @@
#-
from python cimport PyTuple_Check, PyList_Check, PyErr_SetString, Py_INCREF
+from numpy cimport import_array, NPY_UINT16, NPY_UINT32, NPY_UINT64, \
+ npy_intp, PyArray_SimpleNew, PyArray_ContiguousFromAny, \
+ PyArray_FROM_OTF, NPY_CONTIGUOUS, NPY_NOTSWAPPED, \
+ NPY_FORCECAST
+import_array()
+
+cdef object create_numpy_hsize(int rank, hsize_t* dims):
+ # Create an empty Numpy array which can hold HDF5 hsize_t entries
+
+ cdef int typecode
+ cdef npy_intp* dims_npy
+ cdef ndarray arr
+ cdef int i
+
+ if sizeof(hsize_t) == 2:
+ typecode = NPY_UINT16
+ elif sizeof(hsize_t) == 4:
+ typecode = NPY_UINT32
+ elif sizeof(hsize_t) == 8:
+ typecode = NPY_UINT64
+ else:
+ raise RuntimeError("Can't map hsize_t %d to Numpy typecode" % sizeof(hsize_t))
+
+ dims_npy = <npy_intp*>emalloc(sizeof(npy_intp)*rank)
+
+ try:
+ for i from 0<=i<rank:
+ dims_npy[i] = dims[i]
+ arr = PyArray_SimpleNew(rank, dims_npy, typecode)
+ finally:
+ efree(dims_npy)
+
+ return arr
+
+cdef object create_hsize_array(object arr):
+ # Create a NumPy array of hsize_t uints initialized to an existing array
+
+ cdef int typecode
+ cdef ndarray outarr
+
+ if sizeof(hsize_t) == 2:
+ typecode = NPY_UINT16
+ elif sizeof(hsize_t) == 4:
+ typecode = NPY_UINT32
+ elif sizeof(hsize_t) == 8:
+ typecode = NPY_UINT64
+ else:
+ raise RuntimeError("Can't map hsize_t %d to Numpy typecode" % sizeof(hsize_t))
+
+ return PyArray_FROM_OTF(arr, typecode, NPY_CONTIGUOUS | NPY_NOTSWAPPED | NPY_FORCECAST)
+
cdef int require_tuple(object tpl, int none_allowed, int size, char* name) except -1:
# Ensure that tpl is in fact a tuple, or None if none_allowed is nonzero.
# If size >= 0, also ensure that the length matches.
diff --git a/h5py/utils_hl.py b/h5py/utils_hl.py
index 83155c6..daef943 100644
--- a/h5py/utils_hl.py
+++ b/h5py/utils_hl.py
@@ -3,6 +3,7 @@
Utility functions for high-level modules.
"""
from __future__ import with_statement
+from h5py import h5s
from posixpath import basename, normpath
import numpy
@@ -15,56 +16,55 @@ def hbasename(name):
bname = '/'
return bname
-def slicer(shape, args):
- """ Parse a tuple containing Numpy-style extended slices.
- Shape should be a Numpy-style shape tuple.
+def slice_select(space, args):
+ """ Perform a selection on the given HDF5 dataspace, using a tuple
+ of Python extended slice objects. The dataspace may be scalar or
+ simple. The slice argument may be:
- Arguments may contain:
+ 0-tuple:
+ Entire dataspace selected (compatible with scalar)
- 1. Integer/long indices
- 2. Slice objects
- 3. Exactly one ellipsis object
- 4. Strings representing field names (zero or more, in any order)
+ 1-tuple:
+ 1. A single Ellipsis: entire dataspace selected
+ 2. A NumPy array: element-wise selection
+ 3. A single integer or slice (row-broadcasting)
- Return is a 4-tuple containing sub-tuples:
- (start, count, stride, names)
+ n-tuple:
+ 1. slice objects
+ 2. Ellipsis objects
+ 3. Integers
- start: tuple with starting indices
- count: how many elements to select along each axis
- stride: selection pitch for each axis
- names: field names (i.e. for compound types)
+ The return value is the appropriate memory dataspace to use.
"""
- rank = len(shape)
+ if len(args) == 0 or (len(args) == 1 and args[0] is Ellipsis):
+ space.select_all()
+ return space.copy()
- if not isinstance(args, tuple):
- args = (args,)
- args = list(args)
+ if len(args) == 1:
+ argval = args[0]
- slices = []
- names = []
+ if isinstance(argval, numpy.ndarray):
+ # Catch element-wise selection
+ indices = numpy.transpose(argval.nonzero())
+ space.select_elements(indices)
+ return h5s.create_simple((len(indices),))
- # Sort arguments
- for entry in args[:]:
- if isinstance(entry, str):
- names.append(entry)
- else:
- slices.append(entry)
+ # Single-index obj[0] access is always equivalent to obj[0,...].
+ # Pack it back up and send it to the hyperslab machinery
+ args = (argval, Ellipsis)
- # If only named fields are provided
- if len(slices) == 0:
- slices = [Ellipsis]
+ # Proceed to hyperslab selection
- # Hack to allow Numpy-style row indexing
- if len(slices) == 1 and slices[0] != Ellipsis:
- slices.append(Ellipsis)
+ shape = space.shape
+ rank = len(shape)
start = []
count = []
stride = []
# Expand integers and ellipsis arguments to slices
- for dim, arg in enumerate(slices):
+ for dim, arg in enumerate(args):
if isinstance(arg, int) or isinstance(arg, long):
if arg < 0:
@@ -110,7 +110,7 @@ def slicer(shape, args):
count.append(cc)
elif arg == Ellipsis:
- nslices = rank-(len(slices)-1)
+ nslices = rank-(len(args)-1)
if nslices <= 0:
continue
for x in range(nslices):
@@ -122,7 +122,8 @@ def slicer(shape, args):
else:
raise ValueError("Bad slice type %s" % repr(arg))
- return (tuple(start), tuple(count), tuple(stride), tuple(names))
+ space.select_hyperslab(tuple(start), tuple(count), tuple(stride))
+ return h5s.create_simple(tuple(count))
def strhdr(line, char='-'):
""" Print a line followed by an ASCII-art underline """
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/h5py.git
More information about the debian-science-commits
mailing list