Commit 97b7f0ea authored by Roel Aaij's avatar Roel Aaij
Browse files

Add a numbering scheme for the dom id, reference or orca

parent 69a28310
......@@ -151,11 +151,11 @@ endif()
if(ENABLE_PYTHON)
include(XTensor)
# pybind11
find_package(pybind11 REQUIRED)
include(pybind11Tools)
pybind11_add_module(k40gen "src/k40gen/bindings.cpp")
target_link_libraries(k40gen PUBLIC generate)
find_package(pybind11)
target_include_directories(k40gen PUBLIC ${PYBIND11_INCLUDE_DIRS})
if (xtp_BUILTIN)
......@@ -168,10 +168,6 @@ if(ENABLE_PYTHON)
target_compile_definitions(k40gen PRIVATE "-DXTENSOR_ENABLE_XSIMD")
endif()
# pybind11 headers
find_package(pybind11 REQUIRED)
target_include_directories(k40gen PRIVATE ${pybind11_INCLUDE_DIRS})
# numpy headers
find_package(NumPy REQUIRED)
target_include_directories(k40gen PRIVATE ${NUMPY_INCLUDE_DIRS})
......
......@@ -60,4 +60,4 @@ float dot_product(const Vc::SimdArray<float, 3>& left, const Vc::SimdArray<float
}
std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long time_end,
Generators& gens);
Generators& gens, dom_fun_t dom_fun);
......@@ -15,12 +15,20 @@
*/
#pragma once
#include <map>
#include <vector>
#include <numeric>
#include <functional>
#include <array>
#include <random>
#include "storage.h"
unsigned int ref_dom(int dom, int mod);
unsigned int orca_dom(int dom, int mod);
using dom_fun_t = std::function<unsigned int(int, int)>;
namespace Constants {
// 2 * pi
......@@ -41,6 +49,10 @@ namespace Constants {
const float tot_sigma = 2.44078f;
const int random_method = 2;
static const std::map<std::string, std::function<unsigned int(int, int)>> dom_funs
{{"reference", ref_dom},
{"orca", orca_dom}};
}
float cross_prob(const float ct);
......@@ -121,4 +133,5 @@ fill_coincidences(storage_t& times, pmts_t& pmts, size_t idx,
Generators& gens);
std::tuple<storage_t, storage_t> generate(const long start, const long end,
Generators& gens, bool use_avx2);
Generators& gens, std::string scheme,
bool use_avx2);
......@@ -58,4 +58,4 @@ float dot_product(const std::array<float, 3>& left, const std::array<float, 3>&
}
std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const long time_end,
Generators& gens);
Generators& gens, dom_fun_t dom_fun);
......@@ -20,6 +20,7 @@
#include <stdexcept>
#include <vector>
#include <tuple>
#include <sstream>
#include <generate.h>
#include <storage.h>
......@@ -32,6 +33,14 @@ namespace {
using std::array;
}
unsigned int ref_dom(int dom, int mod) {
return 100 * (dom + 1) + mod + 1;
}
unsigned int orca_dom(int dom, int mod) {
return 18 * dom + mod + 1;
}
Generators::Generators(const int seed0, const int seed1, const std::array<float, 4> rates)
: m_seeds{seed0, seed1},
m_rates{std::move(rates)},
......@@ -128,15 +137,27 @@ fill_coincidences(storage_t& times, pmts_t& pmts, size_t idx,
return {pmts.size(), n};
}
std::tuple<storage_t, storage_t> generate(const long start, const long end,
Generators& gens, bool use_avx2) {
std::tuple<storage_t, storage_t> generate(const long start, const long stop,
Generators& gens, std::string dom_fun,
bool use_avx2) {
auto it = Constants::dom_funs.find(dom_fun);
std::stringstream funs;
std::for_each(begin(Constants::dom_funs), end(Constants::dom_funs),
[&funs] (const auto& entry) {
funs << " " << entry.first;
});
if (it == end(Constants::dom_funs)) {
throw std::runtime_error{std::string{"wrong dom mapping function, should be one of: "} +
funs.str()};
}
#ifdef USE_AVX2
if (use_avx2) {
return generate_avx2(start, end, gens);
return generate_avx2(start, stop, gens, it->second);
} else {
return generate_scalar(start, end, gens);
return generate_scalar(start, stop, gens, it->second);
}
#else
return generate_scalar(start, end, gens);
return generate_scalar(start, stop, gens, it->second);
#endif
}
......@@ -66,8 +66,7 @@ auto int_v_to_long_v(const int_v& in) -> pair<long_v, long_v>
};
void fill_values_avx2(long idx_start, long idx_end, storage_t& values, Ranvec1& random,
int dom, int mod,
std::function<int_v(size_t)> pmt_fun) {
int dom_id, std::function<int_v(size_t)> pmt_fun) {
// fill values
size_t n = 0;
const long value_end = idx_end + (2 * long_v::size() - (idx_end % 2 * long_v::size()));
......@@ -85,13 +84,13 @@ void fill_values_avx2(long idx_start, long idx_end, storage_t& values, Ranvec1&
z0 = fact * tot_sigma * z0 + tot_mean;
z1 = fact * tot_sigma * z1 + tot_mean;
auto val0 = simd_cast<int_v>(z0) | (pmt1 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
auto val0 = simd_cast<int_v>(z0) | (pmt1 << 8) | dom_id << 13;
auto [val0_first, val0_second] = int_v_to_long_v(val0);
val0_first.store(&values[vidx]);
val0_second.store(&values[vidx + long_v::size()]);
auto val1 = simd_cast<int_v>(z1) | (pmt2 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
auto val1 = simd_cast<int_v>(z1) | (pmt2 << 8) | dom_id << 13;
auto [val1_first, val1_second] = int_v_to_long_v(val1);
val1_first.store(&values[vidx + 2 * long_v::size()]);
val1_second.store(&values[vidx + 3 *long_v::size()]);
......@@ -99,7 +98,7 @@ void fill_values_avx2(long idx_start, long idx_end, storage_t& values, Ranvec1&
}
std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long time_end,
Generators& gens) {
Generators& gens, dom_fun_t dom_fun) {
Ranvec1 random{Constants::random_method};
random.init(gens.seeds()[0], gens.seeds()[1]);
......@@ -121,6 +120,7 @@ std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long
// First generate some data
for (int dom = 0; dom < Constants::n_dom; ++dom) {
for (int mod = 0; mod < Constants::n_mod; ++mod) {
auto dom_id = dom_fun(dom, mod);
for (int pmt = 0; pmt < Constants::n_pmt; ++pmt) {
size_t pmt_start = idx;
long_v offset;
......@@ -157,14 +157,14 @@ std::tuple<storage_t, storage_t> generate_avx2(const long time_start, const long
idx += 2 * long_v::size();
}
fill_values_avx2(pmt_start, idx, values, random, dom, mod,
fill_values_avx2(pmt_start, idx, values, random, dom_id,
[pmt](size_t) { return int_v(pmt); });
}
// Coincidences
auto [n_times, _] = fill_coincidences(times, pmts, idx, time_start, time_end, gens);
fill_values_avx2(idx, idx + n_times, values, random, dom, mod,
fill_values_avx2(idx, idx + n_times, values, random, dom_id,
[&pmts](size_t n) {
return int_v(pmts.data() + n * int_v::size());
});
......
......@@ -42,7 +42,8 @@ float GenScalar::dot_product(const std::array<float, 3>& left, const std::array<
return r;
}
void fill_values_scalar(long idx_start, long idx_end, storage_t& values, Generators& gens, int dom, int mod,
void fill_values_scalar(long idx_start, long idx_end, storage_t& values, Generators& gens,
unsigned int dom_id,
std::function<unsigned int(size_t)> pmt_fun) {
// fill values
auto& mt = gens.mt;
......@@ -63,16 +64,17 @@ void fill_values_scalar(long idx_start, long idx_end, storage_t& values, Generat
z0 = fact * tot_sigma * z0 + tot_mean;
z1 = fact * tot_sigma * z1 + tot_mean;
auto val0 = static_cast<int>(z0) | (pmt1 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
auto val0 = static_cast<int>(z0) | (pmt1 << 8) | dom_id << 13;
values[vidx] = val0;
auto val1 = static_cast<int>(z1) | (pmt2 << 8) | ((100 * (dom + 1) + mod + 1) << 13);
auto val1 = static_cast<int>(z1) | (pmt2 << 8) | dom_id << 13;
values[vidx + 1] = val1;
}
}
std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const long time_end,
Generators& gens) {
Generators& gens,
dom_fun_t dom_fun) {
auto& mt = gens.mt;
auto& flat = gens.flat;
......@@ -89,6 +91,7 @@ std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const lo
// First generate some data
for (int dom = 0; dom < Constants::n_dom; ++dom) {
for (int mod = 0; mod < Constants::n_mod; ++mod) {
auto dom_id = dom_fun(dom, mod);
for (int pmt = 0; pmt < Constants::n_pmt; ++pmt) {
size_t pmt_start = idx;
long last = time_start;
......@@ -98,13 +101,13 @@ std::tuple<storage_t, storage_t> generate_scalar(const long time_start, const lo
last += static_cast<long>(r + 0.5);
times[idx++] = last;
}
fill_values_scalar(pmt_start, idx, values, gens, dom, mod,
fill_values_scalar(pmt_start, idx, values, gens, dom_id,
[pmt](size_t) { return pmt; });
}
// Coincidences
auto [n_times, _] = fill_coincidences(times, pmts, idx, time_start, time_end, gens);
fill_values_scalar(idx, idx + n_times, values, gens, dom, mod,
fill_values_scalar(idx, idx + n_times, values, gens, dom_id,
[&pmts](size_t n) {
assert(n < pmts.size());
return pmts[n];
......
......@@ -24,9 +24,10 @@ namespace {
}
inline xt::pyarray<long>
generate_k40(const long time_start, const long time_end, Generators& gens, bool avx2)
generate_k40(const long time_start, const long time_end, Generators& gens,
std::string scheme, bool avx2)
{
auto r = generate(time_start, time_end, gens, avx2);
auto r = generate(time_start, time_end, gens, scheme, avx2);
auto times = std::get<0>(r);
auto values = std::get<1>(r);
auto n = times.size();
......
......@@ -60,9 +60,8 @@ if (ROOT_FOUND)
add_test(TestK40ROOT test_k40gen_root)
endif()
find_package(Python3 COMPONENTS Interpreter)
if(Python3_Interpreter_FOUND)
if(ENABLE_PYTHON)
add_test(NAME pytest
COMMAND Python3::Interpreter -m pytest -rsv ${PROJECT_SOURCE_DIR}/tests/python/test_extension.py
COMMAND ${PYTHON_EXECUTABLE} -m pytest -rsv ${PROJECT_SOURCE_DIR}/tests/python/test_extension.py
WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
endif()
......@@ -13,7 +13,7 @@ pair<double, double> generate_l0(float l0_rate, long dt, bool use_avx2) {
long time_start = 0, time_end = 0;
time_end += dt;
auto [times, values] = generate(time_start, time_end, gens, use_avx2);
auto [times, values] = generate(time_start, time_end, gens, "reference", use_avx2);
double av = 0.;
const size_t n_times = times.size();
......
def generate():
import k40gen as m
gens = m.Generators(21341, 1245, [7000., 700., 70., 0.])
times = m.generate_k40(0, int(1e8), gens, False)
times = m.generate_k40(0, int(1e8), gens, False, "reference")
assert(times.shape[0] == 4)
......
......@@ -43,7 +43,7 @@ TEST_CASE( "Rates makes sense [ROOT]", "[rates_ROOT]" ) {
long time_start = 0, time_end = time_start + dt;
auto [times, values] = generate(time_start, time_end, gens, use_avx2);
auto [times, values] = generate(time_start, time_end, gens, "reference", use_avx2);
const size_t n_times = times.size();
for (size_t i = 0; i < n_times - 1; ++i) {
if (((values[i + 1]) >> 8) == (values[i] >> 8)) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment