Discontinuous observable

[lindnemi]

This issue is a follow up to our slack discussion. To compute the size of the basin of attraction I am using a discontinuous observable (converged). With the default setting the quadrature does not converge. With reltol and abstol >0.08 it always finds the same (approximately correct) answer and for reltol and abstol <0.07 it doesn't converge at all. maxiters seems to have no effect.

using DiffEqUncertainty, OrdinaryDiffEq, Distributions

# Define a system [Menck, 2012]

function swing!(du,u,p,t)
    du[1] = u[2] # phase
    du[2] = -p[1]*u[2] + p[2] - p[3] * sin(u[1]) # frequency
end

tspan = (0.0,1000.0)
p = [0.1, 1, 8] # α, P, K
u_fix = [asin(p[2] / p[3]); 0.]
prob = ODEProblem(swing!, u_fix, tspan, p)


# Distribution of initial conditions

u0_dist = [Uniform(u_fix[1] - pi, u_fix[1] + pi), Uniform(-100, 100)]

# Observable for basin size

# An initial condition belongs to the basin if its frequency converges to 0
# Here we assume a trajectory is converged if its end point deviates from 0 by less than 0.1
converged(sol) = abs(sol[2,end]) > 0.1 ? 0 : 1

# MonteCarlo experiment to compute basin stability
@time expectation(converged, prob, u0_dist, p, MonteCarlo(), Tsit5(); trajectories = 1_000)  # ~1.5s
@time expectation(converged, prob, u0_dist, p, MonteCarlo(), Tsit5(), EnsembleSerial(); trajectories = 1_000) # ~3.5s

# [Gerlach, 2020] claims the Koopman algorithm is faster, but here it doesn't converge
expectation(converged, prob, u0_dist, p, Koopman(), Tsit5())

The discontinuous observable seems to be at the heart of the problem. When I replace it with a bump function the Koopman algorithm behaves reasonably.

function bump(sol)
    z = abs(sol[2,end])
    if z < .1
        return 1
    elseif z < .2
        return exp(-1 / (1 - (z-.1)^2))
    else
        return 0
    end
end

@time quad  = expectation(bump, prob, u0_dist, p, Koopman(), Tsit5(); ireltol = .05, iabstol = .05)

[agerlach]

Moving some of the slack convo here for future reference/discussion:

Some quadrature algorithms can handle the discontinuity better than others. The default does not converge, but I find switching to CubaDivonne() gives the best performance. The default setting for CubaDivonne() utilizes a quasi-MC method. I have not tried tuning any of the CubaDivonne() specific setting.

This also support batch processing, so you can do EnsembleThreads() which is used by default with MonteCarlo() .

@time koop = expectation(converged, prob, u0_dist, p, Koopman(), Tsit5(), EnsembleThreads(); quadalg = CubaDivonne(), iabstol= 1e-3, batch=10)

solves in 5.6s in my machine w/ ~3400 simulations

[agerlach]

https://arxiv.org/pdf/1803.06372.pdf

[ChrisRackauckas]

Cubature methods are just a bad choice for characteristic functions, and this is a common use case for quadrature in the Koopman setup since characteristic functions give the probability of the outcome, which is a pretty important quantity.

To handle this better, I think we probably want to define some new special quadrature method in Quadrature.jl specifically for this type of problem. My idea is that you somehow want to use something like arclength continuation to trace the boundary and then calculate the volume of the convex hull of the points, and that last part has a pretty well-established solution: https://www.cs.princeton.edu/~chazelle/pubs/ConvexHullAlgorithm.pdf. I asked @stevengj if he knew of any methods and he mentioned Mathematica does a rootfinding problem for the boundaries. I am not seeing a mention of it in here though: https://reference.wolfram.com/language/tutorial/NIntegrateIntegrationRules.html .

This would probably be a fun topic but it'll take time. That said, I think for higher dimensions the adaptive Monte Carlo quadrature methods will likely be pretty good because they adapt the sampling to try and ignore the zero region, so for now CubaDivonne is a very good suggestion. Even for higher dimensional cases though, I wonder if there's something to specialize on, since if you know that the function is constant, it's a waste of time to get a bunch of points in the interior to calculate the percentage of the total area. Instead, even for higher dimensions, it seems like you want to build a sampler that tries to sample around the boundary point, and then compute the percentage of points within the boundary (avoiding convex hull algorithms though since they grow exponentially in d which is what MC methods want to avoid).

[ChrisRackauckas]

From @dpsanders, https://github.com/bachrathyd/MDBM.jl looks like a perfect way to do this.