From 6b8e1ea724fcd47c72028f5115bb5982a3d02205 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Fri, 5 Feb 2021 12:45:43 +0100 Subject: [PATCH 01/36] LFP classes + unit test. --- coreneuron/io/lfp.hpp | 272 +++++++++++++++++++++++++++++++++++++++++ tests/unit/lfp/lfp.cpp | 102 ++++++++++++++++ 2 files changed, 374 insertions(+) create mode 100644 coreneuron/io/lfp.hpp create mode 100644 tests/unit/lfp/lfp.cpp diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp new file mode 100644 index 000000000..6f4a21d64 --- /dev/null +++ b/coreneuron/io/lfp.hpp @@ -0,0 +1,272 @@ +#ifndef AREA_LFP_H +#define AREA_LFP_H + +#include +#include +#include +#include +#include + +#include +#include + +namespace coreneuron { + + template + using Array3 = std::array; + + template + F dot(const Point3D& p1, const Point3D& p2) { + return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; + } + + template + F norm(const Point3D& p1) { + return std::sqrt(dot(p1, p1)); + } + + template + Array3 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])}; + } + + template + Array3 axp(const Point3D& p1, const F alpha, const Point3D& p2) { + return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]}; + } + + /** + * + * \tparam Point3D Type of 3D point + * \tparam F Floating point type + * \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. + */ + template + F point_source_lfp_factor( + const Point3D& e_pos, + const Point3D& seg_pos, + const F radius, + const F f + ) + { + nrn_assert(radius >= 0.0); + Array3 es = axp(e_pos, -1.0, seg_pos); + return f / std::max(norm(es),radius); + } + + /** + * + * \tparam Point3D Type of 3D point + * \tparam F Floating point type + * \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. + */ + template + F line_source_lfp_factor( + const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, + const F radius, + const F f + ) + { + nrn_assert(radius >= 0.0); + Array3 dx = axp(seg_1, -1.0, seg_0); + Array3 de = axp(e_pos, -1.0, seg_0); + F dx2(dot,F>(dx, dx)); + F dxn(std::sqrt(dx2)); + if (dxn < 1.0e-20) { + return point_source_lfp_factor(e_pos, seg_0, radius, f); + } + F de2(dot,F>(de, de)); + F mu(dot,F>(dx, de) / dx2); + Array3 de_star(axp(de, -mu, dx)); + F de_star2(dot,F>(de_star, de_star)); + F q2(de_star2 / dx2); + + F delta(mu * mu - (de2 - radius * radius) / dx2); + F one_m_mu(1.0 - mu); + auto log_integral = [&q2, &dxn](F a, F b) { + if (std::abs(q2) < std::numeric_limits::epsilon()) { + if (a * b <= 0) { + throw std::invalid_argument( + "Log integral: invalid arguments " + std::to_string(b) + " " + std::to_string(a) + + ". Likely electrode exactly on the segment and no flooring is present."); + } + 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 { + F sqr_delta(std::sqrt(delta)); + F d1(mu - sqr_delta); + F d2(mu + sqr_delta); + F parts = 0.0; + if (d1 > 0.0) { + F b(std::min(d1, 1.0) - mu); + parts += log_integral(-mu, b); + } + if (d2 < 1.0) { + F b(std::max(d2, 0.0) - mu); + parts += log_integral(b, one_m_mu); + }; + // complement + F 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; + }; + + } + + enum LFPCalculatorType { + LineSource, + PointSource + }; + + /** + * \brief LFPCalculator allows calculation of LFP given membrane currents. + */ + template + struct LFPCalculator { + + /** + * LFP Calculator constructor + * \tparam Point3Ds A vector of 3D points type + * \tparam Vector A vector of floats type + * \param comm MPI communicator + * \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 + */ + template + LFPCalculator(MPI_Comm comm, + const Point3Ds& seg_start, + const Point3Ds& seg_end, + const Vector& radius, + const std::vector& segment_ids, + const Point3Ds& electrodes, + double extra_cellular_conductivity) + : comm_(comm) + , segment_ids_(segment_ids) + { + if (seg_start.size() != seg_end.size()) { + throw std::logic_error("Wrong number of segment starts or ends."); + } + if (seg_start.size() != radius.size()) { + throw std::logic_error("Wrong number of radius size."); + } + double f(1.0 / (extra_cellular_conductivity * 4.0 * boost::math::constants::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 + void 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]]; + } + } + std::vector res_reduced(res.size()); + int err = MPI_Allreduce(res.data(), + res_reduced.data(), + res.size(), + MPI_DOUBLE, + MPI_SUM, + comm_ + ); + lfp_values = res_reduced; + } + + std::vector lfp_values; + + private: + + template + inline F getFactor(const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, + const F radius, + const F f) const; + + std::vector > m; + MPI_Comm comm_; + const std::vector& segment_ids_; + }; + + template <> + template + F LFPCalculator::getFactor( + const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, + const F radius, + const F f) const { + return line_source_lfp_factor( + e_pos, + seg_0, + seg_1, + radius, + f + ); + } + + template <> + template + F LFPCalculator::getFactor( + const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, + const F radius, + const F f) const { + return point_source_lfp_factor( + e_pos, + barycenter(seg_0, + seg_1), + radius, + f + ); + } + + +}; + +#endif //AREA_LFP_H + diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp new file mode 100644 index 000000000..6f35120eb --- /dev/null +++ b/tests/unit/lfp/lfp.cpp @@ -0,0 +1,102 @@ +#define BOOST_TEST_MODULE LFPTest +#define BOOST_TEST_MAIN + +#include +#include + +#include +#include +#include "coreneuron/io/lfp.hpp" + +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; +} + +std::array operator * (const double x, const std::array& y) { + return {x * y[0], x * y[1], x * y[2]}; +} + +std::array operator + (const std::array& x, const std::array& y) { + return {x[0] + y[0], x[1] + y[1], x[2] + y[2]}; +} + +BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { + MPI_Init(nullptr, nullptr); + using namespace coreneuron; + double segment_length{1.0e-6}; + double segment_start{1.0e-6}; + std::array segment_start = {0.0, 0.0, segment_start}; + std::array segment_end = {0.0, 0.0, segment_start + segment_length}; + double floor{1.0e-6}; + double 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 = segment_middle + + {0.0, 1.0e-5 - k * 1.0e-6, 0.0}; + std::array circling_elec = segment_middle + + {0.0, circling_radius * std::cos(2.0 * pi * k / 10), + circling_radius * std::sin(2.0 * pi * k / 10)}; + + double l = line_source_lfp_factor( + approaching_elec, + segment_start, + segment_end, + floor, + medium_resistivity_fac); + double l2 = line_source_lfp_factor( + circling_elec, + segment_start, + segment_end, + floor, + medium_resistivity_fac); + double l3 = integral ([&](double x) { + return 1.0 / std::max(floor, + norm,double>( + circling_elec + + (-x) * segment_start + + (x - 1) * segment_end)); + }, 0.0, 1.0, 10000); + // TEST of analytic vs numerical integration + BOOST_REQUIRE_CLOSE(l2, l3, 1.0e-6); + // TEST of LFP Flooring + BOOST_REQUIRE(e[1] < 0.866e-6 || l != 1.0e6); + vals[k] = l2; + } + // TEST of SYMMETRY of LFP FORMULA + for (size_t k = 0; k < 5; k++) { + BOOST_REQUIRE(std::abs((vals[k] - vals[k + 5]) * 0.5 / (vals[k] + vals[k + 5])) + < 1.0e-12); + } + std::vector > start = {{0., 0., 1.}, + {0., 0., 0.5}, + {0.0, 0.0, 0.0}, + {0.0, 0.0, -0.5}}; + std::vector > end = {{0., 0., 0.}, + {0., 0., 1.}, + {0., 0., 0.5}, + {0.0, 0.0, 0.0}}; + std::vector radius{0.1, 0.1, 0.1, 0.1}; + std::vector > elect = {{0.0, 0.3, 0.0}, + {0.0, 0.7, 0.8}}; + + std::vector indices = {0, 1, 2, 3}; + LFPCalculator lfp(MPI_COMM_WORLD, start, end, radius, indices, elect, 1.0); + std::vector res_line_source = lfp.template lfp>({0.0, 1.0, 2.0, 3.0}); + LFPCalculator lfpp(MPI_COMM_WORLD, start, end, radius, indices, elect, 1.0); + std::vector res_point_source = lfpp.template lfp>({0.0, 1.0, 2.0, 3.0}); + BOOST_REQUIRE_CLOSE(res_line_source[0], res_point_source[0], 1.0e-3); + BOOST_REQUIRE_CLOSE(res_line_source[1], res_point_source[1], 1.0e-3); + MPI_Finalize(); + +} + From d7f8b722ff008eb6bac725d78df11de19e9b64a5 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Fri, 5 Feb 2021 13:59:30 +0100 Subject: [PATCH 02/36] Cosmetic. --- coreneuron/io/lfp.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 6f4a21d64..5175e6df0 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -1,5 +1,5 @@ -#ifndef AREA_LFP_H -#define AREA_LFP_H +#ifndef LFP_H +#define LFP_H #include #include @@ -212,6 +212,9 @@ namespace coreneuron { MPI_SUM, comm_ ); + if (err != MPI_SUCCESS) { + MPI_Abort(comm_, err); + } lfp_values = res_reduced; } From c62ec29824cdd13180976f73dec94a9ea501b2b9 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Fri, 5 Feb 2021 14:10:51 +0100 Subject: [PATCH 03/36] Cosmetic. --- coreneuron/io/lfp.hpp | 4 ++-- tests/unit/lfp/lfp.cpp | 48 +++++++++++++++++++++++------------------- 2 files changed, 28 insertions(+), 24 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 5175e6df0..19e8afa17 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -172,10 +172,10 @@ namespace coreneuron { , segment_ids_(segment_ids) { if (seg_start.size() != seg_end.size()) { - throw std::logic_error("Wrong number of segment starts or ends."); + throw std::invalid_argument("Different number of segment starts and ends."); } if (seg_start.size() != radius.size()) { - throw std::logic_error("Wrong number of radius size."); + throw std::invalid_argument("Different number of segments and radii."); } double f(1.0 / (extra_cellular_conductivity * 4.0 * boost::math::constants::pi())); diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 6f35120eb..0b3d83e22 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -32,7 +32,8 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { double segment_length{1.0e-6}; double segment_start{1.0e-6}; std::array segment_start = {0.0, 0.0, segment_start}; - std::array segment_end = {0.0, 0.0, segment_start + segment_length}; + std::array segment_end = + segment_start + {0.0, 0.0, segment_length}; double floor{1.0e-6}; double pi{3.141592653589}; @@ -47,19 +48,19 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { {0.0, circling_radius * std::cos(2.0 * pi * k / 10), circling_radius * std::sin(2.0 * pi * k / 10)}; - double l = line_source_lfp_factor( + double analytic_approaching_lfp = line_source_lfp_factor( approaching_elec, segment_start, segment_end, floor, medium_resistivity_fac); - double l2 = line_source_lfp_factor( + double analytic_circling_lfp = line_source_lfp_factor( circling_elec, segment_start, segment_end, floor, medium_resistivity_fac); - double l3 = integral ([&](double x) { + double numeric_circling_lfp = integral ([&](double x) { return 1.0 / std::max(floor, norm,double>( circling_elec @@ -67,36 +68,39 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { + (x - 1) * segment_end)); }, 0.0, 1.0, 10000); // TEST of analytic vs numerical integration - BOOST_REQUIRE_CLOSE(l2, l3, 1.0e-6); + BOOST_REQUIRE_CLOSE(analytic_circling_lfp, + numeric_circling_lfp, 1.0e-6); // TEST of LFP Flooring - BOOST_REQUIRE(e[1] < 0.866e-6 || l != 1.0e6); - vals[k] = l2; + 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]) * 0.5 / (vals[k] + vals[k + 5])) < 1.0e-12); } - std::vector > start = {{0., 0., 1.}, - {0., 0., 0.5}, - {0.0, 0.0, 0.0}, - {0.0, 0.0, -0.5}}; - std::vector > end = {{0., 0., 0.}, - {0., 0., 1.}, - {0., 0., 0.5}, - {0.0, 0.0, 0.0}}; - std::vector radius{0.1, 0.1, 0.1, 0.1}; - std::vector > elect = {{0.0, 0.3, 0.0}, - {0.0, 0.7, 0.8}}; - + 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(MPI_COMM_WORLD, start, end, radius, indices, elect, 1.0); + LFPCalculator lfp(MPI_COMM_WORLD, segments_starts, segments_ends, radii, indices, electrodes, 1.0); std::vector res_line_source = lfp.template lfp>({0.0, 1.0, 2.0, 3.0}); - LFPCalculator lfpp(MPI_COMM_WORLD, start, end, radius, indices, elect, 1.0); + LFPCalculator lfpp(MPI_COMM_WORLD, segments_starts, segments_ends, radii, indices, electrodes, 1.0); std::vector res_point_source = lfpp.template lfp>({0.0, 1.0, 2.0, 3.0}); BOOST_REQUIRE_CLOSE(res_line_source[0], res_point_source[0], 1.0e-3); BOOST_REQUIRE_CLOSE(res_line_source[1], res_point_source[1], 1.0e-3); MPI_Finalize(); - } From de64826bfe31e0a5bf18a1b1cb73708b092cb883 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Fri, 5 Feb 2021 14:56:59 +0100 Subject: [PATCH 04/36] Polishing test. Make it compile. --- coreneuron/io/lfp.hpp | 3 ++- tests/CMakeLists.txt | 1 + tests/unit/lfp/lfp.cpp | 30 +++++++++++++++++------------- 3 files changed, 20 insertions(+), 14 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 19e8afa17..3c623aa72 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -1,7 +1,6 @@ #ifndef LFP_H #define LFP_H -#include #include #include #include @@ -10,6 +9,8 @@ #include #include +#include "coreneuron/utils/nrn_assert.h" + namespace coreneuron { template 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/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 0b3d83e22..2c816a5a7 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -30,10 +30,10 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { MPI_Init(nullptr, nullptr); using namespace coreneuron; double segment_length{1.0e-6}; - double segment_start{1.0e-6}; - std::array segment_start = {0.0, 0.0, segment_start}; + double segment_start_val{1.0e-6}; + std::array segment_start = std::array{0.0, 0.0, segment_start_val}; std::array segment_end = - segment_start + {0.0, 0.0, segment_length}; + segment_start + std::array{0.0, 0.0, segment_length}; double floor{1.0e-6}; double pi{3.141592653589}; @@ -43,9 +43,9 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { double medium_resistivity_fac{1.0}; for (auto k = 0; k < 10; k++) { std::array approaching_elec = segment_middle + - {0.0, 1.0e-5 - k * 1.0e-6, 0.0}; + std::array{0.0, 1.0e-5 - k * 1.0e-6, 0.0}; std::array circling_elec = segment_middle + - {0.0, circling_radius * std::cos(2.0 * pi * k / 10), + 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( @@ -61,14 +61,16 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { floor, medium_resistivity_fac); double numeric_circling_lfp = integral ([&](double x) { - return 1.0 / std::max(floor, + return 1.0 / std::max(floor, norm,double>( - circling_elec - + (-x) * segment_start + circling_elec + + (-x) * segment_start + (x - 1) * segment_end)); }, 0.0, 1.0, 10000); // TEST of analytic vs numerical integration - BOOST_REQUIRE_CLOSE(analytic_circling_lfp, + 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 @@ -96,11 +98,13 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { {0.0, 0.7, 0.8}}; std::vector indices = {0, 1, 2, 3}; LFPCalculator lfp(MPI_COMM_WORLD, segments_starts, segments_ends, radii, indices, electrodes, 1.0); - std::vector res_line_source = lfp.template lfp>({0.0, 1.0, 2.0, 3.0}); + lfp.template lfp>({0.0, 1.0, 2.0, 3.0}); + std::vector res_line_source = lfp.lfp_values; LFPCalculator lfpp(MPI_COMM_WORLD, segments_starts, segments_ends, radii, indices, electrodes, 1.0); - std::vector res_point_source = lfpp.template lfp>({0.0, 1.0, 2.0, 3.0}); - BOOST_REQUIRE_CLOSE(res_line_source[0], res_point_source[0], 1.0e-3); - BOOST_REQUIRE_CLOSE(res_line_source[1], res_point_source[1], 1.0e-3); + 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); MPI_Finalize(); } From 8f3bd482e7fb3eebb5063796e09a7ad367980782 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Fri, 5 Feb 2021 15:46:00 +0100 Subject: [PATCH 05/36] Unit test CMakeLists.txt. --- tests/unit/lfp/CMakeLists.txt | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 tests/unit/lfp/CMakeLists.txt diff --git a/tests/unit/lfp/CMakeLists.txt b/tests/unit/lfp/CMakeLists.txt new file mode 100644 index 000000000..905a84c5f --- /dev/null +++ b/tests/unit/lfp/CMakeLists.txt @@ -0,0 +1,19 @@ + +# ============================================================================= +# 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 ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${MPI_CXX_LIBRARIES}) + +if(CORENRN_ENABLE_GPU) + target_link_libraries(lfp_test_bin ${link_cudacoreneuron} ${CUDA_LIBRARIES}) +endif() + +add_test(NAME lfp_test COMMAND ${TEST_EXEC_PREFIX} + ${CMAKE_CURRENT_BINARY_DIR}/lfp_test_bin) From 64ae68666776091ba86e291e7d5575afa03969b1 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 8 Feb 2021 14:42:17 +0100 Subject: [PATCH 06/36] Rework LFP using MPI wrappers. --- coreneuron/io/lfp.hpp | 22 ++++++---------------- tests/unit/lfp/CMakeLists.txt | 7 ++++++- tests/unit/lfp/lfp.cpp | 11 ++++++----- 3 files changed, 18 insertions(+), 22 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 3c623aa72..146eb83f9 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -6,10 +6,11 @@ #include #include -#include + #include #include "coreneuron/utils/nrn_assert.h" +#include "coreneuron/mpi/nrnmpidec.h" namespace coreneuron { @@ -162,15 +163,13 @@ namespace coreneuron { * \param extra_cellular_conductivity conductivity of the extra-cellular medium */ template - LFPCalculator(MPI_Comm comm, - const Point3Ds& seg_start, + LFPCalculator(const Point3Ds& seg_start, const Point3Ds& seg_end, const Vector& radius, const std::vector& segment_ids, const Point3Ds& electrodes, double extra_cellular_conductivity) - : comm_(comm) - , segment_ids_(segment_ids) + :segment_ids_(segment_ids) { if (seg_start.size() != seg_end.size()) { throw std::invalid_argument("Different number of segment starts and ends."); @@ -206,16 +205,8 @@ namespace coreneuron { } } std::vector res_reduced(res.size()); - int err = MPI_Allreduce(res.data(), - res_reduced.data(), - res.size(), - MPI_DOUBLE, - MPI_SUM, - comm_ - ); - if (err != MPI_SUCCESS) { - MPI_Abort(comm_, err); - } + int mpi_sum{1}; + nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); lfp_values = res_reduced; } @@ -231,7 +222,6 @@ namespace coreneuron { const F f) const; std::vector > m; - MPI_Comm comm_; const std::vector& segment_ids_; }; diff --git a/tests/unit/lfp/CMakeLists.txt b/tests/unit/lfp/CMakeLists.txt index 905a84c5f..7db0c825e 100644 --- a/tests/unit/lfp/CMakeLists.txt +++ b/tests/unit/lfp/CMakeLists.txt @@ -9,7 +9,12 @@ 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 ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${MPI_CXX_LIBRARIES}) +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}) diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 2c816a5a7..91ca150af 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -2,11 +2,12 @@ #define BOOST_TEST_MAIN #include -#include #include #include + #include "coreneuron/io/lfp.hpp" +#include "coreneuron/mpi/nrnmpi.h" template double integral(F f, double a, double b, int n) { @@ -27,8 +28,8 @@ std::array operator + (const std::array& x, const std::array } BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { - MPI_Init(nullptr, nullptr); using namespace coreneuron; + nrnmpi_init(nullptr, nullptr); 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}; @@ -97,14 +98,14 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { {{0.0, 0.3, 0.0}, {0.0, 0.7, 0.8}}; std::vector indices = {0, 1, 2, 3}; - LFPCalculator lfp(MPI_COMM_WORLD, segments_starts, segments_ends, radii, indices, electrodes, 1.0); + 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(MPI_COMM_WORLD, segments_starts, segments_ends, radii, indices, electrodes, 1.0); + 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); - MPI_Finalize(); + nrnmpi_finalize(); } From 83ce66063007a73f630e4c1e59e9fb802a23a1ec Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 8 Feb 2021 14:57:21 +0100 Subject: [PATCH 07/36] Clang format. --- coreneuron/io/lfp.hpp | 435 +++++++++++++++++++---------------------- tests/unit/lfp/lfp.cpp | 176 ++++++++--------- 2 files changed, 282 insertions(+), 329 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 146eb83f9..2dd66fcfd 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -6,7 +6,6 @@ #include #include - #include #include "coreneuron/utils/nrn_assert.h" @@ -14,253 +13,211 @@ namespace coreneuron { - template - using Array3 = std::array; - - template - F dot(const Point3D& p1, const Point3D& p2) { - return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; +template using Array3 = std::array; + +template +F dot(const Point3D &p1, const Point3D &p2) { + return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; +} + +template F norm(const Point3D &p1) { + return std::sqrt(dot(p1, p1)); +} + +template +Array3 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]) }; +} + +template +Array3 axp(const Point3D &p1, const F alpha, const Point3D &p2) { + return { p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], + p1[2] + alpha * p2[2] }; +} + +/** + * + * \tparam Point3D Type of 3D point + * \tparam F Floating point type + * \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. + */ +template +F point_source_lfp_factor(const Point3D &e_pos, const Point3D &seg_pos, + const F radius, const F f) { + nrn_assert(radius >= 0.0); + Array3 es = axp(e_pos, -1.0, seg_pos); + return f / std::max(norm(es), radius); +} + +/** + * + * \tparam Point3D Type of 3D point + * \tparam F Floating point type + * \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. + */ +template +F line_source_lfp_factor(const Point3D &e_pos, const Point3D &seg_0, + const Point3D &seg_1, const F radius, const F f) { + nrn_assert(radius >= 0.0); + Array3 dx = axp(seg_1, -1.0, seg_0); + Array3 de = axp(e_pos, -1.0, seg_0); + F dx2(dot, F>(dx, dx)); + F dxn(std::sqrt(dx2)); + if (dxn < 1.0e-20) { + return point_source_lfp_factor(e_pos, seg_0, radius, f); + } + F de2(dot, F>(de, de)); + F mu(dot, F>(dx, de) / dx2); + Array3 de_star(axp(de, -mu, dx)); + F de_star2(dot, F>(de_star, de_star)); + F q2(de_star2 / dx2); + + F delta(mu * mu - (de2 - radius * radius) / dx2); + F one_m_mu(1.0 - mu); + auto log_integral = [&q2, &dxn](F a, F b) { + if (std::abs(q2) < std::numeric_limits::epsilon()) { + if (a * b <= 0) { + throw std::invalid_argument( + "Log integral: invalid arguments " + std::to_string(b) + " " + + std::to_string(a) + ". Likely electrode exactly on the segment and " + "no flooring is present."); + } + 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; } - - template - F norm(const Point3D& p1) { - return std::sqrt(dot(p1, p1)); + }; + if (delta <= 0.0) { + return f * log_integral(-mu, one_m_mu); + } else { + F sqr_delta(std::sqrt(delta)); + F d1(mu - sqr_delta); + F d2(mu + sqr_delta); + F parts = 0.0; + if (d1 > 0.0) { + F b(std::min(d1, 1.0) - mu); + parts += log_integral(-mu, b); } - - template - Array3 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])}; + if (d2 < 1.0) { + F b(std::max(d2, 0.0) - mu); + parts += log_integral(b, one_m_mu); + }; + // complement + F 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; + }; +} - template - Array3 axp(const Point3D& p1, const F alpha, const Point3D& p2) { - return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]}; - } +enum LFPCalculatorType { + LineSource, + PointSource +}; - /** - * - * \tparam Point3D Type of 3D point - * \tparam F Floating point type - * \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. - */ - template - F point_source_lfp_factor( - const Point3D& e_pos, - const Point3D& seg_pos, - const F radius, - const F f - ) - { - nrn_assert(radius >= 0.0); - Array3 es = axp(e_pos, -1.0, seg_pos); - return f / std::max(norm(es),radius); +/** + * \brief LFPCalculator allows calculation of LFP given membrane currents. + */ +template +struct LFPCalculator { + + /** + * LFP Calculator constructor + * \tparam Point3Ds A vector of 3D points type + * \tparam Vector A vector of floats type + * \param comm MPI communicator + * \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 + */ + template + LFPCalculator(const Point3Ds &seg_start, const Point3Ds &seg_end, + const 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."); } - - /** - * - * \tparam Point3D Type of 3D point - * \tparam F Floating point type - * \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. - */ - template - F line_source_lfp_factor( - const Point3D& e_pos, - const Point3D& seg_0, - const Point3D& seg_1, - const F radius, - const F f - ) - { - nrn_assert(radius >= 0.0); - Array3 dx = axp(seg_1, -1.0, seg_0); - Array3 de = axp(e_pos, -1.0, seg_0); - F dx2(dot,F>(dx, dx)); - F dxn(std::sqrt(dx2)); - if (dxn < 1.0e-20) { - return point_source_lfp_factor(e_pos, seg_0, radius, f); - } - F de2(dot,F>(de, de)); - F mu(dot,F>(dx, de) / dx2); - Array3 de_star(axp(de, -mu, dx)); - F de_star2(dot,F>(de_star, de_star)); - F q2(de_star2 / dx2); - - F delta(mu * mu - (de2 - radius * radius) / dx2); - F one_m_mu(1.0 - mu); - auto log_integral = [&q2, &dxn](F a, F b) { - if (std::abs(q2) < std::numeric_limits::epsilon()) { - if (a * b <= 0) { - throw std::invalid_argument( - "Log integral: invalid arguments " + std::to_string(b) + " " + std::to_string(a) + - ". Likely electrode exactly on the segment and no flooring is present."); - } - 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 { - F sqr_delta(std::sqrt(delta)); - F d1(mu - sqr_delta); - F d2(mu + sqr_delta); - F parts = 0.0; - if (d1 > 0.0) { - F b(std::min(d1, 1.0) - mu); - parts += log_integral(-mu, b); - } - if (d2 < 1.0) { - F b(std::max(d2, 0.0) - mu); - parts += log_integral(b, one_m_mu); - }; - // complement - F 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; - }; - + if (seg_start.size() != radius.size()) { + throw std::invalid_argument("Different number of segments and radii."); } - - enum LFPCalculatorType { - LineSource, - PointSource - }; - - /** - * \brief LFPCalculator allows calculation of LFP given membrane currents. - */ - template - struct LFPCalculator { - - /** - * LFP Calculator constructor - * \tparam Point3Ds A vector of 3D points type - * \tparam Vector A vector of floats type - * \param comm MPI communicator - * \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 - */ - template - LFPCalculator(const Point3Ds& seg_start, - const Point3Ds& seg_end, - const 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 * boost::math::constants::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 - void 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]]; - } - } - std::vector res_reduced(res.size()); - int mpi_sum{1}; - nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); - lfp_values = res_reduced; - } - - std::vector lfp_values; - - private: - - template - inline F getFactor(const Point3D& e_pos, - const Point3D& seg_0, - const Point3D& seg_1, - const F radius, - const F f) const; - - std::vector > m; - const std::vector& segment_ids_; - }; - - template <> - template - F LFPCalculator::getFactor( - const Point3D& e_pos, - const Point3D& seg_0, - const Point3D& seg_1, - const F radius, - const F f) const { - return line_source_lfp_factor( - e_pos, - seg_0, - seg_1, - radius, - f - ); + double f(1.0 / (extra_cellular_conductivity * 4.0 * + boost::math::constants::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 - F LFPCalculator::getFactor( - const Point3D& e_pos, - const Point3D& seg_0, - const Point3D& seg_1, - const F radius, - const F f) const { - return point_source_lfp_factor( - e_pos, - barycenter(seg_0, - seg_1), - radius, - f - ); + } + + template void 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]]; + } } - - + std::vector res_reduced(res.size()); + int mpi_sum{ 1 }; + nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), + mpi_sum); + lfp_values = res_reduced; + } + + std::vector lfp_values; + +private: + template + inline F getFactor(const Point3D &e_pos, const Point3D &seg_0, + const Point3D &seg_1, const F radius, const F f) const; + + std::vector > m; + const std::vector &segment_ids_; }; -#endif //AREA_LFP_H +template <> +template +F LFPCalculator::getFactor(const Point3D &e_pos, + const Point3D &seg_0, + const Point3D &seg_1, const F radius, + const F f) const { + return line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); +} + +template <> +template +F LFPCalculator::getFactor(const Point3D &e_pos, + const Point3D &seg_0, + const Point3D &seg_1, const F radius, + const F f) const { + return point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), + radius, f); +} +}; +#endif // AREA_LFP_H diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 91ca150af..93629fe5f 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -9,103 +9,99 @@ #include "coreneuron/io/lfp.hpp" #include "coreneuron/mpi/nrnmpi.h" -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; +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; } -std::array operator * (const double x, const std::array& y) { - return {x * y[0], x * y[1], x * y[2]}; +std::array operator*(const double x, + const std::array &y) { + return { x * y[0], x * y[1], x * y[2] }; } -std::array operator + (const std::array& x, const std::array& y) { - return {x[0] + y[0], x[1] + y[1], x[2] + y[2]}; +std::array operator+(const std::array &x, + const std::array &y) { + return { x[0] + y[0], x[1] + y[1], x[2] + y[2] }; } BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { - using namespace coreneuron; - nrnmpi_init(nullptr, nullptr); - 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 = - segment_start + std::array{0.0, 0.0, segment_length}; - double floor{1.0e-6}; - double pi{3.141592653589}; + using namespace coreneuron; + nrnmpi_init(nullptr, nullptr); + 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 = + segment_start + std::array{ 0.0, 0.0, segment_length }; + double floor{ 1.0e-6 }; + double 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 = segment_middle + - std::array{0.0, 1.0e-5 - k * 1.0e-6, 0.0}; - std::array circling_elec = segment_middle + - std::array{0.0, circling_radius * std::cos(2.0 * pi * k / 10), - circling_radius * std::sin(2.0 * pi * k / 10)}; + 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 = + segment_middle + std::array{ 0.0, 1.0e-5 - k * 1.0e-6, 0.0 }; + std::array circling_elec = + segment_middle + + 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,double>( - circling_elec - + (-x) * segment_start - + (x - 1) * 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]) * 0.5 / (vals[k] + 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); - nrnmpi_finalize(); + 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, double>( + circling_elec + (-x) * segment_start + + (x - 1) * 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]) * 0.5 / + (vals[k] + 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); + nrnmpi_finalize(); } - From 9b2204ee4a2057e199cc55dc52c84084bdf7c41a Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 8 Feb 2021 15:21:33 +0100 Subject: [PATCH 08/36] Only use MPI reduction if MPI is available.. --- coreneuron/io/lfp.hpp | 4 ++++ tests/unit/lfp/lfp.cpp | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 2dd66fcfd..a4ab64b92 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -183,9 +183,13 @@ struct LFPCalculator { } } std::vector res_reduced(res.size()); +#if NRNMPI int mpi_sum{ 1 }; nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); +#else + res_reduced = res; +#endif lfp_values = res_reduced; } diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 93629fe5f..c6a9fcb48 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -30,7 +30,9 @@ std::array operator+(const std::array &x, BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { using namespace coreneuron; +#ifdef NRNMPI nrnmpi_init(nullptr, nullptr); +#endif double segment_length{ 1.0e-6 }; double segment_start_val{ 1.0e-6 }; std::array segment_start = @@ -103,5 +105,7 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { 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 } From 5056df6136c930345b1224a6815b4344e562e494 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 8 Feb 2021 15:30:41 +0100 Subject: [PATCH 09/36] Fix previous commit. --- tests/unit/lfp/lfp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index c6a9fcb48..d17a5da64 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -30,7 +30,7 @@ std::array operator+(const std::array &x, BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { using namespace coreneuron; -#ifdef NRNMPI +#if NRNMPI nrnmpi_init(nullptr, nullptr); #endif double segment_length{ 1.0e-6 }; From 59a98b8d7f5f0f5f42201ee00605f5675a32a402 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 8 Feb 2021 16:03:10 +0100 Subject: [PATCH 10/36] Clang format. --- coreneuron/io/lfp.hpp | 297 ++++++++++++++++++++++-------------------- 1 file changed, 155 insertions(+), 142 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index a4ab64b92..fbe29fff8 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -13,27 +13,31 @@ namespace coreneuron { -template using Array3 = std::array; +template +using Array3 = std::array; template -F dot(const Point3D &p1, const Point3D &p2) { - return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; +F dot(const Point3D &p1, const Point3D &p2) +{ + return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; } -template F norm(const Point3D &p1) { - return std::sqrt(dot(p1, p1)); +template +F norm(const Point3D &p1) +{ + return std::sqrt(dot(p1, p1)); } template -Array3 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]) }; +Array3 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])}; } template -Array3 axp(const Point3D &p1, const F alpha, const Point3D &p2) { - return { p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], - p1[2] + alpha * p2[2] }; +Array3 axp(const Point3D &p1, const F alpha, const Point3D &p2) +{ + return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]}; } /** @@ -47,11 +51,11 @@ Array3 axp(const Point3D &p1, const F alpha, const Point3D &p2) { * \return Resistance of the medium from the segment to the electrode. */ template -F point_source_lfp_factor(const Point3D &e_pos, const Point3D &seg_pos, - const F radius, const F f) { - nrn_assert(radius >= 0.0); - Array3 es = axp(e_pos, -1.0, seg_pos); - return f / std::max(norm(es), radius); +F point_source_lfp_factor(const Point3D &e_pos, const Point3D &seg_pos, const F radius, const F f) +{ + nrn_assert(radius >= 0.0); + Array3 es = axp(e_pos, -1.0, seg_pos); + return f / std::max(norm(es), radius); } /** @@ -65,66 +69,68 @@ F point_source_lfp_factor(const Point3D &e_pos, const Point3D &seg_pos, * \return Resistance of the medium from the segment to the electrode. */ template -F line_source_lfp_factor(const Point3D &e_pos, const Point3D &seg_0, - const Point3D &seg_1, const F radius, const F f) { - nrn_assert(radius >= 0.0); - Array3 dx = axp(seg_1, -1.0, seg_0); - Array3 de = axp(e_pos, -1.0, seg_0); - F dx2(dot, F>(dx, dx)); - F dxn(std::sqrt(dx2)); - if (dxn < 1.0e-20) { - return point_source_lfp_factor(e_pos, seg_0, radius, f); - } - F de2(dot, F>(de, de)); - F mu(dot, F>(dx, de) / dx2); - Array3 de_star(axp(de, -mu, dx)); - F de_star2(dot, F>(de_star, de_star)); - F q2(de_star2 / dx2); - - F delta(mu * mu - (de2 - radius * radius) / dx2); - F one_m_mu(1.0 - mu); - auto log_integral = [&q2, &dxn](F a, F b) { - if (std::abs(q2) < std::numeric_limits::epsilon()) { - if (a * b <= 0) { - throw std::invalid_argument( - "Log integral: invalid arguments " + std::to_string(b) + " " + - std::to_string(a) + ". Likely electrode exactly on the segment and " - "no flooring is present."); - } - 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 { - F sqr_delta(std::sqrt(delta)); - F d1(mu - sqr_delta); - F d2(mu + sqr_delta); - F parts = 0.0; - if (d1 > 0.0) { - F b(std::min(d1, 1.0) - mu); - parts += log_integral(-mu, b); +F line_source_lfp_factor(const Point3D &e_pos, + const Point3D &seg_0, + const Point3D &seg_1, + const F radius, + const F f) +{ + nrn_assert(radius >= 0.0); + Array3 dx = axp(seg_1, -1.0, seg_0); + Array3 de = axp(e_pos, -1.0, seg_0); + F dx2(dot, F>(dx, dx)); + F dxn(std::sqrt(dx2)); + if (dxn < 1.0e-20) { + return point_source_lfp_factor(e_pos, seg_0, radius, f); } - if (d2 < 1.0) { - F b(std::max(d2, 0.0) - mu); - parts += log_integral(b, one_m_mu); + F de2(dot, F>(de, de)); + F mu(dot, F>(dx, de) / dx2); + Array3 de_star(axp(de, -mu, dx)); + F de_star2(dot, F>(de_star, de_star)); + F q2(de_star2 / dx2); + + F delta(mu * mu - (de2 - radius * radius) / dx2); + F one_m_mu(1.0 - mu); + auto log_integral = [&q2, &dxn](F a, F b) { + if (std::abs(q2) < std::numeric_limits::epsilon()) { + if (a * b <= 0) { + throw std::invalid_argument("Log integral: invalid arguments " + std::to_string(b) + + " " + std::to_string(a) + + ". Likely electrode exactly on the segment and " + "no flooring is present."); + } + 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 { + F sqr_delta(std::sqrt(delta)); + F d1(mu - sqr_delta); + F d2(mu + sqr_delta); + F parts = 0.0; + if (d1 > 0.0) { + F b(std::min(d1, 1.0) - mu); + parts += log_integral(-mu, b); + } + if (d2 < 1.0) { + F b(std::max(d2, 0.0) - mu); + parts += log_integral(b, one_m_mu); + }; + // complement + F 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; }; - // complement - F 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; - }; } enum LFPCalculatorType { - LineSource, - PointSource + LineSource, + PointSource }; /** @@ -133,95 +139,102 @@ enum LFPCalculatorType { template struct LFPCalculator { - /** - * LFP Calculator constructor - * \tparam Point3Ds A vector of 3D points type - * \tparam Vector A vector of floats type - * \param comm MPI communicator - * \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 - */ - template - LFPCalculator(const Point3Ds &seg_start, const Point3Ds &seg_end, - const 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."); + /** + * LFP Calculator constructor + * \tparam Point3Ds A vector of 3D points type + * \tparam Vector A vector of floats type + * \param comm MPI communicator + * \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 + */ + template + LFPCalculator(const Point3Ds &seg_start, + const Point3Ds &seg_end, + const 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 * boost::math::constants::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); + } + } } - 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 * - boost::math::constants::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 void 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]]; - } - } - std::vector res_reduced(res.size()); + + template + void 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]]; + } + } + std::vector res_reduced(res.size()); #if NRNMPI - int mpi_sum{ 1 }; - nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), - mpi_sum); -#else - res_reduced = res; + int mpi_sum{1}; + nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); +#else + res_reduced = res; #endif - lfp_values = res_reduced; - } + lfp_values = res_reduced; + } - std::vector lfp_values; + std::vector lfp_values; -private: - template - inline F getFactor(const Point3D &e_pos, const Point3D &seg_0, - const Point3D &seg_1, const F radius, const F f) const; + private: + template + inline F getFactor(const Point3D &e_pos, + const Point3D &seg_0, + const Point3D &seg_1, + const F radius, + const F f) const; - std::vector > m; - const std::vector &segment_ids_; + std::vector > m; + const std::vector &segment_ids_; }; template <> template F LFPCalculator::getFactor(const Point3D &e_pos, const Point3D &seg_0, - const Point3D &seg_1, const F radius, - const F f) const { - return line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); + const Point3D &seg_1, + const F radius, + const F f) const +{ + return line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); } template <> template F LFPCalculator::getFactor(const Point3D &e_pos, const Point3D &seg_0, - const Point3D &seg_1, const F radius, - const F f) const { - return point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), - radius, f); + const Point3D &seg_1, + const F radius, + const F f) const +{ + return point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), radius, f); } }; -#endif // AREA_LFP_H +#endif // AREA_LFP_H From 8cd40fa686ebe8bd5615da812025eb44a14a4819 Mon Sep 17 00:00:00 2001 From: Pramod S Kumbhar Date: Tue, 9 Feb 2021 16:27:52 +0100 Subject: [PATCH 11/36] Apply clang-format --- coreneuron/io/lfp.hpp | 80 +++++++++----------- tests/unit/lfp/lfp.cpp | 165 ++++++++++++++++++++--------------------- 2 files changed, 112 insertions(+), 133 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index fbe29fff8..8f99071b8 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -17,26 +17,22 @@ template using Array3 = std::array; template -F dot(const Point3D &p1, const Point3D &p2) -{ +F dot(const Point3D& p1, const Point3D& p2) { return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; } template -F norm(const Point3D &p1) -{ +F norm(const Point3D& p1) { return std::sqrt(dot(p1, p1)); } template -Array3 barycenter(const Point3D &p1, const Point3D &p2) -{ +Array3 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])}; } template -Array3 axp(const Point3D &p1, const F alpha, const Point3D &p2) -{ +Array3 axp(const Point3D& p1, const F alpha, const Point3D& p2) { return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]}; } @@ -51,8 +47,7 @@ Array3 axp(const Point3D &p1, const F alpha, const Point3D &p2) * \return Resistance of the medium from the segment to the electrode. */ template -F point_source_lfp_factor(const Point3D &e_pos, const Point3D &seg_pos, const F radius, const F f) -{ +F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, const F radius, const F f) { nrn_assert(radius >= 0.0); Array3 es = axp(e_pos, -1.0, seg_pos); return f / std::max(norm(es), radius); @@ -69,12 +64,11 @@ F point_source_lfp_factor(const Point3D &e_pos, const Point3D &seg_pos, const F * \return Resistance of the medium from the segment to the electrode. */ template -F line_source_lfp_factor(const Point3D &e_pos, - const Point3D &seg_0, - const Point3D &seg_1, +F line_source_lfp_factor(const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, const F radius, - const F f) -{ + const F f) { nrn_assert(radius >= 0.0); Array3 dx = axp(seg_1, -1.0, seg_0); Array3 de = axp(e_pos, -1.0, seg_0); @@ -128,17 +122,13 @@ F line_source_lfp_factor(const Point3D &e_pos, }; } -enum LFPCalculatorType { - LineSource, - PointSource -}; +enum LFPCalculatorType { LineSource, PointSource }; /** * \brief LFPCalculator allows calculation of LFP given membrane currents. */ template struct LFPCalculator { - /** * LFP Calculator constructor * \tparam Point3Ds A vector of 3D points type @@ -153,14 +143,13 @@ struct LFPCalculator { * medium */ template - LFPCalculator(const Point3Ds &seg_start, - const Point3Ds &seg_end, - const Vector &radius, - const std::vector &segment_ids, - const Point3Ds &electrodes, + LFPCalculator(const Point3Ds& seg_start, + const Point3Ds& seg_end, + const Vector& radius, + const std::vector& segment_ids, + const Point3Ds& electrodes, double extra_cellular_conductivity) - : segment_ids_(segment_ids) - { + : segment_ids_(segment_ids) { if (seg_start.size() != seg_end.size()) { throw std::invalid_argument("Different number of segment starts and ends."); } @@ -171,7 +160,7 @@ struct LFPCalculator { m.resize(electrodes.size()); for (size_t k = 0; k < electrodes.size(); ++k) { - auto &ms = m[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); @@ -180,12 +169,11 @@ struct LFPCalculator { } template - void lfp(const Vector &membrane_current) - { + void 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]; + auto& ms = m[k]; for (size_t l = 0; l < ms.size(); l++) { res[k] += ms[l] * membrane_current[segment_ids_[l]]; } @@ -204,37 +192,35 @@ struct LFPCalculator { private: template - inline F getFactor(const Point3D &e_pos, - const Point3D &seg_0, - const Point3D &seg_1, + inline F getFactor(const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, const F radius, const F f) const; - std::vector > m; - const std::vector &segment_ids_; + std::vector> m; + const std::vector& segment_ids_; }; template <> template -F LFPCalculator::getFactor(const Point3D &e_pos, - const Point3D &seg_0, - const Point3D &seg_1, +F LFPCalculator::getFactor(const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, const F radius, - const F f) const -{ + const F f) const { return line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); } template <> template -F LFPCalculator::getFactor(const Point3D &e_pos, - const Point3D &seg_0, - const Point3D &seg_1, +F LFPCalculator::getFactor(const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, const F radius, - const F f) const -{ + const F f) const { return point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), radius, f); } -}; +}; // namespace coreneuron #endif // AREA_LFP_H diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index d17a5da64..7ed689efb 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -9,103 +9,96 @@ #include "coreneuron/io/lfp.hpp" #include "coreneuron/mpi/nrnmpi.h" -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; +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; } -std::array operator*(const double x, - const std::array &y) { - return { x * y[0], x * y[1], x * y[2] }; +std::array operator*(const double x, const std::array& y) { + return {x * y[0], x * y[1], x * y[2]}; } -std::array operator+(const std::array &x, - const std::array &y) { - return { x[0] + y[0], x[1] + y[1], x[2] + y[2] }; +std::array operator+(const std::array& x, const std::array& y) { + return {x[0] + y[0], x[1] + y[1], x[2] + y[2]}; } BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { - using namespace coreneuron; + using namespace coreneuron; #if NRNMPI - nrnmpi_init(nullptr, nullptr); + 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 = - segment_start + std::array{ 0.0, 0.0, segment_length }; - double floor{ 1.0e-6 }; - double pi{ 3.141592653589 }; + 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 = segment_start + + std::array{0.0, 0.0, segment_length}; + double floor{1.0e-6}; + double 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 = - segment_middle + std::array{ 0.0, 1.0e-5 - k * 1.0e-6, 0.0 }; - std::array circling_elec = - segment_middle + - std::array{ 0.0, - circling_radius * std::cos(2.0 * pi * k / 10), - circling_radius * std::sin(2.0 * pi * k / 10) }; + 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 = + segment_middle + std::array{0.0, 1.0e-5 - k * 1.0e-6, 0.0}; + std::array circling_elec = + segment_middle + 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, double>( - circling_elec + (-x) * segment_start + - (x - 1) * 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]) * 0.5 / - (vals[k] + 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); + 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, double>( + circling_elec + (-x) * segment_start + (x - 1) * 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]) * 0.5 / (vals[k] + 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(); + nrnmpi_finalize(); #endif } From 58fe6f1e3c89d95fac30f976630e9fc46c36732f Mon Sep 17 00:00:00 2001 From: Pramod S Kumbhar Date: Tue, 9 Feb 2021 16:28:36 +0100 Subject: [PATCH 12/36] Apply cmake-format --- tests/unit/lfp/CMakeLists.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/unit/lfp/CMakeLists.txt b/tests/unit/lfp/CMakeLists.txt index 7db0c825e..d7c964be8 100644 --- a/tests/unit/lfp/CMakeLists.txt +++ b/tests/unit/lfp/CMakeLists.txt @@ -1,4 +1,3 @@ - # ============================================================================= # Copyright (C) 2016-2020 Blue Brain Project # @@ -9,7 +8,9 @@ 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} +target_link_libraries( + lfp_test_bin + ${MPI_CXX_LIBRARIES} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} coreneuron ${corenrn_mech_lib} @@ -20,5 +21,4 @@ if(CORENRN_ENABLE_GPU) target_link_libraries(lfp_test_bin ${link_cudacoreneuron} ${CUDA_LIBRARIES}) endif() -add_test(NAME lfp_test COMMAND ${TEST_EXEC_PREFIX} - ${CMAKE_CURRENT_BINARY_DIR}/lfp_test_bin) +add_test(NAME lfp_test COMMAND ${TEST_EXEC_PREFIX} ${CMAKE_CURRENT_BINARY_DIR}/lfp_test_bin) From 06c1519c51eb3e380f55bc81eb1b7db0b93d0f35 Mon Sep 17 00:00:00 2001 From: Pramod S Kumbhar Date: Tue, 9 Feb 2021 17:09:06 +0100 Subject: [PATCH 13/36] Fix the dependency issue : lfp should be compiled after corenrnmech library is built --- tests/unit/lfp/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/unit/lfp/CMakeLists.txt b/tests/unit/lfp/CMakeLists.txt index d7c964be8..c20c61a3a 100644 --- a/tests/unit/lfp/CMakeLists.txt +++ b/tests/unit/lfp/CMakeLists.txt @@ -21,4 +21,6 @@ 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) From 60bfc661fa72bc33fea301d07c73eb6ed16de6a2 Mon Sep 17 00:00:00 2001 From: Ioannis Magkanaris Date: Wed, 10 Feb 2021 17:40:55 +0100 Subject: [PATCH 14/36] Fixes related to PI value read by NEURON - Created new pi global parameter which is read by NEURON globals.dat or directly NEURON in online mode - Use this pi in lfp calculation and test - Removed a deprecated boost test header file --- coreneuron/io/global_vars.cpp | 1 + coreneuron/io/lfp.hpp | 5 ++--- coreneuron/mechanism/register_mech.cpp | 3 ++- coreneuron/nrnconf.h | 1 + tests/unit/alignment/alignment.cpp | 1 - tests/unit/lfp/lfp.cpp | 6 +++--- tests/unit/queueing/test_queueing.cpp | 1 - 7 files changed, 9 insertions(+), 9 deletions(-) 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.hpp b/coreneuron/io/lfp.hpp index 8f99071b8..c34c5b788 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -6,8 +6,7 @@ #include #include -#include - +#include "coreneuron/nrnconf.h" #include "coreneuron/utils/nrn_assert.h" #include "coreneuron/mpi/nrnmpidec.h" @@ -156,7 +155,7 @@ struct LFPCalculator { 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 * boost::math::constants::pi())); + double f(1.0 / (extra_cellular_conductivity * 4.0 * pi)); m.resize(electrodes.size()); for (size_t k = 0; k < electrodes.size(); ++k) { 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/unit/alignment/alignment.cpp b/tests/unit/alignment/alignment.cpp index 0bbf2ec09..4d5e0acea 100644 --- a/tests/unit/alignment/alignment.cpp +++ b/tests/unit/alignment/alignment.cpp @@ -14,7 +14,6 @@ #include #include -#include #include #include "coreneuron/utils/memory.h" diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 7ed689efb..3d73479ff 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -4,11 +4,12 @@ #include #include -#include #include "coreneuron/io/lfp.hpp" #include "coreneuron/mpi/nrnmpi.h" +using namespace coreneuron; + template double integral(F f, double a, double b, int n) { double step = (b - a) / n; // width of each small rectangle @@ -28,7 +29,6 @@ std::array operator+(const std::array& x, const std::array } BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { - using namespace coreneuron; #if NRNMPI nrnmpi_init(nullptr, nullptr); #endif @@ -38,7 +38,7 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { std::array segment_end = segment_start + std::array{0.0, 0.0, segment_length}; double floor{1.0e-6}; - double pi{3.141592653589}; + pi = 3.141592653589; std::array vals; double circling_radius{1.0e-6}; diff --git a/tests/unit/queueing/test_queueing.cpp b/tests/unit/queueing/test_queueing.cpp index 21ae3297f..830deebd0 100644 --- a/tests/unit/queueing/test_queueing.cpp +++ b/tests/unit/queueing/test_queueing.cpp @@ -10,7 +10,6 @@ #define TYPE T::cont #include -#include #include #include From 91cf705a66ef862bf93916376caba49d7928155e Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 15 Feb 2021 09:58:55 +0100 Subject: [PATCH 15/36] Refactoring following review. --- tests/unit/lfp/lfp.cpp | 142 ++++++++++++++++++++--------------------- 1 file changed, 71 insertions(+), 71 deletions(-) diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 3d73479ff..daaf4d345 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -20,84 +20,84 @@ double integral(F f, double a, double b, int n) { return area; } -std::array operator*(const double x, const std::array& y) { - return {x * y[0], x * y[1], x * y[2]}; -} - -std::array operator+(const std::array& x, const std::array& y) { - return {x[0] + y[0], x[1] + y[1], x[2] + y[2]}; -} 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 = segment_start + - std::array{0.0, 0.0, segment_length}; - double floor{1.0e-6}; - pi = 3.141592653589; + 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 = + paxp(segment_start, 1.0, std::array{ 0.0, 0.0, segment_length }); + double floor{ 1.0e-6 }; + double 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 = - segment_middle + std::array{0.0, 1.0e-5 - k * 1.0e-6, 0.0}; - std::array circling_elec = - segment_middle + std::array{0.0, - circling_radius * std::cos(2.0 * pi * k / 10), - circling_radius * std::sin(2.0 * pi * k / 10)}; + 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 = + paxp(segment_middle, 1.0, std::array{ 0.0, 1.0e-5 - k * 1.0e-6, 0.0 }); + std::array circling_elec = + paxp(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, double>( - circling_elec + (-x) * segment_start + (x - 1) * 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]) * 0.5 / (vals[k] + 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); + 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( + paxp(circling_elec, -1.0, + paxp(segment_end, x, paxp(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]) * 0.5 / + (vals[k] + 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 From fb91048701831085b3058bef15220a4dd41ec92a Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 15 Feb 2021 10:47:05 +0100 Subject: [PATCH 16/36] Modify the relative criterion. --- tests/unit/lfp/lfp.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index daaf4d345..36aa5b02c 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -75,8 +75,8 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { } // TEST of SYMMETRY of LFP FORMULA for (size_t k = 0; k < 5; k++) { - BOOST_REQUIRE(std::abs((vals[k] - vals[k + 5]) * 0.5 / - (vals[k] + vals[k + 5])) < 1.0e-12); + 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 } From 5c6e6429d2ec76e17f1941caf46031f721bbd505 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 15 Feb 2021 11:30:47 +0100 Subject: [PATCH 17/36] Refactoring following review. --- coreneuron/io/lfp.hpp | 44 ++++++++++++++++++------------------------- 1 file changed, 18 insertions(+), 26 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index c34c5b788..33309f4c7 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -12,26 +12,23 @@ namespace coreneuron { -template -using Array3 = std::array; +using F = double; +using Point3D = std::array; +using Point3Ds = std::vector; -template F dot(const Point3D& p1, const Point3D& p2) { return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; } -template F norm(const Point3D& p1) { - return std::sqrt(dot(p1, p1)); + return std::sqrt(dot(p1, p1)); } -template -Array3 barycenter(const Point3D& p1, const Point3D& p2) { +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])}; } -template -Array3 axp(const Point3D& p1, const F alpha, const Point3D& p2) { +Point3D paxp(const Point3D& p1, const F alpha, const Point3D& p2) { return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]}; } @@ -45,11 +42,10 @@ Array3 axp(const Point3D& p1, const F alpha, const Point3D& p2) { * \param f conductivity factor 1/([4 pi] * [conductivity]) * \return Resistance of the medium from the segment to the electrode. */ -template F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, const F radius, const F f) { nrn_assert(radius >= 0.0); - Array3 es = axp(e_pos, -1.0, seg_pos); - return f / std::max(norm(es), radius); + Point3D es = paxp(e_pos, -1.0, seg_pos); + return f / std::max(norm(es), radius); } /** @@ -62,24 +58,23 @@ F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, const F * \param f conductivity factor 1/([4 pi] * [conductivity]) * \return Resistance of the medium from the segment to the electrode. */ -template F line_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, const F radius, const F f) { nrn_assert(radius >= 0.0); - Array3 dx = axp(seg_1, -1.0, seg_0); - Array3 de = axp(e_pos, -1.0, seg_0); - F dx2(dot, F>(dx, dx)); + Point3D dx = paxp(seg_1, -1.0, seg_0); + Point3D de = paxp(e_pos, -1.0, seg_0); + F dx2(dot(dx, dx)); F dxn(std::sqrt(dx2)); if (dxn < 1.0e-20) { return point_source_lfp_factor(e_pos, seg_0, radius, f); } - F de2(dot, F>(de, de)); - F mu(dot, F>(dx, de) / dx2); - Array3 de_star(axp(de, -mu, dx)); - F de_star2(dot, F>(de_star, de_star)); + F de2(dot(de, de)); + F mu(dot(dx, de) / dx2); + Point3D de_star(paxp(de, -mu, dx)); + F de_star2(dot(de_star, de_star)); F q2(de_star2 / dx2); F delta(mu * mu - (de2 - radius * radius) / dx2); @@ -141,10 +136,9 @@ struct LFPCalculator { * \param extra_cellular_conductivity conductivity of the extra-cellular * medium */ - template LFPCalculator(const Point3Ds& seg_start, const Point3Ds& seg_end, - const Vector& radius, + const std::vector& radius, const std::vector& segment_ids, const Point3Ds& electrodes, double extra_cellular_conductivity) @@ -190,7 +184,6 @@ struct LFPCalculator { std::vector lfp_values; private: - template inline F getFactor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, @@ -202,7 +195,6 @@ struct LFPCalculator { }; template <> -template F LFPCalculator::getFactor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, @@ -212,14 +204,14 @@ F LFPCalculator::getFactor(const Point3D& e_pos, } template <> -template F LFPCalculator::getFactor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, const F radius, const F f) const { - return point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), radius, f); + return point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), radius, f); } + }; // namespace coreneuron #endif // AREA_LFP_H From 38741008b75616adc50bc60b575d0f6183fbb5c0 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 15 Feb 2021 12:22:14 +0100 Subject: [PATCH 18/36] Unnecessary abs. --- coreneuron/io/lfp.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 33309f4c7..3b1fcabaa 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -80,7 +80,7 @@ F line_source_lfp_factor(const Point3D& e_pos, F delta(mu * mu - (de2 - radius * radius) / dx2); F one_m_mu(1.0 - mu); auto log_integral = [&q2, &dxn](F a, F b) { - if (std::abs(q2) < std::numeric_limits::epsilon()) { + if (q2 < std::numeric_limits::epsilon()) { if (a * b <= 0) { throw std::invalid_argument("Log integral: invalid arguments " + std::to_string(b) + " " + std::to_string(a) + From b05d041ea2803a50306a5e6904ff00db4a8057bd Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 15 Feb 2021 13:19:58 +0100 Subject: [PATCH 19/36] Renaming paxp to paxpy. --- coreneuron/io/lfp.hpp | 10 +++++----- tests/unit/lfp/lfp.cpp | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 3b1fcabaa..014c22a9a 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -28,7 +28,7 @@ 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])}; } -Point3D paxp(const Point3D& p1, const F alpha, const Point3D& p2) { +Point3D paxpy(const Point3D& p1, const F alpha, const Point3D& p2) { return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]}; } @@ -44,7 +44,7 @@ Point3D paxp(const Point3D& p1, const F alpha, const Point3D& p2) { */ F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, const F radius, const F f) { nrn_assert(radius >= 0.0); - Point3D es = paxp(e_pos, -1.0, seg_pos); + Point3D es = paxpy(e_pos, -1.0, seg_pos); return f / std::max(norm(es), radius); } @@ -64,8 +64,8 @@ F line_source_lfp_factor(const Point3D& e_pos, const F radius, const F f) { nrn_assert(radius >= 0.0); - Point3D dx = paxp(seg_1, -1.0, seg_0); - Point3D de = paxp(e_pos, -1.0, seg_0); + Point3D dx = paxpy(seg_1, -1.0, seg_0); + Point3D de = paxpy(e_pos, -1.0, seg_0); F dx2(dot(dx, dx)); F dxn(std::sqrt(dx2)); if (dxn < 1.0e-20) { @@ -73,7 +73,7 @@ F line_source_lfp_factor(const Point3D& e_pos, } F de2(dot(de, de)); F mu(dot(dx, de) / dx2); - Point3D de_star(paxp(de, -mu, dx)); + Point3D de_star(paxpy(de, -mu, dx)); F de_star2(dot(de_star, de_star)); F q2(de_star2 / dx2); diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 36aa5b02c..1d7280c3a 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -30,7 +30,7 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { std::array segment_start = std::array{ 0.0, 0.0, segment_start_val }; std::array segment_end = - paxp(segment_start, 1.0, std::array{ 0.0, 0.0, segment_length }); + paxpy(segment_start, 1.0, std::array{ 0.0, 0.0, segment_length }); double floor{ 1.0e-6 }; double pi = 3.141592653589; @@ -40,9 +40,9 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { double medium_resistivity_fac{ 1.0 }; for (auto k = 0; k < 10; k++) { std::array approaching_elec = - paxp(segment_middle, 1.0, std::array{ 0.0, 1.0e-5 - k * 1.0e-6, 0.0 }); + paxpy(segment_middle, 1.0, std::array{ 0.0, 1.0e-5 - k * 1.0e-6, 0.0 }); std::array circling_elec = - paxp(segment_middle, 1.0, + 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) }); @@ -58,8 +58,8 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { return 1.0 / std::max(floor, norm( - paxp(circling_elec, -1.0, - paxp(segment_end, x, paxp(segment_start, + paxpy(circling_elec, -1.0, + paxpy(segment_end, x, paxpy(segment_start, -1.0, segment_end))))); }, 0.0, 1.0, 10000); From c4240e8e30590999e87c4a1b0c86187c7852efa8 Mon Sep 17 00:00:00 2001 From: Pramod S Kumbhar Date: Mon, 15 Feb 2021 15:03:04 +0100 Subject: [PATCH 20/36] do cmake & clang format --- tests/unit/lfp/lfp.cpp | 140 ++++++++++++++++++++--------------------- 1 file changed, 69 insertions(+), 71 deletions(-) diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index 1d7280c3a..a9df0386b 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -25,79 +25,77 @@ 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 }; - double pi = 3.141592653589; + 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}; + double 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) }); + 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); + 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 From ccd536f20378d4fb9210f34a0e058a1e464d8356 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 15 Feb 2021 16:33:51 +0100 Subject: [PATCH 21/36] More changes following review. --- coreneuron/io/lfp.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 014c22a9a..021d27726 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -68,7 +68,7 @@ F line_source_lfp_factor(const Point3D& e_pos, Point3D de = paxpy(e_pos, -1.0, seg_0); F dx2(dot(dx, dx)); F dxn(std::sqrt(dx2)); - if (dxn < 1.0e-20) { + if (dxn < std::numeric_limits::epsilon()) { return point_source_lfp_factor(e_pos, seg_0, radius, f); } F de2(dot(de, de)); @@ -149,7 +149,7 @@ struct LFPCalculator { 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)); + double f(1.0 / (extra_cellular_conductivity * 4.0 * 3.141592653589)); m.resize(electrodes.size()); for (size_t k = 0; k < electrodes.size(); ++k) { From 1795ca305ba7583c0b34ff361367666209772b4a Mon Sep 17 00:00:00 2001 From: Baudouin Del Marmol Date: Mon, 22 Feb 2021 11:58:04 +0100 Subject: [PATCH 22/36] Refactoring following review - part 1. --- coreneuron/io/lfp.cpp | 122 +++++++++++++++++++++++++++++++++++++++ coreneuron/io/lfp.hpp | 126 ++++++----------------------------------- tests/unit/lfp/lfp.cpp | 1 + 3 files changed, 141 insertions(+), 108 deletions(-) create mode 100644 coreneuron/io/lfp.cpp diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp new file mode 100644 index 000000000..f667be091 --- /dev/null +++ b/coreneuron/io/lfp.cpp @@ -0,0 +1,122 @@ +#include "coreneuron/io/lfp.hpp" + + +namespace coreneuron { + + namespace utils { + + F line_source_lfp_factor(const Point3D& e_pos, + const Point3D& seg_0, + const Point3D& seg_1, + const F radius, + const F f) { + nrn_assert(radius >= 0.0); + Point3D dx = paxpy(seg_1, -1.0, seg_0); + Point3D de = paxpy(e_pos, -1.0, seg_0); + F dx2(dot(dx, dx)); + F dxn(std::sqrt(dx2)); + if (dxn < std::numeric_limits::epsilon()) { + return point_source_lfp_factor(e_pos, seg_0, radius, f); + } + F de2(dot(de, de)); + F mu(dot(dx, de) / dx2); + Point3D de_star(paxpy(de, -mu, dx)); + F de_star2(dot(de_star, de_star)); + F q2(de_star2 / dx2); + + F delta(mu * mu - (de2 - radius * radius) / dx2); + F one_m_mu(1.0 - mu); + auto log_integral = [&q2, &dxn](F a, F b) { + if (q2 < std::numeric_limits::epsilon()) { + if (a * b <= 0) { + throw std::invalid_argument("Log integral: invalid arguments " + std::to_string(b) + + " " + std::to_string(a) + + ". Likely electrode exactly on the segment and " + "no flooring is present."); + } + 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 { + F sqr_delta(std::sqrt(delta)); + F d1(mu - sqr_delta); + F d2(mu + sqr_delta); + F parts = 0.0; + if (d1 > 0.0) { + F b(std::min(d1, 1.0) - mu); + parts += log_integral(-mu, b); + } + if (d2 < 1.0) { + F b(std::max(d2, 0.0) - mu); + parts += log_integral(b, one_m_mu); + }; + // complement + F 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 utils + + + using namespace utils; + + 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 * 3.141592653589)); + + 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]]; + } + } + std::vector res_reduced(res.size()); +#if NRNMPI + int mpi_sum{1}; + nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); +#else + res_reduced = res; +#endif + lfp_values = res_reduced; + } + + +template void LFPCalculator::lfp(const double*& membrane_current); +template void LFPCalculator::lfp(const double*& membrane_current); + + +} // namespace coreneuron diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 021d27726..79137e9c1 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -1,9 +1,7 @@ -#ifndef LFP_H -#define LFP_H +#pragma once #include #include -#include #include #include "coreneuron/nrnconf.h" @@ -12,37 +10,37 @@ namespace coreneuron { +namespace utils { + using F = double; using Point3D = std::array; using Point3Ds = std::vector; -F dot(const Point3D& p1, const Point3D& p2) { +inline F dot(const Point3D& p1, const Point3D& p2) { return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; } -F norm(const Point3D& p1) { +inline F norm(const Point3D& p1) { return std::sqrt(dot(p1, p1)); } -Point3D barycenter(const Point3D& p1, const Point3D& p2) { +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])}; } -Point3D paxpy(const Point3D& p1, const F alpha, const Point3D& p2) { +inline Point3D paxpy(const Point3D& p1, const F alpha, const Point3D& p2) { return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]}; } /** * - * \tparam Point3D Type of 3D point - * \tparam F Floating point type * \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. */ -F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, const F radius, const F f) { +inline F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, const F radius, const F f) { nrn_assert(radius >= 0.0); Point3D es = paxpy(e_pos, -1.0, seg_pos); return f / std::max(norm(es), radius); @@ -50,71 +48,18 @@ F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, const F /** * - * \tparam Point3D Type of 3D point - * \tparam F Floating point type * \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. */ -F line_source_lfp_factor(const Point3D& e_pos, +inline F line_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, const F radius, - const F f) { - nrn_assert(radius >= 0.0); - Point3D dx = paxpy(seg_1, -1.0, seg_0); - Point3D de = paxpy(e_pos, -1.0, seg_0); - F dx2(dot(dx, dx)); - F dxn(std::sqrt(dx2)); - if (dxn < std::numeric_limits::epsilon()) { - return point_source_lfp_factor(e_pos, seg_0, radius, f); - } - F de2(dot(de, de)); - F mu(dot(dx, de) / dx2); - Point3D de_star(paxpy(de, -mu, dx)); - F de_star2(dot(de_star, de_star)); - F q2(de_star2 / dx2); - - F delta(mu * mu - (de2 - radius * radius) / dx2); - F one_m_mu(1.0 - mu); - auto log_integral = [&q2, &dxn](F a, F b) { - if (q2 < std::numeric_limits::epsilon()) { - if (a * b <= 0) { - throw std::invalid_argument("Log integral: invalid arguments " + std::to_string(b) + - " " + std::to_string(a) + - ". Likely electrode exactly on the segment and " - "no flooring is present."); - } - 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 { - F sqr_delta(std::sqrt(delta)); - F d1(mu - sqr_delta); - F d2(mu + sqr_delta); - F parts = 0.0; - if (d1 > 0.0) { - F b(std::min(d1, 1.0) - mu); - parts += log_integral(-mu, b); - } - if (d2 < 1.0) { - F b(std::max(d2, 0.0) - mu); - parts += log_integral(b, one_m_mu); - }; - // complement - F 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; - }; -} + const F f); +} // namespace utils enum LFPCalculatorType { LineSource, PointSource }; @@ -125,8 +70,6 @@ template struct LFPCalculator { /** * LFP Calculator constructor - * \tparam Point3Ds A vector of 3D points type - * \tparam Vector A vector of floats type * \param comm MPI communicator * \param seg_start all segments start owned by the proc * \param seg_end all segments end owned by the proc @@ -141,45 +84,10 @@ struct LFPCalculator { 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 * 3.141592653589)); - - 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); - } - } - } + double extra_cellular_conductivity); template - void 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]]; - } - } - std::vector res_reduced(res.size()); -#if NRNMPI - int mpi_sum{1}; - nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); -#else - res_reduced = res; -#endif - lfp_values = res_reduced; - } + inline void lfp(const Vector& membrane_current); std::vector lfp_values; @@ -200,7 +108,7 @@ F LFPCalculator::getFactor(const Point3D& e_pos, const Point3D& seg_1, const F radius, const F f) const { - return line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); + return utils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); } template <> @@ -209,9 +117,11 @@ F LFPCalculator::getFactor(const Point3D& e_pos, const Point3D& seg_1, const F radius, const F f) const { - return point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), radius, f); + return utils::point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), radius, f); } +template void LFPCalculator::lfp(const double*& membrane_current); +template void LFPCalculator::lfp(const double*& membrane_current); + }; // namespace coreneuron -#endif // AREA_LFP_H diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index a9df0386b..ddf76b8bd 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -9,6 +9,7 @@ #include "coreneuron/mpi/nrnmpi.h" using namespace coreneuron; +using namespace coreneuron::utils; template double integral(F f, double a, double b, int n) { From 833c22c977096f6e8b5f5a2f0ca932753ae3dcba Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 22 Feb 2021 14:12:46 +0100 Subject: [PATCH 23/36] Splitting definition-declaration + changes following review. --- coreneuron/io/lfp.cpp | 36 +++++++++++++++++++---- coreneuron/io/lfp.hpp | 67 ++++++++++++++++++++++++------------------ tests/unit/lfp/lfp.cpp | 2 +- 3 files changed, 70 insertions(+), 35 deletions(-) diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp index f667be091..c6c7ee024 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -1,9 +1,13 @@ #include "coreneuron/io/lfp.hpp" +#include +#include + + namespace coreneuron { - namespace utils { + namespace lfputils { F line_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_0, @@ -64,8 +68,7 @@ namespace coreneuron { } } // namespace utils - - using namespace utils; + using namespace lfputils; template LFPCalculator::LFPCalculator(const Point3Ds& seg_start, @@ -114,9 +117,30 @@ template lfp_values = res_reduced; } +template <> +lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, + const lfputils::Point3D& seg_0, + const lfputils::Point3D& seg_1, + const lfputils::F radius, + const lfputils::F f) const { + return lfputils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); +} -template void LFPCalculator::lfp(const double*& membrane_current); -template void LFPCalculator::lfp(const double*& membrane_current); - +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 Vec& membrane_current); +template void LFPCalculator::lfp(const Vec& 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 index 79137e9c1..9a78a6ae7 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -1,6 +1,5 @@ #pragma once -#include #include #include @@ -10,11 +9,12 @@ namespace coreneuron { -namespace utils { +namespace lfputils { using F = double; using Point3D = std::array; using Point3Ds = std::vector; +using Vec = double*; inline F dot(const Point3D& p1, const Point3D& p2) { return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; @@ -54,7 +54,7 @@ inline F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, c * \param f conductivity factor 1/([4 pi] * [conductivity]) * \return Resistance of the medium from the segment to the electrode. */ -inline F line_source_lfp_factor(const Point3D& e_pos, +F line_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, const F radius, @@ -79,49 +79,60 @@ struct LFPCalculator { * \param extra_cellular_conductivity conductivity of the extra-cellular * medium */ - LFPCalculator(const Point3Ds& seg_start, - const Point3Ds& seg_end, + LFPCalculator(const lfputils::Point3Ds& seg_start, + const lfputils::Point3Ds& seg_end, const std::vector& radius, const std::vector& segment_ids, - const Point3Ds& electrodes, + const lfputils::Point3Ds& electrodes, double extra_cellular_conductivity); template - inline void lfp(const Vector& membrane_current); + void lfp(const Vector& membrane_current); std::vector lfp_values; private: - inline F getFactor(const Point3D& e_pos, - const Point3D& seg_0, - const Point3D& seg_1, - const F radius, - const F f) const; + inline lfputils::F getFactor(const lfputils::Point3D& e_pos, + const lfputils::Point3D& seg_0, + const lfputils::Point3D& seg_1, + const lfputils::F radius, + const lfputils::F f) const; std::vector> m; const std::vector& segment_ids_; }; template <> -F LFPCalculator::getFactor(const Point3D& e_pos, - const Point3D& seg_0, - const Point3D& seg_1, - const F radius, - const F f) const { - return utils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); -} +lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, + const lfputils::Point3D& seg_0, + const lfputils::Point3D& seg_1, + const lfputils::F radius, + const lfputils::F f) const; template <> -F LFPCalculator::getFactor(const Point3D& e_pos, - const Point3D& seg_0, - const Point3D& seg_1, - const F radius, - const F f) const { - return utils::point_source_lfp_factor(e_pos, barycenter(seg_0, seg_1), radius, f); +lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, + const lfputils::Point3D& seg_0, + const lfputils::Point3D& seg_1, + const lfputils::F radius, + const lfputils::F f) const { + return lfputils::point_source_lfp_factor(e_pos, lfputils::barycenter(seg_0, seg_1), radius, f); } -template void LFPCalculator::lfp(const double*& membrane_current); -template void LFPCalculator::lfp(const double*& membrane_current); - +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::Vec& membrane_current); +extern template void LFPCalculator::lfp(const lfputils::Vec& 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/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index ddf76b8bd..d51a64d44 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -9,7 +9,7 @@ #include "coreneuron/mpi/nrnmpi.h" using namespace coreneuron; -using namespace coreneuron::utils; +using namespace coreneuron::lfputils; template double integral(F f, double a, double b, int n) { From 75d902b2495bc14efb60d760e8d615eafc370a9c Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Mon, 22 Feb 2021 14:24:38 +0100 Subject: [PATCH 24/36] More amendments following review. --- coreneuron/io/lfp.cpp | 18 ++++++------------ coreneuron/io/lfp.hpp | 4 +++- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp index c6c7ee024..fc10ca5ef 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -2,7 +2,7 @@ #include #include - +#include namespace coreneuron { @@ -33,10 +33,12 @@ namespace coreneuron { auto log_integral = [&q2, &dxn](F a, F b) { if (q2 < std::numeric_limits::epsilon()) { if (a * b <= 0) { - throw std::invalid_argument("Log integral: invalid arguments " + std::to_string(b) + - " " + std::to_string(a) + + std::ostringstream s; + s << "Log integral: invalid arguments " << b << + " " << a << ". Likely electrode exactly on the segment and " - "no flooring is present."); + << "no flooring is present."; + throw std::invalid_argument(s.str()); } return std::abs(std::log(b / a)) / dxn; } else { @@ -117,14 +119,6 @@ template lfp_values = res_reduced; } -template <> -lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, - const lfputils::Point3D& seg_0, - const lfputils::Point3D& seg_1, - const lfputils::F radius, - const lfputils::F f) const { - return lfputils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); -} template LFPCalculator::LFPCalculator(const lfputils::Point3Ds& seg_start, const lfputils::Point3Ds& seg_end, diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 9a78a6ae7..4ae3ea0cd 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -107,7 +107,9 @@ lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, const lfputils::Point3D& seg_0, const lfputils::Point3D& seg_1, const lfputils::F radius, - const lfputils::F f) const; + const lfputils::F f) const { + return lfputils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); +} template <> lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, From 4361c727dd6767fdbf8647fa51fb860492bc70bc Mon Sep 17 00:00:00 2001 From: Pramod S Kumbhar Date: Mon, 22 Feb 2021 15:21:57 +0100 Subject: [PATCH 25/36] clang-format --- coreneuron/io/lfp.cpp | 121 +++++++++++++++++++++--------------------- coreneuron/io/lfp.hpp | 54 ++++++++++--------- 2 files changed, 88 insertions(+), 87 deletions(-) diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp index fc10ca5ef..30bcd3996 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -7,9 +7,9 @@ namespace coreneuron { - namespace lfputils { +namespace lfputils { - F line_source_lfp_factor(const Point3D& e_pos, +F line_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, const F radius, @@ -34,10 +34,9 @@ namespace coreneuron { 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."; + 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; @@ -67,74 +66,74 @@ namespace coreneuron { } return f * parts; }; - } - } // namespace utils - - using namespace lfputils; +} +} // 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 * 3.141592653589)); +using namespace lfputils; - 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 +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 * 3.141592653589)); + + 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]]; - } +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]]; } - std::vector res_reduced(res.size()); + } + std::vector res_reduced(res.size()); #if NRNMPI - int mpi_sum{1}; - nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); + int mpi_sum{1}; + nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); #else - res_reduced = res; + res_reduced = res; #endif - lfp_values = res_reduced; - } + lfp_values = res_reduced; +} 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); + 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 Vec& membrane_current); -template void LFPCalculator::lfp(const Vec& membrane_current); -template void LFPCalculator::lfp(const std::vector& membrane_current); -template void LFPCalculator::lfp(const std::vector& membrane_current); + 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 Vec& membrane_current); +template void LFPCalculator::lfp(const Vec& membrane_current); +template void LFPCalculator::lfp(const std::vector& membrane_current); +template void LFPCalculator::lfp(const std::vector& membrane_current); -} // namespace coreneuron +} // namespace coreneuron diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 4ae3ea0cd..a09794221 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -40,7 +40,10 @@ inline Point3D paxpy(const Point3D& p1, const F alpha, const Point3D& p2) { * \param f conductivity factor 1/([4 pi] * [conductivity]) * \return Resistance of the medium from the segment to the electrode. */ -inline F point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, const F radius, const F f) { +inline F point_source_lfp_factor(const Point3D& e_pos, + const Point3D& seg_pos, + const F radius, + const F f) { nrn_assert(radius >= 0.0); Point3D es = paxpy(e_pos, -1.0, seg_pos); return f / std::max(norm(es), radius); @@ -59,7 +62,7 @@ F line_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_1, const F radius, const F f); -} // namespace utils +} // namespace lfputils enum LFPCalculatorType { LineSource, PointSource }; @@ -87,16 +90,16 @@ struct LFPCalculator { double extra_cellular_conductivity); template - void lfp(const Vector& membrane_current); + void lfp(const Vector& membrane_current); std::vector lfp_values; private: inline lfputils::F getFactor(const lfputils::Point3D& e_pos, - const lfputils::Point3D& seg_0, - const lfputils::Point3D& seg_1, - const lfputils::F radius, - const lfputils::F f) const; + const lfputils::Point3D& seg_0, + const lfputils::Point3D& seg_1, + const lfputils::F radius, + const lfputils::F f) const; std::vector> m; const std::vector& segment_ids_; @@ -104,37 +107,36 @@ struct LFPCalculator { template <> lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, - const lfputils::Point3D& seg_0, - const lfputils::Point3D& seg_1, - const lfputils::F radius, - const lfputils::F f) const { + const lfputils::Point3D& seg_0, + const lfputils::Point3D& seg_1, + const lfputils::F radius, + const lfputils::F f) const { return lfputils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); } template <> lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, - const lfputils::Point3D& seg_0, - const lfputils::Point3D& seg_1, - const lfputils::F radius, - const lfputils::F f) const { + const lfputils::Point3D& seg_0, + const lfputils::Point3D& seg_1, + const lfputils::F radius, + const lfputils::F 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); + 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); + 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::Vec& membrane_current); extern template void LFPCalculator::lfp(const lfputils::Vec& membrane_current); extern template void LFPCalculator::lfp(const std::vector& membrane_current); extern template void LFPCalculator::lfp(const std::vector& membrane_current); }; // namespace coreneuron - From 72d7abed4cf445458111fb13cc2d5c8e78aaa025 Mon Sep 17 00:00:00 2001 From: Baudouin Date: Wed, 24 Feb 2021 08:52:40 +0100 Subject: [PATCH 26/36] Remove T type definition. --- coreneuron/io/lfp.cpp | 46 +++++++++++++++++++++--------------------- coreneuron/io/lfp.hpp | 47 +++++++++++++++++++++---------------------- 2 files changed, 46 insertions(+), 47 deletions(-) diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp index 30bcd3996..00716fda0 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -9,29 +9,29 @@ namespace coreneuron { namespace lfputils { -F line_source_lfp_factor(const Point3D& e_pos, +double line_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, - const F radius, - const F f) { + 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); - F dx2(dot(dx, dx)); - F dxn(std::sqrt(dx2)); - if (dxn < std::numeric_limits::epsilon()) { + 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); } - F de2(dot(de, de)); - F mu(dot(dx, de) / dx2); + double de2(dot(de, de)); + double mu(dot(dx, de) / dx2); Point3D de_star(paxpy(de, -mu, dx)); - F de_star2(dot(de_star, de_star)); - F q2(de_star2 / dx2); + double de_star2(dot(de_star, de_star)); + double q2(de_star2 / dx2); - F delta(mu * mu - (de2 - radius * radius) / dx2); - F one_m_mu(1.0 - mu); - auto log_integral = [&q2, &dxn](F a, F b) { - if (q2 < std::numeric_limits::epsilon()) { + 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 @@ -47,20 +47,20 @@ F line_source_lfp_factor(const Point3D& e_pos, if (delta <= 0.0) { return f * log_integral(-mu, one_m_mu); } else { - F sqr_delta(std::sqrt(delta)); - F d1(mu - sqr_delta); - F d2(mu + sqr_delta); - F parts = 0.0; + double sqr_delta(std::sqrt(delta)); + double d1(mu - sqr_delta); + double d2(mu + sqr_delta); + double parts = 0.0; if (d1 > 0.0) { - F b(std::min(d1, 1.0) - mu); + double b(std::min(d1, 1.0) - mu); parts += log_integral(-mu, b); } if (d2 < 1.0) { - F b(std::max(d2, 0.0) - mu); + double b(std::max(d2, 0.0) - mu); parts += log_integral(b, one_m_mu); }; // complement - F maxd1_0(std::max(d1, 0.0)), mind2_1(std::min(d2, 1.0)); + 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); } @@ -131,8 +131,8 @@ template LFPCalculator::LFPCalculator(const lfputils::Point3Ds& seg const std::vector& segment_ids, const lfputils::Point3Ds& electrodes, double extra_cellular_conductivity); -template void LFPCalculator::lfp(const Vec& membrane_current); -template void LFPCalculator::lfp(const Vec& membrane_current); +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); diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index a09794221..23030cf76 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -11,16 +11,15 @@ namespace coreneuron { namespace lfputils { -using F = double; -using Point3D = std::array; +using Point3D = std::array; using Point3Ds = std::vector; -using Vec = double*; +using DoublePtr = double*; -inline F dot(const Point3D& p1, const Point3D& p2) { +inline double dot(const Point3D& p1, const Point3D& p2) { return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2]; } -inline F norm(const Point3D& p1) { +inline double norm(const Point3D& p1) { return std::sqrt(dot(p1, p1)); } @@ -28,7 +27,7 @@ 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 F alpha, const Point3D& p2) { +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]}; } @@ -37,13 +36,13 @@ inline Point3D paxpy(const Point3D& p1, const F alpha, const Point3D& p2) { * \param e_pos electrode position * \param seg_pos segment position * \param radius segment radius - * \param f conductivity factor 1/([4 pi] * [conductivity]) + * \param double conductivity factor 1/([4 pi] * [conductivity]) * \return Resistance of the medium from the segment to the electrode. */ -inline F point_source_lfp_factor(const Point3D& e_pos, +inline double point_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_pos, - const F radius, - const F f) { + 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); @@ -57,11 +56,11 @@ inline F point_source_lfp_factor(const Point3D& e_pos, * \param f conductivity factor 1/([4 pi] * [conductivity]) * \return Resistance of the medium from the segment to the electrode. */ -F line_source_lfp_factor(const Point3D& e_pos, +double line_source_lfp_factor(const Point3D& e_pos, const Point3D& seg_0, const Point3D& seg_1, - const F radius, - const F f); + const double radius, + const double f); } // namespace lfputils enum LFPCalculatorType { LineSource, PointSource }; @@ -95,31 +94,31 @@ struct LFPCalculator { std::vector lfp_values; private: - inline lfputils::F getFactor(const lfputils::Point3D& e_pos, + inline double getFactor(const lfputils::Point3D& e_pos, const lfputils::Point3D& seg_0, const lfputils::Point3D& seg_1, - const lfputils::F radius, - const lfputils::F f) const; + const double radius, + const double f) const; std::vector> m; const std::vector& segment_ids_; }; template <> -lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, +double LFPCalculator::getFactor(const lfputils::Point3D& e_pos, const lfputils::Point3D& seg_0, const lfputils::Point3D& seg_1, - const lfputils::F radius, - const lfputils::F f) const { + const double radius, + const double f) const { return lfputils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f); } template <> -lfputils::F LFPCalculator::getFactor(const lfputils::Point3D& e_pos, +double LFPCalculator::getFactor(const lfputils::Point3D& e_pos, const lfputils::Point3D& seg_0, const lfputils::Point3D& seg_1, - const lfputils::F radius, - const lfputils::F f) const { + const double radius, + const double f) const { return lfputils::point_source_lfp_factor(e_pos, lfputils::barycenter(seg_0, seg_1), radius, f); } @@ -135,8 +134,8 @@ extern template LFPCalculator::LFPCalculator(const lfputils::Point3 const std::vector& segment_ids, const lfputils::Point3Ds& electrodes, double extra_cellular_conductivity); -extern template void LFPCalculator::lfp(const lfputils::Vec& membrane_current); -extern template void LFPCalculator::lfp(const lfputils::Vec& membrane_current); +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 From f80818a53320a190bcac4a0e625d33995832d876 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Wed, 24 Feb 2021 15:36:04 +0100 Subject: [PATCH 27/36] Code changes following review. --- coreneuron/io/lfp.cpp | 7 +++---- coreneuron/io/lfp.hpp | 6 ++++-- tests/unit/lfp/lfp.cpp | 4 ++-- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp index 00716fda0..d89307ee6 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -108,14 +108,13 @@ inline void LFPCalculator::lfp(const Vector& membrane_current res[k] += ms[l] * membrane_current[segment_ids_[l]]; } } - std::vector res_reduced(res.size()); #if NRNMPI + lfp_values_.resize(res.size()); int mpi_sum{1}; - nrnmpi_dbl_allreduce_vec(res.data(), res_reduced.data(), res.size(), mpi_sum); + nrnmpi_dbl_allreduce_vec(res.data(), lfp_values_.data(), res.size(), mpi_sum); #else - res_reduced = res; + std::swap(res, lfp_values_); #endif - lfp_values = res_reduced; } diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 23030cf76..aa140812b 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -91,7 +91,9 @@ struct LFPCalculator { template void lfp(const Vector& membrane_current); - std::vector lfp_values; + const std::vector lfp_values() const noexcept { + return lfp_values_; + } private: inline double getFactor(const lfputils::Point3D& e_pos, @@ -99,7 +101,7 @@ struct LFPCalculator { const lfputils::Point3D& seg_1, const double radius, const double f) const; - + std::vector lfp_values_; std::vector> m; const std::vector& segment_ids_; }; diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index d51a64d44..c1180dd0f 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -90,11 +90,11 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { 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; + 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; + 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 From 299f50c869df5126611cc8d4a3ebb327201f5b44 Mon Sep 17 00:00:00 2001 From: Del Marmol Baudouin Date: Wed, 24 Feb 2021 15:36:36 +0100 Subject: [PATCH 28/36] Sort includes. --- tests/unit/alignment/alignment.cpp | 2 +- tests/unit/queueing/test_queueing.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/alignment/alignment.cpp b/tests/unit/alignment/alignment.cpp index 4d5e0acea..f654d27a6 100644 --- a/tests/unit/alignment/alignment.cpp +++ b/tests/unit/alignment/alignment.cpp @@ -13,8 +13,8 @@ #include #include -#include #include +#include #include "coreneuron/utils/memory.h" diff --git a/tests/unit/queueing/test_queueing.cpp b/tests/unit/queueing/test_queueing.cpp index 830deebd0..072c8163c 100644 --- a/tests/unit/queueing/test_queueing.cpp +++ b/tests/unit/queueing/test_queueing.cpp @@ -9,8 +9,8 @@ #define BOOST_TEST_MODULE QueueingTest #define TYPE T::cont -#include #include +#include #include #include From 0db2067bb8ae48eaba0d9fee9e966d41d2fb73fc Mon Sep 17 00:00:00 2001 From: Tristan Carel Date: Wed, 24 Feb 2021 15:44:06 +0100 Subject: [PATCH 29/36] remove semi-colon after closing namespace --- coreneuron/io/lfp.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index aa140812b..f74eab108 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -140,4 +140,4 @@ extern template void LFPCalculator::lfp(const lfputils::DoublePtr& m 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 +} // namespace coreneuron From d747e66bd41ea7a0f23ae70844dc5106c14b6fd4 Mon Sep 17 00:00:00 2001 From: bdelmarm <57898639+bdelmarm@users.noreply.github.com> Date: Wed, 24 Feb 2021 15:45:24 +0100 Subject: [PATCH 30/36] Update coreneuron/io/lfp.hpp Co-authored-by: Tristan Carel --- coreneuron/io/lfp.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index f74eab108..ea6a54dc6 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -91,7 +91,7 @@ struct LFPCalculator { template void lfp(const Vector& membrane_current); - const std::vector lfp_values() const noexcept { + const std::vector& lfp_values() const noexcept { return lfp_values_; } From a2e4e82ff9a94ffa719fff77e604806660acf6b4 Mon Sep 17 00:00:00 2001 From: Tristan Carel Date: Wed, 24 Feb 2021 15:51:56 +0100 Subject: [PATCH 31/36] Fix formatting of lfp header and compilation unit --- coreneuron/io/lfp.cpp | 10 +++++----- coreneuron/io/lfp.hpp | 40 ++++++++++++++++++++-------------------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp index d89307ee6..7677880da 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -10,10 +10,10 @@ namespace coreneuron { 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) { + 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); @@ -111,7 +111,7 @@ inline void LFPCalculator::lfp(const Vector& membrane_current #if NRNMPI lfp_values_.resize(res.size()); int mpi_sum{1}; - nrnmpi_dbl_allreduce_vec(res.data(), lfp_values_.data(), res.size(), mpi_sum); + nrnmpi_dbl_allreduce_vec(res.data(), lfp_values_.data(), res.size(), mpi_sum); #else std::swap(res, lfp_values_); #endif diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index ea6a54dc6..2af81cdec 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -3,9 +3,9 @@ #include #include +#include "coreneuron/mpi/nrnmpidec.h" #include "coreneuron/nrnconf.h" #include "coreneuron/utils/nrn_assert.h" -#include "coreneuron/mpi/nrnmpidec.h" namespace coreneuron { @@ -40,9 +40,9 @@ inline Point3D paxpy(const Point3D& p1, const double alpha, const Point3D& p2) { * \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) { + 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); @@ -57,10 +57,10 @@ inline double point_source_lfp_factor(const Point3D& e_pos, * \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); + const Point3D& seg_0, + const Point3D& seg_1, + const double radius, + const double f); } // namespace lfputils enum LFPCalculatorType { LineSource, PointSource }; @@ -97,10 +97,10 @@ struct LFPCalculator { 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; + 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_; @@ -108,19 +108,19 @@ struct LFPCalculator { 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 { + 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 { + 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); } From e90f8229a6a1c9ab60a7a10c6a4af1061cf6c16e Mon Sep 17 00:00:00 2001 From: Tristan Carel Date: Wed, 24 Feb 2021 16:11:03 +0100 Subject: [PATCH 32/36] fixup! Fix formatting of lfp header and compilation unit --- coreneuron/io/lfp.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 2af81cdec..9cec255dc 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -92,7 +92,7 @@ struct LFPCalculator { void lfp(const Vector& membrane_current); const std::vector& lfp_values() const noexcept { - return lfp_values_; + return lfp_values_; } private: From 19d374e27bbea5be94ce53905fa9d545e438ace6 Mon Sep 17 00:00:00 2001 From: Ioannis Magkanaris Date: Wed, 24 Feb 2021 19:03:04 +0100 Subject: [PATCH 33/36] Fix passing value of pi --- coreneuron/io/lfp.cpp | 2 +- tests/unit/lfp/lfp.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp index 7677880da..d702f3fe2 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -85,7 +85,7 @@ LFPCalculator::LFPCalculator(const Point3Ds& seg_start, 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 * 3.141592653589)); + double f(1.0 / (extra_cellular_conductivity * 4.0 * pi)); m.resize(electrodes.size()); for (size_t k = 0; k < electrodes.size(); ++k) { diff --git a/tests/unit/lfp/lfp.cpp b/tests/unit/lfp/lfp.cpp index c1180dd0f..b65a378a0 100644 --- a/tests/unit/lfp/lfp.cpp +++ b/tests/unit/lfp/lfp.cpp @@ -32,7 +32,7 @@ BOOST_AUTO_TEST_CASE(LFP_PointSource_LineSource) { std::array segment_end = paxpy(segment_start, 1.0, std::array{0.0, 0.0, segment_length}); double floor{1.0e-6}; - double pi = 3.141592653589; + pi = 3.141592653589; std::array vals; double circling_radius{1.0e-6}; From 1983ca84dcae7dc8734b16ab6d33f510b2b4c65f Mon Sep 17 00:00:00 2001 From: Ioannis Magkanaris Date: Thu, 25 Feb 2021 10:58:11 +0100 Subject: [PATCH 34/36] Add pi in pragma update --- coreneuron/apps/main1.cpp | 2 +- coreneuron/io/lfp.cpp | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) 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/lfp.cpp b/coreneuron/io/lfp.cpp index d702f3fe2..18e879416 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -7,6 +7,9 @@ namespace coreneuron { +// extern variables require acc declare +#pragma acc declare create(celsius) + namespace lfputils { double line_source_lfp_factor(const Point3D& e_pos, From 18ca214908127c7845e76f64417d5c65db6d8013 Mon Sep 17 00:00:00 2001 From: Ioannis Magkanaris Date: Thu, 25 Feb 2021 11:02:22 +0100 Subject: [PATCH 35/36] Fix in pragma declare --- coreneuron/io/lfp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/coreneuron/io/lfp.cpp b/coreneuron/io/lfp.cpp index 18e879416..64e7393c5 100644 --- a/coreneuron/io/lfp.cpp +++ b/coreneuron/io/lfp.cpp @@ -8,7 +8,7 @@ namespace coreneuron { // extern variables require acc declare -#pragma acc declare create(celsius) +#pragma acc declare create(pi) namespace lfputils { From cfc7000939f3347a6a8248967550a69e021f5147 Mon Sep 17 00:00:00 2001 From: bdelmarm <57898639+bdelmarm@users.noreply.github.com> Date: Thu, 25 Feb 2021 12:24:21 +0100 Subject: [PATCH 36/36] Update coreneuron/io/lfp.hpp Co-authored-by: Ioannis Magkanaris --- coreneuron/io/lfp.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/coreneuron/io/lfp.hpp b/coreneuron/io/lfp.hpp index 9cec255dc..e47da3c7d 100644 --- a/coreneuron/io/lfp.hpp +++ b/coreneuron/io/lfp.hpp @@ -72,7 +72,6 @@ template struct LFPCalculator { /** * LFP Calculator constructor - * \param comm MPI communicator * \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