Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve UPO storage in periodicorbits.jl #325

Merged
merged 10 commits into from
Mar 18, 2024

Conversation

JonasKoziorek
Copy link
Contributor

I noticed that while using periodicorbits algorithm with parameter roundtol that is high (eg. 8 ), it detects many duplicate UPOs.
MVE:

using ChaosTools
using PredefinedDynamicalSystems: henon
begin
    ds = henon()
    xs = LinRange(-3.0, 3.0, 60)
    ys = LinRange(-10.0, 10.0, 60)
    seeds = [SVector{2}(x,y) for x in xs for y in ys]
    o = 5
    indss, signss = lambdaperms(dimension(ds))
    λ = 0.005
    results1 = sort!(periodicorbits(ds, o, seeds, λ, indss, signss; roundtol = 8))
end

Output:

10-element Vector{SVector{2, Float64}}:
 [-1.13135448, -0.33940634]
 [-1.13135448, -0.33940633]
 [-1.13135447, -0.33940636]
 [-1.13135447, -0.33940635]
 [0.63135447, 0.18940633]
 [0.63135447, 0.18940634]
 [0.63135448, 0.18940634]
 [0.63135448, 0.18940635]
 [0.63135449, 0.18940636]
 [0.63135449, 0.18940636]

My proposed change solves this:

using ChaosTools
using PredefinedDynamicalSystems: henon
begin
    ds = henon()
    xs = LinRange(-3.0, 3.0, 60)
    ys = LinRange(-10.0, 10.0, 60)
    seeds = [SVector{2}(x,y) for x in xs for y in ys]
    o = 5
    indss, signss = lambdaperms(dimension(ds))
    λ = 0.005
    results2 = sort!(periodicorbits(ds, o, seeds, λ, indss, signss; spacetol = 1e-7))
end

Output:

julia> results2
2-element Vector{SVector{2, Float64}}:
 [-1.1313544826122608, -0.33940632522716674]
 [0.6313544679502884, 0.18940632554538597]

I have implemented FP as a binary tree. This way, new fixed points are inserted into FP only if they are spacetol away from the remaining fixed points in FP. Thus duplicates are avoided for a good choice of spacetol.

Additionally, binary tree removes the need for linear search through FP even though this doesn't seem to be a bottleneck.

Let me know if there is something to fix. This is my first PR :).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This file has confused me a bit... So Datastructures provides a binary tree. What is the DummyStructure supposed to do? WHy do we need this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typically binary tree will store value to the left child if it is less than root and to the right side if it is more. I want to store value to the left only if it at least DymmyStructure.eps away from the root. Hence duplicate periodic points that are close closer than DummyStructure.eps won't be stored. DummyStructure is needed for == and < overload so that it can be used in the binary tree from Datastructures. Perhaps there is a more elegant way to do it. I agree that it is confusing.

roundtol::Int = 4
spacetol::Real = 1e-8,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, notice that this is a breaking change right now. You have removed a keyword, so users that were using this code will be getting an error. They should be getting a warning instead.

Do roundtol = nothing
and inside the function have a warning code:

if !isnothing(roundtol)
warn("`roundtol` keyword has been removed in favor of `spacetol`")
end

by the way, why the name spacetol? This is an absolute tolerance right? perhaps abstol is better?

Copy link
Contributor Author

@JonasKoziorek JonasKoziorek Mar 16, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right. I added warning. abstol would stand for absolute tolerance. I renamed spacetol as absintol for absolute intolerance, which is what it actually is: dropping results in absintol neighborhood.

Comment on lines 3 to 4

struct DummyStructure{T}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, please add comments here explaining why this exists (copy paste the explanation you gave in the GitHub discussion)

*if* they are not already contained in `FP`.
This is done so that `FP` doesn't contain duplicate fixed points (notice
that this has nothing to do with `disttol`).
* `absintol = 1e-8`: A detected fixed point isn't stored if it is in `absintol`
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright, then rename this to abstol. it is the same as abstol in the isapprox function so it will make it easier to approach for a user.


FP = typeof(current_state(ds))[]
type = typeof(current_state(ds))
FP = RBTree{DummyStructure{type}}()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here in the FP write a couple of lines of code why you use the tree. The same explanation you gave to me.

Copy link
Member

@Datseris Datseris left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All in all good, I have some standind comments for adding comments in the source code to increase clarity for the non obvious additions.

Copy link
Member

@Datseris Datseris left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perfect. Please rename DummyStructure to NumberWithEpsRadius, which is a much more indicative name of the data structure. Then I'll merge. Would be nice to add a changelog entry and bumb the minor version in project toml. We will release this along with your other PR.

@JonasKoziorek
Copy link
Contributor Author

I renamed it as VectorWithEpsRadius since it's used with SVectors.

Should the changelog include info about both PRs? I will make another PR this week adding the Davidchack-Lai algorithm. It will internally also make use of VectorWithEpsRadius. Perhaps all three could be released together.

@Datseris Datseris merged commit 9631b33 into JuliaDynamics:main Mar 18, 2024
1 of 2 checks passed
@Datseris
Copy link
Member

we can add the changelog in your 3rd PR, it's fine.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants