JuliaDynamics / Attractors.jl

Find attractors (and their basins) of dynamical systems. Perform global continuation. Study global stability (a.k.a. non-local, or resilience). Also tipping points functionality.
MIT License
32 stars 6 forks source link

Make `AttractorsViaRecurrences` more simple when finding and storing attractors #103

Closed Datseris closed 11 months ago

Datseris commented 12 months ago

Closes #102.

EDIT: Updated.

This PR brings many positive changes. It originally was just a simple "fix #102" but this requires much more work... Details:

Datseris commented 12 months ago

Wow, only 1 test has failed! That is amazing.

(well, 3 tests have failed, but 2 of them are irregular grid tests that check how many points were stored in the attractor. the tests themselves don't make much sense, but I've also changed the default storage amount here anyways)

Datseris commented 12 months ago

Extensive tests bring more failures. However, I will first wait to see if @awage agrees with this change and finds it a step forwards or not, before investing any time to fix tests.

awage commented 12 months ago

I am digging into the test results. I will try to see what is wrong with the henon example.

Datseris commented 12 months ago

Sure, but please don't push something on this branch as I am working on it. If you find the solution simply report it as a comment.

Datseris commented 12 months ago

By the way, this change is great even for us, not just for users. This new error clause I've added is very helpful. It detects very well when the algorithm has problems, either due to slowing down (like the problems @KalelR reported), or because the cells are too large and a cell may contain two attractors. I've tested this with the fast-slow systems, as well as with the Lorenz84 where the attractors intersect. Now, the algorithm will error eagerly instead of allowing the iteration to continue and the outcome to be just wrong (and hence not useful).

Note that this will affect the continuation: if the grid or the setup is bad, continuation will also error. But I think this is a great quality of life change because it informs users about bad algorithm behavior early, and this way users don't have to trouble shoot or imagine wjhy the algorithm didn't behave well!

We should focus on fixing the tests soon so that I merge this and tag this, because I am giving a workshop about Attractors.jl in 10 days.

Datseris commented 12 months ago

If the error occurs, tjhat's the error message we get:

 During the phase of locating a new attractor, found via sufficient recurrences,
we encountered a cell of a previously-found attractor. This means that two
attractors intersect in the grid, or that the precision with which we find
and store attractors is not high enough. Either decrease the grid spacing,
or increase `mx_chk_fnd_att` (or both).

Index of cell that this occured at: CartesianIndex(1, 1).
awage commented 12 months ago

Nice! It is really helpful.

I found the bug. Due to the new method the lines between 181 and 186 should be:

            current_basin = bsn_nfo.current_att_label + 1
            cleanup_visited_cells!(bsn_nfo)
            bsn_nfo.visited_cell_label += 2
            bsn_nfo.current_att_label += 2
            reset_basins_counters!(bsn_nfo)
            # We return the label corresponding to the *basin* of the attractor
            return current_basin

If you return ic_label you may have unexpected behavior .

awage commented 12 months ago

Seems that all test passing except for the irregular grids!

Datseris commented 12 months ago

@awage check out lines 281-288 of the source code:


    if next_state != current_state
        # reset counter except in lost state (the counter freezes in this case)
        if current_state == :lost
            bsn_nfo.consecutive_lost = 1
        else
            bsn_nfo.consecutive_match = 1
        end
    end

We have the comment that we "freeze" the lost counter, but in fact, we do the opposite. The counter is reset!!! What are we supposed to do? Freeze it? in this case we need to delete the bsn_nfo.consecutive_lost = 1 line!

awage commented 12 months ago

The comment is misleading. When you switch state the current_state and next_state are different. The internal counter consecutive_match is set to zero.

The exception is when you move out of the :lost state: you set consecutive_lost to one but you do not reset consecutive_match. The algorithm pick up where it left and resume its activity.

I am not satisfied with this part of the code. But I do not find anything better for now. The comment should be updated.

awage commented 12 months ago

Maybe the way is to treat the :lost state separately. I will draft a PR.

awage commented 12 months ago

see #105

awage commented 12 months ago

I have done minimal changes but it is more readable now! I think the behavior of the system is the same.

I am thinking about including a test specifically for the :lost state. What do you think?

awage commented 12 months ago

This is a simple test for the lost state, part of the attractor lies outside the defined grid. What is amazing is that the algorithm finds its way out and identifies the attractor:

    henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
    henon() = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3])
    ds = henon()
    xg = yg = range(-1.0, 1.0; length=100)
    grid = (xg, yg)
    mapper = AttractorsViaRecurrences(ds, grid)
    @show mapper([0. ,0.])
Datseris commented 11 months ago

@awage or @KalelR this PR is now finished. I've updated the top level comment of the PR with all the changes this PR does: https://github.com/JuliaDynamics/Attractors.jl/pull/103#issue-1931705212 .

Please have a read, and give a review. If you give an accepting review, I'll merge.

(I can't merge without a passing review, I've updated the repository settings so that all PRs need a passing review. It's good to ensure I don't just casually push stuff to main that I may regret later.)

awage commented 11 months ago

All good for me. I haven't checked the irregular grid stuff yet. I will do it another day. But the rest is good to merge!