Add summation reports (#514)

* Support i_membrane in summation
* Clang format nitpicks
* Error when trying to use summation report on soma target

Author: jorblancoa

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<SummationReportMapping>();
 }
 
 static size_t memb_list_size(NrnThreadMembList* tml) {

coreneuron/io/reports/nrnreport.hpp

@@ -17,44 +17,82 @@
 #include <string>
 #include <vector>
 #include <set>
+#include <unordered_map>
 
 #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<double> summation_ = {};
+    // Map containing the pointers of the currents and its scaling factor for every segment_id
+    std::unordered_map<int, std::vector<std::pair<double*, int>>> currents_;
+};
+
+struct SummationReportMapping {
+    // Map containing an SummationReport object per report
+    std::unordered_map<std::string, SummationReport> 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"
 
 /// 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<int> 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<std::string> mech_names;  // mechanism names
+    std::vector<std::string> var_names;   // variable names
+    std::vector<int> 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<int> target;                 // list of gids for this report
 };
 
 void setup_report_engine(double dt_report, double mindelay);

coreneuron/io/reports/report_configuration_parser.cpp

@@ -23,52 +23,54 @@
 
 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<std::string> 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<ReportConfiguration> create_report_configurations(const std::string& conf_file,
                                                               const std::string& output_dir,
                                                               std::string& spikes_population_name) {
     std::vector<ReportConfiguration> 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<TargetType>(target);
         std::transform(report.type_str.begin(),
                        report.type_str.end(),
                        report.type_str.begin(),
@@ -79,7 +81,7 @@ std::vector<ReportConfiguration> create_report_configurations(const std::string&
                 nrn_use_fast_imem = true;
                 report.type = IMembraneReport;
             } else {
-                switch (static_cast<TargetType>(target_type)) {
+                switch (report.target_type) {
                     case TargetType::Soma:
                         report.type = SomaReport;
                         break;
@@ -123,11 +125,13 @@ std::vector<ReportConfiguration> 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) {

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,25 +23,48 @@ 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<int>(report_dt / dt);
     gids_to_report.reserve(filtered_gids.size());
     for (const auto& gid: filtered_gids) {
         gids_to_report.push_back(gid.first);
     }
     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<int>(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());

coreneuron/io/reports/report_event.hpp

@@ -32,16 +32,23 @@ using VarsToReport = std::unordered_map<int, std::vector<VarWithMapping>>;
 
 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<int> gids_to_report;
     double tstart;
 };