Add lfp report type

Author: jorblancoa

coreneuron/io/nrn_filehandler.hpp

@@ -134,6 +134,9 @@ class FileHandler {
             read_array<int>(&seg[0], nseg);
             read_array<double>(&lfp_factors[0], nseg);
 
+            std::cout << "=====> NEW CoreNEURON!" << std::endl;
+            std::cout << "nseg = " << nseg << std::endl;
+
             for (int i = 0; i < nseg; i++) {
                 mapinfo->add_segment(sec[i], seg[i]);
                 ntmapping->add_segment_id(seg[i]);

coreneuron/io/nrnsection_mapping.hpp

@@ -161,6 +161,8 @@ struct NrnThreadMappingInfo {
     /** map containing segment ids an its respective lfp factors */
     std::unordered_map<int, double> lfp_factors;
 
+    std::vector<double> _lfp;
+
     /** @brief number of cells */
     size_t size() const {
         return mappingvec.size();

coreneuron/io/reports/nrnreport.hpp

@@ -71,7 +71,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

@@ -140,6 +140,9 @@ std::vector<ReportConfiguration> create_report_configurations(
             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
@@ -72,12 +73,32 @@ 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 sum = 0.0;
+        double* fast_imem_rhs = nt->nrn_fast_imem->nrn_sav_rhs;
+        for (const auto& kv: mapinfo->lfp_factors) {
+                int segment_id = kv.first;
+                double factor = kv.second;
+                if(std::isnan(factor)) {
+                    factor = 0.0;
+                }
+                std::cout << segment_id << " - " << factor << std::endl;
+                sum += fast_imem_rhs[segment_id] * factor;
+        }
+        mapinfo->_lfp[0] = sum;
+    }mdi
+}
+
 /** 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);
+        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

@@ -42,6 +42,7 @@ class ReportEvent: public DiscreteEvent {
     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;

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,13 @@ 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->_lfp.resize(12);
+                vars_to_report = get_lfp_vars_to_report(nt, m_report_config, mapinfo->_lfp.data());
+                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);
@@ -341,6 +349,34 @@ 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 {
+    VarsToReport vars_to_report;
+    /*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);
+    }*/
+    for (int i = 0; i < nt.ncell; i++) {
+        int gid = nt.presyns[i].gid_;
+        if (report.target.find(gid) == report.target.end()) {
+            continue;
+        }
+
+        std::vector<VarWithMapping> to_report;
+        // Add all electrodes to the first gid for now
+        std::vector<int> electrode_ids = {0};
+        for (const auto& electrode_id : electrode_ids) {
+            double* variable = report_variable + electrode_id;
+            to_report.push_back(VarWithMapping(electrode_id, 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,9 @@ 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> map_gids(const NrnThread& nt) const;
 #endif  // defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS)
   protected:

coreneuron/io/reports/sonata_report_handler.cpp

@@ -41,6 +41,7 @@ void SonataReportHandler::register_custom_report(const NrnThread& nt,
 void SonataReportHandler::register_report(const NrnThread& nt,
                                           ReportConfiguration& config,
                                           const VarsToReport& vars_to_report) {
+    std::cout << "Registering report " << config.output_path.data() << std::endl;
     sonata_create_report(config.output_path.data(),
                          config.start,
                          config.stop,