writing measurements to h5 seems to work now

MersenneTwister.h

@@ -127,8 +127,8 @@ public:
 	uint32 get_seed() {return myoneseed;}
 	
 	// Saving and loading generator state
-	void save(std::vector<uint32>& saveArray) const;  // to array of size SAVE
-	void load(const std::vector<uint32>& loadArray);  // from such array
+	void save(std::vector<uint32>& saveArray) const;
+	void load(const std::vector<uint32>& loadArray);
 	friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
 	friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
 
@@ -372,12 +372,16 @@ inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
 inline void MTRand::save(std::vector<uint32>& saveArray) const
 {
 	saveArray = std::vector<uint32>{state, state+N};
+	saveArray.push_back(left);
 }
 
 
 inline void MTRand::load(const std::vector<uint32>& loadArray) {
-	assert(loadArray.size() == N);
+	left = loadArray.back();
+	assert(loadArray.size() == N+1);
 	memcpy(state, loadArray.data(), sizeof(uint32)*N);
+
+	pNext = &state[N-left];
 }
 
 

dump.cpp

 #include "dump.h"
 #include <typeinfo>
+#include <unistd.h>
+#include <sstream>
+#include <iostream>
 
 static bool filter_available(H5Z_filter_t filter) {
 	htri_t avail = H5Zfilter_avail(filter);
 	if(not avail) {
 		return false;
-	}
+        }
 
 	unsigned int filter_info;
 	herr_t status = H5Zget_filter_info(filter, &filter_info);
@@ -16,33 +19,77 @@ static bool filter_available(H5Z_filter_t filter) {
 	return true;
 }
 
-iodump::group_members::group_members(hid_t group)
-	: group_{group} {
+static herr_t H5Ewalk_cb(unsigned int n, const H5E_error2_t *err_desc, void *client_data) {
+	std::stringstream& s = *reinterpret_cast<std::stringstream *>(client_data);
+
+	char *min_str = H5Eget_minor(err_desc->min_num);
+	char *maj_str = H5Eget_major(err_desc->maj_num);
+	s << fmt::format("#{}: {}:{} in {}(): {}\n", n, err_desc->file_name, err_desc->line, err_desc->func_name, err_desc->desc);
+	s << fmt::format("    {}: {}\n", err_desc->maj_num, maj_str);
+	s << fmt::format("    {}: {}\n\n", err_desc->min_num, min_str);
+
+	free(min_str);
+	free(maj_str);
+
+	return 0;
+}
+
+iodump_exception::iodump_exception(const std::string& message)  {
+	std::stringstream s;
+	H5Ewalk(H5E_DEFAULT, H5E_WALK_DOWNWARD, H5Ewalk_cb, &s);
+
+	s << "Error triggered: " << message;
+	message_ = s.str();
+}
+
+const char *iodump_exception::what() const noexcept {
+	return message_.c_str();
 }
 
-iodump::group_members::iterator iodump::group_members::begin() {
-	return iodump::group_members::iterator{group_, 0};
+iodump::group::group(hid_t parent, const std::string& path) {
+	group_ = H5Gopen(parent, path.c_str(), H5P_DEFAULT);
+	if(group_ < 0) {
+		group_ = H5Gcreate2(parent, path.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+		if(group_ < 0) {
+			throw iodump_exception{"H5Gcreate"};
+		}
+	}
 }
 
-iodump::group_members::iterator iodump::group_members::end() {
+iodump::group::~group() {
+	if(group_ < 0)
+		return;
+	herr_t status = H5Gclose(group_);
+	if(status < 0) {
+		std::cerr << iodump_exception{"H5Gclose"}.what();
+		std::cerr << group_ << " …something went wrong in the destructor.\n";
+		assert(false);
+	}
+}
+
+iodump::group::iterator iodump::group::begin() const {
+	return iodump::group::iterator{group_, 0};
+}
+
+iodump::group::iterator iodump::group::end() const {
 	H5G_info_t info;
 	herr_t status = H5Gget_info(group_, &info);
 	if(status < 0)
 		throw iodump_exception{"H5Gget_info"};
 	
-	return iodump::group_members::iterator{group_, info.nlinks};
+	return iodump::group::iterator{group_, info.nlinks};
 }
 
-iodump::group_members::iterator::iterator(hid_t group, uint64_t idx)
+iodump::group::iterator::iterator(hid_t group, uint64_t idx)
 	: group_(group), idx_(idx) {
 }
 
-std::string iodump::group_members::iterator::operator*() {
+std::string iodump::group::iterator::operator*() {
 	ssize_t name_size = H5Lget_name_by_idx(group_, "/", H5_INDEX_NAME, H5_ITER_INC, idx_, nullptr, 0, H5P_DEFAULT);
 	if(name_size < 0)
 		throw iodump_exception{"H5Lget_name_by_idx"};
 
-	std::vector<char> buf(name_size);
+	std::vector<char> buf(name_size+1);
 	name_size = H5Lget_name_by_idx(group_, "/", H5_INDEX_NAME, H5_ITER_INC, idx_, buf.data(), buf.size(), H5P_DEFAULT);
 	if(name_size < 0)
 		throw iodump_exception{"H5Lget_name_by_idx"};
@@ -50,20 +97,17 @@ std::string iodump::group_members::iterator::operator*() {
 	return std::string(buf.data());
 }
 
-iodump::group_members::iterator iodump::group_members::iterator::operator++() {
+iodump::group::iterator iodump::group::iterator::operator++() {
 	idx_++;
 	return *this;
 }
 	
-bool iodump::group_members::iterator::operator!=(const iterator& b) {
+bool iodump::group::iterator::operator!=(const iterator& b) {
 	return idx_ != b.idx_;
 }
 
-iodump::group_members iodump::list() {
-	return group_members{group_};
-}
-	
 iodump iodump::create(const std::string& filename) {
+	H5Eset_auto(H5E_DEFAULT, NULL, NULL);
 	hid_t file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
 
 	if(file < 0)
@@ -73,6 +117,7 @@ iodump iodump::create(const std::string& filename) {
 }
 
 iodump iodump::open_readonly(const std::string& filename) {
+	H5Eset_auto(H5E_DEFAULT, NULL, NULL);
 	hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
 	if(file < 0)
 		throw iodump_exception{"H5Fopen"};
@@ -80,6 +125,11 @@ iodump iodump::open_readonly(const std::string& filename) {
 }
 
 iodump iodump::open_readwrite(const std::string& filename) {
+	H5Eset_auto(H5E_DEFAULT, NULL, NULL);
+	if(access(filename.c_str(), R_OK) != F_OK) {
+		create(filename);
+	}
+
 	hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
 	if(file < 0)
 		throw iodump_exception{"H5Fopen"};
@@ -87,32 +137,78 @@ iodump iodump::open_readwrite(const std::string& filename) {
 }
 
 iodump::iodump(hid_t h5_file) : h5_file_{h5_file} {
-	if(not filter_available(compression_filter_))
+	if(compression_filter_ != 0 and not filter_available(compression_filter_))
 		throw iodump_exception{"H5Filter not available."};
-
-	group_ = H5Gopen(h5_file_, "/", H5P_DEFAULT);
-	if(group_ < 0)
-		throw iodump_exception{"H5Gopen"};
 }
 
 iodump::~iodump() {
-	H5Gclose(group_);
 	H5Fclose(h5_file_);
 }
 
-void iodump::change_group(const std::string& path) {
-	hid_t tmp = H5Gopen(group_, path.c_str(), H5P_DEFAULT);
-	if(tmp < 0)
-		throw iodump_exception{"H5Gopen"};
-	herr_t status = H5Gclose(group_);
-	if(status < 0)
-		throw iodump_exception{"H5Gclose"};
-	group_ = tmp;
-}
+iodump::group iodump::get_root() {
+	return group{h5_file_, "/"};
+}
+
+hid_t iodump::group::create_dataset(const std::string& name, hid_t datatype, hsize_t size, hsize_t chunk_size, H5Z_filter_t compression_filter, bool unlimited) const {
+	herr_t status;
+
+	hid_t dataset;
+	if(exists(name)) {
+		dataset = H5Dopen2(group_, name.c_str(), H5P_DEFAULT);
+		if(dataset < 0)
+			throw iodump_exception{"H5Dopen2"};
+		hid_t dataspace = H5Dget_space(dataset);
+		if(dataspace < 0)
+			throw iodump_exception{"H5Dget_space"};
+
+		hsize_t oldsize = H5Sget_simple_extent_npoints(dataspace);
+		if(oldsize < 0)
+			throw iodump_exception{"H5Sget_simple_extent_npoints"};
+		if(oldsize != size)
+			throw std::runtime_error{"iodump: tried to write into an existing dataset with different dimensions!"};
+	} else {
+		hid_t dataspace;
+		if(not unlimited) {
+			dataspace = H5Screate_simple(1, &size, nullptr);
+		} else {
+			hsize_t maxdim = H5S_UNLIMITED;
+			dataspace = H5Screate_simple(1, &size, &maxdim);
+		}
+
+		hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
+
+		if(chunk_size > 1) { // do not use compression on small datasets
+			status = H5Pset_chunk(plist, 1, &chunk_size);
+			if(status < 0)
+				throw iodump_exception{"H5Pset_chunk"};
+			
+			if(compression_filter == H5Z_FILTER_DEFLATE) {
+				status = H5Pset_deflate(plist, 6);
+				if(status < 0)
+					throw iodump_exception{"H5Pset_deflate"};
+			}
+		}
+
+		dataset = H5Dcreate2(group_, name.c_str(), datatype, dataspace, H5P_DEFAULT, plist, H5P_DEFAULT);
+		if(dataset < 0)
+			throw iodump_exception{"H5Dcreate2"};
+		status = H5Pclose(plist);
+		if(status < 0)
+			throw iodump_exception{"H5Pclose"};
+		status = H5Sclose(dataspace);
+		if(status < 0)
+			throw iodump_exception{"H5Sclose"};
+	}
+	
 
+	return dataset;
+}
 
+iodump::group iodump::group::open_group(const std::string& path) const {
+	return group{group_, path};
+}
 
-size_t iodump::get_extent(const std::string& name) {
+size_t iodump::group::get_extent(const std::string& name) const {
 	hid_t dataset = H5Dopen2(group_, name.c_str(), H5P_DEFAULT);
 	if(dataset < 0)
 		throw iodump_exception{"H5Dopen2"};
@@ -127,3 +223,14 @@ size_t iodump::get_extent(const std::string& name) {
 
 	return size;
 }
+
+bool iodump::group::exists(const std::string &path) const {
+	htri_t exists = H5Lexists(group_, path.c_str(), H5P_DEFAULT);
+	if(exists == 0) {
+		return false;
+	} else if(exists < 0) {
+		throw iodump_exception{"H5Lexists"};
+	}
+
+	return true;
+}

dump.h

@@ -5,10 +5,13 @@
 #include <string>
 #include <fmt/format.h>
 
-struct iodump_exception : public std::exception {
-	std::string message;
-	iodump_exception(const std::string& msg) : message(msg) {}
-	const char *what() { return (std::string("iodump: ")+message).c_str(); }
+class iodump_exception : public std::exception {
+private:
+	std::string message_;
+public:
+	iodump_exception(const std::string& msg);
+
+	const char *what() const noexcept;
 };
 
 
@@ -17,8 +20,9 @@ struct iodump_exception : public std::exception {
 // there will be an error probably.
 class iodump {
 public:
-	// helper for comfortable h5 entry listings.
-	class group_members {
+	// Wrapper around the concept of a HDF5 group.
+	// You can list over the group elements by using the iterators.
+	class group {
 	public:
 		struct iterator {
 			iterator(hid_t group, uint64_t idx);
@@ -30,11 +34,38 @@ public:
 			uint64_t idx_;
 		};
 			
-		group_members(hid_t group);
-		iterator begin();
-		iterator end();
+		group(hid_t parent, const std::string& path);
+		group(const group&) = delete; // did you know this is a thing? Very handy in preventing errors.
+		~group();
+		iterator begin() const;
+		iterator end() const;
+		
+		template <class T>
+		void write(const std::string& name, const std::vector<T>& data) const;
+		template <class T>
+		void write(const std::string& name, const T& value) const; // this is meant for atomic datatypes like int and double
+
+		// insert_back inserts data at the end of the dataset given by name, extending it if necessary.
+		// This only works in read/write mode.
+		template <class T>
+		void insert_back(const std::string& name, const std::vector<T>& data) const;
+		
+		template <class T>
+		void read(const std::string& name, std::vector<T>& data) const;
+		template <class T>
+		void read(const std::string& name, T& value) const;
+		
+		size_t get_extent(const std::string& name) const;
+
+		group open_group(const std::string &path) const; // this works like the cd command
+
+		bool exists(const std::string &path) const; // checks whether an object in the dump file exists
 	private:
-		const hid_t group_;
+		hid_t group_;
+
+		// chunk_size == 0 means contiguous storage
+		// if the dataset already exists, we try to overwrite it. However it must have the same extent for that to work.
+		hid_t create_dataset(const std::string& name, hid_t datatype, hsize_t size, hsize_t chunk_size, H5Z_filter_t compression_filter, bool unlimited) const;
 	};
 	
 	// delete what was there and create a new file for writing
@@ -43,47 +74,27 @@ public:
 	static iodump open_readonly(const std::string& filename);
 	static iodump open_readwrite(const std::string& filename);
 
-	// returns a range expression over the entries of the current group
-	// which you can iterate over using a range-for loop.
-	group_members list();
+	group get_root();
 
+	iodump(iodump&) = delete;
 	~iodump();
-	template <class T>
-	void write(const std::string& name, const std::vector<T>& data);
-	template <class T>
-	void write(const std::string& name, T value); // this is meant for atomic datatypes like int and double
-
-	// insert_back inserts data at the end of the dataset given by name, extending it if necessary.
-	// This only works in read/write mode.
-	template <class T>
-	void insert_back(const std::string& name, const std::vector<T>& data);
-	
-	template <class T>
-	void read(const std::string& name, std::vector<T>& data);
-	template <class T>
-	void read(const std::string& name, T& value);
-	size_t get_extent(const std::string& name);
-
-	void change_group(const std::string &path); // this works like the cd command
 private:
 	iodump(hid_t h5_file);
-
-	template<typename T>
-	static hid_t h5_datatype();
-	template<typename T>
-	static hid_t create_dataset(hid_t group, const std::string& name, hsize_t size, hsize_t chunk_size, H5Z_filter_t compression_filter, bool unlimited);
 	
 	const hid_t h5_file_;
-	hid_t group_;
 
 	// TODO: make these variable if necessary
-	const H5Z_filter_t compression_filter_ = H5Z_FILTER_DEFLATE;
-	const size_t chunk_size_ = 1000;
+	static const H5Z_filter_t compression_filter_ = H5Z_FILTER_DEFLATE;
+	static const size_t chunk_size_ = 1000;
+		
+	template<typename T>
+	constexpr static hid_t h5_datatype();
 };
 
-
 template<typename T>
-hid_t iodump::h5_datatype() {
+constexpr hid_t iodump::h5_datatype() {
+	if(typeid(T) == typeid(char))
+		return H5T_NATIVE_CHAR;
 	if(typeid(T) == typeid(int))
 		return H5T_NATIVE_INT;
 	if(typeid(T) == typeid(short))
@@ -104,63 +115,29 @@ hid_t iodump::h5_datatype() {
 		return H5T_NATIVE_FLOAT;
 	if(typeid(T) == typeid(double))
 		return H5T_NATIVE_DOUBLE;
+		
 
-	throw iodump_exception{fmt::format("unsupported datatype: {}",typeid(T).name())};
+	throw std::runtime_error{fmt::format("unsupported datatype: {}",typeid(T).name())};
+	// If you run into this error, you probably tried to write a non-primitive datatype
+	// to a dump file. see runner_task for a minimalistic example of what to do.
+	// ... or it is a native datatype I forgot to add. Then add it.
 }
 
-// chunk_size == 0 means contiguous storage
-template<class T>
-hid_t iodump::create_dataset(hid_t group, const std::string& name, hsize_t size, hsize_t chunk_size, H5Z_filter_t compression_filter, bool unlimited) {
-	herr_t status;
-
-	hid_t dataspace;
-	if(not unlimited) {
-		dataspace = H5Screate_simple(1, &size, nullptr);
-	} else {
-		hsize_t maxdim = H5S_UNLIMITED;
-		dataspace = H5Screate_simple(1, &size, &maxdim);
-	}
-
-	hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
 
-	if(chunk_size > 1) { // do not use compression on small datasets
-		status = H5Pset_chunk(plist, 1, &chunk_size);
-		if(status < 0)
-			throw iodump_exception{"H5Pset_chunk"};
-		
-		if(compression_filter== H5Z_FILTER_DEFLATE) {
-			status = H5Pset_deflate(plist, 6);
-			if(status < 0)
-				throw iodump_exception{"H5Pset_deflate"};
-		}
-	}
-
-	hid_t dataset = H5Dcreate2(group, name.c_str(), h5_datatype<T>(), dataspace, H5P_DEFAULT, plist, H5P_DEFAULT);
-	if(dataset < 0)
-		throw iodump_exception{"H5Dcreate2"};
-
-	status = H5Pclose(plist);
-	if(status < 0)
-		throw iodump_exception{"H5Pclose"};
-	status = H5Sclose(dataspace);
-	if(status < 0)
-		throw iodump_exception{"H5Sclose"};
-	return dataset;
-}
 	
 
 template<class T>
-void iodump::write(const std::string& name, const std::vector<T>& data) {
+void iodump::group::write(const std::string& name, const std::vector<T>& data) const {
 	int chunk_size = 0;
-	H5Z_filter_t compression_filter= 0;
+	H5Z_filter_t compression_filter = 0;
 	
 	// no compression and chunking unless dataset is big enough
 	if(data.size() >= chunk_size_) { 
-		chunk_size = chunk_size_;
-		compression_filter= compression_filter_;
+		chunk_size = iodump::chunk_size_;
+		compression_filter= iodump::compression_filter_;
 	}
 
-	hid_t dataset = create_dataset<T>(group_, name, data.size(), chunk_size, compression_filter, false);
+	hid_t dataset = create_dataset(name, h5_datatype<T>(), data.size(), chunk_size, compression_filter, false);
 	herr_t status = H5Dwrite(dataset, h5_datatype<T>(), H5S_ALL, H5S_ALL, H5P_DEFAULT, data.data());
 	if(status < 0)
 		throw iodump_exception{"H5Dwrite"};
@@ -170,19 +147,22 @@ void iodump::write(const std::string& name, const std::vector<T>& data) {
 }
 
 template<class T>
-void iodump::write(const std::string& name, T value) {
+void iodump::group::write(const std::string& name, const T& value) const {
 	// I hope nobody copies a lot of small values...
 	write(name, std::vector<T>{value});
 }
 
+template<>
+inline void iodump::group::write(const std::string& name, const std::string& value) const {
+	write(name, std::vector<char>{value.begin(), value.end()});
+}
+
+
 template <class T>
-void iodump::insert_back(const std::string& name, const std::vector<T>& data) {
+void iodump::group::insert_back(const std::string& name, const std::vector<T>& data) const {
 	// If the dataset does not exist, we create a new unlimited one with 0 extent.
-	htri_t exists = H5Lexists(group_, name.c_str(), H5P_DEFAULT);
-	if(exists == 0) {
-		create_dataset<T>(group_, name, 0, chunk_size_, compression_filter_, true);
-	} else if(exists < 0) {
-		throw iodump_exception{"H5Lexists"};
+	if(!exists(name)) {
+		create_dataset(name, h5_datatype<T>(), 0, chunk_size_, compression_filter_, true);
 	}
 
 	hid_t dataset = H5Dopen2(group_, name.c_str(), H5P_DEFAULT);
@@ -193,6 +173,9 @@ void iodump::insert_back(const std::string& name, const std::vector<T>& data) {
 	if(dataspace < 0)
 		throw iodump_exception{"H5Dget_space"};
 
+	hsize_t mem_size = data.size();
+	hid_t memspace = H5Screate_simple(1, &mem_size, nullptr);
+
 	int size = H5Sget_simple_extent_npoints(dataspace);
 	if(size < 0)
 		throw iodump_exception{"H5Sget_simple_extent_npoints"};
@@ -217,10 +200,13 @@ void iodump::insert_back(const std::string& name, const std::vector<T>& data) {
 	if(status < 0)
 		throw iodump_exception{"H5Sselect_hyperslap"};
 
-	status = H5Dwrite(dataset, h5_datatype<T>(), H5S_ALL, dataspace, H5P_DEFAULT, data.data());
+	status = H5Dwrite(dataset, h5_datatype<T>(), memspace, dataspace, H5P_DEFAULT, data.data());
 	if(status < 0)
 		throw iodump_exception{"H5Dwrite"};
 	status = H5Sclose(dataspace);
+	if(status < 0)
+		throw iodump_exception{"H5Sclose"};
+	status = H5Sclose(memspace);
 	if(status < 0)
 		throw iodump_exception{"H5Sclose"};
 	status = H5Dclose(dataset);
@@ -229,7 +215,7 @@ void iodump::insert_back(const std::string& name, const std::vector<T>& data) {
 }
 
 template<class T>
-void iodump::read(const std::string& name, std::vector<T>& data) {
+void iodump::group::read(const std::string& name, std::vector<T>& data) const {
 	hid_t dataset = H5Dopen2(group_, name.c_str(), H5P_DEFAULT);
 	if(dataset < 0)
 		throw iodump_exception{"H5Dopen2"};
@@ -254,9 +240,17 @@ void iodump::read(const std::string& name, std::vector<T>& data) {
 		throw iodump_exception{"H5Dclose"};
 }
 
+template<>
+inline void iodump::group::read(const std::string& name, std::string& value) const {
+	std::vector<char> buf;
+	read(name, buf);
+	value = std::string{buf.begin(), buf.end()};
+}
+
 template<class T>
-void iodump::read(const std::string& name, T& value) {
+void iodump::group::read(const std::string& name, T& value) const {
 	std::vector<T> buf;
 	read(name, buf);
 	value = buf.at(0);
 }
+

mc.cpp

@@ -16,15 +16,6 @@ void abstract_mc::random_init() {
 		rng.reset(new randomnumbergenerator());
 }
 
-void abstract_mc::random_write(iodump& d) {
-	rng->checkpoint_write(d);
-}
-
-void abstract_mc::random_read(iodump& d) {
-	rng.reset(new randomnumbergenerator());
-	rng->checkpoint_read(d);
-}
-
 double abstract_mc::random01() {
 	return rng->d();
 }
@@ -45,27 +36,28 @@ void abstract_mc::_do_update() {
 
 void abstract_mc::_write(const std::string& dir) {
 	iodump meas_file = iodump::open_readwrite(dir+".meas.h5");
-	measure.samples_write(meas_file);
-	
-	iodump dump_file = iodump::create(dir+".dump.h5");
-	measure.checkpoint_write(dump_file);
-	random_write(dump_file);
+	measure.samples_write(meas_file.get_root());
 
-	dump_file.write("sweeps", sweep_);
+	iodump dump_file = iodump::create(dir+".dump.h5");
+	auto g = dump_file.get_root();
 	
-	dump_file.write("thermalization_sweeps", std::min(therm_,sweep_)); // only for convenience
-	checkpoint_write(dump_file);
+	measure.checkpoint_write(g.open_group("measurements"));
+	rng->checkpoint_write(g.open_group("random_number_generator"));
+	checkpoint_write(g.open_group("simulation"));
+
+	g.write("sweeps", sweep_);
+	g.write("thermalization_sweeps", std::min(therm_,sweep_)); // only for convenience
 }
 
 bool abstract_mc::_read(const std::string& dir) {
 	try {
-		iodump dump_file = iodump::open_readonly(dir+"dump.h5");
-		measure.checkpoint_read(dump_file);
-		random_read(dump_file);
-		dump_file.read("sweeps", sweep_);
-
-		checkpoint_read(dump_file);