diff --git a/tests/validation/convergence_test.hpp b/tests/validation/convergence_test.hpp new file mode 100644 index 0000000000000000000000000000000000000000..71700ea7a733c9bc8154c6981dc623ad3590831c --- /dev/null +++ b/tests/validation/convergence_test.hpp @@ -0,0 +1,120 @@ +#pragma once + +#include <util/filter.hpp> +#include <util/rangeutil.hpp> + +#include <json/json.hpp> + +#include "gtest.h" + +#include "trace_analysis.hpp" +#include "validation_data.hpp" + +namespace nest { +namespace mc { + +struct sampler_info { + const char* label; + cell_member_type probe; + simple_sampler sampler; +}; + +/* Common functionality for testing convergence towards + * a reference data set as some parameter of the model + * is changed. + * + * Type parameter Param is the type of the parameter that + * is changed between runs. + */ + +template <typename Param> +class convergence_test_runner { +private: + std::string param_name_; + bool run_validation_; + nlohmann::json meta_; + std::vector<sampler_info> cell_samplers_; + std::map<std::string, trace_data> ref_data_; + std::vector<conv_entry<Param>> conv_tbl_; + +public: + template <typename SamplerInfoSeq> + convergence_test_runner( + const std::string& param_name, + const SamplerInfoSeq& samplers, + const nlohmann::json meta + ): + param_name_(param_name), + run_validation_(false), + meta_(meta) + { + util::assign(cell_samplers_, samplers); + } + + // allow free access to JSON meta data attached to saved traces + nlohmann::json& metadata() { return meta_; } + + void load_reference_data(const util::path& ref_path) { + run_validation_ = false; + try { + ref_data_ = g_trace_io.load_traces(ref_path); + + run_validation_ = util::all_of(cell_samplers_, + [&](const sampler_info& se) { return ref_data_.count(se.label)>0; }); + + EXPECT_TRUE(run_validation_); + } + catch (std::runtime_error&) { + ADD_FAILURE() << "failure loading reference data: " << ref_path; + } + } + + template <typename Model> + void run(Model& m, Param p, float t_end, float dt) { + // reset samplers and attach to probe locations + for (auto& se: cell_samplers_) { + se.sampler.reset(); + m.attach_sampler(se.probe, se.sampler.template sampler<>()); + } + + m.run(t_end, dt); + + for (auto& se: cell_samplers_) { + std::string label = se.label; + const auto& trace = se.sampler.trace; + + // save trace + nlohmann::json trace_meta{meta_}; + trace_meta[param_name_] = p; + + g_trace_io.save_trace(label, trace, trace_meta); + + // compute metrics + if (run_validation_) { + double linf = linf_distance(trace, ref_data_[label]); + auto pd = peak_delta(trace, ref_data_[label]); + + conv_tbl_.push_back({label, p, linf, pd}); + } + } + } + + void report() { + if (run_validation_ && g_trace_io.verbose()) { + // reorder to group by id + util::stable_sort_by(conv_tbl_, [](const conv_entry<Param>& e) { return e.id; }); + report_conv_table(std::cout, conv_tbl_, param_name_); + } + } + + void assert_all_convergence() const { + for (const sampler_info& se: cell_samplers_) { + SCOPED_TRACE(se.label); + assert_convergence(util::filter(conv_tbl_, + [&](const conv_entry<Param>& e) { return e.id==se.label; })); + } + } +}; + +} // namespace mc +} // namespace nest diff --git a/tests/validation/trace_analysis.hpp b/tests/validation/trace_analysis.hpp index 36dc6a91d84e3ea78ab3bb49e599a6f0434a1369..6d9c921d97375c6d48f67dbd3b0e599b926593b1 100644 --- a/tests/validation/trace_analysis.hpp +++ b/tests/validation/trace_analysis.hpp @@ -5,6 +5,7 @@ #include "gtest.h" #include <simple_sampler.hpp> +#include <math.hpp> #include <util/optional.hpp> #include <util/path.hpp> #include <util/rangeutil.hpp> @@ -42,13 +43,15 @@ std::vector<trace_peak> local_maxima(const trace_data& u); util::optional<trace_peak> peak_delta(const trace_data& a, const trace_data& b); // Record for error data for convergence testing. +// Only linf and peak_delta are used for convergence testing below; +// if and param are for record keeping in the validation test itself. template <typename Param> struct conv_entry { std::string id; Param param; double linf; - util::optional<trace_peak> pd; + util::optional<trace_peak> peak_delta; }; template <typename Param> @@ -56,48 +59,45 @@ using conv_data = std::vector<conv_entry<Param>>; // Assert error convergence (gtest). -template <typename Param> -void assert_convergence(const conv_data<Param>& cs) { - if (cs.size()<2) return; +template <typename ConvEntrySeq> +void assert_convergence(const ConvEntrySeq& cs) { + if (util::empty(cs)) return; auto tbound = [](trace_peak p) { return std::abs(p.t)+p.t_err; }; - auto smallest_pd = cs[0].pd; + float peak_dt_bound = math::infinity<>(); - for (unsigned i = 1; i<cs.size(); ++i) { - const auto& p = cs[i-1]; - const auto& c = cs[i]; + for (auto pi = std::begin(cs); std::next(pi)!=std::end(cs); ++pi) { + const auto& p = *pi; + const auto& c = *std::next(pi); EXPECT_LE(c.linf, p.linf) << "L∞ error increase"; - EXPECT_TRUE(c.pd || (!c.pd && !p.pd)) << "divergence in peak count"; - if (c.pd && smallest_pd) { - double t = std::abs(c.pd->t); - EXPECT_LE(t, c.pd->t_err+tbound(*smallest_pd)) << "divergence in max peak displacement"; + if (!c.peak_delta) { + EXPECT_FALSE(p.peak_delta) << "divergence in peak count"; } + else { + double t = std::abs(c.peak_delta->t); + double t_limit = c.peak_delta->t_err+peak_dt_bound; + + EXPECT_LE(t, t_limit) << "divergence in max peak displacement"; - if (c.pd && (!smallest_pd || tbound(*c.pd)<tbound(*smallest_pd))) { - smallest_pd = c.pd; + peak_dt_bound = std::min(peak_dt_bound, tbound(*c.peak_delta)); } } } // Report table of convergence results. -// (Takes collection with pair<string, conv_data> -// entries.) - -template <typename Map> -void report_conv_table(std::ostream& out, const Map& tbl, const std::string& param_name) { - out << "location," << param_name << ",linf,peak_dt,peak_dt_err\n"; - for (const auto& p: tbl) { - const auto& location = p.first; - for (const auto& c: p.second) { - out << location << "," << c.param << "," << c.linf << ","; - if (c.pd) { - out << c.pd->t << "," << c.pd->t_err << "\n"; - } - else { - out << "NA,NA\n"; - } + +template <typename ConvEntrySeq> +void report_conv_table(std::ostream& out, const ConvEntrySeq& tbl, const std::string& param_name) { + out << "id," << param_name << ",linf,peak_dt,peak_dt_err\n"; + for (const auto& c: tbl) { + out << c.id << "," << c.param << "," << c.linf << ","; + if (c.peak_delta) { + out << c.peak_delta->t << "," << c.peak_delta->t_err << "\n"; + } + else { + out << "NA,NA\n"; } } } diff --git a/tests/validation/validate_ball_and_stick.cpp b/tests/validation/validate_ball_and_stick.cpp index e01eb4cf9c6f95a2631e1bbf370762967df71487..f3ca4c949033364c5250ba3c8e92f23f88ac21e5 100644 --- a/tests/validation/validate_ball_and_stick.cpp +++ b/tests/validation/validate_ball_and_stick.cpp @@ -1,70 +1,44 @@ -#include <fstream> -#include <utility> - #include <json/json.hpp> -#include <common_types.hpp> #include <cell.hpp> +#include <common_types.hpp> #include <fvm_multicell.hpp> #include <model.hpp> #include <recipe.hpp> #include <simple_sampler.hpp> -#include <util/rangeutil.hpp> -#include <util/transform.hpp> +#include <util/path.hpp> #include "gtest.h" #include "../test_common_cells.hpp" -#include "../test_util.hpp" +#include "convergence_test.hpp" #include "trace_analysis.hpp" #include "validation_data.hpp" using namespace nest::mc; -struct sampler_info { - const char* label; - cell_member_type probe; - simple_sampler sampler; -}; - template < typename lowered_cell, typename SamplerInfoSeq > void run_ncomp_convergence_test( const char* model_name, - const char* ref_data_path, + const util::path& ref_data_path, const cell& c, SamplerInfoSeq& samplers, float t_end=100.f, float dt=0.001) { - using nlohmann::json; - - SCOPED_TRACE(model_name); - - auto& V = g_trace_io; - bool verbose = V.verbose(); - int max_ncomp = V.max_ncomp(); - - auto keys = util::transform_view(samplers, - [](const sampler_info& se) { return se.label; }); - - bool run_validation = false; - std::map<std::string, trace_data> ref_data; - try { - ref_data = V.load_traces(ref_data_path); - - run_validation = std::all_of(keys.begin(), keys.end(), - [&](const char* key) { return ref_data.count(key)>0; }); - - EXPECT_TRUE(run_validation); - } - catch (std::runtime_error&) { - ADD_FAILURE() << "failure loading reference data: " << ref_data_path; - } + auto max_ncomp = g_trace_io.max_ncomp(); + nlohmann::json meta = { + {"name", "membrane voltage"}, + {"model", model_name}, + {"sim", "nestmc"}, + {"units", "mV"} + }; - std::map<std::string, std::vector<conv_entry<int>>> conv_results; + convergence_test_runner<int> R("ncomp", samplers, meta); + R.load_reference_data(ref_data_path); for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) { for (auto& seg: c.segments()) { @@ -74,48 +48,10 @@ void run_ncomp_convergence_test( } model<lowered_cell> m(singleton_recipe{c}); - // reset samplers and attach to probe locations - for (auto& se: samplers) { - se.sampler.reset(); - m.attach_sampler(se.probe, se.sampler.template sampler<>()); - } - - m.run(t_end, dt); - - for (auto& se: samplers) { - std::string key = se.label; - const simple_sampler& s = se.sampler; - - // save trace - json meta = { - {"name", "membrane voltage"}, - {"model", model_name}, - {"sim", "nestmc"}, - {"ncomp", ncomp}, - {"units", "mV"}}; - - V.save_trace(key, s.trace, meta); - - // compute metrics - if (run_validation) { - double linf = linf_distance(s.trace, ref_data[key]); - auto pd = peak_delta(s.trace, ref_data[key]); - - conv_results[key].push_back({key, ncomp, linf, pd}); - } - } - } - - if (verbose && run_validation) { - report_conv_table(std::cout, conv_results, "ncomp"); - } - - for (auto key: keys) { - SCOPED_TRACE(key); - - const auto& results = conv_results[key]; - assert_convergence(results); + R.run(m, ncomp, t_end, dt); } + R.report(); + R.assert_all_convergence(); } TEST(ball_and_stick, neuron_ref) { @@ -138,26 +74,6 @@ TEST(ball_and_stick, neuron_ref) { samplers); } -TEST(ball_and_taper, neuron_ref) { - using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>; - - cell c = make_cell_ball_and_taper(); - add_common_voltage_probes(c); - - float sample_dt = 0.025f; - sampler_info samplers[] = { - {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)}, - {"taper.mid", {0u, 1u}, simple_sampler(sample_dt)}, - {"taper.end", {0u, 2u}, simple_sampler(sample_dt)} - }; - - run_ncomp_convergence_test<lowered_cell>( - "ball_and_taper", - "neuron_ball_and_taper.json", - c, - samplers); -} - TEST(ball_and_3stick, neuron_ref) { using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>; @@ -207,3 +123,23 @@ TEST(rallpack1, numeric_ref) { 250.f); } +TEST(ball_and_taper, neuron_ref) { + using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>; + + cell c = make_cell_ball_and_taper(); + add_common_voltage_probes(c); + + float sample_dt = 0.025f; + sampler_info samplers[] = { + {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)}, + {"taper.mid", {0u, 1u}, simple_sampler(sample_dt)}, + {"taper.end", {0u, 2u}, simple_sampler(sample_dt)} + }; + + run_ncomp_convergence_test<lowered_cell>( + "ball_and_taper", + "neuron_ball_and_taper.json", + c, + samplers); +} + diff --git a/tests/validation/validate_soma.cpp b/tests/validation/validate_soma.cpp index 31eb784add73f5b675e8b2943e8faa453a0afb26..9f341a502481e4a2f4a3b94fdb51e48f5cf4c9b2 100644 --- a/tests/validation/validate_soma.cpp +++ b/tests/validation/validate_soma.cpp @@ -13,81 +13,38 @@ #include "gtest.h" -#include "../test_util.hpp" #include "../test_common_cells.hpp" +#include "convergence_test.hpp" #include "trace_analysis.hpp" #include "validation_data.hpp" using namespace nest::mc; TEST(soma, numeric_ref) { - // compare voltages against reference data produced from - // nrn/ball_and_taper.py - - using namespace nlohmann; - using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>; - auto& V = g_trace_io; - - bool verbose = V.verbose(); - - // load validation data - - bool run_validation = false; - std::map<std::string, trace_data> ref_data; - const char* key = "soma.mid"; - - const char* ref_data_path = "numeric_soma.json"; - try { - ref_data = V.load_traces(ref_data_path); - run_validation = ref_data.count(key); - - EXPECT_TRUE(run_validation); - } - catch (std::runtime_error&) { - ADD_FAILURE() << "failure loading reference data: " << ref_data_path; - } - // generate test data cell c = make_cell_soma_only(); add_common_voltage_probes(c); + model<lowered_cell> m(singleton_recipe{c}); - float sample_dt = .025; - simple_sampler sampler(sample_dt); - conv_data<float> convs; + float sample_dt = .025f; + sampler_info samplers[] = {{"soma.mid", {0u, 0u}, simple_sampler(sample_dt)}}; - for (auto dt: {0.05f, 0.02f, 0.01f, 0.005f, 0.001f}) { - sampler.reset(); - model<lowered_cell> m(singleton_recipe{c}); - - m.attach_sampler({0u, 0u}, sampler.sampler<>()); - m.run(100, dt); - - // save trace - auto& trace = sampler.trace; - json meta = { - {"name", "membrane voltage"}, - {"model", "soma"}, - {"sim", "nestmc"}, - {"dt", dt}, - {"units", "mV"}}; - - V.save_trace(key, trace, meta); + nlohmann::json meta = { + {"name", "membrane voltage"}, + {"model", "soma"}, + {"sim", "nestmc"}, + {"units", "mV"} + }; - // compute metrics - if (run_validation) { - double linf = linf_distance(trace, ref_data[key]); - auto pd = peak_delta(trace, ref_data[key]); + convergence_test_runner<float> R("dt", samplers, meta); + R.load_reference_data("numeric_soma.json"); - convs.push_back({key, dt, linf, pd}); - } - } - - if (verbose && run_validation) { - std::map<std::string, std::vector<conv_entry<float>>> conv_results = {{key, convs}}; - report_conv_table(std::cout, conv_results, "dt"); + float t_end = 100.f; + for (auto dt: {0.05f, 0.02f, 0.01f, 0.005f, 0.001f}) { + m.reset(); + R.run(m, dt, t_end, dt); } - - SCOPED_TRACE("soma.mid"); - assert_convergence(convs); + R.report(); + R.assert_all_convergence(); } diff --git a/tests/validation/validate_synapses.cpp b/tests/validation/validate_synapses.cpp index 65c99ba345cc6ae7ed52b0703974386c2929bba5..9460d0c3cad2506cc8e6b7104f8cde262d3935c9 100644 --- a/tests/validation/validate_synapses.cpp +++ b/tests/validation/validate_synapses.cpp @@ -1,45 +1,37 @@ -#include <fstream> -#include <utility> - #include <json/json.hpp> -#include <common_types.hpp> #include <cell.hpp> #include <fvm_multicell.hpp> #include <model.hpp> #include <recipe.hpp> #include <simple_sampler.hpp> -#include <util/rangeutil.hpp> -#include <util/transform.hpp> +#include <util/path.hpp> #include "gtest.h" #include "../test_common_cells.hpp" -#include "../test_util.hpp" +#include "convergence_test.hpp" #include "trace_analysis.hpp" #include "validation_data.hpp" using namespace nest::mc; -void run_synapse_test(const char* syn_type, const char* ref_file) { - using namespace nlohmann; - +void run_synapse_test( + const char* syn_type, + const util::path& ref_data_path, + float t_end=70.f, + float dt=0.001) +{ using lowered_cell = fvm::fvm_multicell<double, cell_local_size_type>; - auto& V = g_trace_io; - - bool verbose = V.verbose(); - int max_ncomp = V.max_ncomp(); - // load validation data - auto ref_data = V.load_traces(ref_file); - bool run_validation = - ref_data.count("soma.mid") && - ref_data.count("dend.mid") && - ref_data.count("dend.end"); - - EXPECT_TRUE(run_validation); + auto max_ncomp = g_trace_io.max_ncomp(); + nlohmann::json meta = { + {"name", "membrane voltage"}, + {"model", syn_type}, + {"sim", "nestmc"}, + {"units", "mV"} + }; - // generate test data cell c = make_cell_ball_and_stick(false); // no stimuli parameter_list syn_default(syn_type); c.add_synapse({1, 0.5}, syn_default); @@ -52,66 +44,25 @@ void run_synapse_test(const char* syn_type, const char* ref_file) { {{0u, 0u}, 40.0, 0.04} }; - float sample_dt = .025; - std::pair<const char *, simple_sampler> samplers[] = { - {"soma.mid", simple_sampler(sample_dt)}, - {"dend.mid", simple_sampler(sample_dt)}, - {"dend.end", simple_sampler(sample_dt)} + float sample_dt = 0.025f; + sampler_info samplers[] = { + {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)}, + {"dend.mid", {0u, 1u}, simple_sampler(sample_dt)}, + {"dend.end", {0u, 2u}, simple_sampler(sample_dt)} }; - std::map<std::string, std::vector<conv_entry<int>>> conv_results; + convergence_test_runner<int> R("ncomp", samplers, meta); + R.load_reference_data(ref_data_path); for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) { - for (auto& se: samplers) { - se.second.reset(); - } c.cable(1)->set_compartments(ncomp); model<lowered_cell> m(singleton_recipe{c}); - - // the ball-and-stick-cell (should) have three voltage probes: - // centre of soma, centre of dendrite, end of dendrite. - - m.attach_sampler({0u, 0u}, samplers[0].second.sampler<>()); - m.attach_sampler({0u, 1u}, samplers[1].second.sampler<>()); - m.attach_sampler({0u, 2u}, samplers[2].second.sampler<>()); - m.group(0).enqueue_events(synthetic_events); - m.run(70, 0.001); - - for (auto& se: samplers) { - std::string key = se.first; - const simple_sampler& s = se.second; - - // save trace - json meta = { - {"name", "membrane voltage"}, - {"model", syn_type}, - {"sim", "nestmc"}, - {"ncomp", ncomp}, - {"units", "mV"}}; - - V.save_trace(key, s.trace, meta); - - // compute metrics - if (run_validation) { - double linf = linf_distance(s.trace, ref_data[key]); - auto pd = peak_delta(s.trace, ref_data[key]); - - conv_results[key].push_back({key, ncomp, linf, pd}); - } - } - } - - if (verbose && run_validation) { - report_conv_table(std::cout, conv_results, "ncomp"); - } - - for (auto key: util::transform_view(samplers, util::first)) { - SCOPED_TRACE(key); - const auto& results = conv_results[key]; - assert_convergence(results); + R.run(m, ncomp, t_end, dt); } + R.report(); + R.assert_all_convergence(); } TEST(simple_synapse, expsyn_neuron_ref) {