[C,C++] HDF入門記録，読み書きの基本操作

環境

MacOS, Ubuntu
gcc-9.x
hdf5 1.12.0

目的

wikipedia によると，HDFは，

Hierarchical Data Format（階層的データ形式、略称：HDF）は、大量のデータを格納および構造化するために設計された一連のファイル形式（HDF4、HDF5）。

です．

一つのファイルの中に，グループ(ディレクトリのようなもの)とデータセット(ファイルのようなもの)を再帰的に作成することができ，linkや部分圧縮などの高度な機能も持ちます．

また，HDFはマルチプラットフォームであるとともに，多くの言語でサポートされており，例えばPythonだとライブラリh5pyが開発されています．

このHDF5をC,C++から利用したかったのですが，いざ公式のドキュメントを見ると覚えることが多く，初心者には利用までのハードルが高いと感じました．

そこでひとまず，基本機能であるファイルオープン・クローズ，グループの作成，データセットの作成，データセットの(部分)更新機能まで，自分なりに勉強した記録をまとめておくことにします．．

やること

C,C++からHDFを扱う方法を練習し，手軽に使えるようにまとめます．

公式のexampleはこちら．

といっても，C++はmpi対応していないなどの問題があるらしいので，Cのapiを勉強してラップしていくことにします．

取扱い

利用を始める

プログラムの中で以下のヘッダーをincludeします．

# include <hdf5.h>

コンパイル時にはlibraryのリンクを行います

g++ main.cpp -L$(YOUR_HDF_LIBRARY_PATH) -lhdf

中身を確認するときは，

h5dump data.h5

です．HDFViewというGUIも提供されています．

ファイルのオープン，クローズ

H5Fcreateでファイルが作成できます．H5F...とつくものはファイルハンドルに関連するものです．

char* fname;
hid_t fid = H5Fcreate(fname,H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);

ここで，H5F_ACC_TRUNCを与えています．これは，open modeで，

H5F_ACC_TRUNC	ファイルが存在しても強制的に上書き
H5F_ACC_EXCL	ファイルが存在する場合，失敗

となっています．H5P_DEFAULTのところを変えるとPropertyを変更できるようですが今回はSkip．

開くときはH5Fopenを使います．

char* fname;
hid_t fid = H5Fopen(fname,H5F_ACC_RDONLY,H5P_DEFAULT);

こちらのmodeは

H5F_ACC_RDONLY	読み込みonly
H5F_ACC_RDWR	読み書き可能

です．

開いたファイルは，H5Fcloseで閉じるのを忘れないようにしないといけません．

H5Fclose(fid);

グループの作成

H5Gcreateでグループを作成できます．GはGroupのGです．

char* gname;
H5Gcreate(fid,gname,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);

データセットの作成

データセットを作成するためには，まずH5ScreateでdataSpaceを作成する必要があります．

hsize_t rank = N;
hsize_t dims[N] = {...};
hid_t dataspace = H5Screate_simple(rank,dims,NULL)

ここで，rankとdimsはそれぞれ$N$次元配列の$N$と，それぞれの要素数($Nx,Ny,Nz$,...)に対応しています．

その後H5DcreateでDatasetを作成します．

char* dname;
hid_t dataset = H5Dcreate(fid,dname,H5T_STD_I32LE,dataspace,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);

ここでは4byte(32bit)のint Typeを保存するためのdatasetを作成するためにH5T_STD_I32LEを与えています．

他には例えば，以下があります．

Int	H5T_STD_I32LE
char	H5T_STD_I8LE
float	H5T_IEEE_F32LE
double	H5T_IEEE_F64LE

作成したデータセットにデータを書き込みます．

int* data; // = [Nx*Ny*Nz*...]
H5Dwrite(dataset,H5T_NATIVE_INT,H5S_ALL,H5S_ALL,H5P_DEFAULT,data);

H5T_NATIVE_INT, (H5T_NATIVE_CHAR, H5T_NATIVE_FLOAT, H5T_NATIVE_DOUBLE,...)は書き込みデータのtype指定．H5S_ALLはdataspace全てを使うという意味です．部分的な変更については後で書きます．最後の引数にdataの先頭のポインタを渡します．

dataset，dataspaceも使ったあとは忘れずにcloseします．

H5Sclose(dataspace);
H5Dclose(dataset);

データセットのpadding

データセットを何か特定の値で埋めることもできます．これにはH5P(Property)を用います．

int fillval = -1;
hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
H5Pset_fill_value(plist,H5T_NATIVE_INT,&fillval);
hid_t dataset = H5Dcreate(fid,dname,H5T_STD_I32LE,dataspace,H5P_DEFAULT,plist,H5P_DEFAULT);
H5Pclose(plist);

dataset_createに用いるproperty plistを定義し，fill_valueを与え，H5Dcreateの引数に用います．使用したあとはcloseします．

データセットの(部分)更新

あらかじめ確保されたdatasetの(部分)更新を行う場合，流れとしては，datasetのopen，書き込み用の領域をoffsetなど指定して書き込むという感じになります．他にもstrideやblockも指定できるようですが一旦割愛します，

int* data;
char* dname;
hid_t dataset = H5Dopen(fid,dname,H5P_DEFAULT);
hid_t dataspace = H5Dget_space(dataset);
hsize_t rank;
hsize_t dims[rank];
hsize_t offsets[rank];
hid_t dataspace_buf = H5Screate_simple(rank,dims,NULL);
H5Sselect_hyperslab(dataspace,H5S_SELECT_SET,offsets,NULL,dims,NULL);
H5Dwrite(dataset,H5T_NATIVE_INT,dataspace_buf,dataspace,H5P_DEFAULT,data);

データセットの読み込み

作成したdatasetをopenし，dataspaceを取得します．とりあえずこれをreadに渡すとdataset全体を読み出せます．

hid_t dataset = H5Dopen(fid,dname,H5P_DEFAULT);
hid_t dataspace = H5Dget_space(dataset);
H5Dread(dataset,H5T_NATIVE_INT, dataspace, dataspace, H5P_DEFAULT, data);

以上で，

ファイルのopen close
グループの作成
データセットを配列データから作成，または，ある値で初期化する．
データセットを更新する
データセットを読み込む

という操作ができるようになりました．本当はparallel hdfまで勉強したかったのですがまた今度...

機能をまとめてみる

個人的にはh5pyの使用感が好きなので，少し意識しつつ，以上の機能を構造体にまとめてみました．

# include<hdf5.h>
# include<cassert>
# include<iostream>
# include<vector>
# include<string>
# include<type_traits>

template<typename T>
constexpr bool false_v = false;

template <typename T>
hid_t get_predtype()
{
  if constexpr (std::is_same<T,char>{}){
    return H5T_NATIVE_CHAR;
  }else if constexpr (std::is_same<T,int>{}){
    return H5T_NATIVE_INT;  
  }else if constexpr (std::is_same<T,float>{}){
    return H5T_NATIVE_FLOAT;  
  }else if constexpr (std::is_same<T,double>{}){
    return H5T_NATIVE_DOUBLE;  
  }else{
    static_assert(false_v<T>,"cannot find predtype");
  }
  return 0;
}

template <typename T>
hid_t get_dspace_type()
{
  if constexpr (std::is_same<T,char>{}){
    return H5T_STD_I8LE;
  }else if constexpr (std::is_same<T,int>{}){
    return H5T_STD_I32LE;  
  }else if constexpr (std::is_same<T,float>{}){
    return H5T_IEEE_F32LE;  
  }else if constexpr (std::is_same<T,double>{}){
    return H5T_IEEE_F64LE;  
  }else{
    static_assert(false_v<T>,"cannot find dspacetype");
  }
  return 0;
}



struct h5fp
{
private:
  hid_t fid;

public:
  h5fp(std::string fname, char mode)
  { open(fname,mode); };

  ~h5fp()
  { close(); }

  void open(std::string fname, char mode)
  {
    switch (mode)
    {
    case 'w':
      fid = H5Fcreate(fname.c_str(),H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);
      break;
    case 'r':
      fid = H5Fopen(fname.c_str(),H5F_ACC_RDONLY,H5P_DEFAULT);      
      break;    
    default:
      assert(!"[error] invalid mode from h5f::open"); 
      break;
    }
  }

  void close()
  {
    H5Fclose(fid);    
    return;
  }

  void create_group(std::string gname)
  {
    hid_t gid = H5Gcreate(fid,gname.c_str(),H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
    H5Gclose(gid);
    return;
  }

  template<typename T>
  void create_dataset(std::string dname,std::vector<int> shape,T data[])
  {
    size_t sz = shape.size();
    hsize_t rank = (hsize_t)sz;
    hsize_t dims[sz];
    for (size_t i = 0; i < sz; i++) dims[i] = (hsize_t)shape[i];
    hid_t dataspace = H5Screate_simple(rank,dims,NULL);
    hid_t dataset = H5Dcreate(fid,dname.c_str(),get_dspace_type<T>(),dataspace,H5P_DEFAULT,H5P_DEFAULT,H5P_DEFAULT);
    H5Dwrite(dataset,get_predtype<T>(),H5S_ALL,H5S_ALL,H5P_DEFAULT,data);
    H5Dclose(dataset);
    H5Sclose(dataspace);
    return;
  }

  template<typename T>
  void update(std::string dname,std::vector<int> shape,std::vector<int> offset,T data[])
  {
    hid_t dataset = H5Dopen(fid,dname.c_str(),H5P_DEFAULT);
    hid_t dataspace = H5Dget_space(dataset);
    size_t sz = shape.size();
    hsize_t rank = (hsize_t)sz;
    hsize_t dims[sz];
    hsize_t offsets[sz];
    for (size_t i = 0; i < sz; i++){ 
      dims[i] = (hsize_t)shape[i];
      offsets[i] = (hsize_t)offset[i];
    }
    hid_t dataspace_buf = H5Screate_simple(rank,dims,NULL);
    H5Sselect_hyperslab(dataspace,H5S_SELECT_SET,offsets,NULL,dims,NULL);
    H5Dwrite(dataset,get_predtype<T>(),dataspace_buf,dataspace,H5P_DEFAULT,data);
    H5Dclose(dataset);
    H5Sclose(dataspace);
    return;
  }

  template<typename T>
  void create_dataset(std::string dname,std::vector<int> shape,T fillval)
  {
    size_t sz = shape.size();
    hsize_t rank = (hsize_t)sz;
    hsize_t dims[sz];
    for (size_t i = 0; i < sz; i++) dims[i] = (hsize_t)shape[i];
    hid_t dataspace = H5Screate_simple(rank,dims,NULL);
    hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
    H5Pset_fill_value(plist,get_predtype<T>(),&fillval);
    hid_t dataset = H5Dcreate(fid,dname.c_str(),get_dspace_type<T>(),dataspace,H5P_DEFAULT,plist,H5P_DEFAULT);
    H5Dclose(dataset);
    H5Sclose(dataspace);
    H5Pclose(plist);
    return;
  }

  template<typename T>
  void read_dataset(std::string dname, T data[] )
  {
    hid_t dataset = H5Dopen(fid,dname.c_str(),H5P_DEFAULT);
    hid_t dataspace = H5Dget_space(dataset);
    H5Dread(dataset, get_predtype<T>(), dataspace, dataspace, H5P_DEFAULT, data);
    H5Sclose(dataspace);
    H5Dclose(dataset);
    return;
  }

  std::vector<hsize_t> get_shape(std::string dname)
  {
    hid_t dataset = H5Dopen(fid,dname.c_str(),H5P_DEFAULT);
    hid_t dataspace = H5Dget_space(dataset);
    hsize_t rank = H5Sget_simple_extent_ndims(dataspace);
    std::vector<hsize_t> dims(rank);
    H5Sget_simple_extent_dims(dataspace,dims.data(),NULL);
    H5Dclose(dataset);
    H5Sclose(dataspace);    
    return dims;
  }


};

使って遊んでみました．

  float a[5][8];
  for (int i = 0; i < 5; i++)
  {
    for (int j = 0; j < 8; j++)
    {
      a[i][j] = (i+1)*(j+1);
    }    
  }
  double c[6]={0};
  for (int i = 0; i < 6; i++)
  {
    c[i] = 10.0;
  }
    
  {
    h5fp fp("data.h5",'w');
    fp.create_group("/g");
    fp.create_group("/g/gg1");
    fp.create_group("/g/gg2");
    std::cout<<"created groups"<<std::endl;
    fp.create_dataset("/g/gg1/a",{5,8},&a[0][0]);
    fp.create_dataset("/g/gg2/c",{10,6},0.0);
    for (int i = 0; i < 8; i++)
    {
      for (int j = 0; j < 6; j++)
      {
        c[j] = i;
      }      
      fp.update("/g/gg2/c",{1,6},{i,0},c);
    }    
    std::cout<<"wrote a"<<std::endl;
  }


  int b[5][8];
  for (int i = 0; i < 5; i++)
  {
    for (int j = 0; j < 8; j++)
    {
      b[i][j] = -100;
    }    
  }

  {
    h5fp fp("data.h5",'r');
    std::cout<<"reopened"<<std::endl;
    fp.read_dataset("/g/gg1/a",&b[0][0]);
    std::cout<<"read /g/gg1/a"<<std::endl;
    for (size_t i = 0; i < 5; i++)
    {
      for (size_t j = 0; j < 8; j++)
      {
        std::cout<<b[i][j]<<" ";
      }
      std::cout<<std::endl;
    }
  }

created groups
wrote areopened
read /g/gg1/a1 
2 3 4 5 6 7 8 2 4 6 8 10 12 14 16 3 6 9 12 15 18 21 24 4 8 12 16 20 24 28 32 5 10 15 20 25 30 35 40

シンプルですが自分で使う分には悪くないと思います．もっと機能が欲しいですが，その辺は必要に応じて勉強して増やす感じですかね

追記: parallel HDFへの対応

利用の準備

buildする場合，--enable-parallelが必要です．

./configure --prefix=$(YOUR_INSTALL_PATH) --enable-parallel

Homebrewの場合，

brew install hdf5-mpi

です．

ファイルのオープン

ファイルアクセスハンドラを渡す必要があります．

hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
                 H5Pset_fapl_mpio(plist_id,MPI_COMM_WORLD,MPI_INFO_NULL);
fid = H5Fcreate(fname.c_str(),H5F_ACC_TRUNC,H5P_DEFAULT,plist_id);
H5Pclose(plist_id);

後は普通のhdfと同じように扱います．offsetなどをmpi_rankによって変えながら同一のデータセットに書き出せば良いだけです．

上記でまとめた構造体を使うなら，以下のような変更を加えます．

  void open(std::string fname, char mode)
  {
    hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
                     H5Pset_fapl_mpio(plist_id,MPI_COMM_WORLD,MPI_INFO_NULL);

    switch (mode)
    {
    case 'w':
      fid = H5Fcreate(fname.c_str(),H5F_ACC_TRUNC,H5P_DEFAULT,plist_id);
      break;
    case 'r':
      fid = H5Fopen(fname.c_str(),H5F_ACC_RDONLY,plist_id);      
      break;    
    default:
      assert(!"[error] invalid mode from h5f::open"); 
      break;
    }
    H5Pclose(plist_id);
    return;
  }

  void close()
  {
    int already_finalized;
    MPI_Finalized(&already_finalized);
    if( !already_finalized ){
      H5Fclose(fid);    
    }else{
      H5Fclose(fid);    
      assert(!"coding error! you must destroy ph5fp before MPI_FINALIZE");
    }
    return;
  }

以下はお試しで作った拡散方程式を計算するコードの一部です．
呼び出す際，以下のようにoffsetをrankに依存させれば，独立の場所に書き出せます．

fp.create_dataset(dname,{p.Nx,p.Ny},0.0);
for (int i = 1; i < p.Nx_local+1; i++){
  fp.update(dname,{1,p.Ny},{p.Nx_local*p.mpi_rank+(i-1),0},&v[i][1]);
}

後書き

1年くらい前にc++ apiを一瞬使っていたのですが記憶が吹っ飛んでいました．日頃，まとめるのが大事なんだとは思っているのですが実行するのが...