JuliaDynamics / RecurrenceAnalysis.jl

Recurrence Quantification Analysis in Julia
https://juliadynamics.github.io/RecurrenceAnalysis.jl/stable/
MIT License
45 stars 12 forks source link

Extracting pieces of a AbstractRecurrenceMatrix #118

Open Lucas-Froguel opened 3 years ago

Lucas-Froguel commented 3 years ago

Hello there. I am working in a new way to calculate some RQA quantifiers and, for this purpose, I need to extract some smaller sets of the Recurrence Matrix (RM). I have implemented this in my own code, however I would like to integrate it with DynamicalSystems to use its high speed. In my code, I have defined a RecurrencePlot (RP) as a Array{Bool, 2} and, hence, extracting a smaller block reduces to RP[1:3, 1:3] (for example). Thus, I was expecting that something like RM = RecurrenceMatrix(mapp, e) micro = RM[1:3, 1:3] would return micro as the same type of RM. However, when I do this, I get: typeof(micro) >> SparseArrays.SparseMatrixCSC{Bool,Int64} And this does not allow me to use determinism(micro), for the RQA quantifiers only accept ::AbstractRecurrenceMatrix. Hence, I looked for the documentation of ::AbstractRecurrenceMatrix to see whether I could somehow get my micro to be the same type as RM, but this search was unsuccessful. Is there a way to do this or could this feature somehow be implemented?

It has occured to me that I could get the variable micro the way I want if I do the following: micro = CrossRecurrenceMatrix(mapp[1:3], mapp[1:3], e) However this implies that the whole calculation of the matrix will be made every single time and I was looking for a simpler way to do this, such as the one I mentioned or, at least, not as expensive computationally. I also realize the @windowed macro does something of the fashion, but not quite what I need.

Datseris commented 3 years ago

Hi there,

it is really simple to fix this. You can do a Pull Request to RecurrenceAnalysis.jl that implements the following method:

Base.getindex(R::RecurrenceMatrix, args...) = RecurrenceMatrix(getindex(R.data, args...))

this will transform the sub-selection into a type RecurrenceMatrix instead of SparseMatrix and will allow other quantities to be calculated from it.

Datseris commented 3 years ago

Naturally we should make the function a bit more detailed, so that it checks that your sub-selection is still a square matrix, but other than that, that's about it.

Datseris commented 3 years ago

@heliosdrm actually, we could consider the possibility of allowing recurrence matrices to store views of the data, similarly with how we did for Dataset.

heliosdrm commented 3 years ago

Agreed, we can consider it, although there are some decisions to make about it:

  1. I see no problem in creating a CrossRecurrenceMatrix out of any kind of view of an AbstractRecurrenceMatrix; but in my eyes, a RecurrenceMatrix only makes sense if the view is of another RecurrenceMatrix, and the subset of indices is the same for columns and rows. (Likewise for JointRecurrenceMatrix.) So: do we allow only CrossRecurrenceMatrix to take views, or also to RecurrenceMatrix and JointRecurrenceMatrix with those restrictions?
  2. (Cross/Joint)-recurrence matrices receive the parameter WithinRange or NeighborNumber, which sets how neighbors are defined. Which should be the parameter of view(R, rows, cols)? Without giving much thought to it, I'd say that the simplest would be just inheriting the parameter of R - although that parameter would not mean exactly the same when the data is a view.
Datseris commented 3 years ago
  1. Yeap, it should just inherit the parameter of the parent array. Although I can see how this will make some "problems" if one uses the view with downstream rqa functions, because they implicitly assume some kind of normalization that stems from the parameter type.
  2. I think we allow everything, but explicitly check during the construction of a view of a Joint/RecurrenceMatrix that indeed the indices of x and y dimensions are identical.
heliosdrm commented 3 years ago

I can see how this will make some "problems" if one uses the view with downstream rqa functions, because they implicitly assume some kind of normalization that stems from the parameter type.

The inconsistency would only happen for NeighborNumber, since the columns of the sub-matrix won't have the same number of neighbors anymore. In the case of WithinRange, neighbors will still be within the same fixed range. So, maybe we can trigger a warning if view is used with an AbstractRecurrenceMatrix{NeighborNumber}.

I think we allow everything, but explicitly check during the construction of a view of a Joint/RecurrenceMatrix that indeed the indices of x and y dimensions are identical.

And in addition: provide a cheap way to make a CrossRecurrenceMatrix from any other AbstractRecurrenceMatrix:

CrossRecurrenceMatrix(R::AbstractRecurrenceMatrix{T}) where T = CrossRecurrenceMatrix{T}(R.data)

Thus, if I want a view of a RecurrenceMatrix that wouldn't qualify as RecurrenceMatrix itself, but could be a CrossRecurrenceMatrix, I can do this without penalty:

@view CrossRecurrenceMatrix(R)[rows, columns]
Datseris commented 3 years ago

Yeap all great ideas :+1:

heliosdrm commented 3 years ago

Unfortunately, not that easy. The data inside recurrence matrices are sparse matrices, and some methods use specific functions for such matrices, like nnz, rowvals, etc. But those functions do not work with views of sparse matrices. So we cannot just make a subtype of AbstractRecurrenceMatrix containing a view instead of the original matrix, and expect it to work with the current generic methods. Also, specific methods for views would probably be much slower :disappointed:

Related discussions in Discourse: https://discourse.julialang.org/t/slow-arithmetic-on-views-of-sparse-matrices/3644 https://discourse.julialang.org/t/how-to-get-a-sub-matrix-of-a-large-sparse-matrix-array-efficiently/36741

So, maybe using getindex is the way to go.

Datseris commented 3 years ago

thanks for having a look. Alright, it's not that bad honestly. The copy operation should probably be insignificant in performance compared to actually computing RQA measures.