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

Algorithms to find basin boundary #65

Open
awage opened this issue Apr 25, 2023 · 3 comments
Open

Algorithms to find basin boundary #65

awage opened this issue Apr 25, 2023 · 3 comments
Labels
basin boundaries Related with the boundaries of basins of attraction

Comments

@awage
Copy link
Contributor

awage commented Apr 25, 2023

There is no automated algorithm in DynamicalSystems.jl to explore the boundary between two basins. Here are some ideas from simple to difficult:

  1. Compute the basins on a grid and with some filtering look for cells that have a neighbor from a different basin. This is computationally expensive especially in dimension larger than 3.
  2. Find the attractors with some initial random sampling. With the bisection method find a starting point on the boundary between the basin. Then launch a random search algorithm that explores the boundary. The idea would be to sample a small box of the phase space. If the box is on the boundary then we accept the point and if not it is rejected with some probability. There are method from the Monte Carlo Markov Chain algorithm family that may fit.
  3. Find the saddle on the boundary and iterates backward the dynamical system such that we can explore the stable manifold of the saddle.

There is an enormous amount of works in stochastic search (MCMC, Metropolis-Hasting, hit-and-run ...). And there very nice Julia packages that covers the subject.

I post here a naive attempt with the hit and run algorithm for the Duffing oscillator.

Captura de pantalla_2023-04-25_14-03-43

Smooth boundaries are harder to find.
Captura de pantalla_2023-04-25_14-10-23

@Datseris Datseris added the basin boundaries Related with the boundaries of basins of attraction label Apr 25, 2023
@Datseris
Copy link
Member

What's the code for the second figure? And which of the three options does it implement? Option (2) ?

@awage
Copy link
Contributor Author

awage commented Apr 25, 2023

I post the code here, the option (2) with the hit and run algorithm. I am not sure I implemented it correctly. It is a first attempt. It is not optimal since the purpose of the code was to evaluate the basin entropy of the boundary. But it is the same principle for the exploration of the boundary:

using Plots
using OrdinaryDiffEq
using Attractors
using Statistics


@inline @inbounds function duffing(u, p, t)
    d = p[1]; F = p[2]; omega = p[3]
    du1 = u[2]
    du2 = -d*u[2] + u[1] - u[1]^3 + F*sin(omega*t)
    return SVector{2}(du1, du2)
end

function get_basins_nfo(F, ω, grid)
    ds = StroboscopicMap(2π/ω, duffing, rand(2), [0.15, F, ω], diffeq = (;reltol = 1e-8, alg = Vern9()))
    mapper = AttractorsViaRecurrences(ds, grid; sparse = true, 
        store_once_per_cell = true,
        mx_chk_loc_att = 10000, mx_chk_safety = Int(1e7), show_progress = true)
    return mapper
end

function multi_threaded_basins(mapper, u0s)
	nthreads = Threads.nthreads()
	map_v = [deepcopy(mapper) for i in 1:nthreads]
	bas = zeros(Int16, length(u0s))
	Threads.@threads for j in 1:length(u0s)
		bas[j] = map_v[Threads.threadid()](u0s[j])
	end
	return bas
end

function hit_and_run_random_walk(F, ω, mapper, Nw = 1000)
	# Dummy grid for computations
	xg = yg = range(-3,3,length=10)
	Mx = maximum(xg)
	mx = minimum(xg)
	My = maximum(yg)
	my = minimum(yg)
	# Generate N random samples in a radius ϵ centered at a random point u1 in [xg, yg]
	# First uniform sampler:
	ε = 0.1; N = 100; #Na = length(unique(basins))
	s_v = zeros(Nw)
	# initial condition.
	u1 = [mx + (Mx - mx)*rand(), my + (My - my)*rand()]
	uv = Vector{typeof(u1)}()
	for k in 1:Nw
		uc = _hit_and_run(4*ε, u1, mx, Mx, my, My)
		s  = _get_ball_entropy(ε, N, uc, mapper)
		u1 = _rejection_filter(s, uc, u1)
                if s > 0
                    push!(uv, u1)
                end
		s_v[k] = s
	end
	return s_v, uv
end

function _rejection_filter(s, uc, u1)
	# Rejection filter
	if s == 0.
		# if we are in the interior, move to uc with probability p (arbitrary for now)
		return (rand() < 0.2 ? uc : u1)
	else
		# accept candidate if s > 0
		return uc
	end
end


function _hit_and_run(λ, u1, mx, Mx, my, My)
	reject = true;
	while reject
		# choice of a random direction on the circle:
		θ = rand()*2π
		rd = [cos(θ), sin(θ)]
		uc = u1 .+ rd*λ
		# check bounds
		if (mx < uc[1] < Mx) && (my < uc[2] < My)
			reject = false
			return uc
		end
	end
end

function _get_ball_entropy(ε, N, u1, mapper)
	u2s = [u1 + ε*[rand()-0.5, rand()-0.5] for k in 1:N]
	bsn = multi_threaded_basins(mapper, u2s)
	@show p = _box_entropy(bsn)
	return p
end


function _box_entropy(box_values)
    h = 0.
    for (k,v) in enumerate(unique(box_values))
        p = count( x -> (x == v), box_values)/length(box_values)
        h += p*log(1/p)
    end
    return h
end

# F = 0.128; ω = 1.106
F = 0.1; ω = 0.2
xg = yg = range(-3,3,length = 150)
mapper = get_basins_nfo(F, ω, (xg, yg));

@time b1, v = hit_and_run_random_walk(F, ω, mapper, 1000)

v = Dataset(v)
b,a = basins_of_attraction(mapper, (xg,yg))
heatmap(xg,yg, b')
plot!(v[:,1],v[:,2], seriestype = :scatter)

@Datseris
Copy link
Member

im guessing we could actually make a paper if we come up with an algorithm to find the basin. However, isn't this also what the saddle straddle does in #118 ? Or it is not generic enough to apply to "any" basin boundary?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
basin boundaries Related with the boundaries of basins of attraction
Projects
None yet
Development

No branches or pull requests

2 participants