Hdf5DataReader.cpp

/*

Copyright (c) 2005-2012, University of Oxford.
All rights reserved.

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.
 * Neither the name of the University of Oxford nor the names of its
   contributors may be used to endorse or promote products derived from this
   software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

#include "Hdf5DataReader.hpp"
#include "Exception.hpp"
#include "OutputFileHandler.hpp"
#include "PetscTools.hpp"

#include <cassert>
#include <algorithm>

Hdf5DataReader::Hdf5DataReader(const std::string& rDirectory,
                               const std::string& rBaseName,
                               bool makeAbsolute)
    : mBaseName(rBaseName),
      mIsUnlimitedDimensionSet(false),
      mNumberTimesteps(1),
      mIsDataComplete(true),
      mClosed(false)
{
    RelativeTo::Value relative_to;
    if (makeAbsolute)
    {
        relative_to = RelativeTo::ChasteTestOutput;
    }
    else
    {
        relative_to = RelativeTo::Absolute;
    }
    FileFinder directory(rDirectory, relative_to);
    CommonConstructor(directory, rBaseName);
}

Hdf5DataReader::Hdf5DataReader(const FileFinder& rDirectory,
                               const std::string& rBaseName)
    : mBaseName(rBaseName),
      mIsUnlimitedDimensionSet(false),
      mNumberTimesteps(1),
      mIsDataComplete(true),
      mClosed(false)
{
    CommonConstructor(rDirectory, rBaseName);
}

void Hdf5DataReader::CommonConstructor(const FileFinder& rDirectory, const std::string& rBaseName)
{
    std::string results_dir = rDirectory.GetAbsolutePath();
    if (!rDirectory.IsDir() || !rDirectory.Exists())
    {
        EXCEPTION("Directory does not exist: " + results_dir);
    }
    mDirectory = results_dir;
    assert(*(mDirectory.end()-1) == '/'); // paranoia

    std::string file_name = results_dir + mBaseName + ".h5";
    FileFinder h5_file(file_name,RelativeTo::Absolute);

    if (!h5_file.Exists())
    {
        EXCEPTION("Hdf5DataReader could not open " + file_name + " , as it does not exist.");
    }

    // Open the file and the main dataset
    mFileId = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);

    if (mFileId <= 0)
    {
        EXCEPTION("Hdf5DataReader could not open " << file_name <<
                  " , H5Fopen error code = " << mFileId);
    }

    mVariablesDatasetId = H5Dopen(mFileId, "Data");

    hid_t variables_dataspace = H5Dget_space(mVariablesDatasetId);
    mVariablesDatasetRank = H5Sget_simple_extent_ndims(variables_dataspace);

    // Get the dataset/dataspace dimensions
    hsize_t dataset_max_sizes[MAX_DATASET_RANK];
    H5Sget_simple_extent_dims(variables_dataspace, mVariablesDatasetSizes, dataset_max_sizes);

    for (unsigned i=1; i<MAX_DATASET_RANK; i++)  // Zero is excluded since it may be unlimited
    {
        assert(mVariablesDatasetSizes[i] == dataset_max_sizes[i]);
    }

    // Check if an unlimited dimension has been defined
    if (dataset_max_sizes[0] == H5S_UNLIMITED)
    {
        mIsUnlimitedDimensionSet = true;
        mTimeDatasetId = H5Dopen(mFileId, "Time");

        hid_t timestep_dataspace = H5Dget_space(mTimeDatasetId);

        // Get the dataset/dataspace dimensions
        H5Sget_simple_extent_dims(timestep_dataspace, &mNumberTimesteps, NULL);

    }

    // Get the attribute where the name of the variables are stored
    hid_t attribute_id = H5Aopen_name(mVariablesDatasetId, "Variable Details");

    // Get attribute datatype, dataspace, rank, and dimensions
    hid_t attribute_type  = H5Aget_type(attribute_id);
    hid_t attribute_space = H5Aget_space(attribute_id);

    hsize_t attr_dataspace_dim;
    H5Sget_simple_extent_dims(attribute_space, &attr_dataspace_dim, NULL);

    unsigned num_columns = H5Sget_simple_extent_npoints(attribute_space);
    char* string_array = (char *)malloc(sizeof(char)*MAX_STRING_SIZE*(int)num_columns);
    H5Aread(attribute_id, attribute_type, string_array);

    // Loop over column names and store them.
    for (unsigned index=0; index < num_columns; index++)
    {
        // Read the string from the raw vector
        std::string column_name_unit(&string_array[MAX_STRING_SIZE*index]);

        // Find beginning of unit definition.
        size_t name_length = column_name_unit.find('(');
        size_t unit_length = column_name_unit.find(')') - name_length - 1;

        std::string column_name = column_name_unit.substr(0, name_length);
        std::string column_unit = column_name_unit.substr(name_length+1, unit_length);

        mVariableNames.push_back(column_name);
        mVariableToColumnIndex[column_name] = index;
        mVariableToUnit[column_name] = column_unit;
    }

    // Release all the identifiers
    H5Tclose(attribute_type);
    H5Sclose(attribute_space);
    H5Aclose(attribute_id);

    // Free allocated memory
    free(string_array);

    // Find out if it's incomplete data
    attribute_id = H5Aopen_name(mVariablesDatasetId, "IsDataComplete");
    if (attribute_id < 0)
    {
        // This is in the old format (before we added the IsDataComplete attribute).
        // Just quit (leaving a nasty hdf5 error).
        return;
    }

    attribute_type  = H5Aget_type(attribute_id);
    attribute_space = H5Aget_space(attribute_id);
    unsigned is_data_complete;
    H5Aread(attribute_id, H5T_NATIVE_UINT, &is_data_complete);

    // Release all the identifiers
    H5Tclose(attribute_type);
    H5Sclose(attribute_space);
    H5Aclose(attribute_id);
    mIsDataComplete = (is_data_complete == 1) ? true : false;

    if (is_data_complete)
    {
        return;
    }

    // Incomplete data
    // Read the vector thing
    attribute_id = H5Aopen_name(mVariablesDatasetId, "NodeMap");
    attribute_type  = H5Aget_type(attribute_id);
    attribute_space = H5Aget_space(attribute_id);

    // Get the dataset/dataspace dimensions
    unsigned num_node_indices = H5Sget_simple_extent_npoints(attribute_space);

    // Read data from hyperslab in the file into the hyperslab in memory
    mIncompleteNodeIndices.clear();
    mIncompleteNodeIndices.resize(num_node_indices);
    H5Aread(attribute_id, H5T_NATIVE_UINT, &mIncompleteNodeIndices[0]);

    H5Tclose(attribute_type);
    H5Sclose(attribute_space);
    H5Aclose(attribute_id);
}

std::vector<double> Hdf5DataReader::GetVariableOverTime(const std::string& rVariableName,
                                                        unsigned nodeIndex)
{
    if (!mIsUnlimitedDimensionSet)
    {
        EXCEPTION("The file does not contain time dependent data");
    }

    unsigned actual_node_index = nodeIndex;
    if (!mIsDataComplete)
    {
        unsigned node_index = 0;
        for (node_index=0; node_index<mIncompleteNodeIndices.size(); node_index++)
        {
            if (mIncompleteNodeIndices[node_index]==nodeIndex)
            {
                actual_node_index = node_index;
                break;
            }
        }
        if ( node_index == mIncompleteNodeIndices.size())
        {
            EXCEPTION("The incomplete file does not contain info of node " << nodeIndex);
        }
    }
    if (actual_node_index >= mVariablesDatasetSizes[1])
    {
        EXCEPTION("The file doesn't contain info of node " << actual_node_index);
    }

    std::map<std::string, unsigned>::iterator col_iter = mVariableToColumnIndex.find(rVariableName);
    if (col_iter == mVariableToColumnIndex.end())
    {
        EXCEPTION("The file doesn't contain data for variable " + rVariableName);
    }
    int column_index = (*col_iter).second;

    // Define hyperslab in the dataset.
    hsize_t offset[3] = {0, actual_node_index, column_index};
    hsize_t count[3]  = {mVariablesDatasetSizes[0], 1, 1};
    hid_t variables_dataspace = H5Dget_space(mVariablesDatasetId);
    H5Sselect_hyperslab(variables_dataspace, H5S_SELECT_SET, offset, NULL, count, NULL);

    // Define a simple memory dataspace
    hid_t memspace = H5Screate_simple(1, &mVariablesDatasetSizes[0] ,NULL);

    // Data buffer to return
    std::vector<double> ret(mVariablesDatasetSizes[0]);

    // Read data from hyperslab in the file into the hyperslab in memory
    H5Dread(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, variables_dataspace, H5P_DEFAULT, &ret[0]);

    H5Sclose(variables_dataspace);
    H5Sclose(memspace);

    return ret;
}

std::vector<std::vector<double> > Hdf5DataReader::GetVariableOverTimeOverMultipleNodes(const std::string& rVariableName,
                                                                                       unsigned lowerIndex,
                                                                                       unsigned upperIndex)
{
    if (!mIsUnlimitedDimensionSet)
    {
        EXCEPTION("The file does not contain time dependent data");
    }

    if (!mIsDataComplete)
    {
        EXCEPTION("GetVariableOverTimeOverMultipleNodes() cannot be called using incomplete data sets (those for which data was only written for certain nodes)");
    }

    if (upperIndex > mVariablesDatasetSizes[1])
    {
       EXCEPTION("The file doesn't contain info for node " << upperIndex-1);
    }

    std::map<std::string, unsigned>::iterator col_iter = mVariableToColumnIndex.find(rVariableName);
    if (col_iter == mVariableToColumnIndex.end())
    {
        EXCEPTION("The file doesn't contain data for variable " + rVariableName);
    }
    int column_index = (*col_iter).second;

    // Define hyperslab in the dataset.
    hsize_t offset[3] = {0, lowerIndex, column_index};
    hsize_t count[3]  = {mVariablesDatasetSizes[0], upperIndex-lowerIndex, 1};
    hid_t variables_dataspace = H5Dget_space(mVariablesDatasetId);
    H5Sselect_hyperslab(variables_dataspace, H5S_SELECT_SET, offset, NULL, count, NULL);

    // Define a simple memory dataspace
    hsize_t data_dimensions[2];
    data_dimensions[0] = mVariablesDatasetSizes[0];
    data_dimensions[1] = upperIndex-lowerIndex;
    hid_t memspace = H5Screate_simple(2, data_dimensions, NULL);

    double* data_read = new double[mVariablesDatasetSizes[0]*(upperIndex-lowerIndex)];

    // Read data from hyperslab in the file into the hyperslab in memory
    H5Dread(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, variables_dataspace, H5P_DEFAULT, data_read);

    H5Sclose(variables_dataspace);
    H5Sclose(memspace);

    // Data buffer to return
    unsigned num_nodes_read = upperIndex-lowerIndex;
    unsigned num_timesteps = mVariablesDatasetSizes[0];

    std::vector<std::vector<double> > ret(num_nodes_read);

    for (unsigned node_num=0; node_num<num_nodes_read; node_num++)
    {
        ret[node_num].resize(num_timesteps);
        for (unsigned time_num=0; time_num<num_timesteps; time_num++)
        {
            ret[node_num][time_num] = data_read[num_nodes_read*time_num + node_num];
        }
    }

    delete[] data_read;

    return ret;
}

void Hdf5DataReader::GetVariableOverNodes(Vec data,
                                          const std::string& rVariableName,
                                          unsigned timestep)
{
    if (!mIsDataComplete)
    {
        EXCEPTION("You can only get a vector for complete data");
    }
    if (!mIsUnlimitedDimensionSet && timestep!=0)
    {
        EXCEPTION("The file does not contain time dependent data");
    }

    std::map<std::string, unsigned>::iterator col_iter = mVariableToColumnIndex.find(rVariableName);
    if (col_iter == mVariableToColumnIndex.end())
    {
        EXCEPTION("The file does not contain data for variable " + rVariableName);
    }
    int column_index = (*col_iter).second;

    // Check for valid timestep
    if (timestep >= mNumberTimesteps)
    {
        EXCEPTION("The file does not contain data for timestep number " << timestep);
    }

    int lo, hi, size;
    VecGetSize(data, &size);
    if ((unsigned)size != mVariablesDatasetSizes[1])
    {
        EXCEPTION("Could not read data because Vec is the wrong size");
    }
    // Get range owned by each processor
    VecGetOwnershipRange(data, &lo, &hi);

    if (hi > lo) // i.e. we own some...
    {
        // Define a dataset in memory for this process
        hsize_t v_size[1] = {hi-lo};
        hid_t memspace = H5Screate_simple(1, v_size, NULL);

        // Select hyperslab in the file.
        hsize_t offset[3] = {timestep, lo, column_index};
        hsize_t count[3]  = {1, hi-lo, 1};
        hid_t hyperslab_space = H5Dget_space(mVariablesDatasetId);
        H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, offset, NULL, count, NULL);

        double* p_petsc_vector;
        VecGetArray(data, &p_petsc_vector);
        herr_t err;
        err=H5Dread(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, H5P_DEFAULT, p_petsc_vector);
        assert(err==0);
        VecRestoreArray(data, &p_petsc_vector);

        H5Sclose(hyperslab_space);
        H5Sclose(memspace);
    }
}

std::vector<double> Hdf5DataReader::GetUnlimitedDimensionValues()
{
    // Data buffer to return
    std::vector<double> ret(mNumberTimesteps);

    if (!mIsUnlimitedDimensionSet)
    {
        // Fake it
        assert(mNumberTimesteps==1);
        ret[0] = 0.0;
        return ret;
    }
    // Define hyperslab in the dataset
    hid_t time_dataspace = H5Dget_space(mTimeDatasetId);

    // Define a simple memory dataspace
    hid_t memspace = H5Screate_simple(1, &mNumberTimesteps ,NULL);

    // Read data from hyperslab in the file into the hyperslab in memory
    H5Dread(mTimeDatasetId, H5T_NATIVE_DOUBLE, memspace, time_dataspace, H5P_DEFAULT, &ret[0]);

    H5Sclose(time_dataspace);
    H5Sclose(memspace);

    return ret;
}

void Hdf5DataReader::Close()
{
    if (!mClosed)
    {
        H5Dclose(mVariablesDatasetId);
        if (mIsUnlimitedDimensionSet)
        {
            H5Dclose(mTimeDatasetId);
        }
        H5Fclose(mFileId);
        mClosed = true;
    }
}

Hdf5DataReader::~Hdf5DataReader()
{
    Close();
}

unsigned Hdf5DataReader::GetNumberOfRows()
{
    return mVariablesDatasetSizes[1];
}

std::vector<std::string> Hdf5DataReader::GetVariableNames()
{
    return mVariableNames;
}

std::string Hdf5DataReader::GetUnit(const std::string& rVariableName)
{
    return mVariableToUnit[rVariableName];
}

bool Hdf5DataReader::IsDataComplete()
{
    return mIsDataComplete;
}

std::vector<unsigned> Hdf5DataReader::GetIncompleteNodeMap()
{
    return mIncompleteNodeIndices;
}