Calculate 1 lfp value per gid/timestep

Author: jorblancoa

coreneuron/io/nrn_filehandler.hpp

@@ -111,7 +111,7 @@ class FileHandler {
      * Read count no of mappings for section to segment
      */
     template <typename T>
-    int read_mapping_info(T* mapinfo, NrnThreadMappingInfo* ntmapping) {
+    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];
 
@@ -133,13 +133,10 @@ 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]);
-                ntmapping->add_segment_lfp_factor(seg[i], lfp_factors[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, ntmapping);
+            F.read_mapping_info(smap, ntmapping, cmap);
             cmap->add_sec_map(smap);
         }
 

coreneuron/io/nrnsection_mapping.hpp

@@ -74,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) {}
 
@@ -138,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];
@@ -157,9 +165,6 @@ struct NrnThreadMappingInfo {
     /** list of segment ids */
     std::vector<int> segment_ids;
 
-    /** map containing segment ids an its respective lfp factors */
-    std::unordered_map<int, double> lfp_factors;
-
     std::vector<double> _lfp;
 
     /** @brief number of cells */
@@ -195,10 +200,5 @@ struct NrnThreadMappingInfo {
     void add_segment_id(const int segment_id) {
         segment_ids.push_back(segment_id);
     }
-
-    /** @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});
-    }
 };
 }  // namespace coreneuron

coreneuron/io/reports/report_event.cpp

@@ -25,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;
@@ -77,19 +79,27 @@ 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;
+
+        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;
+                    }
+                    sum += fast_imem_rhs[segment_id] * factor;
+                    count++;
+            }
+            *(to_report.front().var_value) = sum;
         }
-        mapinfo->_lfp[0] = sum;
-    }mdi
+    }
 }
 
 /** on deliver, call ReportingLib and setup next event */
@@ -98,7 +108,9 @@ void ReportEvent::deliver(double t, NetCvode* nc, NrnThread* nt) {
 #pragma omp critical
     {
         summation_alu(nt);
-        lfp_calc(nt);
+        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,7 +36,8 @@ 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;
@@ -53,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

@@ -60,7 +60,8 @@ void ReportHandler::create_report(double dt, double tstop, double delay) {
                 register_custom_report(nt, m_report_config, vars_to_report);
                 break;
             case LFPReport:
-                mapinfo->_lfp.resize(12);
+                // 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());
                 is_soma_target = m_report_config.section_type == SectionType::Soma ||
                                  m_report_config.section_type == SectionType::Cell;
@@ -75,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));
         }
@@ -353,25 +355,15 @@ 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));
-        }
+        double* variable = report_variable + i;
+        to_report.push_back(VarWithMapping(i, variable));
         vars_to_report[gid] = to_report;
     }
     return vars_to_report;