[MRG] MAINT: remove move_to_pos #314

jasmainak · 2021-04-08T04:42:49Z

closes #202

jasmainak · 2021-04-08T04:44:53Z

There is a slight mismatch with original "ground truth" dipole if we do this but I would argue this is more correct. Any thoughts @rythorpe @cjayb @ntolley ?

we could update the old "ground truth" measurement now that we've stress-tested HNN-core and we know it can reproduce the old HNN outputs faithfully

cjayb · 2021-04-08T18:49:28Z

It would be super useful to learn something about the reason that was implemented in the first place! It seems utterly bizarre to not treat x and y coords symmetrically. My concern with just removing it is that this seems like yet another tweak that has a specific and thus "intended" purpose, even though it will most likely be wrong for most other simulations with other goals. perhaps something to do with CSD analysis?

jasmainak · 2021-04-08T19:47:42Z

I agree, we should not merge this until we are sure what was the original intention. My impression was that it was meant for visualization. However, this visualization does not exist any more either in the GUI and there is an equivalent viz in hnn-core that does not depend on this.

My impression was that the outputs change only slightly because of numerical precision issues when using h.pt3dchange. This would reflects how floating points are treated differently in Python vs C. The translation operation should not change anything such as _pardistance and dipole computation as it does not change relative distances? I tried doing h.pt3dchange(x, y, z) to test this hypothesis and that failed the tests too. So, not 100% sure what's going on but I'm still betting my money on rounding issues ...

@samnemo and @stephanie-r-jones any opinions?

samnemo · 2021-04-09T14:48:05Z

You're right that those pt3dchange calls were only written to move the neurons vertically for visualization purposes & shouldn't influence the dipole calculations or circuit connectivity.

…

________________________________ IMPORTANT NOTICE: This e-mail is meant only for the use of the intended recipient. It may contain confidential information which is legally privileged or otherwise protected by law. If you received this e-mail in error or from someone who was not authorized to send it to you, you are strictly prohibited from reviewing, using, disseminating, distributing or copying the e-mail. PLEASE NOTIFY US IMMEDIATELY OF THE ERROR BY RETURN E-MAIL AND DELETE THIS MESSAGE FROM YOUR SYSTEM. Thank you for your cooperation.

rythorpe · 2021-04-23T19:25:22Z

I'm fairly certain @jasmainak is right about the discrepancy being caused by rounding error. To demonstrate this, I've printed self.soma.L before and after calling h.pt3dchange() in the default example in master. As network_builder iterates over all cells and moves them into their final positions, I get the following output snippet. Note that L2Pyr cells should have somas that are 22.1 um long and L5Pyr and baskets cells should have somas that are 39.0 um long.

Before/after:

22.1 22.0999755859375
22.1 22.0999755859375
22.1 22.0999755859375
22.1 22.0999755859375
22.1 22.0999755859375
39.0 39.0
39.0 39.0
39.0 39.0
39.0 39.0
39.0 39.0

Now, if we only move the root compartment (i.e., the soma) of each neuron instead of all sections individually, NEURON should move the whole tree without introducing the rounding error caused by converting floats (except for any rounding error introduced into the soma itself).

    def move_to_pos(self):
        """Move cell to position."""
        x0 = self.soma.x3d(0)
        y0 = self.soma.y3d(0)
        z0 = self.soma.z3d(0)
        dx = self.pos[0] * 100 - x0
        dy = self.pos[2] - y0
        dz = self.pos[1] * 100 - z0
        
        root_sec = self.soma
        for i in range(root_sec.n3d()):
            h.pt3dchange(i, root_sec.x3d(i) + dx, root_sec.y3d(i) + dy,
                         root_sec.z3d(i) + dz, root_sec.diam3d(i),
                         sec=root_sec)

As expected, this throws a similar error as this PR, albeit with fewer mismatched elements. I think the difference here is caused by the fact that we still introduced a small amount of rounding error into the somas themselves while avoiding it in all other segments.

E       AssertionError: 
E       Arrays are not equal
E       
E       Mismatched elements: 13 / 6801 (0.191%)
E       Max absolute difference: 0.0001
E       Max relative difference: 0.00011343
E        x: array([-4.0000e-04, -4.0000e-04, -4.0000e-04, ...,  9.0972e+00,
E               9.0651e+00,  9.0332e+00])
E        y: array([-4.0000e-04, -4.0000e-04, -4.0000e-04, ...,  9.0972e+00,
E               9.0651e+00,  9.0332e+00])

tests/test_parallel_backends.py:167: AssertionError

jasmainak · 2021-04-23T23:33:22Z

@rythorpe @cjayb I actually found a neat way to demo it's a numerical precision issue. Check the last commit 77d4b31. Tests pass on this

What do you think?

codecov-commenter · 2021-04-24T00:04:31Z

Codecov Report

Merging #314 (c4fbd81) into master (5c1d2a2) will decrease coverage by 0.04%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #314      +/-   ##
==========================================
- Coverage   89.85%   89.80%   -0.05%     
==========================================
  Files          13       13              
  Lines        2473     2462      -11     
==========================================
- Hits         2222     2211      -11     
  Misses        251      251

Impacted Files	Coverage Δ
hnn_core/network_builder.py	`92.04% <ø> (-0.12%)`	⬇️
hnn_core/cell.py	`97.19% <100.00%> (-0.11%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 5c1d2a2...c4fbd81. Read the comment docs.

rythorpe · 2021-04-26T00:01:34Z

@rythorpe @cjayb I actually found a neat way to demo it's a numerical precision issue. Check the last commit 77d4b31. Tests pass on this

What do you think?

Sure enough, this makes sense to me. My vote is to move forward with removing the move_to_pos() method and update our ground-truth dipole accordingly.

cjayb · 2021-04-26T15:07:29Z

Isn't the correct fix to remove move_to_pos entirely, though?

So the conclusion is that move_to_pos has lead to slightly wrong currents due to rounding errors? That is, that even though the function was intended to only influence visualisation, it actually also changed the output dipole moment (due to rounding errors)?

This will affect the LFP calculation too.

jasmainak · 2021-04-26T15:16:33Z

yep, I'll remove move_to_pos entirely and update the "ground truth". Just wanted to demo it first and be sure that is indeed the issue

cjayb · 2021-04-26T16:02:40Z

Great! This means I'll take a look at removing the legacy-flag in Network! Now that we'll diverge from legacy GUI output, we'll fly free ;)

jasmainak · 2021-05-05T04:03:14Z

so it looks like Neuron 8.0 also produces some slight differences in results. Not sure it would be straightforward or worth it to track this since none of us are Neuron developers and/or have development version installed. I propose that we accept this slight difference and update the ground truth. What do you guys think?

Question is -- do we update the "ground truth" and require NEURON 8.0 from next release of hnn-core? If we still want to support older NEURON version, we should have two different "ground truths" to compare against. Another possibility is that we could have a "tolerance" parameter, but it's possible that some small bugs pass unnoticed since a tiny change in one missing parameter does not always have a big impact on the simulation output.

All these issues point to the problems of comparing the outputs rather than testing the internals using unit tests. Should we plan to have a "test-a-thon" (hackathon where we write tests) and try to phase out or relax the "ground truth" criteria after that?

jasmainak · 2021-05-05T04:18:21Z

BTW, I just remembered that we used to actually have the dev version of Neuron on Travis in the past. These are the steps you have to follow. And then git bisect the heck out of it. Anybody wants to give it a shot?

rythorpe · 2021-05-18T14:55:44Z

I think we should resolve this before implementing the updated calcium version. Would it be worth attempting to position each cell to its final destination upon construction (i.e., rather than creating the cell and then moving it), or is the goal to completely get rid of actually positioning the Neuron cell object?

jasmainak · 2021-05-18T16:57:46Z

How does this work with the LFP calculation? Are the positions in Neuron used for LFP calculation @cjayb ?

cjayb · 2021-05-18T18:15:15Z

Yes, "real" cell positions are needed for LFP: they are used to calculate the distances to the electrode tips.

jasmainak · 2021-05-18T19:52:08Z

could we make cell spacing an argument to Network class so folks can tune it what they think is realistic for LFP simulation? We can throw an error if there are LFP electrodes and the spacing is too large to be realistic.

Instead of changing it through move_to_pos we do it in Python itself here

cjayb · 2021-05-19T10:08:04Z

How about a Cell.set_rotation-method like LFPy? We'd then do a rotation by -np.pi to get them into the XZ-plane.

On the other hand, why not just decide that in hnn_core, the pyramidal cell apical axis defines the z-axis? We'd need to change _secs_LxPyr for L5 and L2 pyramidals, that's it. Am I missing something there?

My point: if there's no good reason for our cells to lie flat in the XY-plane, why not just re-align them and be done with it?

LFP calculation is affected by the distances, but also geometry: our cells are effectively 2D, as there are only 2 basal dendrites, and the 'oblique' apical stump is in-plane with them.

rythorpe · 2021-05-19T17:39:09Z

could we make cell spacing an argument to Network class so folks can tune it what they think is realistic for LFP simulation? We can throw an error if there are LFP electrodes and the spacing is too large to be realistic.

Instead of changing it through move_to_pos we do it in Python itself here

I'm working on a different branch at the moment that refactors how cells, their positions, and their gids are assigned. Let's save implementing a cell spacing argument at the network-level for that PR and just resolve the move_to_pos dilemma.

What do you think of just changing Cell.move_to_pos() to modify Cell.p_secs in the constructor so that there is a 1-to-1 comparison between the section points specified in p_secs (as accessed by the user in Network.cell_types) and what is actually implemented in the model?

rythorpe · 2021-05-19T17:41:03Z

+1 for redefining _secs_LxPyr so that the apical dendrite extends along the z-axis!

jasmainak · 2021-05-19T19:24:23Z

On the other hand, why not just decide that in hnn_core, the pyramidal cell apical axis defines the z-axis? We'd need to change _secs_LxPyr for L5 and L2 pyramidals, that's it. Am I missing something there?

If we load a different neuron at some point in the future, we'd have to make sure that the conversion from Neuron y axis to our z-axis is done ...

jasmainak · 2021-05-19T19:26:29Z

What do you think of just changing Cell.move_to_pos() to modify Cell.p_secs in the constructor so that there is a 1-to-1 comparison between the section points specified in p_secs (as accessed by the user in Network.cell_types) and what is actually implemented in the model?

+1 for this

cjayb · 2021-05-28T13:20:12Z

If we load a different neuron at some point in the future, we'd have to make sure that the conversion from Neuron y axis to our z-axis is done ...

I'm not sure there really is a "Neuron convention"? In any case, a morphology-reader should not assume anything, but be told:

morph_L5 = read_cell_morphology(filename_L5, cell_type='pyramidal', apical_orientation=(0, 1, 0))
net = Network(..., morphology={'L5Pyr': morph_L5, 'L2Pyr': ...})

It'll be a while before this becomes relevant, but the above would be pretty transparent, no?

jasmainak · 2021-05-28T15:56:50Z

I suppose Neuron convention is whatever is specified here: https://www.neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/topology/geometry.html#pt3dadd. If we decide not to use it in the future, we can fly free!

jasmainak · 2021-05-28T16:01:49Z

@rythorpe in terms of plan of action, I think step 1 is to get this PR into mergeable step so we can remove pt3dadd. This would involve updating / adding data to the data branch and making the CIs pass. For step 2, we can make another PR to change the coordinates.

rythorpe · 2021-05-31T02:57:54Z

I think step 1 is to get this PR into mergeable step so we can remove pt3dadd. This would involve updating / adding data to the data branch and making the CIs pass. For step 2, we can make another PR to change the coordinates.

@jasmainak you mean remove h.pt3dchange right? We'll need to still initialize the points to create Neuron morphology somehow, we just will initialize them to their final position upon creation so that we're not moving them around in Neuron-space.

Am I clear to rebase and push?

jasmainak · 2021-05-31T03:00:03Z

yep you can rebase and push. I am not touching this PR.

remember that the positions in space don't affect the dipole calculation. They only affect the LFP calculation ...

rythorpe · 2021-05-31T03:03:56Z

Sure, but the plan is still to maintain consistent positioning between what is specified in Cell.pos and the corresponding Neuron sections right?

jasmainak · 2021-05-31T03:11:40Z

right, I suppose you'd just have to modify this line:

hnn-core/hnn_core/cell.py

Line 259 in 5c1d2a2

h.pt3dadd(pt[0], pt[1], pt[2], 1, sec=sec)

to:

h.pt3dadd(self.pos[0] + pt[0], 
          self.pos[1] + pt[2],  # Cell z-coordinate is matplotlib y-coordinate / or change cell definition
          self.pos[2] + pt[1], 1, sec=sec)

or some such thing

rythorpe · 2021-05-31T04:27:08Z

I was originally thinking of updating self.p_secs, but I think your method above will be cleaner especially for visualization.

What's the best way to update the data? PR to the "test_data" branch?

jasmainak · 2021-05-31T05:24:21Z

Yep, PR to "test_data" branch would be the best way!

jasmainak · 2021-05-31T19:49:17Z

can you rebase this PR @rythorpe ?

jasmainak · 2021-05-31T19:50:53Z

anything else to do for the PR to be MRG?

rythorpe · 2021-05-31T19:56:08Z

I rebased last night, have we merged to master since then?

jasmainak · 2021-05-31T19:57:27Z

oh indeed, I restarted the CI to redownload the new test data

rythorpe · 2021-05-31T19:58:53Z

Once CI's turn green, I think we're good to go.

rythorpe · 2021-05-31T20:10:36Z

I'm going to go ahead and merge since tests are passing and this PR is something we've all weighed in on.

jasmainak · 2021-05-31T20:11:57Z

Perfect, thanks @rythorpe for taking it to the finish line!

rythorpe · 2021-05-31T20:16:13Z

You did almost all of the real work here @jasmainak! Stayed tuned for a follow up PR re: cell spacing.

cjayb · 2021-06-01T07:31:18Z

hnn_core/cell.py

+        dx = self.pos[0] - self.p_secs['soma']['sec_pts'][0][0]
+        dy = self.pos[2] - self.p_secs['soma']['sec_pts'][0][1]
+        dz = self.pos[1] - self.p_secs['soma']['sec_pts'][0][2]
+
        for sec_name in p_secs:
            sec = h.Section(name=f'{self.name}_{sec_name}')
            self.sections[sec_name] = sec

            h.pt3dclear(sec=sec)
            for pt in p_secs[sec_name]['sec_pts']:
-                h.pt3dadd(pt[0], pt[1], pt[2], 1, sec=sec)
+                h.pt3dadd(pt[0] + dx,
+                          pt[1] + dy,
+                          pt[2] + dz, 1, sec=sec)


The consequences of these lines are very confusing to me (still: I realise nothing has changed in that respect in this PR). Are we reorienting the apical sections to align with the Z-axis here?

In any case, we now have a situation where net.pos_dict, Cell.pos and the Section end points are inconsistent, right?

@rythorpe are you planning to address this in the follow-up you mention? I'm working on some qualified opinions about the spacing-issue too, based on "reduced" model literature.

yes, we should be reorienting the apical sections to align with the Z-axis ... but maybe it's not the case?

The Section end points are specified relative to the Cell.pos or net.pos_dict ...

now that I think of it, to reorient along the Z-axis, you'd need to do:

dy = self.pos[1] - self.p_secs['soma']['sec_pts'][0][2]

what might be more clear though is to apply a trans matrix first:

trans = [[1, 0, 0], [0, 0, 1], [0, 1, 0]] self.p_secs['soma']['sec_pts'] = np.dot(trans, self.p_secs['soma']['sec_pts'])

there is a wealth of functions we can copy from MNE :) There is even a function to align along a target axis:

https://github.com/mne-tools/mne-python/blob/b3f4b05b3933a3b73458205c7e786166f90ad1f7/mne/transforms.py#L298

that we can use to our advantage

I'm not sure I understand your question @cjayb. Are you asking about the implementation in this PR or what we plan to do in later PRs?

If we change the section endpoints in params_default.py, I think we should do that in a separate PR.

jasmainak changed the title ~~MAINT: remove move_to_pos~~ [WIP] MAINT: remove move_to_pos Apr 8, 2021

jasmainak force-pushed the move_to_pos branch from 061617c to f9bf01a Compare April 8, 2021 04:54

jasmainak force-pushed the move_to_pos branch 2 times, most recently from 3fd8db8 to 77d4b31 Compare April 23, 2021 23:32

jasmainak mentioned this pull request May 5, 2021

[MRG] FIX tests for NEURON v8 #324

Merged

This was referenced May 18, 2021

[MAINT][WIP] Simulate dipole enh #315

Closed

[MRG] add method for updating cell biophysics from Network instance #321

Merged

cjayb mentioned this pull request May 19, 2021

[MRG] Improved LFP calculation #329

Merged

12 tasks

jasmainak added 4 commits May 30, 2021 23:05

MAINT: demo what move_to_pos does

f39fa5d

DOC: update whats new

5a1a432

DEMO its numerical precision

5372338

ENH: dummy comment to demo issue with NEURON 8.0

76c1ab2

set position of sections upon creation

c418697

rythorpe force-pushed the move_to_pos branch from 9eedea6 to c418697 Compare May 31, 2021 04:23

update comment and whats_new.rst

c4fbd81

rythorpe mentioned this pull request May 31, 2021

update default dipole data #344

Merged

jasmainak changed the title ~~[WIP] MAINT: remove move_to_pos~~ [MRG] MAINT: remove move_to_pos May 31, 2021

rythorpe merged commit 0f96ef6 into jonescompneurolab:master May 31, 2021

cjayb reviewed Jun 1, 2021

View reviewed changes

jasmainak mentioned this pull request Jun 10, 2021

[MRG] fix layers 5 and 2 swapped depth dimension #352

Merged

rythorpe mentioned this pull request Jul 19, 2021

[MRG] attach drives with a user-defined number of source (artificial) cells #383

Merged

[MRG] MAINT: remove move_to_pos #314

[MRG] MAINT: remove move_to_pos #314

Conversation

jasmainak commented Apr 8, 2021 • edited Loading

jasmainak commented Apr 8, 2021 • edited Loading

cjayb commented Apr 8, 2021

jasmainak commented Apr 8, 2021

samnemo commented Apr 9, 2021 via email

rythorpe commented Apr 23, 2021

jasmainak commented Apr 23, 2021 • edited Loading

codecov-commenter commented Apr 24, 2021 • edited Loading

Codecov Report

rythorpe commented Apr 26, 2021

cjayb commented Apr 26, 2021

jasmainak commented Apr 26, 2021

cjayb commented Apr 26, 2021

jasmainak commented May 5, 2021

jasmainak commented May 5, 2021

rythorpe commented May 18, 2021

jasmainak commented May 18, 2021

cjayb commented May 18, 2021

jasmainak commented May 18, 2021

cjayb commented May 19, 2021

rythorpe commented May 19, 2021 • edited Loading

rythorpe commented May 19, 2021

jasmainak commented May 19, 2021

jasmainak commented May 19, 2021

cjayb commented May 28, 2021

jasmainak commented May 28, 2021

jasmainak commented May 28, 2021

rythorpe commented May 31, 2021

jasmainak commented May 31, 2021

rythorpe commented May 31, 2021

jasmainak commented May 31, 2021 • edited Loading

rythorpe commented May 31, 2021

jasmainak commented May 31, 2021

jasmainak commented May 31, 2021

jasmainak commented May 31, 2021

rythorpe commented May 31, 2021

jasmainak commented May 31, 2021

rythorpe commented May 31, 2021

rythorpe commented May 31, 2021

jasmainak commented May 31, 2021

rythorpe commented May 31, 2021

cjayb Jun 1, 2021

Choose a reason for hiding this comment

jasmainak Jun 1, 2021 • edited Loading

Choose a reason for hiding this comment

jasmainak Jun 1, 2021

Choose a reason for hiding this comment

rythorpe Jun 2, 2021

Choose a reason for hiding this comment

jasmainak commented Apr 8, 2021 •

edited

Loading

jasmainak commented Apr 8, 2021 •

edited

Loading

jasmainak commented Apr 23, 2021 •

edited

Loading

codecov-commenter commented Apr 24, 2021 •

edited

Loading

rythorpe commented May 19, 2021 •

edited

Loading

jasmainak commented May 31, 2021 •

edited

Loading

jasmainak Jun 1, 2021 •

edited

Loading