From efa8a80ca769fb4a41da1a2e39758f10a640a74a Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso Date: Wed, 27 Oct 2021 16:07:58 +0200 Subject: [PATCH 1/6] Read lfp factors from neurodamus --- coreneuron/io/nrn_filehandler.hpp | 8 +++++++- coreneuron/io/nrn_setup.cpp | 2 +- coreneuron/io/nrnsection_mapping.hpp | 17 +++++++++++++++++ 3 files changed, 25 insertions(+), 2 deletions(-) diff --git a/coreneuron/io/nrn_filehandler.hpp b/coreneuron/io/nrn_filehandler.hpp index e7d280023..4783f5ce4 100644 --- a/coreneuron/io/nrn_filehandler.hpp +++ b/coreneuron/io/nrn_filehandler.hpp @@ -14,6 +14,7 @@ #include #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 - int read_mapping_info(T* mapinfo) { + int read_mapping_info(T* mapinfo, NrnThreadMappingInfo* ntmapping) { 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 sec, seg; + std::vector lfp_factors; sec.reserve(nseg); seg.reserve(nseg); + lfp_factors.reserve(nseg); read_array(&sec[0], nseg); read_array(&seg[0], nseg); + read_array(&lfp_factors[0], nseg); 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]); } } return nseg; diff --git a/coreneuron/io/nrn_setup.cpp b/coreneuron/io/nrn_setup.cpp index 63a3c86d2..0b6da62f4 100644 --- a/coreneuron/io/nrn_setup.cpp +++ b/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->add_sec_map(smap); } diff --git a/coreneuron/io/nrnsection_mapping.hpp b/coreneuron/io/nrnsection_mapping.hpp index d7524fc42..9756fca8d 100644 --- a/coreneuron/io/nrnsection_mapping.hpp +++ b/coreneuron/io/nrnsection_mapping.hpp @@ -14,6 +14,7 @@ #include #include #include +#include namespace coreneuron { @@ -153,6 +154,12 @@ struct NrnThreadMappingInfo { /** list of cells mapping */ std::vector mappingvec; + /** list of segment ids */ + std::vector segment_ids; + + /** map containing segment ids an its respective lfp factors */ + std::unordered_map lfp_factors; + /** @brief number of cells */ size_t size() const { return mappingvec.size(); @@ -181,5 +188,15 @@ 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 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 From d349f5d1f40adbd6aeed784ac711c60debffa687 Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso Date: Fri, 5 Nov 2021 17:32:51 +0100 Subject: [PATCH 2/6] Add lfp report type --- coreneuron/io/nrn_filehandler.hpp | 3 ++ coreneuron/io/nrnsection_mapping.hpp | 2 ++ coreneuron/io/reports/nrnreport.hpp | 3 +- .../reports/report_configuration_parser.cpp | 3 ++ coreneuron/io/reports/report_event.cpp | 21 +++++++++++ coreneuron/io/reports/report_event.hpp | 1 + coreneuron/io/reports/report_handler.cpp | 36 +++++++++++++++++++ coreneuron/io/reports/report_handler.hpp | 3 ++ .../io/reports/sonata_report_handler.cpp | 1 + 9 files changed, 72 insertions(+), 1 deletion(-) diff --git a/coreneuron/io/nrn_filehandler.hpp b/coreneuron/io/nrn_filehandler.hpp index 4783f5ce4..8f83d1e02 100644 --- a/coreneuron/io/nrn_filehandler.hpp +++ b/coreneuron/io/nrn_filehandler.hpp @@ -133,6 +133,9 @@ class FileHandler { read_array(&seg[0], nseg); read_array(&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]); diff --git a/coreneuron/io/nrnsection_mapping.hpp b/coreneuron/io/nrnsection_mapping.hpp index 9756fca8d..29836a24b 100644 --- a/coreneuron/io/nrnsection_mapping.hpp +++ b/coreneuron/io/nrnsection_mapping.hpp @@ -160,6 +160,8 @@ struct NrnThreadMappingInfo { /** map containing segment ids an its respective lfp factors */ std::unordered_map lfp_factors; + std::vector _lfp; + /** @brief number of cells */ size_t size() const { return mappingvec.size(); diff --git a/coreneuron/io/reports/nrnreport.hpp b/coreneuron/io/reports/nrnreport.hpp index 7cdc05806..ac85d75be 100644 --- a/coreneuron/io/reports/nrnreport.hpp +++ b/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 diff --git a/coreneuron/io/reports/report_configuration_parser.cpp b/coreneuron/io/reports/report_configuration_parser.cpp index 6c69662b7..528c269a6 100644 --- a/coreneuron/io/reports/report_configuration_parser.cpp +++ b/coreneuron/io/reports/report_configuration_parser.cpp @@ -138,6 +138,9 @@ std::vector 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); diff --git a/coreneuron/io/reports/report_event.cpp b/coreneuron/io/reports/report_event.cpp index 78eb6b268..139fc88b0 100644 --- a/coreneuron/io/reports/report_event.cpp +++ b/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(step) % reporting_period) == 0) { + auto* mapinfo = static_cast(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()); diff --git a/coreneuron/io/reports/report_event.hpp b/coreneuron/io/reports/report_event.hpp index 20d325614..cf2d7bd19 100644 --- a/coreneuron/io/reports/report_event.hpp +++ b/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; diff --git a/coreneuron/io/reports/report_handler.cpp b/coreneuron/io/reports/report_handler.cpp index 84341e7a4..fd385cd68 100644 --- a/coreneuron/io/reports/report_handler.cpp +++ b/coreneuron/io/reports/report_handler.cpp @@ -35,6 +35,7 @@ void ReportHandler::create_report(double dt, double tstop, double delay) { continue; } const std::vector& nodes_to_gid = map_gids(nt); + auto* mapinfo = static_cast(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(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 to_report; + // Add all electrodes to the first gid for now + std::vector 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 ReportHandler::map_gids(const NrnThread& nt) const { std::vector nodes_gid(nt.end, -1); diff --git a/coreneuron/io/reports/report_handler.hpp b/coreneuron/io/reports/report_handler.hpp index 084886c96..db62bf993 100644 --- a/coreneuron/io/reports/report_handler.hpp +++ b/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& nodes_to_gids) const; + VarsToReport get_lfp_vars_to_report(const NrnThread& nt, + ReportConfiguration& report, + double* report_variable) const; std::vector map_gids(const NrnThread& nt) const; #endif // defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS) protected: diff --git a/coreneuron/io/reports/sonata_report_handler.cpp b/coreneuron/io/reports/sonata_report_handler.cpp index 99f5709d5..bc38ef27b 100644 --- a/coreneuron/io/reports/sonata_report_handler.cpp +++ b/coreneuron/io/reports/sonata_report_handler.cpp @@ -57,6 +57,7 @@ std::pair SonataReportHandler::get_population_info(int gid) { 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, From ab5c2c345353201d1ff01b0d4ad21cf554368b8a Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso Date: Fri, 3 Dec 2021 11:52:28 +0100 Subject: [PATCH 3/6] Calculate 1 lfp value per gid/timestep --- coreneuron/io/nrn_filehandler.hpp | 7 ++--- coreneuron/io/nrn_setup.cpp | 2 +- coreneuron/io/nrnsection_mapping.hpp | 16 +++++----- coreneuron/io/reports/report_event.cpp | 38 ++++++++++++++++-------- coreneuron/io/reports/report_event.hpp | 4 ++- coreneuron/io/reports/report_handler.cpp | 20 ++++--------- 6 files changed, 45 insertions(+), 42 deletions(-) diff --git a/coreneuron/io/nrn_filehandler.hpp b/coreneuron/io/nrn_filehandler.hpp index 8f83d1e02..99908550a 100644 --- a/coreneuron/io/nrn_filehandler.hpp +++ b/coreneuron/io/nrn_filehandler.hpp @@ -111,7 +111,7 @@ class FileHandler { * Read count no of mappings for section to segment */ template - 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(&seg[0], nseg); read_array(&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; diff --git a/coreneuron/io/nrn_setup.cpp b/coreneuron/io/nrn_setup.cpp index 0b6da62f4..b078be1ef 100644 --- a/coreneuron/io/nrn_setup.cpp +++ b/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); } diff --git a/coreneuron/io/nrnsection_mapping.hpp b/coreneuron/io/nrnsection_mapping.hpp index 29836a24b..5b9c7bd57 100644 --- a/coreneuron/io/nrnsection_mapping.hpp +++ b/coreneuron/io/nrnsection_mapping.hpp @@ -74,6 +74,9 @@ struct CellMapping { /** list of section lists (like soma, axon, apic) */ std::vector secmapvec; + /** map containing segment ids an its respective lfp factors */ + std::unordered_map 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 segment_ids; - /** map containing segment ids an its respective lfp factors */ - std::unordered_map lfp_factors; - std::vector _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 diff --git a/coreneuron/io/reports/report_event.cpp b/coreneuron/io/reports/report_event.cpp index 139fc88b0..9deaff30c 100644 --- a/coreneuron/io/reports/report_event.cpp +++ b/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(step) % reporting_period) == 0) { auto* mapinfo = static_cast(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()); diff --git a/coreneuron/io/reports/report_event.hpp b/coreneuron/io/reports/report_event.hpp index cf2d7bd19..0f1a07358 100644 --- a/coreneuron/io/reports/report_event.hpp +++ b/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 gids_to_report; double tstart; VarsToReport vars_to_report; + ReportType report_type; }; #endif // defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS) diff --git a/coreneuron/io/reports/report_handler.cpp b/coreneuron/io/reports/report_handler.cpp index fd385cd68..1b69cd6ff 100644 --- a/coreneuron/io/reports/report_handler.cpp +++ b/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,12 +355,6 @@ 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(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()) { @@ -366,12 +362,8 @@ VarsToReport ReportHandler::get_lfp_vars_to_report(const NrnThread& nt, } std::vector to_report; - // Add all electrodes to the first gid for now - std::vector 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; From 84cdf5c67916805d8db92656fd6286a33ebc5b62 Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso Date: Thu, 6 Jan 2022 18:05:14 +0100 Subject: [PATCH 4/6] Remove cout message --- coreneuron/io/reports/sonata_report_handler.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/coreneuron/io/reports/sonata_report_handler.cpp b/coreneuron/io/reports/sonata_report_handler.cpp index bc38ef27b..99f5709d5 100644 --- a/coreneuron/io/reports/sonata_report_handler.cpp +++ b/coreneuron/io/reports/sonata_report_handler.cpp @@ -57,7 +57,6 @@ std::pair SonataReportHandler::get_population_info(int gid) { 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, From b6c8c2e78b541d7709208b70db57adad4f5d19bb Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso Date: Tue, 5 Apr 2022 18:32:03 +0200 Subject: [PATCH 5/6] Take into account IClamp for the lfp calculation --- coreneuron/io/reports/report_event.cpp | 16 ++++++++++++---- coreneuron/io/reports/report_handler.cpp | 17 ++++++++++++++--- coreneuron/io/reports/report_handler.hpp | 3 ++- 3 files changed, 28 insertions(+), 8 deletions(-) diff --git a/coreneuron/io/reports/report_event.cpp b/coreneuron/io/reports/report_event.cpp index 9deaff30c..5a943bf02 100644 --- a/coreneuron/io/reports/report_event.cpp +++ b/coreneuron/io/reports/report_event.cpp @@ -80,7 +80,7 @@ void ReportEvent::lfp_calc(NrnThread* nt) { if (step > 0 && (static_cast(step) % reporting_period) == 0) { auto* mapinfo = static_cast(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; @@ -94,7 +94,13 @@ void ReportEvent::lfp_calc(NrnThread* nt) { if(std::isnan(factor)) { factor = 0.0; } - sum += fast_imem_rhs[segment_id] * factor; + 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; @@ -107,8 +113,10 @@ void ReportEvent::deliver(double t, NetCvode* nc, NrnThread* nt) { /* reportinglib is not thread safe */ #pragma omp critical { - summation_alu(nt); - if (report_type == ReportType::LFPReport) { + 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 diff --git a/coreneuron/io/reports/report_handler.cpp b/coreneuron/io/reports/report_handler.cpp index 1b69cd6ff..53efb137d 100644 --- a/coreneuron/io/reports/report_handler.cpp +++ b/coreneuron/io/reports/report_handler.cpp @@ -62,7 +62,7 @@ void ReportHandler::create_report(double dt, double tstop, double delay) { 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()); + 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); @@ -353,14 +353,25 @@ VarsToReport ReportHandler::get_synapse_vars_to_report( VarsToReport ReportHandler::get_lfp_vars_to_report(const NrnThread& nt, ReportConfiguration& report, - double* report_variable) const { + double* report_variable, + const std::vector& 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 to_report; double* variable = report_variable + i; to_report.push_back(VarWithMapping(i, variable)); diff --git a/coreneuron/io/reports/report_handler.hpp b/coreneuron/io/reports/report_handler.hpp index db62bf993..746edbaee 100644 --- a/coreneuron/io/reports/report_handler.hpp +++ b/coreneuron/io/reports/report_handler.hpp @@ -46,7 +46,8 @@ class ReportHandler { const std::vector& nodes_to_gids) const; VarsToReport get_lfp_vars_to_report(const NrnThread& nt, ReportConfiguration& report, - double* report_variable) const; + double* report_variable, + const std::vector& nodes_to_gids) const; std::vector map_gids(const NrnThread& nt) const; #endif // defined(ENABLE_BIN_REPORTS) || defined(ENABLE_SONATA_REPORTS) protected: From 8936692dadbc4a4ddc1de1c63f0b9d9fb507df71 Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso Date: Wed, 6 Apr 2022 11:22:18 +0200 Subject: [PATCH 6/6] Clang format --- coreneuron/io/reports/report_event.cpp | 29 ++++++++++++------------ coreneuron/io/reports/report_handler.cpp | 3 ++- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/coreneuron/io/reports/report_event.cpp b/coreneuron/io/reports/report_event.cpp index 5a943bf02..6189c4526 100644 --- a/coreneuron/io/reports/report_event.cpp +++ b/coreneuron/io/reports/report_event.cpp @@ -89,19 +89,19 @@ void ReportEvent::lfp_calc(NrnThread* nt) { 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++; + 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; } @@ -115,8 +115,7 @@ void ReportEvent::deliver(double t, NetCvode* nc, NrnThread* nt) { { if (report_type == ReportType::SummationReport) { summation_alu(nt); - } - else if (report_type == ReportType::LFPReport) { + } else if (report_type == ReportType::LFPReport) { lfp_calc(nt); } // each thread needs to know its own step diff --git a/coreneuron/io/reports/report_handler.cpp b/coreneuron/io/reports/report_handler.cpp index 53efb137d..e720f331e 100644 --- a/coreneuron/io/reports/report_handler.cpp +++ b/coreneuron/io/reports/report_handler.cpp @@ -62,7 +62,8 @@ void ReportHandler::create_report(double dt, double tstop, double delay) { 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); + 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);