Merge pull request #796 from JuliaReach/schillic/examples

Better script output (#!jl) in examples and unused packages removed from test/Project.toml

examples/Brusselator/Brusselator.jl

@@ -25,7 +25,7 @@
 # The numerical values for the model's constants (in their respective units) are
 # ``A = 1, B = 1.5``.
 
-using ReachabilityAnalysis  #!jl
+using ReachabilityAnalysis
 
 @taylorize function brusselator!(du, u, p, t)
     local A = 1.0
@@ -49,19 +49,19 @@ end
 # The initial set is ``x ∈ [0.8, 1]``, ``y ∈ [0, 0.2]``. The time horizon is
 # ``18``. These settings are taken from [^CAS13].
 
-U₀ = (0.8 .. 1.0) × (0.0 .. 0.2)  #!jl
-prob = @ivp(u' = brusselator!(u), u(0) ∈ U₀, dim:2)  #!jl
+U₀ = (0.8 .. 1.0) × (0.0 .. 0.2)
+prob = @ivp(u' = brusselator!(u), u(0) ∈ U₀, dim:2)
 
-T = 18.0;  #!jl
+T = 18.0;
 
 # ## Analysis
 
 # We use the `TMJets` algorithm with sixth-order expansion in time and
 # second-order expansion in the spatial variables. The version `TMJets20` is
 # more efficient here.
 
-alg = TMJets20(; orderT=6, orderQ=2)  #!jl
-sol = solve(prob; T=T, alg=alg);  #!jl
+alg = TMJets20(; orderT=6, orderQ=2)
+sol = solve(prob; T=T, alg=alg);
 
 # ## Results
 
@@ -97,7 +97,7 @@ bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim:2);
 
 # First we solve for ``r = 0.01``:
 
-sol_01 = solve(bruss(0.01); T=30.0, alg=alg);  #!jl
+sol_01 = solve(bruss(0.01); T=30.0, alg=alg);
 
 #-
 
@@ -115,7 +115,7 @@ plot!(fig, U0(0.01); color=:orange, lab="Uo", xlims=(0.6, 1.3))  #!jl
 # flowpipe zoomed to the last portion and compare ``r = 0.01`` with a larger set
 # of initial states, ``r = 0.1``.
 
-sol_1 = solve(bruss(0.1); T=30.0, alg=alg);  #!jl
+sol_1 = solve(bruss(0.1); T=30.0, alg=alg);
 
 #-
 

examples/Building/Building.jl

@@ -45,6 +45,7 @@
 
 using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath
 using ReachabilityBase.Arrays: SingleEntryVector
+using ReachabilityAnalysis: add_dimension
 
 const x25 = SingleEntryVector(25, 48, 1.0)
 const x25e = SingleEntryVector(25, 49, 1.0)
@@ -71,8 +72,6 @@ function building_BLDF01()
     return prob_BLDF01
 end
 
-using ReachabilityAnalysis: add_dimension
-
 function building_BLDC01()
     A, B, X0, U = building_common()
     n = size(A, 1)
@@ -127,42 +126,42 @@ prob_BLDF01 = building_BLDF01();
 sol_BLDF01_dense = solve(prob_BLDF01; T=20.0,
                          alg=LGG09(; δ=0.004, vars=(25), n=48))
 
-fig = plot(sol_BLDF01_dense; vars=(0, 25), linecolor=:blue, color=:blue,
-           alpha=0.8, lw=1.0, xlab="t", ylab="x25")
+fig = plot(sol_BLDF01_dense; vars=(0, 25), linecolor=:blue, color=:blue,  #!jl
+           alpha=0.8, lw=1.0, xlab="t", ylab="x25")  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 # Safety properties:
 
 @assert ρ(x25, sol_BLDF01_dense) <= 5.1e-3 "the property should be proven"  # BLDF01 - BDS01
 @assert !(ρ(x25, sol_BLDF01_dense) <= 4e-3) "the property should not be proven"  # BLDF01 - BDU01
-ρ(x25, sol_BLDF01_dense)
+ρ(x25, sol_BLDF01_dense)  #!jl
 
 #-
 
 @assert !(ρ(x25, sol_BLDF01_dense(20.0)) <= -0.78e-3) "the property should not be proven"  # BLDF01 - BDU02
-ρ(x25, sol_BLDF01_dense(20.0))
+ρ(x25, sol_BLDF01_dense(20.0))  #!jl
 
 # #### Discrete time
 
 sol_BLDF01_discrete = solve(prob_BLDF01; T=20.0,
                             alg=LGG09(; δ=0.01, vars=(25), n=48, approx_model=NoBloating()));
 
-fig = plot(sol_BLDF01_discrete; vars=(0, 25), linecolor=:blue, color=:blue,
-           alpha=0.8, lw=1.0, xlab="t", ylab="x25")
+fig = plot(sol_BLDF01_discrete; vars=(0, 25), linecolor=:blue, color=:blue,  #!jl
+           alpha=0.8, lw=1.0, xlab="t", ylab="x25")  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 # Safety properties:
 
 @assert ρ(x25, sol_BLDF01_discrete) <= 5.1e-3 "the property should be proven"  # BLDF01 - BDS01
 @assert !(ρ(x25, sol_BLDF01_discrete) <= 4e-3) "the property should not be proven"  # BLDF01 - BDU01
-ρ(x25, sol_BLDF01_discrete)
+ρ(x25, sol_BLDF01_discrete)  #!jl
 
 #-
 
 @assert !(ρ(x25, sol_BLDF01_discrete(20.0)) <= -0.78e-3) "the property should not be proven"  # BLDF01 - BDU02
-ρ(x25, sol_BLDF01_discrete(20.0))
+ρ(x25, sol_BLDF01_discrete(20.0))  #!jl
 
 # ### BLDC01
 
@@ -172,42 +171,42 @@ prob_BLDC01 = building_BLDC01();
 
 sol_BLDC01_dense = solve(prob_BLDC01; T=20.0, alg=LGG09(; δ=0.005, vars=(25), n=49))
 
-fig = plot(sol_BLDC01_dense; vars=(0, 25), linecolor=:blue, color=:blue,
-           alpha=0.8, lw=1.0, xlab="t", ylab="x25")
+fig = plot(sol_BLDC01_dense; vars=(0, 25), linecolor=:blue, color=:blue,  #!jl
+           alpha=0.8, lw=1.0, xlab="t", ylab="x25")  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 # Safety properties
 
 @assert ρ(x25e, sol_BLDC01_dense) <= 5.1e-3 "the property should be proven"  # BLDC01 - BDS01
 @assert !(ρ(x25e, sol_BLDC01_dense) <= 4e-3) "the property should not be proven"  # BLDC01 - BDU01
-ρ(x25e, sol_BLDC01_dense)
+ρ(x25e, sol_BLDC01_dense)  #!jl
 
 #-
 
 @assert !(ρ(x25e, sol_BLDC01_dense(20.0)) <= -0.78e-3) "the property should not be proven"  # BLDC01 - BDU02
-ρ(x25, sol_BLDF01_discrete(20.0))
+ρ(x25, sol_BLDF01_discrete(20.0))  #!jl
 
 # #### Discrete time
 
 sol_BLDC01_discrete = solve(prob_BLDC01; T=20.0,
                             alg=LGG09(; δ=0.01, vars=(25), n=49, approx_model=NoBloating()))
 
-fig = plot(sol_BLDC01_discrete; vars=(0, 25), linecolor=:blue, color=:blue,
-           alpha=0.8, lw=1.0, xlab="t", ylab="x25")
+fig = plot(sol_BLDC01_discrete; vars=(0, 25), linecolor=:blue, color=:blue,  #!jl
+           alpha=0.8, lw=1.0, xlab="t", ylab="x25")  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 # Safety properties
 
 @assert ρ(x25e, sol_BLDC01_discrete) <= 5.1e-3 "the property should be proven" # BLDC01 - BDS01
 @assert !(ρ(x25e, sol_BLDC01_discrete) <= 4e-3) "the property should not be proven" # BLDC01 - BDU01
-ρ(x25e, sol_BLDC01_discrete)
+ρ(x25e, sol_BLDC01_discrete)  #!jl
 
 #-
 
 @assert !(ρ(x25e, sol_BLDC01_discrete(20.0)) <= -0.78e-3) "the property should not be proven" # BLDC01 - BDU02
-ρ(x25e, sol_BLDC01_discrete(20.0))
+ρ(x25e, sol_BLDC01_discrete(20.0))  #!jl
 
 # ## References
 

examples/DuffingOscillator/DuffingOscillator.jl

@@ -7,7 +7,7 @@
 
 # ## Model description
 
-using ReachabilityAnalysis  #!jl
+using ReachabilityAnalysis
 
 const ω = 1.2
 
@@ -42,18 +42,18 @@ sol = solve(prob; tspan=(0.0, 20 * T), alg=TMJets21a());
 using Plots  #!jl
 #!jl import DisplayAs  #hide
 
-fig = plot(sol; vars=(0, 1), xlab="t", ylab="x", lw=0.2, color=:blue)
+fig = plot(sol; vars=(0, 1), xlab="t", ylab="x", lw=0.2, color=:blue)  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 #-
 
-fig = plot(sol; vars=(0, 2), xlab="t", ylab="v", lw=0.2, color=:blue)
+fig = plot(sol; vars=(0, 2), xlab="t", ylab="v", lw=0.2, color=:blue)  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 #-
 
-fig = plot(sol; vars=(1, 2), xlab="x", ylab="v", lw=0.5, color=:red)
+fig = plot(sol; vars=(1, 2), xlab="x", ylab="v", lw=0.5, color=:red)  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide

examples/EMBrake/EMBrake.jl

@@ -33,7 +33,7 @@
 # - **Constant inputs**: The inputs are only uncertain in the initial value, and
 #    constant over time: ``u(0) ∈ \mathcal{U}``, ``\dot{u}(t) = 0.``
 
-using ReachabilityAnalysis, SparseArrays  #!jl
+using ReachabilityAnalysis, SparseArrays
 
 # ## Common functionality between model variants
 

examples/Heat3D/Heat3D.jl

@@ -24,6 +24,6 @@
 # - **Constant inputs**: The inputs are only uncertain in the initial value, and
 #    constant over time: ``u(0) ∈ \mathcal{U}``, ``\dot{u}(t) = 0.``
 
-using ReachabilityAnalysis, SparseArrays, JLD2, ReachabilityBase.CurrentPath  #!jl
+using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath
 
 path = @current_path("Heat3D", "HEAT01.jld2")

examples/ISS/ISS.jl

@@ -35,7 +35,7 @@
 # - **ISSC01** (constant inputs): The inputs are only uncertain in their initial
 #   value, and constant over time: ``u(0) ∈ \mathcal{U}``, ``\dot{u}(t) = 0``.
 
-using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath  #!jl
+using ReachabilityAnalysis, JLD2, ReachabilityBase.CurrentPath
 using ReachabilityAnalysis: add_dimension
 
 path = @current_path("ISS", "ISS.jld2")
@@ -121,17 +121,17 @@ sol_ISSF01 = solve(prob_ISSF01; T=20.0,
 # The solution `sol_ISSF01` is a 270-dimensional set that (only) contains
 # template reach sets for the linear combination `C_3 x(t)`.
 
-dim(sol_ISSF01)
+dim(sol_ISSF01)  #!jl
 
 # For visualization, it is necessary to specify that we want to plot "time" vs.
 # ``y_3(t)``. We can transform the flowpipe on the output ``y_3(t)`` by
 # "flattening" the flowpipe along directions `dirs`.
 
-πsol_ISSF01 = flatten(sol_ISSF01);
+πsol_ISSF01 = flatten(sol_ISSF01);  #!jl
 
 # Now `πsol_ISSF01` is a one-dimensional flowpipe.
 
-dim(πsol_ISSF01)
+dim(πsol_ISSF01)  #!jl
 
 #md # !!! tip "Performance tip"
 #md #     Note that projecting the solution along direction ``C_3`` corresponds
@@ -149,11 +149,11 @@ LazySets.set_ztol(Float64, 1e-8);  # use higher precision for the plots #!jl
 
 #-
 
-fig = Plots.plot(πsol_ISSF01[1:10:end]; vars=(0, 1), linecolor=:blue, color=:blue,
-                 alpha=0.8, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0],
-                 ytick=[-0.00075, -0.0005, -0.00025, 0, 0.00025, 0.0005, 0.00075],
-                 xlims=(0.0, 20.0), ylims=(-0.00075, 0.00075),
-                 bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm)
+fig = plot(πsol_ISSF01[1:10:end]; vars=(0, 1), linecolor=:blue, color=:blue,  #!jl
+           alpha=0.8, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0],  #!jl
+           ytick=[-0.00075, -0.0005, -0.00025, 0, 0.00025, 0.0005, 0.00075],  #!jl
+           xlims=(0.0, 20.0), ylims=(-0.00075, 0.00075),  #!jl
+           bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm)  #!jl
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 # ### ISSC01
@@ -165,15 +165,15 @@ sol_ISSC01 = solve(prob_ISSC01; T=20.0,
 
 # We can flatten the flowpipe to the output ``y_3(t)`` as before:
 
-πsol_ISSC01 = flatten(sol_ISSC01);
+πsol_ISSC01 = flatten(sol_ISSC01);  #!jl
 
 #-
 
-fig = Plots.plot(πsol_ISSC01; vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8,
-                 lw=1.0, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0],
-                 ytick=[-0.0002, -0.0001, 0.0, 0.0001, 0.0002],
-                 xlims=(0.0, 20.0), ylims=(-0.0002, 0.0002),
-                 bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm)
+fig = plot(πsol_ISSC01; vars=(0, 1), linecolor=:blue, color=:blue, alpha=0.8,  #!jl
+           lw=1.0, xlab=L"t", ylab=L"y_{3}", xtick=[0, 5, 10, 15, 20.0],  #!jl
+           ytick=[-0.0002, -0.0001, 0.0, 0.0001, 0.0002],  #!jl
+           xlims=(0.0, 20.0), ylims=(-0.0002, 0.0002),  #!jl
+           bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm)  #!jl
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide
 
 #-

examples/LaubLoomis/LaubLoomis.jl

@@ -11,7 +11,7 @@
 # enzymatic activities. The dynamics can be defined by the following
 # seven-dimensional system of differential equations.
 
-using ReachabilityAnalysis  #!jl
+using ReachabilityAnalysis
 
 @taylorize function laubloomis!(dx, x, p, t)
     dx[1] = 1.4 * x[3] - 0.9 * x[1]
@@ -68,11 +68,11 @@ sol_1z = overapproximate(sol_1, Zonotope);
 # We verify that the specification holds:
 
 @assert ρ(e4, sol_1z) < 4.5 "the property should be proven"
-ρ(e4, sol_1z)
+ρ(e4, sol_1z)  #!jl
 
 # To compute the width of the final box, we use the support function:
 
-ρ(e4, sol_1z[end]) + ρ(-e4, sol_1z[end])
+ρ(e4, sol_1z[end]) + ρ(-e4, sol_1z[end])  #!jl
 
 # ### Case 2: Intermediate initial set
 
@@ -87,11 +87,11 @@ sol_2z = overapproximate(sol_2, Zonotope);
 # We verify that the specification holds:
 
 @assert ρ(e4, sol_2z) < 4.5 "the property should be proven"
-ρ(e4, sol_2z)
+ρ(e4, sol_2z)  #!jl
 
 # Width of final box:
 
-ρ(e4, sol_2z[end]) + ρ(-e4, sol_2z[end])
+ρ(e4, sol_2z[end]) + ρ(-e4, sol_2z[end])  #!jl
 
 # ### Case 3: Large initial set
 
@@ -106,30 +106,27 @@ sol_3z = overapproximate(sol_3, Zonotope);
 # We verify that the specification holds:
 
 @assert ρ(e4, sol_3z) < 5.0 "the property should be proven"
-ρ(e4, sol_3z)
+ρ(e4, sol_3z)  #!jl
 
 # Width of final box:
 
-ρ(e4, sol_3z[end]) + ρ(-e4, sol_3z[end])
+ρ(e4, sol_3z[end]) + ρ(-e4, sol_3z[end])  #!jl
 
 # ## Results
 
 using Plots, Plots.PlotMeasures, LaTeXStrings  #!jl
 #!jl import DisplayAs  #hide
 
-fig = plot(sol_3z; vars=(0, 4), linecolor="green", color=:green, alpha=0.8, lab="W = 0.1")
-plot!(fig, sol_2z; vars=(0, 4), linecolor="blue", color=:blue, alpha=0.8, lab="W = 0.05")
-plot!(fig, sol_1z; vars=(0, 4), linecolor="yellow", color=:yellow, alpha=0.8, lab="W = 0.01",
-      tickfont=font(10, "Times"), guidefontsize=15,
-      xlab=L"t",
-      ylab=L"x_4",
-      xtick=[0.0, 5.0, 10.0, 15.0, 20.0], ytick=[1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0],
-      xlims=(0.0, 20.0), ylims=(1.5, 5.02),
-      bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
-      size=(600, 600))
-
-plot!(fig, x -> x, x -> 4.5, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="")
-plot!(fig, x -> x, x -> 5.0, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="")
+fig = plot(sol_3z; vars=(0, 4), linecolor="green", color=:green, alpha=0.8, lab="W = 0.1")  #!jl
+plot!(fig, sol_2z; vars=(0, 4), linecolor="blue", color=:blue, alpha=0.8, lab="W = 0.05")  #!jl
+plot!(fig, sol_1z; vars=(0, 4), linecolor="yellow", color=:yellow, alpha=0.8, lab="W = 0.01",  #!jl
+      tickfont=font(10, "Times"), guidefontsize=15, xlab=L"t", ylab=L"x_4",  #!jl
+      xtick=[0.0, 5.0, 10.0, 15.0, 20.0], ytick=[1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0],  #!jl
+      xlims=(0.0, 20.0), ylims=(1.5, 5.02),  #!jl
+      bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,  #!jl
+      size=(600, 600))  #!jl
+plot!(fig, x -> x, x -> 4.5, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="")  #!jl
+plot!(fig, x -> x, x -> 5.0, 0.0, 20.0; line=2, color="red", linestyle=:dash, lab="")  #!jl
 
 #!jl DisplayAs.Text(DisplayAs.PNG(fig))  #hide