juliadynamics / recurrenceanalysis.jl Goto Github PK

View Code? Open in Web Editor NEW

44.0 10.0 11.0 1.41 MB

Recurrence Quantification Analysis in Julia

Home Page: https://juliadynamics.github.io/RecurrenceAnalysis.jl/stable/

License: MIT License

Julia 97.21% R 0.69% MATLAB 0.95% Python 1.16%

recurrence-quantification-analysis nonlinear-dynamics julia rqa recurence nonlinear-timeseries-analysis hacktoberfest

recurrenceanalysis.jl's Introduction

RecurrenceAnalysis.jl

A Julia package that offers tools for computing Recurrence Plots and exploring them within the framework of Recurrence Quantification Analysis and Recurrence Network Analysis. It can be used as a standalone package, or as part of DynamicalSystems.jl.

To install it, run import Pkg; Pkg.add("RecurrenceAnalysis").

All further information is provided in the documentation, which you can either find online or build locally by running the docs/make.jl file.

recurrenceanalysis.jl's People

Contributors

Stargazers

Watchers

Forkers

tkelman felixcremer asinghvi17 ginkulv t-bltg rs7q5 battyone isgasho pagesys allomulder pucicu

recurrenceanalysis.jl's Issues

1.0.0 release

For stability reasons, a 1.0 release is necessary, so that we can start following SymVer as far as breaking changes are concerned. @heliosdrm is there anything missing or anything you would like to do before we do this?

Warning: implicit `dims=2` argument now has to be passed explicitly to specify that distances between columns should be computed

During testing with pkg > test RecurrenceAnalysis, the test threw this warning:

┌ Warning: implicit `dims=2` argument now has to be passed explicitly to specify that distances between columns should be computed
│   caller = ip:0x0
└ @ Core :-1

Do you want to join JuliaDynamics and DynamicalSystems.jl?

JuliaDynamics

Hi!

I was not aware of this package, and I can see that it is very interesting and helpful! In addition I can also see that there is a significant overlap with DynamicalSystems.jl and what you are trying to do here!

This is an invitation to join the JuliaDynamics organization, which I believe this package is a perfect fit. Unfortunately joining the organization means that you would have to share owner privilliges with user @Datseris, but I can assure you that this is nothing to worry about.

Joining DynamicalSystems.jl

This is also a second invitation to join DynamicalSystems.jl.
The way DynamicalSystems.jl works is that it re-exports existing packages (currently DynamicalSystemsBase, ChaosTools and in the near future GAIO) and host a unified documentation for them.

I have seen a significant amount of overlap already with what you have here. For example we also have:

timeseries embedding (supporting multi-time and multi-series)
renyi entropies
mutual information (I need to compare ours with yours, because your method may be better)
Methods to estimate optimal embedding parameters (but not the False nearest neighbors one!)

Should you wish to join DynamicalSystems.jl there are some effort that needs to be put into this (from both parties). We need to eradicate the duplicated functionality, e.g. 1 embedding method should be used, which ever is best. The same goes for all other overlapping. In addition, documentation strings should be enriched with descriptions and citations of what is going on.

Why?

DynamicalSystems.jl is an award-winning software. Besides this, it aims to be a trusty companion to both students and scientists working on nonlinear dynamics and chaos.

DynamicalSystems.jl will benefit significantly from having your package there, because I believe it is not only a job well done, but also very useful!

RecurrenceAnalysis.jl will also benefit from joining DynamicalSystems.jl / JuliaDynamics : it will get much more exposure and more testing. This will help the package mature into something robust and trustworthy.

In addition, for users, it is much more helpful to have things in a unified and "respected" (i.e. with real people and real packages) organization (in this case JuliaDynamics) that has active members and inspires trust for the users. It also helps finding packages. For example even though I am very active in nonlinear dynamics in Julia, I had no idea about your package until today.

Run femtocleaner

[DOC] Enrich RQA docstrings with comments for usefulness / understanding

I suggest that little by little we start enriching the docstrings of all RQA functions by copy/pasting the relevant text from Norbert's book. This will really make the RQA page a "library" page (which is my life-long dream with DynamicalSystems.jl).

For example, in the determinism we can copy paste:

Systems possessing deterministic (rule-obeying) dynamics are characterized by diagonal lines indicating repeating recurrences within a state. For periodic signals the diagonal lines are long. For chaotic signals the diagonal lines are short. For stochastic signals the diagonal lines are absent, save for chance recurrences forming very short lines. DET can be interpreted as the predictability of the system more so for periodic behaviors than chaotic processes. But it must be realized that DET does not have the same strict meaning of the determinism of a process.

And so on and so forth.

Specialized type for Boolean sparse arrays

Recurrence matrices are wrappers of SparseMatrixCSC{Bool,Int}, which are useful to make some operations. However it is inefficient to create them (cf. #76), and also take more memory than what is needed (all nonzero values are true, so they don't add any information at all).

It would be useful to have a special type of sparse matrices for these cases, which only store row indices and the last index of each column - but not the values, with optimized methods. This might come in handy even for other applications, not only for recurrence analysis, so in my opinion it is worth to be in a separate package (outside JuliaDynamics).

I feel like investigating if such a package exists, and otherwise draft it. When it is working, we could use it to improve performance in this package.

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

Conversion of inputs in distancematrix

TL;DR: I think we should revert that change introduced in a1d0404

The code of distancematrix was changed in a1d0404 and 86d7356 from:

if sx[2] > _maxdimension
    return _distancematrix(Matrix(x), Matrix(y), metric)
else
    return _distancematrix(Dataset(x), Dataset(y), metric)
end

to:

if sx[2] < MAXDIM # convert to Dataset here, it is significantly faster
    return _distancematrix(Dataset(x), Dataset(y), metric)
else
    return _distancematrix(x, y, metric)
end

Let aside the name of the constant (previously _maxdimension, changed to MAXDIM) and the sign of the comparison, the main difference is that when the dimension of the embedding space is greater equal than MAXDIM, x and y were converted into matrices to gain performance (cf. https://gist.github.com/Datseris/3ca61b32dea2033aeae5c4dfe3e0b78c), and now they are left as they were.

But in many (perhaps most) cases, if x or y been made by embed or reconstruct, they won't be matrices, but Datasets, so even although their dimension is greater than MAXDIM, the two branches of the condition will run the same method.

Or worse (although that's a far-fetched rare case): if one of them is a Matrix and the other a Dataset, and their dimension is greater than MAXDIM, no appropriate method of _distancematrix will be found and the function will fail.

Calculation of RQA parameters for "white vertical structures"

Create the dedicated functions, and extend rqa to support the new related keys.

Specific code for `recurrencematrix`

recurrencematrix calls crossrecurrencematrix with the same x and y arguments, thus duplicating the number of comparisons and the size of the nonzero indices of the resulting matrix.

If we have different RecurrenceMatrix and CrossRecurrenceMatrix types, the size and number of calculations for the former may be reduced. Type-specific methods for RQA and for the creation of the plottable matrix can take into account the symmetry of recurrence matrices. (And also for joint recurrence matrices.)

deprecation of `crossrecurrencematrix`

recurrencematrix, crossrecurrencematrix and jointrecurrencematrix are deprecated in favor of the type constructors RecurrenceMatrix, etc. But crossrecurrencematrix is still used internally to create the data field (i.e. the actual matrix).

This is a bit confusing. Should we change the name of that internal function?

On the plot of a recurrence matrix

As I am preparing some material for the upcoming Julia tutorial I have noticed some difficulties when actually making the plots of recurrence matrices.

In my tutorial I will be use the Lorenz system as a case study. To explain a recurrence matrix I will be comparing a parameter where the lorenz system is periodic and a parameter where it is chaotic. Here is an example that shows both a recurrence plot and a trajectory for the two cases:

lor = Systems.lorenz()
figure(figsize = (12,8))

for (i, ρ) in enumerate((160.0, 28.0))
    set_parameter!(lor, 2, ρ)
    t, dt = i == 1 ? (50.0, 0.01) : (50.0, 0.01)
    tr = trajectory(lor, t; dt = dt, Ttr = 2000.0)
    
    x = tr[:, 2]
    subplot2grid((3,2), (0, i-1))
    plot((0:dt:t)[1:2000], x[1:2000], "k")
    title("ρ = $ρ, " * (ρ != 160.0 ? "not periodic" : "periodic"))
    
    subplot2grid((3,2), (1, i-1), rowspan = 2)
    R = RecurrenceMatrix(tr, 5.0)
    imshow(recurrenceplot(R), cmap = "binary_r",  extent = (0, t, 0, t));
    xlabel("X"); ylabel("Z"); 
end

which gives

What I noticed is that the bottom left plot is extremely misleading. For this parameter value ρ=160 the Lorenz system is periodic. In fact examining the Poincare surface of section reveals that the system has "period 3" in the Poincare map. In addition, by looking at the timeseries plot one can see that the period of the trajectory is around 1 unit of time. This means that there should be a diagonal line spanning the recurrence plot, starting from time 1 (notice that the plots above are in correct time units not index units).

This is not what I see. I investigated further, by comparing the recurrenceplot function with the function scatterdata, which gives exactly the same plot data but in the form of a scatter plot. Here is the result:

lor = Systems.lorenz()
set_parameter!(lor, 2, 160.0)
tr = trajectory(lor, 50.0; dt = 0.01, Ttr = 2000.0)

R = RecurrenceMatrix(tr, 5)
L = size(R)[1]

figure(figsize = (10, 5))
subplot(121)
plot(RecurrenceAnalysis.scatterdata(R)..., linestyle = "None", marker = "o", ms = 0.1, color = "black")
gca().set_aspect("equal"); xlim(1, L); ylim(1, L)
subplot(122)
imshow(recurrenceplot(R), cmap = "binary_r",  extent = (0, L, 0, L));

which gives

In the above the axis are in index units, divide by 100 to get real time. Anyway, there is a massive visual difference between the two plots. My first concern with this issue is that it should be mentioned in the docs that imshow-like plots using recurrenceplot can be misleading. We should also export scatterdata and probably find a better name like scatterrecurrenceplot (becomes too long).

I now have a second, conceptual concern. In the above plot I will now zoom in the bottom left corner of both plots, to come closer to time=1 (i.e. index=100) which is the period of the orbit:

Now you can see that the imshow plot actually has the same information as the scatter plot. (I never doubted that, but a new user can be greatly misled).

Anyway, my second conceptual concern is: why are the off diagonals not continuous? Why do they have white interruptions? Since the orbit is periodic exact diagonal should exist right? And they should be uninterrupted. I am asking so I can give a proper explanation of the event during the tutorial.

EDIT: I Think I have resolved the second concern. Doing

figure()
plot3D(columns(tr)..., lw = 0.1)

shows that the orbit thickens when the Z variable reaches the maximum value. Must be some kind of weird quasi-periodicity.

Then the width may actually exceed the threshold of distance of 5 I had in the original computation. Using ε = 10 instead resolves the problem.

SECOND EDIT: The above case is probably a higher multiplicity of period 3, like when the periodicity has just splitted and you require a lot of time to return to original state.

`rqa` functions should also take in SparseArrays

This makes it easier to calculate metrics on segments of a recurrence matrix.

names of RQA functions

After some discussion in #24 I propose the following changes in the names of the functions for RQA. It means deprecating three of the older names, and adding a couple of parameters that I have not actually seen reported in papers, like the entropy of vertical lines or the longest recurrence time. But the advantage is that the set of functions and their names follow a more systematic pattern, easier to remember. This is what pyunicorn and PyRQA also do, so it facilitates understanding for users of both Python and Julia.

(Each function name is accompanied by the key of the corresponding dictionary item returned by rqa in square brackets.)

These conventional names remain without changes:

recurrencerate ["RR"]
determinism ["DET"]
divergence ["DIV"]
laminarity ["LAM"]
trappingtime ["TT"]
trend ["TND"]

In addition, diagonal lines, vertical lines and recurrence times will have functions to calculate their mean and maximum value, and their (Shannon) entropy, with the following common pattern (this may be nicely meta-programmed):

For diagonal lines:

dl_average ["L"] (deprecates avgdiag)
dl_max ["Lmax"] (deprecates maxdiag)
dl_entropy ["ENT"] (deprecates rqaentropy)

For vertical lines:

vl_average (equivalent to trappingtime, which will remain)
vl_max ["Vmax"] (deprecates maxdiag)
vl_entropy ["VENT"]

For recurrence times (all new)

rt_average ["MRT"] (with alias meanrecurrencetime", as reported in some papers)
rt_max ["RTmax"]
rt_entropy ["RTENT"]

Also add nmprt (with key "NMPRT"] for the number of the most probable recurrence time.

There are some keys that vary across sources. Namely, Norbert's paper (2007) uses "ENTR" and "TREND", where the RQA book of 2015 uses "ENT" and "TND". I have followed the latter, which are shorter, but we are now on time to change it.

(I'm also closing #30, which is redundant with this)

Classical RQA measures for periodic and chaotic orbits

Here are two parameters for the Lorenz system:

One is clearly periodic, the other is clearly chaotic. I now want to demonstrate the use of the classical RQA measures in these two cases:

lor = Systems.lorenz()
for (i, ρ) in enumerate((69.75, 28.0))
    println("\nρ = $ρ ")
    set_parameter!(lor, 2, ρ)
    ε = 2.0
    t, dt = 20.0, 0.005
    tr = trajectory(lor, t; dt = dt, Ttr = 200.0)
    R = RecurrenceMatrix(tr, ε)

    println(determinism(R))
    println(dl_max(R))
    println(dl_average(R))
    println(dl_entropy(R))
end


ρ = 69.75 
0.9996688958227877
2598
83.35871404399323
4.626554662807655

ρ = 28.0 
0.9997196293101414
4000
74.48904006046862
4.938120518819617

Here are my questions:

Almost all parameters are identical between the measures. How can this be / How can I understand this? I expected big differences at least in the entropies.
The maximum diagonal for the_chaotic_ case is higher than the regular case. Conceptually I cannot imagine this being true. Why does it happen?
Can I actually use the classical RQA measures to point out differences between periodic and (simple) chaotic trajectories, or should I go about a different approach in the video tutorial?

Don't make a copy for equal length data in JointRecurrenceMatrix

Coming from #31

[Doc] Theiler = 0

Does the value of 0 for the theiler window exclude the LOI or it includes it?

Functions for analysis of (cross-)recurrence networks

This includes the calculation of:

Adjacency matrix (based on recurrencematrix)
Edge density (based on recurrencerate)
Global clustering coefficient or network transitivity
Average path length
Global efficiency

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

Extension of RQA and CRQA to trajectories

@heliosdrm Just wanted to bring to your attention this extension to trajectories (in the DynamicalSystems sense) i.e. multidimensional timeseries of RecurrenceAnalysis. (https://www.frontiersin.org/articles/10.3389/fpsyg.2016.01835/full)

The illustration with Lorenz63 demonstrates the usefulness of the trajectory approach versus the timeseries approach.

Deprecation warnings for Julia 0.6

See https://travis-ci.org/heliosdrm/RecurrenceAnalysis.jl/jobs/219044031
Mostly related to element-wise operations with arrays.

Plot recipes

Any plans for Plots.jl? Would be happy to help.

Citation of RecurrenceAnalysis

What is the best way to cite RecurrenceAnalysis.jl?
I have seen the Citation of DynamicalSystems.jl, but is this also applicable for this package?
Could you provide a CITATION.bib file for RecurrenceAnalysis.jl?

`lmin` in `recurrencestructures`

I acknowledge a bit of inconsistency here: the default lmin in RQA functions is 2. But when I wrote recurrencestructures I set lmin=1 there. Why?

For the diagonal and vertical histograms, lmin=n makes zero the values of the (n-1) first bins. The rest is the same. Thus, I initially thought that the default lmin=2 was not very useful, it just removed some information.

However, the story is different when recurrence times are also counted: ruling out lines shorter than lmin means that sets of recurrence times between isolated points are merged in a smaller number of longer recurrence times. So for recurrence times changing lmin means more than just "deleting" or ignoring some part of the output.

All in all, I agree that a consistent lmin=2 everywhere is preferable.

An alternative would be to give always the "full" diagonal and vertical line histograms, and just use lmin for the calculation of recurrence times. The disadvantage is that this means calculating the vertical histograms twice.

I'm debating myself, don't know what is the best.

Performance regression due to new code of `trend`

Part of the big improvement in speed that I saw in v1.0.0-alpha has been lost with the changes in trend (instead of a gain in time > 50% it is around 40%) . Changes are necessary because the calculations were wrong, but I want to look for a smarter, faster way of calculating it.

Error tagging new release

The URL of this package does not match that stored in METADATA.jl.
cc: @Datseris

divergence in perpendicular recurrence plot : estimating lyapunov exponents

This paper : https://www.sciencedirect.com/science/article/abs/pii/S0375960199007513?via%3Dihub

discusses of a way to compute the maximum Lyapunov exponent from an RP using what the authors call "divergence in perpendicular recurrence plot ".

Seems to be useful and nicely ties with the other method we have to estimate λ from timeseries in DynamicalSystems.jl

Documentation structure

Hey @heliosdrm ,

I will be writing the documentation page of RecurrenceAnalysis this week, and finalize #9 ! But since I am not super familiar with all the methods, I was wondering whether you could suggest a structure for the doc pages.

I intend to use the material of the README file and simply add the documentation strings as well as runnable examples. Can you suggest some thematic structure, as in which pages will contain what, etc.? For example ChaosTools is structured like this: https://juliadynamics.github.io/DynamicalSystems.jl/latest/chaos/overview/

Document return values of empty histograms

Continuing from #41 , will do it soon.

Compat bounds

In this moment all compat bounds in Project.toml are lower bounds. But the recommended policy for the General Registry is to include also upper bounds.

If we want to follow that recommendation, we should check what are the compatible lower and upper bounds of dependencies for the next release. According to the last successful CI check, the following package versions are compatible (if we consider the PR #84):

DelayEmbeddings 1.2.0
Distances 0.8.2
StaticArrays 0.12.1
UnicodePlots 1.1.0

The corresponding minor versions are safe upper bounds. Moreover, Distances 0.7 will not compatible after merging #84. So, the line for that package should be `Distances = "0.8".

It remains to be checked what lower bounds might be given for the other packages, if we want to provide a wider range of compatible versions.

(This also applies to the Julia versions: it works for any 1.x version, but I've not checked if it would still be compatlble with 0.7 - if we still want it to be compatible with it.)

Clean up source a bit

It is possible to clean up and re-organize the source code a bit more for readability and consiseness . A lot of development was put into releasing 1.0 but since then I didn't really put effort to "clean things up".

[FR] Estimations of `ε` [$100]

Would be great if we had a function with various methods to estimate ε , the threshold of a recurrence plot.

This could work similarly with estimate_delay which gives delay times for timeseries.

Different methods are described in the book from N. Marwan and citations within, see chapter 1 (we cite the book in most docs). Some of these methods are actually very easy and straightforward and would be a great issue for newcomer contributor.

There is a $100 open bounty on this issue. Add to the bounty at Bountysource.

nomenclature of entropy functions

Following up #25, particularly the issue about the entropy measures derived from recurrence plots.

In summary, there are three entropy-related parameters, which are the Shannon entropies of (1) the "black" diagonal lines (the usual RQA entropy, which is now calculated by rqaentropy), (2) the black vertical lines (Leonardi's paper to be published), and (3) the white vertical lines (recurrence time entropy).

Besides, recurrence plots can be used to estimate the Kolmogorov entropy, but that is another story I think.

I don't think that eventually removing rqaentropy due to the wealth of entropy-related parameters is a good idea. It is one of the parameters that I have seen reported most frequently, after the recurrence rate and determinism, so its calculation should be in the same level. Yes, it can be calculated from the histogram with some other generic entropy function, but all RQA parameters can eventually be calculated from the histogram and some simple function (actually simpler than the entropy), and moreover there are options like the minimum length that should be taken into account.

pyunicorn has different entropy parameters corresponding to the three that Norbert mentioned, cf. http://www.pik-potsdam.de/~donges/pyunicorn/api/timeseries/recurrence_plot.html : diag_entropy, vert_entropy, and white_vert_entropy. Those names are clear and meaningful, but not really the style of Julia functions.

@pucicu: any suggestion for the names?

Kolmogorov Sinai entropy from a RP

Chapter two of the book N. Marwan & C.L. Webber, "Recurrence Quantification Analysis. Theory and Best Practices" discusses how the Kolmogorov entropy can be connected to and computed from Recurrence matrices. See 2.2 and 2.3.3 .

A method for estimating KS entropy has long been requested from ChaosTools, see JuliaDynamics/ChaosTools.jl#6 .

Would be cool to implement the algorithm in this package!

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

Error tagging new release

The tag name "v1.0.0-alpha" is not of the appropriate SemVer form (vX.Y.Z).
cc: @Datseris

tau_recurrence ?

The function

function tau_recurrence(x::AbstractMatrix{Bool})
    n = minimum(size(x))
    [count(!iszero, diag(x,d))/(n-d) for d in (0:n-1)]
end

is in the source, but it is not exported or used anywhere else.

What's wrong with Windows x86?

Appveyor fails for the x86 versions. No idea of what's happening there.

Return Type of rqa()

Why is the return type of rqa a Dict?
With a NamedTuple you still have the clear connection between the values and the name of the metric.
Furthermore, the return values are ordered, and it is a little bit cheaper because it needs less allocations.
See https://discourse.julialang.org/t/namedtuples-vs-dict/13119
I made some tests with NamedTuples, with the code on my branch.
For a 100 x 100 RecurrenceMatrix the number of allocations is reduced from 257 to 209 and the running time is reduced from 112,5 microseconds to 107,8 microseconds.
For me this is important, because I want to apply rqa on a large number of pixels (at least a million) with around 100 timesteps.

Ambiguous definition of `Dataset`

When I run some examples of recent issues, where I create Datasets from DynamicalSystemsBase (e.g. a Lorenz system) and then try to analyze them with RecurrenceMatrix etc., I get an error because the type of Dataset that the functions of this package recognize is RecurrenceAnalysis.Dataset, but they are DynamicalSystemsBase.Dataset.

I don't know if this is something that is wrong in my configuration...

Denominator of recurrence rate

Now it is the size of the entire matrix. This is the simplest choice, and it is also coherent with the formula given in Norbert's 2007 paper, which is the most popular reference for RQA (Phys. Reports 438:237-329).

However, when a Theiler window is applied, it might also make sense to reduce the denominator, since there are less points that may be detected as recurrences. Actually the formula in the first chapter of the RQA book from 2015 has N(N-1) in the denominator, accounting for the exclusion of the N points of the LOI.

Matlab's CRP Toolbox does adjust the denominator according to the Theiler window. Other packages that I know don't, e.g. crqa or PyRQA. (I'm not counting pyunicorn, because that package just doesn't have the option of setting a Theiler window).

Specialized method for `_recurrence_matrix(::Vector, ::Vector)`

Currently all recurrence_matrix calls transform data to Datasets in all cases. This is always faster for Matrices but I do not believe that this does not induce a performance regression for vectors.

A method should be added to the low level function _recurrence_matrix that does not do this transformation and also does not care about the metric and instead uses abs(xi - yi)

(notice: to get the gist of this issue is straightforward, but requires you to have a look at the source file: https://github.com/JuliaDynamics/RecurrenceAnalysis.jl/blob/master/src/matrices.jl

Metric parameter

In #9 @Datseris suggested:

[Optional/Suggestion] : Rework the handling of the metrics: instead of asking for a string ask for the Metric instance. E.g. "euclidean" => Euclidean(). This is important for education as it teaches the Julia way. But this is just my humble opinion.

If a funcion asks for a Metric instance, the user must be able to create it (e.g. calling Euclidean()). But this means that RecurrenceAnalysis.jl should export the Metric types, although they may have no other use than creating the parameters needed by the functions of this module. Is it not tidier to keep the Metric types unexported, and let the functions use them internally?

That is the reason why functions like distancematrix, etc. asked for strings, and a Dict{String,Metric} was used internally. If strings are bad style in Julia, perhaps symbols would be better (e.g. :euclidean, :max, etc.)?

In the notebook you also noted a risk of type instability with the current approach based on strings. Type instability pitfalls are still a bit obscure for me, so I'm probably missing something, but I don't understand why it may happen in this case: the variables that receive the strings are not used for the Metric objects in any part of the code. The string is just passed as a key of a dictionary of metrics, and the corresponding item of the dictionary is copied in a different variable.

Fail tests for Julia 0.7

Check https://pkg.julialang.org/logs/RecurrenceAnalysis_0.7.log

[DOC] Example (with plots) of using the windowed functionality

I think it will be easy to create an example where a big recurrence matrix is plotted and then "subsections" of it are created using the @windowed functionality. Something very similar with how it is displayed in Norbert's papers, but not as a static image: as a Julia code.

If someone provides an example that using @windowed makes sense (scientifically I mean) I will take care of the plotting part.

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

Parameters of "empty" histograms

This may happen more frequently in vertical structures, if there are no laminar states (no vertical lines, specially if a minimum line length is defined). But it might happen also with some signals for diagonal lines and recurrence times.

In such case the calculation of various parameters (averages, entropies, etc.) gives NaN. But is that the "right" solution? The condition is easy to check, so we may set them to zero in those cases if it's preferable. (Perhaps with a warning?)

Error tagging new release

The tag name "v0.3.0-beta" is not of the appropriate SemVer form (vX.Y.Z).
cc: @heliosdrm

RecurrenceMatrix for irregular time series

Ozken et al describe a possibility to use recurrence plots also for irregularly sampled data.
Is this in the scope of this package?

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

Dedicated types for Recurrence matrices

as discussed in #20 having dedicated types RecurrenceMatrix, CrossRecurrenceMatrix etc, all of which are subtypes of AbstractRecurrenceMatrix has a lot of benefits. Let's discuss the implementation here.

What type parameters should be there, besides the ones of the data field?
Any other fields besides the data field?

To-dos:

Create types
Ensure recurrencematrix and co return the types
Enforce that the user is passing a recurrence matrix to e.g. laminarity, because the dispatch happens on ::AbstractRecurrenceMatrix (so far there is no guarantee that x is not a unfitting matrix)
Pretty printing (which I super love). So when the user sees x they will get info that it is a recurrence matrix, etc. This is really cool for front end users!

Joining DynamicalSystems.jl !

Following from here: https://gist.github.com/Datseris/3ca61b32dea2033aeae5c4dfe3e0b78c and from the discussion in #8

Here is what I believe are the best steps forward to make this really nice package an "official" part of JuliaDynamics:

Once these are done, and if you agree of course, we make RA an official part of DS by

Make sure the source code is clear. I will read the source and wherever there is something unclear I'll either clarify myself or ask for your help in claryfying.
re-export RA from DS.
adding it to its documentation page as a set of dedicated pages (much like ChaosTools). For this to happen we will also have to write a proper documentation (i.e. separate / expand the current README to documenter acceptable files that expand the docstrings. I will take care of that).
Change metadata links and tag a new release!!!

You might think number 1 is irrelevant or useless or unecessary but the philoshopy of JuliaDynamics is that the source code is clear and understandable.

Error: function recurrence_matrix does not accept keyword arguments

This

using DynamicalSystems, RecurrenceAnalysis

ro = Systems.roessler(ones(3), a=0.15, b=0.20, c=10.0)
N = 2000; dt = 0.05
tr = trajectory(ro, N*dt; dt = dt, Ttr = 10.0)

R = RecurrenceMatrix(tr, 5.0; metric = "manhattan", scale=maximum)

gives

function recurrence_matrix does not accept keyword arguments

Stacktrace:
 [1] kwfunc(::Any) at .\boot.jl:332
 [2] #RecurrenceMatrix#1(::String, ::Base.Iterators.Pairs{Symbol,typeof(maximum),Tuple{Symbol},NamedTuple{(:scale,),Tuple{typeof(maximum)}}}, ::Type{RecurrenceMatrix}, ::Dataset{3,Float64}, ::Float64) at C:\Users\Yuxi\.julia\dev\RecurrenceAnalysis\src\matrices.jl:173
 [3] (::Core.var"#kw#Type")(::NamedTuple{(:metric, :scale),Tuple{String,typeof(maximum)}}, ::Type{RecurrenceMatrix}, ::Dataset{3,Float64}, ::Float64) at .\none:0
 [4] top-level scope at In[32]:4

Performance of computation of recurrence matrices

Initially asked by @asinghvi17 , is it possible to parallelize or put on the GPU the calculation of Recurrence matrices? I have no knowledge on the GPU, so I will only talk about parallelization. The core function that does the number crunching is:

function _recurrence_matrix(xx::Dataset, yy::Dataset, ε, metric::Metric)
    x = xx.data
    y = yy.data
    rowvals = Vector{Int}()
    colvals = Vector{Int}()
    for j in 1:length(y)
        nzcol = 0
        for i in 1:length(x)
            @inbounds if evaluate(metric, x[i], y[j]) ≤ ε
                push!(rowvals, i)
                nzcol += 1
            end
        end
        append!(colvals, fill(j, (nzcol,)))
    end
    nzvals = fill(true, (length(rowvals),))
    return sparse(rowvals, colvals, nzvals, length(x), length(y))
end

You can see that the algorithm operates via push! instead of pre-allocating. It may be a tiny bit trick, but I don't think it should be impossible to parallelize this.

@heliosdrm do you have any take on this?

Improve performance with `fixedrate=true`

In the calculation of recurrence matrices with fixed recurrence rate there is the issue commented in #46 (indirectly calling distancematrix, which is expensive):

RecurrenceAnalysis.jl/src/matrices.jl

Lines 220 to 222 in 2961a32

    
           if fixedrate 
        
               sfun = (m) -> quantile(m[:], ε) 
        
               return recurrence_matrix(x, y, 1; scale=sfun, fixedrate=false, metric=metric)

A function that gives the same result as quantile(distancematrix(x,y), ε) without calling distancematrix would help.

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

	if fixedrate
	sfun = (m) -> quantile(m[:], ε)
	return recurrence_matrix(x, y, 1; scale=sfun, fixedrate=false, metric=metric)

juliadynamics / recurrenceanalysis.jl Goto Github PK

recurrenceanalysis.jl's Introduction

RecurrenceAnalysis.jl

recurrenceanalysis.jl's People

Contributors

Stargazers

Watchers

Forkers

recurrenceanalysis.jl's Issues

JuliaDynamics

Joining DynamicalSystems.jl

Why?

Recommend Projects

Recommend Topics

Recommend Org