int main(int argc, char *argv[])
{
  float *two_theta = NULL;
  int *counts = NULL,  rank, i;
  hid_t fid, dataid, fapl;
  hsize_t *dim = NULL;
  hid_t datatype, dataspace, memdataspace;

  /*
   * Open file, thereby enforcing proper file close
   * semantics
   */
  fapl = H5Pcreate(H5P_FILE_ACCESS);
  H5Pset_fclose_degree(fapl,H5F_CLOSE_STRONG);
  fid = H5Fopen("NXfile.h5", H5F_ACC_RDONLY,fapl);  
  H5Pclose(fapl);

  /*
   * open and read the counts dataset
   */
  dataid = H5Dopen(fid,"/scan/data/counts");
  dataspace = H5Dget_space(dataid);
  rank = H5Sget_simple_extent_ndims(dataspace);
  dim = malloc(rank*sizeof(hsize_t));
  H5Sget_simple_extent_dims(dataspace, dim, NULL);
  counts = malloc(dim[0]*sizeof(int));
  memdataspace = H5Tcopy(H5T_NATIVE_INT32);
  H5Dread(dataid,memdataspace,H5S_ALL, H5S_ALL,H5P_DEFAULT, counts);
  H5Dclose(dataid);
  H5Sclose(dataspace);
  H5Tclose(memdataspace);

  /*
   * open and read the two_theta data set
   */
  dataid = H5Dopen(fid,"/scan/data/two_theta");
  dataspace = H5Dget_space(dataid);
  rank = H5Sget_simple_extent_ndims(dataspace);
  dim = malloc(rank*sizeof(hsize_t));
  H5Sget_simple_extent_dims(dataspace, dim, NULL);
  two_theta = malloc(dim[0]*sizeof(float));
  memdataspace = H5Tcopy(H5T_NATIVE_FLOAT);
  H5Dread(dataid,memdataspace,H5S_ALL, H5S_ALL,H5P_DEFAULT, two_theta);
  H5Dclose(dataid);
  H5Sclose(dataspace);
  H5Tclose(memdataspace);



  H5Fclose(fid);

  for(i = 0; i < dim[0]; i++){
    printf("%8.2f %10d\n", two_theta[i], counts[i]);
  }
  
}

int_f
nh5sselect_elements_c ( hid_t_f *space_id , int_f *op, size_t_f *nelements,  hsize_t_f *coord)
{
  int ret_value = -1;
  hid_t c_space_id;
  H5S_seloper_t c_op;
  herr_t  status;
  int rank;
  int i, j;
  hsize_t *c_coord;
  size_t c_nelements;
/*
  if (*op != H5S_SELECT_SET_F) return ret_value;
*/
  if (*op != H5S_SELECT_SET) return ret_value;
  c_op =  H5S_SELECT_SET;

  c_space_id = *space_id;
  rank = H5Sget_simple_extent_ndims(c_space_id);

  c_coord = malloc(sizeof(hsize_t)*rank*(*nelements));
  if(!c_coord) return ret_value;
  for (i=0; i< *nelements; i++) {
      for (j = 0; j < rank; j++) {
          c_coord[j+i*rank] = (hsize_t)coord[j + i*rank];
      }
  }

  c_nelements = *nelements;
  status = H5Sselect_elements(c_space_id, c_op, c_nelements, (const hsize_t **)c_coord);
  if ( status >= 0  ) ret_value = 0;
  HDfree(c_coord);
  return ret_value;
}

int_f
nh5sget_select_elem_pointlist_c( hid_t_f *space_id ,hsize_t_f * startpoint,
                                  hsize_t_f * numpoints, hsize_t_f * buf)
{
  int ret_value = -1;
  hid_t c_space_id;
  hsize_t c_num_points;
  hsize_t c_startpoint,* c_buf;
  int i, rank;

  c_space_id = *space_id;
  c_num_points = (hsize_t)* numpoints;

  rank = H5Sget_simple_extent_ndims(c_space_id);
  if (rank < 0 ) return ret_value;

  c_startpoint = (hsize_t)*startpoint;
  c_buf = (hsize_t*)malloc(sizeof(hsize_t)*(size_t)(c_num_points*rank));
  if (!c_buf) return ret_value;
  ret_value = H5Sget_select_elem_pointlist(c_space_id, c_startpoint,
                                            c_num_points, c_buf);
  for (i = c_num_points*rank-1; i >= 0; i--) {
      buf[i] = (hsize_t_f)(c_buf[i]+1);
  }

  if (ret_value  >= 0  ) ret_value = 0;

  HDfree(c_buf);

  return ret_value;
}

int_f
nh5sget_select_hyper_blocklist_c( hid_t_f *space_id ,hsize_t_f * startblock,
                                  hsize_t_f * num_blocks, hsize_t_f * buf)
{
  int ret_value = -1;
  hid_t c_space_id;
  hsize_t c_num_blocks;

  hsize_t i;
  int rank;
  hsize_t c_startblock, *c_buf;

  c_space_id = *space_id;
  c_num_blocks = * num_blocks;

  rank = H5Sget_simple_extent_ndims(c_space_id);
  if (rank < 0 ) return ret_value;
  c_startblock = (hsize_t)*startblock;

  c_buf = (hsize_t*)malloc(sizeof(hsize_t)*(size_t)(c_num_blocks*2*rank));
  if (!c_buf) return ret_value;

  ret_value = H5Sget_select_hyper_blocklist(c_space_id, c_startblock,
                                            c_num_blocks, c_buf);
  for(i = 0; i < c_num_blocks*2*rank; i++)
  {
      buf[i] = (hsize_t_f)c_buf[i] +1;
  }
  HDfree(c_buf);
  if (ret_value  >= 0  ) ret_value = 0;
  return ret_value;
}

//-----------------------------------------------------------------------------
int MGL_EXPORT mgl_datac_read_hdf(HADT d,const char *fname,const char *data)
{
	hid_t hf,hd,hs;
	hsize_t dims[4];
	long rank, res = H5Fis_hdf5(fname);
	if(res<=0)	{	return false;	}
	hf = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
	if(hf<0)	return false;
	hd = H5Dopen(hf,data);
	if(hd<0)	return false;
	hs = H5Dget_space(hd);
	rank = H5Sget_simple_extent_ndims(hs);
	if(rank>0 && rank<=3)
	{
		H5Sget_simple_extent_dims(hs,dims,0);
		if(rank==2)			{	dims[2]=dims[0];	dims[0]=dims[1]=1;	}
		else if(rank==3)	{	dims[2]=dims[1];	dims[1]=dims[0];	dims[0]=1;	}
//		else if(rank>3)		continue;
		mgl_datac_create(d,dims[2],dims[1],dims[0]);
#if MGL_USE_DOUBLE
		H5Dread(hd, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, d->a);
#else
		H5Dread(hd, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, d->a);
#endif
	}
	H5Sclose(hs);	H5Dclose(hd);	H5Fclose(hf);	return true;
}

int read_vector(hid_t loc_id,const char *name,T *value,int len)
{
  hid_t space_id,dset_id,dti_file;
  int hdferr,ndims;
  hsize_t dims[2], maxdims[2];

  printf("Attempting to read vector '%s'\n",name);

  dset_id= H5Dopen2(loc_id,name,H5P_DEFAULT);

  space_id= H5Dget_space(dset_id);
  ndims= H5Sget_simple_extent_ndims(space_id);
  if(ndims!=1) return -1;
  H5Sget_simple_extent_dims(space_id,dims,maxdims);

  if(dims[0]>len) return -1;

  hdferr= H5Dread(dset_id,TypeTraits<T>::hdf5Type(),H5S_ALL,H5S_ALL,H5P_DEFAULT,value);
  if(hdferr<0) return -1;

  H5Sclose(space_id);
  H5Dclose(dset_id);

  return dims[0];
}

cv::Mat hdf5file::load(const std::string& name)
{
	hdf5dataset hdataset(*this, name);
	hid_t dataset_id = hdataset.handle();

	hid_t space_id = H5Dget_space(dataset_id);
	hid_t type_id = hdataset.type();
	//hid_t type_id = H5Tget_native_type(type_id_nnative, H5T_DIR_DEFAULT);

	hsize_t rank = H5Sget_simple_extent_ndims(space_id);

	std::vector<hsize_t> dims_out(rank);
	H5Sget_simple_extent_dims(space_id, &(dims_out[0]), NULL);

	//due to different types, the dimension sizes are copied in a vector of different type
	std::vector<int> dims_out_i(rank);
	std::copy(dims_out.begin(), dims_out.end(), dims_out_i.begin());

	int channels = 1;

	if(rank == 3 && dims_out[rank-1] <= 4)
	{
		channels = dims_out[rank-1];
		--rank;
	}

	cv::Mat dataset = cv::Mat(rank, &(dims_out_i[0]), CV_MAKETYPE(h5types::typeid_to_ocv(type_id), channels));

	hdataset.read(dataset.data);

	H5Sclose(space_id);

	return dataset;
}

int_f
h5dset_extent_c ( hid_t_f *dset_id , hsize_t_f *dims)
/******/
{
  hid_t c_space_id;
  hsize_t c_dims[H5S_MAX_RANK];
  int rank;
  int i;
  int status;
  int ret_value = -1;

  if((c_space_id = H5Dget_space((hid_t)*dset_id)) < 0) return ret_value;

  rank = H5Sget_simple_extent_ndims(c_space_id);
  H5Sclose(c_space_id);
  if(rank < 0 ) return ret_value;


  /*
   * Reverse dimensions due to C-FORTRAN storage order.
   */
  for(i = 0; i < rank; i++)
      c_dims[i] = (hsize_t)dims[rank - i - 1];

  status = H5Dset_extent((hid_t)*dset_id, c_dims);

  if(status >= 0)
      ret_value = 0;
  return ret_value;
}

// JRC: This fat interface may not scale?  What about
// scalar attributes?
herr_t VsH5Attribute::getDoubleVectorValue(std::vector<double>* dvals) {
  herr_t err = 0;
  size_t npoints;
  hid_t atype = H5Aget_type(getId());
  H5T_class_t type = H5Tget_class(atype);
  hid_t aspace = H5Aget_space(getId());
  size_t rank = H5Sget_simple_extent_ndims(aspace);

  if (type != H5T_FLOAT) {
    VsLog::warningLog() <<"VsH5Attribute::getDoubleVectorValue() - Requested attribute " <<getShortName()
         <<" is not a floating point vector." <<std::endl;
    dvals->resize(0);
    return -1;
  }
  
  if (rank == 0) {
    dvals->resize(1);
    double v;
    err = H5Aread(getId(), H5T_NATIVE_DOUBLE, &v);
    (*dvals)[0] = v;
    return err;
  }
  
  // rank>0
  npoints = H5Sget_simple_extent_npoints(aspace);
  double* v = new double[npoints];
  err = H5Aread(getId(), H5T_NATIVE_DOUBLE, v);
  dvals->resize(npoints);
  for (size_t i = 0; i<npoints; ++i) {
    (*dvals)[i] = v[i];
  }
  delete [] v;
  
  return err;
}

int luaC_h5_get_ndims(lua_State *L)
{
  const char *dsetnm = luaL_checkstring(L, 1);

  if (PresentFile < 0) {
    luaL_error(L, "no open file");
  }
  hid_t dset = H5Dopen(PresentFile, dsetnm, H5P_DEFAULT);
  if (dset < 0) {
    luaL_error(L, "no data set named %s", dsetnm);
  }
  hid_t fspc = H5Dget_space(dset);
  size_t Nd  = H5Sget_simple_extent_ndims(fspc);
  hsize_t *dims = (hsize_t*) malloc(Nd*sizeof(hsize_t));
  int *dims_int =     (int*) malloc(Nd*sizeof(int));

  H5Sget_simple_extent_dims(fspc, dims, NULL);

  int i;
  for (i=0; i<Nd; ++i) {
    dims_int[i] = dims[i];
  }
  _pusharray_i(L, dims_int, Nd);

  H5Sclose(fspc);
  H5Dclose(dset);
  free(dims);
  free(dims_int);

  return 1;
}

void h5_value_doubles(hid_t file_id, char *group, char *name, double **values,
	int *numValues)
{
  int kk;
  hid_t group_id = H5Gopen(file_id, group, H5P_DEFAULT);
  hid_t data_id = H5Dopen(group_id, name, H5P_DEFAULT);
  hid_t data_space = H5Dget_space(data_id);
  int rank = H5Sget_simple_extent_ndims(data_space);
  hsize_t dims[H5S_MAX_RANK], maxdim[H5S_MAX_RANK];
  H5Sget_simple_extent_dims(data_space, dims, maxdim);
  hid_t data_type = H5Dget_type(data_id);
  H5T_class_t data_class = H5Tget_native_type(data_type, H5T_DIR_DEFAULT);
  size_t data_size = H5Tget_size(data_class);  
  hsize_t elements = 1;
  for (kk=0; kk<rank; kk++)
    elements *= dims[kk];
  void *buf = (void *) MALLOC((size_t)(elements*data_size));
  H5Dread(data_id, data_class, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
  *values = buf;
  *numValues = elements;
  H5Tclose(data_type);
  H5Tclose(data_class);
  H5Sclose(data_space);
  H5Dclose(data_id);
  H5Gclose(group_id); 
}

int_f
nh5sget_simple_extent_dims_c ( hid_t_f *space_id , hsize_t_f *dims, hsize_t_f *maxdims)
{
  int ret_value = -1;
  hid_t c_space_id;
  hsize_t *c_dims;
  hsize_t *c_maxdims;
  int status;
  int rank;
  int i;

  c_space_id = *space_id;
  rank = H5Sget_simple_extent_ndims(c_space_id);
  if (rank < 0) return ret_value;

  c_dims = malloc(sizeof(hsize_t)*rank);
  if (!c_dims) return ret_value;

  c_maxdims = malloc(sizeof(hsize_t)*rank);
  if (!c_maxdims) return ret_value;

  status = H5Sget_simple_extent_dims(c_space_id, c_dims, c_maxdims);
  /*
   * Reverse dimensions due to C-FORTRAN storage order.
   */
  for (i=0; i < rank; i++) {
      dims[rank - i - 1] = (hsize_t_f)c_dims[i];
      maxdims[rank - i - 1] = (hsize_t_f)c_maxdims[i];
  }

  if ( status >= 0  ) ret_value = rank;
  HDfree(c_dims);
  HDfree(c_maxdims);
  return ret_value;
}

int_f
nh5sget_select_bounds_c( hid_t_f *space_id , hsize_t_f * start, hsize_t_f * end)
{
  int ret_value = -1;
  hid_t c_space_id;
  hsize_t* c_start, *c_end;
  int i, rank;

  c_space_id = *space_id;
  rank = H5Sget_simple_extent_ndims(c_space_id);
  if (rank < 0 ) return ret_value;

  c_start =(hsize_t*) malloc(sizeof(hsize_t)*rank);
  if (!c_start) return ret_value;

  c_end = (hsize_t*)malloc(sizeof(hsize_t)*rank);
  if(!c_end) return ret_value;

  ret_value = H5Sget_select_bounds(c_space_id, c_start, c_end);
  for(i = 0; i < rank; i++)
  {
    start[i] = (hsize_t_f)(c_start[rank-i-1]+1);
    end[i] = (hsize_t_f)(c_end[rank-i-1]+1);
  }
  if (ret_value  >= 0  ) ret_value = 0;

  HDfree(c_start);
  HDfree(c_end);

  return ret_value;
}

//hyperslab write
void hdf5IoDataModel::writeHyperslab(const QString &dataset_name, DataType type, quint64 *offset, quint64 *stride, quint64 *count, quint64 *block, quint64 *values_shape, void *values)
{
    hid_t dataset_id = d->datasetId(dataset_name);

    //the selection within the file dataset's dataspace
    hid_t file_dataspace = H5Dget_space(dataset_id);
    if(H5Sselect_hyperslab(file_dataspace, H5S_SELECT_SET, offset, stride, count, block)<0) {
        dtkError() << "ERROR selecting hyperslab" << dataset_name;
    }

    //set the dimensions of values. memory dataspace and the selection within it
    hid_t values_dataspace = H5Screate_simple(H5Sget_simple_extent_ndims(file_dataspace), values_shape, NULL);

    switch(type) {
    case dtkIoDataModel::Int:
        // TODO put d->prop_list_id instead of H5P_DEFAULT ????????
        d->status = H5Dwrite(dataset_id, H5T_NATIVE_INT, values_dataspace, file_dataspace, H5P_DEFAULT, values);
        break;
    case dtkIoDataModel::LongLongInt:
        d->status = H5Dwrite(dataset_id, H5T_NATIVE_LLONG, values_dataspace, file_dataspace, H5P_DEFAULT, values);
        break;
    case dtkIoDataModel::Double:
        d->status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, values_dataspace, file_dataspace, H5P_DEFAULT, values);
        break;
    default:
        dtkError() << "write method: Datatype not supported";
    };
    if(d->status<0) {
        dtkError() << "error writing hyperslab" << dataset_name;
    }

    H5Sclose(file_dataspace);
    H5Sclose(values_dataspace);

}

//-*****************************************************************************
// Get the dimensions directly off of the dataspace on the dataset
// This isn't suitable for string and wstring
void
ReadDataSetDimensions( hid_t iParent,
                       const std::string &iName,
                       hsize_t iExtent,
                       Dimensions &oDims )
{
    // Open the data set.
    hid_t dsetId = H5Dopen( iParent, iName.c_str(), H5P_DEFAULT );
    ABCA_ASSERT( dsetId >= 0, "Cannot open dataset: " << iName );
    DsetCloser dsetCloser( dsetId );

    // Read the data space.
    hid_t dspaceId = H5Dget_space( dsetId );
    ABCA_ASSERT( dspaceId >= 0, "Could not get dataspace for dataSet: "
                 << iName );
    DspaceCloser dspaceCloser( dspaceId );

    H5S_class_t dspaceClass = H5Sget_simple_extent_type( dspaceId );
    if ( dspaceClass == H5S_SIMPLE )
    {
        // Get the dimensions
        int rank = H5Sget_simple_extent_ndims( dspaceId );
        ABCA_ASSERT( rank == 1, "H5Sget_simple_extent_ndims() must be 1." );

        hsize_t hdim = 0;
        rank = H5Sget_simple_extent_dims( dspaceId, &hdim, NULL );
        oDims.setRank(1);
        oDims[0] = hdim / iExtent;
    }
    else
    {
        oDims.setRank(1);
        oDims[0] = 0;
    }
}

inline size_t DataSpace::getNumberDimensions() const{
    const int ndim = H5Sget_simple_extent_ndims(_hid);
    if(ndim < 0){
        HDF5ErrMapper::ToException<DataSetException>("Unable to get dataspace number of dimensions");
    }
    return size_t(ndim);
}

/* read attribute value */
char* H5Dataset_read_attr_value(hid_t aid)
{
    hid_t asid=-1, atid=-1;
    int i, rank;
    hsize_t *dims;
    size_t size=1;
    char *attr_buf=NULL;
  
    asid = H5Aget_space(aid);
    atid = H5Aget_type(aid);
    rank = H5Sget_simple_extent_ndims(asid);

    if (rank > 0) {
        dims = (hsize_t *)malloc(rank * sizeof(hsize_t));
        H5Sget_simple_extent_dims(asid, dims, NULL);
        for (i=0; i<rank; i++) {
            size *= (size_t)dims[i];
            free(dims);
        }
        size *= H5Tget_size(atid);
        attr_buf = (char *)malloc(size);

        if (H5Aread( aid, atid, attr_buf) < 0) {
            free(attr_buf);
            attr_buf = NULL;
        }

    }
    
    if( atid > 0) H5Tclose(atid);
    if (asid > 0) H5Sclose(asid);
 
    return attr_buf;
}

bool Hdf5Dataset::getSphereRI(MapVecDouble &mvec)
{
  hsize_t dims_out[2], count[2], offset[2], dimsm[2];
  hid_t dataspace = H5Dget_space(this->sphere_dataset_); // dataspace handle
  int rank = H5Sget_simple_extent_ndims(dataspace);
  herr_t status_n = H5Sget_simple_extent_dims(dataspace, dims_out, NULL);
  herr_t status;
  offset[0] = 0;
  offset[1] = 0;
  count[0] = dims_out[0];
  count[1] = 4;
  double data_out[count[0]][count[1]];
  status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL, count, NULL);
  dimsm[0] = count[0];
  dimsm[1] = count[1];
  hid_t memspace;
  memspace = H5Screate_simple(RANK_OUT, dimsm, NULL);
  status = H5Dread(this->sphere_dataset_, H5T_NATIVE_DOUBLE, memspace, dataspace, H5P_DEFAULT, data_out);
  for (int i = 0; i < count[0]; i++)
  {
    std::vector< double > sphere_center(3);
    double ri;
    for (int j = 0; j < 3; j++)
    {
      sphere_center[j] = data_out[i][j];
    }
    for (int k = 3; k < 4; k++)
    {
      ri = data_out[i][k];
    }
    mvec.insert(std::pair< std::vector< double >, double >(sphere_center, ri));
   }
  return 0;
}

static void
ufo_hdf5_reader_open (UfoReader *reader,
                      const gchar *filename)
{
    UfoHdf5ReaderPrivate *priv;
    gchar *h5_filename;
    gchar *h5_dataset;
    gchar **components;

    priv = UFO_HDF5_READER_GET_PRIVATE (reader);
    components = g_strsplit (filename, ":", 2);

    if (components[1] == NULL) {
        g_warning ("hdf5: must specify dataset name after color");
        return;
    }

    h5_filename = components[0];
    h5_dataset = components[1];

    priv->file_id = H5Fopen (h5_filename, H5F_ACC_RDWR, H5P_DEFAULT);
    priv->dataset_id = H5Dopen (priv->file_id, h5_dataset, H5P_DEFAULT);
    priv->src_dataspace_id = H5Dget_space (priv->dataset_id);
    priv->n_dims = H5Sget_simple_extent_ndims (priv->src_dataspace_id);

    if (priv->n_dims > 3)
        g_error ("read:hdf5: no support for four-dimensional data");

    H5Sget_simple_extent_dims (priv->src_dataspace_id, priv->dims, NULL);

    priv->current = 0;
    g_strfreev (components);
}

ls2_hdf5_read_anchors(hid_t file, vector2** anchors, size_t *no_anchors)
{
    hid_t dataset, dataspace, memspace;
    hsize_t dims[2];
    int rank;

    dataset = H5Dopen2(file, "/Anchors", H5P_DEFAULT);
    dataspace = H5Dget_space(dataset);
    rank = H5Sget_simple_extent_ndims(dataspace);
    if (rank != 2) {
        fprintf(stderr, "/Anchors wrong rank %d\n", (int) rank);
        return -1;
    }
    H5Sget_simple_extent_dims(dataspace, dims, NULL);
    if (dims[1] != 2) {
        fprintf(stderr, "/Anchors dimension (%d, %d)\n",
                (int) dims[0], (int) dims[1]);
        return -1;
    }
    *no_anchors = (size_t) dims[0];
    *anchors = (vector2 *) calloc(*no_anchors, sizeof(vector2));
    memspace = H5Screate_simple(2, dims, NULL);
    H5Dread(dataset, H5T_NATIVE_FLOAT, memspace, dataspace, H5P_DEFAULT,
            *anchors);
    H5Dclose(dataset);
    H5Sclose(dataspace);
    H5Sclose(memspace);

    return 0;
}