Commit e7957243 authored by Roel Aaij's avatar Roel Aaij
Browse files

Reorganise code, cleanup build, add some tests.

parent 726f9965
......@@ -5,11 +5,17 @@ list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake) # for find_package
option(ENABLE_PYTHON "Enable python bindings" TRUE)
option(ENABLE_VC "Enable usage of VC library" TRUE)
option(ENABLE_TESTS "Enable tests" TRUE)
include(ExternalProject)
include(FindPackageHandleStandardArgs)
if (ENABLE_VC)
find_package(Vc)
if (NOT Vc_FOUND or Vc_VERSION VERSION_LESS "1.3.2")
if (NOT Vc_FOUND OR "${Vc_VERSION}" VERSION_LESS "1.3.2")
message(STATUS "Using builtin Vc")
set(EXTERNAL_DIR "external")
set(Vc_VERSION "1.4.1")
......@@ -19,8 +25,6 @@ set(Vc_DESTDIR "${CMAKE_BINARY_DIR}/${EXTERNAL_DIR}")
set(Vc_LIBNAME "${CMAKE_STATIC_LIBRARY_PREFIX}Vc${CMAKE_STATIC_LIBRARY_SUFFIX}")
set(Vc_LIBRARY "${Vc_DESTDIR}/lib/${Vc_LIBNAME}")
include(ExternalProject)
include(FindPackageHandleStandardArgs)
ExternalProject_Add(VC
URL ${Vc_SRC_URI}
......@@ -73,8 +77,7 @@ ExternalProject_Add(vectorclass
SOURCE_DIR vectorclass
INSTALL_DIR ${vectorclass_ROOTDIR}
CONFIGURE_COMMAND ""
UPDATE_COMMAND unzip -o special.zip
BUILD_COMMAND ""
BUILD_COMMAND unzip -d <SOURCE_DIR> -o <SOURCE_DIR>/special.zip
INSTALL_COMMAND ${CMAKE_COMMAND} -E copy_directory <SOURCE_DIR> <INSTALL_DIR>
STEP_TARGETS install
)
......@@ -97,7 +100,9 @@ if(ENABLE_PYTHON)
# pybind11
find_package(pybind11)
if (NOT pybind11_FOUND or pybind11_VERSION VERSION_LESS "2.1.0")
if (NOT pybind11_FOUND OR pybind11_VERSION VERSION_LESS "2.1.0")
message(STATUS "Using builtin pybind11")
set(BUILTIN_PYBIND TRUE)
set(PYBIND11_VERSION "2.2.4")
set(PYBIND11_DESTDIR "${CMAKE_BINARY_DIR}/${EXTERNAL_DIR}")
......@@ -130,7 +135,8 @@ find_package_handle_standard_args(pybind11
VERSION_VAR PYBIND11_VERSION)
# install(DIRECTORY ${PYBIND11_DESTDIR}/ DESTINATION ".")
else()
set(BUILTIN_PYBIND FALSE)
endif()
else()
set(pybind11_FOUND FALSE)
......@@ -147,14 +153,11 @@ set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG")
set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g -DDEBUG")
set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lpthread")
find_package(Threads REQUIRED)
SET(CMAKE_CXX_FLAGS "${CMAKE_C_FLAGS} ${INCLUDES} -W -O3 -DNDEBUG")
SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} " )
include(VcMacros)
vc_set_preferred_compiler_flags()
message(STATUS "${Vc_INCLUDE_DIR}")
if (USE_AVX2)
include_directories(${Vc_INCLUDE_DIR})
......@@ -177,16 +180,27 @@ endif()
include_directories(${CMAKE_SOURCE_DIR}/lib)
include_directories(${CMAKE_SOURCE_DIR}/lib/generate)
add_library(generate SHARED ${generate_SOURCES})
if (USE_AVX2)
target_link_libraries(generate PUBLIC ${Vc_LIBRARIES})
target_include_directories(generate
PUBLIC ${CMAKE_BINARY_DIR}/include/vectorclass)
add_dependencies(generate vectorclass-install)
target_compile_definitions(generate PUBLIC "-DVc_IMPL=AVX2")
endif()
if(pybind11_FOUND)
include(pybind11Tools)
pybind11_add_module(k40gen ${generate_SOURCES} "src/k40gen/bindings.cpp")
add_dependencies(k40gen pybind11)
pybind11_add_module(k40gen "src/k40gen/bindings.cpp")
target_include_directories(k40gen PRIVATE ${PYTHON_SITE_PACKAGES}/numpy/core/include)
target_link_libraries(k40gen PUBLIC generate)
if(BUILTIN_PYBIND)
add_dependencies(k40gen pybind11)
endif()
endif()
if (USE_AVX2)
target_link_libraries(k40gen PUBLIC ${Vc_LIBRARIES})
target_include_directories(k40gen
PUBLIC ${CMAKE_BINARY_DIR}/include/vectorclass)
add_dependencies(k40gen vectorclass)
target_compile_definitions(k40gen PUBLIC "-DVc_IMPL=AVX2")
if(ENABLE_TESTS)
enable_testing()
add_subdirectory(tests)
endif()
......@@ -39,9 +39,6 @@ namespace Constants {
const float tot_sigma = 2.44078f;
const int random_method = 2;
const std::array<float, 4> rates{7000., 700, 70., 7.};
const float tau_l0 = 1.0e9f / rates[0];
}
float cross_prob(const float ct);
......@@ -56,7 +53,7 @@ public:
Generators(const int seed0, const int seed1, const std::array<float, 4> rates);
std::mt19937_64 mt;
double coincidence_rate;
const double coincidence_rate = 0.;
std::geometric_distribution<long> coincidence;
std::uniform_real_distribution<double> flat;
std::array<float, 31> prob1D;
......@@ -71,6 +68,20 @@ public:
return m_seeds;
}
float tau_l0() const {
return 1.0e9f / m_rates[0];
}
size_t n_expect(const long interval) const {
float total_rate = std::accumulate(begin(m_rates), end(m_rates), 0.f);
// Fudge with 30% extra space to avoid reallocating
double n_per_pmt = 1.3 * total_rate * (double)(interval) / 1e9;
if (n_per_pmt > std::numeric_limits<float>::max()) {
throw std::domain_error{"rate of " + std::to_string(n_per_pmt) + " is too large"};
}
return std::lround(storage::n_dom * storage::n_mod * (float)n_per_pmt);
}
};
template<typename Container>
......@@ -99,7 +110,8 @@ unsigned int random_index(const Container& buffer, const double random) {
}
}
size_t fill_coincidences(storage_t& times, size_t idx, const long time_start, const long time_end,
size_t fill_coincidences(storage_t& times, size_t& idx,
const long time_start, const long time_end,
Generators& gens);
std::tuple<storage_t, storage_t> generate(const long start, const long end,
......
......@@ -33,7 +33,7 @@ namespace storage {
#ifdef HAVE_CUDA
using storage_t = thrust::host_vector<long, thrust::cuda::experimental::pinned_allocator<int>>;
#else
using storage_t = std::vector<long>;
using storage_t = std::vector<long, aligned_allocator<long>>;
#endif
using queue_t = std::vector<std::tuple<storage_t, storage_t>>;
......@@ -21,8 +21,6 @@
#include <vector>
#include <tuple>
#include <optional>
#include <generate.h>
#include <storage.h>
......@@ -49,8 +47,9 @@ Generators::Generators(const int seed0, const int seed1, const std::array<float,
m_rates{std::move(rates)},
m_seed_seq{{seed0, seed1}},
mt{m_seed_seq},
coincidence_rate{std::accumulate(std::next(begin(Constants::rates)),
end(Constants::rates), 0.0)},
coincidence_rate{std::accumulate(std::next(begin(m_rates)),
end(m_rates), 0.0)},
// p = 1 - exp(-lambda)
coincidence{1. - exp(-coincidence_rate / 1e9)},
flat{0., 1.}
{
......@@ -60,6 +59,10 @@ Generators::Generators(const int seed0, const int seed1, const std::array<float,
using namespace GenScalar;
#endif
if (m_rates[0] < 1.f) {
throw std::domain_error{"L0 rate must be > 1 Hz"};
}
// probabilities for coincidences
for (size_t i = 0; i < prob1D.size(); ++i) {
prob2D[i][i] = 0.0;
......@@ -81,7 +84,8 @@ float cross_prob(const float ct) {
return std::exp(ct * (Constants::p2 + ct * (Constants::p3 + ct * Constants::p4)));
}
size_t fill_coincidences(storage_t& times, size_t idx, const long time_start, const long time_end,
size_t fill_coincidences(storage_t& times, size_t& idx,
const long time_start, const long time_end,
Generators& gen) {
const auto& prob1D = gen.prob1D;
const auto& prob2D = gen.prob2D;
......@@ -89,8 +93,14 @@ size_t fill_coincidences(storage_t& times, size_t idx, const long time_start, co
auto& mt = gen.mt;
auto& flat = gen.flat;
if (gen.coincidence_rate < 0.001) {
return idx;
}
// Fill coincidences
size_t n = 0;
for (long t1 = time_start ; t1 < time_end; t1 += gen.coincidence(mt)) {
++n;
// generate two-fold coincidence
const unsigned int pmt1 = random_index(prob1D, flat(mt));
const unsigned int pmt2 = random_index(prob2D[pmt1], flat(mt));
......@@ -119,18 +129,21 @@ size_t fill_coincidences(storage_t& times, size_t idx, const long time_start, co
}
catch (const std::domain_error&) {}
}
return idx;
return n;
}
std::tuple<storage_t, storage_t> generate(const long start, const long end,
Generators& gens, bool use_avx2) {
#ifdef USE_AVX2
if (use_avx2) {
std::cout << "Generating AVX2" << std::endl;
return generate_avx2(start, end, gens);
} else {
std::cout << "Generating scalar" << std::endl;
return generate_scalar(start, end, gens);
}
#else
std::cout << "Generating scalar" << std::endl;
return generate_scalar(start, end, gens);
#endif
}
......@@ -67,12 +67,13 @@ std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long
// Assume times.size() is multiple of 8 here
assert(storage::n_hits % 16 == 0);
storage_t times; times.resize(storage::n_hits);
storage_t values; values.resize(storage::n_hits);
times[0] = 0l;
const size_t n_expect = gens.n_expect(time_end - time_start);
const float tau_l0 = gens.tau_l0();
size_t storage_size = n_expect + long_v::size() - n_expect % long_v::size();
storage_t mod_times; mod_times.resize(10000);
storage_t mod_values; mod_values.resize(10000);
storage_t times; times.resize(storage_size);
storage_t values; values.resize(storage_size);
times[0] = 0l;
auto int_v_to_long_v = [] (const int_v& in) -> pair<long_v, long_v>
{
......@@ -97,7 +98,7 @@ std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long
// We have to rescale our drawn random numbers here to make sure
// we have the right coverage
r = -1.f * Constants::tau_l0 * log(r);
r = -1.f * tau_l0 * log(r);
auto tmp_i = simd_cast<int_v>(r + 0.5f);
// Prefix sum from stackoverflow:
......@@ -123,7 +124,7 @@ std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long
}
// Coincidences
idx = fill_coincidences(times, idx, mod_start, time_end, gens);
fill_coincidences(times, idx, mod_start, time_end, gens);
// fill values
const size_t value_end = idx + (2 * long_v::size() - (idx % 2 * long_v::size()));
......
......@@ -14,6 +14,7 @@
* limitations under the License.
*/
#include <cassert>
#include <stdexcept>
#include <tuple>
#include <storage.h>
......@@ -24,6 +25,10 @@ namespace {
// from fit to Gaussian part of ToT distribution
const float tot_mean = Constants::tot_mean;
const float tot_sigma = Constants::tot_sigma;
using namespace storage;
using std::cout;
using std::endl;
}
float GenScalar::dot_product(const std::array<float, 3>& left, const std::array<float, 3>& right) {
......@@ -41,16 +46,13 @@ std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const lo
auto& mt = gens.mt;
auto& flat = gens.flat;
// Assume times.size() is multiple of 8 here
assert(storage::n_hits % 16 == 0);
const size_t n_expect = gens.n_expect(time_end - time_start);
const float tau_l0 = gens.tau_l0();
storage_t times; times.resize(storage::n_hits);
storage_t values; values.resize(storage::n_hits);
storage_t times; times.resize(n_expect);
storage_t values; values.resize(n_expect);
times[0] = 0l;
storage_t mod_times; mod_times.resize(10000);
storage_t mod_values; mod_values.resize(10000);
size_t idx = 0;
// First generate some data
......@@ -62,16 +64,17 @@ std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const lo
long last = 0l;
while(last < time_end && idx < times.size() - 2) {
// Generate times
float r = -1.f * Constants::tau_l0 * log(flat(mt));
float r = -1.f * tau_l0 * log(flat(mt));
last += static_cast<long>(r + 0.5);
times[idx++] = last;
}
// Coincidences
idx = fill_coincidences(times, idx, mod_start, time_end, gens);
fill_coincidences(times, idx, time_start, time_end, gens);
// fill values
const size_t value_end = idx + (2 - (idx % 2));
for (size_t vidx = mod_start; vidx < value_end; vidx += 2) {
long pmt1 = flat_pmt(mt);
long pmt2 = flat_pmt(mt);
......
#include "pybind11/pybind11.h"
#include "pybind11/stl.h"
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/iostream.h>
#include "xtensor/xmath.hpp"
#include "xtensor/xarray.hpp"
......@@ -50,5 +51,8 @@ PYBIND11_MODULE(k40gen, m)
generate_k40
)pbdoc";
m.def("generate_k40", generate_k40, "Generate K40 background");
m.def("generate_k40", &generate_k40,
py::call_guard<py::scoped_ostream_redirect,
py::scoped_estream_redirect>{},
"Generate K40 background");
}
find_package(Catch2)
if(NOT Catch2_FOUND)
# catch
set(catch_VERSION "2.6.1")
set(catch_SRC_URI "https://github.com/catchorg/Catch2/releases/download/v${catch_VERSION}/catch.hpp")
set(catch_DESTDIR "${CMAKE_BINARY_DIR}/include")
set(catch_ROOTDIR "${catch_DESTDIR}/catch2")
ExternalProject_Add(catch
URL ${catch_SRC_URI}
URL_HASH SHA256=a53ef31cf0bd4c1038cdaf8a38d6c88ee762e6aba181f9534a026363bc73e430
SOURCE_DIR catch2
INSTALL_DIR ${catch_ROOTDIR}
CONFIGURE_COMMAND ""
INSTALL_COMMAND ${CMAKE_COMMAND} -E copy_directory <SOURCE_DIR> <INSTALL_DIR>
STEP_TARGETS install)
set(Catch2_FOUND TRUE)
set(builtin_catch TRUE)
else()
set(builtin_catch FALSE)
endif()
include_directories(include)
add_library(test_functions STATIC
lib/coincidences.cpp
lib/generate.cpp)
add_executable(test_k40gen
src/test_coincidences.cpp
src/test_rates.cpp
src/test_main.cpp)
target_link_libraries(test_k40gen PRIVATE generate test_functions)
add_test(TestCoincidences test_coincidences)
find_package(ROOT)
#include <array>
#include <numeric>
#include <random>
#include <iostream>
#include <generate.h>
#include <storage.h>
#include <range/v3/view/zip.hpp>
#include <range/v3/algorithm/sort.hpp>
#define CATCH_CONFIG_MAIN
#include <catch2/catch.hpp>
using namespace std;
template <std::size_t N>
struct get_n {
template <typename T>
auto operator()(T&& t) const ->
decltype(std::get<N>(std::forward<T>(t))) {
return std::get<N>(std::forward<T>(t));
}
};
pair<double, double> coincidence_rate(array<float, 4> rates) {
Generators gens{1052, 9523, rates};
long dt = std::lround(1e9);
const size_t n_expect = gens.n_expect(dt);
storage_t times; times.resize(n_expect);
storage_t values; values.resize(n_expect);
size_t idx = 0;
long time_start = 0, time_end = 0;
double av = 0.;
const int n_iter = 10000;
for (int i = 0; i < n_iter; ++i) {
time_end += dt;
auto n_coincidence = fill_coincidences(times, idx, time_start, time_end, gens);
idx = 0;
av += n_coincidence / double{n_iter};
time_start = time_end;
}
return {gens.coincidence_rate, av};
}
pair<double, double> generate_times(array<float, 4> rates, bool use_avx2) {
Generators gens{1052, 9523, rates};
long dt = std::lround(1e8);
long time_start = 0, time_end = 0;
time_end += dt;
auto [times, values] = generate(time_start, time_end, gens, use_avx2);
// Zip the ids and hits together to sort them together
auto zipped = ranges::view::zip(times, values);
// Sort the hits
ranges::sort(zipped, std::less<long>{}, get_n<0>{});
double av = 0.;
const size_t n_times = times.size();
for (size_t i = 0; i < n_times - 1; ++i) {
av += static_cast<double>(times[i]) / n_times;
}
return {av, times.size()};
}
TEST_CASE( "Rates make sense", "[rates]" ) {
auto check_rates = [](array<float, 4> rates, bool use_avx2) {
const auto [av_rate, av_n] = generate_times(rates, use_avx2);
auto coincidence_rate = std::accumulate(begin(rates), end(rates), 0.);
REQUIRE(std::abs(coincidence_rate - av_rate) / coincidence_rate < 1e-4);
cout << av_rate << " " << av_n << endl;
};
check_rates({7000.}, false);
check_rates({1000., 100.}, false);
check_rates({7000., 700., 70., 7.}, false);
}
TEST_CASE( "Coincidences make sense", "[coincidence]" ) {
auto check_coincidence = [](array<float, 4> rates) {
const auto [rate, av] = coincidence_rate(rates);
auto coincidence_rate = std::accumulate(std::next(begin(rates)), end(rates), 0.);
REQUIRE(rate == coincidence_rate);
if (rate != 0.) {
REQUIRE(std::abs(rate - av) / rate < 1e-4);
}
};
check_coincidence({7000., 700., 70., 7.});
check_coincidence({1000., 100.});
check_coincidence({7000.});
}
#pragma once
#include <array>
std::pair<double, double> generate_times(std::array<float, 4> rates, bool use_avx2);
std::pair<double, double> coincidence_rate(std::array<float, 4> rates);
#include <array>
#include <numeric>
#include <random>
#include <iostream>
#include <generate.h>
#include <storage.h>
using namespace std;
pair<double, double> coincidence_rate(array<float, 4> rates) {
Generators gens{1052, 9523, rates};
long dt = std::lround(1e9);
const size_t n_expect = gens.n_expect(dt);
storage_t times; times.resize(n_expect);
storage_t values; values.resize(n_expect);
size_t idx = 0;
long time_start = 0, time_end = 0;
double av = 0.;
const int n_iter = 10000;
for (int i = 0; i < n_iter; ++i) {
time_end += dt;
auto n_coincidence = fill_coincidences(times, idx, time_start, time_end, gens);
idx = 0;
av += n_coincidence / double{n_iter};
time_start = time_end;
}
return {gens.coincidence_rate, av};
}
#include <generate.h>
#include <storage.h>
#include <range/v3/view/zip.hpp>
#include <range/v3/algorithm/sort.hpp>
template <std::size_t N>
struct get_n {
template <typename T>
auto operator()(T&& t) const ->
decltype(std::get<N>(std::forward<T>(t))) {
return std::get<N>(std::forward<T>(t));
}
};
using namespace std;
pair<double, double> generate_times(array<float, 4> rates, bool use_avx2) {
Generators gens{1052, 9523, rates};
long dt = std::lround(1e8);
long time_start = 0, time_end = 0;
time_end += dt;
auto [times, values] = generate(time_start, time_end, gens, use_avx2);
// Zip the ids and hits together to sort them together
auto zipped = ranges::view::zip(times, values);
// Sort the hits
ranges::sort(zipped, std::less<long>{}, get_n<0>{});
double av = 0.;
const size_t n_times = times.size();
for (size_t i = 0; i < n_times - 1; ++i) {
av += static_cast<double>(times[i]) / n_times;
}
return {av, times.size()};
}
#include <numeric>
#include <catch2/catch.hpp>
#include <test_functions.h>
using namespace std;
TEST_CASE( "Coincidences make sense", "[coincidence]" ) {
auto check_coincidence = [](array<float, 4> rates) {
const auto [rate, av] = coincidence_rate(rates);
auto coincidence_rate = std::accumulate(std::next(begin(rates)), end(rates), 0.);
REQUIRE(rate == coincidence_rate);
if (rate != 0.) {
REQUIRE(std::abs(rate - av) / rate < 1e-4);
}
};
check_coincidence({7000., 700., 70., 7.});