Mutable Parameters

CMakeLists.txt

@@ -357,8 +357,7 @@ if (BUILD_TESTING)
                   OPTIONS "BENCHMARK_ENABLE_TESTING OFF" "CMAKE_BUILD_TYPE release" "BUILD_SHARED_LIBS OFF")
     CPMFindPackage(NAME googletest
                   GITHUB_REPOSITORY google/googletest
-                  GIT_TAG release-1.12.1
-                  VERSION 1.12.1
+                  VERSION 1.17.0
                   OPTIONS "INSTALL_GTEST OFF" "BUILD_GMOCK OFF")
     if(benchmark_ADDED)
         target_link_libraries(ext-bench INTERFACE benchmark)

arbor/backends/multicore/shared_state.cpp

@@ -23,7 +23,6 @@
 
 #include "multicore_common.hpp"
 #include "shared_state.hpp"
-#include "fvm.hpp"
 
 namespace arb {
 namespace multicore {
@@ -183,7 +182,6 @@ void istim_state::add_current(const arb_value_type time, array& current_density)
 }
 
 // shared_state methods:
-
 shared_state::shared_state(task_system_handle,    // ignored in mc backend
                            arb_size_type n_cell,
                            arb_size_type n_cv_,
@@ -400,10 +398,8 @@ void shared_state::instantiate(arb::mechanism& m,
     bool peer_indices = !pos_data.peer_cv.empty();
 
     // store indices for random number generation
-    if (m.mech_.n_random_variables) {
-        store.gid_ = pos_data.gid;
-        store.idx_ = pos_data.idx;
-    }
+    store.gid_ = pos_data.gid;
+    if (m.mech_.n_random_variables) store.idx_ = pos_data.idx;
 
     // Allocate view pointers (except globals!)
     store.state_vars_.resize(m.mech_.n_state_vars); m.ppack_.state_vars = store.state_vars_.data();
@@ -533,5 +529,14 @@ void shared_state::instantiate(arb::mechanism& m,
     }
 }
 
+void shared_state::update_density_data(cell_gid_type lid, arb_mechanism_ppack& ppack, cell_gid_type pid, arb_value_type val) {
+    for (auto idx = 0ul; idx < ppack.width; ++idx) {
+        if (lid != ppack.vec_ci[ppack.node_index[idx]]) continue;
+        std::cerr << "old=" << ppack.parameters[pid][idx];
+        ppack.parameters[pid][idx] = val;
+        std::cerr << " new=" << ppack.parameters[pid][idx] << '\n';
+    }
+}
+
 } // namespace multicore
 } // namespace arb

arbor/backends/multicore/shared_state.hpp

@@ -239,6 +239,8 @@ struct ARB_ARBOR_API shared_state:
         sample_time_host = util::range_pointer_view(sample_time);
         sample_value_host = util::range_pointer_view(sample_value);
     }
+
+    void update_density_data(cell_gid_type gid, arb_mechanism_ppack& ppack, cell_gid_type pid, arb_value_type val);
 };
 
 // For debugging only:

arbor/backends/shared_state_base.hpp

@@ -42,8 +42,8 @@ struct shared_state_base {
         // samples
         auto n_samples = util::sum_by(samples, [] (const auto& s) {return s.size();});
         if (d->sample_time.size() < n_samples) {
-            d->sample_time = array(n_samples);
-            d->sample_value = array(n_samples);
+            d->sample_time.resize(n_samples);
+            d->sample_value.resize(n_samples);
         }
         initialize(samples, d->sample_events);
         // thresholds
@@ -80,6 +80,11 @@ struct shared_state_base {
         }
     }
 
+    void update_density_data(cell_gid_type gid, arb_mechanism_ppack& ppack, cell_gid_type pid, arb_value_type val) {
+        auto d = static_cast<D*>(this);
+        d->update_density_data(gid, ppack, pid, val);
+    }
+
     arb_value_type* mechanism_state_data(const mechanism& m,
                                          const std::string& key) {
         auto d = static_cast<D*>(this);

arbor/benchmark_cell_group.cpp

@@ -11,6 +11,7 @@
 #include "profile/profiler_macro.hpp"
 
 #include "util/span.hpp"
+#include "util/maputil.hpp"
 
 template<typename K>
 void serialize(arb::serializer& s, const K& k, const arb::benchmark_cell_group&);
@@ -23,8 +24,7 @@ benchmark_cell_group::benchmark_cell_group(const std::vector<cell_gid_type>& gid
                                            const recipe& rec,
                                            cell_label_range& cg_sources,
                                            cell_label_range& cg_targets):
-    gids_(gids)
-{
+    gids_(gids) {
     for (auto gid: gids_) {
         if (!rec.get_probes(gid).empty()) {
             throw bad_cell_probe(cell_kind::benchmark, gid);
@@ -47,28 +47,39 @@ benchmark_cell_group::benchmark_cell_group(const std::vector<cell_gid_type>& gid
 }
 
 void benchmark_cell_group::reset() {
-    for (auto& c: cells_) {
-        c.time_sequence.reset();
-    }
-
+    for (auto& c: cells_) c.time_sequence.reset();
     clear_spikes();
 }
 
-void benchmark_cell_group::t_serialize(serializer& ser, const std::string& k) const {
-    serialize(ser, k, *this);
-}
-void benchmark_cell_group::t_deserialize(serializer& ser, const std::string& k) {
-    deserialize(ser, k, *this);
+void
+benchmark_cell_group::edit_cell(cell_gid_type gid, std::any cell_edit) {
+    try {
+        auto bench_edit = std::any_cast<benchmark_cell_editor>(cell_edit);
+        auto lid = util::binary_search_index(gids_, gid);
+        if (!lid) throw arb::arbor_internal_error{"gid " + std::to_string(gid) + " erroneuosly dispatched to cell group."};
+        benchmark_cell& lowered = cells_[*lid];
+        auto tmp = benchmark_cell{.source=lowered.source, .target=lowered.target, .time_sequence=std::move(lowered.time_sequence), .realtime_ratio=lowered.realtime_ratio};
+        bench_edit(tmp);
+        if (tmp.source != lowered.source) throw bad_cell_edit(gid, "Source is not editable.");
+        if (tmp.target != lowered.target) throw bad_cell_edit(gid, "Target is not editable.");
+        // Write back
+        lowered.time_sequence = std::move(tmp.time_sequence);
+        lowered.realtime_ratio = tmp.realtime_ratio;
+    }
+    catch (const std::bad_any_cast&) {
+        throw bad_cell_edit(gid, "Not a Benchmark editor (C++ type-id: '" + std::string{cell_edit.type().name()} + "')");
+    }
 }
 
-cell_kind benchmark_cell_group::get_cell_kind() const {
-    return cell_kind::benchmark;
-}
+void benchmark_cell_group::t_serialize(serializer& ser, const std::string& k) const { serialize(ser, k, *this); }
+
+void benchmark_cell_group::t_deserialize(serializer& ser, const std::string& k) { deserialize(ser, k, *this); }
+
+cell_kind benchmark_cell_group::get_cell_kind() const { return cell_kind::benchmark; }
 
 void benchmark_cell_group::advance(epoch ep,
                                    time_type dt,
-                                   const event_lane_subrange& event_lanes)
-{
+                                   const event_lane_subrange& event_lanes) {
     using std::chrono::high_resolution_clock;
     using duration_type = std::chrono::duration<double, std::micro>;
 
@@ -97,13 +108,9 @@ void benchmark_cell_group::advance(epoch ep,
     PL();
 };
 
-const std::vector<spike>& benchmark_cell_group::spikes() const {
-    return spikes_;
-}
+const std::vector<spike>& benchmark_cell_group::spikes() const { return spikes_; }
 
-void benchmark_cell_group::clear_spikes() {
-    spikes_.clear();
-}
+void benchmark_cell_group::clear_spikes() { spikes_.clear(); }
 
 void benchmark_cell_group::add_sampler(sampler_association_handle h,
                                    cell_member_predicate probeset_ids,

arbor/benchmark_cell_group.hpp

@@ -32,6 +32,8 @@ class benchmark_cell_group: public cell_group {
 
     void remove_all_samplers() override {}
 
+    void edit_cell(cell_gid_type gid, std::any edit) override;
+
     ARB_SERDES_ENABLE(benchmark_cell_group, cells_, spikes_, gids_);
 
     void t_serialize(serializer& ser, const std::string& k) const override;

arbor/cable_cell_group.cpp

@@ -18,6 +18,7 @@
 #include "util/partition.hpp"
 #include "util/range.hpp"
 #include "util/span.hpp"
+#include "util/maputil.hpp"
 
 namespace arb {
 
@@ -478,6 +479,20 @@ void cable_cell_group::remove_all_samplers() {
     sampler_map_.clear();
 }
 
+void
+cable_cell_group::edit_cell(cell_gid_type gid, std::any cell_edit) {
+    auto lid = util::binary_search_index(gids_, gid);
+    if (!lid) throw arb::arbor_internal_error{"gid " + std::to_string(gid) + " erroneuosly dispatched to cell group."};
+
+    if (auto cable_edit = std::any_cast<cable_cell_density_editor>(&cell_edit); cable_edit) {
+        lowered_->edit_density_parameter(gid, *lid, *cable_edit);
+    }
+    else {
+        throw bad_cell_edit(gid, "Not a Cable Cell editor (C++ type-id: '" + std::string{cell_edit.type().name()} + "')");
+    }
+
+}
+
 std::vector<probe_metadata> cable_cell_group::get_probe_metadata(const cell_address_type& probeset_id) const {
     // SAFETY: Probe associations are fixed after construction, so we do not
     //         need to grab the mutex.

arbor/cable_cell_group.hpp

@@ -44,6 +44,8 @@ struct ARB_ARBOR_API cable_cell_group: public cell_group {
 
     void remove_all_samplers() override;
 
+    void edit_cell(cell_gid_type gid, std::any edit) override;
+
     std::vector<probe_metadata> get_probe_metadata(const cell_address_type&) const override;
 
     ARB_SERDES_ENABLE(cable_cell_group, gids_, spikes_, lowered_);

arbor/cell_group.hpp

@@ -8,6 +8,7 @@
 #include <arbor/schedule.hpp>
 #include <arbor/spike.hpp>
 #include <arbor/serdes.hpp>
+#include <arbor/arbexcept.hpp>
 
 #include "epoch.hpp"
 #include "event_lane.hpp"
@@ -20,8 +21,7 @@
 //   ranges are needed to map (gid, label) pairs to their corresponding lid sets.
 namespace arb {
 
-class cell_group {
-public:
+struct cell_group {
     virtual ~cell_group() = default;
 
     virtual cell_kind get_cell_kind() const = 0;
@@ -32,20 +32,23 @@ class cell_group {
     virtual const std::vector<spike>& spikes() const = 0;
     virtual void clear_spikes() = 0;
 
-    // Sampler association methods below should be thread-safe, as they might be invoked
-    // from a sampler call back called from a different cell group running on a different thread.
-
+    // Sampler association methods below should be thread-safe, as they might be
+    // invoked from a sampler call back called from a different cell group
+    // running on a different thread.
     virtual void add_sampler(sampler_association_handle, cell_member_predicate, schedule, sampler_function) = 0;
     virtual void remove_sampler(sampler_association_handle) = 0;
     virtual void remove_all_samplers() = 0;
 
-    // Probe metadata queries might also be called while a simulation is running, and so should
-    // also be thread-safe.
+    // allow editing of certain cell properties
+    virtual void edit_cell(cell_gid_type gid, std::any edit) = 0;
 
+    // Probe metadata queries might also be called while a simulation is
+    // running, and so should also be thread-safe.
     virtual std::vector<probe_metadata> get_probe_metadata(const cell_address_type&) const { return {}; }
+
     // trampolines for serialization
     virtual void t_serialize(serializer& s, const std::string&) const = 0;
-    virtual void t_deserialize(serializer& s, const std::string&)  = 0;
+    virtual void t_deserialize(serializer& s, const std::string&) = 0;
 };
 
 using cell_group_ptr = std::unique_ptr<cell_group>;

arbor/fvm_lowered_cell.hpp

@@ -231,12 +231,15 @@ struct fvm_initialization_data {
 struct fvm_lowered_cell {
     virtual void reset() = 0;
 
-    virtual fvm_initialization_data initialize(const std::vector<cell_gid_type>& gids, const recipe& rec) = 0;
+    virtual fvm_initialization_data initialize(const std::vector<cell_gid_type>& gids,
+                                               const recipe& rec) = 0;
 
     virtual fvm_integration_result integrate(const timestep_range& dts,
                                              const event_lane_subrange& event_lanes,
                                              const std::vector<std::vector<sample_event>>& staged_samples) = 0;
 
+    virtual void edit_density_parameter(cell_gid_type gid, cell_gid_type lid, const cable_cell_density_editor& edit) = 0;
+
     virtual arb_value_type time() const = 0;
 
     virtual ~fvm_lowered_cell() {}

arbor/fvm_lowered_cell_impl.hpp

@@ -12,6 +12,8 @@
 #include <utility>
 #include <vector>
 
+#include <iostream>
+
 #include <arbor/assert.hpp>
 #include <arbor/common_types.hpp>
 #include <arbor/cable_cell.hpp>
@@ -55,8 +57,7 @@ struct fvm_lowered_cell_impl: public fvm_lowered_cell {
 
     value_type time() const override { return state_->time; }
 
-    //Exposed for testing purposes
-    std::vector<mechanism_ptr>& mechanisms() { return mechanisms_; }
+    void edit_density_parameter(cell_gid_type, cell_gid_type, const cable_cell_density_editor&) override;
 
     ARB_SERDES_ENABLE(fvm_lowered_cell_impl<Backend>, seed_, state_);
 
@@ -317,6 +318,24 @@ fvm_lowered_cell_impl<Backend>::add_probes(const std::vector<cell_gid_type>& gid
     }
 }
 
+template <typename Backend> void
+fvm_lowered_cell_impl<Backend>::edit_density_parameter(cell_gid_type gid,
+                                                       cell_gid_type lid,
+                                                       const cable_cell_density_editor& edit) {
+    std::cerr << "gid=" << gid << " lid=" << lid << '\n';
+    auto mech = std::find_if(mechanisms_.begin(), mechanisms_.end(), [&edit](const auto& m) { return m->internal_name() == edit.mechanism; });
+    if (mech == mechanisms_.end()) throw bad_cell_edit{gid, "no such mechanism: " + edit.mechanism};
+    auto ptr = mech->get();
+    if (ptr->kind() != arb_mechanism_kind_density) throw bad_cell_edit{gid, "not a density mechanism: " + edit.mechanism};
+    auto params = util::make_range(ptr->mech_.parameters, ptr->mech_.parameters + ptr->mech_.n_parameters);
+    for (const auto& [key, val]: edit.values) {
+        auto param = std::find_if(params.begin(), params.end(), [&key](const auto& p) { return p.name == key; });
+        if (params.end() == param) throw bad_cell_edit{gid, "no paramter " + key + " in mechanism: " + edit.mechanism};
+        auto pid = param - params.begin();
+        state_->update_density_data(lid, ptr->ppack_, pid, val);
+    }
+}
+
 template <typename Backend> fvm_initialization_data
 fvm_lowered_cell_impl<Backend>::initialize(const std::vector<cell_gid_type>& gids,
                                            const recipe& rec) {
@@ -402,18 +421,19 @@ fvm_lowered_cell_impl<Backend>::initialize(const std::vector<cell_gid_type>& gid
         auto it = util::max_element_by(fvm_info.num_sources, [](auto elem) {return util::second(elem);});
         max_detector = it->second;
     }
-    std::vector<arb_index_type> src_to_spike, cv_to_cell;
 
+    std::vector<arb_index_type> src_to_spike;
     if (post_events_) {
         for (auto cell_idx: make_span(ncell)) {
             for (auto lid: make_span(fvm_info.num_sources[gids[cell_idx]])) {
                 src_to_spike.push_back(cell_idx * max_detector + lid);
             }
         }
         src_to_spike.shrink_to_fit();
-        cv_to_cell = D.geometry.cv_to_cell;
     }
 
+    auto cv_to_cell = D.geometry.cv_to_cell;
+
     // map control volume ids to global cell ids
     std::vector<arb_index_type> cv_to_gid(D.geometry.cv_to_cell.size());
     std::transform(D.geometry.cv_to_cell.begin(), D.geometry.cv_to_cell.end(),
@@ -439,9 +459,9 @@ fvm_lowered_cell_impl<Backend>::initialize(const std::vector<cell_gid_type>& gid
                                             data_alignment? data_alignment: 1u,
                                             seed_);
 
+    target_handles_.resize(mech_data.n_target);
     // Keep track of mechanisms by name for probe lookup.
     std::unordered_map<std::string, mechanism*> mechptr_by_name;
-    target_handles_.resize(mech_data.n_target);
 
     unsigned mech_id = 0;
     for (const auto& [name, config]: mech_data.mechanisms) {
@@ -535,7 +555,7 @@ fvm_lowered_cell_impl<Backend>::initialize(const std::vector<cell_gid_type>& gid
                 voltage_mechanisms_.emplace_back(mech.release());
                 break;
             }
-            default:;
+            default:
                 throw invalid_mechanism_kind(config.kind);
         }
     }

arbor/include/arbor/arbexcept.hpp

@@ -47,6 +47,12 @@ struct ARB_SYMBOL_VISIBLE resolution_disabled: arbor_exception {
     {}
 };
 
+struct ARB_SYMBOL_VISIBLE bad_cell_edit: arbor_exception {
+    bad_cell_edit(cell_gid_type gid, std::string why):
+        arbor_exception("Cannot edit cell gid=" + std::to_string(gid) + ": " + why)
+    {}
+    std::string why;
+};
 
 struct ARB_SYMBOL_VISIBLE dup_cell_probe: arbor_exception {
     dup_cell_probe(cell_kind kind, cell_gid_type gid, cell_tag_type tag);