using ValidatedNumerics # pulls in all packages using Plots using IntervalConstraintProgramming, IntervalArithmetic, IntervalRootFinding, Interact @manipulate for a in 0.0:0.1:1.0 C = @constraint sin(y * exp(x) - 1) - $a <= 0.0 X = (-5..5) × (-5..5) # box @time paving = pave(C, X, 0.1) plot(paving.inner, lw=1, label=""); plot!(paving.boundary, lw=0, label="") end C = @constraint (y-5)*cos(4*sqrt((x-4)^2+y^2)) - x*sin(2*sqrt(x^2 + y^2)) >= 0 X = IntervalBox(-5..5, 2) @time paving = pave(C, X, 0.05) plot(paving.inner, lw=0, aspect_ratio=1, leg=false) plot!(paving.boundary, lw=0) using Makie #plotlyjs() f1(x, y) = sin(y*exp(x) - 1) r = -3:0.01:3 Makie.surface(r, r, f1) #surface!(r, r, (x,y)->0.0, alpha=0.5, cbar=false) #zlims!(-1, 1) gr() Plots.contour(r, r, (x,y)->sin(y*exp(x)-1), levels=[0.0, 0.001], cbar=false) Plots.contour!(r, r, (x,y)->(x + 3) * (y^3 - 7) + 18, levels=[0.0, 0.001], cbar=false, color = cgrad([:blue, :blue])) using IntervalConstraintProgramming, Plots, IntervalArithmetic, Interact @manipulate for a in 1:10 C1 = @constraint (x + 3)*(y^3 - $a) + 18 == 0 C2 = @constraint sin(y*exp(x) - 1) == 0 X = IntervalBox(-5..5, 2) p1 = pave(C1, X, 0.1); @time p2 = pave(C2, X, 0.1) plot(p1.boundary, lw=0); plot!(p2.boundary, lw=0, leg=false) @time p3 = pave(C1 ∩ C2, X, 0.1); plot!(p3.boundary, lw=3, label="") end using IntervalRootFinding, StaticArrays X = IntervalBox(-5..5, 2) f(x, y) = SVector( (x + 3) * (y^3 - 7) + 18, sin(y * exp(x) - 1) ) f(xx) = f(xx[1], xx[2]) @time rts = roots(f, X, Newton, 1e-20) display(rts[1:5]); display(length(rts)) all(isunique, rts) gr() r = -1:0.01:5 Plots.contour(r, r, (x,y)->sin(y*exp(x)-1), levels=[0.0, 0.001], cbar=false) Plots.contour!(r, r, (x,y)->(x + 3) * (y^3 - 7) + 18, levels=[0.0, 0.001], cbar=false, color = cgrad([:blue, :blue])) mids = mid.(interval.(rts)) Plots.scatter!(first.(mids), last.(mids), ms=1, aspect_ratio=1) G(X) = 1 + sum(X.^2) / 4000 - prod( cos(X[i] / √i) for i in 1:length(X) ) using IntervalRootFinding @time rts = roots(∇(G), IntervalBox(-10..10, 2)) using IntervalArithmetic, StaticArrays, IntervalRootFinding using Makie # new Julia 2D and 3D plotting package G(x, y) = G([x, y]) mids = mid.(interval.(rts)) vals = G.(mids) N = 10; r = -N:0.1:N Makie.surface( r, r, G, colormap = [(:black, 0.1), (:red, 0.5), (:green, 0.8)] ) Makie.scatter!(first.(mids), last.(mids), G.(mids), markersize=0.8) using Plots, IntervalConstraintProgramming, IntervalRootFinding, IntervalArithmetic using ForwardDiff, LinearAlgebra stability(f, X) = count(eigvals(ForwardDiff.hessian(f, X)) .> 0) # can make rigorous mids = mid.(interval.(rts)) stabilities = stability.(G, mids) mins = mids[stabilities .== 2]; maxs = mids[stabilities .== 0]; stats = mids[stabilities .== 1] N = 50; r = -N:0.1:N Makie.surface( r, r, G, colormap = [(:black, 0.1), (:red, 0.5), (:green, 0.8)], ) Makie.scatter!(first.(mins), last.(mins), G.(mins), markersize=0.8, color=:blue) Makie.scatter!(first.(maxs), last.(maxs), G.(maxs), markersize=0.8, color=:red) Makie.scatter!(first.(stats), last.(stats), G.(stats), markersize=0.8, color=:green) using Plots using StaticArrays, DifferentialEquations, IntervalRootFinding plotlyjs() function parameterized_lorenz(du, u, p, t) du[1] = p[1] * (u[2] - u[1]) du[2] = u[1] * (p[2] - u[3]) - u[2] du[3] = u[1] * u[2] - p[3] * u[3] end lorenz2(u, p) = SVector( p[1]*(u[2]-u[1]), u[1]*(p[2]-u[3]) - u[2], u[1]*u[2] - p[3]*u[3] ) function draw_lorenz(ρ, T) p = [10.0, ρ, 8/3] J = ForwardDiff.jacobian(u->lorenz2(u, p), SVector(0, 0, 0)) λs, vs = eig(Array(J)) which = indmax(λs) # unstable is biggest unstab_direction = vs[:, which] u0 = 1e-5 * unstab_direction tspan = (0.0, T) prob = ODEProblem(parameterized_lorenz, u0, tspan, p) sol = solve(prob) plot(sol, vars=(1, 2, 3), c=:blue, alpha=0.8, lw=1) u0 = -1e-5 * unstab_direction tspan = (0.0, T) prob = ODEProblem(parameterized_lorenz, u0, tspan, p) sol = solve(prob) plot!(sol, vars=(1, 2, 3), c=:red, alpha=0.8, lw=1) X = IntervalBox(-50..51, 3) rts = roots(u->lorenz2(u, p), X) mids = mid.([rt.interval for rt in rts]) scatter!(getindex.(mids, 1), getindex.(mids, 2), getindex.(mids, 3), alpha=1.0, ms=5) xlims!(-30, 30) ylims!(-30, 30) zlims!(-1, 50) end using Interact, IntervalArithmetic # plotlyjs() # @manipulate throttle=0.5 for ρ in slider(0:0.1:100, value=28.0), T in 1:0.1:100 # draw_lorenz(ρ, T) # end using IntervalOptimisation, IntervalArithmetic @manipulate for n in Interact.slider(1:100, value=1) @time global_min, minimizers = minimise(G, IntervalBox(-600..600, n)) global_min, minimizers end using IntervalOptimisation, IntervalArithmetic f(x, y) = (1.5 - x * (1 - y))^2 + (2.25 - x * (1 - y^2))^2 + (2.625 - x * (1 - y^3))^2 f(X) = f(X...) #X = (-1e6..1e6) × (-1e6..1e6) X = minimise_icp(f, X) using IntervalArithmetic, IntervalOptimisation, Interact, Plots, IntervalConstraintProgramming gr() X = IntervalBox(-2..2, 2) C = @constraint x^2 + y^2 <= 1 paving1 = pave(C, X, 0.05) @manipulate for a in -2:0.01:2 C2 = @constraint x + y - $a == 0 paving2 = pave(C2, X, 0.01) plot(paving1.inner, lw=0, aspect_ratio=1, leg=false); plot!(paving1.boundary, lw=0, aspect_ratio=1) plot!(paving2.boundary) xlims!(-2, 2); ylims!(-2, 2) end f(X) = (X[1] + X[2]) X = (-∞..∞) × (-∞..∞) constraints = [IntervalOptimisation.Constraint(x->(x[1]^2 + x[2]^2), -∞..1)] @time minimise_icp_constrained(f, X, constraints, abstol=1e-5) using Pkg # built-in package manager Pkg.add("IntervalArithmetic") # install - do only once using IntervalArithmetic, BenchmarkTools X = interval(0.1, 0.3) # checked interval constructor typeof(X) # type is inferred X = 0.1..0.3 # automatic outward rounding showfull(stdout, X) f(x) = x^2 - 2 # one-line mathematical function definition f(X) X, 0 ∈ f(X) # Unicode operators via LaTeX completions: \in ∈ ∈ X = 3..∞ # \infty ∞ Y = f(X) (6..8) ⊆ Y (6..8) ∩ Y (-1..2)^2 sin(3..4 ) X = big(π)..big(π) typeof(X) showfull(stdout, X) diam(f(X)) using ErrorfreeArithmetic # developed by Jeffrey Sarnoff two_sum(0.1, 0.3) # returns pair of `Float64`s that sum exactly to `0.1 + 0.3` @edit two_sum(0.1, 0.3) +(0.1, 0.3, RoundDown), +(0.1, 0.3, RoundUp) X, Y = (0.1..0.2), (0.3..0.4) @btime $X + $Y # use $ interpolation to avoid performance hit from global variables setrounding(Interval, :none) # redefines all functions @btime $X + $Y setrounding(Interval, :accurate) @btime $X + $Y setrounding(Interval, :slow) # change rounding mode four times @btime $X + $Y setrounding(Interval, :tight) # reset for rest of talk! using IntervalArithmetic using Base.Threads, BenchmarkTools function f!(y, x) # f1!(y, x) = ( y .= x ./ 2 ) for i in 1:length(x) y[i] = x[i] / 2 end end function f_threaded!(y, x) @threads for i in 1:length(x) y[i] = x[i] / 2 end end x = [ (i..(i+1)) for i in 1:1000] # array comprehension y = similar(x); y2 = similar(x); @btime f!($y, $x) @btime f_threaded!($y2, $x) y == y2 12.5 / 3.78 setformat(decorations=true) x = -1..2 sqrt(x) x = DecoratedInterval(-1..2) typeof(x) sqrt(x) function mince(x::Interval, n) nodes = range(x.lo, stop=x.hi, length=n+1) return [Interval(nodes[i], nodes[i+1]) for i in 1:length(nodes) - 1] end f(x) = cos(x) + 0.5*sin(2*x) function image(f, X) II = f(X) if II.lo == -Inf II = -100..(II.hi) end return II end intervals = mince(-5..5, 10) images = image.(f, intervals) using Interact, Plots, IntervalArithmetic; gr() function bound_function(f, X, N ,ylims=(-2,2)) @manipulate for N in Interact.slider(1:200, value=1) intervals = mince(X, N) images = image.(f, intervals) boxes = [IntervalBox(intervals[i], images[i]) for i in 1:length(intervals)] contains_zero = [0 ∈ images[i] for i in 1:length(images)] if length(boxes[contains_zero]) > 0 plot(boxes[contains_zero], ylim=ylims) # uses "plot recipe" for IntervalBoxes end if length(boxes[.!(contains_zero)]) > 0 plot!(boxes[.!(contains_zero)], ylim=ylims) end plot!(f, c=:green, lw=3) end end f(x) = cos(x) + 0.5*sin(2*x) bound_function(f, -5..5, 50, (-2,2)) f(x) = (1/50)*log(abs(3*(1-x) + 1)) + x^2 + 1 bound_function(f, 0.8..2.0, 100, (1,4)) bound_function(x->sin(1/x), 0.01..0.2, (-1,1)) function integrate(f, X, N) intervals = mince(X, N) return sum(diam.(intervals) .* f.(intervals)) end integrate(x->x^2, 0..1, 100000) X = (1..2) × (3..4) # \times typeof(X) IntervalBox(-10..10, 50) Y = bisect(X, 0.5) typeof(Y) X1, X2 = Y X1 X ∩ Y[1] # \cap X ∪ Y[1] # \cup setdiff(X, Y[1]) struct MyInterval{T <: AbstractFloat} # type parameter T must be subtype of AbstractFloat lo::T hi::T end import Base: + +(x::MyInterval{T}, y::MyInterval{T}) where {T} = MyInterval(x.lo + y.lo, x.hi + y.hi) x = MyInterval(0.1, 0.3) y = MyInterval(0.5, 0.6) # automatically *infers* correct type T x + y @code_native x + y x = MyInterval(big"1.2", big"3.4") x + x X = (1..2) × (3..4) # \times dump(X) using ForwardDiff # also ReverseDiff ForwardDiff.derivative(X -> X^2 - 2X, 4..5) ff(XX) = ( (X, Y) = XX; [X^2 + Y^2, X - Y] ) J = ForwardDiff.jacobian(ff, [big(3..4), 5..6]) # J^-1 # inv(J) J, inv(Interval{Float64}.(J)) methods(roots) G(X) = 1 + sum(abs2, X) / 4000 - prod( cos(X[i] / √i) for i in 1:length(X) ) using IntervalArithmetic, IntervalRootFinding @time rts = roots(∇(G), IntervalBox(-600..600, 2)); length(rts) all(isunique, rts) function simple_minimize(f, X, ϵ=1e-3) L = [X]; R = typeof(X)[ ]; m = ∞ while !isempty(L) X = popfirst!(L) # first element: queue f(X).lo > m && continue # short-circuiting if-then if diam(X) < ϵ push!(R, X) continue end current_max = f(Interval.(mid.(X))) # no intermediate temporary is created: "fusion" # pointwise function application ("broadcast") current_max < m && (m = current_max) append!(L, bisect(X)) end return m, R end @time simple_minimize(G, IntervalBox(-600..600, 50)) @time simple_minimize(G, IntervalBox(-600..600, 2)) struct QuadraticForm A::Matrix{Float64} end (Q::QuadraticForm)(x) = dot(x, Q.A * x) Q = QuadraticForm([2 0; 0 3]) Q([3, 4]) import Base: * *(M::Matrix, X::IntervalBox{2,T}) where {T} = *(M, X.v) Q(IntervalBox(3, 4)) @time simple_minimize(Q, IntervalBox(-1000..1000, 2)) @time simple_minimize(G, -600..600) @time minval, minimisers = minimise(G, IntervalBox(-600..600, 50)) using IntervalConstraintProgramming, IntervalArithmetic, Plots, Interact C = @constraint 2 <= x^3 + y^3 <= 5 X = (-5..5) × (-5..5) # \times @manipulate for n in Interact.slider(0.5:0.5:8, value=0.5) ϵ = 2.0^(-n) paving = pave(C, X, ϵ) plot(paving.inner, aspect_ratio=1, leg=false) plot!(paving.boundary, ann=(0, 0, "$(Vol(paving))")) end function bisection(f, X, ϵ=1e-3) L = [X]; R = [ ] while !isempty(L) X = pop!(L) if diam(X) < ϵ push!(R, X) continue end if 0 ∈ f(X) X1, X2 = bisect(X) push!(L, X1, X2) end end return R end bisection(x->x^2 - 2, -∞..∞) function minimize(f, X, ϵ=1e-3) L = [X]; R = typeof(X)[ ]; m = ∞ while !isempty(L) X = shift!(L) # first element: queue if f(X).lo > m continue end if diam(X) < ϵ push!(R, X) continue end current_max = f(Interval(mid(X))) if current_max < m m = current_max end push!(L, bisect(X)...) end return m, R end using IntervalArithmetic, IntervalRootFinding, Plots plotlyjs() m, R = minimize(x->(x^2 - 2)^2 - 2, -10..11) carrom(x, y) = - (exp(2*abs(1 - sqrt(x^2 + y^2) / pi)) * cos(x)^2 * cos(y)^2) / 30 carrom(X) = ( (x, y) = X; carrom(x, y) ) plotlyjs() surface(-5:0.1:5, -5:0.1:5, carrom) using IntervalArithmetic, IntervalOptimisation, Plots @time globalmin, minimisers = minimise(carrom, (-10..10) × (-10..10), 2e-3) gr() contour(-10:0.05:10, -10:0.05:10, carrom, aspect_ratio=1) plot!(minimisers) Meta.show_sexpr(:(x^2 + y^2)) using IntervalArithmetic x = y = -∞..∞ # intervals a = x^2 b = y^2 c = a + b c = c ∩ (-∞..1) c a, (c - b) # curre ntly inconsistent a = a ∩ (c - b) # make consistent x = √a ∪ (-√a) # x^2 = a using IntervalConstraintProgramming C = @constraint x^2 + y^2 <= 1 #using MacroTools C.contractor.backward.code using Interact C = @constraint 2 <= x^3 + y^2 <= 5 X = (-5..5) × (-5..5) # \times @manipulate for n in slider(0.5:0.5:8, value=0.5) eps = 2.0^(-n) paving = pave(C, X, eps) plot(paving.inner, aspect_ratio=:equal) plot!(paving.boundary, ann=(0, 0, "$(Vol(paving))")) end straight(x0, y0, x, m) = y0 + m * (x - x0) function interval_newton_full(f, X0, ymin=-10, ymax=10) @manipulate for n in slider(1:20, value=1), α in linspace(0, 1, 50) X = X0 # draw graph of function over interval X xx = X.lo:0.0001:X.hi p = plot(xx, map(f, xx), c="blue", lw=3, xlim=(X.lo, X.hi), ylim=(ymin, ymax), legend=:false) hline!([0], color="magenta", lw=3, linestyle=:dash) Xs = [X] new_Xs = [] N1 = ∅ N2 = ∅ for i in 1:n-1 for X in Xs x0 = (1-α)*X.lo + α*X.hi # mid(X) deriv = ForwardDiff.derivative(f, X) if 0 ∈ deriv N1 = x0 - f(@interval(x0)) / @interval(deriv.lo, -0.0) N2 = x0 - f(@interval(x0)) / @interval(0.0, deriv.hi) N1 = N1 ∩ X N2 = N2 ∩ X if !(isempty(N1)) push!(new_Xs, N1) end if !(isempty(N2)) push!(new_Xs, N2) end else N1 = x0 - f(@interval(x0)) / deriv N1 = N1 ∩ X if !(isempty(N1)) push!(new_Xs, N1) end end end Xs = new_Xs new_Xs = [] end for X in Xs #if n > 1 plot!([X.lo, X.hi], [0,0], c="cyan", linewidth=4, alpha=0.3) m = (1-α)*X.lo + α*X.hi scatter!([m], [0], c="green") scatter!([m], [f(m)], c="red") plot!([m, m], [0, f(m)], c="green", ls=:dash) #end x0 = (1-α)*X.lo + α*X.hi # mid(X) deriv = ForwardDiff.derivative(f, X) # draw initial point y0 = f(x0) scatter!([x0], [y0], c="red") # draw cone for m in linspace(deriv.lo, deriv.hi, 100) plot!([X.lo, X.hi], [ y0 + m*(x-x0) for x in [X.lo, X.hi]], color="gray", alpha=0.2) end if 0 ∈ deriv N1 = x0 - f(@interval(x0)) / @interval(deriv.lo, -0.0) N2 = x0 - f(@interval(x0)) / @interval(0.0, deriv.hi) N1 = N1 ∩ X N2 = N2 ∩ X plot!([N1.lo, N1.hi], [0,0], c="red", linewidth=4, alpha=0.8) plot!([N2.lo, N2.hi], [0,0], c="red", linewidth=4, alpha=0.8) else N1 = x0 - f(@interval(x0)) / deriv N1 = N1 ∩ X plot!([N1.lo, N1.hi], [0,0], c="red", linewidth=4, alpha=0.8) end end #text(0, 5, "$(length(Xs))") p end end xx = -1:0.01:10 f(x) = x^2 - 2 x0 = 8 interval_newton_full(x->f(x), 0..4, -10, 26) interval_newton_full(x->x^2 - 2, -3..5, -10, 10) gr() @manipulate for a in -2:0.01:2 C1 = @constraint x + 2y - $a == 0 X = IntervalBox(0..5, -5..1) paving1 = pave(C1, X, 0.01) C2 = @constraint x + y <= 1 paving2 = pave(C2, X, 0.01) plot(paving1.boundary, aspect_ratio=1) plot!(paving2.inner, lw=0) xlims!(-2, 8) ylims!(-10, 8) end