Add support for on-line LFP calculations

CMakeLists.txt

@@ -10,7 +10,7 @@ cmake_minimum_required(VERSION 3.15 FATAL_ERROR)
 # standalone.
 project(
   coreneuron
-  VERSION 9.0.0
+  VERSION 8.2.1
   LANGUAGES CXX)
 
 # ~~~

coreneuron/io/nrn_filehandler.hpp

@@ -14,6 +14,7 @@
 #include <sys/stat.h>
 
 #include "coreneuron/utils/nrn_assert.h"
+#include "coreneuron/io/nrnsection_mapping.hpp"
 
 namespace coreneuron {
 /** Encapsulate low-level reading of coreneuron input data files.
@@ -110,27 +111,45 @@ class FileHandler {
      * Read count no of mappings for section to segment
      */
     template <typename T>
-    int read_mapping_info(T* mapinfo) {
+    int read_mapping_info(T* mapinfo, NrnThreadMappingInfo* ntmapping, CellMapping* cmap) {
         int nsec, nseg, n_scan;
+        size_t total_lfp_factors;
+        int num_electrodes;
         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 %d", name, &nsec, &nseg, &total_lfp_factors, &num_electrodes);
 
-        nrn_assert(n_scan == 3);
+        nrn_assert(n_scan == 5);
 
         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(total_lfp_factors);
 
             read_array<int>(&sec[0], nseg);
             read_array<int>(&seg[0], nseg);
+            if (total_lfp_factors > 0) {
+                read_array<double>(&lfp_factors[0], total_lfp_factors);
+            }
 
+            int factor_offset = 0;
             for (int i = 0; i < nseg; i++) {
                 mapinfo->add_segment(sec[i], seg[i]);
+                ntmapping->add_segment_id(seg[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/nrn_setup.cpp

@@ -959,7 +959,7 @@ void read_phase3(NrnThread& nt, UserParams& userParams) {
         // read section-segment mapping for every section list
         for (int j = 0; j < nseclist; j++) {
             SecMapping* smap = new SecMapping();
-            F.read_mapping_info(smap);
+            F.read_mapping_info(smap, ntmapping, cmap);
             cmap->add_sec_map(smap);
         }
 

coreneuron/io/nrnsection_mapping.hpp

@@ -14,6 +14,7 @@
 #include <vector>
 #include <map>
 #include <iostream>
+#include <unordered_map>
 
 namespace coreneuron {
 
@@ -73,6 +74,9 @@ struct CellMapping {
     /** list of section lists (like soma, axon, apic) */
     std::vector<SecMapping*> secmapvec;
 
+    /** map containing segment ids an its respective lfp factors */
+    std::unordered_map<int, std::vector<double>> lfp_factors;
+
     CellMapping(int g)
         : gid(g) {}
 
@@ -96,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();
@@ -137,6 +150,11 @@ struct CellMapping {
         return count;
     }
 
+    /** @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() {
         for (size_t i = 0; i < secmapvec.size(); i++) {
             delete secmapvec[i];
@@ -153,6 +171,11 @@ struct NrnThreadMappingInfo {
     /** list of cells mapping */
     std::vector<CellMapping*> mappingvec;
 
+    /** list of segment ids */
+    std::vector<int> segment_ids;
+
+    std::vector<double> _lfp;
+
     /** @brief number of cells */
     size_t size() const {
         return mappingvec.size();
@@ -181,5 +204,19 @@ struct NrnThreadMappingInfo {
     void add_cell_mapping(CellMapping* c) {
         mappingvec.push_back(c);
     }
+
+    /** @brief add a new segment */
+    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/nrnreport.hpp

@@ -76,7 +76,8 @@ enum ReportType {
     SynapseReport,
     IMembraneReport,
     SectionReport,
-    SummationReport
+    SummationReport,
+    LFPReport
 };
 
 // enumerate that defines the section type for a Section report

coreneuron/io/reports/report_configuration_parser.cpp

@@ -138,6 +138,9 @@ std::vector<ReportConfiguration> create_report_configurations(const std::string&
             report_type = SynapseReport;
         } else if (report.type_str == "summation") {
             report_type = SummationReport;
+        } else if (report.type_str == "lfp") {
+            nrn_use_fast_imem = true;
+            report_type = LFPReport;
         } else {
             std::cerr << "Report error: unsupported type " << report.type_str << std::endl;
             nrn_abort(1);

coreneuron/io/reports/report_event.cpp

@@ -10,6 +10,7 @@
 #include "coreneuron/sim/multicore.hpp"
 #include "coreneuron/io/reports/nrnreport.hpp"
 #include "coreneuron/utils/nrn_assert.h"
+#include "coreneuron/io/nrnsection_mapping.hpp"
 #ifdef ENABLE_BIN_REPORTS
 #include "reportinglib/Records.h"
 #endif  // ENABLE_BIN_REPORTS
@@ -24,11 +25,13 @@ ReportEvent::ReportEvent(double dt,
                          double tstart,
                          const VarsToReport& filtered_gids,
                          const char* name,
-                         double report_dt)
+                         double report_dt,
+                         ReportType type)
     : dt(dt)
     , tstart(tstart)
     , report_path(name)
     , report_dt(report_dt)
+    , report_type(type)
     , vars_to_report(filtered_gids) {
     nrn_assert(filtered_gids.size());
     step = tstart / dt;
@@ -72,12 +75,53 @@ void ReportEvent::summation_alu(NrnThread* nt) {
     }
 }
 
+void ReportEvent::lfp_calc(NrnThread* nt) {
+    // Calculate lfp only on reporting steps
+    if (step > 0 && (static_cast<int>(step) % reporting_period) == 0) {
+        auto* mapinfo = static_cast<NrnThreadMappingInfo*>(nt->mapping);
+        double* fast_imem_rhs = nt->nrn_fast_imem->nrn_sav_rhs;
+        auto& summation_report = nt->summation_report_handler_->summation_reports_[report_path];
+        for (const auto& kv: vars_to_report) {
+            int gid = kv.first;
+            const auto& to_report = kv.second;
+            const auto& cell_mapping = mapinfo->get_cell_mapping(gid);
+            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;
+                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++;
+                }
+            }
+            for (int i = 0; i < to_report.size(); i++) {
+                *(to_report[i].var_value) = lfp_values[i];
+            }
+        }
+    }
+}
+
 /** on deliver, call ReportingLib and setup next event */
 void ReportEvent::deliver(double t, NetCvode* nc, NrnThread* nt) {
 /* reportinglib is not thread safe */
 #pragma omp critical
     {
-        summation_alu(nt);
+        if (report_type == ReportType::SummationReport) {
+            summation_alu(nt);
+        } else if (report_type == ReportType::LFPReport) {
+            lfp_calc(nt);
+        }
         // each thread needs to know its own step
 #ifdef ENABLE_BIN_REPORTS
         records_nrec(step, gids_to_report.size(), gids_to_report.data(), report_path.data());

coreneuron/io/reports/report_event.hpp

@@ -36,12 +36,14 @@ class ReportEvent: public DiscreteEvent {
                 double tstart,
                 const VarsToReport& filtered_gids,
                 const char* name,
-                double report_dt);
+                double report_dt,
+                ReportType type);
 
     /** on deliver, call ReportingLib and setup next event */
     void deliver(double t, NetCvode* nc, NrnThread* nt) override;
     bool require_checkpoint() override;
     void summation_alu(NrnThread* nt);
+    void lfp_calc(NrnThread* nt);
 
   private:
     double dt;
@@ -52,6 +54,7 @@ class ReportEvent: public DiscreteEvent {
     std::vector<int> gids_to_report;
     double tstart;
     VarsToReport vars_to_report;
+    ReportType report_type;
 };
 #endif  // defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS)
 

coreneuron/io/reports/report_handler.cpp

@@ -35,6 +35,7 @@ void ReportHandler::create_report(double dt, double tstop, double delay) {
             continue;
         }
         const std::vector<int>& nodes_to_gid = map_gids(nt);
+        auto* mapinfo = static_cast<NrnThreadMappingInfo*>(nt.mapping);
         VarsToReport vars_to_report;
         bool is_soma_target;
         switch (m_report_config.type) {
@@ -58,6 +59,14 @@ void ReportHandler::create_report(double dt, double tstop, double delay) {
                                                               nodes_to_gid);
                 register_custom_report(nt, m_report_config, vars_to_report);
                 break;
+            case LFPReport:
+                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 ||
+                                 m_report_config.section_type == SectionType::Cell;
+                register_section_report(nt, m_report_config, vars_to_report, is_soma_target);
+                break;
             default:
                 vars_to_report = get_synapse_vars_to_report(nt, m_report_config, nodes_to_gid);
                 register_custom_report(nt, m_report_config, vars_to_report);
@@ -67,7 +76,8 @@ void ReportHandler::create_report(double dt, double tstop, double delay) {
                                                               t,
                                                               vars_to_report,
                                                               m_report_config.output_path.data(),
-                                                              m_report_config.report_dt);
+                                                              m_report_config.report_dt,
+                                                              m_report_config.type);
             report_event->send(t, net_cvode_instance, &nt);
             m_report_events.push_back(std::move(report_event));
         }
@@ -341,6 +351,55 @@ VarsToReport ReportHandler::get_synapse_vars_to_report(
     return vars_to_report;
 }
 
+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()) {
+            continue;
+        }
+        // IClamp is needed for the LFP calculation
+        auto mech_id = nrn_get_mechtype("IClamp");
+        Memb_list* ml = nt._ml_list[mech_id];
+        if (ml) {
+            for (int j = 0; j < ml->nodecount; j++) {
+                auto segment_id = ml->nodeindices[j];
+                if ((nodes_to_gids[segment_id] == gid)) {
+                    double* var_value = get_var_location_from_var_name(mech_id, "i", ml, j);
+                    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;
+        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;
+}
+
 // map GIDs of every compartment, it consist in a backward sweep then forward sweep algorithm
 std::vector<int> ReportHandler::map_gids(const NrnThread& nt) const {
     std::vector<int> nodes_gid(nt.end, -1);

coreneuron/io/reports/report_handler.hpp

@@ -44,6 +44,10 @@ class ReportHandler {
     VarsToReport get_synapse_vars_to_report(const NrnThread& nt,
                                             ReportConfiguration& report,
                                             const std::vector<int>& nodes_to_gids) const;
+    VarsToReport get_lfp_vars_to_report(const NrnThread& nt,
+                                        ReportConfiguration& report,
+                                        double* report_variable,
+                                        const std::vector<int>& nodes_to_gids) const;
     std::vector<int> map_gids(const NrnThread& nt) const;
 #endif  // defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS)
   protected: