LoginSignup
10
10

More than 1 year has passed since last update.

@Miya42

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

Last updated at Posted at 2021-11-20

[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でdata*S*paceを作成する必要があります．

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を一瞬使っていたのですが記憶が吹っ飛んでいました．日頃，まとめるのが大事なんだとは思っているのですが実行するのが...

10
10
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
10
10