nelpy / nelpy Goto Github PK

If we make a SpikeTrain object, we assume a single unit, whereas if we make a SpikeTrainArray object we assume multiple units. However, we may want to merge both into one class, namely SpikeTrain that can have multiple units.

Input validation becomes tricky though, since we support taking in lists, arrays, and generally any array-like object.

How can we (robustly) differentiate between st1 = [[2,3,4],[3,4],[5,7,8]] and st2 = [3,6,8] ? In both cases, len(st) == 3 and doing things like np.squeeze(st) doesn't get us much further, since then

st1.ndim == 1 and st2.ndim == 1 and st1.size = 3 and st2.size == 3 and even st1.shape == (3, ) and st2.shape == (3, ).

@shayokdutta any ideas?

epocharray = EpochArray([[3,4],[5,8],[10,12], [16,20], [22,23]])

epocharray[::2].time works

but

epocharray[::3].time

returns

array([[ 3, 16],
       [ 4, 20]])

instead of

array([[ 3, 4],
       [ 16, 20]])

as expected.

I suspect this is due to some of Emily's legacy code in my EpochArray.__init__ method.

Similar to intersection, but an EpochArray's complement within the current EpochArray.

One way to do this is to take the complement in [0, end] and then to intersect with another EpochArray.

The complement can be found efficiently by a few manipulations to starts and stops. However, if epochs are overlapping, then how should the complement be defined? Perhaps the complement of the final merge() should be returned?

If we eventually fix the determine_num_units() method, then test with

sta = SpikeTrainArray([[1],[2],[3],[4],[5]])

In other words, 5 units each with one spike. Do we get the correct result out?

If fs is None, then it can later be set; but the first time it's set, samples and time are left unaltered. So we may have a situation where time == samples, yet there is a non-unity fs. Clearly the expected behavior should be to warn the user that time is changing, and then to enforce the time = samples / fs relationship. This just makes the initial setting of fs a little more tricky, but it should still be a relatively easy fix.

See annotate() here: http://neo.readthedocs.io/en/0.4.1/core.html#example-spiketrain

I can't figure out how to access (or intercept) the __repr__ from the parent class (hmmlearn.hmm.PoissonHmm)?

See https://docs.python.org/3/library/stdtypes.html#typememoryview

Also, in neuropipes, we will often want to work with views of the curated data, I think.

I really like to return empty objects when indexing out of bounds, or similar, but in order for the iterator functionality to work, I need to raise an IndexError, which prevents me from returning an empty object.

How can I rewrite the code, or return something AND raise an exception for the iterator to work?

See http://notes.brazen.ca/python-iterators.html for more on iterators.

Currently, we have code of the form

if support is None:
    self.support = EpochArray(np.array([0, xdata[-1]]), fs=fs)
else:
    #restrict signal to only those within the support

but this is not the desired behavior. If no support is given to AnalogSignal, it should figure out which contiguous segments the signal is defined on by looking at xdata. If no xdata is given, AND no support is given, then the above is appropriate. If no xdata is given, but an explicit support is given, well, then we'll have to think very carefully about what that means, and if it should even be allowed / supported.

Many things wold be simpler (where, logical indexing, mean, std, len, dim, shape, scalar multiplication, scalar addition and subtraction, vector addition and subtraction, etc.) but other things in the interface may not make sense, like transpose, vector multiplication, squeezing, adding new dimensions, etc.

Essentially, for a wrapper, we'll have to re-implement or link MANY methods, but with inheritance we'll have to police and restrict many operations that don't make sense.

Which approach is better?

last_st = np.array([unit[-1] for unit in sampleArray]).max() raises an IndexError when a spiketrain (unit) is empty (which can easily happen).

We currently use endpt = (0, foo), but we should be able to check whether the axis is reversed, in which case we should form it differently, like endpt = (foo, 0) or something.

I successfully implemented a hack to properly center the y-label when an (x+y)-scalebar is used. However, the code is now a little clunky, and should be improved. In addition, I should implement ec=None in the fix, so that I don't draw multiple lines over each other.

See here:

def __getitem__(self, i):
        '''
        Get the item or slice :attr:`i`.
        '''
        obj = super(AnalogSignal, self).__getitem__(i)
        if isinstance(i, int):  # a single point in time across all channels
            obj = pq.Quantity(obj.magnitude, units=obj.units)
        elif isinstance(i, tuple):
            j, k = i
            if isinstance(j, int):  # extract a quantity array
                obj = pq.Quantity(obj.magnitude, units=obj.units)
            else:
                if isinstance(j, slice):
                    if j.start:
                        obj.t_start = (self.t_start +
                                       j.start * self.sampling_period)
                    if j.step:
                        obj.sampling_period *= j.step
                elif isinstance(j, np.ndarray):
                    raise NotImplementedError("Arrays not yet supported")
                    # in the general case, would need to return IrregularlySampledSignal(Array)
                else:
                    raise TypeError("%s not supported" % type(j))
                if isinstance(k, int):
                    obj = obj.reshape(-1, 1)
        elif isinstance(i, slice):
            if i.start:
                obj.t_start = self.t_start + i.start * self.sampling_period
        else:
            raise IndexError("index should be an integer, tuple or slice")
        return obj

sampleArray = np.array([np.array(st) for st in samples]) assumes that there are multiple units, and will break if not.

• Line expansion
• Plot subset of SpikeTrain
• Get labels from SpikeTrain object
• Receive argument to re-index
• Clean up docstring

sta = nel.SpikeTrainArray([[2,3,4],[3],[5],[1,2,3,4,5]]) works
sta = nel.SpikeTrainArray([[2,3,4],[3],5,[1,2,3,4,5]]) fails

We need to somehow coerce each element into the expected type. However, to know the proper thing to do, we need to know if it's a single or multi-unit object. See #46

Check when an axis is reversed, and adjust default behavior. Also build in a way to explicitly configure orientation that is axis-orientation-aware.

extend AnalogSignal to have an EpochArray support, and to support indexing, slicing, .mean(), .std(), .max(), .min(), .trim([from, to], .clip(min, max), and .smooth()

Questions: are all of these sensible for a generic AnalogSignal class, or are they more specialized to the LFP instantiation? How about the position class that Emily has created as a subclass of AnalogSignal, for example? Worth a second look, and thinking a little before we commit.

Currently only implemented in EpochArray.

Depending on where the data comes from, (esp if we read it in from .mat files) we could maybe face NaNs of some sort.

Raster should work on SpikeTrain, SpikeTrainArray, BinnedSpikeTrain, BinnedSpikeTrainArray, list, numpy array, and numpy n-d arrays.

Also make sure that even SpikeTrainArray and SpikeTrain has .n_unit properties that work as expected.

Emily sorted epoch start times in her EpochArray class. But this is inconsistent with other behavior, and returns unexpected results when you do epocharray[::-1] for example, in which case the user hopes that the order of the epochs is reversed.

However, I have not checked thoroughly why she sorted? Was this an implicit assumption in the merge() method? That needs to be fixed / improved anyway, so maybe it's not critical, but I'm removing the sorting now, and need to confirm (through writing tests, most likely) that it doesn't break anything.

Also consider adding an issorted property as with some of our other objects.

We currently do sort starting times in EpochArray, which doesn't necessarily make sense, but we do NOT sort in SpikeTrain, where it almost certainly DOES make sense.

It may be very useful to be able to do something like

analogsignal.where(foo > 5) or analogsignal[analogsignal < 3] or even analgsignal[analogsignal > analogsignal2] and so on, also supporting ==.

This means that we'll have to support indexing by boolean arrays, and we may need to check for other things, too. This goes against duck typing, so a better alternative might be to inherit from numpy array, but this warrants its own discussion. See #37 .

Also, the behavior should be specialized, but natural, for each type of object. For example, epocharray[epocharray > 100] should return only those epochs in epocharray whose (starts? stops? centers?) are greater than 100. Of course we can already do something very similar like so

epochs_greater_than_100 = EpochArray([100, infty])
epocharray[epochs_greater_than_100]

but the point is to make the API as easy, quick, elegant, and intuitive to use as possible.

Should be a method if numpy's is a method, and a property otherwise.

There are times when we need to specify an Epoch without finite bounds, like this [100, infty). We may want to use this, for example, if we want to take the intersection with another EpochArray or other object such that the new object (after the intersection or restriction) is only defined for 100 seconds and later.

I think Emily had something like this in her time_slice and time_slices methods, where None was taken to be unbounded. I think it makes sense to support None in this context, since there is no upper bound in our example above, but also to support np.infty.

If we support None, then we'll have to change that to np.infty pretty early on in the __init__ method to avoid trouble with input validation where we check that stop > start, for example, as well as other places where we compute the epoch duration, etc.

np.infty should work just fine with duration, simply returning np.infty again. But this case (with one unbounded, and also with both unbounded edges) should definitely be written up as a test case!

Easy install works locally though? Not sure what's going on...

Remember to implement .shape()

And which ones do I currently return? It is wholly possible that I am not consistent in behavior in this regard. This should be figured out and fixed before it becomes too unwieldy.

It is probably a good idea to have a .copy() method for each object as well.

Consider this example:

stdata = np.array([1, 2, 3, 4, 5, 6, 6.5, 7, 8, 10, 18, 18.5, 23, 25])
ep2 = EpochArray(np.array([[1,2],[4,6],[7,10]]), fs=1)
st = SpikeTrain(stdata, fs=1, support=ep2)
print(st.support.time)
bst = BinnedSpikeTrain(st, ds=2.8)

print(bst.bins)

Here, since ds=2.8 we have overlapping binned support for adjacent epochs. How should this be handled? Bin edges can not be assumed to be monotone, even if epochs in EpochArray are monotone and non-overlapping.

Neo uses explicit units, which can be a little messy and annoying, but also very elegant and powerful.

In contrast, so far we assume Hz for sampling rate, and seconds for time. (Probably cm for position). Should we consider adopting a more powerful approach a la Neo?

We have time and samples. Should we rather use times and samples? I'm not sure. For an AnalogSignal time is more appropriate for one of the dimensions. But for a spike train, "times" seems more appropriate.

Also add in .trim() for time, similar to .clip in numpy.

• fix iterator end
• add support for steps
• change to epoch-based support
• add syntax to index spike times directly

Currently we can slice a SpikeTrain object like so st[:5] which will return the first 5 spikes. This is cool, since we can use it to see how long we need to wait for the first five spikes (st[:5].support.duration) or we can use it to estimate the average firing rate over the first 100 spikes, for example, as opposed to some fixed length of time.

However, this type of slicing will likely be of relatively limited use in practice, and the epoch-based indexing should see much more frequent and powerful use.

Now, an alternative idea is to use the slicing syntax to slice by Epochs. For example, if an EpochArray has 100 Epochs, then st[:5] will return a SpikeTrain whose support is the first five epochs in st.support (which is an EpochArray).

If the spiketrain has only one continuous support, then st[:5] will return the entire spiketrain, just like numpy array slicing. But in this way, say we have an EpochArray corresponding to SWR events, then we can get the spiketrain during the first SWR event with st[0] and the first five SWR events with st[:5 ] for example.

A useful thing to add might then be to add a __repr__ string for how many Epochs are in the support of a SpikeTrain object.

Need to think about this a little more, but here are some ideas:

• Epoch + Epoch ==> EpochArray

• EpochArray + EpochArray ==> merge()

• EpochArray - Eppoch(Array) ==> set removal

• EpochArray * 1.1 ==> expand

• EpochArray * 0.9 ==> shrink / contract

• SpikeTrain + SpikeTrain ==> add on EpochArray supports

• SpikeTrain * 5 ==> ? expand?

• AnalogSignal + AnalogSignal ==> ? It could make sense to actually ADD analog signals as opposed to merging them; same with scaling them. What should the operators do in these cases? Is consistency more important, or readability? AnalogSignal * 5 really seems to suggest that we simply scale the signal, whereas SpikeTrain * 5 does not.

This is especially useful if there are many cells; see Ting's NN paper figure for an example

I think it is best to enforce a merge().final on ANY support object in nelpy.core

Otherwise, how should we interpret a spiketrain with overlapping epochs? How should we then bin such a spiketrain, and how should we inherit its support?

Certainly the spikes don't happen multiple times simply because they fall within multiple epochs, and the total duration of the spiketrain should not be the simple sum of all the epochs, UNLESS the support is the merge().final EpochArray.

add spiketrain.units[list] indexing --- how? requires units to be its own class; okay; but how do we then restrict the parent object's units?

add spiketrain.units.labels (defaults to enumerate + 1)

Events are similar to Epochs, but have no duration so to speak. It could also be cool to have the option to have a label list, so that each event has an optional associated label. See NEO for some pointers.

I want to be able to call ep[999] to return an empty EpochArray when it has less than 999 epochs, and this works already. However, when I do

for ep in epochArray:
    print(ep)

I enter an infinite loop. How can I fix this iterator in EpochArray, and all other similar locations?