From cbeb991328a97026d61e2d6f85a8b781b917a514 Mon Sep 17 00:00:00 2001 From: Jorge Blanco Alonso Date: Tue, 5 Apr 2022 18:32:03 +0200 Subject: [PATCH] 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: