Add support for LFP reports

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,7 +111,7 @@ 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;
         char line_buf[max_line_length], name[max_line_length];
 
@@ -123,14 +124,19 @@ class FileHandler {
 
         if (nseg) {
             std::vector<int> sec, seg;
+            std::vector<double> lfp_factors;
             sec.reserve(nseg);
             seg.reserve(nseg);
+            lfp_factors.reserve(nseg);
 
             read_array<int>(&sec[0], nseg);
             read_array<int>(&seg[0], nseg);
+            read_array<double>(&lfp_factors[0], nseg);
 
             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]);
             }
         }
         return nseg;

coreneuron/io/nrn_setup.cpp

@@ -953,7 +953,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, double> lfp_factors;
+
     CellMapping(int g)
         : gid(g) {}
 
@@ -137,6 +141,11 @@ 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});
+    }
+
     ~CellMapping() {
         for (size_t i = 0; i < secmapvec.size(); i++) {
             delete secmapvec[i];
@@ -153,6 +162,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 +195,10 @@ 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);
+    }
 };
 }  // 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,49 @@ 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 count = 0;
+            double sum = 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;
+                }
+                sum += (fast_imem_rhs[segment_id] + iclamp) * factor;
+                count++;
+            }
+            *(to_report.front().var_value) = sum;
+        }
+    }
+}
+
 /** 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,15 @@ 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:
+                // 1 lfp value per gid
+                mapinfo->_lfp.resize(nt.ncell);
+                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 +77,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 +352,35 @@ 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 {
+    auto& summation_report = nt.summation_report_handler_->summation_reports_[report.output_path];
+    VarsToReport vars_to_report;
+    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];
+        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));
+            }
+        }
+        std::vector<VarWithMapping> to_report;
+        double* variable = report_variable + i;
+        to_report.push_back(VarWithMapping(i, variable));
+        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: