Skip to content

Conversation

@jjnurminen
Copy link
Contributor

@jjnurminen jjnurminen commented May 3, 2016

Discussed in #3097

Summary: Elekta/Neuromag DACQ (data acquisition) supports rather flexible event and
averaging logic that is currently not implemented in mne-python. It also stores all averaging
parameters in the fiff file, so raw data can be easily reaveraged in postprocessing.
The purpose of this PR is

  1. extract the relevant info from the fiff file
  2. implement support for averaging according to DACQ categories (or to modify the categories first)

Implementation: a class that takes all the relevant info from raw.info['acq_pars'].

API:
get_condition() gives averaging info for a given condition, that can be fed to mne.Epochs

Technical details: DACQ supports defining 32 different events which correspond to trigger
transitions. Events support pre- and post-transition bit masking. Based on the events, 32 averaging categories can be defined. Each category defines a reference event that is the zero-time point for collecting the corresponding epochs. Epoch length is defined by the start and end times (given relative to the reference event). A conditional ('required') event is also supported; if defined, it must appear in a specified time window (before or after the reference event) for the epoch to be valid.

mne/event.py Outdated
'eog':self.eogmax, 'ecg':self.ecgmax}

def get_epochs(self, raw, catname, picks=None, reject=None, stim_channel=None, mask=0):
""" Get mne.Epochs instance corresponding to the given category. """
Copy link
Member

Choose a reason for hiding this comment

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

Nest the import here. In other words, put this just below this line:

    from .epochs import Epochs

Then delete the import from the top. This will solve the circular import problem you're having.

@agramfort
Copy link
Member

can you run flake8 tool on your code and follow the pep8 guidelines?
you could put the code in an elekta_events.py file for now.

@coveralls
Copy link

coveralls commented May 4, 2016

Coverage Status

Coverage decreased (-0.2%) to 90.651% when pulling e9a5353 on jjnurminen:elekta_avg into f932004 on mne-tools:master.

@jjnurminen
Copy link
Contributor Author

jjnurminen commented May 9, 2016

I dug a little into the TRIUX/Vectorview issue.
In my current understanding, the code is already compatible with both TRIUX and Vectorview data, as long as it was measured with DACQ version => 3.4 (released in 2005). According to the fiff info, the sample_audvis_raw.fif example data was measured in 2002, so it's pretty old. For testing the averager, I would definitely like to use newer data.
I can try to implement support for DACQ < 3.4 if I can get the specs from Elekta, but I'm not sure how many people will actually be processing such ancient files. ;)

@coveralls
Copy link

Coverage Status

Coverage decreased (-10.8%) to 79.988% when pulling 77fa95e on jjnurminen:elekta_avg into f932004 on mne-tools:master.

@jjnurminen
Copy link
Contributor Author

jjnurminen commented May 9, 2016

I created a fiff file (95 megs) with various event/category definitions and some environmental interference acting as "signal". With that file, the code above now seems to give results identical to the MaxFilter averager (except for some baseline effect, haven't yet figured out where it comes from). I guess this fiff could be used for the unit tests eventually.

@coveralls
Copy link

coveralls commented May 9, 2016

Coverage Status

Coverage decreased (-0.1%) to 90.694% when pulling 31e2494 on jjnurminen:elekta_avg into f932004 on mne-tools:master.

@larsoner
Copy link
Member

larsoner commented May 9, 2016

I guess this fiff could be used for the unit tests eventually.

Sure. Maybe to keep file size down you could do raw.pick_channels(...) and pick a few MEG channels, plus the STI channel(s) you need. That plus the -ave.fif from Elekta's program would make for a good set of unit tests.

@jjnurminen
Copy link
Contributor Author

I still get a DC shift compared to the Elekta averager, so my unit tests don't pass at the moment. I'm waiting for response from Elekta to confirm that their code really doesn't do any baseline or filtering operations.

@larsoner
Copy link
Member

Have you tried putting in a baseline to see if it gets closer? You could also compute the PSD of the two datasets to get some evidence of whether or not filtering has been applied

@jjnurminen
Copy link
Contributor Author

@Eric89GXL If I apply a pre-stimulus baseline to the epochs before averaging, I get almost 100% match with Elekta data, but Elekta claims that they don't use any baseline in their averager.

BTW, do you know why the audvis_raw example data does not have info['acq_pars'] at all? I guess it was somehow modified and rewritten? My goal would be to implement support also for pre-3.4 DACQ versions (older Vectorview systems), so I guess more testing data would be needed for that. I can try to find some old data from our lab.

@larsoner
Copy link
Member

I get almost 100% match with Elekta data, but Elekta claims that they don't use any baseline in their averager.

Our precedent in cases like this is to put some indignant comment like "# XXX this baseline is necessary to get a good match even though the Elekta averager isn't supposed to baseline correct???" and do what is necessary to get a good match :) That way it won't hold up the PR while we wait for clarification.

The sample data are old, and from MGH. Maybe they have a special setup...? I don't think anyone manually removed it. The best thing is probably to find data from your lab.

@jjnurminen
Copy link
Contributor Author

jjnurminen commented May 27, 2016

OK. Just to be sure, can you confirm my understanding that the following sequence of operations should perform no baselining, detrending or filtering whatsoever:

  1. raw = mne.io.raw(filename)
  2. eps = mne.Epochs(raw [...], baseline=None, detrend=None)
  3. avs = eps.average()

@larsoner
Copy link
Member

Yep, that should do it. You can confirm by doing avs.plot(), you'll see the butterfly plot drifts

@jjnurminen
Copy link
Contributor Author

Yep. The funny thing is that the MaxAve data is also shifted, but by a different amount. Applying baseline to both makes them identical. Anyway, evoked responses are usually baselined before analysis, so I guess it's not a big deal. I'll just write the unit tests for the baselined data as you suggested.

@larsoner
Copy link
Member

Don't forget the angry/confused comment, too :) But seriously put in a comment with XXX in it somewhere, it's often how we flag things that are mysterious and can hopefully be fixed later

@jjnurminen
Copy link
Contributor Author

Will do.

@jjnurminen
Copy link
Contributor Author

There's now an unit test that averages all categories from the raw data and compares to Elekta averaged data to within 1 fT(/cm). I truncated the raw file to 3 MEG channels.

@larsoner
Copy link
Member

Can you open a PR to add the data file to the mne-testing-data repository?

@larsoner
Copy link
Member

I guess you'll probably add a raw and an evoked file?

@jjnurminen
Copy link
Contributor Author

jjnurminen commented May 31, 2016

Sure, but I'm a bit confused about where the files should go? The test files under mne/io/tests/data/ don't seem to be in the mne-testing-data repository.

@agramfort
Copy link
Member

agramfort commented May 31, 2016 via email

@jjnurminen
Copy link
Contributor Author

@agramfort OK. Should I put it under /MEG/sample then?

@larsoner
Copy link
Member

larsoner commented May 31, 2016 via email

@jjnurminen
Copy link
Contributor Author

Created PR at mne-tools/mne-testing-data#11

mne/event.py Outdated
return events_out


class Elekta_event(object):
Copy link
Member

Choose a reason for hiding this comment

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

ElektaEvent

Copy link
Member

Choose a reason for hiding this comment

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

Actually, this class is pretty thin. Can you change it just to use a dict?

Copy link
Contributor Author

@jjnurminen jjnurminen May 31, 2016

Choose a reason for hiding this comment

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

Can think about it. I would prefer classes, since there is basically no cost and they are easier to expand later. Also if the user wants to create his own categories/events, sanity checking is much easier inside a class. (I would like to make some API expansions later so that one can easily change averaging parameters or create new events)

Copy link
Member

Choose a reason for hiding this comment

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

I think all these classes should be private if possible -- people should never need to see them if we can avoid it. But maybe I'm missing some future API design plans.

Copy link
Member

Choose a reason for hiding this comment

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

found this from an old tweet by @agramfort : https://www.youtube.com/watch?v=o9pEzgHorH0 :) we shouldn't overfit for future design plans ...

@jjnurminen
Copy link
Contributor Author

jjnurminen commented Sep 20, 2016

Thanks @jasmainak, I understood that python_reference.rst should be ok already

@jjnurminen
Copy link
Contributor Author

Not sure what the coverage error is about - should there be more tests?

@larsoner
Copy link
Member

larsoner commented Sep 20, 2016 via email

mne/event.py Outdated
s += '>'
return s

def __getitem__(self, items):
Copy link
Member

Choose a reason for hiding this comment

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

else:
raise KeyError('No such category')

def __len__(self):
Copy link
Member

Choose a reason for hiding this comment

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

Same here

@larsoner
Copy link
Member

Example looks good

https://2598-1301584-gh.circle-artifacts.com/0//home/ubuntu/mne-python/doc/_build/html/auto_examples/io/plot_elekta_epochs.html#sphx-glr-auto-examples-io-plot-elekta-epochs-py

Can you add your new dataset to _download_all_datasets in datasets/utils.py? It will prevent the progress bar junk in example rendering in master.

If you want to venture into bonus territory you can also modify circle.yml to also download it for PRs if necessary based on a regex like we do for sample and some BS datasets.

But otherwise LGTM, thanks for sticking with so many revisions. The next PR should be easier :₩

@larsoner
Copy link
Member

That was supposed to be a smiley but apparently I mistyped on my phone

@larsoner
Copy link
Member

Looking at coverage, the problems are __getitem__, __repr__, events
property, and four lines related to self.compat, at least the first three
should be easily fixable in tests

On Sep 20, 2016 7:30 AM, "Eric Larson" larson.eric.d@gmail.com wrote:

Take a look at the report and see which lines you missed. You might need
some assert_raises to hit exceptions, for example. Try to get over 90%
coverage. Sometimes it seems not to update the comment or status properly,
so be sure to check the site. And failing builds, if you have them, do not
contribute coverage, so those need to pass first (on mobile so can't see
directly).

@jjnurminen
Copy link
Contributor Author

But otherwise LGTM, thanks for sticking with so many revisions. The next PR should be easier

Thanks to everybody (and especially Eric) for devoting a lot of time and effort to this.
BTW, only now I noticed the 'Allow edits from maintainers' button. I wish I had found it six months ago, then you guys would have fixed everything ;)

@jjnurminen
Copy link
Contributor Author

Now it doesn't seem to run coverage checks, so I'm not sure what's the situation.

@larsoner
Copy link
Member

It reported on your last commit -- it just waits until all Travis runs are done. You have 87% coverage on changed lines, but the misses are mostly in datasets and self.compat so it's fine (unless you can easily touch the self.compat lines)

@jjnurminen
Copy link
Contributor Author

Ok. I'm not sure if compat needs more tests. Of course I could e.g. load the multimodal dataset and see that compat gets switched on, but I wonder if it justifies the hassle and increased testing time.

@larsoner
Copy link
Member

Yeah definitely not using the multimodal dataset, but wouldn't our sample (and therefore testing, which is a truncated version of that) have it, too, since it's old?

@jjnurminen
Copy link
Contributor Author

Yep. Tested in latest commit.

@larsoner larsoner changed the title MRG: Elekta averager MRG+1: Elekta averager Sep 21, 2016
@larsoner
Copy link
Member

LGTM +1 for merge. @agramfort feel free to merge if you're happy

@agramfort agramfort merged commit 577d5fa into mne-tools:master Sep 24, 2016
@agramfort
Copy link
Member

thx @jjnurminen

@jjnurminen
Copy link
Contributor Author

Thanks!

On 24 Sep 2016 23:06, "Alexandre Gramfort" notifications@github.com wrote:

thx @jjnurminen https://github.com/jjnurminen


You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
#3205 (comment),
or mute the thread
https://github.com/notifications/unsubscribe-auth/AIGRObiQI9MZsZi7J9sfpfcIqITde3-Mks5qtYKUgaJpZM4IWNS6
.

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.