[MRG] Improved LFP calculation #329
Conversation
Codecov Report
@@ Coverage Diff @@
## master #329 +/- ##
==========================================
- Coverage 91.00% 89.77% -1.23%
==========================================
Files 13 15 +2
Lines 2557 2886 +329
==========================================
+ Hits 2327 2591 +264
- Misses 230 295 +65
Continue to review full report at Codecov.
|
|
I created this gist to validate the LFP calculations are correct. |
This comment has been minimized.
This comment has been minimized.
|
Looks tantalizing! By the way, your script could use some print statements with |
Really?! Can you run it cell-by-cell? All the laminar analysis should work for sure, I suspect the problem is with trying to reset NEURON ( I've separated the laminar and grid analyses into separate gists (links here and here, respectively). The 20 x 20 grid is super-slow, indicating some work is needed to optimise the calculation. I couldn't get |
This comment has been minimized.
This comment has been minimized.
This comment has been minimized.
This comment has been minimized.
hnn_core/network.py
Outdated
| @@ -849,6 +860,11 @@ def _reset_drives(self): | |||
| for drive_name in self.external_drives.keys(): | |||
| self.external_drives[drive_name]['events'] = list() | |||
|
|
|||
| def _reset_rec_array(self): | |||
There was a problem hiding this comment.
I have thought about these reset methods. I think this is an indication that these objects don't belong as attributes of the Network object. They should be instantiated in simulate_dipole and returned separately
hnn_core/tests/test_extracellular.py
Outdated
| net = default_network(deepcopy(params), add_drives_from_params=True) | ||
|
|
||
| # Test LFP electrodes | ||
| electrode_pos = (2, 400, 2) |
There was a problem hiding this comment.
is this still with the xy axis swap? Probably doesn't matter for test but might be worth switching
There was a problem hiding this comment.
throw a warning/error if all electrodes are too far from any cell?
There was a problem hiding this comment.
Didn't see that! That test function just tests the API, but I still made some clarifications to emphasise what's being tested.
Not sure there's a heuristic for "too far", especially since we want to compare with "macroscopic" recordings.
hnn_core/tests/test_extracellular.py
Outdated
|
|
||
| for method in ['psa', 'lsa']: | ||
| res = _transfer_resistance(sec, elec_pos, conductivity, method) | ||
| assert_allclose(res, target_vals[method]) |
There was a problem hiding this comment.
why allclose and not an exact equal?
There was a problem hiding this comment.
Good question. I guess I just avoid "exactness" when it comes to floating point calculations. Note that the test is an independent implementation of the maths in the calculating function.
I've reduced rtol to 1e-12, doesn't get much more 'exact' than that :)
hnn_core/tests/test_extracellular.py
Outdated
| return (Q * cosT) / (4 * np.pi * R ** 2) | ||
|
|
||
|
|
||
| # require MPI to speed up due to large number of extracellular electrodes |
There was a problem hiding this comment.
but this will not speed up in the CIs ...
There was a problem hiding this comment.
while I like this code, I think this is not really a "unit" test. Could we perhaps maintain this in your example repository?
There was a problem hiding this comment.
I'm not sure I know the definition, but this is clearly the slowest test. I guess I'm OK with cutting it out, will push update soon.
| 't_evprox_1': 7, | ||
| 't_evdist_1': 17}) |
There was a problem hiding this comment.
why not using your own API? :)
There was a problem hiding this comment.
for brevity! I hope to address this in a later proposal. I think the API should be something like:
net = default_network(params)
net.add_drives_from_params(modify_params={'t_evprox_1': 7, 't_evdist_1': 17})WDYT?
There was a problem hiding this comment.
do you really want to be loading drives from files and modifying them on the fly? It feels very opaque ...
There was a problem hiding this comment.
A discussion for another thread ;)
|
Okay I did a full read through of the code ... looks reasonable to me. Let's iterate more after merging. I don't want to block other PRs for fear of rebase issues with this PR. Would you be able to look at the math/code soon @ntolley ? |
|
Honestly I am ok with merging now and then raising an issue in the unlikely case that the math doesn't add up. It's more or less just my own desire to understand the calculations. Everything else I've been able to read carefully and am quite happy with what I've seen! Also I'm less concerned since this is a completely new feature and doesn't risk breaking existing functionality. |
| # To avoid very large values when electrode is placed (anywhere) on | ||
| # the section axis, enforce minimal perpendicular distance | ||
| R2 = np.maximum(R2, min_distance ** 2) |
There was a problem hiding this comment.
Maybe just think about what min_distance is doing before merging @jasmainak and @ntolley
I think this'll come up again when repositioning the cells, so I wouldn't call it a show-stopper for merging this PR. Just something to keep in mind.
There was a problem hiding this comment.
humm ... why change R without changing b etc? Won't the math screw up then?
There was a problem hiding this comment.
Well the math gets "screwed up" either way, but this is to avoid division by zero. The situation is clearer in PSA, where the euclidean distance is restricted to be > min_distance. I noticed that LFPy and NetPyNE use a limit on R, too. We can return to this: it will resurface when talking about physical distances. Right now it "works".
There was a problem hiding this comment.
If the electrode is so close to section axis, wouldn't the voltage simply be the membrane potential? I guess accessible with seg.v ?
| each segment junction as a point extracellular current source. | ||
| ``'lsa'`` (line source approximation) treats each segment as a line | ||
| source of current, which extends from the previous to the next segment | ||
| center point: |---x---|, where x is the current segment flanked by |. |
There was a problem hiding this comment.
what is | in this diagram?
There was a problem hiding this comment.
The 'x' is an element of seg_ctr, and the | are start and end on L118-9.
|
@cjayb are you ready to merge? I think I won't be able to sleep at night with a 1500 line PR waiting to be merged and 300 comments ... :) |
|
I'm ready to let this go :) There'll be plenty of opportunities to scrutinise the details in conjunction with the 0.2 release. Merge away, I say! I'll take the heat if we find something idiotic in there ;) I'm (also) very convinced that eventually this needs to be pulled out of |
|
Indeed, I had plenty to work with from the get-go! |
This PR builds on and supersedes #150 (commits by @ntolley and @jasmainak retained in squash)
sec(0.5); see this gist for details)and CSDcalculation in (weak) gamma examplewhats_new.rst&api.rstcreate demo using single instance of(move to explore response properties of individual HNN network cells #323)Celladd tests for validity of PSA and LSA values close to cells (how?)not needed as long as the above three tests passoptimise for speedfuture PR should improve implementation, possibly getting rid ofextra_gather_scatterCloses #68