Hdf5DataReader.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 #include "Hdf5DataReader.hpp"
00029 #include "Exception.hpp"
00030 #include "OutputFileHandler.hpp"
00031 
00032 #include <cassert>
00033 
00034 
00035 Hdf5DataReader::Hdf5DataReader(const std::string& rDirectory,
00036                                const std::string& rBaseName,
00037                                bool makeAbsolute)
00038     : mBaseName(rBaseName),
00039       mIsUnlimitedDimensionSet(false),
00040       mNumberTimesteps(1),
00041       mIsDataComplete(true)
00042 {
00043     RelativeTo::Value relative_to;
00044     if (makeAbsolute)
00045     {
00046         relative_to = RelativeTo::ChasteTestOutput;
00047     }
00048     else
00049     {
00050         relative_to = RelativeTo::Absolute;
00051     }
00052     FileFinder directory(rDirectory, relative_to);
00053     CommonConstructor(directory, rBaseName);
00054 }
00055 
00056 Hdf5DataReader::Hdf5DataReader(const FileFinder& rDirectory,
00057                                const std::string& rBaseName)
00058     : mBaseName(rBaseName),
00059       mIsUnlimitedDimensionSet(false),
00060       mNumberTimesteps(1),
00061       mIsDataComplete(true)
00062 {
00063     CommonConstructor(rDirectory, rBaseName);
00064 }
00065 
00066 void Hdf5DataReader::CommonConstructor(const FileFinder& rDirectory, const std::string& rBaseName)
00067 {
00068     std::string results_dir = rDirectory.GetAbsolutePath();
00069     if (!rDirectory.IsDir() || !rDirectory.Exists())
00070     {
00071         EXCEPTION("Directory does not exist: " + results_dir);
00072     }
00073     mDirectory = results_dir;
00074     assert(*(mDirectory.end()-1) == '/'); // paranoia
00075 
00076     std::string file_name = results_dir + mBaseName + ".h5";
00077 
00078     // Open the file and the main dataset
00079     mFileId = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
00080 
00081     if (mFileId <= 0)
00082     {
00083         EXCEPTION("Hdf5DataReader could not open " + file_name);
00084     }
00085     mVariablesDatasetId = H5Dopen(mFileId, "Data");
00086 
00087     hid_t variables_dataspace = H5Dget_space(mVariablesDatasetId);
00088     mVariablesDatasetRank = H5Sget_simple_extent_ndims(variables_dataspace);
00089 
00090     // Get the dataset/dataspace dimensions
00091     hsize_t dataset_max_sizes[MAX_DATASET_RANK];
00092     H5Sget_simple_extent_dims(variables_dataspace, mVariablesDatasetSizes, dataset_max_sizes);
00093 
00094     for (unsigned i=1; i<MAX_DATASET_RANK; i++)  // Zero is excluded since it may be unlimited
00095     {
00096         assert(mVariablesDatasetSizes[i] == dataset_max_sizes[i]);
00097     }
00098 
00099     // Check if an unlimited dimension has been defined
00100     if (dataset_max_sizes[0] == H5S_UNLIMITED)
00101     {
00102         mIsUnlimitedDimensionSet = true;
00103         mTimeDatasetId = H5Dopen(mFileId, "Time");
00104 
00105         hid_t timestep_dataspace = H5Dget_space(mTimeDatasetId);
00106 
00107         // Get the dataset/dataspace dimensions
00108         H5Sget_simple_extent_dims(timestep_dataspace, &mNumberTimesteps, NULL);
00109 
00110     }
00111 
00112     // Get the attribute where the name of the variables are stored
00113     hid_t attribute_id = H5Aopen_name(mVariablesDatasetId, "Variable Details");
00114 
00115     // Get attribute datatype, dataspace, rank, and dimensions
00116     hid_t attribute_type  = H5Aget_type(attribute_id);
00117     hid_t attribute_space = H5Aget_space(attribute_id);
00118 
00119     hsize_t attr_dataspace_dim;
00120     H5Sget_simple_extent_dims(attribute_space, &attr_dataspace_dim, NULL);
00121 
00122     unsigned num_columns = H5Sget_simple_extent_npoints(attribute_space);
00123     char* string_array = (char *)malloc(sizeof(char)*MAX_STRING_SIZE*(int)num_columns);
00124     H5Aread(attribute_id, attribute_type, string_array);
00125 
00126     // Loop over column names and store them.
00127     for (unsigned index=0; index < num_columns; index++)
00128     {
00129         // Read the string from the raw vector
00130         std::string column_name_unit(&string_array[MAX_STRING_SIZE*index]);
00131 
00132         // Find beginning of unit definition.
00133         size_t name_length = column_name_unit.find('(');
00134         size_t unit_length = column_name_unit.find(')') - name_length - 1;
00135 
00136         std::string column_name = column_name_unit.substr(0, name_length);
00137         std::string column_unit = column_name_unit.substr(name_length+1, unit_length);
00138 
00139         mVariableNames.push_back(column_name);
00140         mVariableToColumnIndex[column_name] = index;
00141         mVariableToUnit[column_name] = column_unit;
00142     }
00143 
00144     // Release all the identifiers
00145     H5Tclose(attribute_type);
00146     H5Sclose(attribute_space);
00147     H5Aclose(attribute_id);
00148 
00149     // Free allocated memory
00150     free(string_array);
00151 
00152     // Find out if it's incomplete data
00153     attribute_id = H5Aopen_name(mVariablesDatasetId, "IsDataComplete");
00154     if (attribute_id < 0)
00155     {
00156         // This is in the old format (before we added the IsDataComplete attribute).
00157         // Just quit (leaving a nasty hdf5 error).
00158         return;
00159     }
00160 
00161     attribute_type  = H5Aget_type(attribute_id);
00162     attribute_space = H5Aget_space(attribute_id);
00163     unsigned is_data_complete;
00164     H5Aread(attribute_id, H5T_NATIVE_UINT, &is_data_complete);
00165 
00166     // Release all the identifiers
00167     H5Tclose(attribute_type);
00168     H5Sclose(attribute_space);
00169     H5Aclose(attribute_id);
00170     mIsDataComplete = (is_data_complete == 1) ? true : false;
00171 
00172     if (is_data_complete)
00173     {
00174         return;
00175     }
00176 
00177     // Incomplete data
00178     // Read the vector thing
00179     attribute_id = H5Aopen_name(mVariablesDatasetId, "NodeMap");
00180     attribute_type  = H5Aget_type(attribute_id);
00181     attribute_space = H5Aget_space(attribute_id);
00182 
00183     // Get the dataset/dataspace dimensions
00184     unsigned num_node_indices = H5Sget_simple_extent_npoints(attribute_space);
00185 
00186     // Read data from hyperslab in the file into the hyperslab in memory
00187     mIncompleteNodeIndices.clear();
00188     mIncompleteNodeIndices.resize(num_node_indices);
00189     H5Aread(attribute_id, H5T_NATIVE_UINT, &mIncompleteNodeIndices[0]);
00190 
00191     H5Tclose(attribute_type);
00192     H5Sclose(attribute_space);
00193     H5Aclose(attribute_id);
00194 }
00195 
00196 std::vector<double> Hdf5DataReader::GetVariableOverTime(const std::string& rVariableName,
00197                                                         unsigned nodeIndex)
00198 {
00199     if (!mIsUnlimitedDimensionSet)
00200     {
00201         EXCEPTION("The file does not contain time dependent data");
00202     }
00203 
00204     unsigned actual_node_index = nodeIndex;
00205     if (!mIsDataComplete)
00206     {
00207         unsigned node_index = 0;
00208         for (node_index=0; node_index<mIncompleteNodeIndices.size(); node_index++)
00209         {
00210             if (mIncompleteNodeIndices[node_index]==nodeIndex)
00211             {
00212                 actual_node_index = node_index;
00213                 break;
00214             }
00215         }
00216         if ( node_index == mIncompleteNodeIndices.size())
00217         {
00218             std::stringstream ss;
00219             ss << "The incomplete file does not contain info of node " << nodeIndex ;
00220             EXCEPTION(ss.str());
00221         }
00222     }
00223     if (actual_node_index >= mVariablesDatasetSizes[1])
00224     {
00225         std::stringstream ss;
00226         ss << "The file doesn't contain info of node " << actual_node_index ;
00227         EXCEPTION(ss.str());
00228     }
00229 
00230     std::map<std::string, unsigned>::iterator col_iter = mVariableToColumnIndex.find(rVariableName);
00231     if (col_iter == mVariableToColumnIndex.end())
00232     {
00233         EXCEPTION("The file doesn't contain data for variable " + rVariableName);
00234     }
00235     int column_index = (*col_iter).second;
00236 
00237     // Define hyperslab in the dataset.
00238     hsize_t offset[3] = {0, actual_node_index, column_index};
00239     hsize_t count[3]  = {mVariablesDatasetSizes[0], 1, 1};
00240     hid_t variables_dataspace = H5Dget_space(mVariablesDatasetId);
00241     H5Sselect_hyperslab(variables_dataspace, H5S_SELECT_SET, offset, NULL, count, NULL);
00242 
00243     // Define a simple memory dataspace
00244     hid_t memspace = H5Screate_simple(1, &mVariablesDatasetSizes[0] ,NULL);
00245 
00246     // Data buffer to return
00247     std::vector<double> ret(mVariablesDatasetSizes[0]);
00248 
00249     // Read data from hyperslab in the file into the hyperslab in memory
00250     H5Dread(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, variables_dataspace, H5P_DEFAULT, &ret[0]);
00251 
00252     H5Sclose(variables_dataspace);
00253     H5Sclose(memspace);
00254 
00255     return ret;
00256 }
00257 
00258 void Hdf5DataReader::GetVariableOverNodes(Vec data,
00259                                           const std::string& rVariableName,
00260                                           unsigned timestep)
00261 {
00262     if (!mIsDataComplete)
00263     {
00264         EXCEPTION("You can only get a vector for complete data");
00265     }
00266     if (!mIsUnlimitedDimensionSet && timestep!=0)
00267     {
00268         EXCEPTION("The file does not contain time dependent data");
00269     }
00270 
00271     std::map<std::string, unsigned>::iterator col_iter = mVariableToColumnIndex.find(rVariableName);
00272     if (col_iter == mVariableToColumnIndex.end())
00273     {
00274         EXCEPTION("The file does not contain data for variable " + rVariableName);
00275     }
00276     int column_index = (*col_iter).second;
00277 
00278     // Check for valid timestep
00279     if (timestep >= mNumberTimesteps)
00280     {
00281         std::stringstream ss;
00282         ss << "The file does not contain data for timestep number " << timestep;
00283         EXCEPTION(ss.str());
00284     }
00285 
00286     int lo, hi, size;
00287     VecGetSize(data, &size);
00288     if ((unsigned)size != mVariablesDatasetSizes[1])
00289     {
00290         EXCEPTION("Could not read data because Vec is the wrong size");
00291     }
00292     // Get range owned by each processor
00293     VecGetOwnershipRange(data, &lo, &hi);
00294 
00295     if (hi > lo) // i.e. we own some...
00296     {
00297         // Define a dataset in memory for this process
00298         hsize_t v_size[1] = {hi-lo};
00299         hid_t memspace = H5Screate_simple(1, v_size, NULL);
00300 
00301         // Select hyperslab in the file.
00302         hsize_t offset[3] = {timestep, lo, column_index};
00303         hsize_t count[3]  = {1, hi-lo, 1};
00304         hid_t hyperslab_space = H5Dget_space(mVariablesDatasetId);
00305         H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, offset, NULL, count, NULL);
00306 
00307         double* p_petsc_vector;
00308         VecGetArray(data, &p_petsc_vector);
00309         herr_t err;
00310         err=H5Dread(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, H5P_DEFAULT, p_petsc_vector);
00311         assert(err==0);
00312         VecRestoreArray(data, &p_petsc_vector);
00313 
00314         H5Sclose(hyperslab_space);
00315         H5Sclose(memspace);
00316     }
00317 }
00318 
00319 std::vector<double> Hdf5DataReader::GetUnlimitedDimensionValues()
00320 {
00321     // Data buffer to return
00322     std::vector<double> ret(mNumberTimesteps);
00323 
00324     if (!mIsUnlimitedDimensionSet)
00325     {
00326         // Fake it
00327         assert(mNumberTimesteps==1);
00328         ret[0] = 0.0;
00329         return ret;
00330     }
00331     // Define hyperslab in the dataset
00332     hid_t time_dataspace = H5Dget_space(mTimeDatasetId);
00333 
00334     // Define a simple memory dataspace
00335     hid_t memspace = H5Screate_simple(1, &mNumberTimesteps ,NULL);
00336 
00337     // Read data from hyperslab in the file into the hyperslab in memory
00338     H5Dread(mTimeDatasetId, H5T_NATIVE_DOUBLE, memspace, time_dataspace, H5P_DEFAULT, &ret[0]);
00339 
00340     H5Sclose(time_dataspace);
00341     H5Sclose(memspace);
00342 
00343     return ret;
00344 }
00345 
00346 void Hdf5DataReader::Close()
00347 {
00349     H5Dclose(mVariablesDatasetId);
00350 
00351     if (mIsUnlimitedDimensionSet)
00352     {
00353         H5Dclose(mTimeDatasetId);
00354     }
00355 
00356     H5Fclose(mFileId);
00357 }
00358 
00359 Hdf5DataReader::~Hdf5DataReader()
00360 {
00361     Close();
00362 }
00363 
00364 unsigned Hdf5DataReader::GetNumberOfRows()
00365 {
00366     return mVariablesDatasetSizes[1];
00367 }
00368 
00369 std::vector<std::string> Hdf5DataReader::GetVariableNames()
00370 {
00371     return mVariableNames;
00372 }
00373 
00374 std::string Hdf5DataReader::GetUnit(const std::string& rVariableName)
00375 {
00376     return mVariableToUnit[rVariableName];
00377 }
00378 
00379 bool Hdf5DataReader::IsDataComplete()
00380 {
00381     return mIsDataComplete;
00382 }
00383 
00384 std::vector<unsigned> Hdf5DataReader::GetIncompleteNodeMap()
00385 {
00386     return mIncompleteNodeIndices;
00387 }

Generated on Tue May 31 14:31:42 2011 for Chaste by  doxygen 1.5.5