Add stim demultiplex rhs intan

[h-mayorquin]

Should fix #1649

We will be needing test data but this work for reading the data for me:

from spikeinterface.extractors import IntanRecordingExtractor

recording = IntanRecordingExtractor(file_path, stream_name="Stim channel")

recording.get_traces()

If you can give it a try @TommasoLambresa and @FrancescoNegri and check that the outputs matches what you expect it would be great.

@zm711

[h-mayorquin]

No idea why the docs are failing on this one.

[zm711]

No idea why the docs are failing on this one.

We've been having a lot of failures of gin lately. it's intermittent (3.9, 1.26 testing failed too because of that). I'm rerunning now :)

[zm711]

First round of comments.

[TommasoLambresa]

If you can give it a try @TommasoLambresa and @FrancescoNegri and check that the outputs matches what you expect it would >be great.

@h-mayorquin, @zm711
We're still here, sorry for the delay. We moved on this branch pointing to this PR, but running the lines that you sent us nothing changes compared with the previous neo version. We tried to debug, but it never enters in the function that you have modified: when spikeinterface calls the neorawio parse_header function, we see that the funtion read_rhs is called. In this function there is no reference to the function modified in this PR, so we're asking: how can we test the changes?
We tried also to print some random text but nothing happens.

[zm711]

Is your neo install editable?

If you just installed neo during your spikeinterface install it would not be editable. So even if you switched branches it wouldn't change anything. You would need to editable install neo in the environment that you're using for spikeinterface to be able to test changes. Let us know if you need help with that!

[TommasoLambresa]

Yes, we installed neo editable! We are also sure that the neo version in environment is this newer one. @zm711
Here the steps that we adopted:

1. Cloning the repo h-mayorquin/python-neo
2. Checkout on the add_intan_stim branch
3. Create an environment with spikeinterface
4. Moving within the local neo folder and digiting pip install -e .

[zm711]

You would need to request the "Stim channel" stream_name or stream_index 11. Did you do that from spikeinterface. I would recommend trying this from neo directly. EDITED for correct stream name

so,

from neo.rawio import IntanRawIO

reader = IntanRawIO(path/to/data)
reader.parse_header()
sig_chunk = reader.get_analogsignal_chunk(stream_index=11)

[TommasoLambresa]

By trying these lines, it raised this error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[6], line 5
      3 reader = IntanRawIO("Z:/SNN/K-SNN/R24-21/02_stim_240621_132910.rhs")
      4 reader.parse_header()
----> 5 sig_chunk = reader.get_analogsignal_chunk(stream_index=11)
      7 import matplotlib.pyplot as plt
      8 plt.plot(sig_chunk.T)

File ~\Repositories\test_NEOIntan\mayorquin-neo\neo\rawio\baserawio.py:805, in BaseRawIO.get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, stream_index, channel_indexes, channel_names, channel_ids, prefer_slice)
    799     error_message = (
    800         "get_analogsignal_chunk can't be called on a file with no signal streams or channels."
    801         "Double check that your file contains signal streams and channels."
    802     )
    803     raise AttributeError(error_message)
--> 805 stream_index = self._get_stream_index_from_arg(stream_index)
    806 channel_indexes = self._get_channel_indexes(stream_index, channel_indexes, channel_names, channel_ids)
    808 # some check on channel_indexes

File ~\Repositories\test_NEOIntan\mayorquin-neo\neo\rawio\baserawio.py:653, in BaseRawIO._get_stream_index_from_arg(self, stream_index_arg)
    651 else:
    652     if stream_index_arg < 0 or stream_index_arg >= self.header["signal_streams"].size:
--> 653         raise ValueError(f"stream_index must be between 0 and {self.header['signal_streams'].size}")
    654     stream_index = stream_index_arg
    655 return stream_index

ValueError: stream_index must be between 0 and 3

As you can see the functions are called from my local repository, which contains the neo updated version. It seems it doesn't recognize the file as containing a stim stream, but I'm sure it's wrong!

The previous code from spikeinterface was already with the specified stream_name = "Stim channel":

from spikeinterface.extractors import IntanRecordingExtractor

recording = IntanRecordingExtractor("Z:/SNN/K-SNN/R24-21/02_stim_240621_132910.rhs", stream_name="Stim channel")
data = recording.get_traces()

Actually, these lines didn't raised an error but simply returned the correct signals without converting from the binary format.

[zm711]

Sorry I gave you the stream_id instead of the stream index. That was my bad! This is why Heberto and I would prefer if this function accepted the id instead of the index.

Going back to your spikeinterface. You mean the data wasn't demultiplexed. It is still one array instead of multiple arrays? We need a scaling step to switch to current which spikeinterface doesn't do automatically. Let me rewrite my neo script for you :)

[zm711]

from neo.rawio import IntanRawIO

reader = IntanRawIO(path/to/data)
reader.parse_header()
signal_streams = reader.header['signal_streams']
stream_ids = signal_streams['id'].tolist()
stream_index = stream_ids.index('11')
sig_chunk = reader.get_analogsignal_chunk(stream_index=stream_index)
final_stim = reader.rescale_signal_raw_to_float(sig_chunk, stream_index=stream_index)

Sorry kept hitting enter accidentally. Try this. The key is you also need to rescale which doesn't occur with the raw data.

[h-mayorquin]

Hi, @TommasoLambresa

Using this branch, this is what I get for the file that you shared with me:

from spikeinterface.extractors import IntanRecordingExtractor
from matplotlib import pyplot as plt


recording = IntanRecordingExtractor(file_path, stream_name="Stim channel")

traces = recording.get_traces(return_scaled=True)
times = recording.get_times()
num_channels = traces.shape[1]
for i in range(num_channels):
    channel_traces = traces[:, i]
    if channel_traces.sum() > 0:
        plt.plot(times, channel_traces, label=f"Channel {i}")
        
plt.legend()
plt.xlabel("Time (s)")
plt.ylabel("Amperes")

[TommasoLambresa]

Hi @h-mayorquin
The same goes for me, but unfortunately it is not the correct conversion. This is basically the same trace got without the binary conversion (it is just scaled accordingly to some gain and offset). I'm looking what goes wrong in this conversion function that you wrote. The correctly decoded traces contains biphasic anodic-cathodic stimuli of $\pm 60 \mu A$. Debugging within your function it seems the conversion is not made properly. I'm trying to understand what's wrong.

The same goes from the neo side, since the conversion function is the same. @zm711

[h-mayorquin]

Hi @h-mayorquin The same goes for me, but unfortunately it is not the correct conversion. This is basically the same trace got without the binary conversion (it is just scaled accordingly to some gain and offset). I'm looking what goes wrong in this conversion function that you wrote. The correctly decoded traces contains biphasic anodic-cathodic stimuli of ± 60 μ A . Debugging within your function it seems the conversion is not made properly. I'm trying to understand what's wrong.

Thanks a lot, this is very helpful and precisely the reason why I was waiting for someone to implement this. Two questions:

1. Can you reproduce the figure that I have? (this is to see that we can get you the code to make debugging faster in the future)
2. Can you share the code that you use for decoding with me? Maybe that way I can spot the difference and what I did wrong.

[TommasoLambresa]

1. Can you reproduce the figure that I have? (this is to see that we can get you the code to make debugging faster in the future)

Yes, I have the same result:

from spikeinterface.extractors import IntanRecordingExtractor
from matplotlib import pyplot as plt


recording = IntanRecordingExtractor("C:/Users/t.lambresa/Downloads/newintanFile 1.rhs", stream_name="Stim channel")

traces = recording.get_traces(return_scaled=True)
times = recording.get_times()
num_channels = traces.shape[1]
for i in range(num_channels):
    channel_traces = traces[:, i]
    if channel_traces.sum() > 0:
        plt.plot(times, channel_traces, label=f"Channel {i}")
        
plt.legend()
plt.xlabel("Time (s)")
plt.ylabel("Amperes")

2. Can you share the code that you use for decoding with me? Maybe that way I can spot the difference and what I did wrong.

Sure here is the code @h-mayorquin:

def decode_stim_data(orig_data, stim_step_size):
    """Converts the stimulation data from the Intan RHS format to micro-amps.
    Parameters
    ----------
    stim_data : np.ndarray (dtype=uint16)
        One-dimensional array of STIM samples read from the .rhs file.
    stim_step_size : float
        'Stim Step Size' value read from the header (in amperes per step).
    Returns
    -------
    stim_curr : np.ndarray (float)
        Current in micro-amps (positive or negative).
    """
    stim_data = orig_data.copy()
    compliance_limit_data = stim_data >= 2**15
    stim_data = stim_data - (compliance_limit_data * 2**15)
    charge_recovery_data = stim_data >= 2**14
    stim_data = stim_data - (charge_recovery_data * 2**14)
    amp_settle_data = stim_data >= 2**13
    stim_data = stim_data - (amp_settle_data * 2**13)
    stim_polarity = stim_data >= 2**8
    stim_data = stim_data - (stim_polarity * 2**8)
    stim_polarity = 1 - 2 * stim_polarity  # Convert (0 = pos, 1 = neg) to +/-1
    stim_data = stim_data * stim_polarity
    # Convert to amps
    stim_curr = stim_data * stim_step_size
    return stim_curr, compliance_limit_data, charge_recovery_data, amp_settle_data

[h-mayorquin]

Ok, can you confirm this is what I should get for channel 112:

[TommasoLambresa]

Yes, now it works :)

[zm711]

@h-mayorquin wasn't our whole discussion for spikeinterface that we weren't going to allow people to do this through spikeinterface (without going through the scale mechanism)?

[h-mayorquin]

Yeah, indeed. The scaling on the first code that I shown is wrong because it will put the scaling as one. There is another error though, I am fixing it.

[TommasoLambresa]

I think the issue could be that the signal is casted as uint but it is int actually. And the exadecimal to take the first 8 bits is 0x7F10, taking into account that the first 8 bits codify the magnitude. By using 0xFF you consider only the last 8 bits if I'm not wrong. @h-mayorquin

[zm711]

Hey @TommasoLambresa I would buy you're right on the bit issues (I'm terrible on doing hex etc which is why I trust Heberto for that. But it is uint16 :)

Unless you think Intan miswrote.

[TommasoLambresa]

No sorry, you're right! I'm terrible too...

[zm711]

We will get this all sorted :)

Like I said in my note, spikeinterface is moving toward limiting non-voltage scaling in the get_traces function. So there will be a workaround to do what you want. That is why I was having you try to test in neo first. We allow for any units for scaling. Then Heberto will make sure it scales correctly in neuroconv. But spikeinterface will require an extra few lines of code to scale correctly. I'll wait until Heberto posts his fix then I would recommend you both test in neo first :)

[TommasoLambresa]

Perfect! Thank you guys @zm711 @h-mayorquin.

[h-mayorquin]

Here is the gin PR:

https://gin.g-node.org/NeuralEnsemble/ephy_testing_data/pulls/151

[zm711]

gin file is merged and can be added to testing.

[h-mayorquin]

I added a test, here is the data:

And I am testing that the first pulse matches the description that the authors gave us to ensure that the decoding is done properly

[zm711]

We may have discovered an issue with nixio by updating the io... let me look into that for a moment before we merge this.

[zm711]

looks like setuptools just did a release five hours ago. I think that new release is breaking our nixio install. I'll see if I can limit setup tools for now, but I'll be away from my computer for a minute so I'll have to debug this more tomorrow.

[h-mayorquin]

G-Node/nixpy#552

The last version of setuptools fixes it so just set a ceiling for setuptools

https://setuptools.pypa.io/en/stable/history.html?utm_source=chatgpt.com#deprecations-and-removals

(they probably went back after they realize a lot of people just won't upload and got people yelling at them x D)

[zm711]

I would push to comment these two spots out for the moment. We can add this later once we have the test data that way we can get this merged into intan.

[zm711]

@h-mayorquin could you scan through my changes really quick just to make sure I didn't do anything stupid!

[zm711]

@h-mayorquin I added in the dc only test folder, because it has stim file tests and it tests the fact that the end user can choose not to save normal amp files this will allow to test one-file-per-channel for dc amp and for stim file. The other test is nice but not necessary to get this merged :)

[zm711]

@h-mayorquin I marked three comments of what I want you to glance at before we merge this. The fact that tests are passing means that my assumptions of the data org appear to be correct, but if you can think of a clearer/better way to do what I've done feel free to tweak. Once you glance over my additions this is ready for merge.

[zm711]

@h-mayorquin does this make sense or can you think of a better way of doing this?

[h-mayorquin]

I don't understand the context. Why are we doing this? I re-read

        # if using the `header-attached` or `one-file-per-signal` formats
        # the amplifier data is stored in the same place in the header so no matter what the DC amp
        # uses the same info as the RHS amp

There is something strange about it and makes me confused.

[zm711]

This is the only change here (the diff is wrong). Basically if the user shut off saving of files in the file-per-channel we need to skip this section. Does this make sense.

[zm711]

@h-mayorquin you want to push your last change with our patch here? Then I merged. We can add the rhs channel test either now or in a new PR.

[h-mayorquin]

Merge it like this, my comments are a mess.

[zm711]

For this we check if these exist and if they do we analyze. This was my hacky way of doing this (since the user can choose individual channels to stim (ie stim files only exist for a subset of amp channels). Check this for me.

[zm711]

Also fine if you want to do your reorg for this if statement before we merge :)

[h-mayorquin]

Na, same thing, this is working, I am afraid of breaking something and I want to be sure we have some sort of support on the next release.