Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MRG: BEM model and solution in Python #2193

Merged
merged 14 commits into from
Jul 13, 2015
Merged

MRG: BEM model and solution in Python #2193

merged 14 commits into from
Jul 13, 2015

Conversation

larsoner
Copy link
Member

@larsoner larsoner commented Jun 8, 2015

Closes #2190.

I translated all the code but haven't tested it at all. In fact I know it will not work. But the current plan for usage is something like this:

from mne import make_bem_model, make_bem_solution
subject = 'sample'
model = make_bem_model(subject)  # puts surfaces and conductances together
solution = make_bem_solution(model)  # does the potential computations

triple = np.dot(cross, y3)

# XXX refactor with _get_solids?
l1 = np.sqrt((y1 * y1).sum())
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(y1 * y1).sum() -> np.dot(y1, y1)

fix below too.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also you use use functions from math module for floats: sqrt, abs ...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah, the reason I have it in this form is that I am almost certainly going to vectorize the operation, in which case having it in the sum form will be more convenient to start translating to higher dimensions

@larsoner
Copy link
Member Author

larsoner commented Jun 9, 2015

It seems to be working for me now, but it is horribly slow. Hopefully I can speed it up tomorrow.

@dgwakeman
Copy link
Contributor

You are awesome! I am so excited for this!

@larsoner
Copy link
Member Author

FYI I added a couple more files to the testing dataset. BEMs have now been created for ico-2 surfaces (320 instead of 5120 or 1280) which is much faster.

The code at this point seems to work for homog 1-layer, I just need to track down a bug in the 3-layer case. Speed isn't too much worse than C for small shells (<= 1280) and actually faster for standard size (5120) because, while calculating the matrix coefficients is slower in Python, the inversion is much faster:

ico       type    layers  time
2 (320)   C       1       0.0300931930542
2 (320)   Python  1       0.250679969788
2 (320)   C       3       0.347531795502
2 (320)   Python  3       2.41311502457
3 (1280)  C       1       0.621548891068
3 (1280)  Python  1       2.16235113144
3 (1280)  C       3       12.2742669582
3 (1280)  Python  3       21.4681899548
4 (5120)  C       1       24.4842910767
4 (5120)  Python  1       27.2853899002
4 (5120)  C       3       643.578150034
4 (5120)  Python  3       282.933717012

Just need to track down where that bug is, which I probably won't get to until next week.

@larsoner larsoner modified the milestone: 0.10 Jun 16, 2015
@larsoner larsoner changed the title WIP: BEM model and solution in Python MRG: BEM model and solution in Python Jul 7, 2015
@larsoner
Copy link
Member Author

larsoner commented Jul 7, 2015

Okay, 1- and 3-layer model and solution computation works now. Ready for review/merge from my end.

@larsoner
Copy link
Member Author

larsoner commented Jul 9, 2015

Ping @agramfort @dengemann this is ready for review

@dengemann
Copy link
Member

@Eric89GXL amazing! How fast is it compared to the original code?

@larsoner
Copy link
Member Author

larsoner commented Jul 9, 2015

See comment above, faster for normal use case (3-layer 5120)

@larsoner
Copy link
Member Author

larsoner commented Jul 9, 2015

There might be additional ways of optimizing it, too, if there is interest

@dengemann
Copy link
Member

here might be additional ways of optimizing it, too, if there is interest

I'm already very impressed. Eric what do we need here, any testing desired?

@larsoner
Copy link
Member Author

larsoner commented Jul 9, 2015

Hopefully the unit tests are sufficient to make us confident in using it. But sure, if you want to test it out and make comments, it's welcome.

@dengemann
Copy link
Member

But sure, if you want to test it out and make comments, it's welcome.

I will trust you if you're confident ;)

#


def _calc_beta(rk, rk1):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

short comment

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That would require (re-)understanding what it does first :) but seriously, I'll add one. It also reminds me that I need to check for duplication of functionality across the _compute_forward module, I think these are different but I should check.

@larsoner
Copy link
Member Author

larsoner commented Jul 9, 2015

Check out the tests, hopefully those will make us confident :) If they don't, then they can be expanded.

* This is the general multilayer case

"""
from scipy import linalg
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why import here?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we nest certain scipy imports now to save on import speed, see mne/tests/test_import_nesting.py

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

right! thanks

@larsoner
Copy link
Member Author

larsoner commented Jul 9, 2015

@dengemann comments addressed. FYI you asked about speed, and that made me optimize a bit more:

2 (320)  1       C       0.283627033234
2 (320)  1       Python  0.291443109512
2 (320)  3       C       0.465950012207
2 (320)  3       Python  2.45959687233
3 (1280)  1       C       0.660896062851
3 (1280)  1       Python  1.62492084503
3 (1280)  3       C       12.4721641541
3 (1280)  3       Python  11.3650319576
4 (5120)  1       C       24.9352688789
4 (5120)  1       Python  11.2639200687
4 (5120)  3       C       644.138396025
4 (5120)  3       Python  118.89906311

@dengemann
Copy link
Member

Great!!! I will recompute all my BEMs to enjoy the speedy drive :)

On 09 Jul 2015, at 18:31, Eric Larson notifications@github.com wrote:

@dengemann comments addressed. FYI you asked about speed, and that made me optimize a bit more:

2 (320) 1 C 0.283627033234
2 (320) 1 Python 0.291443109512
2 (320) 3 C 0.465950012207
2 (320) 3 Python 2.45959687233
3 (1280) 1 C 0.660896062851
3 (1280) 1 Python 1.62492084503
3 (1280) 3 C 12.4721641541
3 (1280) 3 Python 11.3650319576
4 (5120) 1 C 24.9352688789
4 (5120) 1 Python 11.2639200687
4 (5120) 3 C 644.138396025
4 (5120) 3 Python 118.89906311

Reply to this email directly or view it on GitHub.

v2 = tri_rr[np.newaxis, 1, :] - fros
v3 = tri_rr[np.newaxis, 2, :] - fros
triples = _fast_cross_nd_sum(v1, v2, v3)
l1 = np.sqrt(np.sum(v1 * v1, axis=1))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is it faster than np.linalg.norm( ..., axis=1) ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes sadly (and axis is only available on 1.9+), I tested both. I think it depends on your application. Probably for larger matrices norm(..., axis) would be faster, we'll have to decide on a case-by-case basis.

@agramfort
Copy link
Member

please update what's new, rebase and +1 for merge !

thanks awesome !

@larsoner
Copy link
Member Author

Rebased and updated, I'll merge once Travis comes back happy

@agramfort
Copy link
Member

agramfort commented Jul 13, 2015 via email

larsoner added a commit that referenced this pull request Jul 13, 2015
MRG: BEM model and solution in Python
@larsoner larsoner merged commit f60a18a into mne-tools:master Jul 13, 2015
@larsoner larsoner deleted the bem branch July 13, 2015 20:01
@larsoner
Copy link
Member Author

Thanks for the reviews @dengemann @agramfort

@dengemann
Copy link
Member

Super cool, Eric -- amazing!

On Mon, Jul 13, 2015 at 10:01 PM, Eric Larson notifications@github.com
wrote:

Thanks for the reviews @dengemann https://github.com/dengemann
@agramfort https://github.com/agramfort


Reply to this email directly or view it on GitHub
#2193 (comment)
.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants