Skip to content
This repository was archived by the owner on Mar 20, 2023. It is now read-only.

Commit b234502

Browse files
bdelmarmDel Marmol Baudouinpramodkiomaganaristristan0x
authored
LFP Calculator (#492)
* New LFP calculator class - Allows for LineSource and PointSource approximations * Read value of `pi` from NEURON in online or offline mode Co-authored-by: Del Marmol Baudouin <[email protected]> Co-authored-by: Pramod S Kumbhar <[email protected]> Co-authored-by: Ioannis Magkanaris <[email protected]> Co-authored-by: Tristan Carel <[email protected]> Co-authored-by: Tristan Carel <[email protected]>
1 parent 0219aa6 commit b234502

File tree

11 files changed

+420
-6
lines changed

11 files changed

+420
-6
lines changed

‎coreneuron/apps/main1.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -490,7 +490,7 @@ extern "C" int run_solve_core(int argc, char** argv) {
490490

491491
// clang-format off
492492

493-
#pragma acc update device(celsius, secondorder) if (compute_gpu)
493+
#pragma acc update device(celsius, secondorder, pi) if (compute_gpu)
494494
// clang-format on
495495
{
496496
double v = corenrn_param.voltage;

‎coreneuron/io/global_vars.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,7 @@ void set_globals(const char* path, bool cli_global_seed, int cli_global_seed_val
4848
(*n2v)["celsius"] = PSD(0, &celsius);
4949
(*n2v)["dt"] = PSD(0, &dt);
5050
(*n2v)["t"] = PSD(0, &t);
51+
(*n2v)["PI"] = PSD(0, &pi);
5152

5253
if (corenrn_embedded) { // CoreNEURON embedded, get info direct from NEURON
5354

‎coreneuron/io/lfp.cpp‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
#include "coreneuron/io/lfp.hpp"
2+
3+
#include <cmath>
4+
#include <limits>
5+
#include <sstream>
6+
7+
8+
namespace coreneuron {
9+
10+
// extern variables require acc declare
11+
#pragma acc declare create(pi)
12+
13+
namespace lfputils {
14+
15+
double line_source_lfp_factor(const Point3D& e_pos,
16+
const Point3D& seg_0,
17+
const Point3D& seg_1,
18+
const double radius,
19+
const double f) {
20+
nrn_assert(radius >= 0.0);
21+
Point3D dx = paxpy(seg_1, -1.0, seg_0);
22+
Point3D de = paxpy(e_pos, -1.0, seg_0);
23+
double dx2(dot(dx, dx));
24+
double dxn(std::sqrt(dx2));
25+
if (dxn < std::numeric_limits<double>::epsilon()) {
26+
return point_source_lfp_factor(e_pos, seg_0, radius, f);
27+
}
28+
double de2(dot(de, de));
29+
double mu(dot(dx, de) / dx2);
30+
Point3D de_star(paxpy(de, -mu, dx));
31+
double de_star2(dot(de_star, de_star));
32+
double q2(de_star2 / dx2);
33+
34+
double delta(mu * mu - (de2 - radius * radius) / dx2);
35+
double one_m_mu(1.0 - mu);
36+
auto log_integral = [&q2, &dxn](double a, double b) {
37+
if (q2 < std::numeric_limits<double>::epsilon()) {
38+
if (a * b <= 0) {
39+
std::ostringstream s;
40+
s << "Log integral: invalid arguments " << b << " " << a
41+
<< ". Likely electrode exactly on the segment and "
42+
<< "no flooring is present.";
43+
throw std::invalid_argument(s.str());
44+
}
45+
return std::abs(std::log(b / a)) / dxn;
46+
} else {
47+
return std::log((b + std::sqrt(b * b + q2)) / (a + std::sqrt(a * a + q2))) / dxn;
48+
}
49+
};
50+
if (delta <= 0.0) {
51+
return f * log_integral(-mu, one_m_mu);
52+
} else {
53+
double sqr_delta(std::sqrt(delta));
54+
double d1(mu - sqr_delta);
55+
double d2(mu + sqr_delta);
56+
double parts = 0.0;
57+
if (d1 > 0.0) {
58+
double b(std::min(d1, 1.0) - mu);
59+
parts += log_integral(-mu, b);
60+
}
61+
if (d2 < 1.0) {
62+
double b(std::max(d2, 0.0) - mu);
63+
parts += log_integral(b, one_m_mu);
64+
};
65+
// complement
66+
double maxd1_0(std::max(d1, 0.0)), mind2_1(std::min(d2, 1.0));
67+
if (maxd1_0 < mind2_1) {
68+
parts += 1.0 / radius * (mind2_1 - maxd1_0);
69+
}
70+
return f * parts;
71+
};
72+
}
73+
} // namespace lfputils
74+
75+
using namespace lfputils;
76+
77+
template <LFPCalculatorType Type, typename SegmentIdTy>
78+
LFPCalculator<Type, SegmentIdTy>::LFPCalculator(const Point3Ds& seg_start,
79+
const Point3Ds& seg_end,
80+
const std::vector<double>& radius,
81+
const std::vector<SegmentIdTy>& segment_ids,
82+
const Point3Ds& electrodes,
83+
double extra_cellular_conductivity)
84+
: segment_ids_(segment_ids) {
85+
if (seg_start.size() != seg_end.size()) {
86+
throw std::invalid_argument("Different number of segment starts and ends.");
87+
}
88+
if (seg_start.size() != radius.size()) {
89+
throw std::invalid_argument("Different number of segments and radii.");
90+
}
91+
double f(1.0 / (extra_cellular_conductivity * 4.0 * pi));
92+
93+
m.resize(electrodes.size());
94+
for (size_t k = 0; k < electrodes.size(); ++k) {
95+
auto& ms = m[k];
96+
ms.resize(seg_start.size());
97+
for (size_t l = 0; l < seg_start.size(); l++) {
98+
ms[l] = getFactor(electrodes[k], seg_start[l], seg_end[l], radius[l], f);
99+
}
100+
}
101+
}
102+
103+
template <LFPCalculatorType Type, typename SegmentIdTy>
104+
template <typename Vector>
105+
inline void LFPCalculator<Type, SegmentIdTy>::lfp(const Vector& membrane_current) {
106+
std::vector<double> res(m.size());
107+
for (size_t k = 0; k < m.size(); ++k) {
108+
res[k] = 0.0;
109+
auto& ms = m[k];
110+
for (size_t l = 0; l < ms.size(); l++) {
111+
res[k] += ms[l] * membrane_current[segment_ids_[l]];
112+
}
113+
}
114+
#if NRNMPI
115+
lfp_values_.resize(res.size());
116+
int mpi_sum{1};
117+
nrnmpi_dbl_allreduce_vec(res.data(), lfp_values_.data(), res.size(), mpi_sum);
118+
#else
119+
std::swap(res, lfp_values_);
120+
#endif
121+
}
122+
123+
124+
template LFPCalculator<LineSource>::LFPCalculator(const lfputils::Point3Ds& seg_start,
125+
const lfputils::Point3Ds& seg_end,
126+
const std::vector<double>& radius,
127+
const std::vector<int>& segment_ids,
128+
const lfputils::Point3Ds& electrodes,
129+
double extra_cellular_conductivity);
130+
template LFPCalculator<PointSource>::LFPCalculator(const lfputils::Point3Ds& seg_start,
131+
const lfputils::Point3Ds& seg_end,
132+
const std::vector<double>& radius,
133+
const std::vector<int>& segment_ids,
134+
const lfputils::Point3Ds& electrodes,
135+
double extra_cellular_conductivity);
136+
template void LFPCalculator<LineSource>::lfp(const DoublePtr& membrane_current);
137+
template void LFPCalculator<PointSource>::lfp(const DoublePtr& membrane_current);
138+
template void LFPCalculator<LineSource>::lfp(const std::vector<double>& membrane_current);
139+
template void LFPCalculator<PointSource>::lfp(const std::vector<double>& membrane_current);
140+
141+
} // namespace coreneuron

‎coreneuron/io/lfp.hpp‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
#pragma once
2+
3+
#include <array>
4+
#include <vector>
5+
6+
#include "coreneuron/mpi/nrnmpidec.h"
7+
#include "coreneuron/nrnconf.h"
8+
#include "coreneuron/utils/nrn_assert.h"
9+
10+
namespace coreneuron {
11+
12+
namespace lfputils {
13+
14+
using Point3D = std::array<double, 3>;
15+
using Point3Ds = std::vector<Point3D>;
16+
using DoublePtr = double*;
17+
18+
inline double dot(const Point3D& p1, const Point3D& p2) {
19+
return p1[0] * p2[0] + p1[1] * p2[1] + p1[2] * p2[2];
20+
}
21+
22+
inline double norm(const Point3D& p1) {
23+
return std::sqrt(dot(p1, p1));
24+
}
25+
26+
inline Point3D barycenter(const Point3D& p1, const Point3D& p2) {
27+
return {0.5 * (p1[0] + p2[0]), 0.5 * (p1[1] + p2[1]), 0.5 * (p1[2] + p2[2])};
28+
}
29+
30+
inline Point3D paxpy(const Point3D& p1, const double alpha, const Point3D& p2) {
31+
return {p1[0] + alpha * p2[0], p1[1] + alpha * p2[1], p1[2] + alpha * p2[2]};
32+
}
33+
34+
/**
35+
*
36+
* \param e_pos electrode position
37+
* \param seg_pos segment position
38+
* \param radius segment radius
39+
* \param double conductivity factor 1/([4 pi] * [conductivity])
40+
* \return Resistance of the medium from the segment to the electrode.
41+
*/
42+
inline double point_source_lfp_factor(const Point3D& e_pos,
43+
const Point3D& seg_pos,
44+
const double radius,
45+
const double f) {
46+
nrn_assert(radius >= 0.0);
47+
Point3D es = paxpy(e_pos, -1.0, seg_pos);
48+
return f / std::max(norm(es), radius);
49+
}
50+
51+
/**
52+
*
53+
* \param e_pos electrode position
54+
* \param seg_pos segment position
55+
* \param radius segment radius
56+
* \param f conductivity factor 1/([4 pi] * [conductivity])
57+
* \return Resistance of the medium from the segment to the electrode.
58+
*/
59+
double line_source_lfp_factor(const Point3D& e_pos,
60+
const Point3D& seg_0,
61+
const Point3D& seg_1,
62+
const double radius,
63+
const double f);
64+
} // namespace lfputils
65+
66+
enum LFPCalculatorType { LineSource, PointSource };
67+
68+
/**
69+
* \brief LFPCalculator allows calculation of LFP given membrane currents.
70+
*/
71+
template <LFPCalculatorType Ty, typename SegmentIdTy = int>
72+
struct LFPCalculator {
73+
/**
74+
* LFP Calculator constructor
75+
* \param seg_start all segments start owned by the proc
76+
* \param seg_end all segments end owned by the proc
77+
* \param radius fence around the segment. Ensures electrode cannot be
78+
* arbitrarily close to the segment
79+
* \param electrodes positions of the electrodes
80+
* \param extra_cellular_conductivity conductivity of the extra-cellular
81+
* medium
82+
*/
83+
LFPCalculator(const lfputils::Point3Ds& seg_start,
84+
const lfputils::Point3Ds& seg_end,
85+
const std::vector<double>& radius,
86+
const std::vector<SegmentIdTy>& segment_ids,
87+
const lfputils::Point3Ds& electrodes,
88+
double extra_cellular_conductivity);
89+
90+
template <typename Vector>
91+
void lfp(const Vector& membrane_current);
92+
93+
const std::vector<double>& lfp_values() const noexcept {
94+
return lfp_values_;
95+
}
96+
97+
private:
98+
inline double getFactor(const lfputils::Point3D& e_pos,
99+
const lfputils::Point3D& seg_0,
100+
const lfputils::Point3D& seg_1,
101+
const double radius,
102+
const double f) const;
103+
std::vector<double> lfp_values_;
104+
std::vector<std::vector<double>> m;
105+
const std::vector<SegmentIdTy>& segment_ids_;
106+
};
107+
108+
template <>
109+
double LFPCalculator<LineSource>::getFactor(const lfputils::Point3D& e_pos,
110+
const lfputils::Point3D& seg_0,
111+
const lfputils::Point3D& seg_1,
112+
const double radius,
113+
const double f) const {
114+
return lfputils::line_source_lfp_factor(e_pos, seg_0, seg_1, radius, f);
115+
}
116+
117+
template <>
118+
double LFPCalculator<PointSource>::getFactor(const lfputils::Point3D& e_pos,
119+
const lfputils::Point3D& seg_0,
120+
const lfputils::Point3D& seg_1,
121+
const double radius,
122+
const double f) const {
123+
return lfputils::point_source_lfp_factor(e_pos, lfputils::barycenter(seg_0, seg_1), radius, f);
124+
}
125+
126+
extern template LFPCalculator<LineSource>::LFPCalculator(const lfputils::Point3Ds& seg_start,
127+
const lfputils::Point3Ds& seg_end,
128+
const std::vector<double>& radius,
129+
const std::vector<int>& segment_ids,
130+
const lfputils::Point3Ds& electrodes,
131+
double extra_cellular_conductivity);
132+
extern template LFPCalculator<PointSource>::LFPCalculator(const lfputils::Point3Ds& seg_start,
133+
const lfputils::Point3Ds& seg_end,
134+
const std::vector<double>& radius,
135+
const std::vector<int>& segment_ids,
136+
const lfputils::Point3Ds& electrodes,
137+
double extra_cellular_conductivity);
138+
extern template void LFPCalculator<LineSource>::lfp(const lfputils::DoublePtr& membrane_current);
139+
extern template void LFPCalculator<PointSource>::lfp(const lfputils::DoublePtr& membrane_current);
140+
extern template void LFPCalculator<LineSource>::lfp(const std::vector<double>& membrane_current);
141+
extern template void LFPCalculator<PointSource>::lfp(const std::vector<double>& membrane_current);
142+
} // namespace coreneuron

‎coreneuron/mechanism/register_mech.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,10 +19,11 @@
1919

2020
namespace coreneuron {
2121
int secondorder = 0;
22-
double t, dt, celsius;
22+
double t, dt, celsius, pi;
2323
// declare copyin required for correct initialization
2424
#pragma acc declare copyin(secondorder)
2525
#pragma acc declare copyin(celsius)
26+
#pragma acc declare copyin(pi)
2627
int rev_dt;
2728

2829
using Pfrv = void (*)();

‎coreneuron/nrnconf.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@ using Symbol = char;
3232
#define VECTORIZE 1
3333

3434
extern double celsius;
35+
extern double pi;
3536

3637
extern double t, dt;
3738
extern int rev_dt;

‎tests/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ if(Boost_FOUND)
2929
add_subdirectory(unit/interleave_info)
3030
add_subdirectory(unit/alignment)
3131
add_subdirectory(unit/queueing)
32+
add_subdirectory(unit/lfp)
3233
endif()
3334
message(STATUS "Boost found, unit tests enabled")
3435
else()

‎tests/unit/alignment/alignment.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,8 @@
1313
#include <cstring>
1414
#include <stdint.h>
1515

16-
#include <boost/test/unit_test.hpp>
17-
#include <boost/test/test_case_template.hpp>
1816
#include <boost/mpl/list.hpp>
17+
#include <boost/test/unit_test.hpp>
1918

2019
#include "coreneuron/utils/memory.h"
2120

‎tests/unit/lfp/CMakeLists.txt‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
# =============================================================================
2+
# Copyright (C) 2016-2020 Blue Brain Project
3+
#
4+
# See top-level LICENSE file for details.
5+
# =============================================================================
6+
7+
include_directories(${CMAKE_SOURCE_DIR}/coreneuron ${Boost_INCLUDE_DIRS})
8+
file(GLOB lfp_test_src "*.cpp")
9+
10+
add_executable(lfp_test_bin ${lfp_test_src})
11+
target_link_libraries(
12+
lfp_test_bin
13+
${MPI_CXX_LIBRARIES}
14+
${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
15+
coreneuron
16+
${corenrn_mech_lib}
17+
${reportinglib_LIBRARY}
18+
${sonatareport_LIBRARY})
19+
20+
if(CORENRN_ENABLE_GPU)
21+
target_link_libraries(lfp_test_bin ${link_cudacoreneuron} ${CUDA_LIBRARIES})
22+
endif()
23+
24+
add_dependencies(lfp_test_bin nrniv-core)
25+
26+
add_test(NAME lfp_test COMMAND ${TEST_EXEC_PREFIX} ${CMAKE_CURRENT_BINARY_DIR}/lfp_test_bin)

0 commit comments

Comments
 (0)