|
7 | 7 | */ |
8 | 8 |
|
9 | 9 | #include "report_event.hpp" |
| 10 | + |
| 11 | +#include <numeric> |
| 12 | + |
10 | 13 | #include "coreneuron/sim/multicore.hpp" |
11 | 14 | #include "coreneuron/io/reports/nrnreport.hpp" |
12 | 15 | #include "coreneuron/utils/nrn_assert.h" |
@@ -77,33 +80,44 @@ void ReportEvent::summation_alu(NrnThread* nt) { |
77 | 80 | } |
78 | 81 | } |
79 | 82 |
|
| 83 | +/** @brief Compute Local Field Potentials (LFP) for each cell. |
| 84 | + * |
| 85 | + * For every segment of a cell, computes the total membrane current (imem + iclamp) |
| 86 | + * and accumulates the weighted contribution to each electrode using precomputed |
| 87 | + * transfer factors. Results are written into the report output buffers. |
| 88 | + */ |
80 | 89 | void ReportEvent::lfp_calc(NrnThread* nt) { |
81 | 90 | auto* mapinfo = static_cast<NrnThreadMappingInfo*>(nt->mapping); |
82 | 91 | double* fast_imem_rhs = nt->nrn_fast_imem->nrn_sav_rhs; |
83 | 92 | auto& summation_report = nt->summation_report_handler_->summation_reports_[report_path]; |
84 | 93 | for (const auto& kv: vars_to_report) { |
85 | 94 | int gid = kv.first; |
86 | | - const auto& to_report = kv.second; |
| 95 | + const auto& electrode_outputs = kv.second; |
87 | 96 | const auto& cell_mapping = mapinfo->get_cell_mapping(gid); |
88 | | - int num_electrodes = cell_mapping->num_electrodes(); |
89 | | - std::vector<double> lfp_values(num_electrodes, 0.0); |
90 | | - for (const auto& kv: cell_mapping->lfp_factors) { |
91 | | - int segment_id = kv.first; |
92 | | - const auto& factors = kv.second; |
93 | | - int electrode_id = 0; |
94 | | - for (const auto& factor: factors) { |
95 | | - double iclamp = 0.0; |
96 | | - for (const auto& value: summation_report.currents_[segment_id]) { |
97 | | - double current_value = *value.first; |
98 | | - int scale = value.second; |
99 | | - iclamp += current_value * scale; |
100 | | - } |
101 | | - lfp_values[electrode_id] += (fast_imem_rhs[segment_id] + iclamp) * factor; |
102 | | - electrode_id++; |
| 97 | + const auto n_electrodes = cell_mapping->num_electrodes(); |
| 98 | + const auto n_segments = cell_mapping->lfp_segment_ids.size(); |
| 99 | + std::vector<double> lfp_values(n_electrodes, 0.0); |
| 100 | + for (size_t i = 0; i < n_segments; i++) { |
| 101 | + const auto segment_id = cell_mapping->lfp_segment_ids[i]; |
| 102 | + |
| 103 | + // compute imem + iclamp |
| 104 | + const double imem = std::accumulate(summation_report.currents_[segment_id].begin(), |
| 105 | + summation_report.currents_[segment_id].end(), |
| 106 | + fast_imem_rhs[segment_id], |
| 107 | + [](double sum, const auto& value) { |
| 108 | + return sum + *value.first * value.second; |
| 109 | + }); |
| 110 | + |
| 111 | + // dot product with the factors |
| 112 | + const double* factors = &cell_mapping->lfp_factors_flat[i * n_electrodes]; |
| 113 | + for (size_t e = 0; e < n_electrodes; e++) { |
| 114 | + lfp_values[e] += imem * factors[e]; |
103 | 115 | } |
104 | 116 | } |
105 | | - for (int i = 0; i < to_report.size(); i++) { |
106 | | - *(to_report[i].var_value) = lfp_values[i]; |
| 117 | + |
| 118 | + // write LFP values to report output buffers |
| 119 | + for (size_t e = 0; e < electrode_outputs.size(); e++) { |
| 120 | + *(electrode_outputs[e].var_value) = lfp_values[e]; |
107 | 121 | } |
108 | 122 | } |
109 | 123 | } |
|
0 commit comments