LFP Calculator

coreneuron/apps/main1.cpp

@@ -490,7 +490,7 @@ extern "C" int run_solve_core(int argc, char** argv) {
 
     // clang-format off
 
-    #pragma acc update device(celsius, secondorder) if (compute_gpu)
+    #pragma acc update device(celsius, secondorder, pi) if (compute_gpu)
     // clang-format on
     {
         double v = corenrn_param.voltage;

coreneuron/io/global_vars.cpp

@@ -48,6 +48,7 @@ void set_globals(const char* path, bool cli_global_seed, int cli_global_seed_val
     (*n2v)["celsius"] = PSD(0, &celsius);
     (*n2v)["dt"] = PSD(0, &dt);
     (*n2v)["t"] = PSD(0, &t);
+    (*n2v)["PI"] = PSD(0, &pi);
 
     if (corenrn_embedded) {  // CoreNEURON embedded, get info direct from NEURON
 

coreneuron/io/lfp.cpp

@@ -0,0 +1,141 @@
+#include "coreneuron/io/lfp.hpp"
+
+#include <cmath>
+#include <limits>
+#include <sstream>
+
+
+namespace coreneuron {
+
+// extern variables require acc declare
+#pragma acc declare create(pi)
+
+namespace lfputils {
+
+double line_source_lfp_factor(const Point3D& e_pos,
+                              const Point3D& seg_0,
+                              const Point3D& seg_1,
+                              const double radius,
+                              const double f) {
+    nrn_assert(radius >= 0.0);
+    Point3D dx = paxpy(seg_1, -1.0, seg_0);
+    Point3D de = paxpy(e_pos, -1.0, seg_0);
+    double dx2(dot(dx, dx));
+    double dxn(std::sqrt(dx2));
+    if (dxn < std::numeric_limits<double>::epsilon()) {
+        return point_source_lfp_factor(e_pos, seg_0, radius, f);
+    }
+    double de2(dot(de, de));
+    double mu(dot(dx, de) / dx2);
+    Point3D de_star(paxpy(de, -mu, dx));
+    double de_star2(dot(de_star, de_star));
+    double q2(de_star2 / dx2);
+
+    double delta(mu * mu - (de2 - radius * radius) / dx2);
+    double one_m_mu(1.0 - mu);
+    auto log_integral = [&q2, &dxn](double a, double b) {
+        if (q2 < std::numeric_limits<double>::epsilon()) {
+            if (a * b <= 0) {
+                std::ostringstream s;
+                s << "Log integral: invalid arguments " << b << " " << a
+                  << ". Likely electrode exactly on the segment and "
+                  << "no flooring is present.";
+                throw std::invalid_argument(s.str());
+            }
+            return std::abs(std::log(b / a)) / dxn;
+        } else {
+            return std::log((b + std::sqrt(b * b + q2)) / (a + std::sqrt(a * a + q2))) / dxn;
+        }
+    };
+    if (delta <= 0.0) {
+        return f * log_integral(-mu, one_m_mu);
+    } else {
+        double sqr_delta(std::sqrt(delta));
+        double d1(mu - sqr_delta);
+        double d2(mu + sqr_delta);
+        double parts = 0.0;
+        if (d1 > 0.0) {
+            double b(std::min(d1, 1.0) - mu);
+            parts += log_integral(-mu, b);
+        }
+        if (d2 < 1.0) {
+            double b(std::max(d2, 0.0) - mu);
+            parts += log_integral(b, one_m_mu);
+        };
+        // complement
+        double maxd1_0(std::max(d1, 0.0)), mind2_1(std::min(d2, 1.0));
+        if (maxd1_0 < mind2_1) {
+            parts += 1.0 / radius * (mind2_1 - maxd1_0);
+        }
+        return f * parts;
+    };
+}
+}  // namespace lfputils
+
+using namespace lfputils;
+
+template <LFPCalculatorType Type, typename SegmentIdTy>
+LFPCalculator<Type, SegmentIdTy>::LFPCalculator(const Point3Ds& seg_start,
+                                                const Point3Ds& seg_end,
+                                                const std::vector<double>& radius,
+                                                const std::vector<SegmentIdTy>& segment_ids,
+                                                const Point3Ds& electrodes,
+                                                double extra_cellular_conductivity)
+    : segment_ids_(segment_ids) {
+    if (seg_start.size() != seg_end.size()) {
+        throw std::invalid_argument("Different number of segment starts and ends.");
+    }
+    if (seg_start.size() != radius.size()) {
+        throw std::invalid_argument("Different number of segments and radii.");
+    }
+    double f(1.0 / (extra_cellular_conductivity * 4.0 * pi));
+
+    m.resize(electrodes.size());
+    for (size_t k = 0; k < electrodes.size(); ++k) {
+        auto& ms = m[k];
+        ms.resize(seg_start.size());
+        for (size_t l = 0; l < seg_start.size(); l++) {
+            ms[l] = getFactor(electrodes[k], seg_start[l], seg_end[l], radius[l], f);
+        }
+    }
+}
+
+template <LFPCalculatorType Type, typename SegmentIdTy>
+template <typename Vector>
+inline void LFPCalculator<Type, SegmentIdTy>::lfp(const Vector& membrane_current) {
+    std::vector<double> res(m.size());
+    for (size_t k = 0; k < m.size(); ++k) {
+        res[k] = 0.0;
+        auto& ms = m[k];
+        for (size_t l = 0; l < ms.size(); l++) {
+            res[k] += ms[l] * membrane_current[segment_ids_[l]];
+        }
+    }
+#if NRNMPI
+    lfp_values_.resize(res.size());
+    int mpi_sum{1};
+    nrnmpi_dbl_allreduce_vec(res.data(), lfp_values_.data(), res.size(), mpi_sum);
+#else
+    std::swap(res, lfp_values_);
+#endif
+}
+
+
+template LFPCalculator<LineSource>::LFPCalculator(const lfputils::Point3Ds& seg_start,
+                                                  const lfputils::Point3Ds& seg_end,
+                                                  const std::vector<double>& radius,
+                                                  const std::vector<int>& segment_ids,
+                                                  const lfputils::Point3Ds& electrodes,
+                                                  double extra_cellular_conductivity);
+template LFPCalculator<PointSource>::LFPCalculator(const lfputils::Point3Ds& seg_start,
+                                                   const lfputils::Point3Ds& seg_end,
+                                                   const std::vector<double>& radius,
+                                                   const std::vector<int>& segment_ids,
+                                                   const lfputils::Point3Ds& electrodes,
+                                                   double extra_cellular_conductivity);
+template void LFPCalculator<LineSource>::lfp(const DoublePtr& membrane_current);
+template void LFPCalculator<PointSource>::lfp(const DoublePtr& membrane_current);
+template void LFPCalculator<LineSource>::lfp(const std::vector<double>& membrane_current);
+template void LFPCalculator<PointSource>::lfp(const std::vector<double>& membrane_current);
+
+}  // namespace coreneuron

coreneuron/io/lfp.hpp

@@ -0,0 +1,142 @@
+#pragma once
+
+#include <array>
+#include <vector>
+
+#include "coreneuron/mpi/nrnmpidec.h"
+#include "coreneuron/nrnconf.h"
+#include "coreneuron/utils/nrn_assert.h"
+
+namespace coreneuron {
+
+namespace lfputils {
+
+using Point3D = std::array<double, 3>;
+using Point3Ds = std::vector<Point3D>;
+using DoublePtr = double*;
+
+inline double dot(const Point3D& p1, const Point3D& p2) {
+    return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2];
+}
+
+inline double norm(const Point3D& p1) {
+    return std::sqrt(dot(p1, p1));
+}
+
+inline Point3D barycenter(const Point3D& p1, const Point3D& p2) {
+    return {0.5 * (p1[0] + p2[0]), 0.5 * (p1[1] + p2[1]), 0.5 * (p1[2] + p2[2])};
+}
+
+inline Point3D paxpy(const Point3D& p1, const double alpha, const Point3D& p2) {
+    return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]};
+}
+
+/**
+ *
+ * \param e_pos electrode position
+ * \param seg_pos segment position
+ * \param radius segment radius
+ * \param double conductivity factor 1/([4 pi] * [conductivity])
+ * \return Resistance of the medium from the segment to the electrode.
+ */
+inline double point_source_lfp_factor(const Point3D& e_pos,
+                                      const Point3D& seg_pos,
+                                      const double radius,
+                                      const double f) {
+    nrn_assert(radius >= 0.0);
+    Point3D es = paxpy(e_pos, -1.0, seg_pos);
+    return f / std::max(norm(es), radius);
+}
+
+/**
+ *
+ * \param e_pos electrode position
+ * \param seg_pos segment position
+ * \param radius segment radius
+ * \param f conductivity factor 1/([4 pi] * [conductivity])
+ * \return Resistance of the medium from the segment to the electrode.
+ */
+double line_source_lfp_factor(const Point3D& e_pos,
+                              const Point3D& seg_0,
+                              const Point3D& seg_1,
+                              const double radius,
+                              const double f);
+}  // namespace lfputils
+
+enum LFPCalculatorType { LineSource, PointSource };
+
+/**
+ * \brief LFPCalculator allows calculation of LFP given membrane currents.
+ */
+template <LFPCalculatorType Ty, typename SegmentIdTy = int>
+struct LFPCalculator {
+    /**
+     * LFP Calculator constructor
+     * \param seg_start all segments start owned by the proc
+     * \param seg_end all segments end owned by the proc
+     * \param radius fence around the segment. Ensures electrode cannot be
+     * arbitrarily close to the segment
+     * \param electrodes positions of the electrodes
+     * \param extra_cellular_conductivity conductivity of the extra-cellular
+     * medium
+     */
+    LFPCalculator(const lfputils::Point3Ds& seg_start,
+                  const lfputils::Point3Ds& seg_end,
+                  const std::vector<double>& radius,
+                  const std::vector<SegmentIdTy>& segment_ids,
+                  const lfputils::Point3Ds& electrodes,
+                  double extra_cellular_conductivity);
+
+    template <typename Vector>
+    void lfp(const Vector& membrane_current);
+
+    const std::vector<double>& lfp_values() const noexcept {
+        return lfp_values_;
+    }
+
+  private:
+    inline double getFactor(const lfputils::Point3D& e_pos,
+                            const lfputils::Point3D& seg_0,
+                            const lfputils::Point3D& seg_1,
+                            const double radius,
+                            const double f) const;
+    std::vector<double> lfp_values_;
+    std::vector<std::vector<double>> m;
+    const std::vector<SegmentIdTy>& segment_ids_;
+};
+
+template <>
+double LFPCalculator<LineSource>::getFactor(const lfputils::Point3D& e_pos,
+                                            const lfputils::Point3D& seg_0,
+                                            const lfputils::Point3D& seg_1,
+                                            const double radius,
+                                            const double f) const {
+    return lfputils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f);
+}
+
+template <>
+double LFPCalculator<PointSource>::getFactor(const lfputils::Point3D& e_pos,
+                                             const lfputils::Point3D& seg_0,
+                                             const lfputils::Point3D& seg_1,
+                                             const double radius,
+                                             const double f) const {
+    return lfputils::point_source_lfp_factor(e_pos, lfputils::barycenter(seg_0, seg_1), radius, f);
+}
+
+extern template LFPCalculator<LineSource>::LFPCalculator(const lfputils::Point3Ds& seg_start,
+                                                         const lfputils::Point3Ds& seg_end,
+                                                         const std::vector<double>& radius,
+                                                         const std::vector<int>& segment_ids,
+                                                         const lfputils::Point3Ds& electrodes,
+                                                         double extra_cellular_conductivity);
+extern template LFPCalculator<PointSource>::LFPCalculator(const lfputils::Point3Ds& seg_start,
+                                                          const lfputils::Point3Ds& seg_end,
+                                                          const std::vector<double>& radius,
+                                                          const std::vector<int>& segment_ids,
+                                                          const lfputils::Point3Ds& electrodes,
+                                                          double extra_cellular_conductivity);
+extern template void LFPCalculator<LineSource>::lfp(const lfputils::DoublePtr& membrane_current);
+extern template void LFPCalculator<PointSource>::lfp(const lfputils::DoublePtr& membrane_current);
+extern template void LFPCalculator<LineSource>::lfp(const std::vector<double>& membrane_current);
+extern template void LFPCalculator<PointSource>::lfp(const std::vector<double>& membrane_current);
+}  // namespace coreneuron

coreneuron/mechanism/register_mech.cpp

@@ -19,10 +19,11 @@
 
 namespace coreneuron {
 int secondorder = 0;
-double t, dt, celsius;
+double t, dt, celsius, pi;
 // declare copyin required for correct initialization
 #pragma acc declare copyin(secondorder)
 #pragma acc declare copyin(celsius)
+#pragma acc declare copyin(pi)
 int rev_dt;
 
 using Pfrv = void (*)();

coreneuron/nrnconf.h

@@ -32,6 +32,7 @@ using Symbol = char;
 #define VECTORIZE   1
 
 extern double celsius;
+extern double pi;
 
 extern double t, dt;
 extern int rev_dt;

tests/CMakeLists.txt

@@ -29,6 +29,7 @@ if(Boost_FOUND)
     add_subdirectory(unit/interleave_info)
     add_subdirectory(unit/alignment)
     add_subdirectory(unit/queueing)
+    add_subdirectory(unit/lfp)
   endif()
   message(STATUS "Boost found, unit tests enabled")
 else()

tests/unit/alignment/alignment.cpp

@@ -13,9 +13,8 @@
 #include <cstring>
 #include <stdint.h>
 
-#include <boost/test/unit_test.hpp>
-#include <boost/test/test_case_template.hpp>
 #include <boost/mpl/list.hpp>
+#include <boost/test/unit_test.hpp>
 
 #include "coreneuron/utils/memory.h"
 

tests/unit/lfp/CMakeLists.txt

@@ -0,0 +1,26 @@
+# =============================================================================
+# Copyright (C) 2016-2020 Blue Brain Project
+#
+# See top-level LICENSE file for details.
+# =============================================================================
+
+include_directories(${CMAKE_SOURCE_DIR}/coreneuron ${Boost_INCLUDE_DIRS})
+file(GLOB lfp_test_src "*.cpp")
+
+add_executable(lfp_test_bin ${lfp_test_src})
+target_link_libraries(
+  lfp_test_bin
+  ${MPI_CXX_LIBRARIES}
+  ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
+  coreneuron
+  ${corenrn_mech_lib}
+  ${reportinglib_LIBRARY}
+  ${sonatareport_LIBRARY})
+
+if(CORENRN_ENABLE_GPU)
+  target_link_libraries(lfp_test_bin ${link_cudacoreneuron} ${CUDA_LIBRARIES})
+endif()
+
+add_dependencies(lfp_test_bin nrniv-core)
+
+add_test(NAME lfp_test COMMAND ${TEST_EXEC_PREFIX} ${CMAKE_CURRENT_BINARY_DIR}/lfp_test_bin)