serialize.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <iganet/core/core.hpp>
#include <iganet/utils/tensorarray.hpp>
 
#include <nlohmann/json.hpp>
#include <pugixml.hpp>
 
namespace iganet::utils {
 
struct Serializable {
  virtual ~Serializable() = default;
 
  virtual nlohmann::json to_json() const = 0;
 
  virtual void pretty_print(std::ostream &os = Log(log::info)) const = 0;
};
 
template <typename T, std::size_t N>
inline auto to_json(const torch::TensorAccessor<T, N> &accessor) {
  auto json = nlohmann::json::array();
 
  if constexpr (N == 1) {
    for (int64_t i = 0; i < accessor.size(0); ++i)
      json.push_back(accessor[i]);
  } else if constexpr (N == 2) {
    for (int64_t i = 0; i < accessor.size(0); ++i)
      for (int64_t j = 0; j < accessor.size(1); ++j)
        json.push_back(accessor[i][j]);
  } else if constexpr (N == 3) {
    for (int64_t i = 0; i < accessor.size(0); ++i)
      for (int64_t j = 0; j < accessor.size(1); ++j)
        for (int64_t k = 0; k < accessor.size(2); ++k)
          json.push_back(accessor[i][j][k]);
  } else if constexpr (N == 4) {
    for (int64_t i = 0; i < accessor.size(0); ++i)
      for (int64_t j = 0; j < accessor.size(1); ++j)
        for (int64_t k = 0; k < accessor.size(2); ++k)
          for (int64_t l = 0; l < accessor.size(3); ++l)
            json.push_back(accessor[i][j][k][l]);
  }
 
  return json;
}
 
template <typename T, std::size_t N>
inline auto to_json(const torch::Tensor &tensor) {
  if (tensor.is_cuda()) {
    auto [tensor_cpu, accessor] = to_tensorAccessor<T, N>(tensor, torch::kCPU);
    return to_json(accessor);
  } else {
    auto accessor = to_tensorAccessor<T, N>(tensor);
    return to_json(accessor);
  }
}
 
template <typename T, std::size_t N, std::size_t M>
inline auto to_json(const utils::TensorArray<M> &tensors) {
  auto json = nlohmann::json::array();
 
  for (std::size_t i = 0; i < M; ++i) {
    if (tensors[i].is_cuda()) {
      auto [tensor_cpu, accessor] =
          to_tensorAccessor<T, N>(tensors[i], torch::kCPU);
      json.push_back(to_json<T, N>(accessor));
    } else {
      auto accessor = to_tensorAccessor<T, N>(tensors[i]);
      json.push_back(to_json<T, N>(accessor));
    }
  }
 
  return json;
}
 
#ifdef IGANET_WITH_GISMO
template <typename T, int Rows, int Cols, int Options>
inline auto to_json(const gismo::gsMatrix<T, Rows, Cols, Options> &matrix,
                    bool flatten = false, bool transpose = false) {
  auto json = nlohmann::json::array();
 
  if constexpr (Options == gismo::RowMajor) {
    if (flatten) {
      if (transpose) {
        for (std::size_t j = 0; j < matrix.cols(); ++j)
          for (std::size_t i = 0; i < matrix.rows(); ++i)
            json.push_back(matrix(i, j));
      } else {
        for (std::size_t i = 0; i < matrix.rows(); ++i)
          for (std::size_t j = 0; j < matrix.cols(); ++j)
            json.push_back(matrix(i, j));
      }
    } else {
      if (transpose) {
        for (std::size_t j = 0; j < matrix.cols(); ++j) {
          auto data = nlohmann::json::array();
          for (std::size_t i = 0; i < matrix.rows(); ++i) {
            data.push_back(matrix(i, j));
          }
          json.emplace_back(data);
        }
      } else {
        for (std::size_t i = 0; i < matrix.rows(); ++i) {
          auto data = nlohmann::json::array();
          for (std::size_t j = 0; j < matrix.cols(); ++j) {
            data.push_back(matrix(i, j));
          }
          json.emplace_back(data);
        }
      }
    }
 
  } else if constexpr (Options == gismo::ColMajor) {
    if (flatten) {
      if (transpose) {
        for (std::size_t i = 0; i < matrix.rows(); ++i)
          for (std::size_t j = 0; j < matrix.cols(); ++j)
            json.push_back(matrix(i, j));
      } else {
        for (std::size_t j = 0; j < matrix.cols(); ++j)
          for (std::size_t i = 0; i < matrix.rows(); ++i)
            json.push_back(matrix(i, j));
      }
    } else {
      if (transpose) {
        for (std::size_t i = 0; i < matrix.rows(); ++i) {
          auto data = nlohmann::json::array();
          for (std::size_t j = 0; j < matrix.cols(); ++j) {
            data.push_back(matrix(i, j));
          }
          json.emplace_back(data);
        }
      } else {
        for (std::size_t j = 0; j < matrix.cols(); ++j) {
          auto data = nlohmann::json::array();
          for (std::size_t i = 0; i < matrix.rows(); ++i) {
            data.push_back(matrix(i, j));
          }
          json.emplace_back(data);
        }
      }
    }
 
  } else
    throw std::runtime_error("Invalid matrix options");
 
  return json;
}
 
template <typename T> inline auto to_json(const gismo::gsBSpline<T> &bspline) {
  auto json = nlohmann::json();
 
  json["degrees"] = nlohmann::json::array();
 
  for (std::size_t i = 0; i < bspline.parDim(); ++i)
    json["degrees"].push_back(bspline.degree(i));
 
  json["geoDim"] = bspline.geoDim();
  json["parDim"] = bspline.parDim();
 
  json["ncoeffs"] = nlohmann::json::array();
  for (std::size_t i = 0; i < bspline.parDim(); ++i)
    json["ncoeffs"].push_back(bspline.basis().size(i));
 
  json["coeffs"] = to_json(bspline.coefs());
 
  json["nknots"] = nlohmann::json::array();
  for (std::size_t i = 0; i < bspline.parDim(); ++i)
    json["nknots"].push_back(bspline.knots(i).size());
 
  json["knots"] = nlohmann::json::array();
  for (std::size_t i = 0; i < bspline.parDim(); ++i)
    json["knots"].push_back(bspline.knots(i));
 
  return json;
}
 
template <int d, typename T>
inline auto to_json(const gismo::gsTensorBSpline<d, T> &bspline) {
  auto json = nlohmann::json();
 
  json["degrees"] = nlohmann::json::array();
 
  for (std::size_t i = 0; i < bspline.parDim(); ++i)
    json["degrees"].push_back(bspline.degree(i));
 
  json["geoDim"] = bspline.geoDim();
  json["parDim"] = bspline.parDim();
 
  json["ncoeffs"] = nlohmann::json::array();
  for (std::size_t i = 0; i < bspline.parDim(); ++i)
    json["ncoeffs"].push_back(bspline.basis().size(i));
 
  json["coeffs"] = to_json(bspline.coefs());
 
  json["nknots"] = nlohmann::json::array();
  for (std::size_t i = 0; i < bspline.parDim(); ++i)
    json["nknots"].push_back(bspline.knots(i).size());
 
  json["knots"] = nlohmann::json::array();
  for (std::size_t i = 0; i < bspline.parDim(); ++i)
    json["knots"].push_back(bspline.knots(i));
 
  return json;
}
 
template <typename T>
inline auto to_json(const gismo::gsGeometry<T> &geometry) {
 
  if (auto patch = dynamic_cast<const gismo::gsBSpline<T> *>(&geometry))
    return to_json(*patch);
  else if (auto patch =
               dynamic_cast<const gismo::gsTensorBSpline<2, T> *>(&geometry))
    return to_json(*patch);
  else if (auto patch =
               dynamic_cast<const gismo::gsTensorBSpline<3, T> *>(&geometry))
    return to_json(*patch);
  else if (auto patch =
               dynamic_cast<const gismo::gsTensorBSpline<4, T> *>(&geometry))
    return to_json(*patch);
  else
    return nlohmann::json("{ Invalid patch type }");
}
 
template <typename T>
inline auto
to_json(const typename gismo::gsMultiPatch<T>::ifContainer &interfaces) {
 
  auto json = nlohmann::json::array();
 
  for (auto const &interface : interfaces) {
    auto interface_json = nlohmann::json();
 
    interface_json["patches"] = {interface.first().patchIndex(),
                                 interface.second().patchIndex()};
    interface_json["sides"] = {interface.first().side().index(),
                               interface.second().side().index()};
    interface_json["direction"] = "NOT IMPLEMENTED YET";
    interface_json["orientation"] = "NOT IMPLEMENTED YET";
 
    json.push_back(interface_json);
  }
 
  return json;
}
 
template <typename T>
inline auto
to_json(const typename gismo::gsMultiPatch<T>::bContainer &boundaries) {
 
  auto json = nlohmann::json::array();
 
  for (auto const &boundary : boundaries) {
    auto boundary_json = nlohmann::json();
 
    boundary_json["patch"] = boundary.patchIndex();
    boundary_json["side"] = boundary.side().index();
 
    json.push_back(boundary_json);
  }
 
  return json;
}
 
template <typename T>
inline auto to_json(const gismo::gsMultiPatch<T> &mp, bool verbose = false) {
 
  auto json = nlohmann::json();
 
  // Create list of patch indices
  auto patches_json = nlohmann::json::array();
  for (std::size_t i = 0; i < mp.nPatches(); ++i)
    patches_json.push_back(i);
 
  json["patches"] = patches_json;
  json["interfaces"] = to_json<T>(mp.interfaces());
  json["boundaries"] = to_json<T>(mp.boundaries());
 
  if (verbose) {
    auto patches_json = nlohmann::json::array();
 
    for (std::size_t i = 0; i < mp.nPatches(); ++i)
      patches_json.push_back(to_json(mp.patch(i)));
 
    json["patches"] = patches_json;
  }
 
  return json;
}
#endif
 
template <typename T, std::size_t N>
inline pugi::xml_node &to_xml(const torch::TensorAccessor<T, N> &accessor,
                              torch::IntArrayRef sizes, pugi::xml_node &root,
                              std::string tag = "Matrix", int id = 0,
                              const std::string &label = "", int index = -1);
 
template <typename T, std::size_t N>
inline pugi::xml_document to_xml(const torch::TensorAccessor<T, N> &accessor,
                                 torch::IntArrayRef sizes,
                                 std::string tag = "Matrix", int id = 0,
                                 const std::string &label = "", int index = -1) {
  pugi::xml_document doc;
  pugi::xml_node root = doc.append_child("xml");
  to_xml(accessor, sizes, root, tag, id, label, index);
 
  return doc;
}
 
template <typename T, std::size_t N>
inline pugi::xml_node &to_xml(const torch::TensorAccessor<T, N> &accessor,
                              torch::IntArrayRef sizes, pugi::xml_node &root,
                              std::string tag, int id, const std::string &label,
                              int index) {
 
  // add node
  pugi::xml_node node = root.append_child(tag.c_str());
 
  if (id >= 0)
    node.append_attribute("id") = id;
 
  if (index >= 0)
    node.append_attribute("index") = index;
 
  if (!label.empty())
    node.append_attribute("label") = label.c_str();
 
  // add rows/cols or dimensions
  if (tag == "Matrix") {
    if constexpr (N == 1) {
      node.append_attribute("rows") = sizes[0];
      node.append_attribute("cols") = 1;
 
      std::stringstream ss;
      for (std::size_t i = 0; i < sizes[0]; ++i)
        ss << std::to_string(accessor[i]) << (i < sizes[0] - 1 ? " " : "");
      node.append_child(pugi::node_pcdata).set_value(ss.str().c_str());
    } else if constexpr (N == 2) {
      node.append_attribute("rows") = sizes[0];
      node.append_attribute("cols") = sizes[1];
 
      std::stringstream ss;
      for (std::size_t i = 0; i < sizes[0]; ++i)
        for (std::size_t j = 0; j < sizes[1]; ++j)
          ss << std::to_string(accessor[i][j])
             << (j < sizes[1] - 1 ? " " : (i < sizes[0] - 1 ? " " : ""));
      node.append_child(pugi::node_pcdata).set_value(ss.str().c_str());
    } else
      throw std::runtime_error(
          "Tag \"Matrix\" only supports 1- and 2-dimensional tensors");
  } else {
    std::stringstream ss;
    for (const auto &size : sizes)
      ss << std::to_string(size) << " ";
 
    pugi::xml_node dims = node.append_child("Dimensions");
    dims.append_child(pugi::node_pcdata).set_value(ss.str().c_str());
 
    ss.str("");
    if constexpr (N == 1) {
      for (std::size_t i = 0; i < sizes[0]; ++i)
        ss << std::to_string(accessor[i]) << " ";
    } else if constexpr (N == 2) {
      for (std::size_t i = 0; i < sizes[0]; ++i)
        for (std::size_t j = 0; j < sizes[1]; ++j)
          ss << std::to_string(accessor[i][j]) << " ";
    } else if constexpr (N == 3) {
      for (std::size_t i = 0; i < sizes[0]; ++i)
        for (std::size_t j = 0; j < sizes[1]; ++j)
          for (std::size_t k = 0; k < sizes[2]; ++k)
            ss << std::to_string(accessor[i][j][k]) << " ";
    } else if constexpr (N == 4) {
      for (std::size_t i = 0; i < sizes[0]; ++i)
        for (std::size_t j = 0; j < sizes[1]; ++j)
          for (std::size_t k = 0; k < sizes[2]; ++k)
            for (std::size_t l = 0; l < sizes[3]; ++l)
              ss << std::to_string(accessor[i][j][k][l]) << " ";
 
    } else if constexpr (N == 5) {
      for (std::size_t i = 0; i < sizes[0]; ++i)
        for (std::size_t j = 0; j < sizes[1]; ++j)
          for (std::size_t k = 0; k < sizes[2]; ++k)
            for (std::size_t l = 0; l < sizes[3]; ++l)
              for (std::size_t m = 0; m < sizes[4]; ++m)
                ss << std::to_string(accessor[i][j][k][l][m]) << " ";
    } else if constexpr (N == 6) {
      for (std::size_t i = 0; i < sizes[0]; ++i)
        for (std::size_t j = 0; j < sizes[1]; ++j)
          for (std::size_t k = 0; k < sizes[2]; ++k)
            for (std::size_t l = 0; l < sizes[3]; ++l)
              for (std::size_t m = 0; m < sizes[4]; ++m)
                for (std::size_t n = 0; n < sizes[5]; ++n)
                  ss << std::to_string(accessor[i][j][k][l][m][n]) << " ";
    } else
      throw std::runtime_error(
          "Dimensions higher than 6 are not implemented yet");
 
    pugi::xml_node data = node.append_child("Data");
    data.append_child(pugi::node_pcdata).set_value(ss.str().c_str());
  }
 
  return root;
}
 
template <typename T, std::size_t N>
inline pugi::xml_node &to_xml(const torch::Tensor &tensor,
                              pugi::xml_node &root,
                              std::string tag = "Matrix", int id = 0,
                              const std::string &label = "", int index = -1);
 
template <typename T, std::size_t N>
inline pugi::xml_document to_xml(const torch::Tensor &tensor,
                                 std::string tag = "Matrix", int id = 0,
                                 const std::string &label = "", int index = -1) {
  pugi::xml_document doc;
  pugi::xml_node root = doc.append_child("xml");
  to_xml<T, N>(tensor, root, tag, id, label, index);
 
  return doc;
}
 
template <typename T, std::size_t N>
inline pugi::xml_node &to_xml(const torch::Tensor &tensor, pugi::xml_node &root,
                              std::string tag, int id, const std::string &label,
                              int index) {
 
  if (tensor.is_cuda()) {
    auto [tensor_cpu, accessor] = to_tensorAccessor<T, N>(tensor, torch::kCPU);
    return to_xml(accessor, tensor.sizes(), root, tag, id, label, index);
  } else {
    auto accessor = to_tensorAccessor<T, N>(tensor);
    return to_xml(accessor, tensor.sizes(), root, tag, id, label, index);
  }
}
 
template <typename T, std::size_t N, std::size_t M>
inline pugi::xml_node &to_xml(const utils::TensorArray<M> &tensors,
                              pugi::xml_node &root,
                              std::string tag = "Matrix", int id = 0,
                              const std::string &label = "");
 
template <typename T, std::size_t N, std::size_t M>
inline pugi::xml_document to_xml(const utils::TensorArray<M> &tensors,
                                 std::string tag = "Matrix", int id = 0,
                                 const std::string &label = "", int index = -1) {
  pugi::xml_document doc;
  pugi::xml_node root = doc.append_child("xml");
  to_xml<T, N>(tensors, root, tag, id, label);
 
  return doc;
}
 
template <typename T, std::size_t N, std::size_t M>
inline pugi::xml_node &to_xml(const utils::TensorArray<M> &tensors,
                              pugi::xml_node &root, std::string tag, int id,
                              const std::string &label) {
 
  for (std::size_t i = 0; i < M; ++i) {
    if (tensors[i].is_cuda()) {
      auto [tensor_cpu, accessor] =
          to_tensorAccessor<T, N>(tensors[i], torch::kCPU);
      to_xml(accessor, tensors[i].sizes(), root, tag, id, label, i);
    } else {
      auto accessor = to_tensorAccessor<T, N>(tensors[i]);
      to_xml(accessor, tensors[i].sizes(), root, tag, id, label, i);
    }
  }
 
  return root;
}
 
template <typename T, std::size_t N>
inline torch::TensorAccessor<T, N> &
from_xml(const pugi::xml_document &doc, torch::TensorAccessor<T, N> &accessor,
         torch::IntArrayRef sizes, std::string tag = "Matrix", int id = 0,
         const std::string &label = "", int index = -1) {
  return from_xml(doc.child("xml"), accessor, sizes, tag, id, label, index);
}
 
template <typename T, std::size_t N>
inline torch::TensorAccessor<T, N> &
from_xml(const pugi::xml_node &root, torch::TensorAccessor<T, N> &accessor,
         torch::IntArrayRef sizes, std::string tag = "Matrix", int id = 0,
         const std::string &label = "", int index = -1) {
 
  return accessor;
}
 
template <typename T, std::size_t N>
inline torch::Tensor &from_xml(const pugi::xml_node &root,
                               torch::Tensor &tensor,
                               std::string tag = "Matrix", int id = 0,
                               const std::string &label = "", bool alloc = true,
                               int index = -1);
 
template <typename T, std::size_t N>
inline torch::Tensor &
from_xml(const pugi::xml_document &doc, torch::Tensor &tensor,
         std::string tag = "Matrix", int id = 0,
         const std::string &label = "",
         bool alloc = true, int index = -1) {
  return from_xml<T, N>(doc.child("xml"), tensor, tag, id, label, alloc,
                        index);
}
 
template <typename T, std::size_t N>
inline torch::Tensor &
from_xml(const pugi::xml_node &root, torch::Tensor &tensor,
         std::string tag, int id, const std::string &label, bool alloc,
         int index) {
 
  // Loop through all nodes
  for (pugi::xml_node node : root.children(tag.c_str())) {
 
    if ((id >= 0 ? node.attribute("id").as_int() == id : true) &&
        (index >= 0 ? node.attribute("index").as_int() == index : true) &&
        (!label.empty() ? node.attribute("label").value() == label : true)) {
 
      if (tag == "Matrix") {
 
        int64_t rows = node.attribute("rows").as_int();
        int64_t cols = node.attribute("cols").as_int();
 
        if constexpr (N == 1) {
          if (cols != 1)
            throw std::runtime_error("Invalid matrix dimensions");
          if (!alloc && (tensor.dim() != 1 || tensor.size(0) != rows))
            throw std::runtime_error("Invalid matrix dimensions");
          if (alloc && (tensor.dim() != 1 || tensor.size(0) != rows))
            tensor = torch::zeros({rows}, tensor.options());
        } else if constexpr (N == 2) {
          if (!alloc &&
              (tensor.dim() != 2 || tensor.size(0) != rows ||
               tensor.size(1) != cols))
            throw std::runtime_error("Invalid matrix dimensions");
          if (alloc &&
              (tensor.dim() != 2 || tensor.size(0) != rows ||
               tensor.size(1) != cols))
            tensor = torch::zeros({rows, cols}, tensor.options());
        } else {
          throw std::runtime_error(
              "Tag \"Matrix\" only supports 1- and 2-dimensional tensors");
        }
 
        std::string values = std::regex_replace(
            node.text().get(), std::regex("[\t\r\n\a]+| +"), " ");
 
        auto [tensor_cpu, accessor] =
            to_tensorAccessor<T, N>(tensor, torch::kCPU);
        auto value = strtok(&values[0], " ");
 
        if constexpr (N == 1) {
          for (int64_t i = 0; i < rows; ++i) {
            if (value == nullptr)
              throw std::runtime_error(
                  "XML object does not provide enough coefficients");
            accessor[i] = static_cast<T>(std::stod(value));
            value = strtok(nullptr, " ");
          }
        } else if constexpr (N == 2) {
          for (int64_t i = 0; i < rows; ++i)
            for (int64_t j = 0; j < cols; ++j) {
              if (value == nullptr)
                throw std::runtime_error(
                    "XML object does not provide enough coefficients");
              accessor[i][j] = static_cast<T>(std::stod(value));
              value = strtok(nullptr, " ");
            }
        }
 
        if (value != nullptr)
          throw std::runtime_error("XML object provides too many coefficients");
 
        if (tensor.device().type() != torch::kCPU)
          tensor = std::move(tensor_cpu);
 
        return tensor;
 
      } else {
 
        // Check for "Dimensions"
        if (pugi::xml_node dims = node.child("Dimensions")) {
          std::vector<int64_t> sizes;
 
          std::string values = std::regex_replace(
              dims.text().get(), std::regex("[\t\r\n\a]+| +"), " ");
          for (auto value = strtok(&values[0], " "); value != nullptr;
               value = strtok(nullptr, " "))
            sizes.push_back(static_cast<std::size_t>(std::stoi(value)));
 
          if (!alloc && (tensor.sizes() != sizes))
            throw std::runtime_error("Invalid tensor dimensions");
 
          else if (alloc && (tensor.sizes() != sizes))
            tensor = torch::zeros(torch::IntArrayRef{sizes}, tensor.options());
 
          if (sizes.size() != N)
            throw std::runtime_error("Invalid tensor dimensions");
 
          // Check for "Data"
          if (pugi::xml_node data = node.child("Data")) {
            std::string values = std::regex_replace(
                data.text().get(), std::regex("[\t\r\n\a]+| +"), " ");
 
            auto [tensor_cpu, accessor] =
                to_tensorAccessor<T, N>(tensor, torch::kCPU);
            auto value = strtok(&values[0], " ");
 
            if constexpr (N == 1) {
              for (int64_t i = 0; i < sizes[0]; ++i) {
                if (value == nullptr)
                  throw std::runtime_error(
                      "XML object does not provide enough coefficients");
 
                accessor[i] = static_cast<T>(std::stod(value));
                value = strtok(nullptr, " ");
              }
            } else if constexpr (N == 2) {
              for (int64_t i = 0; i < sizes[0]; ++i)
                for (int64_t j = 0; j < sizes[1]; ++j) {
                  if (value == nullptr)
                    throw std::runtime_error(
                        "XML object does not provide enough coefficients");
 
                  accessor[i][j] = static_cast<T>(std::stod(value));
                  value = strtok(nullptr, " ");
                }
            } else if constexpr (N == 3) {
              for (int64_t i = 0; i < sizes[0]; ++i)
                for (int64_t j = 0; j < sizes[1]; ++j)
                  for (int64_t k = 0; k < sizes[2]; ++k) {
                    if (value == nullptr)
                      throw std::runtime_error(
                          "XML object does not provide enough coefficients");
 
                    accessor[i][j][k] = static_cast<T>(std::stod(value));
                    value = strtok(nullptr, " ");
                  }
            } else if constexpr (N == 4) {
              for (int64_t i = 0; i < sizes[0]; ++i)
                for (int64_t j = 0; j < sizes[1]; ++j)
                  for (int64_t k = 0; k < sizes[2]; ++k)
                    for (int64_t l = 0; l < sizes[3]; ++l) {
                      if (value == nullptr)
                        throw std::runtime_error(
                            "XML object does not provide enough coefficients");
 
                      accessor[i][j][k][l] = static_cast<T>(std::stod(value));
                      value = strtok(nullptr, " ");
                    }
            } else if constexpr (N == 5) {
              for (int64_t i = 0; i < sizes[0]; ++i)
                for (int64_t j = 0; j < sizes[1]; ++j)
                  for (int64_t k = 0; k < sizes[2]; ++k)
                    for (int64_t l = 0; l < sizes[3]; ++l)
                      for (int64_t m = 0; m < sizes[4]; ++m) {
                        if (value == nullptr)
                          throw std::runtime_error(
                              "XML object does not provide enough "
                              "coefficients");
 
                        accessor[i][j][k][l][m] =
                            static_cast<T>(std::stod(value));
                        value = strtok(nullptr, " ");
                      }
            } else if constexpr (N == 6) {
              for (int64_t i = 0; i < sizes[0]; ++i)
                for (int64_t j = 0; j < sizes[1]; ++j)
                  for (int64_t k = 0; k < sizes[2]; ++k)
                    for (int64_t l = 0; l < sizes[3]; ++l)
                      for (int64_t m = 0; m < sizes[4]; ++m)
                        for (int64_t n = 0; n < sizes[5]; ++n) {
                          if (value == nullptr)
                            throw std::runtime_error(
                                "XML object does not provide enough "
                                "coefficients");
 
                          accessor[i][j][k][l][m][n] =
                              static_cast<T>(std::stod(value));
                          value = strtok(nullptr, " ");
                        }
            }
 
            if (value != nullptr)
              throw std::runtime_error(
                  "XML object provides too many coefficients");
 
            if (tensor.device().type() != torch::kCPU)
              tensor = std::move(tensor_cpu);
 
            return tensor;
          } // "Data"
        }   // "Dimenions"
 
        throw std::runtime_error(
            "XML object does not provide a \"Dimensions\" tag");
 
        return tensor;
      }
 
    } // try next node
  }   // "tag"
 
  throw std::runtime_error(
      "XML object does not provide tag with given id, index, and/or label");
  return tensor;
}
 
template <typename T, std::size_t N, std::size_t M>
inline utils::TensorArray<M> &
from_xml(const pugi::xml_node &root, utils::TensorArray<M> &tensors,
         std::string tag = "Matrix", int id = 0, bool alloc = true,
         const std::string &label = "");
 
template <typename T, std::size_t N, std::size_t M>
inline utils::TensorArray<M> &
from_xml(const pugi::xml_document &doc, utils::TensorArray<M> &tensors,
         std::string tag = "Matrix", int id = 0, bool alloc = true,
         const std::string &label = "") {
 
  return from_xml<T, N>(doc.child("xml"), tensors, tag, id, alloc, label);
}
 
template <typename T, std::size_t N, std::size_t M>
inline utils::TensorArray<M> &
from_xml(const pugi::xml_node &root, utils::TensorArray<M> &tensors,
         std::string tag, int id, bool alloc, const std::string &label) {
 
  for (std::size_t i = 0; i < M; ++i) {
    from_xml<T, N>(root, tensors[i], tag, id, label, alloc, i);
  }
 
  return tensors;
}
 
} // namespace iganet::utils