gistfile1.txt

## Type aliases for standard C types, introduced for better readability

typealias C_int Int32
typealias C_unsigned Uint32
typealias C_char Uint8
typealias C_unsigned_long_long Uint64
typealias C_size_t Uint64

####################
## HDF5 interface ##
####################

## HDF5 types and constants

typealias C_hid_t  C_int
typealias C_herr_t C_int
typealias C_hsize_t C_unsigned_long_long
typealias C_tri_t C_int
typealias C_H5T_sign_t C_int
typealias C_H5T_class_t C_int

hdf5_symbols = {:H5E_DEFAULT    => convert(C_int, 0),
                :H5P_DEFAULT    => convert(C_int, 0),
                # file access modes
                :H5F_ACC_RDONLY => convert(C_unsigned, 0x00),
                :H5F_ACC_RDWR   => convert(C_unsigned, 0x01),
                :H5F_ACC_TRUNC  => convert(C_unsigned, 0x02),
                :H5F_ACC_EXCL   => convert(C_unsigned, 0x04),
                :H5F_ACC_DEBUG  => convert(C_unsigned, 0x08),
                :H5F_ACC_CREAT  => convert(C_unsigned, 0x10),
                # object types (C enum H5Itype_t)
                :H5I_FILE => 1,
                :H5I_GROUP => 2,
                :H5I_DATATYPE => 3,
                :H5I_DATASPACE => 4,
                :H5I_DATASET => 5,
                :H5I_ATTR => 6,
                :H5I_REFERENCE => 7,
                # type classes (C enum H5T_class_t)
                :H5T_INTEGER => 0,
                :H5T_FLOAT => 1,
                :H5T_TIME => 2,
                :H5T_STRING => 3,
                :H5T_BITFIELD => 4,
                :H5T_OPAQUE => 5,
                :H5T_COMPOUND => 6,
                :H5T_REFERENCE => 7,
                :H5T_ENUM => 8,
                :H5T_VLEN => 9,
                :H5T_ARRAY => 10,
                # Sign types (C enum H5T_sign_t)
                :H5T_SGN_NONE => 0, ## unsigned
                :H5T_SGN_2 => 1, ## 2's complement
                }

## Julia types corresponding to the HDF5 base types
hdf5_type_map = {(hdf5_symbols[:H5T_INTEGER],
                  hdf5_symbols[:H5T_SGN_2],
                  1) => Int8,
                 (hdf5_symbols[:H5T_INTEGER],
                  hdf5_symbols[:H5T_SGN_2],
                  2) => Int16,
                 (hdf5_symbols[:H5T_INTEGER],
                  hdf5_symbols[:H5T_SGN_2],
                  4) => Int32,
                 (hdf5_symbols[:H5T_INTEGER],
                  hdf5_symbols[:H5T_SGN_2],
                  8) => Int64,
                 (hdf5_symbols[:H5T_INTEGER],
                  hdf5_symbols[:H5T_SGN_NONE],
                  1) => Uint8,
                 (hdf5_symbols[:H5T_INTEGER],
                  hdf5_symbols[:H5T_SGN_NONE],
                  2) => Uint16,
                 (hdf5_symbols[:H5T_INTEGER],
                  hdf5_symbols[:H5T_SGN_NONE],
                  4) => Uint32,
                 (hdf5_symbols[:H5T_INTEGER],
                  hdf5_symbols[:H5T_SGN_NONE],
                  8) => Uint64,
                 (hdf5_symbols[:H5T_FLOAT],
                  nothing,
                  4) => Float32,
                 (hdf5_symbols[:H5T_FLOAT],
                  nothing,
                  8) => Float64,
                 }

## Load the HDF5 wrapper library and disable automatic error printing.
#hdf5lib = dlopen("hdf5_wrapper")
hdf5lib = dlopen("/opt/local/lib/libhdf5.dylib")
status = ccall(dlsym(hdf5lib, :H5Eset_auto2),
               C_herr_t,
               (C_hid_t, Ptr{Void}, Ptr{Void}),
               hdf5_symbols[:H5E_DEFAULT], C_NULL, C_NULL)
assert(status==0)


## HDF5 uses a plain integer to refer to each file, group, or
## dataset. These are wrapped into special types in order to allow
## method dispatch.

abstract HDF5Object

type HDF5File <: HDF5Object
   id::C_hid_t
   filename::String

   function HDF5File(id, filename)
       f = new(id, filename)
       finalizer(f, close)
       f
   end
end

type HDF5Group <: HDF5Object
   id::C_hid_t
   file::HDF5File

   function HDF5Group(id, file)
       g = new(id, file)
       finalizer(g, close)
       g
   end
end

type HDF5Dataset <: HDF5Object
   id::C_hid_t
   file::HDF5File

   function HDF5Dataset(id, file)
       ds = new(id, file)
       finalizer(ds, close)
       ds
   end
end

## HDF5-specific exception type

type HDF5Exception <: Exception
   msg::String
   object::HDF5Object
end

## HDF5 has a generic "close" function that closes groups,
## datasets, and named types. That allows a single close method
## for everything exept files.

function close(object::HDF5Object)
    status = ccall(dlsym(hdf5lib, :H5Oclose),
                   C_herr_t,
                   (C_hid_t,),
                   object.id)
    if status < 0
        error(HDF5Exception("could not close HDF5 object", object))
    end
end

function close(file::HDF5File)
    status = ccall(dlsym(hdf5lib, :H5Fclose),
                   C_herr_t,
                   (C_hid_t,),
                   file.id)
    if status < 0
        error(HDF5Exception("could not close HDF5 file", file))
    end
end

## h5open is made to resemble open() from io.j. It incorporates
## the functionality of H5Fopen and H5Fcreate.

function h5open(fname::String, rd::Bool, wr::Bool, cr::Bool, tr::Bool)
    if !rd
        error("HDF5 files have no write-only mode")
    end
    if cr != tr
        error("Truncation and creation are identical for HDF5 files")
    end
    if cr && !rd
        error("Can't create a read-only HDF5 file")
    end
    if cr
        fid = ccall(dlsym(hdf5lib, :H5Fcreate),
                    C_hid_t,
                    (Ptr{C_char}, C_unsigned, C_hid_t, C_hid_t),
                    fname,
                    hdf5_symbols[:H5F_ACC_TRUNC],
                    hdf5_symbols[:H5P_DEFAULT],
                    hdf5_symbols[:H5P_DEFAULT])
    else
        fid = ccall(dlsym(hdf5lib, :H5Fopen),
                    C_hid_t,
                    (Ptr{C_char}, C_unsigned, C_hid_t),
                    fname,
                    wr ? hdf5_symbols[:H5F_ACC_RDWR] :
                         hdf5_symbols[:H5F_ACC_RDONLY],
                    hdf5_symbols[:H5P_DEFAULT])
    end
    if fid < 0
        error("could not open HDF5 file ", fname)
    end
    HDF5File(fid, fname)
end

function h5open(fname::String, mode::String)
    mode == "r"  ? h5open(fname, true,  false, false, false) :
    mode == "r+" ? h5open(fname, true,  true , false, false) :
    mode == "w"  ? h5open(fname, true,  true , true , true)  :
    mode == "w+" ? h5open(fname, true,  true , true , true)  :
    error("invalid open mode: ", mode)
end

h5open(fname::String) = h5open(fname, true, false, false, false)

function root_group(h5file::HDF5File)
    gid = ccall(dlsym(hdf5lib, :H5Gopen2),
                C_hid_t,
                (C_hid_t, Ptr{C_char}, C_hid_t),
                h5file.id, "/", hdf5_symbols[:H5P_DEFAULT])
    if gid < 0
        error(HDF5Exception("could not open root group in HDF5 file ", h5file))
    end
    HDF5Group(gid, h5file)
end

function ref(group::HDF5Group, path::ASCIIString)
    obj_id = ccall(dlsym(hdf5lib, :H5Oopen),
                   C_hid_t,
                   (C_hid_t, Ptr{C_char}, C_hid_t),
                   group.id, path, hdf5_symbols[:H5P_DEFAULT])
    if obj_id < 0
        error(HDF5Exception("could not access path $path ", group))
    end
    obj_type = ccall(dlsym(hdf5lib, :H5Iget_type),
                     C_int,
                     (C_hid_t, ),
                     obj_id)
    obj_type == hdf5_symbols[:H5I_GROUP] ? HDF5Group(obj_id, group.file) :
    obj_type == hdf5_symbols[:H5I_DATATYPE] ? HDF5NamedType(obj_id, group.file) :
    obj_type == hdf5_symbols[:H5I_DATASET] ? HDF5Dataset(obj_id, group.file) :
    error(HDF5Exception("invalid object type for path $path ", group))
end

ref(file::HDF5File, path::ASCIIString) = ref(root_group(file), path)


function dataspace_id(ds::HDF5Dataset)
    space_id = ccall(dlsym(hdf5lib, :H5Dget_space),
                     C_hid_t,
                     (C_hid_t,),
                     ds.id)
    if space_id < 0
        error(HDF5Exception("could not open HDF5 dataspace", ds))
    end
    is_simple = ccall(dlsym(hdf5lib, :H5Sis_simple),
                      C_tri_t,
                      (C_hid_t,),
                      space_id)
    if is_simple <= 0   ## failure (negative) or false (0)
        close_dataspace(space_id)
        error(HDF5Exception("can't handle non-simple dataspace", ds))
    end
    space_id
end

function close_dataspace(id::C_hid_t)
    status = ccall(dlsym(hdf5lib, :H5Sclose),
                   C_herr_t,
                   (C_hid_t,),
                   id)
    if status < 0
        error(HDF5Exception("could not close HDF5 dataspace", id))
    end
end

function size(ds::HDF5Dataset)
    space_id = dataspace_id(ds)
    numdims = ccall(dlsym(hdf5lib, :H5Sget_simple_extent_ndims),
                    C_int,
                    (C_hid_t,),
                    space_id)
    dims = Array(C_hsize_t, numdims)
    numdims = ccall(dlsym(hdf5lib, :H5Sget_simple_extent_dims),
                    C_int,
                    (C_hid_t, Ptr{C_hsize_t}, Ptr{C_hsize_t}),
                    space_id, dims, C_NULL)
    close_dataspace(space_id)
    ntuple(numdims, i->dims[end-i+1])
end

size(ds::HDF5Dataset, d) = size(ds)[d]


######################################################################
######################################################################
######################################################################

h5f = h5open("test.h5")
ds = h5f["/universe/atoms"]

## This is the exact code of the function "size". Put at the
## top level, it works fine.
space_id = dataspace_id(ds)
numdims = ccall(dlsym(hdf5lib, :H5Sget_simple_extent_ndims),
                C_int,
                (C_hid_t,),
                space_id)
dims = Array(C_hsize_t, numdims)
numdims = ccall(dlsym(hdf5lib, :H5Sget_simple_extent_dims),
                C_int,
                (C_hid_t, Ptr{C_hsize_t}, Ptr{C_hsize_t}),
                space_id, dims, C_NULL)
close_dataspace(space_id)
println(ntuple(numdims, i->dims[end-i+1]))

## Here I call the function "size" and it crashes with a segmentatio
## fault.
println(size(ds))

## Here's what gdb says about the crash:
##
## Program received signal EXC_BAD_ACCESS, Could not access memory.
## Reason: KERN_INVALID_ADDRESS at address: 0x0000000000000001
## 0x0000000117c3c084 in H5Sget_simple_extent_dims ()

## 
close(h5f)