Mutable Parameters

.editorconfig

@@ -0,0 +1,14 @@
+root = true
+
+# Unix-style newlines with a newline ending every file
+[*]
+end_of_line = lf
+insert_final_newline = true
+
+[*.py]
+charset = utf-8
+
+# 4 space indentation
+[*{.py,.cpp,.cu,.hpp}]
+indent_style = space
+indent_size = 4

.github/workflows/test-matrix.yml

@@ -29,8 +29,8 @@ jobs:
         - {
             name:  "Linux Min Clang",
             os:    "ubuntu-22.04",
-            cc:    "clang-13",
-            cxx:   "clang++-13",
+            cc:    "clang-14",
+            cxx:   "clang++-14",
             py:    "3.10",
             cmake: "4.0.x",
             mpi:   "ON",

arbor/adex_cell_group.cpp

@@ -3,6 +3,7 @@
 #include <arbor/arbexcept.hpp>
 
 #include "arbor/math.hpp"
+#include "util/maputil.hpp"
 #include "util/rangeutil.hpp"
 #include "util/span.hpp"
 #include "label_resolution.hpp"
@@ -99,6 +100,49 @@ void adex_cell_group::reset() {
     spikes_.clear();
 }
 
+void
+adex_cell_group::edit_cell(cell_gid_type gid, std::any cell_edit) {
+    try {
+        auto adex_edit = std::any_cast<adex_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."};
+        auto& lowered = cells_[*lid];
+        auto tmp = adex_cell {
+            .source = lowered.source, // Label of source
+            .target = lowered.target, // Label of target
+            .delta  = lowered.delta * U::mV,
+            .V_th   = lowered.V_th  * U::mV,
+            .C_m    = lowered.C_m   * U::nF,
+            .E_L    = lowered.E_L   * U::mV,
+            .E_R    = lowered.E_R   * U::mV,
+            .V_m    = lowered.V_m   * U::mV,
+            .t_ref  = lowered.t_ref * U::ms,
+            .g      = lowered.g     * U::uS,
+            .tau    = lowered.tau   * U::ms,
+            .w      = lowered.w     * U::pA,
+            .a      = lowered.a     * U::uS,
+            .b      = lowered.b     * U::nA,
+        };
+        adex_edit(tmp);
+        // NOTE: we forbid writing to V_m? Reasons
+        //       * the cell might be in the refractory period which causes semantic issues
+        //         - return to normal or not?
+        //         - what should probes return
+        //       * V_m is the _initial state_ only
+        if (tmp.V_m.value_as(U::mV) != lowered.V_m) throw bad_cell_edit(gid, "Initial voltage is not editable.");
+        if (tmp.w.value_as(U::pA) != lowered.w) throw bad_cell_edit(gid, "Adaption parameter is not editable.");
+        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 = adex_lowered_cell{tmp};
+
+    }
+    catch (const std::bad_any_cast& ){
+        throw bad_cell_edit(gid, "Not an AdEx editor (C++ type-id: '" + std::string{cell_edit.type().name()} + "')");
+    }
+}
+
+
 // integrate a single cell's state from current time `cur` tos final time `end`.
 // Extra parameters
 // * the cell cannot be updated until time `nxt`, which might be in the past or future.
@@ -118,7 +162,7 @@ void integrate_until(adex_lowered_cell& cell, const time_type end, const time_ty
     auto delta = end - cur;
     // membrane potential deviation from resting value
     auto dE = cell.V_m - cell.E_L;
-    // leak current 
+    // leak current
     auto il = cell.g*dE;
     // spike current
     auto is = cell.g*cell.delta*exp((cell.V_m - cell.V_th)/cell.delta);

arbor/adex_cell_group.hpp

@@ -89,6 +89,8 @@ struct ARB_ARBOR_API adex_cell_group: public cell_group {
 
     static bool backend_supported(backend_kind kind) { return kind == backend_kind::multicore; }
 
+    void edit_cell(cell_gid_type gid, std::any edit) override;
+
 private:
     enum class adex_probe_kind { voltage, adaption };
 

arbor/backends/gpu/shared_state.cpp

@@ -381,6 +381,24 @@ std::pair<arb_value_type, arb_value_type> shared_state::voltage_bounds() const {
     return minmax_value_impl(n_cv, voltage.data());
 }
 
+void shared_state::update_range_parameter(unsigned mid,
+                                          arb_index_type lid,
+                                          arb_mechanism_ppack& ppack,
+                                          cell_gid_type pid,
+                                          const std::vector<arb_value_type>& vals) {
+    auto off = 0;
+    auto& store = storage[mid];
+    auto& param_d = store.parameters_[pid];
+    auto param_h = memory::on_host(memory::device_view<arb_value_type>(param_d, ppack.width));
+    for (auto idx = 0ul; idx < ppack.width; ++idx) {
+        if (lid != ppack.vec_ci[ppack.node_index[idx]]) continue;
+        param_h[idx] = vals[off];
+        ++off;
+    }
+    // conservative copy of the _full_ array...
+    memory::copy(param_h.data(), param_d);
+}
+
 void shared_state::take_samples() {
     sample_events.mark();
     if (!sample_events.empty()) {

arbor/backends/gpu/shared_state.hpp

@@ -233,6 +233,12 @@ struct ARB_ARBOR_API shared_state: shared_state_base<shared_state, array, ion_st
         sample_time_host  = memory::on_host(sample_time);
         sample_value_host = memory::on_host(sample_value);
     }
+
+    void update_range_parameter(unsigned mid,                              // mechanism id
+                                arb_index_type cell_id,                    // to filter CVs of the cell we are editing
+                                arb_mechanism_ppack& ppack,                // ppack, used to read from CV -> cell_id
+                                cell_gid_type pid,                         // index of parameter
+                                const std::vector<arb_value_type>& val);   // values we are writing
 };
 
 // For debugging only

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_,
@@ -403,10 +401,8 @@ unsigned 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();
@@ -537,5 +533,19 @@ unsigned shared_state::instantiate(arb::mechanism& m,
     return id;
 }
 
+void shared_state::update_range_parameter(unsigned mid,
+                                          arb_index_type lid,
+                                          arb_mechanism_ppack& ppack,
+                                          cell_gid_type pid,
+                                          const std::vector<arb_value_type>& vals) {
+    auto off = 0;
+    auto param = storage[mid].parameters_[pid];
+    for (auto idx = 0ul; idx < ppack.width; ++idx) {
+        if (lid != ppack.vec_ci[ppack.node_index[idx]]) continue;
+        param[idx] = vals[off];
+        ++off;
+    }
+}
+
 } // namespace multicore
 } // namespace arb

arbor/backends/multicore/shared_state.hpp

@@ -238,6 +238,14 @@ struct ARB_ARBOR_API shared_state:
         sample_time_host = util::range_pointer_view(sample_time);
         sample_value_host = util::range_pointer_view(sample_value);
     }
+
+    bool mechanism_matches();
+
+    void update_range_parameter(unsigned mid,                              // mechanism id
+                                arb_index_type cell_id,                    // to filter CVs of the cell we are editing
+                                arb_mechanism_ppack& ppack,                // ppack, used to read from CV -> cell_id
+                                cell_gid_type pid,                         // index of parameter
+                                const std::vector<arb_value_type>& val);   // values we are writing
 };
 
 // For debugging only:

arbor/backends/shared_state_base.hpp

@@ -42,7 +42,7 @@ 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_time  = array(n_samples);
             d->sample_value = array(n_samples);
         }
         initialize(samples, d->sample_events);
@@ -80,6 +80,12 @@ struct shared_state_base {
         }
     }
 
+    // overwrite a RANGE-type parameter in the mechanism `ppack` is attached to.
+    void update_range_parameter(cell_gid_type lid, arb_mechanism_ppack& ppack, cell_gid_type pid, const std::vector<arb_value_type>& vals) {
+        auto d = static_cast<D*>(this);
+        d->update_range_parameter(lid, ppack, pid, vals);
+    }
+
     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(cell);
 };
 
-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.cpp

@@ -81,6 +81,8 @@ struct cable_cell_impl {
     using index_type = cable_cell::index_type;
     using size_type  = cable_cell::size_type;
 
+    bool is_mut_ = false;
+
     // The label dictionary.
     label_dict dictionary;
 
@@ -117,7 +119,8 @@ struct cable_cell_impl {
 
     // Discretization
     std::optional<cv_policy> discretization_;
-    cable_cell_impl(const arb::morphology& m, const label_dict& labels, const decor& decorations, const std::optional<cv_policy>& cvp):
+    cable_cell_impl(const arb::morphology& m, const label_dict& labels, const decor& decorations, const std::optional<cv_policy>& cvp, cable_cell_mutability is_mut=cable_cell_mutability::disabled):
+        is_mut_(is_mut == cable_cell_mutability::enabled),
         dictionary(labels),
         provider(m, dictionary),
         decorations(decorations),
@@ -239,8 +242,12 @@ void cable_cell_impl::init() {
     }
 }
 
-cable_cell::cable_cell(const arb::morphology& m, const decor& decorations, const label_dict& dictionary, const std::optional<cv_policy>& cvp):
-    impl_(make_impl(new cable_cell_impl(m, dictionary, decorations, cvp)))
+cable_cell::cable_cell(const arb::morphology& m,
+                       const decor& decorations,
+                       const label_dict& dictionary,
+                       const std::optional<cv_policy>& cvp,
+                       cable_cell_mutability is_mut):
+    impl_(make_impl(new cable_cell_impl(m, dictionary, decorations, cvp, is_mut)))
 {}
 
 cable_cell::cable_cell(): impl_(make_impl(new cable_cell_impl())) {}
@@ -249,6 +256,8 @@ cable_cell::cable_cell(const cable_cell& other):
     impl_(make_impl(new cable_cell_impl(*other.impl_)))
 {}
 
+bool cable_cell::is_editable() const { return impl_->is_mut_; }
+
 const label_dict& cable_cell::labels() const { return impl_->dictionary; }
 const concrete_embedding& cable_cell::embedding() const { return impl_->provider.embedding(); }
 const arb::morphology& cable_cell::morphology() const { return impl_->provider.morphology(); }

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 {
 
@@ -491,4 +492,22 @@ std::vector<probe_metadata> cable_cell_group::get_probe_metadata(const cell_addr
     }
     return result;
 }
+
+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."};    
+    try {
+        auto cc_edit = std::any_cast<cable_cell_editor>(cell_edit);
+        lowered_->edit_cell(gid, *lid, cc_edit);
+    }
+    catch (std::bad_any_cast& ex) {
+        throw bad_cell_edit(gid, "Not a Cable Cell editor (C++ type-id: '"
+                                 + std::string{cell_edit.type().name()}
+                                 + " ./. "
+                                 + std::string{typeid(cable_cell_editor).name()}
+                                 + "')");
+    }
+}
+
 } // namespace arb