Commit 4828f95b authored by Lukas Weber's avatar Lukas Weber
Browse files

testing and implementing missing features/tools

parent 1cdff5f2
......@@ -21,6 +21,7 @@ set(SRCs
mc.cpp
runner_single.cpp
runner_task.cpp
results.cpp
)
add_library(load_leveller ${SRCs})
......
......@@ -3,20 +3,17 @@
#include <unistd.h>
#include <sstream>
#include <iostream>
#include <sys/file.h>
static bool filter_available(H5Z_filter_t filter) {
htri_t avail = H5Zfilter_avail(filter);
if(not avail) {
if(avail == 0) {
return false;
}
unsigned int filter_info;
herr_t status = H5Zget_filter_info(filter, &filter_info);
if(status < 0 or not (filter_info & H5Z_FILTER_CONFIG_ENCODE_ENABLED) or not (filter_info & H5Z_FILTER_CONFIG_DECODE_ENABLED)) {
return false;
}
return true;
return status >= 0 && (filter_info & H5Z_FILTER_CONFIG_ENCODE_ENABLED) != 0 && (filter_info & H5Z_FILTER_CONFIG_DECODE_ENABLED) != 0;
}
static herr_t H5Ewalk_cb(unsigned int n, const H5E_error2_t *err_desc, void *client_data) {
......@@ -25,8 +22,8 @@ static herr_t H5Ewalk_cb(unsigned int n, const H5E_error2_t *err_desc, void *cli
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);
s << fmt::format(" {}: {}\n", static_cast<int>(err_desc->maj_num), maj_str);
s << fmt::format(" {}: {}\n\n", static_cast<int>(err_desc->min_num), min_str);
free(min_str);
free(maj_str);
......@@ -57,8 +54,9 @@ iodump::group::group(hid_t parent, const std::string& path) {
}
iodump::group::~group() {
if(group_ < 0)
if(group_ < 0) {
return;
}
herr_t status = H5Gclose(group_);
if(status < 0) {
std::cerr << iodump_exception{"H5Gclose"}.what();
......@@ -72,10 +70,11 @@ iodump::group::iterator iodump::group::begin() const {
}
iodump::group::iterator iodump::group::end() const {
H5G_info_t info;
H5G_info_t info{};
herr_t status = H5Gget_info(group_, &info);
if(status < 0)
if(status < 0) {
throw iodump_exception{"H5Gget_info"};
}
return iodump::group::iterator{group_, info.nlinks};
}
......@@ -86,13 +85,15 @@ iodump::group::iterator::iterator(hid_t group, uint64_t idx)
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)
if(name_size < 0) {
throw iodump_exception{"H5Lget_name_by_idx"};
}
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)
if(name_size < 0) {
throw iodump_exception{"H5Lget_name_by_idx"};
}
return std::string(buf.data());
}
......@@ -107,38 +108,42 @@ bool iodump::group::iterator::operator!=(const iterator& b) {
}
iodump iodump::create(const std::string& filename) {
H5Eset_auto(H5E_DEFAULT, NULL, NULL);
H5Eset_auto(H5E_DEFAULT, nullptr, nullptr);
hid_t file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
if(file < 0)
if(file < 0) {
throw iodump_exception{"H5Fcreate"};
}
return iodump{file};
return iodump{filename, file};
}
iodump iodump::open_readonly(const std::string& filename) {
H5Eset_auto(H5E_DEFAULT, NULL, NULL);
H5Eset_auto(H5E_DEFAULT, nullptr, nullptr);
hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
if(file < 0)
if(file < 0) {
throw iodump_exception{"H5Fopen"};
return iodump{file};
}
return iodump{filename, file};
}
iodump iodump::open_readwrite(const std::string& filename) {
H5Eset_auto(H5E_DEFAULT, NULL, NULL);
H5Eset_auto(H5E_DEFAULT, nullptr, nullptr);
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)
if(file < 0) {
throw iodump_exception{"H5Fopen"};
return iodump{file};
}
return iodump{filename, file};
}
iodump::iodump(hid_t h5_file) : h5_file_{h5_file} {
if(compression_filter_ != 0 and not filter_available(compression_filter_))
iodump::iodump(std::string filename, hid_t h5_file) : filename_{std::move(filename)}, h5_file_{h5_file} {
if(compression_filter_ != 0 && !filter_available(compression_filter_)) {
throw iodump_exception{"H5Filter not available."};
}
}
iodump::~iodump() {
......@@ -149,59 +154,51 @@ 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 {
iodump::h5_handle 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)
h5_handle dataset{H5Dopen2(group_, name.c_str(), H5P_DEFAULT), H5Dclose};
h5_handle dataspace{H5Dget_space(*dataset), H5Sclose};
hsize_t oldsize = H5Sget_simple_extent_npoints(*dataspace);
if(oldsize < 0) {
throw iodump_exception{"H5Sget_simple_extent_npoints"};
if(oldsize != size)
}
if(oldsize != size) {
throw std::runtime_error{"iodump: tried to write into an existing dataset with different dimensions!"};
}
return dataset;
} else {
hid_t dataspace;
hid_t dataspace_h;
if(not unlimited) {
dataspace = H5Screate_simple(1, &size, nullptr);
dataspace_h = H5Screate_simple(1, &size, nullptr);
} else {
hsize_t maxdim = H5S_UNLIMITED;
dataspace = H5Screate_simple(1, &size, &maxdim);
dataspace_h = H5Screate_simple(1, &size, &maxdim);
}
h5_handle dataspace{dataspace_h, H5Sclose};
hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
h5_handle plist{H5Pcreate(H5P_DATASET_CREATE), H5Pclose};
if(chunk_size > 1) { // do not use compression on small datasets
status = H5Pset_chunk(plist, 1, &chunk_size);
if(status < 0)
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)
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 h5_handle{H5Dcreate2(group_, name.c_str(), datatype, *dataspace, H5P_DEFAULT, *plist, H5P_DEFAULT), H5Dclose};
}
return dataset;
}
iodump::group iodump::group::open_group(const std::string& path) const {
......@@ -209,17 +206,13 @@ iodump::group iodump::group::open_group(const std::string& path) const {
}
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"};
hid_t dataspace = H5Dget_space(dataset);
if(dataspace < 0)
throw iodump_exception{"H5Dget_space"};
h5_handle dataset{H5Dopen2(group_, name.c_str(), H5P_DEFAULT), H5Dclose};
h5_handle dataspace{H5Dget_space(*dataset), H5Sclose};
int size = H5Sget_simple_extent_npoints(dataspace); // rank > 1 will get flattened when loaded.
if(size < 0)
int size = H5Sget_simple_extent_npoints(*dataspace); // rank > 1 will get flattened when loaded.
if(size < 0) {
throw iodump_exception{"H5Sget_simple_extent_npoints"};
}
return size;
}
......@@ -228,9 +221,39 @@ 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) {
}
if(exists < 0) {
throw iodump_exception{"H5Lexists"};
}
return true;
}
iodump::h5_handle::h5_handle(hid_t handle, herr_t (*closer)(hid_t))
: closer_{closer}, handle_{handle} {
if(handle < 0) {
throw iodump_exception{"h5_handle"};
}
}
iodump::h5_handle::h5_handle(h5_handle&& h) noexcept
: closer_{h.closer_}, handle_{h.handle_} {
h.handle_ = -1;
}
iodump::h5_handle::~h5_handle() {
if(handle_ < 0) {
return;
}
herr_t status = closer_(handle_);
if(status < 0) {
std::cerr << iodump_exception{"~h5_handle"}.what() << "\n";
std::abort();
}
}
hid_t iodump::h5_handle::operator*() {
return handle_;
}
......@@ -19,7 +19,33 @@ public:
// writing depending on how you open it. If you write to a read only file,
// there will be an error probably.
class iodump {
public:
private:
// helper class to make sure those pesky HDF5 hid_t handles are always closed
class h5_handle {
public:
h5_handle(hid_t handle, herr_t (*closer)(hid_t));
~h5_handle();
h5_handle(h5_handle&) = delete;
h5_handle(h5_handle&&) noexcept;
hid_t operator*();
private:
herr_t (*closer_)(hid_t);
hid_t handle_;
};
iodump(std::string filename, hid_t h5_file);
const std::string filename_;
const hid_t h5_file_;
// TODO: make these variable if necessary
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();
public:
// Wrapper around the concept of a HDF5 group.
// You can list over the group elements by using the iterators.
class group {
......@@ -65,7 +91,7 @@ public:
// 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;
iodump::h5_handle 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
......@@ -78,17 +104,8 @@ public:
iodump(iodump&) = delete;
~iodump();
private:
iodump(hid_t h5_file);
const hid_t h5_file_;
// TODO: make these variable if necessary
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();
friend class group;
};
template<typename T>
......@@ -124,9 +141,6 @@ constexpr hid_t iodump::h5_datatype() {
// ... or it is a native datatype I forgot to add. Then add it.
}
template<class T>
void iodump::group::write(const std::string& name, const std::vector<T>& data) const {
int chunk_size = 0;
......@@ -138,13 +152,10 @@ void iodump::group::write(const std::string& name, const std::vector<T>& data) c
compression_filter= iodump::compression_filter_;
}
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());
h5_handle 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"};
status = H5Dclose(dataset);
if(status < 0)
throw iodump_exception{"H5Dclose"};
}
template<class T>
......@@ -166,79 +177,58 @@ void iodump::group::insert_back(const std::string& name, const std::vector<T>& d
create_dataset(name, h5_datatype<T>(), 0, chunk_size_, compression_filter_, true);
}
hid_t 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"};
h5_handle dataset{H5Dopen2(group_, name.c_str(), H5P_DEFAULT), H5Dclose};
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"};
herr_t status = H5Sclose(dataspace);
if(status < 0)
throw iodump_exception{"H5Sclose"};
hsize_t new_size = size + data.size();
status = H5Dset_extent(dataset, &new_size);
if(status < 0)
throw iodump_exception{"H5Pset_extent"};
h5_handle memspace{H5Screate_simple(1, &mem_size, nullptr), H5Sclose};
int size;
herr_t status;
{ // limit the scope of the dataspace handle
h5_handle dataspace{H5Dget_space(*dataset), H5Sclose};
size = H5Sget_simple_extent_npoints(*dataspace);
if(size < 0) {
throw iodump_exception{"H5Sget_simple_extent_npoints"};
}
if(data.size() > 0) {
hsize_t new_size = size + data.size();
status = H5Dset_extent(*dataset, &new_size);
if(status < 0) {
throw iodump_exception{"H5Pset_extent"};
}
}
} // because it will be reopened after this
dataspace = H5Dget_space(dataset);
if(dataspace < 0)
throw iodump_exception{"H5Dget_space"};
h5_handle dataspace{H5Dget_space(*dataset), H5Sclose};
// select the hyperslab of the extended area
hsize_t pos = size;
hsize_t extent = data.size();
status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, &pos, nullptr, &extent, nullptr);
status = H5Sselect_hyperslab(*dataspace, H5S_SELECT_SET, &pos, nullptr, &extent, nullptr);
if(status < 0)
throw iodump_exception{"H5Sselect_hyperslap"};
status = H5Dwrite(dataset, h5_datatype<T>(), memspace, 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);
if(status < 0)
throw iodump_exception{"H5Dclose"};
}
template<class T>
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"};
hid_t dataspace = H5Dget_space(dataset);
if(dataspace < 0)
throw iodump_exception{"H5Dget_space"};
h5_handle dataset{H5Dopen2(group_, name.c_str(), H5P_DEFAULT), H5Dclose};
h5_handle dataspace{H5Dget_space(*dataset), H5Sclose};
int size = H5Sget_simple_extent_npoints(dataspace); // rank > 1 will get flattened when loaded.
int size = H5Sget_simple_extent_npoints(*dataspace); // rank > 1 will get flattened when loaded.
if(size < 0)
throw iodump_exception{"H5Sget_simple_extent_npoints"};
data.resize(size);
herr_t status = H5Dread(dataset, h5_datatype<T>(), H5S_ALL, H5P_DEFAULT, H5P_DEFAULT, data.data());
herr_t status = H5Dread(*dataset, h5_datatype<T>(), H5S_ALL, H5P_DEFAULT, H5P_DEFAULT, data.data());
if(status < 0)
throw iodump_exception{"H5Dread"};
status = H5Sclose(dataspace);
if(status < 0)
throw iodump_exception{"H5Sclose"};
status = H5Dclose(dataset);
if(status < 0)
throw iodump_exception{"H5Dclose"};
}
template<>
......
#include "evalable.h"
#include <map>
#include <fmt/format.h>
#include "measurements.h"
evalable::evalable(const std::string& name, std::vector<std::string> used_observables, func fun)
: name_{name}, used_observables_{used_observables}, fun_{fun} {
evalable::evalable(std::string name, std::vector<std::string> used_observables, func fun)
: name_{std::move(name)}, used_observables_{std::move(used_observables)}, fun_{std::move(fun)} {
// evalable names also should be valid HDF5 paths
if(not measurements::observable_name_is_legal(name)) {
throw std::runtime_error{fmt::format("illegal evalable name '{}': must not contain . or /", name)};
}
}
const std::string& evalable::name() const {
......@@ -19,19 +24,21 @@ void evalable::jackknife(const results& res, observable_result& obs_res) const {
size_t bin_count = -1; // maximal value
for(const auto& obs_name : used_observables_) {
if(res.observables.count(obs_name) <= 0) {
throw std::runtime_error(fmt::format("evalable '{}': used observable '{}' not found in Monte Carlo results."));
throw std::runtime_error(fmt::format("evalable '{}': used observable '{}' not found in Monte Carlo results.", name_, obs_name));
}
const auto& obs = res.observables.at(obs_name);
if(obs.rebinning_bin_count < bin_count)
if(obs.rebinning_bin_count < bin_count) {
bin_count = obs.rebinning_bin_count;
}
observables.emplace_back(obs);
}
obs_res.rebinning_bin_count = bin_count;
if(bin_count == 0)
if(bin_count == 0) {
return;
}
std::vector<std::vector<double>> jacked_means(observables.size());
std::vector<std::vector<double>> sums(observables.size());
......@@ -60,7 +67,7 @@ void evalable::jackknife(const results& res, observable_result& obs_res) const {
}
std::vector<double> jacked_eval = fun_(jacked_means);
if(jacked_eval_mean.size() == 0)
if(jacked_eval_mean.empty())
jacked_eval_mean.resize(jacked_eval.size(), 0);
if(jacked_eval_mean.size() != jacked_eval.size()) {
throw std::runtime_error(fmt::format("evalable '{}': evalables must not change their dimensions depending on the input"));
......@@ -78,8 +85,9 @@ void evalable::jackknife(const results& res, observable_result& obs_res) const {
// calculate the function estimator from the complete dataset.
for(size_t obs_idx = 0; obs_idx < observables.size(); obs_idx++) {
jacked_means[obs_idx] = sums[obs_idx];
for(size_t i = 0; i < observables[obs_idx].mean.size(); i++)
jacked_means[obs_idx][i] /= bin_count;
for(auto& jacked_mean : jacked_means[obs_idx]) {
jacked_mean /= bin_count;
}
}
std::vector<double> complete_eval = fun_(jacked_means);
......@@ -103,15 +111,16 @@ void evalable::jackknife(const results& res, observable_result& obs_res) const {
}
std::vector<double> jacked_eval = fun_(jacked_means);
if(obs_res.error.size() == 0)
if(obs_res.error.empty()) {
obs_res.error.resize(jacked_eval.size(), 0);
}
for(size_t i = 0; i < obs_res.error.size(); i++) {
obs_res.error[i] += pow(jacked_eval[i]-jacked_eval_mean[i], 2);
}
}
for(size_t i = 0; i < obs_res.error.size(); i++) {
obs_res.error[i] = sqrt((bin_count-1)/bin_count * obs_res.error[i]);
obs_res.error[i] = sqrt((bin_count-1) * obs_res.error[i]/bin_count);
}
}
......
......@@ -16,7 +16,7 @@ public:
//
typedef std::function<std::vector<double>(const std::vector<std::vector<double>>& observables)> func;
evalable(const std::string& name, std::vector<std::string> used_observables, func fun);
evalable(std::string name, std::vector<std::string> used_observables, func fun);
const std::string& name() const;
void jackknife(const results& res, observable_result &out) const;
......
......@@ -6,17 +6,29 @@
#include "mc.h"
namespace load_leveller {
template <class mc_runner>
static int run_mc(mc_factory mccreator, int argc, char **argv) {
inline int merge_only(jobinfo job, const mc_factory& mccreator, int argc, char **argv) {
for(size_t task_id = 0; task_id < job.task_names.size(); task_id++) {
std::vector<evalable> evalables;
std::unique_ptr<abstract_mc> sys{mccreator(job.jobfile_name, job.task_names[task_id])};
sys->register_evalables(evalables);
job.merge_task(task_id, evalables);
std::cout << fmt::format("-- {} merged\n", job.taskdir(task_id));
}
return 0;
}
inline int run_mc(int (*starter)(jobinfo job, const mc_factory&, int argc, char **argv), mc_factory mccreator, int argc, char **argv) {
if(argc < 2) {
std::cerr << fmt::format("{0} JOBFILE\n{0} single JOBFILE\n{0} merge JOBFILE\n\n Without further flags, the MPI scheduler is started. 'single' starts a single core scheduler (useful for debugging), 'merge' merges unmerged measurement results.\n", argv[0]);
return -1;
}
std::string jobfile{argv[1]};
jobinfo job{jobfile};
mc_runner r;
return r.start(jobfile, mccreator, argc, argv);