Hdf5File.cpp

Go to the documentation of this file.

/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
// Code only included on platforms with libhdf5
// Dummy placeholders defined in Hdf5File.h will be used if libhdf5 is
// unavailable on compile time

#include "Hdf5File.h"

#include "ptree2hdf5.h"
#include "bio/AttributeContainer.h"
#include "bio/MeshState.h"
#include "sim/Sim.h"
#include "sim/SimState.h"
#include "util/misc/log_debug.h"
#include "util/misc/XmlWriterSettings.h"

#include <boost/algorithm/string/predicate.hpp>
#include <boost/property_tree/ptree.hpp>
#include <boost/property_tree/xml_parser.hpp>

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>

using namespace std;
using namespace boost::property_tree;
using namespace boost::property_tree::xml_parser;
using namespace SimPT_Sim::Util;

namespace {
    // TODO: Document!
    bool write_array_rank1(hid_t loc_id, string name, void const* data, hid_t type_id, hsize_t dim_1);
    bool write_array_rank1(hid_t loc_id, string name, const vector<string>& data);
    bool write_array_rank2(hid_t loc_id, string name, void const* data, hid_t type_id, hsize_t dim_1, hsize_t dim_2);
    // Not used for now
    // bool write_array_rank3(hid_t loc_id, string name, void const* data, hid_t type_id, hsize_t dim_1, hsize_t dim_2, hsize_t dim_3);

    template <typename T>
    bool read_array_rank1(hid_t loc_id, string name, T** data, hid_t type_id, size_t* dim_1, bool optional = false);

    template <typename T>
    bool read_array_rank2(hid_t loc_id, string name, T** data, hid_t type_id, size_t* dim_1, size_t* dim_2, bool optional = false);

    list<string> find_matching_links(hid_t loc_id, function<bool(string)> match);

}

namespace {
    bool write_array_rank1(hid_t loc_id, string name, void const* data, hid_t type_id, hsize_t dim_1) {
        // No writing if something called 'name' already exists in loc_id
        if (0 == H5Lexists(loc_id, name.c_str(), H5P_DEFAULT)) {
            // Create rank 1 dataspace
            hid_t dspace = H5Screate_simple(1, &dim_1, 0);
            // Create dataset 'name' of provided simple type using the rank 1 dataspace created before
            hid_t dset = H5Dcreate(loc_id, name.c_str(), type_id, dspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            // Write data using the native datatype corresponding to the provided type_id
            hid_t mem_type_id = H5Tget_native_type(type_id, H5T_DIR_ASCEND);
            H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
            // Clean up HDF5 handles
            H5Tclose(mem_type_id);
            H5Dclose(dset);
            H5Sclose(dspace);
            return true;
        } else {
            return false;
        }
    }

    bool write_array_rank1(hid_t loc_id, string name, const vector<string>& data) {
        // No writing if something called 'name' already exists in loc_id
        if (0 == H5Lexists(loc_id, name.c_str(), H5P_DEFAULT)) {
            // Use C-style string datatype
            hid_t dtype = H5Tcopy(H5T_C_S1);
            H5Tset_size(dtype, H5T_VARIABLE);
            // Create rank 1 dataspace
            hsize_t dim_1 = data.size();
            hid_t dspace = H5Screate_simple(1, &dim_1, 0);
            // Create dataset 'name' of provided simple type using the rank 1 dataspace created before
            hid_t dset = H5Dcreate(loc_id, name.c_str(), dtype, dspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            // Write data as an array of pointers to the C-string representation of the strings
            vector<const char*> raw_string_ptrs;
            for (const string & s : data) raw_string_ptrs.push_back(s.c_str());
            H5Dwrite(dset, dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, raw_string_ptrs.data());
            // Clean up HDF5 handles
            H5Dclose(dset);
            H5Sclose(dspace);
            H5Tclose(dtype);
            return true;
        } else {
            return false;
        }
    }

    bool write_array_rank2(hid_t loc_id, string name, void const* data, hid_t type_id, hsize_t dim_1, hsize_t dim_2) {
        // No writing if something called 'name' already exists in loc_id
        if (0 == H5Lexists(loc_id, name.c_str(), H5P_DEFAULT)) {
            // Create rank 2 dataspace
            hsize_t dims[2] = {dim_1, dim_2};
            hid_t dspace = H5Screate_simple(2, dims, 0);
            // Create dataset 'name' of provided simple type using the rank 2 dataspace created before
            hid_t dset = H5Dcreate(loc_id, name.c_str(), type_id, dspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            // Write data
            hid_t mem_type_id = H5Tget_native_type(type_id, H5T_DIR_ASCEND);
            H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
            // Clean up HDF5 handles
            H5Tclose(mem_type_id);
            H5Dclose(dset);
            H5Sclose(dspace);
            return true;
        } else {
            return false;
        }
    }

    // Not used for now
    /*
    bool write_array_rank3(hid_t loc_id, string name, void const* data, hid_t type_id, hsize_t dim_1, hsize_t dim_2, hsize_t dim_3) {
        // No writing if something called 'name' already exists in loc_id
        if (0 == H5Lexists(loc_id, name.c_str(), H5P_DEFAULT)) {
            // Create rank 3 dataspace
            hsize_t dims[3] = {dim_1, dim_2, dim_3};
            hid_t dspace = H5Screate_simple(3, dims, 0);
            // Create dataset 'name' of provided simple type using the rank 3 dataspace created before
            hid_t dset = H5Dcreate(loc_id, name.c_str(), type_id, dspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            // Write data
            hid_t mem_type_id = H5Tget_native_type(type_id, H5T_DIR_ASCEND);
            H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
            // Clean up HDF5 handles
            H5Tclose(mem_type_id);
            H5Dclose(dset);
            H5Sclose(dspace);
            return true;
        } else {
            return false;
        }
    }
    */

    template <typename T>
    bool read_array_rank1(hid_t loc_id, string name, T** data, hid_t mem_type_id, size_t* dim_1, bool optional) {
        // HDF5 error status. Used to check if bad things are happening.
        herr_t h5_status;

        // Quickly check for existence of a link called 'name'
        if (0 == H5Lexists(loc_id, name.c_str(), H5P_DEFAULT)) {
            if (optional) {
                return false;
            } else {
                string errmsg = string("No dataset \'") + name + string("\'");
                throw runtime_error(errmsg);
            }
        }

        // Try to open the dataset
        hid_t dset = H5Dopen(loc_id, name.c_str(), H5P_DEFAULT);
        if (dset < 0) {
            if (optional) {
                return false;
            } else {
                string errmsg = string("Could not open dataset \'") + name + string("\'");
                throw runtime_error(errmsg);
            }
        }

        // Get the associated dataspace
        hid_t dspace = H5Dget_space(dset);
        int ndims = H5Sget_simple_extent_ndims(dspace);
        if (ndims != 1) {
            string errmsg = string("Dataset \'") + name + string("\' must be rank 1");
            throw runtime_error(errmsg);
        }

        // Number of elements in data
        hsize_t _dim_1;
        H5Sget_simple_extent_dims(dspace, &_dim_1, 0);
        *dim_1 = static_cast<size_t>(_dim_1);

        // Allocate space for data
        *data = new T[_dim_1];

        // Read the data
        h5_status = H5Dread(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, *data);
        if (h5_status < 0) {
            string errmsg = string("Could not read dataset \'") + name + string("\'");
            throw runtime_error(errmsg);
        }

        // TODO: is type-checking feasible?

        // Clean up HDF5 handles
        H5Sclose(dspace);
        H5Dclose(dset);
        return true;
    }

    template <typename T>
    bool read_array_rank2(hid_t loc_id, string name, T** data, hid_t mem_type_id, size_t* dim_1, size_t* dim_2, bool optional) {
        // HDF5 error status. Used to check if bad things are happening.
        herr_t h5_status;

        // Quickly check for existence of a link called 'name'
        if (0 == H5Lexists(loc_id, name.c_str(), H5P_DEFAULT)) {
            if (optional) {
                return false;
            } else {
                string errmsg = string("No dataset \'") + name + string("\'");
                throw runtime_error(errmsg);
            }
        }

        // Try to open the dataset
        hid_t dset = H5Dopen(loc_id, name.c_str(), H5P_DEFAULT);
        if (dset < 0) {
            if (optional) {
                return false;
            } else {
                string errmsg = string("Could not open dataset \'") + name + string("\'");
                throw runtime_error(errmsg);
            }
        }

        // Get the associated dataspace
        hid_t dspace = H5Dget_space(dset);
        int ndims = H5Sget_simple_extent_ndims(dspace);
        if (ndims != 2) {
            string errmsg = string("Dataset \'") + name + string("\' must be rank 2");
            throw runtime_error(errmsg);
        }

        // Number of elements in data
        hsize_t _dims[2];
        H5Sget_simple_extent_dims(dspace, _dims, 0);
        *dim_1 = static_cast<size_t>(_dims[0]);
        *dim_2 = static_cast<size_t>(_dims[1]);

        // Allocate space for data
        *data = new T[_dims[0] * _dims[1]];

        // Read the data
        h5_status = H5Dread(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, *data);
        if (h5_status < 0) {
            string errmsg = string("Could not read dataset \'") + name + string("\'");
            throw runtime_error(errmsg);
        }

        // TODO: is type-checking feasible?

        // Clean up HDF5 handles
        H5Sclose(dspace);
        H5Dclose(dset);
        return true;
    }

    //herr_t (*H5L_iterate_t)( hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)
    herr_t collect_link_names(hid_t , const char* name, const H5L_info_t* , void* op_data) {
        list<string>* link_names = static_cast<list<string>*>(op_data);
        link_names->push_back(string(name));
        return 0;
    }

    list<string> find_matching_links(hid_t loc_id, function<bool(string)> match) {
        // Collect all link names at location loc_id
        list<string> link_names;
        H5Literate(loc_id, H5_INDEX_NAME, H5_ITER_NATIVE, 0, &collect_link_names, &link_names);
        // Remove link names that DON'T match the 'match' predicate
        link_names.remove_if([&](string name)->bool{ return !match(name); });
        return link_names;
    }
}

namespace SimPT_Shell {

using namespace std;
using boost::property_tree::ptree;
using SimPT_Sim::Sim;
using SimPT_Sim::SimState;
using SimPT_Sim::MeshState;

// Note that a hid_t identifier -1 for m_file_id represents a closed file
Hdf5File::Hdf5File()
    : m_file_id(-1), m_file_path("") {}

Hdf5File::Hdf5File(const string& file_path)
    : m_file_id(-1), m_file_path("")
{
    // Open resets m_file_id and m_file_path iff successful
    Open(file_path);
}

Hdf5File::~Hdf5File()
{
    if (IsOpen()) {
        Close();
    }
}

bool Hdf5File::Close()
{
    bool status = false;
    if (m_file_id >= 0) {
        H5Fclose(m_file_id);
        m_file_path = "";
        status = true;
    }
    return status;
}

string Hdf5File::GetFileExtension()
{
    return "h5";
}

string Hdf5File::GetFilePath()
{
    return m_file_path;
}

bool Hdf5File::IsOpen()
{
    return (m_file_path != "") && (m_file_id != -1);
}

bool Hdf5File::Open(const string& file_path)
{
    string const here = string(VL_HERE) + "> ";
    // Did we succeed in providing a valid HDF5 file, either by opening an
    // existing one or by creating a new one?
    bool successful = true;
    bool new_file = false;

    // Close the storage if any was open
    if (IsOpen()) Close(); 

    // Open or create the file if file name is not empty
    if (file_path != "") {
        // Try to open the file if it exists and is an HDF5 file
        // Save old error handler and turn off error handling temporarily since this function
        // fail and will clutter up the stdout.
        herr_t (*old_func)(hid_t, void*);
        void *old_client_data;
        H5Eget_auto(H5E_DEFAULT, &old_func, &old_client_data);
        H5Eset_auto(H5E_DEFAULT, NULL, NULL);
        htri_t is_hdf5 = H5Fis_hdf5(file_path.c_str());
        // Restore previous error handler
        H5Eset_auto(H5E_DEFAULT, old_func, old_client_data);
        if (0 < is_hdf5) {
            // File exists and is in HDF5 format, open it
            m_file_id = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
            if (m_file_id < 0) {
                // File exists and is HDF5 but can not be opened, pretty strange: bail out!
                cerr << here << "[Hdf5File] Error: \'" << file_path << "\' seems to be valid HDF5 but could not be opened.\
                    Might be you don\'t have r/w permissions." << endl;
                successful = false;
            } else {
                successful = true;
                new_file = false;
                m_file_path = file_path;
            }
        } else if (0 == is_hdf5) {
            // File exists but is not in HDF5 format: bail out!
            cerr << here << "[Hdf5File] Error: \'" << file_path << "\' is not a valid HDF5 file.\
                I can neither open it, nor will I truncate it." << endl;
            successful = false;
        } else {
            // File does not exist, create a new one
            m_file_id = H5Fcreate(file_path.c_str(), H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT);
            if (m_file_id < 0) {
                // Probably no permissions to create a file
                cerr << here << "[Hdf5File] Error: \'" << file_path << "\' could not be created nor opened.\
                    Check your r/w permissions." << endl;
                successful = false;
            } else {
                successful = true;
                new_file = true;
                m_file_path = file_path;
            }
        }
    }

    // Populate the new HDF5 file with the basic skeleton to keep time-series of simulation steps
    //  - /time_steps
    //  - /time_steps_idx
    if (successful && new_file) {
        // Create two extendable 1D double precision and integer array datasets
        // to hold the simulation time steps in seconds (double precision floats)
        // to hold the simulation time steps indices (integers)
        // This data set must be extendable to allow appending time steps during the simulation
        hsize_t time_steps_dims = 0;                        // Initial size: no elements
        hsize_t time_steps_max_dims = H5S_UNLIMITED;        // Can extend with no bounds

        // Create the data space, first argument 'rank' is 1 (i.e. 1 dimensional array)
        hid_t time_steps_dspace = H5Screate_simple(1, &time_steps_dims, &time_steps_max_dims);

        // Create properties for dataset creation: chunking is necessary for extendable data sets
        hid_t time_steps_dset_props = H5Pcreate(H5P_DATASET_CREATE);
        hsize_t time_steps_chunk_dims = 10;                                 // Rather arbitrary chunk size
        H5Pset_chunk(time_steps_dset_props, 1, &time_steps_chunk_dims);     // 1 is the rank of a 1D array

        // Create data set '/time_steps' (link creation and data access properties are default)
        hid_t time_steps_dset = H5Dcreate(m_file_id, "/time_steps", H5T_IEEE_F64LE,
                time_steps_dspace, H5P_DEFAULT, time_steps_dset_props, H5P_DEFAULT);
        if (time_steps_dset < 0) {
            cerr << here << "[Hdf5File] Error: could not create \'/time_steps\' dataset." << endl;
            successful = false;
        }

        // Create data set '/time_steps_idx' (link creation and data access properties are default)
        hid_t time_steps_idx_dset = H5Dcreate(m_file_id, "/time_steps_idx", H5T_STD_I32LE,
                time_steps_dspace, H5P_DEFAULT, time_steps_dset_props, H5P_DEFAULT);
        if (time_steps_idx_dset < 0) {
            cerr << here << "[Hdf5File] Error: could not create \'/time_steps_idx\' dataset." << endl;
            successful = false;
        }

        // Clean up
        H5Dclose(time_steps_dset);
        H5Dclose(time_steps_idx_dset);
        H5Sclose(time_steps_dspace);
        H5Pclose(time_steps_dset_props);
    }

    return successful;
}

SimState Hdf5File::Read(int step_idx)
{
    if (!IsOpen()) {
        throw runtime_error("[Hdf5File]: File needs to be open prior to reading.");
    }

    // The Sim/MeshState that will be filled up with information form HDF5 and returned
    MeshState mstate;
    SimState sstate;

    // =========================================================================
    // Read information about time steps contained in the file
    double* time_steps_data = nullptr;
    size_t time_steps_dim1 = 0;
    read_array_rank1(m_file_id, "time_steps",  &time_steps_data, H5T_NATIVE_DOUBLE, &time_steps_dim1);

    int* time_steps_idx_data = nullptr;
    size_t time_steps_idx_dim1 = 0;
    read_array_rank1(m_file_id, "time_steps_idx", &time_steps_idx_data, H5T_NATIVE_INT, &time_steps_idx_dim1);

    // Check size consistency of /time_steps and /time_step_idx
    if (time_steps_dim1 != time_steps_idx_dim1) {
        string errmsg = string("Error reading file \'") + m_file_path
                + string("\': The /time_steps and /time_steps_idx datasets must be the same size. File might be corrupted.");
        // Clean up temp storage
        delete[] time_steps_data;
        delete[] time_steps_idx_data;
        throw runtime_error(errmsg);
    }

    if (-1 == step_idx) {
        // Use the last time step if no step was requested
        step_idx = time_steps_idx_dim1 - 1;
    }

    const size_t sim_step = time_steps_idx_data[step_idx];
    const double sim_time = time_steps_data[step_idx];

    // Try to open the requested time step
    string step_grp_name = string("/step_") + to_string(sim_step);
    hid_t step_grp_id = H5Gopen(m_file_id, step_grp_name.c_str(), H5P_DEFAULT);
    if (step_grp_id < 0) {
        string errmsg = string("Error reading file \'") + m_file_path
                + string("\': The group /step_") + to_string(sim_step)
                + string(" does not exist. File might be corrupted.");
        // Clean up temp storage
        delete[] time_steps_data;
        delete[] time_steps_idx_data;
        throw runtime_error(errmsg);
    }

    // Clean up temp storage
    delete[] time_steps_data;
    delete[] time_steps_idx_data;

    // =========================================================================

    ReadNodes(step_grp_id, mstate);
    ReadCells(step_grp_id, mstate);
    ReadBPNodes(step_grp_id, mstate);
    // In the very special case of just one cell (e.g. in the beginning of
    // a simulation, there's no walls. Hence: don't read walls in that case)
    if (1 < mstate.GetNumCells()) {
        ReadWalls(step_grp_id, mstate);
        ReadBPWalls(step_grp_id, mstate);
    }
    //ptree params = ReadParameters(step_grp_id);
    ptree params = ReadPtreeFromStringAttribute(step_grp_id, "parameters");
    ptree re_state = ReadPtreeFromStringAttribute(step_grp_id, "random_engine");

    sstate.SetTimeStep(sim_step);
    sstate.SetTime(sim_time);
    sstate.SetParameters(params);
    sstate.SetRandomEngineState(re_state);
    sstate.SetMeshState(mstate);

    return sstate;
}

void Hdf5File::ReadNodes(hid_t loc_id, MeshState& mstate)
{
    // Reading nodes_id
    int* nodes_id_data = nullptr;
    size_t nodes_id_dim1 = 0;
    read_array_rank1(loc_id, "nodes_id", &nodes_id_data, H5T_NATIVE_INT, &nodes_id_dim1);
    // Will be used as reference size to validate other sizes
    const size_t& num_nodes = nodes_id_dim1;

    // Reading nodes_xy
    double* nodes_xy_data = nullptr;
    size_t nodes_xy_dim1 = 0;
    size_t nodes_xy_dim2 = 0;
    read_array_rank2(loc_id, "nodes_xy", &nodes_xy_data, H5T_NATIVE_DOUBLE, &nodes_xy_dim1, &nodes_xy_dim2);

    // Verify dimensions
    if (2 != nodes_xy_dim2) {
        string errmsg = string("Error reading file \'") + m_file_path
            + string("\': Second dimension of the /step_*/nodes_xy dataset must equal 2. File might be corrupted.");
        // Clean up temp storage
        delete[] nodes_id_data;
        delete[] nodes_xy_data;
        throw runtime_error(errmsg);
    }
    if (num_nodes != nodes_xy_dim1) {
        string errmsg = string("Error reading file \'") + m_file_path
            + string("\': Dimensions of \'nodes_xy\' and \'nodes_id\' datasets don\'t match. File might be corrupted.");
        // Clean up temp storage
        delete[] nodes_id_data;
        delete[] nodes_xy_data;
        throw runtime_error(errmsg);
    }

    // Build nodes from data read before
    mstate.SetNumNodes(num_nodes);
    for (size_t node_idx = 0; node_idx < num_nodes; ++node_idx) {
            assert( node_idx == static_cast<unsigned int>(nodes_id_data[node_idx]) && "Error for node ids");
        mstate.SetNodeID(node_idx, nodes_id_data[node_idx]);
        mstate.SetNodeX(node_idx, nodes_xy_data[2 * node_idx]);
        mstate.SetNodeY(node_idx, nodes_xy_data[2 * node_idx + 1]);
    }

    // Clean up temp storage
    delete[] nodes_id_data;
    delete[] nodes_xy_data;
    
    // Find datasets that are called 'nodes_attr_<name>' (I.e.: node attributes)
    list<string> attr_names = find_matching_links(loc_id, [](string name)->bool{
        return boost::starts_with(name, "nodes_attr_");
    });

    // Try to read all matching node attribute datasets
    for (string long_name : attr_names) {
        // Magic number 11 is the length of the prefix "nodes_attr_" that needs top be stripped away
        string short_name = long_name.substr(11, string::npos);

        // Check if "long_name" refers to a valid rank-1 dataset in loc_id
        H5O_info_t obj_info;
        H5Oget_info_by_name(loc_id, long_name.c_str(), &obj_info, H5P_DEFAULT);
        if (H5O_TYPE_DATASET == obj_info.type) {
            // Open the dataset to get all info
            hid_t dset = H5Dopen(loc_id, long_name.c_str(), H5P_DEFAULT);

            // Make sure dimensions match
            hid_t dspace = H5Dget_space(dset);
            // Check if rank = 1
            const int ndims = H5Sget_simple_extent_ndims(dspace);
            if (1 != ndims) {
                H5Sclose(dspace);
                break;
            }
            // Check if number of attribute values = number of nodes
            hsize_t dims = 0;
            H5Sget_simple_extent_dims(dspace, &dims, 0);
            if (num_nodes != dims) {
                H5Sclose(dspace);
                break;
            }

            // Get the type of the attribute from HDF5 (either int or double)
            hid_t dtype = H5Dget_type(dset);
            H5T_class_t dtype_class = H5Tget_class(dtype);
            H5Tclose(dtype);
            if (H5T_INTEGER == dtype_class) {
                // Read the data directly into a vector
                vector<int> attr_values(num_nodes);
                H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, attr_values.data());
                // Add attribute and its values to the nodes
                mstate.NodeAttributeContainer().Add<int>(short_name);
                mstate.NodeAttributeContainer().SetAll<int>(short_name, attr_values);
            } else if (H5T_FLOAT == dtype_class) {
                // Read the data directly into a vector
                vector<double> attr_values(num_nodes);
                H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, attr_values.data());
                // Add attribute and its values to the nodes
                mstate.NodeAttributeContainer().Add<double>(short_name);
                mstate.NodeAttributeContainer().SetAll<double>(short_name, attr_values);
            } else if (H5T_STRING == dtype_class) {
                // Read the raw string data directly into a vector of C-style strings
                vector<char*> raw_attr_values(num_nodes);
                hid_t ntype = H5Tcopy(H5T_C_S1);
                H5Tset_size(ntype, H5T_VARIABLE);
                H5Dread(dset, ntype, H5S_ALL, H5S_ALL, H5P_DEFAULT, raw_attr_values.data());
                H5Tclose(ntype);

                // Convert the C-style strings to a C++ string vector
                vector<string> attr_values(raw_attr_values.begin(), raw_attr_values.end());

                // Add attribute and its values to the nodes
                mstate.NodeAttributeContainer().Add<string>(short_name);
                mstate.NodeAttributeContainer().SetAll<string>(short_name, attr_values);
            } else {
                // Unsupported datatype: ignore (or generate a warning)
            }

            // Clean up
            H5Dclose(dset);
        }
    }
}

void Hdf5File::ReadCells(hid_t loc_id, MeshState& mstate)
{
        // Reading 'cells_id'
        int* cells_id_data = nullptr;
        size_t cells_id_dim1 = 0;
        read_array_rank1(loc_id, "cells_id", &cells_id_data, H5T_NATIVE_INT, &cells_id_dim1);
        // Will be used as reference size to validate other sizes
        const size_t & num_cells = cells_id_dim1;

        // Reading 'cells_num_nodes'
        int* cells_num_nodes_data = nullptr;
        size_t cells_num_nodes_dim1 = 0;
        read_array_rank1(loc_id, "cells_num_nodes", &cells_num_nodes_data, H5T_NATIVE_INT, &cells_num_nodes_dim1);

        // Reading 'cells_nodes'
        int* cells_nodes_data = nullptr;
        size_t cells_nodes_dim1 = 0, cells_nodes_dim2 = 0;
        read_array_rank2(loc_id, "cells_nodes", &cells_nodes_data, H5T_NATIVE_INT, &cells_nodes_dim1, &cells_nodes_dim2);
        const size_t & max_num_cell_nodes = cells_nodes_dim2;

        // Reading 'cells_num_walls'
        int* cells_num_walls_data = nullptr;
        size_t cells_num_walls_dim1 = 0;
        read_array_rank1(loc_id, "cells_num_walls", &cells_num_walls_data, H5T_NATIVE_INT, &cells_num_walls_dim1);

        // Reading 'cells_walls'
        // (Don't read 'cells_walls' in the special case of one cell, which has no walls)
        int* cells_walls_data = nullptr;
        size_t cells_walls_dim1 = 0, cells_walls_dim2 = 0;
        if (1 != num_cells) {
                read_array_rank2(loc_id, "cells_walls", &cells_walls_data, H5T_NATIVE_INT, &cells_walls_dim1, &cells_walls_dim2);
        }
        const size_t & max_num_cell_walls = cells_walls_dim2;
    
        // Verify dimensions
        if (num_cells != cells_num_nodes_dim1) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': Dimensions of \'cells_num_nodes\' and \'cells_id\' datasets don\'t match. File might be corrupted.");
                // Clean up temp storage
                delete[] cells_id_data;
                delete[] cells_num_nodes_data;
                delete[] cells_nodes_data;
                delete[] cells_num_walls_data;
                delete[] cells_walls_data;
                throw runtime_error(errmsg);
        }
        if (num_cells != cells_nodes_dim1) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': Dimensions of \'cells_nodes\' and \'cells_id\' datasets don\'t match. File might be corrupted.");
                // Clean up temp storage
                delete[] cells_id_data;
                delete[] cells_num_nodes_data;
                delete[] cells_nodes_data;
                delete[] cells_num_walls_data;
                delete[] cells_walls_data;
                throw runtime_error(errmsg);
        }
        if (num_cells != cells_num_walls_dim1) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': Dimensions of \'cells_num_walls\' and \'cells_id\' datasets don\'t match. File might be corrupted.");
                // Clean up temp storage
                delete[] cells_id_data;
                delete[] cells_num_nodes_data;
                delete[] cells_nodes_data;
                delete[] cells_num_walls_data;
                delete[] cells_walls_data;
                throw runtime_error(errmsg);
        }
        // Don't pay attention to anything related to 'cells_walls'
        // in the special case of one cell, which has no walls
        if (1 != num_cells && num_cells != cells_walls_dim1) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': Dimensions of \'cells_walls\' and \'cells_id\' datasets don\'t match. File might be corrupted.");
                // Clean up temp storage
                delete[] cells_id_data;
                delete[] cells_num_nodes_data;
                delete[] cells_nodes_data;
                delete[] cells_num_walls_data;
                delete[] cells_walls_data;
                throw runtime_error(errmsg);
        }

        // Build cells from data read before (new API)
        mstate.SetNumCells(num_cells);
        for (size_t cell_idx = 0; cell_idx < num_cells; ++cell_idx) {
                // ID of the cell
                assert(cell_idx == static_cast<unsigned int>(cells_id_data[cell_idx]) && "Error for cell ids");
                mstate.SetCellID(cell_idx, cells_id_data[cell_idx]);
                // List of node IDs that make up this cell's polygon
                vector<int> cell_nodes;
                for (size_t cell_node_idx = 0; cell_node_idx < static_cast<size_t>(cells_num_nodes_data[cell_idx]);
                        ++cell_node_idx) {
                        // Data from 2D HDF5 dataset '/step_n/cells_nodes' is stored row-wise
                        // in the C++ array called 'cells_nodes_data'. This means that:
                        //  - first dimension is: cell node ids
                        //  - second dimension is: cells
                        // and hence the index is:
                        size_t idx = cell_node_idx + (cell_idx * max_num_cell_nodes);
                        cell_nodes.push_back(cells_nodes_data[idx]);
                }
                // Set the current cell's nodes in MeshState
                mstate.SetCellNodes(cell_idx, cell_nodes);
                // List of wall IDs that connect this cell to its neighbors
                vector<int> cell_walls;
                for (size_t cell_wall_idx = 0; cell_wall_idx < static_cast<size_t>(cells_num_walls_data[cell_idx]);
                        ++cell_wall_idx) {
                        // Data from 2D HDF5 dataset '/step_n/cells_walls' is stored row-wise
                        // in the C++ array called 'cells_walls_data'. This means that:
                        //  - first dimension is: cell wall ids
                        //  - second dimension is: cells
                        // and hence the index is:
                        size_t idx = cell_wall_idx + (cell_idx * max_num_cell_walls);
                        cell_walls.push_back(cells_walls_data[idx]);
                }
                // Set the current cell's walls in MeshState
                mstate.SetCellWalls(cell_idx, cell_walls);
        }

        // Clean up temp storage
        delete[] cells_id_data;
        delete[] cells_num_nodes_data;
        delete[] cells_nodes_data;
        delete[] cells_num_walls_data;
        delete[] cells_walls_data;

        // Find datasets that are called 'cells_attr_<name>' (I.e.: cell attributes)
        list<string> attr_names = find_matching_links(loc_id, [](string name)->bool {
                return boost::starts_with(name, "cells_attr_");
        });

        // Try to read all matching cell attribute datasets
        for (string long_name : attr_names) {
                // Magic number 11 is the length of the prefix "cells_attr_" that needs top be stripped away
                string short_name = long_name.substr(11, string::npos);

                // Check if "long_name" refers to a valid rank-1 dataset in loc_id
                H5O_info_t obj_info;
                H5Oget_info_by_name(loc_id, long_name.c_str(), &obj_info, H5P_DEFAULT);
                if (H5O_TYPE_DATASET == obj_info.type) {
                        // Open the dataset to get all info
                        hid_t dset = H5Dopen(loc_id, long_name.c_str(), H5P_DEFAULT);

                        // Make sure dimensions match
                        hid_t dspace = H5Dget_space(dset);
                        // Check if rank = 1
                        const int ndims = H5Sget_simple_extent_ndims(dspace);
                        if (1 != ndims) {
                                H5Sclose(dspace);
                                break;
                        }
                        // Check if number of attribute values = number of cells
                        hsize_t dims = 0;
                        H5Sget_simple_extent_dims(dspace, &dims, 0);
                        if (num_cells != dims) {
                                H5Sclose(dspace);
                                break;
                        }

                        // Get the type of the attribute from HDF5 (either int or double)
                        hid_t dtype = H5Dget_type(dset);
                        H5T_class_t dtype_class = H5Tget_class(dtype);
                        H5Tclose(dtype);
                        if (H5T_INTEGER == dtype_class) {
                                // Read the data directly into a vector
                                vector<int> attr_values(num_cells);
                                H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, attr_values.data());
                                // Add attribute and its values to the cells
                                mstate.CellAttributeContainer().Add<int>(short_name);
                                mstate.CellAttributeContainer().SetAll<int>(short_name, attr_values);
                        } else if (H5T_FLOAT == dtype_class) {
                                // Read the data directly into a vector
                                vector<double> attr_values(num_cells);
                                H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, attr_values.data());
                                // Add attribute and its values to the cells
                                mstate.CellAttributeContainer().Add<double>(short_name);
                                mstate.CellAttributeContainer().SetAll<double>(short_name, attr_values);
                        } else if (H5T_STRING == dtype_class) {
                                // Read the raw string data directly into a vector of C-style strings
                                vector<char*> raw_attr_values(num_cells);
                                hid_t ntype = H5Tcopy(H5T_C_S1);
                                H5Tset_size(ntype, H5T_VARIABLE);
                                H5Dread(dset, ntype, H5S_ALL, H5S_ALL, H5P_DEFAULT, raw_attr_values.data());
                                H5Tclose(ntype);

                                // Convert the C-style strings to a C++ string vector
                                vector < string > attr_values(raw_attr_values.begin(), raw_attr_values.end());

                                // Add attribute and its values to the cells
                                mstate.CellAttributeContainer().Add<string>(short_name);
                                mstate.CellAttributeContainer().SetAll<string>(short_name, attr_values);
                        } else {
                                // Unsupported datatype: ignore (or generate a warning)
                        }

                        // Clean up
                        H5Dclose(dset);
                }
        }
}

void Hdf5File::ReadBPNodes(hid_t loc_id, MeshState& mstate)
{
        // Reading 'bp_nodes'
        int*   bp_nodes_data = 0;
        size_t bp_nodes_dim1 = 0;
        read_array_rank1(loc_id, "bp_nodes", &bp_nodes_data, H5T_NATIVE_INT, &bp_nodes_dim1);

        // Copy to temporary vectors
        vector<int> bp_nodes(bp_nodes_dim1);
        for (size_t node_idx = 0; node_idx < bp_nodes_dim1; ++node_idx) {
                bp_nodes[node_idx] = bp_nodes_data[node_idx];
        }

        // Use temp vectors to set boundary polygon in the mesh state
        mstate.SetBoundaryPolygonNodes(bp_nodes);

        // Clean up temp storage
        delete[] bp_nodes_data;
}

void Hdf5File::ReadBPWalls(hid_t loc_id, MeshState& mstate)
{
        // Reading 'bp_walls'
        int*   bp_walls_data = nullptr;
        size_t bp_walls_dim1 = 0;
        read_array_rank1(loc_id, "bp_walls", &bp_walls_data, H5T_NATIVE_INT, &bp_walls_dim1);

        // Copy to temporary vectors
        vector<int> bp_walls(bp_walls_dim1);
        for (size_t wall_idx = 0; wall_idx < bp_walls_dim1; ++wall_idx) {
                bp_walls[wall_idx] = bp_walls_data[wall_idx];
        }

        // Use temp vectors to set boundary polygon in the mesh state
        mstate.SetBoundaryPolygonWalls(bp_walls);

        // Clean up temp storage
        delete[] bp_walls_data;
}

void Hdf5File::ReadWalls(hid_t loc_id, MeshState& mstate)
{
        // Reading 'walls_id'
        int* walls_id_data = nullptr;
        size_t walls_id_dim1 = 0;
        read_array_rank1(loc_id, "walls_id", &walls_id_data, H5T_NATIVE_INT, &walls_id_dim1);
        // Will be used as reference size to validate other sizes
        const size_t & num_walls = walls_id_dim1;

        // Reading 'walls_cells'
        int* walls_cells_data = nullptr;
        size_t walls_cells_dim1 = 0, walls_cells_dim2 = 0;
        read_array_rank2(loc_id, "walls_cells", &walls_cells_data, H5T_NATIVE_INT, &walls_cells_dim1,
                &walls_cells_dim2);

        // Reading 'walls_nodes'
        int* walls_nodes_data = nullptr;
        size_t walls_nodes_dim1 = 0, walls_nodes_dim2 = 0;
        read_array_rank2(loc_id, "walls_nodes", &walls_nodes_data, H5T_NATIVE_INT, &walls_nodes_dim1,
                &walls_nodes_dim2);

        // Verify dimensions
        if (2 != walls_cells_dim2) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': Second dimension of the /step_*/walls_cells dataset must equal 2. File might be corrupted.");
                // Clean up temp storage
                delete[] walls_id_data;
                delete[] walls_cells_data;
                delete[] walls_nodes_data;
                throw runtime_error(errmsg);
        }
        if (2 != walls_nodes_dim2) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': Second dimension of the /step_*/walls_nodes dataset must equal 2. File might be corrupted.");
                // Clean up temp storage
                delete[] walls_id_data;
                delete[] walls_cells_data;
                delete[] walls_nodes_data;
                throw runtime_error(errmsg);
        }
        if (num_walls != walls_cells_dim1) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': Dimensions of \'walls_id\' and \'walls_cells\' datasets don\'t match. File might be corrupted.");
                // Clean up temp storage
                delete[] walls_id_data;
                delete[] walls_cells_data;
                delete[] walls_nodes_data;
                throw runtime_error(errmsg);
        }
        if (num_walls != walls_nodes_dim1) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': Dimensions of \'walls_id\' and \'walls_nodes\' datasets don\'t match. File might be corrupted.");
                // Clean up temp storage
                delete[] walls_id_data;
                delete[] walls_cells_data;
                delete[] walls_nodes_data;
                throw runtime_error(errmsg);
        }

        // Build walls from data read before (new API)
        mstate.SetNumWalls(num_walls);
        for (size_t wall_idx = 0; wall_idx < num_walls; ++wall_idx) {
                // ID of the wall
                assert(wall_idx == static_cast<unsigned int>(walls_id_data[wall_idx]) && "Error for wall ids");
                mstate.SetWallID(wall_idx, walls_id_data[wall_idx]);
                // - IDs of the two cells the wall is connecting
                // - IDs of the two nodes that are the endpoints of the wall
                //   (the IDs of cells and nodes are stored row-major in C++:
                //    * first dimension is: IDs of cells / nodes
                //    * second dimension is: walls)
                mstate.SetWallCells(wall_idx,
                        make_pair(walls_cells_data[2 * wall_idx], walls_cells_data[2 * wall_idx + 1]));
                mstate.SetWallNodes(wall_idx,
                        make_pair(walls_nodes_data[2 * wall_idx], walls_nodes_data[2 * wall_idx + 1]));
        }

        // Clean up temp storage
        delete[] walls_id_data;
        delete[] walls_cells_data;
        delete[] walls_nodes_data;

        // Find datasets that are called 'walls_attr_<name>' (I.e.: wall attributes)
        list < string > attr_names = find_matching_links(loc_id, [](string name)->bool {
                return boost::starts_with(name, "walls_attr_");
        });

        // Try to read all matching wall attribute datasets
        for (string long_name : attr_names) {
                // Magic number 11 is the length of the prefix "cells_attr_" that needs top be stripped away
                string short_name = long_name.substr(11, string::npos);

                // Check if "long_name" refers to a valid rank-1 dataset in loc_id
                H5O_info_t obj_info;
                H5Oget_info_by_name(loc_id, long_name.c_str(), &obj_info, H5P_DEFAULT);
                if (H5O_TYPE_DATASET == obj_info.type) {
                        // Open the dataset to get all info
                        hid_t dset = H5Dopen(loc_id, long_name.c_str(), H5P_DEFAULT);

                        // Make sure dimensions match
                        hid_t dspace = H5Dget_space(dset);
                        // Check if rank = 1
                        const int ndims = H5Sget_simple_extent_ndims(dspace);
                        if (1 != ndims) {
                                H5Sclose(dspace);
                                break;
                        }
                        // Check if number of attribute values = number of walls
                        hsize_t dims = 0;
                        H5Sget_simple_extent_dims(dspace, &dims, 0);
                        if (num_walls != dims) {
                                H5Sclose(dspace);
                                break;
                        }

                        // Get the type of the attribute from HDF5 (either int or double)
                        hid_t dtype = H5Dget_type(dset);
                        H5T_class_t dtype_class = H5Tget_class(dtype);
                        H5Tclose(dtype);
                        if (H5T_INTEGER == dtype_class) {
                                // Read the data directly into a vector
                                vector<int> attr_values(num_walls);
                                H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, attr_values.data());
                                // Add attribute and its values to the walls
                                mstate.WallAttributeContainer().Add<int>(short_name);
                                mstate.WallAttributeContainer().SetAll<int>(short_name, attr_values);
                        } else if (H5T_FLOAT == dtype_class) {
                                // Read the data directly into a vector
                                vector<double> attr_values(num_walls);
                                H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, attr_values.data());
                                // Add attribute and its values to the walls
                                mstate.WallAttributeContainer().Add<double>(short_name);
                                mstate.WallAttributeContainer().SetAll<double>(short_name, attr_values);
                        } else if (H5T_STRING == dtype_class) {
                                // Read the raw string data directly into a vector of C-style strings
                                vector<char*> raw_attr_values(num_walls);
                                hid_t ntype = H5Tcopy(H5T_C_S1);
                                H5Tset_size(ntype, H5T_VARIABLE);
                                H5Dread(dset, ntype, H5S_ALL, H5S_ALL, H5P_DEFAULT, raw_attr_values.data());
                                H5Tclose(ntype);

                                // Convert the C-style strings to a C++ string vector
                                vector < string > attr_values(raw_attr_values.begin(), raw_attr_values.end());

                                // Add attribute and its values to the walls
                                mstate.WallAttributeContainer().Add<string>(short_name);
                                mstate.WallAttributeContainer().SetAll<string>(short_name, attr_values);
                        } else {
                                // Unsupported datatype: ignore (or generate a warning)
                        }

                        // Clean up
                        H5Dclose(dset);
                }
        }
}

ptree Hdf5File::ReadParameters(hid_t loc_id)
{
        hid_t param_grp_id = H5Gopen(loc_id, "parameters", H5P_DEFAULT);
        if (param_grp_id < 0) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': The group \'/step_*/parameters\' does not exist. File might be corrupted.");
                throw runtime_error(errmsg);
        }
        ptree params = hdf52ptree(param_grp_id);
        H5Gclose(param_grp_id);
        return params;
}

ptree Hdf5File::ReadPtreeFromStringAttribute(hid_t loc_id, const string& attr_name)
{
        // Empty ptree to be filled with the contents of attr_name
        ptree pt;

        // Try to open the requested attribute at loc_id
        hid_t attr_id = H5Aopen(loc_id, attr_name.c_str(), H5P_DEFAULT);

        // If attribute was not found, an empty ptree will be returned
        if (attr_id > 0) {
                hid_t attr_dtype = H5Aget_type(attr_id);
                hid_t attr_dspace = H5Aget_space(attr_id);

                // Check if 'attr_name' is a string attribute indeed
                if (H5T_STRING == H5Tget_class(attr_dtype)) {
                        // Note that H5T_STRING attributes are scalars so one does not check for ndims and dims.
                        // Just read the attribute as a C string.
                        char* attr_value;
                        hid_t attr_ntype = H5Tcopy(H5T_C_S1);
                        H5Tset_size(attr_ntype, H5T_VARIABLE);
                        H5Aread(attr_id, attr_ntype, &attr_value);
                        H5Tclose(attr_ntype);

                        // Convert the attribute value from a C string into a stream to be translated to a ptree.
                        stringstream pt_strm;
                        pt_strm << attr_value;
                        read_xml(pt_strm, pt, trim_whitespace);

                        // Cleanup
                        H5Dvlen_reclaim(attr_ntype, attr_dspace, H5P_DEFAULT, attr_value);
                }

                // Clean up handles and signal success to the calling function
                H5Sclose(attr_dspace);
                H5Tclose(attr_dtype);
                H5Aclose(attr_id);
        }

        return pt;
}

vector<pair<int, double>> Hdf5File::TimeSteps()
{
        if (!IsOpen()) {
                throw runtime_error("[Hdf5File]: File needs to be open prior to reading.");
        }

        // =========================================================================
        // Read information about time steps contained in the file
        double* time_steps_data = nullptr;
        size_t time_steps_dim1 = 0;
        read_array_rank1(m_file_id, "time_steps", &time_steps_data, H5T_NATIVE_DOUBLE, &time_steps_dim1);

        int* time_steps_idx_data = nullptr;
        size_t time_steps_idx_dim1 = 0;
        read_array_rank1(m_file_id, "time_steps_idx", &time_steps_idx_data, H5T_NATIVE_INT, &time_steps_idx_dim1);

        // Check size consistency of /time_steps and /time_step_idx
        if (time_steps_dim1 != time_steps_idx_dim1) {
                string errmsg = string("Error reading file \'") + m_file_path
                        + string("\': The /time_steps and /time_steps_idx datasets must be the same size. File might be corrupted.");
                // Clean up temp storage
                delete[] time_steps_data;
                delete[] time_steps_idx_data;
                throw runtime_error(errmsg);
        }

        vector<pair<int, double>> time_steps(time_steps_dim1);
        for (size_t i = 0; i < time_steps_dim1; ++i) {
                time_steps[i] = make_pair(time_steps_idx_data[i], time_steps_data[i]);
        }

        // Clean up temp storage
        delete[] time_steps_data;
        delete[] time_steps_idx_data;

        return time_steps;
}

void Hdf5File::Write(shared_ptr<const SimPT_Sim::Sim> sim)
{
        // Just a wrapper around Write() to comply with the interface of FileViewer.
        // (could be replaced with WriteState such that Sim doesn't have to be passed around)
        WriteState(sim->GetState());
}

void Hdf5File::WriteCells(hid_t loc_id, const MeshState& mstate)
{
        const size_t num_cells = mstate.GetNumCells();
        // Cell information is stored in a variety of arrays and matrices that
        // keep geometric information as well as other cell attributes.
        // All arrays have length 'num_cells', unless explicitly specified otherwise.
        // All matrices are stored as row-major C-style arrays.
        //
        // Each cell has a unique ID, these are stored in the array 'cells_id'
        //
        // Geometric information of a cell consists of:
        // - Number of nodes used to define a cell: stored in the 'cells_num_nodes' array.
        // - IDs of the nodes that define the cell's polygon: 'cells_nodes' matrix
        //   Each row in this matrix represents a cell. Elements of each row are
        //   IDs of the nodes (cfr.: 'nodes_id') that make up the cell's polygon.
        //   Since not all cells have the same number of nodes, some storage space is lost.
        //   Although cell's polygons are well-defined by both 'cells_nodes' and
        //   'cells_num_nodes' combined, padding (-1) is used to make the
        //   'cells_nodes' array easier to read and interpret.
        //
        // Attributes of cells might be for instance:
        // - Values of cell-based chemicals.
        // - Area and target area of each cell.
        // - Type of each cell.
        // - And many others ...

        // Cell IDs
        const vector<int> cells_id = mstate.GetCellsID();
        write_array_rank1(loc_id, "cells_id", cells_id.data(), H5T_STD_I32LE, num_cells);

        // Max number of nodes in a cell
        size_t max_cells_num_nodes = 0;
        const vector<size_t> cells_num_nodes = mstate.GetCellsNumNodes();
        for (size_t cell_num_nodes : cells_num_nodes) {
                if (cell_num_nodes > max_cells_num_nodes) {
                        max_cells_num_nodes = cell_num_nodes;
                }
        }

        // Max number of walls in a cell
        size_t max_cells_num_walls = 0;
        const vector<size_t> cells_num_walls = mstate.GetCellsNumWalls();
        for (size_t cell_num_walls : cells_num_walls) {
                if (cell_num_walls > max_cells_num_walls) {
                        max_cells_num_walls = cell_num_walls;
                }
        }

        // Allocate space for geometric data and copy that from the state
        int* cells_num_nodes_data = new int[num_cells];
        int* cells_nodes_data = new int[num_cells * max_cells_num_nodes];
        int* cells_num_walls_data = new int[num_cells];
        int* cells_walls_data = new int[num_cells * max_cells_num_walls];

        for (size_t cell_idx = 0; cell_idx < num_cells; ++cell_idx) {
                // Number of nodes in the cell
                cells_num_nodes_data[cell_idx] = cells_num_nodes[cell_idx];
                // Nodes of the cells (i.e.: polygonal geometry information)
                size_t cell_node_idx = 0;
                for (const int cell_node_id : mstate.GetCellNodes(cell_idx)) {
                        cells_nodes_data[(cell_idx * max_cells_num_nodes) + cell_node_idx] = cell_node_id;
                        ++cell_node_idx;
                }
                // Pad the remaining positions in the row with -1
                while (cell_node_idx < max_cells_num_nodes) {
                        cells_nodes_data[(cell_idx * max_cells_num_nodes) + cell_node_idx] = -1;
                        ++cell_node_idx;
                }

                // Number of walls in the cell
                cells_num_walls_data[cell_idx] = cells_num_walls[cell_idx];
                // Walls of the cells (i.e.: connectivity information)
                size_t cell_wall_idx = 0;
                for (const int cell_wall_id : mstate.GetCellWalls(cell_idx)) {
                        cells_walls_data[(cell_idx * max_cells_num_walls) + cell_wall_idx] = cell_wall_id;
                        ++cell_wall_idx;
                }
                // Pad the remaining positions in the row with -1
                while (cell_wall_idx < max_cells_num_walls) {
                        cells_walls_data[(cell_idx * max_cells_num_walls) + cell_wall_idx] = -1;
                        ++cell_wall_idx;
                }
        }

        // Write data arrays to HDF5
        // (Don't write cells_walls in the special case of one cell, which has no walls)
        write_array_rank1(loc_id, "cells_num_nodes", cells_num_nodes_data, H5T_STD_I32LE, num_cells);
        write_array_rank2(loc_id, "cells_nodes", cells_nodes_data, H5T_STD_I32LE, num_cells, max_cells_num_nodes);
        write_array_rank1(loc_id, "cells_num_walls", cells_num_walls_data, H5T_STD_I32LE, num_cells);
        if (1 != num_cells) {
                write_array_rank2(loc_id, "cells_walls", cells_walls_data, H5T_STD_I32LE, num_cells, max_cells_num_walls);
        }

        // Clean up temporary arrays
        delete[] cells_num_nodes_data;
        delete[] cells_nodes_data;
        delete[] cells_num_walls_data;
        delete[] cells_walls_data;

        // Cell attributes (int)
        for (string attr_name : mstate.CellAttributeContainer().GetNames<int>()) {
                const vector<int> attr_values = mstate.CellAttributeContainer().GetAll<int>(attr_name);
                string attr_name_hdf5 = string("cells_attr_") + attr_name;
                write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values.data(), H5T_STD_I32LE, num_cells);
        }

        // Cell attributes (double)
        for (string attr_name : mstate.CellAttributeContainer().GetNames<double>()) {
                const vector<double> attr_values = mstate.CellAttributeContainer().GetAll<double>(attr_name);
                string attr_name_hdf5 = string("cells_attr_") + attr_name;
                write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values.data(), H5T_IEEE_F64LE, num_cells);
        }

        // Cell attributes (string)
        for (string attr_name : mstate.CellAttributeContainer().GetNames<string>()) {
                const vector<string> attr_values = mstate.CellAttributeContainer().GetAll<string>(attr_name);
                string attr_name_hdf5 = string("cells_attr_") + attr_name;
                write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values);
        }
}

void Hdf5File::WriteBPNodes(hid_t loc_id, const MeshState& mstate)
{
        vector<int> bp_nodes = mstate.GetBoundaryPolygonNodes();
        write_array_rank1(loc_id, "bp_nodes", bp_nodes.data(), H5T_STD_I32LE, bp_nodes.size());
}

void Hdf5File::WriteBPWalls(hid_t loc_id, const MeshState& mstate)
{
        vector<int> bp_walls = mstate.GetBoundaryPolygonWalls();
        write_array_rank1(loc_id, "bp_walls", bp_walls.data(), H5T_STD_I32LE, bp_walls.size());
}

void Hdf5File::WriteMesh(hid_t loc_id, const MeshState& mstate)
{
        string const here = string(VL_HERE) + "> ";
        // Split up in three parts for readability
        WriteNodes(loc_id, mstate);
        WriteCells(loc_id, mstate);
        WriteBPNodes(loc_id, mstate);
        // In the very special case of just one cell (e.g. in the beginning of
        // a simulation, there's no walls. Hence: don't write walls in that case)
        const size_t num_cells = mstate.GetNumCells();
        const size_t num_walls = mstate.GetNumWalls();
        if (num_cells > 1 && num_walls > 0) {
                WriteWalls(loc_id, mstate);
                WriteBPWalls(loc_id, mstate);
        } else if ((num_cells == 1 && num_walls != 0) || (num_cells > 1 && num_walls == 0)) {
                cerr << here
                        << "[Hdf5File] Error: inconsistent number of cells and walls: (#cells > 1) iff (#walls > 0)"
                        << endl;
        }
        // The case (#cells == 1) && (#walls == 0) means that no walls should be written, which is handled silently
}

void Hdf5File::WriteNodes(hid_t loc_id, const MeshState& mstate)
{
        const size_t num_nodes = mstate.GetNumNodes();
        // TODO: Document how nodes are written to HDF5

        // Node IDs
        vector<int> nodes_id = mstate.GetNodesID();
        write_array_rank1(loc_id, "nodes_id", nodes_id.data(), H5T_STD_I32LE, num_nodes);

        // Node xy coordinates
        double* nodes_xy_data = new double[2 * num_nodes];
        for (size_t node_idx = 0; node_idx < num_nodes; ++node_idx) {
                nodes_xy_data[2 * node_idx] = mstate.GetNodeX(node_idx);
                nodes_xy_data[2 * node_idx + 1] = mstate.GetNodeY(node_idx);
        }
        write_array_rank2(loc_id, "nodes_xy", nodes_xy_data, H5T_IEEE_F64LE, num_nodes, 2);
        delete[] nodes_xy_data;

        // Node attributes (int)
        for (string attr_name : mstate.NodeAttributeContainer().GetNames<int>()) {
                vector<int> attr_values = mstate.NodeAttributeContainer().GetAll<int>(attr_name);
                string attr_name_hdf5 = string("nodes_attr_") + attr_name;
                write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values.data(), H5T_STD_I32LE, num_nodes);
        }

        // Node attributes (double)
        for (string attr_name : mstate.NodeAttributeContainer().GetNames<double>()) {
                vector<double> attr_values = mstate.NodeAttributeContainer().GetAll<double>(attr_name);
                string attr_name_hdf5 = string("nodes_attr_") + attr_name;
                write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values.data(), H5T_IEEE_F64LE, num_nodes);
        }

        // Node attributes (string)
        for (string attr_name : mstate.NodeAttributeContainer().GetNames<string>()) {
                vector < string > attr_values = mstate.NodeAttributeContainer().GetAll<string>(attr_name);
                string attr_name_hdf5 = string("nodes_attr_") + attr_name;
                write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values);
        }
}

void Hdf5File::WriteParameters(hid_t loc_id, const ptree& params)
{
        // Create a HDF5 group /step_n/parameters where the subtrees aof param
        // will be stored as subgroups
        hid_t param_grp = H5Gcreate(loc_id, "parameters", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        // Transform the ptree structure to HDF5 subgroups and named attributes
        // with parameter values inside the group [loc_id]/parameters
        ptree2hdf5(params, param_grp);
        // Clean up HDF5 handles
        H5Gclose(param_grp);
}

void Hdf5File::WritePtreeToStringAttribute(hid_t loc_id, const string& attr_name, const ptree& pt)
{
        // Get the ptree in a C style string
        stringstream pt_strm;
        write_xml(pt_strm, pt, XmlWriterSettings::GetTab());
        string pt_str = pt_strm.str();
        const char* pt_cstr = pt_str.c_str();

        // Write pt_str to HDF5 attribute attr_name as a variable-length string
        hid_t attr_dspace = H5Screate(H5S_SCALAR);
        hid_t attr_dtype = H5Tcopy(H5T_C_S1);
        H5Tset_size(attr_dtype, H5T_VARIABLE);
        hid_t attr_id = H5Acreate(loc_id, attr_name.c_str(), attr_dtype, attr_dspace, H5P_DEFAULT, H5P_DEFAULT);
        H5Awrite(attr_id, attr_dtype, &pt_cstr);
        H5Aclose(attr_id);
        H5Sclose(attr_dspace);
        H5Tclose(attr_dtype);
}

void Hdf5File::WriteState(const SimState& sstate)
{
        const string here = string(VL_HERE) + "> ";

        // ========================================================================
        // First do the following checks
        // - Can '/time_steps' be opened & are the dimensions correct?
        //  -> fatal if condition not met
        // - Can '/time_steps_idx' be opened & are the dimensions correct?
        //  -> fatal if condition not met
        // - Is sstate.GetTimeStep() already present in 'time_steps_idx' OR are there any time steps beyond?
        //  -> remove all time steps(_idx) beyond that point, then add new time step
        // - Is the group '/step_n' already present for the current step?
        //  -> remove all /step_n groups starting from that point, then add new time step group
        // ========================================================================

        double* time_steps_data = nullptr;
        size_t time_steps_dim1 = 0;
        read_array_rank1(m_file_id, "time_steps", &time_steps_data, H5T_NATIVE_DOUBLE, &time_steps_dim1);

        int* time_steps_idx_data = nullptr;
        size_t time_steps_idx_dim1 = 0;
        read_array_rank1(m_file_id, "time_steps_idx", &time_steps_idx_data, H5T_NATIVE_INT, &time_steps_idx_dim1);

        // Sanity check: dimensions of 'time_steps' and 'time_steps_idx' must match
        if (time_steps_dim1 != time_steps_idx_dim1) {
                // Clean up temporary structures for time_steps
                delete[] time_steps_data;
                delete[] time_steps_idx_data;
                throw runtime_error("[Hdf5File] Error: dim of time_steps and time_steps_idx do not match.");
        }

        // Shorthand notations
        size_t& num_time_steps = time_steps_dim1;
        const double sim_time = sstate.GetTime();
        const int sim_step = sstate.GetTimeStep();
        string step_grp_name = string("/step_") + to_string(sim_step);

        //cout << "HDF5: trying to write " << step_grp_name << endl;

        // Safe cases where no truncation is done are:
        //   num_time_steps = 0
        //   sim_step > time_steps_idx[end]
        // Hence truncation must be performed in the following case:
        if (num_time_steps != 0 && sim_step <= time_steps_idx_data[num_time_steps - 1]) {
                // cout << "[Hdf5File] File will be truncated." << endl;
                // To find the last step in time_steps_idx BEFORE sim_step loop over
                // the time_steps_idx values backwards until sim_step is reached
                size_t i = num_time_steps - 1;
                // cout << "[Hdf5File] Removing /step_[";
                while (static_cast<int>(i) >= 0 && sim_step <= time_steps_idx_data[i]) {
                        // cout << time_steps_idx_data[i] << ' ';
                        string name = string("step_") + to_string(time_steps_idx_data[i]);
                        // Read: http://www.hdfgroup.org/HDF5/doc/H5.user/Performance.html#Freespace
                        // to see why your data files are getting bigger while being truncated.
                        H5Ldelete(m_file_id, name.c_str(), H5P_DEFAULT);
                        --i;
                }
                //cout << ']' << endl;
                // Now i is the index of the last time_steps_idx BEFORE sim_step
                // cout << "HDF5: i = " << i << " time_steps_idx[i] = " << time_steps_idx_data[i] << endl;
                // Remove elements with index i+1 -> num_time_steps-1 from time_steps(_idx)
                // by resizing both datasets to contain i + 1 elements.
                {
                        hid_t dset = H5Dopen(m_file_id, "time_steps", H5P_DEFAULT);
                        hsize_t dims;
                        hid_t dspace = H5Dget_space(dset);
                        H5Sget_simple_extent_dims(dspace, &dims, 0);
                        dims = i + 1;
                        H5Dset_extent(dset, &dims);
                        H5Sclose(dspace);
                        H5Dclose(dset);
                }
                {
                        hid_t dset = H5Dopen(m_file_id, "time_steps_idx", H5P_DEFAULT);
                        hsize_t dims;
                        hid_t dspace = H5Dget_space(dset);
                        H5Sget_simple_extent_dims(dspace, &dims, 0);
                        dims = i + 1;
                        H5Dset_extent(dset, &dims);
                        H5Sclose(dspace);
                        H5Dclose(dset);
                }
        }

        // Clean up temporary structures for time_steps
        delete[] time_steps_data;
        delete[] time_steps_idx_data;

        // ========================================================================

        // Proceed with saving the current sim state
        // (HDF5 checks are limited to the minimum since they are already caught above)
        // Extending /time_steps with current step =============================
        hid_t time_steps_dset = H5Dopen(m_file_id, "time_steps", H5P_DEFAULT);

        // Get the current size of /time_steps, increase the size by one
        // and extend the /time_steps dataset to accomodate the new time step value
        hsize_t time_steps_dims;
        hid_t time_steps_dspace = H5Dget_space(time_steps_dset);
        H5Sget_simple_extent_dims(time_steps_dspace, &time_steps_dims, 0);
        time_steps_dims++;
        H5Dset_extent(time_steps_dset, &time_steps_dims);

        // To append the current time value to /time_steps properly:
        //  1. create the destination dataspace to append the current time step value
        //     to the newly extended /time_steps
        //     (i.e. select the last point of /time_steps)
        //          H5S_SELECT_SET overwrites any previous selection
        //          1 is the number of elements that are selected
        //          last_time_step_coord is the index of the last /time_steps element after expansion
        //  2. create the source dataspace for the current time step value
        //     (i.e. just a simple dataspace curr_time_step_dspace for one element)
        //  3. Write the current time value with the previously defined source/destination spaces

        // 1.
        hsize_t last_time_step_coord = time_steps_dims - 1;
        H5Sselect_elements(time_steps_dspace, H5S_SELECT_SET, 1, &last_time_step_coord);
        // 2.
        hsize_t curr_time_step_dims[] = { 1 };
        hid_t curr_time_step_dspace = H5Screate_simple(1, curr_time_step_dims, NULL);
        // 3.
        H5Dwrite(time_steps_dset, H5T_IEEE_F64LE, curr_time_step_dspace, time_steps_dspace, H5P_DEFAULT, &sim_time);
        // Clean up HDF5 handles
        H5Sclose(curr_time_step_dspace);
        H5Sclose(time_steps_dspace);
        H5Dclose(time_steps_dset);
        // End of extending /time_steps ========================================

        // Extending /time_steps_idx with current step index ===================
        hid_t time_steps_idx_dset = H5Dopen(m_file_id, "time_steps_idx", H5P_DEFAULT);

        // Get the current size of /time_steps_idx, increase the size by one
        // and extend the /time_steps_idx data set to accommodate the new time step value
        hsize_t time_steps_idx_dims;
        hid_t time_steps_idx_dspace = H5Dget_space(time_steps_idx_dset);
        H5Sget_simple_extent_dims(time_steps_idx_dspace, &time_steps_idx_dims, 0);
        time_steps_idx_dims++;
        H5Dset_extent(time_steps_idx_dset, &time_steps_idx_dims);

        // To append the current step index to /time_steps_idx properly:
        //  1. create the destination data space to append the current time step index
        //     to the newly extended /time_steps_idx
        //     (i.e. select the last point of /time_steps_idx)
        //          H5S_SELECT_SET overwrites any previous selection
        //          1 is the number of elements that are selected
        //          last_time_step_coord is the index of the last /time_steps_idx element after expansion
        //  2. create the source data space for the current step index
        //     (i.e. just a simple data space curr_time_step_idx_dspace for one element)
        //  3. Write the current step index value with the previously defined source/destination spaces

        // 1.
        hsize_t last_time_step_idx_coord = time_steps_idx_dims - 1;
        H5Sselect_elements(time_steps_idx_dspace, H5S_SELECT_SET, 1, &last_time_step_idx_coord);
        // 2.
        hsize_t curr_time_step_idx_dims[] = { 1 };
        hid_t curr_time_step_idx_dspace = H5Screate_simple(1, curr_time_step_idx_dims, NULL);
        // 3.
        H5Dwrite(time_steps_idx_dset, H5T_STD_I32LE, curr_time_step_idx_dspace, time_steps_idx_dspace, H5P_DEFAULT,
                &sim_step);
        // Clean up HDF5 handles
        H5Sclose(curr_time_step_idx_dspace);
        H5Sclose(time_steps_idx_dspace);
        H5Dclose(time_steps_idx_dset);
        // End of extending /time_steps_idx ====================================

        // Creating & writing the '/step_n' group ==============================
        hid_t step_grp = H5Gcreate(m_file_id, step_grp_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        // Save the mesh: Nodes, Cells & Walls to '/step_n' group
        WriteMesh(step_grp, sstate.GetMeshState());
        // Get and write the ptree of current simulation parameters
        // TODO: Naively write the whole parameters ptree as a serialized XML version of the ptree
        // TODO: to an attribute of the current '/step_n' group called 'parameters'
        // TODO: This is a quick fix for something that should be resolved by rethinking
        // TODO: the format for saving data in XML which I believe is flawed atm.
        WritePtreeToStringAttribute(step_grp, "parameters", sstate.GetParameters());
        // Get and write the ptree with the current random engine state
        WritePtreeToStringAttribute(step_grp, "random_engine", sstate.GetRandomEngineState());
        // Close the '/step_n' group
        H5Gclose(step_grp);
        // End of creating & writing the '/step_n' group =======================
}

void Hdf5File::WriteWalls(hid_t loc_id, const MeshState& mstate)
{
    const size_t num_walls = mstate.GetNumWalls();
    // TODO: Document how walls are written to HDF5
    // - 'walls_id'
    //   IDs of the walls. These are being used by cells to refer to their walls
    // - 'walls_cells' array keeps the 'neighbor cells' information
    //   ie, which cells are connected by a specific wall
    //   rows are walls, the two cols are the cells that are connected
    // - 'walls_nodes' array keeps the wall end-nodes information
    //   rows are walls, the two cols are the nodes n1 and n2, the endpoints of the wall

    // Wall IDs
    const vector<int> walls_id = mstate.GetWallsID();
    write_array_rank1(loc_id, "walls_id", walls_id.data(), H5T_STD_I32LE, num_walls);

    // Copy wall properties into temporary arrays used for saving to HDF5
    int* walls_cells_data = new int[2 * num_walls];
    int* walls_nodes_data = new int[2 * num_walls];

    for (size_t wall_idx = 0; wall_idx < num_walls; ++wall_idx) {
        // IDs of the two cells connected by the wall
        pair<int, int> cells_ids = mstate.GetWallCells(wall_idx);
        walls_cells_data[2 * wall_idx]     = cells_ids.first;
        walls_cells_data[2 * wall_idx + 1] = cells_ids.second;
        // IDs of the two endpoints of the wall
        pair<int, int> nodes_ids = mstate.GetWallNodes(wall_idx);
        walls_nodes_data[2 * wall_idx]     = nodes_ids.first;
        walls_nodes_data[2 * wall_idx + 1] = nodes_ids.second;
    }

    // Write data arrays to HDF5
    write_array_rank2(loc_id, "walls_cells", walls_cells_data, H5T_STD_I32LE, num_walls, 2);
    write_array_rank2(loc_id, "walls_nodes", walls_nodes_data, H5T_STD_I32LE, num_walls, 2);

    // Clean up temporary arrays
    delete[] walls_cells_data;
    delete[] walls_nodes_data;

    // Wall attributes (int)
    for (string attr_name : mstate.WallAttributeContainer().GetNames<int>()) {
        const vector<int> attr_values = mstate.WallAttributeContainer().GetAll<int>(attr_name);
        string attr_name_hdf5 = string("walls_attr_") + attr_name;
        write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values.data(), H5T_STD_I32LE, num_walls);
    }

    // Wall attributes (double)
    for (string attr_name : mstate.WallAttributeContainer().GetNames<double>()) {
        const vector<double> attr_values = mstate.WallAttributeContainer().GetAll<double>(attr_name);
        string attr_name_hdf5 = string("walls_attr_") + attr_name;
        write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values.data(), H5T_IEEE_F64LE, num_walls);
    }

    // Wall attributes (string)
    for (string attr_name : mstate.WallAttributeContainer().GetNames<string>()) {
        const vector<string> attr_values = mstate.WallAttributeContainer().GetAll<string>(attr_name);
        string attr_name_hdf5 = string("walls_attr_") + attr_name;
        write_array_rank1(loc_id, attr_name_hdf5.c_str(), attr_values);
    }
}

} // namespace