Draft: HDF5 interface

CMakeLists.txt

@@ -16,6 +16,7 @@ option(docasa         "Enable CASA"                                     off)
 option(docs           "Build documentation"                             off)
 option(coverage       "Build coverage"                                  off)
 option(onnxrt         "Build with ONNXruntime interface"                off)
+option(hdf5           "Build with HDF5 interface"                       off)
 
 if(NOT CMAKE_BUILD_TYPE)
   message(STATUS "Setting build type to 'Release' as none was specified.")

cmake_files/LookUpHighFive.cmake

@@ -0,0 +1,40 @@
+# Installs BlueBrain/HighFive into build directory
+#
+# - GIT_REPOSITORY: defaults to origin BlueBrain/HighFive repo on github
+# - GIT_TAG: defaults to master
+include(ExternalProject)
+if(HighFive_ARGUMENTS)
+  cmake_parse_arguments(HighFive "" "GIT_REPOSITORY;GIT_TAG" ""
+    ${HighFive_ARGUMENTS})
+endif()
+if(NOT HighFive_GIT_REPOSITORY)
+  set(HighFive_GIT_REPOSITORY https://github.com/BlueBrain/HighFive)
+endif()
+if(NOT HighFive_GIT_TAG)
+  set(HighFive_GIT_TAG master)
+endif()
+
+set(HighFive_DIR "${CMAKE_INSTALL_PREFIX}/external")
+ExternalProject_Add(
+  Lookup-HighFive
+  GIT_REPOSITORY ${HighFive_GIT_REPOSITORY}
+  GIT_TAG ${HighFive_GIT_TAG}
+  PREFIX "${CMAKE_BINARY_DIR}/external"
+  #INSTALL_DIR ${EXTERNAL_ROOT}
+  CMAKE_ARGS
+  -DCMAKE_INSTALL_PREFIX=${HighFive_DIR}
+  -DCMAKE_INSTALL_LIBDIR=${CMAKE_SHARED_LIBRARY_PREFIX}
+  -DHIGHFIVE_EXAMPLES=OFF
+  -DHIGHFIVE_USE_BOOST=OFF
+  -DHIGHFIVE_UNIT_TESTS=OFF
+  # Wrap download, configure and build steps in a script to log output
+  UPDATE_COMMAND ""
+  LOG_DOWNLOAD ON
+  LOG_CONFIGURE ON
+  LOG_BUILD ON
+  LOG_INSTALL ON
+)
+set(HighFive_INCLUDE_DIR "${HighFive_DIR}/include")
+#set(HighFive_LIBRARY_DIR "${HighFive_DIR}/${CMAKE_SHARED_LIBRARY_PREFIX}")
+#set(HighFive_LIBRARIES "HighFive")
+

cmake_files/dependencies.cmake

@@ -129,5 +129,19 @@ if(docasa)
   set(PURIFY_CASACORE TRUE)
 endif()
 
+set(PURIFY_H5 FALSE)
+if(hdf5)
+  find_package(HDF5 COMPONENTS CXX REQUIRED)
+  message(STATUS "Found HDF5 include dir: ${HDF5_INCLUDE_DIR}")
+  message(STATUS "Found HDF5 CXX library: ${HDF5_CXX_LIBRARIES}")
+  find_package(HighFive QUIET)
+  if(NOT HighFive_FOUND)
+    message(STATUS "HighFive not found. Attempt to install...")
+    include(LookUpHighFive)
+    set(PURIFY_HIGHFIVE_LOOKUP TRUE)
+  endif()
+  set(PURIFY_H5 TRUE)
+endif()
+
 # Add script to execute to make sure libraries in the build tree can be found
 add_to_ld_path("${EXTERNAL_ROOT}/lib")

cpp/CMakeLists.txt

@@ -23,6 +23,8 @@ add_include_dir(
   ${yaml-cpp_INCLUDE_DIR}
   ${Cubature_INCLUDE_DIR}
   ${CImg_INCLUDE_DIR}
+  ${HDF5_INCLUDE_DIR}
+  ${HighFive_INCLUDE_DIR}
 )
 
 

cpp/purify/CMakeLists.txt

@@ -16,8 +16,8 @@ endfunction()
 configure_file(config.in.h "${PROJECT_BINARY_DIR}/include/purify/config.h")
 
 set(HEADERS
-  logging.h kernels.h pfitsio.h  types.h
-  IndexMapping.h utilities.h operators.h  wproj_utilities.h
+  logging.h kernels.h pfitsio.h types.h
+  IndexMapping.h utilities.h operators.h wproj_utilities.h
   cimg.h uvfits.h convolution.h measurement_operator_factory.h wavelet_operator_factory.h distribute.h
   update_factory.h
   convergence_factory.h
@@ -42,6 +42,10 @@ if(PURIFY_CASACORE)
   list(APPEND HEADERS casacore.h)
 endif()
 
+if(PURIFY_H5)
+  list(APPEND HEADERS h5reader.h)
+endif()
+
 if(PURIFY_MPI)
   list(APPEND HEADERS mpi_utilities.h distribute.h DistributeSparseVector.h
   random_update_factory.h
@@ -70,6 +74,11 @@ if(PURIFY_ARRAYFIRE)
   target_link_libraries(libpurify ${ArrayFire_LIBRARIES})
 endif()
 
+if(PURIFY_H5)
+  target_link_libraries(libpurify ${HDF5_CXX_LIBRARIES})
+  target_include_directories(libpurify SYSTEM PUBLIC ${HDF5_INCLUDE_DIR} ${HighFive_INCLUDE_DIR})
+endif()
+
 target_include_directories(libpurify PUBLIC
   $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/..>
   $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
@@ -111,6 +120,10 @@ if(PURIFY_CASACORE_LOOKUP)
   add_dependencies(libpurify Lookup-CasaCore)
 endif()
 
+if(PURIFY_HIGHFIVE_LOOKUP)
+  add_dependencies(libpurify Lookup-HighFive)
+endif()
+
 install(FILES ${HEADERS} DESTINATION include/purify)
 install(TARGETS libpurify
   EXPORT PurifyTargets

cpp/purify/config.in.h

@@ -28,6 +28,9 @@
 //! Whether PURIFY is using (and SOPT was built with) onnxrt support
 #cmakedefine PURIFY_ONNXRT
 
+//! Whether PURIFY is using HDF5 support
+#cmakedefine PURIFY_H5
+
 #include <string>
 #include <tuple>
 #include <cstdint>

cpp/purify/h5reader.h

@@ -0,0 +1,34 @@
+#ifndef PURIFY_H5READER_H
+#define PURIFY_H5READER_H
+#include "purify/config.h"
+#include "purify/types.h"
+#include "purify/logging.h"
+
+#include "highfive/H5File.hpp"
+
+#include <string>
+#include <vector>
+
+namespace purify {
+
+/// @brief Purify interface class to handle HDF5 input files
+class H5Handler {
+ public:
+  H5Handler() = delete;
+
+  H5Handler(const std::string& filename) : _file(filename) {}
+
+  template <typename T = double>
+  std::vector<T> read(const std::string& label) {
+    auto dataset = _file.getDataSet(label);
+    return dataset.read<std::vector<T>>();
+  }
+
+ private:
+  /// HDF5 file
+  const HighFive::File _file;
+};
+
+}  // namespace purify
+
+#endif

cpp/purify/uvw_utilities.cc

@@ -3,11 +3,19 @@
 #include <fstream>
 #include <random>
 #include <sys/stat.h>
+#ifdef PURIFY_H5
+#include "purify/h5reader.h"
+#endif
 #include "purify/logging.h"
 #include "purify/operators.h"
 
 namespace purify {
 namespace utilities {
+
+bool has_suffix(const std::string &str, const std::string &suff) {
+  return str.size() >= suff.size() && str.compare(str.size() - suff.size(), suff.size(), suff) == 0;
+}
+
 Matrix<t_real> generate_antennas(const t_uint N, const t_real scale) {
   Matrix<t_real> B = Matrix<t_real>::Zero(N, 3);
   const t_real mean = 0;
@@ -170,7 +178,7 @@ t_real streamtoreal(std::ifstream &stream) {
   return std::stod(input);
 }
 
-utilities::vis_params read_visibility(const std::string &vis_name, const bool w_term) {
+utilities::vis_params read_visibility_csv(const std::string &vis_name, const bool w_term) {
   /*
     Reads an csv file with u, v, visibilities and returns the vectors.
 
@@ -220,6 +228,58 @@ utilities::vis_params read_visibility(const std::string &vis_name, const bool w_
   return uv_vis;
 }
 
+#ifdef PURIFY_H5
+utilities::vis_params read_visibility_h5(const std::string &vis_name, const bool w_term) {
+  /*
+    Reads an HDF5 file with u, v, visibilities and returns the vectors.
+
+    vis_name:: name of input HDF5 file containing [u, v, real(V), imag(V)].
+  */
+  H5Handler vis_file(vis_name);
+  utilities::vis_params uv_vis;
+
+  std::vector<t_real> utemp = vis_file.read<t_real>("u");
+  uv_vis.u = Eigen::Map<Vector<t_real>>(utemp.data(), utemp.size(), 1);
+
+  // found that a reflection is needed for the orientation
+  // of the gridded image to be correct
+  std::vector<t_real> vtemp = vis_file.read<t_real>("v");
+  uv_vis.v = -Eigen::Map<Vector<t_real>>(vtemp.data(), vtemp.size(), 1);
+
+  if (w_term) {
+    std::vector<t_real> wtemp = vis_file.read<t_real>("w");
+    uv_vis.w = Eigen::Map<Vector<t_real>>(wtemp.data(), wtemp.size(), 1);
+  } else {
+    uv_vis.w = Vector<t_real>::Zero(utemp.size());
+  }
+
+  std::vector<t_real> retemp = vis_file.read<t_real>("re");
+  std::vector<t_real> imtemp = vis_file.read<t_real>("im");
+  std::vector<t_real> sigma = vis_file.read<t_real>("sigma");
+  assert(retemp.size() == imtemp.size());
+
+  uv_vis.vis = Vector<t_complex>::Zero(retemp.size());
+  uv_vis.weights = Vector<t_complex>::Zero(retemp.size());
+  for (size_t i = 0; i < retemp.size(); ++i) {
+    uv_vis.vis(i) = t_complex(retemp[i], imtemp[i]);
+    uv_vis.weights(i) = 1 / sigma[i];
+  }
+
+  uv_vis.ra = 0;
+  uv_vis.dec = 0;
+  uv_vis.average_frequency = 0;
+
+  return uv_vis;
+}
+#endif
+
+utilities::vis_params read_visibility(const std::string &vis_name, const bool w_term) {
+#ifdef PURIFY_H5
+  if (has_suffix(vis_name, ".h5")) return read_visibility_h5(vis_name, w_term);
+#endif
+  return read_visibility_csv(vis_name, w_term);
+}
+
 void write_visibility(const utilities::vis_params &uv_vis, const std::string &file_name,
                       const bool w_term) {
   /*

cpp/purify/uvw_utilities.h

@@ -189,6 +189,10 @@ utilities::vis_params read_ant_positions_to_coverage(const std::string &pos_name
 }
 //! Reading reals from visibility file (including nan's and inf's)
 t_real streamtoreal(std::ifstream &stream);
+//! Reads in visibility csv file
+utilities::vis_params read_visibility_csv(const std::string &vis_name, const bool w_term = false);
+//! Reads in visibility HDF5 file
+utilities::vis_params read_visibility_h5(const std::string &vis_name, const bool w_term = false);
 //! Reads in visibility file
 utilities::vis_params read_visibility(const std::string &vis_name, const bool w_term = false);
 //! Read visibility files from name of vector

cpp/tests/CMakeLists.txt

@@ -5,7 +5,7 @@ if(Catch2_FOUND)
   target_link_libraries(common_catch_main_object Catch2::Catch2)
 endif()
 if(EIGEN3_INCLUDE_DIR)
-target_include_directories(common_catch_main_object SYSTEM PUBLIC ${EIGEN3_INCLUDE_DIR})
+  target_include_directories(common_catch_main_object SYSTEM PUBLIC ${EIGEN3_INCLUDE_DIR})
 endif()
 if(sopt_FOUND)
   target_include_directories(common_catch_main_object SYSTEM PUBLIC ${sopt_INCLUDE_DIR})
@@ -44,19 +44,23 @@ if(PURIFY_ARRAYFIRE)
   add_catch_test(operators_gpu LIBRARIES libpurify ${ArrayFire_LIBRARIES})
 endif()
 
+if(PURIFY_H5)
+  add_catch_test(purify_h5 LIBRARIES libpurify DEPENDS lookup_dependencies)
+endif()
+
 if(PURIFY_MPI)
   add_library(common_mpi_catch_main_object OBJECT "common_mpi_catch_main.cc")
-  target_include_directories(common_mpi_catch_main_object
-    PUBLIC ${PROJECT_BINARY_DIR}/include ${MPI_CXX_INCLUDE_PATH})
   if(Catch2_FOUND)
     target_link_libraries(common_mpi_catch_main_object Catch2::Catch2)
   endif()
-  if(sopt_FOUND)
-    target_include_directories(common_mpi_catch_main_object SYSTEM PUBLIC ${sopt_INCLUDE_DIR})
-  endif()
   if(EIGEN3_INCLUDE_DIR)
     target_include_directories(common_mpi_catch_main_object SYSTEM PUBLIC ${EIGEN3_INCLUDE_DIR})
   endif()
+  if(sopt_FOUND)
+    target_include_directories(common_mpi_catch_main_object SYSTEM PUBLIC ${sopt_INCLUDE_DIR})
+  endif()
+  target_include_directories(common_mpi_catch_main_object
+    PUBLIC ${PROJECT_BINARY_DIR}/include ${MPI_CXX_INCLUDE_PATH})
 
   function(add_mpi_test testname)
     add_catch_test(${testname} COMMON_MAIN common_mpi_catch_main_object NOTEST ${ARGN})
@@ -80,7 +84,7 @@ if(PURIFY_MPI)
   add_mpi_test(mpi_measurement_factory LIBRARIES libpurify)
   add_mpi_test(mpi_measurement_operator LIBRARIES libpurify)
   add_mpi_test(mpi_read_measurements LIBRARIES libpurify)
-  add_catch_test(mpi_utilities LIBRARIES libpurify)
+  add_mpi_test(mpi_utilities LIBRARIES libpurify)
   add_mpi_test(mpi_wavelet_factory LIBRARIES libpurify)
   add_mpi_test(mpi_wide_field_utilities LIBRARIES libpurify)
   add_mpi_test(parallel_mpi_utilities LIBRARIES libpurify)

cpp/tests/purify_h5.cc

@@ -0,0 +1,36 @@
+#include "catch2/catch_all.hpp"
+
+#include "purify/config.h"
+#include "purify/types.h"
+#include "purify/logging.h"
+
+#include "purify/directories.h"
+#include "purify/h5reader.h"
+
+#include <iostream>
+#include <vector>
+
+using namespace purify;
+using namespace purify::notinstalled;
+
+TEST_CASE("Purify H5", "[HDF5]") {
+  H5Handler f(atca_filename("0332-391.h5"));
+
+  const std::vector<double> u = f.read("u");
+  const std::vector<double> v = f.read("v");
+  const std::vector<double> w = f.read("w");
+  const std::vector<double> re = f.read("re");
+  const std::vector<double> im = f.read("im");
+  const std::vector<double> sigma = f.read("sigma");
+
+  std::cout << "u size = " << u.size() << std::endl;
+  std::cout << "v size = " << v.size() << std::endl;
+  std::cout << "w size = " << w.size() << std::endl;
+  std::cout << "re size = " << re.size() << std::endl;
+  std::cout << "im size = " << im.size() << std::endl;
+  std::cout << "sigma size = " << sigma.size() << std::endl;
+
+  const bool pass = u.size() > 0 && u.size() == v.size() && u.size() == w.size() &&
+                    u.size() == re.size() && u.size() == im.size() && u.size() == sigma.size();
+  CHECK(pass);
+}