Refactor artifacts detection + add sturation detection

[samuelgarcia]

@JoeZiminski @oliche

This refactor the artifact detection and blanking

• with a main multi method function : detect_artifact_periods()
• return periods with the new base_period_dtype instead list of list
• SilencedPeriodsRecording also use this new dtype
• DetectArtifactAndSilentPeriodsRecording (replace SilencedArtifactsRecording) do detect and silence at the same
  time to be used in preprocess dict (@chrishalcrow this is for you)
• refactor the method of Pierre based on envleope
• incorporate the method of @oliche and @JoeZiminski to detect saturation

Could also be done in this PR:

• add a margin_ms option in SilentPeriodsRecording to ewxtend the zeroing with a mragin : do we need this ?

To be done in 2 futures PR:

• handle margin and taper on border of SilentPeriod (more tricky that it appear when period is at border...)
• rewrite RemoveArtifactsRecording based on SilencedPeriodsRecording and the new dtype

[yger]

@samuelgarcia do not forget, in this PR, to change DetectThresholdCrossing(recording, ...) to DetectThresholdCrossing(envelope, ...) in the detec_period_artefacts_by_envelope

[samuelgarcia]

Hi @oliche @JoeZiminski
I copy paste and clean a bit your PR #4301 into this one.
Have a look and tell if this is OK.

I think we do not need to scaled traces at each chunk we could only could the unscaled threhold in the init and then we could avoid the costly convertion to float32 or float64.

Also I put the unit in uV to be more SI friendly.

[samuelgarcia]

@alejoe91 @chrishalcrow ready to review on my side

[samuelgarcia]

add detect_artifact_method detect_artifact_kwargs

[chrishalcrow]

Hello, trying this out on real data with saturation. I get this output

artifact_periods=array([(0, 637721, 638769), (0, 638955, 639206), (0, 923219, 923602),
       (0, 923605, 923606), (0, 923607, 923835)],

Are these periods inclusive? And if they are, should the last pair be merged?

[JoeZiminski]

Hey @samuelgarcia thanks for this! will check this out tomorrow. I accidentally did not push two tidy-up commits to #4301, do you mind if I push to your branch to add some changes? I think it is mostly superficial stuff like docstrings etc.

[samuelgarcia]

Hello, trying this out on real data with saturation. I get this output

artifact_periods=array([(0, 637721, 638769), (0, 638955, 639206), (0, 923219, 923602),
       (0, 923605, 923606), (0, 923607, 923835)],

Are these periods inclusive? And if they are, should the last pair be merged?

Yes, the last pair should be merged.
Let us check. the _collapse_events() is done for this purpose.
@JoeZiminski can you look at this ?

[samuelgarcia]

Hey @samuelgarcia thanks for this! will check this out tomorrow. I accidentally did not push two tidy-up commits to #4301, do you mind if I push to your branch to add some changes? I think it is mostly superficial stuff like docstrings etc.

yes of course. push anything you want

[JoeZiminski]

Hey @samuelgarcia thanks for this, I added the content of those commits I missed. I ran into some issues with the (now renamed) uV_per_second_threshold . This is a threshold for the first derivative in units uV/s so it can be a very small number (e.g. IBL default in volts, converted to uV is 1e-2; for example, with the test data the threshold is around 0.0004 which goes negative after the round trip to int and back. Please correct me if wrong, the dynamic range of NP probes pm 5mV at int16 the bit resolution is 0.153 μV. I think this might be why @chrishalcrow is seeing those missed samples in his data.

If possible I'd prefer to cast and scale the data to avoid such issues, then we always know we are solid whatever the user inputs. I get 0.047298s to cast, offset and scale 10 million samples so I don't think this will be a bottleneck. WDYT?

[samuelgarcia]

@JoeZiminski : why did you put 50 as offset here ?

[JoeZiminski]

Apologies forgot to say, just to ensure it is tested (if it was zero and was accidentally not applied in the function, the test would not flag it)

[JoeZiminski]

(that was just a quick test, see update below)

[samuelgarcia]

Please correct me if wrong, the dynamic range of NP probes pm 5mV at int16 the bit resolution is 0.153 μV. I think this might be why @chrishalcrow is seeing those missed samples in his data.

NP1 is 12bits, so I think that the lsb is 5000 / 2**11 = 2.4uV
Am I wrong ?

[JoeZiminski]

Hey @samuelgarcia I didn't realize the bit depth of the ADC is not 16-bit! I see here ADC resolution is given as 10 bits. In any case, the resolution from the gains/offsets used to scale the data won't be able to represent the values used for the uV_per_second_threshold?

[oliche]

@JoeZiminski the sample to volts conversion is a probe dependent variable. To get a feel of the complexity of the matter, here are the Neuropixel possible configs for all probes released to date: https://docs.google.com/spreadsheets/d/1UZycmxwwaJGCJ4AknyXrHkDXWOUO24TtAnBbe1mQ9u4/edit?gid=0#gid=0

As you can see there are potentially multiple conversion factors depending on the probe vintage.

It seems the SpikeInterface canon is to leave this up to the user. Our design in https://github.com/int-brain-lab/ibl-neuropixel has been to read it from the metadata, making sure we are compatible with the oldest probe versions we observed at IBL.

[JoeZiminski]

Thanks @oliche, in this case @samuelgarcia I think we should definitely scale the data rather than the parameters, as we cannot be sure the gains / offsets used to scale the data will not destroy information when converting float-to-int e.g. uV_per_second_threshold. We could play around with applying the division by fs at another point but it seems a little brittle and not worth the very minor speed-up.

[JoeZiminski]

Hey @samuelgarcia @alejoe91 thanks for the feedback, I'm running into the below issue calculating the derivative on int16 data. I don't think this is an artefact of the test data but would be good to check, @oliche too.

See below for the test data in float and int. When taking the derivative, information is lost due to overflow.

The above is the data converted to int, each timeseries is the data from one channel. In many instances, the data after the saturation is negative in int space, resulting in overflow

e.g. by definition the saturation point is at the max of the int16 range, 32767. Often, when the saturation goes back to baseline, the datapoint is negative in int space e.g. -10000. When you take the diff e.g. -10000 - 32767 you get overflow and the information is lost. It seems like this would happen in real data as well?

[oliche]

@JoeZiminski I think the test data is not realistic. Usually the maximum int recorded on neuropixel is 512 on NP1 and 8192 on NP2 so you wouldn't have the issues above as the max diff would be 8192 * 2 (2 ** 14).

However in general I think it is a poor practice to compute in int16 for this kind of issues and patchy support of the numpy suite. Users may unwittingly create weird situations, exactly like this one: #4175.

Anyways for this specific saturation detection problem, you just have to fix your test data by staying in the maxint ranges above.

[JoeZiminski]

Cheers @oliche! That's good to know for NP probes. For other providers like Nexus, run their ADC at 16 bit-depth (e.g. here), in this case it might be a problem? That being said I don't know how often these probes saturate and whether people are likely to run this function for non-NP1 probes. But yes I think the easiest solution to avoid these worries is to perform computation in float.

[oliche]

If you think you may exceed half the the range of int16 (ie higher than 2 ** 14 or lower than -2 ** 14), then you can check for the max(abs()) values. If they are above or below, you may throw a warning and only perform the absolute saturation detection and skip the one based on the derivative ?

That would prevent too many false positives.

[JoeZiminski]

Thanks @oliche, yes in this case I would definitely suggest @samuelgarcia we perform operations in float here to avoid having to worry about the unpredictable bit-depth of the data.

[alejoe91]

Hi guys, I've debugged this a bit and I think it's much simpler to specify the diff threshold in diff_threshold_uV. This would be the threshold in the uV per sample, for example, if diff_threshold_uV=200 it means that if the next sample has a diff larger than 200uV, the current sample is flagged as an artifact. I think this makes it much easier to grasp, with respect to using uv_per_seconds ot uv_per_ms.

Thoughts?

[JoeZiminski]

Hey @alejoe91 thanks I completely agree that is much more intuitive and easy to use. Thanks for finalising this! My only comments are:

1. The traces should be cast to float before the derivative threshold is applied (to avoid overflow issues as above). As far as I can tell the traces are currently passed as int16 by default to compute

2. I am not 100% convinced that the saturation detection and envelope detection should be merged into a single function. What is the use case of the envelope detection / when should it be applied? Is it similar to the saturation detection in that it should be computed before any preprocessing has been performed and applied after? If the use case are exactly the same and should be applied in exactly the same way, and these are just different methods to achieve the same thing, then I can see the benefit.

   But otherwise, I cannot see the benefit, and it makes things confusing as a user. I am a 'user' interested in saturation detection and now I have spent 20 minutes researching envelope detection to see if I should be using that instead (which I have just done, and I'm still not sure). Passing a dictionary of kwargs that changes based on the value of another argument is also quite confusing, and you loose IDE autocomplete and it makes documentation difficult. I have not seen this pattern before in common frameworks, are there any other examples of it? If all three functions are made publically available and documented its not so bad, but in that case I cannot see the marginal benefit of the detect_artifact_periods function.

[alejoe91]

Hey @alejoe91 thanks I completely agree that is much more intuitive and easy to use. Thanks for finalising this! My only comments are:

1. The traces should be cast to float before the derivative threshold is applied (to avoid overflow issues as above). As far as I can tell the traces are currently passed as int16 by default to compute
2. I am not 100% convinced that the saturation detection and envelope detection should be merged into a single function. What is the use case of the envelope detection / when should it be applied? Is it similar to the saturation detection in that it should be computed before any preprocessing has been performed and applied after? If the use case are exactly the same and should be applied in exactly the same way, and these are just different methods to achieve the same thing, then I can see the benefit.
   But otherwise, I cannot see the benefit, and it makes things confusing as a user. I am a 'user' interested in saturation detection and now I have spent 20 minutes researching envelope detection to see if I should be using that instead (which I have just done, and I'm still not sure). Passing a dictionary of kwargs that changes based on the value of another argument is also quite confusing, and you loose IDE autocomplete and it makes documentation difficult. I have not seen this pattern before in common frameworks, are there any other examples of it? If all three functions are made publically available and documented its not so bad, but in that case I cannot see the marginal benefit of the detect_artifact_periods function.

Hey @alejoe91 thanks I completely agree that is much more intuitive and easy to use. Thanks for finalising this! My only comments are:

1. The traces should be cast to float before the derivative threshold is applied (to avoid overflow issues as above). As far as I can tell the traces are currently passed as int16 by default to compute

Done! Just casting to float at the top of the compute!

2. I am not 100% convinced that the saturation detection and envelope detection should be merged into a single function. What is the use case of the envelope detection / when should it be applied? Is it similar to the saturation detection in that it should be computed before any preprocessing has been performed and applied after? If the use case are exactly the same and should be applied in exactly the same way, and these are just different methods to achieve the same thing, then I can see the benefit.
   But otherwise, I cannot see the benefit, and it makes things confusing as a user. I am a 'user' interested in saturation detection and now I have spent 20 minutes researching envelope detection to see if I should be using that instead (which I have just done, and I'm still not sure). Passing a dictionary of kwargs that changes based on the value of another argument is also quite confusing, and you loose IDE autocomplete and it makes documentation difficult. I have not seen this pattern before in common frameworks, are there any other examples of it? If all three functions are made publically available and documented its not so bad, but in that case I cannot see the marginal benefit of the detect_artifact_periods function.

I think I agree with this. We can keep it simple, have the 2 functions and users can directly choose what they want to do. @samuelgarcia ok to remove detect_artifact_periods?

[alejoe91]

Another note: the last commits also add the option to annotate a recording with saturation_threshold_uV and use that automatically when set. In a follow up, we cna port the logic from ibl-neuropixel to set this automatically based on metadata. We'll have to see how to do it for Open Ephys as well

[alejoe91]

I suggest making this optional and have a check_for_peak_source flag

[alejoe91]

@JoeZiminski can you check this 300uV/sample value?