Record online LFP from multiple electrodes

Author: jorblancoa

coreneuron/io/nrn_filehandler.hpp

@@ -113,30 +113,38 @@ class FileHandler {
     template <typename T>
     int read_mapping_info(T* mapinfo, NrnThreadMappingInfo* ntmapping, CellMapping* cmap) {
         int nsec, nseg, n_scan;
+        size_t total_lfp_factors;
         char line_buf[max_line_length], name[max_line_length];
 
         F.getline(line_buf, sizeof(line_buf));
-        n_scan = sscanf(line_buf, "%s %d %d", name, &nsec, &nseg);
+        n_scan = sscanf(line_buf, "%s %d %d %zd", name, &nsec, &nseg, &total_lfp_factors);
 
-        nrn_assert(n_scan == 3);
+        nrn_assert(n_scan == 4);
 
         mapinfo->name = std::string(name);
 
         if (nseg) {
             std::vector<int> sec, seg;
             std::vector<double> lfp_factors;
+
             sec.reserve(nseg);
             seg.reserve(nseg);
-            lfp_factors.reserve(nseg);
+            lfp_factors.reserve(total_lfp_factors);
 
             read_array<int>(&sec[0], nseg);
             read_array<int>(&seg[0], nseg);
-            read_array<double>(&lfp_factors[0], nseg);
+            read_array<double>(&lfp_factors[0], total_lfp_factors);
+            int num_electrodes = read_int();
 
+            int factor_offset = 0;
             for (int i = 0; i < nseg; i++) {
                 mapinfo->add_segment(sec[i], seg[i]);
                 ntmapping->add_segment_id(seg[i]);
-                cmap->add_segment_lfp_factor(seg[i], lfp_factors[i]);
+                int factor_offset = i * num_electrodes;
+                if (total_lfp_factors > 0) {
+                    std::vector<double> segment_factors(lfp_factors.begin() + factor_offset, lfp_factors.begin() + factor_offset + num_electrodes);
+                    cmap->add_segment_lfp_factor(seg[i], segment_factors);
+                }
             }
         }
         return nseg;

coreneuron/io/nrnsection_mapping.hpp

@@ -75,7 +75,7 @@ struct CellMapping {
     std::vector<SecMapping*> secmapvec;
 
     /** map containing segment ids an its respective lfp factors */
-    std::unordered_map<int, double> lfp_factors;
+    std::unordered_map<int, std::vector<double>> lfp_factors;
 
     CellMapping(int g)
         : gid(g) {}
@@ -100,6 +100,15 @@ struct CellMapping {
                                });
     }
 
+    /** @brief return the number of electrodes in the lfp_factors map **/
+    int num_electrodes() const {
+        int num_electrodes = 0;
+        if (!lfp_factors.empty()) {
+            num_electrodes = lfp_factors.begin()->second.size();
+        }
+        return num_electrodes;
+    }
+
     /** @brief number of section lists */
     size_t size() const noexcept {
         return secmapvec.size();
@@ -141,9 +150,9 @@ struct CellMapping {
         return count;
     }
 
-    /** @brief add the lfp factor of a segment_id */
-    void add_segment_lfp_factor(const int segment_id, double factor) {
-        lfp_factors.insert({segment_id, factor});
+    /** @brief add the lfp electrode factors of a segment_id */
+    void add_segment_lfp_factor(const int segment_id, std::vector<double> factors) {
+        lfp_factors.insert({segment_id, factors});
     }
 
     ~CellMapping() {
@@ -200,5 +209,14 @@ struct NrnThreadMappingInfo {
     void add_segment_id(const int segment_id) {
         segment_ids.push_back(segment_id);
     }
+
+    /** @brief Resize the lfp vector */
+    void prepare_lfp() {
+        size_t lfp_size = 0;
+        for (const auto& mapping: mappingvec) {
+            lfp_size += mapping->num_electrodes();
+        }
+        _lfp.resize(lfp_size);
+    }
 };
 }  // namespace coreneuron

coreneuron/io/reports/report_event.cpp

@@ -85,25 +85,29 @@ void ReportEvent::lfp_calc(NrnThread* nt) {
             int gid = kv.first;
             const auto& to_report = kv.second;
             const auto& cell_mapping = mapinfo->get_cell_mapping(gid);
-
-            int count = 0;
-            double sum = 0.0;
+            int num_electrodes = cell_mapping->num_electrodes();
+            std::vector<double> lfp_values (num_electrodes, 0.0);
             for (const auto& kv: cell_mapping->lfp_factors) {
                 int segment_id = kv.first;
-                double factor = kv.second;
-                if (std::isnan(factor)) {
-                    factor = 0.0;
-                }
-                double iclamp = 0.0;
-                for (const auto& value: summation_report.currents_[segment_id]) {
-                    double current_value = *value.first;
-                    int scale = value.second;
-                    iclamp += current_value * scale;
+                std::vector<double> factors = kv.second;
+                int electrode_id = 0;
+                for (auto& factor: factors) {
+                    if (std::isnan(factor)) {
+                        factor = 0.0;
+                    }
+                    double iclamp = 0.0;
+                    for (const auto& value: summation_report.currents_[segment_id]) {
+                        double current_value = *value.first;
+                        int scale = value.second;
+                        iclamp += current_value * scale;
+                    }
+                    lfp_values[electrode_id] += (fast_imem_rhs[segment_id] + iclamp) * factor;
+                    electrode_id++;
                 }
-                sum += (fast_imem_rhs[segment_id] + iclamp) * factor;
-                count++;
             }
-            *(to_report.front().var_value) = sum;
+            for(int i=0; i < to_report.size(); i++) {
+                *(to_report[i].var_value) = lfp_values[i];
+            }
         }
     }
 }

coreneuron/io/reports/report_handler.cpp

@@ -60,8 +60,7 @@ void ReportHandler::create_report(double dt, double tstop, double delay) {
                 register_custom_report(nt, m_report_config, vars_to_report);
                 break;
             case LFPReport:
-                // 1 lfp value per gid
-                mapinfo->_lfp.resize(nt.ncell);
+                mapinfo->prepare_lfp();
                 vars_to_report =
                     get_lfp_vars_to_report(nt, m_report_config, mapinfo->_lfp.data(), nodes_to_gid);
                 is_soma_target = m_report_config.section_type == SectionType::Soma ||
@@ -356,8 +355,15 @@ VarsToReport ReportHandler::get_lfp_vars_to_report(const NrnThread& nt,
                                                    ReportConfiguration& report,
                                                    double* report_variable,
                                                    const std::vector<int>& nodes_to_gids) const {
+    const auto* mapinfo = static_cast<NrnThreadMappingInfo*>(nt.mapping);
+    if (!mapinfo) {
+        std::cerr << "[LFP] Error : mapping information is missing for a Cell group "
+                  << nt.ncell << '\n';
+        nrn_abort(1);
+    }
     auto& summation_report = nt.summation_report_handler_->summation_reports_[report.output_path];
     VarsToReport vars_to_report;
+    off_t offset_lfp = 0;
     for (int i = 0; i < nt.ncell; i++) {
         int gid = nt.presyns[i].gid_;
         if (report.target.find(gid) == report.target.end()) {
@@ -373,10 +379,22 @@ VarsToReport ReportHandler::get_lfp_vars_to_report(const NrnThread& nt,
                 summation_report.currents_[segment_id].push_back(std::make_pair(var_value, -1));
             }
         }
+        const auto& cell_mapping = mapinfo->get_cell_mapping(gid);
+        if (cell_mapping == nullptr) {
+            std::cerr
+                << "[LFP] Error : Compartment mapping information is missing for gid "
+                << gid << '\n';
+            nrn_abort(1);
+        }
         std::vector<VarWithMapping> to_report;
-        double* variable = report_variable + i;
-        to_report.push_back(VarWithMapping(i, variable));
-        vars_to_report[gid] = to_report;
+        int num_electrodes = cell_mapping->num_electrodes();
+        for (int electrode_id = 0; electrode_id < num_electrodes; electrode_id++) {
+            to_report.emplace_back(VarWithMapping(electrode_id, report_variable + offset_lfp));
+            offset_lfp++;
+        }
+        if (!to_report.empty()) {
+            vars_to_report[gid] = to_report;
+        }
     }
     return vars_to_report;
 }