arrayh5.c

/* Copyright (c) 1999-2008 Massachusetts Institute of Technology
 *
 * 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.
 */

#include <stdlib.h>
#include <string.h>

/* don't use new HDF5 1.8 API (which isn't even fully documented yet, grrr) */
#define H5_USE_16_API 1

#include <hdf5.h>

#include "arrayh5.h"

#define CHECK(cond, msg) { if (!(cond)) { fprintf(stderr, "arrayh5 error: %s\n", msg); exit(EXIT_FAILURE); } }

#define CHK_MALLOC(p, t, n) CHECK(p = (t *) malloc(sizeof(t) * (n)), "out of memory")

/* Normally, HDF5 prints out all sorts of error messages, e.g. if a dataset
   can't be found, in addition to returning an error code.  The following
   macro can be wrapped around code to temporarily suppress error messages. */

#define SUPPRESS_HDF5_ERRORS(statements) { \
     H5E_auto_t xxxxx_err_func; \
     void *xxxxx_err_func_data; \
     H5Eget_auto(&xxxxx_err_func, &xxxxx_err_func_data); \
     H5Eset_auto(NULL, NULL); \
     { statements; } \
     H5Eset_auto(xxxxx_err_func, xxxxx_err_func_data); \
}

arrayh5 arrayh5_create_withdata(int rank, const int *dims, double *data)
{
     arrayh5 a;
     int i;

     CHECK(rank >= 0, "non-positive rank");
     a.rank = rank;

     CHK_MALLOC(a.dims, int, rank);

     a.N = 1;
     for (i = 0; i < rank; ++i) {
	  a.dims[i] = dims[i];
	  a.N *= dims[i];
     }

     if (data)
	  a.data = data;
     else {
	  CHK_MALLOC(a.data, double, a.N);
     }
     return a;
}

arrayh5 arrayh5_create(int rank, const int *dims)
{
     return arrayh5_create_withdata(rank, dims, NULL);
}

arrayh5 arrayh5_clone(arrayh5 a)
{
     arrayh5 b = arrayh5_create(a.rank, a.dims);
     if (a.data) memcpy(b.data, a.data, sizeof(double) * a.N);
     return b;
}

void arrayh5_destroy(arrayh5 a)
{
     free(a.dims);
     free(a.data);
}

int arrayh5_conformant(arrayh5 a, arrayh5 b)
{
     int i;

     if (a.rank != b.rank)
	  return 0;
     for (i = 0; i < a.rank; ++i)
	  if (a.dims[i] != b.dims[i])
	       return 0;
     return 1;
}

static void rtranspose(int curdim, int rank, const int *dims,
		       int curindex, int curindex_t,
		       const double *data, double *data_t)
{
     int prod_before = 1, prod_after = 1;
     int i;

     if (rank == 0) {
	  *data_t = *data;
	  return;
     }

     for (i = 0; i < curdim; ++i)
	  prod_before *= dims[i];
     for (i = curdim + 1; i < rank; ++i)
	  prod_after *= dims[i];

     if (curdim == rank - 1) {
	  for (i = 0; i < dims[curdim]; ++i)
	       data_t[curindex_t + i * prod_before] = data[curindex + i];
     }
     else {
	  for (i = 0; i < dims[curdim]; ++i)
	       rtranspose(curdim + 1, rank, dims,
			  curindex + i * prod_after,
			  curindex_t + i * prod_before,
			  data, data_t);
     }
}

void arrayh5_transpose(arrayh5 *a)
{
     double *data_t;
     int i;

     CHK_MALLOC(data_t, double, a->N);
     rtranspose(0, a->rank, a->dims, 0, 0, a->data, data_t);
     free(a->data);
     a->data = data_t;

     for (i = 0; i < a->rank - 1 - i; ++i) {
	  int dummy = a->dims[i];
	  a->dims[i] = a->dims[a->rank - 1 - i];
	  a->dims[a->rank - 1 - i] = dummy;
     }
}

void arrayh5_getrange(arrayh5 a, double *min, double *max)
{
     int i;

     CHECK(a.N > 0, "no elements in array");
     *min = *max = a.data[0];
     for (i = 1; i < a.N; ++i) {
	  if (a.data[i] < *min)
	       *min = a.data[i];
	  if (a.data[i] > *max)
	       *max = a.data[i];
     }
}

static herr_t find_dataset(hid_t group_id, const char *name, void *d)
{
     char **dname = (char **) d;
     H5G_stat_t info;

     H5Gget_objinfo(group_id, name, 1, &info);
     if (info.type == H5G_DATASET) {
	  CHK_MALLOC(*dname, char, strlen(name) + 1);
	  strcpy(*dname, name);
	  return 1;
     }
     return 0;
}

typedef enum { NO_ERROR = 0, OPEN_FAILED, NO_DATA, READ_FAILED, SLICE_FAILED,
	     INVALID_SLICE, INVALID_RANK, OPEN_DATA_FAILED } arrayh5_err;

const char arrayh5_read_strerror[][100] = {
     "no error",
     "error opening HD5 file",
     "couldn't find data set in HDF5 file",
     "error reading data from HDF5",
     "error reading data slice from HDF5",
     "invalid slice of HDF5 data",
     "non-positive rank in HDF file",
     "error opening data set in HDF file",
};

int arrayh5_read(arrayh5 *a, const char *fname, const char *datapath,
		 char **dataname,
		 int nslicedims_, const int *slicedim_, const int *islice_,
		 const int *center_slice)
{
     hid_t file_id = -1, data_id = -1, space_id = -1;
     char *dname = NULL;
     int err = NO_ERROR;
     hsize_t i, rank, *dims_copy, *maxdims, *slicedim = 0;
     int *islice = 0;
     int *dims = 0;
     hsize_t nslicedims = (hsize_t) nslicedims_;

     CHECK(a, "NULL array passed to arrayh5_read");
     a->dims = NULL;
     a->data = NULL;

     file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
     if (file_id < 0) {
	  err = OPEN_FAILED;
	  goto done;
     }

     if (datapath && datapath[0]) {
	  CHK_MALLOC(dname, char, strlen(datapath) + 1);
	  strcpy(dname, datapath);
     }
     else {
	  if (H5Giterate(file_id, "/", NULL, find_dataset, &dname) <= 0) {
	       err = NO_DATA;
	       goto done;
	  }
     }

     data_id = H5Dopen(file_id, dname);
     if (data_id < 0) {
	  err = OPEN_DATA_FAILED;
	  goto done;
     }

     space_id = H5Dget_space(data_id);
     rank = H5Sget_simple_extent_ndims(space_id);
     if (rank <= 0) {
	  err = INVALID_RANK;
	  goto done;
     }

     CHK_MALLOC(dims, int, rank);
     CHK_MALLOC(dims_copy, hsize_t, rank);
     CHK_MALLOC(maxdims, hsize_t, rank);

     H5Sget_simple_extent_dims(space_id, dims_copy, maxdims);
     for (i = 0; i < rank; ++i)
	  dims[i] = dims_copy[i];

     free(maxdims);
     free(dims_copy);

     for (i = 0; i < nslicedims && slicedim_[i] == NO_SLICE_DIM; ++i)
	  ;

     if (i == nslicedims) { /* no slices */
	  *a = arrayh5_create(rank, dims);

	  if (H5Dread(data_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
		      H5P_DEFAULT, (void *) a->data) < 0) {
	       err = READ_FAILED;
	       goto done;
	  }
     }
     else if (nslicedims > 0) {
	  int j, rank2 = rank;
	  hsize_t *start;
	  hsize_t *count;
	  hid_t mem_space_id;
	  herr_t readerr;

	  CHK_MALLOC(slicedim, hsize_t, nslicedims);
	  CHK_MALLOC(islice, int, nslicedims);

	  for (i = j = 0; i < nslicedims; ++i)
	       if (slicedim_[i] != NO_SLICE_DIM) {
		    if (slicedim_[i] == LAST_SLICE_DIM)
			 slicedim[j] = rank - 1;
		    else
			 slicedim[j] = (hsize_t) slicedim_[i];
		    if (slicedim[j] >= rank) {
			 err = INVALID_SLICE;
			 goto done;
		    }
		    islice[j] = islice_[i];
		    if (center_slice[i])
			 islice[j] += dims[slicedim[j]] / 2;
		    if (islice[j] < 0 || islice[j] >= dims[slicedim[j]]) {
			 err = INVALID_SLICE;
			 goto done;
		    }
		    j++;
	       }
	  nslicedims = j;

	  CHK_MALLOC(start, hsize_t, rank);
	  CHK_MALLOC(count, hsize_t, rank);

	  for (i = 0; i < rank; ++i) {
	       count[i] = dims[i];
	       start[i] = 0;
	  }
	  for (i = 0; i < nslicedims; ++i) {
	       start[slicedim[i]] = islice[i];
	       count[slicedim[i]] = 1;
	  }

	  H5Sselect_hyperslab(space_id, H5S_SELECT_SET,
			      start, NULL, count, NULL);

	  for (i = j = 0; i < rank; ++i)
	       if (count[i] > 1)
		    dims[j++] = count[i];
	  rank2 = j;

	  *a = arrayh5_create(rank2, dims);

	  mem_space_id = H5Screate_simple(rank, count, NULL);
	  H5Sselect_all(mem_space_id);

	  readerr = H5Dread(data_id, H5T_NATIVE_DOUBLE,
			    mem_space_id, space_id,
			    H5P_DEFAULT, (void *) a->data);

	  H5Sclose(mem_space_id);
	  free(count);
	  free(start);

	  if (readerr < 0) {
	       err = SLICE_FAILED;
	       goto done;
	  }
     }
     else {
	  err = INVALID_SLICE;
	  goto done;
     }

 done:
     if (err != NO_ERROR)
	  arrayh5_destroy(*a);
     free(islice);
     free(slicedim);
     free(dims);
     if (space_id >= 0)
	  H5Sclose(space_id);
     if (data_id >= 0)
	  H5Dclose(data_id);
     if (dataname)
	  *dataname = dname;
     else
	  free(dname);
     if (file_id >= 0)
	  H5Fclose(file_id);

     return err;
}

static int dataset_exists(hid_t id, const char *name)
{
     hid_t data_id;
     SUPPRESS_HDF5_ERRORS(data_id = H5Dopen(id, name));
     if (data_id >= 0)
          H5Dclose(data_id);
     return (data_id >= 0);
}

void arrayh5_write(arrayh5 a, char *filename, char *dataname,
		   short append_data)
{
     int i;
     hid_t file_id, space_id, type_id, data_id;
     hsize_t *dims_copy;

     if (append_data)
	  file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
     else
	  file_id = H5Fcreate(filename, H5F_ACC_TRUNC,
			      H5P_DEFAULT, H5P_DEFAULT);
     CHECK(file_id >= 0, "error opening HDF5 output file");

     if (dataset_exists(file_id, dataname))
	  H5Gunlink(file_id, dataname);  /* delete it */

     CHECK(a.rank > 0, "non-positive rank");
     CHK_MALLOC(dims_copy, hsize_t, a.rank);
     for (i = 0; i < a.rank; ++i)
	  dims_copy[i] = a.dims[i];
     space_id = H5Screate_simple(a.rank, dims_copy, NULL);
     free(dims_copy);

     type_id = H5T_NATIVE_DOUBLE;
     data_id = H5Dcreate(file_id, dataname, type_id, space_id, H5P_DEFAULT);
     H5Sclose(space_id);

     H5Dwrite(data_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, a.data);

     H5Dclose(data_id);
     H5Fclose(file_id);
}

int arrayh5_read_rank(const char *fname, const char *datapath, int *rank)
{
     hid_t file_id = -1, data_id = -1, space_id = -1;
     char *dname = NULL;
     int err = NO_ERROR;

     file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
     if (file_id < 0) {
	  err = OPEN_FAILED;
	  goto done;
     }

     if (datapath && datapath[0]) {
	  CHK_MALLOC(dname, char, strlen(datapath) + 1);
	  strcpy(dname, datapath);
     }
     else {
	  if (H5Giterate(file_id, "/", NULL, find_dataset, &dname) <= 0) {
	       err = NO_DATA;
	       goto done;
	  }
     }

     data_id = H5Dopen(file_id, dname);
     if (data_id < 0) {
	  err = OPEN_DATA_FAILED;
	  goto done;
     }

     space_id = H5Dget_space(data_id);
     *rank = H5Sget_simple_extent_ndims(space_id);

 done:
     if (space_id >= 0)
	  H5Sclose(space_id);
     if (data_id >= 0)
	  H5Dclose(data_id);
     free(dname);
     if (file_id >= 0)
	  H5Fclose(file_id);

     return err;
}