From dff0a86f80a4f660054e23c2afb05f1b9611d780 Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Wed, 29 Sep 2021 10:07:17 +0200
Subject: [PATCH] Segfault on instantiating mechanisms on empty regions (#1657)

More, subtler fallout of the ABI patchset: When mechanisms are added to empty regions, various, layered segfaults
appear.
1. `set_parameter` will receive a `nullptr` upon linear search, throwing `no such parameter`
2. Guarding this for finite `width` uncovers further need in the ABI interface methods
3. This is caused by early exit in `instantiate` on `width == 0`
4. Getting rid of this is better, as we can treat vacuous mechanism painting like all others...
5. ...but causes another segfault, since we default index arrays to `cv_pos.back()` but `cv_pos` is empty.

This patch fixes all of the above and in addition removes mechanisms without support before reifying them
on the FVM cells. Relevant tests are added.
---
 arbor/backends/multicore/shared_state.cpp |  15 ++--
 arbor/fvm_layout.cpp                      |  10 +--
 test/unit/test_abi.cpp                    | 100 +++++++++++++++++++++-
 test/unit/test_fvm_lowered.cpp            |  27 ++++++
 4 files changed, 139 insertions(+), 13 deletions(-)

diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp
index 6936f56d..34a9853d 100644
--- a/arbor/backends/multicore/shared_state.cpp
+++ b/arbor/backends/multicore/shared_state.cpp
@@ -430,17 +430,20 @@ std::size_t extend_width(const arb::mechanism& mech, std::size_t width) {
 void shared_state::set_parameter(mechanism& m, const std::string& key, const std::vector<arb_value_type>& values) {
     if (values.size()!=m.ppack_.width) throw arbor_internal_error("mechanism field size mismatch");
 
+    bool found = false;
     arb_value_type* data = nullptr;
     for (arb_size_type i = 0; i<m.mech_.n_parameters; ++i) {
         if (key==m.mech_.parameters[i].name) {
             data = m.ppack_.parameters[i];
+            found = true;
             break;
         }
     }
 
-    if (!data) throw arbor_internal_error(util::pprintf("no such parameter '{}'", key));
+    if (!found) throw arbor_internal_error(util::pprintf("no such parameter '{}'", key));
 
     if (!m.ppack_.width) return;
+
     auto width_padded = extend_width<arb_value_type>(m, m.ppack_.width);
     copy_extend(values, util::range_n(data, width_padded), values.back());
 }
@@ -519,9 +522,6 @@ void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_o
         m.ppack_.ion_states[idx] = { oion->iX_.data(), oion->eX_.data(), oion->Xi_.data(), oion->Xo_.data(), oion->charge.data() };
     }
 
-    // If there are no sites (is this ever meaningful?) there is nothing more to do.
-    if (m.ppack_.width==0) return;
-
     // Initialize state and parameter vectors with default values.
     {
         // Allocate bulk storage
@@ -567,8 +567,11 @@ void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_o
         store.indices_ = iarray(count*index_width_padded, 0, pad);
         chunk_writer writer(store.indices_.data(), index_width_padded);
         // Setup node indices
-        m.ppack_.node_index = writer.append(pos_data.cv, pos_data.cv.back());
-
+        //   We usually insert cv.size() == width elements into node index (length: width_padded >= width)
+        //   and pad by the last element of cv. If width == 0 we must choose a different pad, that will not
+        //   really be used, as width == width_padded == 0. Nevertheless, we need to pass it.
+        auto pad_val = pos_data.cv.empty() ? 0 : pos_data.cv.back();
+        m.ppack_.node_index = writer.append(pos_data.cv, pad_val);
         auto node_index = util::range_n(m.ppack_.node_index, index_width_padded);
         // Make SIMD index constraints and set the view
         store.constraints_ = make_constraint_partition(node_index, m.ppack_.width, m.iface_.partition_width);
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index d2ee99ce..12ed2a34 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -943,7 +943,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
         }
 
         update_ion_support(info, config.cv);
-        M.mechanisms[name] = std::move(config);
+        if (!config.cv.empty()) M.mechanisms[name] = std::move(config);
     }
 
     // Synapses:
@@ -1085,7 +1085,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
         update_ion_support(info, config.cv);
 
         M.n_target += config.target.size();
-        M.mechanisms[name] = std::move(config);
+        if (!config.cv.empty()) M.mechanisms[name] = std::move(config);
     }
     M.post_events = post_events;
 
@@ -1136,7 +1136,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
         std::unique_copy(config.cv.begin(), config.cv.end(), std::back_inserter(config.cv_unique));
         config.cv_unique.shrink_to_fit();
 
-        M.stimuli = std::move(config);
+        if (!config.cv.empty()) M.stimuli = std::move(config);
     }
 
     // Ions:
@@ -1205,7 +1205,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
             config.init_econc[i] *= oo_cv_area;
         }
 
-        M.ions[ion] = std::move(config);
+        if (!config.cv.empty()) M.ions[ion] = std::move(config);
     }
 
     std::unordered_map<std::string, mechanism_desc> revpot_tbl;
@@ -1276,7 +1276,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
                         config.param_values.emplace_back(kv.first, std::vector<value_type>(config.cv.size(), kv.second));
                     }
 
-                    M.mechanisms[revpot.name()] = std::move(config);
+                    if (!config.cv.empty()) M.mechanisms[revpot.name()] = std::move(config);
                 }
             }
         }
diff --git a/test/unit/test_abi.cpp b/test/unit/test_abi.cpp
index 5db2510b..5cc736e1 100644
--- a/test/unit/test_abi.cpp
+++ b/test/unit/test_abi.cpp
@@ -59,7 +59,7 @@ TEST(abi, multicore_initialisation) {
     layout.weight.assign(ncv, 1.);
     for (arb_size_type i = 0; i<ncv; ++i) layout.cv.push_back(i);
 
-    shared_state.instantiate(mech, 42, {}, layout);
+    EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout));
 
     {
         ASSERT_EQ(globals.size(), mech.mech_.n_globals);
@@ -93,6 +93,53 @@ TEST(abi, multicore_initialisation) {
     }
 }
 
+TEST(abi, multicore_null) {
+    std::vector<arb_field_info> globals = {{ "G0", "kg",  123.0,     0.0, 2000.0},
+                                           { "G1", "lb",  456.0,     0.0, 2000.0},
+                                           { "G2", "gr",  789.0,     0.0, 2000.0}};
+    std::vector<arb_field_info> states  = {{ "S0", "nA",      0.123, 0.0, 2000.0},
+                                           { "S1", "mV",      0.456, 0.0, 2000.0}};
+    std::vector<arb_field_info> params  = {{ "P0", "lm", -123.0,     0.0, 2000.0}};
+
+    arb_mechanism_type type{};
+    type.abi_version = ARB_MECH_ABI_VERSION;
+    type.globals    = globals.data(); type.n_globals    = globals.size();
+    type.parameters = params.data();  type.n_parameters = params.size();
+    type.state_vars = states.data();  type.n_state_vars = states.size();
+
+    arb_mechanism_interface iface { arb_backend_kind_cpu,
+                                    1,
+                                    1,
+                                    nullptr,
+                                    nullptr,
+                                    nullptr,
+                                    nullptr,
+                                    nullptr,
+                                    nullptr };
+
+    auto mech = arb::mechanism(type, iface);
+
+    arb_size_type ncell = 1;
+    arb_size_type ncv = 0;
+    std::vector<arb_index_type> cv_to_intdom(ncv, 0);
+    std::vector<arb_value_type> temp(ncv, 23);
+    std::vector<arb_value_type> diam(ncv, 1.);
+    std::vector<arb_value_type> vinit(ncv, -65);
+    std::vector<arb::fvm_gap_junction> gj = {};
+    std::vector<arb_index_type> src_to_spike = {};
+
+    arb::multicore::shared_state shared_state(ncell, ncell, 0,
+                                              cv_to_intdom, cv_to_intdom,
+                                              gj, vinit, temp, diam, src_to_spike,
+                                              mech.data_alignment());
+
+    arb::mechanism_layout layout;
+    layout.weight.assign(ncv, 1.);
+    for (arb_size_type i = 0; i<ncv; ++i) layout.cv.push_back(i);
+
+    EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout));
+}
+
 #ifdef ARB_GPU_ENABLED
 
 namespace {
@@ -155,7 +202,7 @@ TEST(abi, gpu_initialisation) {
     layout.weight.assign(ncv, 1.);
     for (arb_size_type i = 0; i<ncv; ++i) layout.cv.push_back(i);
 
-    shared_state.instantiate(mech, 42, {}, layout);
+    EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout));
 
     {
         ASSERT_EQ(globals.size(), mech.mech_.n_globals);
@@ -187,4 +234,53 @@ TEST(abi, gpu_initialisation) {
         }
     }
 }
+
+TEST(abi, gpu_null) {
+    std::vector<arb_field_info> globals = {{ "G0", "kg",  123.0,     0.0, 2000.0},
+                                           { "G1", "lb",  456.0,     0.0, 2000.0},
+                                           { "G2", "gr",  789.0,     0.0, 2000.0}};
+    std::vector<arb_field_info> states  = {{ "S0", "nA",      0.123, 0.0, 2000.0},
+                                           { "S1", "mV",      0.456, 0.0, 2000.0}};
+    std::vector<arb_field_info> params  = {{ "P0", "lm", -123.0,     0.0, 2000.0}};
+
+    arb_mechanism_type type{};
+    type.abi_version = ARB_MECH_ABI_VERSION;
+    type.globals    = globals.data(); type.n_globals    = globals.size();
+    type.parameters = params.data();  type.n_parameters = params.size();
+    type.state_vars = states.data();  type.n_state_vars = states.size();
+
+    arb_mechanism_interface iface { arb_backend_kind_gpu,
+                                    1,
+                                    1,
+                                    nullptr,
+                                    nullptr,
+                                    nullptr,
+                                    nullptr,
+                                    nullptr,
+                                    nullptr };
+
+    auto mech = arb::mechanism(type, iface);
+
+    arb_size_type ncell = 1;
+    arb_size_type ncv = 0;
+    std::vector<arb_index_type> cv_to_intdom(ncv, 0);
+    std::vector<arb_value_type> temp(ncv, 23);
+    std::vector<arb_value_type> diam(ncv, 1.);
+    std::vector<arb_value_type> vinit(ncv, -65);
+    std::vector<arb::fvm_gap_junction> gj = {};
+    std::vector<arb_index_type> src_to_spike = {};
+
+    arb::gpu::shared_state shared_state(ncell, ncell, 0,
+                                        cv_to_intdom, cv_to_intdom,
+                                        gj, vinit, temp, diam, src_to_spike,
+                                        1);
+
+    arb::mechanism_layout layout;
+    layout.weight.assign(ncv, 1.);
+    for (arb_size_type i = 0; i<ncv; ++i) layout.cv.push_back(i);
+
+    EXPECT_NO_THROW(shared_state.instantiate(mech, 42, {}, layout));
+}
+
+
 #endif
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index fcaa3dd9..d8295388 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -535,6 +535,33 @@ TEST(fvm_lowered, derived_mechs) {
     }
 }
 
+TEST(fvm_lowered, null_region) {
+    arb::proc_allocation resources;
+    if (auto nt = arbenv::get_env_num_threads()) {
+        resources.num_threads = nt;
+    }
+    else {
+        resources.num_threads = arbenv::thread_concurrency();
+    }
+
+    soma_cell_builder builder(6);
+    builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
+    auto cell = builder.make_cell();
+
+    cell.decorations.paint(reg::nil(), "test_kin1");
+    cell.decorations.paint(reg::nil(), "custom_kin1");
+
+    cable1d_recipe rec(cable_cell{cell});
+    rec.catalogue() = make_unit_test_catalogue();
+    rec.catalogue().derive("custom_kin1", "test_kin1", {{"tau", 20.0}});
+
+    auto ctx = make_context(resources);
+    auto decomp = partition_load_balance(rec, ctx);
+    simulation sim(rec, decomp, ctx);
+    EXPECT_NO_THROW(sim.run(30.0, 1.f/1024));
+}
+
+
 // Test that ion charge is propagated into mechanism variable.
 
 TEST(fvm_lowered, read_valence) {
-- 
GitLab