gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
hdf5_combine_tools.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
19 
20 #ifndef __hdf5_combine_tools_hpp__
21 #define __hdf5_combine_tools_hpp__
22 
23 #include <vector>
24 #include <sstream>
25 #include <unordered_set>
26 #include <unordered_map>
27 #include <hdf5.h>
28 
34 
35 namespace Gambit
36 {
37  namespace Printers
38  {
39  namespace HDF5
40  {
41  template <typename T>
42  inline T type_ret(){return T();}
43 
44  template <class U, typename... T>
45  void Enter_HDF5(hid_t dataset, T&... params);
46 
47  struct read_hdf5
48  {
49  template <typename U, typename T>
50  static void run(U, hid_t &dataset, std::vector <T> &vec)
51  {
52  hid_t dataspace = H5Dget_space(dataset);
53  hssize_t dim_t = H5Sget_simple_extent_npoints(dataspace);
54  std::vector<U> data(dim_t);
55  H5Dread( dataset, get_hdf5_data_type<U>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&data[0]);
56  vec = std::vector<T>(data.begin(), data.end());
57  H5Sclose(dataspace);
58  }
59  };
60 
61  struct copy_hdf5
62  {
63  template <typename U>
64  static void run(U, hid_t &dataset_out, std::vector<hid_t> &datasets, unsigned long long &size_tot, std::vector<unsigned long long> &sizes, hid_t &old_dataset, size_t offset)
65  {
66  std::vector<U> data(size_tot);
67  unsigned long long j = 0;
68 
69  if (old_dataset >= 0)
70  {
71  hid_t space = H5Dget_space(old_dataset);
72  hsize_t dim_t = H5Sget_simple_extent_npoints(space);
73  H5Sclose(space);
74  data.resize(dim_t + size_tot);
75  H5Dread(old_dataset, get_hdf5_data_type<U>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&data[0]);
76  j = dim_t;
77  //std::cout <<"Read "<<dim_t<<" points from old dataset (into position 0)!"<<std::endl;
78  //for(unsigned long long i=0; i<dim_t; i++) // debug
79  //{
80  // std::cout<<" "<<data[i]<<std::endl;
81  //}
82  }
83  else
84  {
85  //std::cout << "No old dataset found!" << std::endl;
86  }
87 
88  for (int i = 0, end = datasets.size(); i < end; i++)
89  {
90  hsize_t dim_t;
91  if(datasets[i] >= 0)
92  {
93  // Check that we allocated enough space for this dataset!
94  hid_t space = H5Dget_space(datasets[i]);
95  dim_t = H5Sget_simple_extent_npoints(space);
96  H5Sclose(space);
97  if(j+dim_t > data.size())
98  {
99  if( (data.size() - j) < sizes[i]) // was !=, but I think it is fine so long as there isn't too *little* space. TODO: might be an odds measuring bug though since I think these should indeed be == if everything has happened correctly.
100  {
101  // This is how much space we expect to have left in the buffer
102  // at the end. If not, something bad has happened!
103  std::ostringstream errmsg;
104  errmsg << "Error copying parameter. Dataset in input file " << i
105  << " was larger than previously measured and doesn't fit in the allocated buffer! "
106  << "This is a bug in the combination tools, please report it." <<std::endl;
107  errmsg << "(j="<<j<<", dim_t="<<dim_t<<", data.size()="<<data.size()<<", sizes["<<i<<"]="<<sizes[i]<<")"<<std::endl;
108  printer_error().raise(LOCAL_INFO, errmsg.str());
109  }
110  // Just some junk buffer points at the end of the last data set.
111  // Need to do a more careful dataset selection to chop these out so
112  // that the data fits in our buffer (and in the space already allocated
113  // in the output file)
114  hid_t memspace_id, dspace_id;
115  std::pair<hid_t,hid_t> chunk_ids = HDF5::selectChunk(datasets[i], 0, sizes[i]);
116  memspace_id = chunk_ids.first;
117  dspace_id = chunk_ids.second;
118 
119  H5Dread(datasets[i], get_hdf5_data_type<U>::type(), memspace_id, dspace_id, H5P_DEFAULT, (void *)&data[j]);
120  H5Sclose(memspace_id);
121  H5Sclose(dspace_id);
122  //std::cout <<"Read "<<sizes[i]<<" points from dataset "<<i<<" into position "<<j<<std::endl;
123  //std::ostringstream errmsg;
124  //errmsg << "Error copying parameter. Dataset in input file " << i
125  // << " was larger than previously measured and doesn't fit in the allocated buffer! "
126  // << "This is a bug in the combination tools, please report it." <<std::endl;
127  //errmsg << "(j="<<j<<", dim_t="<<dim_t<<", data.size()="<<data.size()<<")"<<std::endl;
128  //printer_error().raise(LOCAL_INFO, errmsg.str());
129  }
130  else
131  {
132  // Read the whole dataset in to buffer (faster than selecting, I should think)
133  H5Dread(datasets[i], get_hdf5_data_type<U>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&data[j]);
134  //std::cout <<"Read "<<dim_t<<" points from dataset "<<i<<" into position "<<j<<std::endl;
135  //for(unsigned long long k=j; k<j+dim_t; k++) // debug
136  //{
137  // std::cout<<" "<<data[k]<<std::endl;
138  //}
139  }
140  }
141  else
142  {
143  // dataset didn't exist for this file, skip it.
144  dim_t = 0;
145  }
146 
147  // Check size consistency
148  if (dim_t >= sizes[i])
149  {
150  // Data had expected size, no problem
151  // NOTE: Why can dim_t be larger than expected? Possible padding at end or something?
152  j += sizes[i];
153  }
154  else if(dim_t==0)
155  {
156  // Data was missing, but also probably fine, just skip it
157  j += sizes[i];
158  }
159  else
160  {
161  // Data has some other unexpected size, error!
162  std::ostringstream errmsg;
163  errmsg << "Error copying parameter "". Dataset in input file " << i << " did not have the expected size" <<std::endl;
164  errmsg << "(sizes["<<i<<"] = "<<sizes[i]<<" was less than dim_t = "<<dim_t<<")";
165  printer_error().raise(LOCAL_INFO, errmsg.str());
166  }
167  }
168  // Select dataspace for writing in the output file. Partial write now possible, so we
169  // have to select the hyperslab that we want to write into
170  hid_t memspace_id, dspace_id;
171  std::pair<hid_t,hid_t> chunk_ids = HDF5::selectChunk(dataset_out, offset, data.size());
172  memspace_id = chunk_ids.first;
173  dspace_id = chunk_ids.second;
174  herr_t err = H5Dwrite( dataset_out, get_hdf5_data_type<U>::type(), memspace_id, dspace_id, H5P_DEFAULT, (void *)&data[0]);
175  if(err<0)
176  {
177  std::ostringstream errmsg;
178  errmsg << "Error copying parameter. HD5write failed." <<std::endl;
179  printer_error().raise(LOCAL_INFO, errmsg.str());
180  }
181  //std::cout<<"Wrote "<<data.size()<<" points into output dataset at position "<<offset<<std::endl;
182  H5Sclose(memspace_id);
183  H5Sclose(dspace_id);
184 
185  // Was:
186  //H5Dwrite( dataset_out, get_hdf5_data_type<U>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&data[0]);
187  }
188  };
189 
191  {
192  template <typename U>
193  static void run (U, hid_t &dataset_out, hid_t &dataset2_out, std::vector<hid_t> &datasets, std::vector<hid_t> &datasets2, const unsigned long long size, const std::unordered_map<PPIDpair, unsigned long long, PPIDHash, PPIDEqual>& RA_write_hash, const std::vector<std::vector <unsigned long long> > &pointid, const std::vector<std::vector <unsigned long long> > &rank, const std::vector<unsigned long long> &aux_sizes, hid_t &/*old_dataset*/, hid_t &/*old_dataset2*/)
194  {
195  std::vector<U> output(size, 0);
196  std::vector<int> valids(size, 0);
197 
198  // Should no longer need the old datasets, they should have already been copied during "copy_hdf5"
199  // if (old_dataset >= 0 && old_dataset2 >= 0)
200  // {
201  // hid_t space = HDF5::getSpace(old_dataset);
202  // hid_t space2 = HDF5::getSpace(old_dataset2);
203  // hsize_t old_size = H5Sget_simple_extent_npoints(space);
204  // hsize_t old_size2 = H5Sget_simple_extent_npoints(space2);
205  // HDF5::closeSpace(space);
206  // HDF5::closeSpace(space2);
207  // if(old_size > size or old_size2 > size)
208  // {
209  // std::ostringstream errmsg;
210  // errmsg << "Error copying old dataset into buffer! The old dataset has a larger size than has been allocated for new data! (old_size="<<old_size<<", old_size2="<<old_size2<<", new_size="<<size<<")";
211  // printer_error().raise(LOCAL_INFO, errmsg.str());
212  // }
213  // H5Dread(old_dataset, get_hdf5_data_type<U>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&output[0]);
214  // H5Dread(old_dataset2, get_hdf5_data_type<int>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&valids[0]);
215  // }
216 
217  // Instead we just need to read in the recently-copied primary points and replace
218  // entries as needed (TODO: would use much less memory if we just write individual replacements
219  // straight to disk)
220  if (dataset_out >= 0 && dataset2_out >= 0)
221  {
222  hid_t space = HDF5::getSpace(dataset_out);
223  hid_t space2 = HDF5::getSpace(dataset2_out);
224  hsize_t out_size = HDF5::getSimpleExtentNpoints(space);
225  hsize_t out_size2 = HDF5::getSimpleExtentNpoints(space2);
226  HDF5::closeSpace(space);
227  HDF5::closeSpace(space2);
228  if(out_size > size or out_size2 > size)
229  {
230  std::ostringstream errmsg;
231  errmsg << "Error copying dataset into buffer for RA replacements! The dataset has a larger size than has been allocated for new data! (out_size="<<out_size<<", out_size2="<<out_size2<<", expected_size="<<size<<")";
232  printer_error().raise(LOCAL_INFO, errmsg.str());
233  }
234  H5Dread(dataset_out, get_hdf5_data_type<U>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&output[0]);
235  H5Dread(dataset2_out, get_hdf5_data_type<int>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&valids[0]);
236  }
237  else
238  {
239  std::ostringstream errmsg;
240  errmsg << "Error copying datasets into buffer for RA replacements! Could not open datasets (were they created properly during 'copy_hdf5' operation?).";
241  printer_error().raise(LOCAL_INFO, errmsg.str());
242  }
243 
244  // Check that all the input given is consistent in length
245  size_t ndsets = datasets.size();
246  #define DSET_SIZE_CHECK(VEC) \
247  if(VEC.size()!=ndsets) \
248  { \
249  std::ostringstream errmsg; \
250  errmsg << STRINGIFY(VEC) << " vector has inconsistent size! ("<<VEC.size()<<", should be "<<ndsets<<"). This is a bug, please report it."; \
251  printer_error().raise(LOCAL_INFO, errmsg.str()); \
252  }
253  DSET_SIZE_CHECK(datasets2)
254  DSET_SIZE_CHECK(aux_sizes)
255  DSET_SIZE_CHECK(pointid)
257  #undef DSET_SIZE_CHECK
258 
259  // Additional iterators to be iterated in sync with dataset iteration
260  // We are actually only copying over one 'parameter' in this function,
261  // but are looping over matching datasets from the temp files from
262  // each rank. So each iteration goes to the next file, but remains
263  // targeting the same parameter.
264  auto st = aux_sizes.begin();
265  auto pt = pointid.begin();
266  auto ra = rank.begin();
267  auto itv = datasets2.begin();
268  for (auto it = datasets.begin(); it != datasets.end();
269  ++it, ++itv, ++pt, ++ra, ++st)
270  {
271  if(*it < 0)
272  {
273  // Dataset wasn't opened, probably some RA parameter just happened to not
274  // exist in a certain temporary file. Skip this dataset. Though check that
275  // Dataset2 also wasn't opened.
276  if(*itv >= 0)
277  {
278  std::ostringstream errmsg;
279  errmsg << "dataset2 iterator ('isvalid') points to an open dataset, while dataset iterator (main dataset) does not. This is inconsistent and indicates either a bug in this combine code, or in the code which generated the datasets, please report it.";
280  printer_error().raise(LOCAL_INFO, errmsg.str());
281  }
282  }
283  else
284  {
285  hid_t space = H5Dget_space(*it);
286  hssize_t dim_t = H5Sget_simple_extent_npoints(space);
287  std::vector<U> data(dim_t);
288  std::vector<bool> valid;
289  Enter_HDF5<read_hdf5> (*itv, valid);
290  H5Dread(*it, get_hdf5_data_type<U>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&data[0]);
291 
292  if((unsigned long long)dim_t < *st)
293  {
294  std::ostringstream errmsg;
295  errmsg << "Error copying aux parameter. Input file smaller than required.";
296  printer_error().raise(LOCAL_INFO, errmsg.str());
297  }
298 
299  for (int i = 0, end = *st; i < end; i++)
300  {
301  if (valid[i])
302  {
303  // Look up target for write in hash map
304  std::unordered_map<PPIDpair, unsigned long long, PPIDHash, PPIDEqual>::const_iterator ihash = RA_write_hash.find(PPIDpair((*pt)[i],(*ra)[i]));
305  if(ihash != RA_write_hash.end())
306  {
307  // found hash key, copy data
308  unsigned long long temp = ihash->second;
309  if(temp > size)
310  {
311  std::ostringstream errmsg;
312  errmsg << "Error copying random access parameter. The hash entry for "
313  << "pt number " << (*pt)[i] << " of rank " << (*ra)[i]
314  << " targets the point outside the size of the output dataset ("<<temp<<" > "<<size<<")."
315  << "This indicates"
316  << " a bug in the hash generation, please report it.";
317  printer_error().raise(LOCAL_INFO, errmsg.str());
318  }
319  output[temp] = data[i];
320  valids[temp] = 1;
321  }
322  else
323  {
324  std::ostringstream errmsg;
325  errmsg << "Error copying random access parameter. Could not find "
326  << "pt number " << (*pt)[i] << " of rank " << (*ra)[i]
327  << " in the output dataset (hash entry was not found).";
328  printer_error().raise(LOCAL_INFO, errmsg.str());
329  }
330  }
331  }
332  } // end if
333  }
334 
335  H5Dwrite( dataset_out, get_hdf5_data_type<U>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&output[0]);
336  H5Dwrite( dataset2_out, get_hdf5_data_type<int>::type(), H5S_ALL, H5S_ALL, H5P_DEFAULT, (void *)&valids[0]);
337  }
338  };
339 
340  template <class U, typename... T>
341  inline void Enter_HDF5(hid_t dataset, T&... params)
342  {
343  if(dataset<0)
344  {
345  std::ostringstream errmsg;
346  errmsg << "Invalid dataset supplied to Enter_HDF5 routine!";
347  printer_error().raise(LOCAL_INFO, errmsg.str());
348  }
349 
350  hid_t dtype = H5Dget_type(dataset);
351  if(dtype<0)
352  {
353  std::ostringstream errmsg;
354  errmsg << "Failed to detect type for dataset provides as argument for Enter_HDF5 routine!";
355  printer_error().raise(LOCAL_INFO, errmsg.str());
356  }
357 
358  //H5T_class_t cl = H5Tget_class(dtype);
359  hid_t type = H5Tget_native_type(dtype, H5T_DIR_DESCEND);
360  if(type<0)
361  {
362  std::ostringstream errmsg;
363  errmsg << "Failed to detect native type for dataset provides as argument for Enter_HDF5 routine!";
364  printer_error().raise(LOCAL_INFO, errmsg.str());
365  }
366 
367  if (H5Tequal(type, get_hdf5_data_type<float>::type()))
368  {
369  U::run(float(), dataset, params...);
370  }
371  else if (H5Tequal(type, get_hdf5_data_type<double>::type()))
372  {
373  U::run(double(), dataset, params...);
374  }
375  else if (H5Tequal(type, get_hdf5_data_type<long double>::type()))
376  {
377  U::run(type_ret<long double>(), dataset, params...);
378  }
379  else if (H5Tequal(type, get_hdf5_data_type<char>::type()))
380  {
381  U::run(char(), dataset, params...);
382  }
383  else if (H5Tequal(type, get_hdf5_data_type<short>::type()))
384  {
385  U::run(short(), dataset, params...);
386  }
387  else if (H5Tequal(type, get_hdf5_data_type<int>::type()))
388  {
389  U::run(int(), dataset, params...);
390  }
391  else if (H5Tequal(type, get_hdf5_data_type<long>::type()))
392  {
393  U::run(long(), dataset, params...);
394  }
395  else if (H5Tequal(type, get_hdf5_data_type<long long>::type()))
396  {
397  U::run(type_ret<long long>(), dataset, params...);
398  }
399  else if (H5Tequal(type, get_hdf5_data_type<unsigned char>::type()))
400  {
401  U::run(type_ret<unsigned char>(), dataset, params...);
402  }
403  else if (H5Tequal(type, get_hdf5_data_type<unsigned short>::type()))
404  {
405  U::run(type_ret<unsigned short>(), dataset, params...);
406  }
407  else if (H5Tequal(type, get_hdf5_data_type<unsigned int>::type()))
408  {
409  U::run(type_ret<unsigned int>(), dataset, params...);
410  }
411  else if (H5Tequal(type, get_hdf5_data_type<unsigned long>::type()))
412  {
413  U::run(type_ret<unsigned long>(), dataset, params...);
414  }
415  else if (H5Tequal(type, get_hdf5_data_type<unsigned long long>::type()))
416  {
417  U::run(type_ret<unsigned long long>(), dataset, params...);
418  }
419  else
420  {
421  std::ostringstream errmsg;
422  errmsg << "Could not deduce input hdf5 parameter type.";
423  printer_error().raise(LOCAL_INFO, errmsg.str());
424  }
425 
426  H5Tclose(dtype);
427  }
428 
430  {
431  private:
432  std::string group_name;
433  std::vector<std::string> param_names, aux_param_names;
434  std::unordered_set<std::string> param_set, aux_param_set; // for easier finding
435  std::vector<unsigned long long> cum_sizes;
436  std::vector<unsigned long long> sizes;
437  unsigned long long size_tot;
438  unsigned long long size_tot_l;
439  std::string root_file_name;
440  bool do_cleanup; // If true, delete old temporary files after successful combination
441  std::string get_fname(const size_t); // get name of ith temp file
442  bool skip_unreadable; // If true, ignores temp files that fail to open
443  std::vector<hid_t> files;
444  std::vector<hid_t> groups;
445  std::vector<hid_t> aux_groups;
446  std::vector<std::string> file_names; // Names of temp files to combine
447  bool custom_mode; // Running in mode which allows 'custom' filenames (for combining output from multiple runs)
448 
449  public:
450  hdf5_stuff(const std::string &base_file_name, const std::string &output_file, const std::string &group_name, const size_t num, const bool cleanup, const bool skip, const std::vector<std::string>& input_files);
451  ~hdf5_stuff(); // close files on destruction
452  void Enter_Aux_Parameters(const std::string &output_file, bool resume = false);
453  };
454 
455  inline void combine_hdf5_files(const std::string output_file, const std::string &base_file_name, const std::string &group, const size_t num, const bool resume, const bool cleanup, const bool skip, const std::vector<std::string> input_files = std::vector<std::string>())
456  {
457  hdf5_stuff stuff(base_file_name, output_file, group, num, cleanup, skip, input_files);
458 
459  stuff.Enter_Aux_Parameters(output_file, resume);
460  }
461 
462  // Helper function to compute target point hash for RA combination
463  std::unordered_map<PPIDpair, unsigned long long, PPIDHash, PPIDEqual> get_RA_write_hash(hid_t, std::unordered_set<PPIDpair,PPIDHash,PPIDEqual>&);
464 
466  std::pair<std::vector<std::string>,std::vector<size_t>> find_temporary_files(const std::string& finalfile);
467  std::pair<std::vector<std::string>,std::vector<size_t>> find_temporary_files(const std::string& finalfile, size_t& max_i);
468 
469  }
470  }
471 }
472 
473 #endif
474 
greatScanData data
Definition: great.cpp:38
dataset H5Sget_simple_extent_npoints
Definition: hdf5tools.cpp:720
EXPORT_SYMBOLS error & printer_error()
Printer errors.
static void run(U, hid_t &dataset_out, hid_t &dataset2_out, std::vector< hid_t > &datasets, std::vector< hid_t > &datasets2, const unsigned long long size, const std::unordered_map< PPIDpair, unsigned long long, PPIDHash, PPIDEqual > &RA_write_hash, const std::vector< std::vector< unsigned long long > > &pointid, const std::vector< std::vector< unsigned long long > > &rank, const std::vector< unsigned long long > &aux_sizes, hid_t &, hid_t &)
LOCAL_INFO macro.
void Enter_HDF5(hid_t dataset, T &... params)
#define LOCAL_INFO
Definition: local_info.hpp:34
void run(bool, bool, bool, int, double, double, int, int, int, int, int, double, char[], int, int[], bool, bool, bool, bool, double, int, double(*)(double *, int, int, void *), void(*)(int, int, int, double *, double *, double *, double, double, double, void *), void *)
dataset getSimpleExtentNpoints
Definition: hdf5tools.cpp:720
std::vector< unsigned long long > sizes
std::vector< std::string > param_names
void Enter_Aux_Parameters(const std::string &output_file, bool resume=false)
std::vector< unsigned long long > cum_sizes
Definitions of new MPI datatypes needed by printers.
std::unordered_map< PPIDpair, unsigned long long, PPIDHash, PPIDEqual > get_RA_write_hash(hid_t, std::unordered_set< PPIDpair, PPIDHash, PPIDEqual > &)
Helper function to create output hash map for RA points note: left_to_match points will be erased as ...
static void run(U, hid_t &dataset, std::vector< T > &vec)
void combine_hdf5_files(const std::string output_file, const std::string &base_file_name, const std::string &group, const size_t num, const bool resume, const bool cleanup, const bool skip, const std::vector< std::string > input_files=std::vector< std::string >())
A simple C++ wrapper for the MPI C bindings.
#define DSET_SIZE_CHECK(VEC)
std::unordered_set< std::string > param_set
Exception objects required for standalone compilation.
hid_t closeSpace(hid_t space_id)
Close dataspace.
static void run(U, hid_t &dataset_out, std::vector< hid_t > &datasets, unsigned long long &size_tot, std::vector< unsigned long long > &sizes, hid_t &old_dataset, size_t offset)
std::vector< T > vec(std::vector< T > vector)
Definition: daFunk.hpp:142
A collection of tools for interacting with HDF5 databases.
std::vector< std::string > file_names
Base template is left undefined in order to raise a compile error if specialisation doesn&#39;t exist...
Definition: hdf5tools.hpp:77
pointID / process number pair Used to identify a single parameter space point
TODO: see if we can use this one:
Definition: Analysis.hpp:33
std::pair< std::vector< std::string >, std::vector< size_t > > find_temporary_files(const std::string &finalfile)
Search for temporary files to be combined.
std::pair< hid_t, hid_t > selectChunk(const hid_t dset_id, std::size_t offset, std::size_t length)
Select a simple hyperslab in a 1D dataset.
Definition: hdf5tools.cpp:733