#ifndef PARALLEL_CHECKER_HPP
#define PARALLEL_CHECKER_HPP

#include <utils/pugs_config.hpp>
#ifdef PUGS_HAS_HDF5
#include <highfive/H5File.hpp>
#endif   // PUGS_HAS_HDF5

#include <mesh/Connectivity.hpp>
#include <mesh/ItemArrayVariant.hpp>
#include <mesh/ItemValueVariant.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <utils/Demangle.hpp>
#include <utils/Filesystem.hpp>
#include <utils/Messenger.hpp>
#include <utils/SourceLocation.hpp>

#include <fstream>

template <typename DataType, ItemType item_type, typename ConnectivityPtr>
void parallel_check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
                    const std::string& name,
                    const SourceLocation& source_location = SourceLocation{});

template <typename DataType, ItemType item_type, typename ConnectivityPtr>
void parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value,
                    const std::string& name,
                    const SourceLocation& source_location = SourceLocation{});

class ParallelChecker
{
 public:
  enum class Mode
  {
    automatic,   // write in sequential, read in parallel
    read,
    write
  };

 private:
  static ParallelChecker* m_instance;

  Mode m_mode  = Mode::automatic;
  size_t m_tag = 0;

  std::string m_filename = "parallel_checker.h5";

  ParallelChecker() = default;

 public:
  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
  friend void parallel_check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
                             const std::string& name,
                             const SourceLocation& source_location);

  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
  friend void parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
                             const std::string& name,
                             const SourceLocation& source_location);

#ifdef PUGS_HAS_HDF5
 private:
  template <typename T>
  struct TinyVectorDataType;

  template <size_t Dimension, typename DataT>
  struct TinyVectorDataType<TinyVector<Dimension, DataT>> : public HighFive::DataType
  {
    TinyVectorDataType()
    {
      hsize_t dim[]     = {Dimension};
      auto h5_data_type = HighFive::create_datatype<DataT>();
      _hid              = H5Tarray_create(h5_data_type.getId(), 1, dim);
    }
  };

  template <typename T>
  struct TinyMatrixDataType;

  template <size_t M, size_t N, typename DataT>
  struct TinyMatrixDataType<TinyMatrix<M, N, DataT>> : public HighFive::DataType
  {
    TinyMatrixDataType()
    {
      hsize_t dim[]     = {M, N};
      auto h5_data_type = HighFive::create_datatype<DataT>();
      _hid              = H5Tarray_create(h5_data_type.getId(), 2, dim);
    }
  };

  HighFive::File
  _createOrOpenFileRW() const
  {
    if (m_tag == 0) {
      createDirectoryIfNeeded(m_filename);
      return HighFive::File{m_filename, HighFive::File::Truncate};
    } else {
      return HighFive::File{m_filename, HighFive::File::ReadWrite};
    }
  }

  void
  _printHeader(const std::string& name, const SourceLocation& source_location) const
  {
    std::cout << rang::fg::cyan << " - " << rang::fgB::cyan << "parallel checker" << rang::fg::cyan << " for \""
              << rang::fgB::magenta << name << rang::fg::cyan << "\" tag " << rang::fgB::blue << m_tag
              << rang::fg::reset << '\n';
    std::cout << rang::fg::cyan << " | from " << rang::fgB::blue << source_location.filename() << rang::fg::reset << ':'
              << rang::style::bold << source_location.line() << rang::style::reset << '\n';
  }

  template <typename DataType>
  void
  _writeArray(HighFive::Group& group, const std::string& name, const Array<DataType>& array) const
  {
    using data_type = std::remove_const_t<DataType>;
    if constexpr (is_tiny_vector_v<data_type>) {
      auto dataset = group.createDataSet(name, HighFive::DataSpace{std::vector<size_t>{array.size()}},
                                         TinyVectorDataType<data_type>{});

      dataset.template write_raw<typename data_type::data_type>(&(array[0][0]), TinyVectorDataType<data_type>{});
    } else if constexpr (is_tiny_matrix_v<data_type>) {
      auto dataset = group.createDataSet(name, HighFive::DataSpace{std::vector<size_t>{array.size()}},
                                         TinyMatrixDataType<data_type>{});

      dataset.template write_raw<typename data_type::data_type>(&(array[0](0, 0)), TinyMatrixDataType<data_type>{});
    } else {
      auto dataset = group.createDataSet<data_type>(name, HighFive::DataSpace{std::vector<size_t>{array.size()}});
      dataset.template write_raw<data_type>(&(array[0]));
    }
  }

  template <typename DataType>
  void
  _writeTable(HighFive::Group& group, const std::string& name, const Table<DataType>& table) const
  {
    using data_type = std::remove_const_t<DataType>;
    if constexpr (is_tiny_vector_v<data_type>) {
      auto dataset =
        group.createDataSet(name,
                            HighFive::DataSpace{std::vector<size_t>{table.numberOfRows(), table.numberOfColumns()}},
                            TinyVectorDataType<data_type>{});

      dataset.template write_raw<typename data_type::data_type>(&(table(0, 0)[0]), TinyVectorDataType<data_type>{});
    } else if constexpr (is_tiny_matrix_v<data_type>) {
      auto dataset =
        group.createDataSet(name,
                            HighFive::DataSpace{std::vector<size_t>{table.numberOfRows(), table.numberOfColumns()}},
                            TinyMatrixDataType<data_type>{});

      dataset.template write_raw<typename data_type::data_type>(&(table(0, 0)(0, 0)), TinyMatrixDataType<data_type>{});
    } else {
      auto dataset =
        group.createDataSet<data_type>(name, HighFive::DataSpace{
                                               std::vector<size_t>{table.numberOfRows(), table.numberOfColumns()}});
      dataset.template write_raw<data_type>(&(table(0, 0)));
    }
  }

  template <typename DataType>
  Array<std::remove_const_t<DataType>>
  _readArray(HighFive::Group& group, const std::string& name) const
  {
    using data_type = std::remove_const_t<DataType>;

    auto dataset = group.getDataSet(name);
    Array<data_type> array(dataset.getElementCount());

    if constexpr (is_tiny_vector_v<data_type>) {
      dataset.read<data_type>(&(array[0]), TinyVectorDataType<data_type>{});
    } else if constexpr (is_tiny_matrix_v<data_type>) {
      dataset.read<data_type>(&(array[0]), TinyMatrixDataType<data_type>{});
    } else {
      dataset.read<data_type>(&(array[0]));
    }
    return array;
  }

  template <typename DataType>
  Table<std::remove_const_t<DataType>>
  _readTable(HighFive::Group& group, const std::string& name) const
  {
    using data_type = std::remove_const_t<DataType>;

    auto dataset = group.getDataSet(name);
    Table<data_type> table(dataset.getDimensions()[0], dataset.getDimensions()[1]);

    if constexpr (is_tiny_vector_v<data_type>) {
      dataset.read<data_type>(&(table(0, 0)), TinyVectorDataType<data_type>{});
    } else if constexpr (is_tiny_matrix_v<data_type>) {
      dataset.read<data_type>(&(table(0, 0)), TinyMatrixDataType<data_type>{});
    } else {
      dataset.read<data_type>(&(table(0, 0)));
    }
    return table;
  }

  size_t
  _getConnectivityId(const std::shared_ptr<const IConnectivity>& i_connectivity) const
  {
    switch (i_connectivity->dimension()) {
    case 1: {
      return dynamic_cast<const Connectivity<1>&>(*i_connectivity).id();
    }
    case 2: {
      return dynamic_cast<const Connectivity<2>&>(*i_connectivity).id();
    }
    case 3: {
      return dynamic_cast<const Connectivity<3>&>(*i_connectivity).id();
    }
    default: {
      throw UnexpectedError("unexpected connectivity dimension");
    }
    }
  }

  template <ItemType item_type>
  Array<const int>
  _getItemNumber(const std::shared_ptr<const IConnectivity>& i_connectivity) const
  {
    switch (i_connectivity->dimension()) {
    case 1: {
      const Connectivity<1>& connectivity = dynamic_cast<const Connectivity<1>&>(*i_connectivity);
      return connectivity.number<item_type>().arrayView();
    }
    case 2: {
      const Connectivity<2>& connectivity = dynamic_cast<const Connectivity<2>&>(*i_connectivity);
      return connectivity.number<item_type>().arrayView();
    }
    case 3: {
      const Connectivity<3>& connectivity = dynamic_cast<const Connectivity<3>&>(*i_connectivity);
      return connectivity.number<item_type>().arrayView();
    }
    default: {
      throw UnexpectedError("unexpected connectivity dimension");
    }
    }
  }

  template <ItemType item_type>
  Array<const int>
  _getItemOwner(const std::shared_ptr<const IConnectivity>& i_connectivity) const
  {
    switch (i_connectivity->dimension()) {
    case 1: {
      const Connectivity<1>& connectivity = dynamic_cast<const Connectivity<1>&>(*i_connectivity);
      return connectivity.owner<item_type>().arrayView();
    }
    case 2: {
      const Connectivity<2>& connectivity = dynamic_cast<const Connectivity<2>&>(*i_connectivity);
      return connectivity.owner<item_type>().arrayView();
    }
    case 3: {
      const Connectivity<3>& connectivity = dynamic_cast<const Connectivity<3>&>(*i_connectivity);
      return connectivity.owner<item_type>().arrayView();
    }
    default: {
      throw UnexpectedError("unexpected connectivity dimension");
    }
    }
  }

  template <ItemType item_type>
  void
  _writeItemNumbers(const std::shared_ptr<const IConnectivity> i_connectivity,
                    HighFive::File file,
                    HighFive::Group group) const
  {
    std::string item_number_path = "/connectivities/" + std::to_string(this->_getConnectivityId(i_connectivity)) + '/' +
                                   std::string{itemName(item_type)};

    if (not file.exist(item_number_path)) {
      HighFive::Group item_number_group = file.createGroup(item_number_path);
      this->_writeArray(item_number_group, "numbers", this->_getItemNumber<item_type>(i_connectivity));
    }

    HighFive::DataSet item_numbers = file.getDataSet(item_number_path + "/numbers");
    group.createHardLink("numbers", item_numbers);
  }

  template <typename DataType, ItemType item_type>
  void
  _checkIsComparable(const std::string& name,
                     const SourceLocation& source_location,
                     const std::vector<size_t> data_shape,
                     const std::shared_ptr<const IConnectivity>& i_connectivity,
                     HighFive::Group group) const
  {
    const std::string reference_name          = group.getAttribute("name").read<std::string>();
    const std::string reference_file_name     = group.getAttribute("filename").read<std::string>();
    const std::string reference_function_name = group.getAttribute("function").read<std::string>();
    const size_t reference_line_number        = group.getAttribute("line").read<size_t>();
    const size_t reference_dimension          = group.getAttribute("dimension").read<size_t>();
    const std::string reference_item_type     = group.getAttribute("item_type").read<std::string>();
    const std::string reference_data_type     = group.getAttribute("data_type").read<std::string>();

    bool is_comparable = true;
    if (i_connectivity->dimension() != reference_dimension) {
      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different support dimensions: reference ("
                << rang::fgB::yellow << reference_dimension << rang::fgB::red << ") / target (" << rang::fgB::yellow
                << i_connectivity->dimension() << rang::fg::reset << ")\n";
      is_comparable = false;
    }
    if (itemName(item_type) != reference_item_type) {
      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different item types: reference (" << rang::fgB::yellow
                << reference_item_type << rang::fgB::red << ") / target (" << rang::fgB::yellow << itemName(item_type)
                << rang::fg::reset << ")\n";
      is_comparable = false;
    }
    if (demangle<DataType>() != reference_data_type) {
      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different data types: reference (" << rang::fgB::yellow
                << reference_data_type << rang::fgB::red << ") / target (" << rang::fgB::yellow << demangle<DataType>()
                << rang::fg::reset << ")\n";
      is_comparable = false;
    }
    std::vector reference_data_shape = group.getDataSet(reference_name).getSpace().getDimensions();
    if (reference_data_shape.size() != data_shape.size()) {
      std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different data shape kind: reference ("
                << rang::fgB::yellow << reference_data_shape.size() << "d array" << rang::fgB::red << ") / target ("
                << rang::fgB::yellow << data_shape.size() << "d array" << rang::fg::reset << ")\n";
      is_comparable = false;
    }
    {
      bool same_shapes = true;
      for (size_t i = 1; i < reference_data_shape.size(); ++i) {
        same_shapes &= (reference_data_shape[i] == data_shape[i]);
      }
      if (not same_shapes) {
        std::cout << rang::fg::cyan << " | " << rang::fgB::red << "different data shape: reference ("
                  << rang::fgB::yellow << "*";
        for (size_t i = 1; i < reference_data_shape.size(); ++i) {
          std::cout << ":" << reference_data_shape[i];
        }
        std::cout << rang::fgB::red << ") / target (" << rang::fgB::yellow << "*";
        for (size_t i = 1; i < reference_data_shape.size(); ++i) {
          std::cout << ":" << data_shape[i];
        }
        std::cout << rang::fg::reset << ")\n";
        is_comparable = false;
      }
    }
    if (name != reference_name) {
      // Just warn for different labels (maybe useful for some kind of
      // debugging...)
      std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "different names: reference (" << rang::fgB::yellow
                << reference_name << rang::fgB::magenta << ") / target (" << rang::fgB::yellow << name
                << rang::fg::reset << ")\n";
      std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "reference from " << rang::fgB::blue
                << reference_file_name << rang::fg::reset << ':' << rang::style::bold << reference_line_number
                << rang::style::reset << '\n';
      if ((reference_function_name.size() > 0) or (source_location.function().size() > 0)) {
        std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "reference function " << rang::fgB::blue
                  << reference_function_name << rang::fg::reset << '\n';
        std::cout << rang::fg::cyan << " | " << rang::fgB::magenta << "target function " << rang::fgB::blue
                  << source_location.function() << rang::fg::reset << '\n';
      }
    }

    if (not parallel::allReduceAnd(is_comparable)) {
      throw NormalError("cannot compare data");
    }
  }

 private:
  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
  void
  write(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
        const std::string& name,
        const SourceLocation& source_location)
  {
    HighFive::SilenceHDF5 m_silence_hdf5{true};
    this->_printHeader(name, source_location);

    try {
      HighFive::File file = this->_createOrOpenFileRW();

      auto group = file.createGroup("values/" + std::to_string(m_tag));

      group.createAttribute("filename", std::string{source_location.filename()});
      group.createAttribute("function", std::string{source_location.function()});
      group.createAttribute("line", static_cast<size_t>(source_location.line()));
      group.createAttribute("name", name);

      std::shared_ptr<const IConnectivity> i_connectivity = item_value.connectivity_ptr();
      group.createAttribute("dimension", static_cast<size_t>(i_connectivity->dimension()));
      group.createAttribute("item_type", std::string{itemName(item_type)});
      group.createAttribute("data_type", demangle<DataType>());

      this->_writeArray(group, name, item_value.arrayView());

      this->_writeItemNumbers<item_type>(i_connectivity, file, group);

      ++m_tag;

      std::cout << rang::fg::cyan << " | writing " << rang::fgB::green << "success" << rang::fg::reset << '\n';
    }
    catch (HighFive::Exception& e) {
      throw NormalError(e.what());
    }
  }

  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
  void
  write(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
        const std::string& name,
        const SourceLocation& source_location)
  {
    HighFive::SilenceHDF5 m_silence_hdf5{true};
    this->_printHeader(name, source_location);

    try {
      HighFive::File file = this->_createOrOpenFileRW();

      auto group = file.createGroup("values/" + std::to_string(m_tag));

      group.createAttribute("filename", std::string{source_location.filename()});
      group.createAttribute("function", std::string{source_location.function()});
      group.createAttribute("line", static_cast<size_t>(source_location.line()));
      group.createAttribute("name", name);

      std::shared_ptr<const IConnectivity> i_connectivity = item_array.connectivity_ptr();
      group.createAttribute("dimension", static_cast<size_t>(i_connectivity->dimension()));
      group.createAttribute("item_type", std::string{itemName(item_type)});
      group.createAttribute("data_type", demangle<DataType>());

      this->_writeTable(group, name, item_array.tableView());

      this->_writeItemNumbers<item_type>(i_connectivity, file, group);

      ++m_tag;

      std::cout << rang::fg::cyan << " | writing " << rang::fgB::green << "success" << rang::fg::reset << '\n';
    }
    catch (HighFive::Exception& e) {
      throw NormalError(e.what());
    }
  }

  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
  void
  compare(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
          const std::string& name,
          const SourceLocation& source_location)
  {
    HighFive::SilenceHDF5 m_silence_hdf5{true};
    this->_printHeader(name, source_location);

    try {
      HighFive::File file{m_filename, HighFive::File::ReadOnly};

      auto group = file.getGroup("values/" + std::to_string(m_tag));

      const std::string reference_name = group.getAttribute("name").read<std::string>();

      std::shared_ptr<const IConnectivity> i_connectivity = item_value.connectivity_ptr();

      this->_checkIsComparable<DataType, item_type>(name, source_location,
                                                    std::vector<size_t>{item_value.numberOfItems()}, i_connectivity,
                                                    group);

      Array<const int> reference_item_numbers = this->_readArray<int>(group, "numbers");

      Array<const DataType> reference_item_value = this->_readArray<DataType>(group, reference_name);

      Array<const int> item_numbers = this->_getItemNumber<item_type>(i_connectivity);

      using ItemId = ItemIdT<item_type>;

      std::unordered_map<int, ItemId> item_number_to_item_id_map;

      for (ItemId item_id = 0; item_id < item_numbers.size(); ++item_id) {
        const auto& [iterator, success] =
          item_number_to_item_id_map.insert(std::make_pair(item_numbers[item_id], item_id));

        if (not success) {
          throw UnexpectedError("item numbers have duplicate values");
        }
      }

      Assert(item_number_to_item_id_map.size() == item_numbers.size());

      Array<int> index_in_reference(item_numbers.size());
      index_in_reference.fill(-1);
      for (size_t i = 0; i < reference_item_numbers.size(); ++i) {
        const auto& i_number_to_item_id = item_number_to_item_id_map.find(reference_item_numbers[i]);
        if (i_number_to_item_id != item_number_to_item_id_map.end()) {
          index_in_reference[i_number_to_item_id->second] = i;
        }
      }

      if (parallel::allReduceMin(min(index_in_reference)) < 0) {
        throw NormalError("some item numbers are not defined in reference");
      }

      Array<const int> owner = this->_getItemOwner<item_type>(i_connectivity);

      bool has_own_differences = false;
      bool is_same             = true;

      for (ItemId item_id = 0; item_id < item_value.numberOfItems(); ++item_id) {
        if (reference_item_value[index_in_reference[item_id]] != item_value[item_id]) {
          is_same = false;
          if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
            has_own_differences = true;
          }
        }
      }

      is_same             = parallel::allReduceAnd(is_same);
      has_own_differences = parallel::allReduceOr(has_own_differences);

      if (is_same) {
        std::cout << rang::fg::cyan << " | compare: " << rang::fgB::green << "success" << rang::fg::reset << '\n';
      } else {
        if (has_own_differences) {
          std::cout << rang::fg::cyan << " | compare: " << rang::fgB::red << "failed!" << rang::fg::reset;
        } else {
          std::cout << rang::fg::cyan << " | compare: " << rang::fgB::yellow << "not synchronized" << rang::fg::reset;
        }
        std::cout << rang::fg::cyan << " [see \"" << rang::fgB::blue << "parallel_differences_" << m_tag << "_*"
                  << rang::fg::cyan << "\" files for details]" << rang::fg::reset << '\n';

        {
          std::ofstream fout(std::string{"parallel_differences_"} + stringify(m_tag) + std::string{"_"} +
                             stringify(parallel::rank()));

          fout.precision(15);
          for (ItemId item_id = 0; item_id < item_value.numberOfItems(); ++item_id) {
            if (reference_item_value[index_in_reference[item_id]] != item_value[item_id]) {
              const bool is_own_difference = (parallel::rank() == static_cast<size_t>(owner[item_id]));
              if (is_own_difference) {
                fout << rang::fgB::red << "[ own ]" << rang::fg::reset;
              } else {
                fout << rang::fgB::yellow << "[ghost]" << rang::fg::reset;
              }
              fout << " rank=" << parallel::rank() << " owner=" << owner[item_id] << " item_id=" << item_id
                   << " number=" << item_numbers[item_id]
                   << " reference=" << reference_item_value[index_in_reference[item_id]]
                   << " target=" << item_value[item_id]
                   << " difference=" << reference_item_value[index_in_reference[item_id]] - item_value[item_id] << '\n';
              if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
                has_own_differences = true;
              }
            }
          }
        }

        if (parallel::allReduceAnd(has_own_differences)) {
          throw NormalError("calculations differ!");
        }
      }

      ++m_tag;
    }
    catch (HighFive::Exception& e) {
      throw NormalError(e.what());
    }
  }

  template <typename DataType, ItemType item_type, typename ConnectivityPtr>
  void
  compare(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
          const std::string& name,
          const SourceLocation& source_location)
  {
    HighFive::SilenceHDF5 m_silence_hdf5{true};
    this->_printHeader(name, source_location);

    try {
      HighFive::File file{m_filename, HighFive::File::ReadOnly};

      auto group = file.getGroup("values/" + std::to_string(m_tag));

      const std::string reference_name = group.getAttribute("name").read<std::string>();

      std::shared_ptr<const IConnectivity> i_connectivity = item_array.connectivity_ptr();

      this->_checkIsComparable<DataType, item_type>(name, source_location,
                                                    std::vector<size_t>{item_array.numberOfItems(),
                                                                        item_array.sizeOfArrays()},
                                                    i_connectivity, group);

      Array<const int> reference_item_numbers    = this->_readArray<int>(group, "numbers");
      Table<const DataType> reference_item_array = this->_readTable<DataType>(group, reference_name);

      Array<const int> item_numbers = this->_getItemNumber<item_type>(i_connectivity);

      using ItemId = ItemIdT<item_type>;

      std::unordered_map<int, ItemId> item_number_to_item_id_map;

      for (ItemId item_id = 0; item_id < item_numbers.size(); ++item_id) {
        const auto& [iterator, success] =
          item_number_to_item_id_map.insert(std::make_pair(item_numbers[item_id], item_id));

        if (not success) {
          throw UnexpectedError("item numbers have duplicate values");
        }
      }

      Assert(item_number_to_item_id_map.size() == item_numbers.size());

      Array<int> index_in_reference(item_numbers.size());
      index_in_reference.fill(-1);
      for (size_t i = 0; i < reference_item_numbers.size(); ++i) {
        const auto& i_number_to_item_id = item_number_to_item_id_map.find(reference_item_numbers[i]);
        if (i_number_to_item_id != item_number_to_item_id_map.end()) {
          index_in_reference[i_number_to_item_id->second] = i;
        }
      }

      if (parallel::allReduceMin(min(index_in_reference)) < 0) {
        throw NormalError("some item numbers are not defined in reference");
      }

      Array<const int> owner = this->_getItemOwner<item_type>(i_connectivity);

      bool has_own_differences = false;
      bool is_same             = true;

      for (ItemId item_id = 0; item_id < item_array.numberOfItems(); ++item_id) {
        for (size_t i = 0; i < reference_item_array.numberOfColumns(); ++i) {
          if (reference_item_array[index_in_reference[item_id]][i] != item_array[item_id][i]) {
            is_same = false;
            if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
              has_own_differences = true;
            }
          }
        }
      }

      is_same             = parallel::allReduceAnd(is_same);
      has_own_differences = parallel::allReduceOr(has_own_differences);

      if (is_same) {
        std::cout << rang::fg::cyan << " | compare: " << rang::fgB::green << "success" << rang::fg::reset << '\n';
      } else {
        if (has_own_differences) {
          std::cout << rang::fg::cyan << " | compare: " << rang::fgB::red << "failed!" << rang::fg::reset;
        } else {
          std::cout << rang::fg::cyan << " | compare: " << rang::fgB::yellow << "not synchronized" << rang::fg::reset;
        }
        std::cout << rang::fg::cyan << " [see \"" << rang::fgB::blue << "parallel_differences_" << m_tag << "_*"
                  << rang::fg::cyan << "\" files for details]" << rang::fg::reset << '\n';

        {
          std::ofstream fout(std::string{"parallel_differences_"} + stringify(m_tag) + std::string{"_"} +
                             stringify(parallel::rank()));

          fout.precision(15);
          for (ItemId item_id = 0; item_id < item_array.numberOfItems(); ++item_id) {
            for (size_t i = 0; i < item_array.sizeOfArrays(); ++i) {
              if (reference_item_array[index_in_reference[item_id]][i] != item_array[item_id][i]) {
                const bool is_own_difference = (parallel::rank() == static_cast<size_t>(owner[item_id]));
                if (is_own_difference) {
                  fout << rang::fgB::red << "[ own ]" << rang::fg::reset;
                } else {
                  fout << rang::fgB::yellow << "[ghost]" << rang::fg::reset;
                }
                fout << " rank=" << parallel::rank() << " owner=" << owner[item_id] << " item_id=" << item_id
                     << " column=" << i << " number=" << item_numbers[item_id]
                     << " reference=" << reference_item_array[index_in_reference[item_id]][i]
                     << " target=" << item_array[item_id][i]
                     << " difference=" << reference_item_array[index_in_reference[item_id]][i] - item_array[item_id][i]
                     << '\n';
                if (static_cast<size_t>(owner[item_id]) == parallel::rank()) {
                  has_own_differences = true;
                }
              }
            }
          }
        }

        if (parallel::allReduceAnd(has_own_differences)) {
          throw NormalError("calculations differ!");
        }
      }

      ++m_tag;
    }
    catch (HighFive::Exception& e) {
      throw NormalError(e.what());
    }
  }

#else    // PUGS_HAS_HDF5

  template <typename T>
  void
  write(const T&, const std::string&, const SourceLocation&)
  {
    throw UnexpectedError("parallel checker cannot be used without HDF5 support");
  }

  template <typename T>
  void
  compare(const T&, const std::string&, const SourceLocation&)
  {
    throw UnexpectedError("parallel checker cannot be used without HDF5 support");
  }
#endif   // PUGS_HAS_HDF5

 public:
  static void create();
  static void destroy();

  static ParallelChecker&
  instance()
  {
    return *m_instance;
  }

  Mode
  mode() const
  {
    return m_mode;
  }

  void
  setMode(const Mode& mode)
  {
    if (m_tag != 0) {
      throw UnexpectedError("Cannot modify parallel checker mode if it was already used");
    }
    m_mode = mode;
  }

  const std::string&
  filename() const
  {
    return m_filename;
  }

  void
  setFilename(const std::string& filename)
  {
    if (m_tag != 0) {
      throw UnexpectedError("Cannot modify parallel checker file if it was already used");
    }
    m_filename = filename;
  }

  bool
  isWriting() const
  {
    bool is_writting = false;
    switch (m_mode) {
    case Mode::automatic: {
      is_writting = (parallel::size() == 1);
      break;
    }
    case Mode::write: {
      is_writting = true;
      break;
    }
    case Mode::read: {
      is_writting = false;
      break;
    }
    }

    if ((is_writting) and (parallel::size() > 1)) {
      throw NotImplementedError("parallel check write in parallel");
    }

    return is_writting;
  }
};

template <typename DataType, ItemType item_type, typename ConnectivityPtr>
void
parallel_check(const ItemArray<DataType, item_type, ConnectivityPtr>& item_array,
               const std::string& name,
               const SourceLocation& source_location)
{
  if (ParallelChecker::instance().isWriting()) {
    ParallelChecker::instance().write(item_array, name, source_location);
  } else {
    ParallelChecker::instance().compare(item_array, name, source_location);
  }
}

template <typename DataType, ItemType item_type, typename ConnectivityPtr>
void
parallel_check(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value,
               const std::string& name,
               const SourceLocation& source_location)
{
  if (ParallelChecker::instance().isWriting()) {
    ParallelChecker::instance().write(item_value, name, source_location);
  } else {
    ParallelChecker::instance().compare(item_value, name, source_location);
  }
}

template <size_t Dimension, typename DataType>
void PUGS_INLINE
parallel_check(const DiscreteFunctionP0<Dimension, DataType>& discrete_function,
               const std::string& name,
               const SourceLocation& source_location = SourceLocation{})
{
  parallel_check(discrete_function.cellValues(), name, source_location);
}

template <size_t Dimension, typename DataType>
void PUGS_INLINE
parallel_check(const DiscreteFunctionP0Vector<Dimension, DataType>& discrete_function,
               const std::string& name,
               const SourceLocation& source_location = SourceLocation{})
{
  parallel_check(discrete_function.cellArrays(), name, source_location);
}

void parallel_check(const ItemValueVariant& item_value_variant,
                    const std::string& name,
                    const SourceLocation& source_location = SourceLocation{});

void parallel_check(const DiscreteFunctionVariant& discrete_function_variant,
                    const std::string& name,
                    const SourceLocation& source_location = SourceLocation{});

#endif   // PARALLEL_CHECKER_HPP