diff --git a/coreneuron/io/nrn_setup.cpp b/coreneuron/io/nrn_setup.cpp index 8fa8afa19..4e1a91a88 100644 --- a/coreneuron/io/nrn_setup.cpp +++ b/coreneuron/io/nrn_setup.cpp @@ -36,6 +36,7 @@ #include "coreneuron/io/mech_report.h" #include "coreneuron/apps/corenrn_parameters.hpp" #include "coreneuron/io/nrn_setup.hpp" +#include "coreneuron/io/reports/nrnreport.hpp" // callbacks into nrn/src/nrniv/nrnbbcore_write.cpp #include "coreneuron/sim/fast_imem.hpp" @@ -951,6 +952,7 @@ void read_phase3(NrnThread& nt, UserParams& userParams) { // set pointer in NrnThread nt.mapping = (void*) ntmapping; + nt.summation_report_handler_ = std::make_unique(); } static size_t memb_list_size(NrnThreadMembList* tml) { diff --git a/coreneuron/io/reports/nrnreport.hpp b/coreneuron/io/reports/nrnreport.hpp index f08209778..89e15fba2 100644 --- a/coreneuron/io/reports/nrnreport.hpp +++ b/coreneuron/io/reports/nrnreport.hpp @@ -17,11 +17,25 @@ #include #include #include +#include #define REPORT_MAX_NAME_LEN 256 #define REPORT_MAX_FILEPATH_LEN 4096 namespace coreneuron { + +struct SummationReport { + // Contains the values of the summation with index == segment_id + std::vector summation_ = {}; + // Map containing the pointers of the currents and its scaling factor for every segment_id + std::unordered_map>> currents_; +}; + +struct SummationReportMapping { + // Map containing an SummationReport object per report + std::unordered_map summation_reports_; +}; + // name of the variable in mod file that is used to indicate which synapse // is enabled or disable for reporting #define SELECTED_VAR_MOD_NAME "selected_for_report" @@ -29,32 +43,56 @@ namespace coreneuron { /// name of the variable in mod file used for setting synapse id #define SYNAPSE_ID_MOD_NAME "synapseID" +/* + * Defines the type of target, as per the following syntax: + * 0=Compartment, 1=Cell/Soma, Section { 2=Axon, 3=Dendrite, 4=Apical } + * The "Comp" variations are compartment-based (all segments, not middle only) + */ +enum class TargetType { + Compartment = 0, + Soma = 1, + Axon = 2, + Dendrite = 3, + Apical = 4, + AxonComp = 5, + DendriteComp = 6, + ApicalComp = 7, +}; + // enumerate that defines the type of target report requested -enum ReportType { SomaReport, CompartmentReport, SynapseReport, IMembraneReport, SectionReport }; +enum ReportType { + SomaReport, + CompartmentReport, + SynapseReport, + IMembraneReport, + SectionReport, + SummationReport +}; // enumerate that defines the section type for a Section report enum SectionType { Axon, Dendrite, Apical }; struct ReportConfiguration { - std::string name; // name of the report - std::string output_path; // full path of the report - std::string target_name; // target of the report - std::string mech_name; // mechanism name - std::string var_name; // variable name - std::string unit; // unit of the report - std::string format; // format of the report (Bin, hdf5, SONATA) - std::string type_str; // type of report string - std::string population_name; // population name of the report - ReportType type; // type of the report - SectionType section_type; // type of section report - bool section_all_compartments; // flag for section report (all values) - int mech_id; // mechanism - double report_dt; // reporting timestep - double start; // start time of report - double stop; // stop time of report - int num_gids; // total number of gids - int buffer_size; // hint on buffer size used for this report - std::set target; // list of gids for this report + std::string name; // name of the report + std::string output_path; // full path of the report + std::string target_name; // target of the report + std::vector mech_names; // mechanism names + std::vector var_names; // variable names + std::vector mech_ids; // mechanisms + std::string unit; // unit of the report + std::string format; // format of the report (Bin, hdf5, SONATA) + std::string type_str; // type of report string + std::string population_name; // population name of the report + TargetType target_type; // type of the target + ReportType type; // type of the report + SectionType section_type; // type of section report + bool section_all_compartments; // flag for section report (all values) + double report_dt; // reporting timestep + double start; // start time of report + double stop; // stop time of report + int num_gids; // total number of gids + int buffer_size; // hint on buffer size used for this report + std::set target; // list of gids for this report }; void setup_report_engine(double dt_report, double mindelay); diff --git a/coreneuron/io/reports/report_configuration_parser.cpp b/coreneuron/io/reports/report_configuration_parser.cpp index c8b56852d..a0367264a 100644 --- a/coreneuron/io/reports/report_configuration_parser.cpp +++ b/coreneuron/io/reports/report_configuration_parser.cpp @@ -23,30 +23,33 @@ namespace coreneuron { -/* - * Defines the type of target, as per the following syntax: - * 0=Compartment, 1=Cell/Soma, Section { 2=Axon, 3=Dendrite, 4=Apical } - * The "Comp" variations are compartment-based (all segments, not middle only) - */ -enum class TargetType { - Compartment = 0, - Soma = 1, - Axon = 2, - Dendrite = 3, - Apical = 4, - AxonComp = 5, - DendriteComp = 6, - ApicalComp = 7, -}; /* - * Split filter string ("mech.var_name") into mech_id and var_name + * Split filter comma separated strings ("mech.var_name") into mech_name and var_name */ void parse_filter_string(const std::string& filter, ReportConfiguration& config) { - std::istringstream iss(filter); - std::string token; - std::getline(iss, config.mech_name, '.'); - std::getline(iss, config.var_name, '.'); + std::vector mechanisms; + std::stringstream ss(filter); + std::string mechanism; + // Multiple report variables are separated by `,` + while (getline(ss, mechanism, ',')) { + mechanisms.push_back(mechanism); + + // Split mechanism name and corresponding reporting variable + std::string mech_name; + std::string var_name; + std::istringstream iss(mechanism); + std::getline(iss, mech_name, '.'); + std::getline(iss, var_name, '.'); + if (var_name.empty()) { + var_name = "i"; + } + config.mech_names.emplace_back(mech_name); + config.var_names.emplace_back(var_name); + if (mech_name == "i_membrane") { + nrn_use_fast_imem = true; + } + } } std::vector create_report_configurations(const std::string& conf_file, @@ -54,21 +57,20 @@ std::vector create_report_configurations(const std::string& std::string& spikes_population_name) { std::vector reports; std::string report_on; - int target_type; + int target; std::ifstream report_conf(conf_file); int num_reports = 0; report_conf >> num_reports; for (int i = 0; i < num_reports; i++) { ReportConfiguration report; - // mechansim id registered in coreneuron - report.mech_id = -1; report.buffer_size = 4; // default size to 4 Mb report_conf >> report.name >> report.target_name >> report.type_str >> report_on >> - report.unit >> report.format >> target_type >> report.report_dt >> report.start >> + report.unit >> report.format >> target >> report.report_dt >> report.start >> report.stop >> report.num_gids >> report.buffer_size >> report.population_name; + report.target_type = static_cast(target); std::transform(report.type_str.begin(), report.type_str.end(), report.type_str.begin(), @@ -79,7 +81,7 @@ std::vector create_report_configurations(const std::string& nrn_use_fast_imem = true; report.type = IMembraneReport; } else { - switch (static_cast(target_type)) { + switch (report.target_type) { case TargetType::Soma: report.type = SomaReport; break; @@ -123,11 +125,13 @@ std::vector create_report_configurations(const std::string& } } else if (report.type_str == "synapse") { report.type = SynapseReport; + } else if (report.type_str == "summation") { + report.type = SummationReport; } else { std::cerr << "Report error: unsupported type " << report.type_str << std::endl; nrn_abort(1); } - if (report.type == SynapseReport) { + if (report.type == SynapseReport || report.type == SummationReport) { parse_filter_string(report_on, report); } if (report.num_gids) { diff --git a/coreneuron/io/reports/report_event.cpp b/coreneuron/io/reports/report_event.cpp index 36b208a5c..6bc1b9364 100644 --- a/coreneuron/io/reports/report_event.cpp +++ b/coreneuron/io/reports/report_event.cpp @@ -8,6 +8,7 @@ #include "report_event.hpp" #include "coreneuron/sim/multicore.hpp" +#include "coreneuron/io/reports/nrnreport.hpp" #include "coreneuron/utils/nrn_assert.h" #ifdef ENABLE_BIN_REPORTS #include "reportinglib/Records.h" @@ -22,13 +23,16 @@ namespace coreneuron { ReportEvent::ReportEvent(double dt, double tstart, const VarsToReport& filtered_gids, - const char* name) + const char* name, + double report_dt) : dt(dt) , tstart(tstart) - , report_path(name) { + , report_path(name) + , report_dt(report_dt) { VarsToReport::iterator it; nrn_assert(filtered_gids.size()); step = tstart / dt; + reporting_period = static_cast(report_dt / dt); gids_to_report.reserve(filtered_gids.size()); for (const auto& gid: filtered_gids) { gids_to_report.push_back(gid.first); @@ -36,11 +40,31 @@ ReportEvent::ReportEvent(double dt, std::sort(gids_to_report.begin(), gids_to_report.end()); } +void ReportEvent::summation_alu(NrnThread* nt) { + // Sum currents only on reporting steps + if (static_cast(step) % reporting_period == 0) { + auto& summation_report = nt->summation_report_handler_->summation_reports_[report_path]; + // Add currents of all variables in each segment + double sum = 0.0; + for (const auto& kv: summation_report.currents_) { + int segment_id = kv.first; + for (const auto& value: kv.second) { + double current_value = *value.first; + int scale = value.second; + sum += current_value * scale; + } + summation_report.summation_[segment_id] = sum; + sum = 0.0; + } + } +} + /** 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); // 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 0d4f18917..8c5da5b95 100644 --- a/coreneuron/io/reports/report_event.hpp +++ b/coreneuron/io/reports/report_event.hpp @@ -32,16 +32,23 @@ using VarsToReport = std::unordered_map>; class ReportEvent: public DiscreteEvent { public: - ReportEvent(double dt, double tstart, const VarsToReport& filtered_gids, const char* name); + ReportEvent(double dt, + double tstart, + const VarsToReport& filtered_gids, + const char* name, + double report_dt); /** 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); private: double dt; double step; std::string report_path; + double report_dt; + int reporting_period; std::vector gids_to_report; double tstart; }; diff --git a/coreneuron/io/reports/report_handler.cpp b/coreneuron/io/reports/report_handler.cpp index 14c33fecb..d7ceb6498 100644 --- a/coreneuron/io/reports/report_handler.cpp +++ b/coreneuron/io/reports/report_handler.cpp @@ -19,9 +19,11 @@ void ReportHandler::create_report(double dt, double tstop, double delay) { } m_report_config.stop = std::min(m_report_config.stop, tstop); - m_report_config.mech_id = nrn_get_mechtype(m_report_config.mech_name.data()); - if (m_report_config.type == SynapseReport && m_report_config.mech_id == -1) { - std::cerr << "[ERROR] mechanism to report: " << m_report_config.mech_name + for (const auto& mech: m_report_config.mech_names) { + m_report_config.mech_ids.emplace_back(nrn_get_mechtype(mech.data())); + } + if (m_report_config.type == SynapseReport && m_report_config.mech_ids.empty()) { + std::cerr << "[ERROR] mechanism to report: " << m_report_config.mech_names[0] << " is not mapped in this simulation, cannot report on it \n"; nrn_abort(1); } @@ -57,6 +59,13 @@ void ReportHandler::create_report(double dt, double tstop, double delay) { m_report_config.section_all_compartments); register_compartment_report(nt, m_report_config, vars_to_report); break; + case SummationReport: + vars_to_report = get_summation_vars_to_report(nt, + m_report_config.target, + m_report_config, + nodes_to_gid); + register_custom_report(nt, m_report_config, vars_to_report); + break; default: vars_to_report = get_custom_vars_to_report(nt, m_report_config, nodes_to_gid); register_custom_report(nt, m_report_config, vars_to_report); @@ -65,7 +74,8 @@ void ReportHandler::create_report(double dt, double tstop, double delay) { auto report_event = std::make_unique(dt, t, vars_to_report, - m_report_config.output_path.data()); + m_report_config.output_path.data(), + m_report_config.report_dt); report_event->send(t, net_cvode_instance, &nt); m_report_events.push_back(std::move(report_event)); } @@ -253,6 +263,93 @@ VarsToReport ReportHandler::get_section_vars_to_report(const NrnThread& nt, return vars_to_report; } +VarsToReport ReportHandler::get_summation_vars_to_report( + const NrnThread& nt, + const std::set& target, + ReportConfiguration& report, + const std::vector& nodes_to_gids) const { + VarsToReport vars_to_report; + const auto* mapinfo = static_cast(nt.mapping); + auto& summation_report = nt.summation_report_handler_->summation_reports_[report.output_path]; + if (!mapinfo) { + std::cerr << "[COMPARTMENTS] Error : mapping information is missing for a Cell group " + << nt.ncell << '\n'; + nrn_abort(1); + } + if (report.target_type == TargetType::Soma) { + std::cerr << "[SUMMATION] Error with report: '" << report.name + << "' Soma target not supported with summation reports" << std::endl; + 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; + } + bool has_imembrane = false; + // In case we need convertion of units + int scale = 1; + for (auto i = 0; i < report.mech_ids.size(); ++i) { + auto mech_id = report.mech_ids[i]; + auto var_name = report.var_names[i]; + auto mech_name = report.mech_names[i]; + if (mech_name != "i_membrane") { + // need special handling for Clamp processes to flip the current value + if (mech_name == "IClamp" || mech_name == "SEClamp") { + scale = -1; + } + Memb_list* ml = nt._ml_list[mech_id]; + if (!ml) { + continue; + } + + for (int j = 0; j < ml->nodecount; j++) { + auto segment_id = ml->nodeindices[j]; + if ((nodes_to_gids[ml->nodeindices[j]] == gid)) { + double* var_value = + get_var_location_from_var_name(mech_id, var_name.data(), ml, j); + summation_report.currents_[segment_id].push_back( + std::make_pair(var_value, scale)); + } + } + } else { + has_imembrane = true; + } + } + if (target.find(gid) != target.end()) { + const auto& cell_mapping = mapinfo->get_cell_mapping(gid); + if (cell_mapping == nullptr) { + std::cerr + << "[SUMMATION] Error : Compartment mapping information is missing for gid " + << gid << '\n'; + nrn_abort(1); + } + std::vector to_report; + to_report.reserve(cell_mapping->size()); + summation_report.summation_.resize(nt.end); + const auto& section_mapping = cell_mapping->secmapvec; + for (const auto& sections: section_mapping) { + for (auto& section: sections->secmap) { + // compartment_id + int section_id = section.first; + auto& segment_ids = section.second; + for (const auto& segment_id: segment_ids) { + /** corresponding voltage in coreneuron voltage array */ + if (has_imembrane) { + summation_report.currents_[segment_id].push_back( + std::make_pair(nt.nrn_fast_imem->nrn_sav_rhs + segment_id, 1)); + } + double* variable = summation_report.summation_.data() + segment_id; + to_report.emplace_back(VarWithMapping(section_id, variable)); + } + } + } + vars_to_report[gid] = to_report; + } + } + return vars_to_report; +} + VarsToReport ReportHandler::get_custom_vars_to_report(const NrnThread& nt, ReportConfiguration& report, const std::vector& nodes_to_gids) const { @@ -262,7 +359,11 @@ VarsToReport ReportHandler::get_custom_vars_to_report(const NrnThread& nt, if (report.target.find(gid) == report.target.end()) { continue; } - Memb_list* ml = nt._ml_list[report.mech_id]; + // There can only be 1 mechanism + nrn_assert(report.mech_ids.size() == 1); + auto mech_id = report.mech_ids[0]; + auto var_name = report.var_names[0]; + Memb_list* ml = nt._ml_list[mech_id]; if (!ml) { continue; } @@ -271,7 +372,7 @@ VarsToReport ReportHandler::get_custom_vars_to_report(const NrnThread& nt, for (int j = 0; j < ml->nodecount; j++) { double* is_selected = - get_var_location_from_var_name(report.mech_id, SELECTED_VAR_MOD_NAME, ml, j); + get_var_location_from_var_name(mech_id, SELECTED_VAR_MOD_NAME, ml, j); bool report_variable = false; /// if there is no variable in mod file then report on every compartment @@ -282,10 +383,9 @@ VarsToReport ReportHandler::get_custom_vars_to_report(const NrnThread& nt, report_variable = *is_selected != 0.; } if ((nodes_to_gids[ml->nodeindices[j]] == gid) && report_variable) { - double* var_value = - get_var_location_from_var_name(report.mech_id, report.var_name.data(), ml, j); + double* var_value = get_var_location_from_var_name(mech_id, var_name.data(), ml, j); double* synapse_id = - get_var_location_from_var_name(report.mech_id, SYNAPSE_ID_MOD_NAME, ml, j); + get_var_location_from_var_name(mech_id, SYNAPSE_ID_MOD_NAME, ml, j); nrn_assert(synapse_id && var_value); to_report.emplace_back(static_cast(*synapse_id), var_value); } diff --git a/coreneuron/io/reports/report_handler.hpp b/coreneuron/io/reports/report_handler.hpp index cab03c981..403511eb1 100644 --- a/coreneuron/io/reports/report_handler.hpp +++ b/coreneuron/io/reports/report_handler.hpp @@ -46,6 +46,10 @@ class ReportHandler { double* report_variable, SectionType section_type, bool all_compartments) const; + VarsToReport get_summation_vars_to_report(const NrnThread& nt, + const std::set& target, + ReportConfiguration& report, + const std::vector& nodes_to_gids) const; VarsToReport get_custom_vars_to_report(const NrnThread& nt, ReportConfiguration& report, const std::vector& nodes_to_gids) const; diff --git a/coreneuron/sim/multicore.hpp b/coreneuron/sim/multicore.hpp index 820c7d79a..036ac3432 100644 --- a/coreneuron/sim/multicore.hpp +++ b/coreneuron/sim/multicore.hpp @@ -12,7 +12,9 @@ #include "coreneuron/mechanism/membfunc.hpp" #include "coreneuron/utils/memory.h" #include "coreneuron/mpi/nrnmpi.h" +#include "coreneuron/io/reports/nrnreport.hpp" #include +#include namespace coreneuron { class NetCon; @@ -129,8 +131,10 @@ struct NrnThread: public MemoryManaged { int _net_send_buffer_cnt = 0; int* _net_send_buffer = nullptr; - int* _watch_types = nullptr; /* nullptr or 0 terminated array of integers */ - void* mapping = nullptr; /* section to segment mapping information */ + int* _watch_types = nullptr; /* nullptr or 0 terminated array of integers */ + void* mapping = nullptr; /* section to segment mapping information */ + std::unique_ptr summation_report_handler_; /* report to ALU (values of + the current summation */ TrajectoryRequests* trajec_requests = nullptr; /* per time step values returned to NEURON */ /* Needed in case there are FOR_NETCON statements in use. */