DynamicalSystemsBase.jl icon indicating copy to clipboard operation
DynamicalSystemsBase.jl copied to clipboard

Is there a specific reason why Float64 was enforced in PoincareMap?

Open hibiki-kato opened this issue 3 months ago • 1 comments

BigFloat support in PoincareMap

Right now, PoincareMap is hard-wired to Float64 in several fields:

mutable struct PoincareMap{I<:ContinuousTimeDynamicalSystem, F, P, R, V} <: DiscreteTimeDynamicalSystem
	ds::I
	plane_distance::F
 	planecrossing::P
	Tmax::Float64 # here
	rootkw::R
	state_on_plane::V
    tcross::Float64 # here
	t::Int
    # These two fields are for setting the state of the pmap from the plane
    # (i.e., given a D-1 dimensional state, create the full D-dimensional state)
    dummy::Vector{Float64} # here
    diffidxs::Vector{Int}
    state_initial::V
end

and initializer is

function PoincareMap(
		ds::DS, plane;
        Tmax = 1e3,
	    direction = -1, u0 = nothing,
	    rootkw = (xrtol = 1e-6, atol = 1e-8)
	) where {DS<:ContinuousTimeDynamicalSystem}

    reinit!(ds, u0)
    D = dimension(ds)
	check_hyperplane_match(plane, D)
	planecrossing = PlaneCrossing(plane, direction > 0)
	plane_distance = (t) -> planecrossing(ds(t))
    v = recursivecopy(current_state(ds))
    dummy = zeros(D)
    diffidxs = _indices_on_poincare_plane(plane, D)
	pmap = PoincareMap(
        ds, plane_distance, planecrossing, Tmax, rootkw,
        v, current_time(ds), 0, dummy, diffidxs, recursivecopy(v), # problem occurs here!
    )
    step!(pmap) # this ensures that the state is on the hyperplane
    pmap.state_initial = recursivecopy(current_state(pmap))
    pmap.t = 0 # first step is 0
    return pmap
end

However, current_time(ds) adheres to the type of the initial state, u0. So if u0 is BigFloat, the integrator’s time is also BigFloat, which the struct doesn't accept

Proposal: Make PoincaréMap parametric over the time type and state element type, e.g.

mutable struct PoincareMap{I<:ContinuousTimeDynamicalSystem,Tv<:AbstractFloat}
    ds::I
    ...
    Tmax::Tv
    tcross::Tv
    dummy::Vector{Tv}
    ...
end

and convert literals like 1e3 into Tv(1e3) inside the constructor. Is there a reason Float64 was fixed originally, or would a PR making this change be welcome?

hibiki-kato avatar Oct 07 '25 22:10 hibiki-kato

no reason whatsoever, please feel free to make a PR that generalizes this to a parameterized type as you suggest!

Datseris avatar Oct 08 '25 10:10 Datseris