diff --git a/CMakeLists.txt b/CMakeLists.txt index 183639480..0a02ec93d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,7 +10,7 @@ cmake_minimum_required(VERSION 3.15 FATAL_ERROR) # standalone. project( coreneuron - VERSION 9.0.0 + VERSION 8.2.1 LANGUAGES CXX) # ~~~ diff --git a/coreneuron/io/nrn_filehandler.hpp b/coreneuron/io/nrn_filehandler.hpp index 65cab4b97..9476a6c53 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,27 +111,45 @@ 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, CellMapping* cmap) { int nsec, nseg, n_scan; + size_t total_lfp_factors; + int num_electrodes; char line_buf[max_line_length], name[max_line_length]; F.getline(line_buf, sizeof(line_buf)); - n_scan = sscanf(line_buf, "%s %d %d", name, &nsec, &nseg); + n_scan = sscanf( + line_buf, "%s %d %d %zd %d", name, &nsec, &nseg, &total_lfp_factors, &num_electrodes); - nrn_assert(n_scan == 3); + nrn_assert(n_scan == 5); mapinfo->name = std::string(name); if (nseg) { std::vector sec, seg; + std::vector lfp_factors; + sec.reserve(nseg); seg.reserve(nseg); + lfp_factors.reserve(total_lfp_factors); read_array(&sec[0], nseg); read_array(&seg[0], nseg); + if (total_lfp_factors > 0) { + read_array(&lfp_factors[0], total_lfp_factors); + } + int factor_offset = 0; for (int i = 0; i < nseg; i++) { mapinfo->add_segment(sec[i], seg[i]); + ntmapping->add_segment_id(seg[i]); + int factor_offset = i * num_electrodes; + if (total_lfp_factors > 0) { + std::vector segment_factors(lfp_factors.begin() + factor_offset, + lfp_factors.begin() + factor_offset + + num_electrodes); + cmap->add_segment_lfp_factor(seg[i], segment_factors); + } } } return nseg; diff --git a/coreneuron/io/nrn_setup.cpp b/coreneuron/io/nrn_setup.cpp index 703e853d8..19d5b5297 100644 --- a/coreneuron/io/nrn_setup.cpp +++ b/coreneuron/io/nrn_setup.cpp @@ -959,7 +959,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); } diff --git a/coreneuron/io/nrnsection_mapping.hpp b/coreneuron/io/nrnsection_mapping.hpp index d7524fc42..11e79d086 100644 --- a/coreneuron/io/nrnsection_mapping.hpp +++ b/coreneuron/io/nrnsection_mapping.hpp @@ -14,6 +14,7 @@ #include #include #include +#include namespace coreneuron { @@ -73,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) {} @@ -96,6 +100,15 @@ struct CellMapping { }); } + /** @brief return the number of electrodes in the lfp_factors map **/ + int num_electrodes() const { + int num_electrodes = 0; + if (!lfp_factors.empty()) { + num_electrodes = lfp_factors.begin()->second.size(); + } + return num_electrodes; + } + /** @brief number of section lists */ size_t size() const noexcept { return secmapvec.size(); @@ -137,6 +150,11 @@ struct CellMapping { return count; } + /** @brief add the lfp electrode factors of a segment_id */ + void add_segment_lfp_factor(const int segment_id, std::vector factors) { + lfp_factors.insert({segment_id, factors}); + } + ~CellMapping() { for (size_t i = 0; i < secmapvec.size(); i++) { delete secmapvec[i]; @@ -153,6 +171,11 @@ struct NrnThreadMappingInfo { /** list of cells mapping */ std::vector mappingvec; + /** list of segment ids */ + std::vector segment_ids; + + std::vector _lfp; + /** @brief number of cells */ size_t size() const { return mappingvec.size(); @@ -181,5 +204,19 @@ 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 Resize the lfp vector */ + void prepare_lfp() { + size_t lfp_size = 0; + for (const auto& mapping: mappingvec) { + lfp_size += mapping->num_electrodes(); + } + _lfp.resize(lfp_size); + } }; } // namespace coreneuron 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..6182efd0c 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 @@ -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,53 @@ 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* 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 num_electrodes = cell_mapping->num_electrodes(); + std::vector lfp_values(num_electrodes, 0.0); + for (const auto& kv: cell_mapping->lfp_factors) { + int segment_id = kv.first; + std::vector factors = kv.second; + int electrode_id = 0; + for (auto& factor: factors) { + 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; + } + lfp_values[electrode_id] += (fast_imem_rhs[segment_id] + iclamp) * factor; + electrode_id++; + } + } + for (int i = 0; i < to_report.size(); i++) { + *(to_report[i].var_value) = lfp_values[i]; + } + } + } +} + /** 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()); diff --git a/coreneuron/io/reports/report_event.hpp b/coreneuron/io/reports/report_event.hpp index 20d325614..0f1a07358 100644 --- a/coreneuron/io/reports/report_event.hpp +++ b/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 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 84341e7a4..760ed0749 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,14 @@ 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->prepare_lfp(); + 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 +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)); } @@ -341,6 +351,55 @@ 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& nodes_to_gids) const { + 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); + } + auto& summation_report = nt.summation_report_handler_->summation_reports_[report.output_path]; + VarsToReport vars_to_report; + off_t offset_lfp = 0; + 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]; + if (ml) { + 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)); + } + } + } + const auto& cell_mapping = mapinfo->get_cell_mapping(gid); + if (cell_mapping == nullptr) { + std::cerr << "[LFP] Error : Compartment mapping information is missing for gid " << gid + << '\n'; + nrn_abort(1); + } + std::vector to_report; + int num_electrodes = cell_mapping->num_electrodes(); + for (int electrode_id = 0; electrode_id < num_electrodes; electrode_id++) { + to_report.emplace_back(VarWithMapping(electrode_id, report_variable + offset_lfp)); + offset_lfp++; + } + if (!to_report.empty()) { + 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..746edbaee 100644 --- a/coreneuron/io/reports/report_handler.hpp +++ b/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& nodes_to_gids) const; + VarsToReport get_lfp_vars_to_report(const NrnThread& nt, + ReportConfiguration& report, + 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: diff --git a/tests/integration/ring/0_3.dat b/tests/integration/ring/0_3.dat index b79f96e2e..7e68f2cdb 100644 Binary files a/tests/integration/ring/0_3.dat and b/tests/integration/ring/0_3.dat differ diff --git a/tests/integration/ring/1_3.dat b/tests/integration/ring/1_3.dat index ae6b27a25..a4579fd1a 100644 Binary files a/tests/integration/ring/1_3.dat and b/tests/integration/ring/1_3.dat differ