diff --git a/coreneuron/apps/main1.cpp b/coreneuron/apps/main1.cpp index fa67d8fe7..262dcb7b3 100644 --- a/coreneuron/apps/main1.cpp +++ b/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; diff --git a/coreneuron/io/global_vars.cpp b/coreneuron/io/global_vars.cpp index f69f6256c..c9890cc4e 100644 --- a/coreneuron/io/global_vars.cpp +++ b/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 diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp new file mode 100644 index 000000000..64e7393c5 --- /dev/null +++ b/coreneuron/io/lfp.cpp @@ -0,0 +1,141 @@ +#include "coreneuron/io/lfp.hpp" + +#include +#include +#include + + +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::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::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 +LFPCalculator::LFPCalculator(const Point3Ds& seg_start, + const Point3Ds& seg_end, + const std::vector& radius, + const std::vector& 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 +template +inline void LFPCalculator::lfp(const Vector& membrane_current) { + std::vector 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::LFPCalculator(const lfputils::Point3Ds& seg_start, + const lfputils::Point3Ds& seg_end, + const std::vector& radius, + const std::vector& segment_ids, + const lfputils::Point3Ds& electrodes, + double extra_cellular_conductivity); +template LFPCalculator::LFPCalculator(const lfputils::Point3Ds& seg_start, + const lfputils::Point3Ds& seg_end, + const std::vector& radius, + const std::vector& segment_ids, + const lfputils::Point3Ds& electrodes, + double extra_cellular_conductivity); +template void LFPCalculator::lfp(const DoublePtr& membrane_current); +template void LFPCalculator::lfp(const DoublePtr& membrane_current); +template void LFPCalculator::lfp(const std::vector& membrane_current); +template void LFPCalculator::lfp(const std::vector& membrane_current); + +} // namespace coreneuron diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp new file mode 100644 index 000000000..e47da3c7d --- /dev/null +++ b/coreneuron/io/lfp.hpp @@ -0,0 +1,142 @@ +#pragma once + +#include +#include + +#include "coreneuron/mpi/nrnmpidec.h" +#include "coreneuron/nrnconf.h" +#include "coreneuron/utils/nrn_assert.h" + +namespace coreneuron { + +namespace lfputils { + +using Point3D = std::array; +using Point3Ds = std::vector; +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 +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& radius, + const std::vector& segment_ids, + const lfputils::Point3Ds& electrodes, + double extra_cellular_conductivity); + + template + void lfp(const Vector& membrane_current); + + const std::vector& 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 lfp_values_; + std::vector> m; + const std::vector& segment_ids_; +}; + +template <> +double LFPCalculator::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::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::LFPCalculator(const lfputils::Point3Ds& seg_start, + const lfputils::Point3Ds& seg_end, + const std::vector& radius, + const std::vector& segment_ids, + const lfputils::Point3Ds& electrodes, + double extra_cellular_conductivity); +extern template LFPCalculator::LFPCalculator(const lfputils::Point3Ds& seg_start, + const lfputils::Point3Ds& seg_end, + const std::vector& radius, + const std::vector& segment_ids, + const lfputils::Point3Ds& electrodes, + double extra_cellular_conductivity); +extern template void LFPCalculator::lfp(const lfputils::DoublePtr& membrane_current); +extern template void LFPCalculator::lfp(const lfputils::DoublePtr& membrane_current); +extern template void LFPCalculator::lfp(const std::vector& membrane_current); +extern template void LFPCalculator::lfp(const std::vector& membrane_current); +} // namespace coreneuron diff --git a/coreneuron/mechanism/register_mech.cpp b/coreneuron/mechanism/register_mech.cpp index 05f8dbe20..3e5b452cc 100644 --- a/coreneuron/mechanism/register_mech.cpp +++ b/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 (*)(); diff --git a/coreneuron/nrnconf.h b/coreneuron/nrnconf.h index b3fb173a3..064cc1547 100644 --- a/coreneuron/nrnconf.h +++ b/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; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4d2ec88d1..0e9d3b965 100644 --- a/tests/CMakeLists.txt +++ b/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() diff --git a/tests/unit/alignment/alignment.cpp b/tests/unit/alignment/alignment.cpp index 0bbf2ec09..f654d27a6 100644 --- a/tests/unit/alignment/alignment.cpp +++ b/tests/unit/alignment/alignment.cpp @@ -13,9 +13,8 @@ #include #include -#include -#include #include +#include #include "coreneuron/utils/memory.h" diff --git a/tests/unit/lfp/CMakeLists.txt b/tests/unit/lfp/CMakeLists.txt new file mode 100644 index 000000000..c20c61a3a --- /dev/null +++ b/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) diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp new file mode 100644 index 000000000..b65a378a0 --- /dev/null +++ b/tests/unit/lfp/lfp.cpp @@ -0,0 +1,103 @@ +#define BOOST_TEST_MODULE LFPTest +#define BOOST_TEST_MAIN + +#include + +#include + +#include "coreneuron/io/lfp.hpp" +#include "coreneuron/mpi/nrnmpi.h" + +using namespace coreneuron; +using namespace coreneuron::lfputils; + +template +double integral(F f, double a, double b, int n) { + double step = (b - a) / n; // width of each small rectangle + double area = 0.0; // signed area + for (int i = 0; i < n; i++) { + area += f(a + (i + 0.5) * step) * step; // sum up each small rectangle + } + return area; +} + + +BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { +#if NRNMPI + nrnmpi_init(nullptr, nullptr); +#endif + double segment_length{1.0e-6}; + double segment_start_val{1.0e-6}; + std::array segment_start = std::array{0.0, 0.0, segment_start_val}; + std::array segment_end = + paxpy(segment_start, 1.0, std::array{0.0, 0.0, segment_length}); + double floor{1.0e-6}; + pi = 3.141592653589; + + std::array vals; + double circling_radius{1.0e-6}; + std::array segment_middle{0.0, 0.0, 1.5e-6}; + double medium_resistivity_fac{1.0}; + for (auto k = 0; k < 10; k++) { + std::array approaching_elec = + paxpy(segment_middle, 1.0, std::array{0.0, 1.0e-5 - k * 1.0e-6, 0.0}); + std::array circling_elec = + paxpy(segment_middle, + 1.0, + std::array{0.0, + circling_radius * std::cos(2.0 * pi * k / 10), + circling_radius * std::sin(2.0 * pi * k / 10)}); + + double analytic_approaching_lfp = line_source_lfp_factor( + approaching_elec, segment_start, segment_end, floor, medium_resistivity_fac); + double analytic_circling_lfp = line_source_lfp_factor( + circling_elec, segment_start, segment_end, floor, medium_resistivity_fac); + double numeric_circling_lfp = integral( + [&](double x) { + return 1.0 / std::max(floor, + norm(paxpy(circling_elec, + -1.0, + paxpy(segment_end, + x, + paxpy(segment_start, -1.0, segment_end))))); + }, + 0.0, + 1.0, + 10000); + // TEST of analytic vs numerical integration + std::clog << "ANALYTIC line source " << analytic_circling_lfp + << " vs NUMERIC line source LFP " << numeric_circling_lfp << "\n"; + BOOST_REQUIRE_CLOSE(analytic_circling_lfp, numeric_circling_lfp, 1.0e-6); + // TEST of LFP Flooring + BOOST_REQUIRE((approaching_elec[1] < 0.866e-6) ? analytic_approaching_lfp == 1.0e6 : true); + vals[k] = analytic_circling_lfp; + } + // TEST of SYMMETRY of LFP FORMULA + for (size_t k = 0; k < 5; k++) { + BOOST_REQUIRE(std::abs((vals[k] - vals[k + 5]) / + std::max(std::abs(vals[k]), std::abs(vals[k + 5]))) < 1.0e-12); + } + std::vector> segments_starts = {{0., 0., 1.}, + {0., 0., 0.5}, + {0.0, 0.0, 0.0}, + {0.0, 0.0, -0.5}}; + std::vector> segments_ends = {{0., 0., 0.}, + {0., 0., 1.}, + {0., 0., 0.5}, + {0.0, 0.0, 0.0}}; + std::vector radii{0.1, 0.1, 0.1, 0.1}; + std::vector> electrodes = {{0.0, 0.3, 0.0}, {0.0, 0.7, 0.8}}; + std::vector indices = {0, 1, 2, 3}; + LFPCalculator lfp(segments_starts, segments_ends, radii, indices, electrodes, 1.0); + lfp.template lfp>({0.0, 1.0, 2.0, 3.0}); + std::vector res_line_source = lfp.lfp_values(); + LFPCalculator lfpp( + segments_starts, segments_ends, radii, indices, electrodes, 1.0); + lfpp.template lfp>({0.0, 1.0, 2.0, 3.0}); + std::vector res_point_source = lfpp.lfp_values(); + BOOST_REQUIRE_CLOSE(res_line_source[0], res_point_source[0], 1.0); + BOOST_REQUIRE_CLOSE(res_line_source[1], res_point_source[1], 1.0); +#if NRNMPI + nrnmpi_finalize(); +#endif +} diff --git a/tests/unit/queueing/test_queueing.cpp b/tests/unit/queueing/test_queueing.cpp index 21ae3297f..072c8163c 100644 --- a/tests/unit/queueing/test_queueing.cpp +++ b/tests/unit/queueing/test_queueing.cpp @@ -9,9 +9,8 @@ #define BOOST_TEST_MODULE QueueingTest #define TYPE T::cont -#include -#include #include +#include #include #include