diff --git a/Project.toml b/Project.toml index 7aa00bd34..ec61d8032 100644 --- a/Project.toml +++ b/Project.toml @@ -29,6 +29,7 @@ Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/demo/plotting/Visualise_random_slicing.jl b/demo/plotting/Visualise_random_slicing.jl new file mode 100644 index 000000000..8823090c8 --- /dev/null +++ b/demo/plotting/Visualise_random_slicing.jl @@ -0,0 +1,80 @@ +using UncertaintyQuantification, DataFrames, Plots, ColorSchemes +using UncertaintyQuantification: lo, hi + +X1 = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(0, 2), :σ => 4)), :X1) +X2 = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(0, 2), :σ => Interval(5, 5.5))), :X2) + +function limitstate(df) + return 300 .- 2* df.X2 .^2 + df.X1 .+ df.X2 .^3 - 3 *df.X1.^2 .+ df.X1 .*df.X2 +end + +model = Model(limitstate, :g) + +inputs =[X1; X2] + +@time samples = sample(inputs, SobolSampling(1000)) + +propagate_intervals!(model, samples) + +gs = samples.g + +safe = lo.(gs) .>= 0 +fail = hi.(gs) .<= 0 +unsure = lo.(gs) .<= 0 .&& 0 .<= hi.(gs) + +N_grid = 1000 + +xs_physical = range(-15, 12, length = N_grid) +ys_physical = range(-10, 10, length = N_grid) + +XY_physical = [[x; y] for x in xs_physical, y in ys_physical] +XY_physical = reduce(hcat, XY_physical[:]) + +samples_physical = DataFrame(:X1 => XY_physical[1,:], :X2 => XY_physical[2,:]) + +evaluate!(model, samples_physical) +contour(xs_physical, ys_physical, samples_physical.g .<= 0, color = :red, linewidth = 2, levels = 1, label = "g", colorbar =false) + + +plot(samples.X1[safe], samples.X2[safe], alpha = 0.2, color = "blue", label = "safe") +plot!(samples.X1[fail], samples.X2[fail], alpha = 0.2, color = "red", label = "failed") +plot!(samples.X1[unsure], samples.X2[unsure], alpha = 0.2, color = "yellow", label = "indetermined") +contour!(xs_physical, ys_physical, samples_physical.g .<= 0, color = :red, linewidth = 2, levels = 1, label = "g", colorbar =false) + + +### Subset viz, with tighter dists + +X1_ = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(0, 2), :σ => 2)), :X1) +X2_ = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(-1, 1), :σ => 1.5)), :X2) + +inputs2 = [X1_, X2_] + +ss = SubSetInfinity(1000, 0.1, 10, 0.5) +@time pf_ss, outputs_ss1, outputs_ss2 = probability_of_failure(limitstate, inputs2, RandomSlicing(ss)) + +samples_lo = outputs_ss1[2] +N_levels = maximum(samples_lo.level) + +colour = colormap("RdBu", N_levels) + +p1 = plot(samples_lo.X1[samples_lo.level .== 1], samples_lo.X2[samples_lo.level .==1], color = colour[1], label ="level1", alpha = 0.2) +for i = 2:N_levels + plot!(p1, samples_lo.X1[samples_lo.level .==i], samples_lo.X2[samples_lo.level .==i], color = colour[i], label ="level$i", alpha = 0.2) +end +contour!(p1, xs_physical, ys_physical, samples_physical.g .<= 0, color = :red, linewidth = 2, levels = 1, label = true, colorbar =false, title = "lower bound") + +## Lower bound + +samples_lo = outputs_ss2[2] +N_levels = maximum(samples_lo.level) + +# colour = colormap("RdBu", N_levels) + +p2 = plot(samples_lo.X1[samples_lo.level .== 1], samples_lo.X2[samples_lo.level .==1], color = colour[1], label ="level1", alpha = 0.2) +for i = 2:N_levels + plot!(p2, samples_lo.X1[samples_lo.level .==i], samples_lo.X2[samples_lo.level .==i], color = colour[i], label ="level$i", alpha = 0.2) +end +contour!(p2, xs_physical, ys_physical, samples_physical.g .<= 0, color = :red, linewidth = 2, levels = 1, label = true, colorbar =false, title = "upper bound") + +# Side by side +p = plot(p1, p2, layout = (1, 2), size = (1200, 500)) \ No newline at end of file diff --git a/demo/plotting/plot_inputs.jl b/demo/plotting/plot_inputs.jl new file mode 100644 index 000000000..1168c0de9 --- /dev/null +++ b/demo/plotting/plot_inputs.jl @@ -0,0 +1,24 @@ +using UncertaintyQuantification, Plots + +X1 = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(-1, 2), :σ => 1)), :X1) +X2 = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(-1, 2), :σ => 2)), :X2) +X3 = IntervalVariable(-1, 2, :s) +h = RandomVariable(Normal(0.24, 0.01), :h) # height + +plot(X1) # p-box +plot(X1, color = "red") # red p-box +plot(X3) # Interval +plot(h) # distribution + +inputs = [X1, X2, X3, h] + +plot(inputs) # Everything together + +samples = sample(inputs, 200) + +plot(X1) +plot!(samples.X1) # Samples ecdf samples of X1 + +plot(samples.X1[1], samples.X2[1]) # Plots 2D box + +plot(samples.X1, samples.X2) # Plots bivariate random sets of X1 and X2 diff --git a/docs/src/manual/gettingstarted.md b/docs/src/manual/gettingstarted.md index c2da09cb3..c3da62343 100644 --- a/docs/src/manual/gettingstarted.md +++ b/docs/src/manual/gettingstarted.md @@ -47,6 +47,23 @@ to_standard_normal_space!(x, samples) to_physical_space!(x, samples) ``` +Examples of some distributions available from *Distributions.jl* are shown below. + +```@example rv +using Plots #hide + +inputs = [ #hide + RandomVariable(Normal(0, 1), :normal), #hide + RandomVariable(Beta(2, 5), :beta), #hide + RandomVariable(Gamma(2, 1), :gamma), #hide + RandomVariable(Exponential(1), :exponential), #hide + RandomVariable(Uniform(0, 0.3), :uniform), #hide + RandomVariable(LogNormal(0, 1), :lognormal), #hide +]#hide + +plot(inputs)#hide +``` + ### Imprecise Probabilities To represent purely epistemic uncertainty we provide the [`IntervalVariable`](@ref) type. @@ -73,6 +90,21 @@ Truncated p-boxes can be created by passing optional lower and upper bounds to t x = RandomVariable(pbox, :x) ``` +Examples of some parametric p-boxes are shown below. + +```@example rv +inputs = [ #hide + RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(-1, 2), :σ => Interval(0.5, 2))), :normal), #hide + RandomVariable(ProbabilityBox{Beta}(Dict(:α => Interval(1, 3), :β => Interval(2, 5))), :beta), #hide + RandomVariable(ProbabilityBox{Gamma}(Dict(:α => Interval(1, 3), :θ => Interval(0.5, 2))), :gamma), #hide + RandomVariable(ProbabilityBox{Exponential}(Dict(:θ => Interval(0.5, 2))), :exponential), #hide + RandomVariable(ProbabilityBox{Uniform}(Dict(:a => Interval(-1, 0), :b => Interval(1, 2))), :uniform), #hide + RandomVariable(ProbabilityBox{LogNormal}(Dict(:μ => Interval(-1, 1), :σ => Interval(0.5, 2))), :lognormal), #hide +]#hide + +plot(inputs)#hide + +``` ## Dependencies Modelling of dependent variables is exposed through the [`JointDistribution`](@ref). In its simplest form this is a wrapper around any implementation of the `MultivariateDistribution` interface defined by *Distributions.jl*. To construct a [`JointDistribution`](@ref) pass the `MultivariateDistribution` and a `vector{`Symbol`}` of matching dimension. diff --git a/src/UncertaintyQuantification.jl b/src/UncertaintyQuantification.jl index b77dc2b36..00b8700fb 100644 --- a/src/UncertaintyQuantification.jl +++ b/src/UncertaintyQuantification.jl @@ -22,6 +22,7 @@ using Random using Reexport using Roots using StatsBase +using RecipesBase @reexport using Distributions @@ -227,4 +228,11 @@ include("reliability/probabilityoffailure.jl") include("reliability/probabilityoffailure_imprecise.jl") include("sensitivity/sobolindices.jl") +include("plotting/plot_recipes.jl") + +include("util/fourier-transform.jl") +include("util/wrap.jl") +include("util/imprecise.jl") +include("util/kde.jl") + end diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl index 10ec31c34..88a146fdd 100644 --- a/src/inputs/imprecise/interval.jl +++ b/src/inputs/imprecise/interval.jl @@ -44,6 +44,11 @@ function bounds(i::Interval) return i.lb, i.ub end +hi(X::Interval) = X.ub +lo(X::Interval) = X.lb +hi(X::Real) = X +lo(X::Real) = X + """ IntervalVariable(lb::Real, ub::Real, name::Symbol) diff --git a/src/plotting/plot_recipes.jl b/src/plotting/plot_recipes.jl new file mode 100644 index 000000000..eeb3ac1d4 --- /dev/null +++ b/src/plotting/plot_recipes.jl @@ -0,0 +1,292 @@ +DEFAULT_ALPHA = 0.2 +DEFAULT_LABEL = "" +DEFAULT_GRID = false +DEFAULT_LEGEND = true +DEFAULT_CDF = false +DEFAULT_FILL_DISTRIBUTION=true +DEFAULT_FILL_IMPRECISE=true +DEFAULT_DISTRIBUTION = :pdf +DEFAULT_FILL = :gray +DEFAULT_COLOUR_PDF = :blue +DEFAULT_COLOUR_UPPER = :red +DEFAULT_COLOUR_LOWER = :black + +DEFAULT_INTERVAL_WIDTH=1.5 +DEFAULT_INTERVAL_EDGE_ALPHA=1 + +DEFAULT_DISTRIBUTION_WIDTH=2 + +DEFAULT_PLOT_RANGE_EXTEND_DENSITY = 0.2 +DEFAULT_PLOT_RANGE_EXTEND = 0.2 +DEFAULT_PLOT_RANGE_INTERVAL = 0.4 +DEFAULT_PLOT_GRID_NUMBER = 500 +DEFAULT_FONT_SIZE = 18 +DEFAULT_TICK_SIZE = 12 + +### +# Plots for UQInputs +### +@recipe function _plot( + x::RandomVariable{T}; cdf_on=DEFAULT_CDF, shade=DEFAULT_FILL_DISTRIBUTION +) where {T<:UnivariateDistribution} + grid --> DEFAULT_GRID + legend --> DEFAULT_LEGEND + label --> String(x.name) + (cdf_on ? ylabel --> "cdf" : ylabel --> "pdf") + seriescolor --> :auto + + lo = quantile(x, 0.001) + hi = quantile(x, 0.999) + w = hi - lo + lo -= abs(w * DEFAULT_PLOT_RANGE_EXTEND_DENSITY) + hi += abs(w * DEFAULT_PLOT_RANGE_EXTEND_DENSITY) + + xs = range(lo, hi, DEFAULT_PLOT_GRID_NUMBER) + ys = cdf_on ? cdf.(Ref(x), xs) : pdf.(Ref(x), xs) + + # Primary line: either cdf or pdf + @series begin + seriestype := :path + fill := nothing + alpha --> 1 + linewidth --> DEFAULT_DISTRIBUTION_WIDTH + label --> String(x.name) + xs, ys + end + + # Optional fill for PDF only: reuse color, don't advance palette + if !cdf_on && shade + @series begin + primary := false # <-- reuse the color from the primary line + seriestype := :path + fillrange := 0 # fill down to baseline + fillcolor := :match + fillalpha --> DEFAULT_ALPHA + linewidth --> 0 # fill only; no extra line + label := "" + xs, ys + end + end +end + +@recipe function _plot(x::IntervalVariable; shade=DEFAULT_FILL_IMPRECISE) + # --- plot-level defaults (soft) --- + grid --> DEFAULT_GRID + legend --> DEFAULT_LEGEND + ylabel --> "cdf" + label --> String(x.name) # single legend entry + + seriescolor --> :auto + + lo_grid = x.lb + hi_grid = x.ub + + width = hi_grid - lo_grid + + plot_lo = lo_grid - abs(width * DEFAULT_PLOT_RANGE_INTERVAL) # For adding a slight width to the left and right side of the interval plot + plot_hi = hi_grid + abs(width * DEFAULT_PLOT_RANGE_INTERVAL) + + x_grid = range(lo_grid, hi_grid, DEFAULT_PLOT_GRID_NUMBER) + + cdf_lo = x_grid .>= x.ub + cdf_hi = x_grid .> x.lb + + # Plot upper cdf (primary, inherits colour and label) + @series begin + seriestype := :path + alpha --> 1 + linewidth --> DEFAULT_DISTRIBUTION_WIDTH + x_grid, cdf_hi + end + + # Plot lower cdf + @series begin + primary := false + seriestype := :path + alpha --> 1 + label := "" + linewidth --> DEFAULT_DISTRIBUTION_WIDTH + x_grid, cdf_lo + end + + # Add invisible lines so interval bounds don't touch plot boundaries + @series begin + primary := false + seriestype := :path + alpha --> 0 + label := "" + linewidth --> 0 + range(plot_lo, plot_hi, 100), zeros(100) + end + + # Plot fill + if shade + @series begin + primary := false + seriestype := :path + fillcolor := :match # match this series' line color + fillrange := cdf_hi + color := DEFAULT_FILL + fillalpha := DEFAULT_ALPHA + linewidth --> 0 # draw fill only here + label := "" + x_grid, cdf_lo + end + end +end + +@recipe function _plot( + x::RandomVariable{T}; shade=DEFAULT_FILL_IMPRECISE +) where {T<:ProbabilityBox} + # --- plot-level defaults (soft) --- + grid --> DEFAULT_GRID + legend --> DEFAULT_LEGEND + ylabel --> "cdf" + label --> String(x.name) # single legend entry + + seriescolor --> :auto + + lo_grid = quantile(x, 0.001).lb + hi_grid = quantile(x, 0.999).ub + width = hi_grid - lo_grid + lo_grid = lo_grid - abs(width * DEFAULT_PLOT_RANGE_EXTEND) + hi_grid = hi_grid + abs(width * DEFAULT_PLOT_RANGE_EXTEND) + + x_grid = range(lo_grid, hi_grid, DEFAULT_PLOT_GRID_NUMBER) + cdf_evals = cdf.(Ref(x), x_grid) + + # Plot upper cdf (primary, inherits colour and label) + @series begin + seriestype := :path + alpha --> 1 + linewidth --> DEFAULT_DISTRIBUTION_WIDTH + x_grid, hi.(cdf_evals) + end + + # Plot lower cdf + @series begin + primary := false + seriestype := :path + alpha --> 1 + linewidth --> DEFAULT_DISTRIBUTION_WIDTH + label := "" + x_grid, lo.(cdf_evals) + end + + # Plot fill + if shade + @series begin + primary := false + seriestype := :path + fillcolor := :match # match this series' line color + fillrange := hi.(cdf_evals) + fillalpha --> DEFAULT_ALPHA + linewidth --> 0 # draw fill only here + label := "" + x_grid, lo.(cdf_evals) + end + end +end + +@recipe function _plot(x::Vector{T}) where {T<:UQInput} + # Filter out Parameter objects + x_no_params = filter(xi -> !isa(xi, Parameter), x) + + N = length(x_no_params) + cols = ceil(Int, sqrt(N)) + rows = ceil(Int, N / cols) + layout --> (rows, cols) + + # Choose a grid palette once (users can still override via plot(...; palette=...)) + # palette --> :default # or :default, :Dark2_8, etc. + + for i in 1:N + @series begin + subplot := i + seriescolor --> i # <-- panel i uses the i-th color from the current palette + x_no_params[i] + end + end +end + +### +# This code is a modified version of the plot recipe from IntervalArithmetic.jl +# https://github.com/JuliaIntervals/IntervalArithmetic.jl +### + +# Plot a 2D IntervalBox: +@recipe function _plot(x::Interval, y::Interval) + seriesalpha --> DEFAULT_ALPHA + seriestype := :shape + + label := false + + linecolor --> :black # Explicitly set edge color + linewidth --> DEFAULT_INTERVAL_WIDTH # Make edges more visible + linealpha --> DEFAULT_INTERVAL_EDGE_ALPHA + + x = [x.lb, x.ub, x.ub, x.lb] + y = [y.lb, y.lb, y.ub, y.ub] + + x, y +end + +# Plot a vector of 2D IntervalBoxes: +@recipe function _plot(xx::Vector{T}, yy::Vector{T}) where {T<:Interval} + seriesalpha --> DEFAULT_ALPHA + seriestype := :shape + + label := false + + linecolor := :black # Explicitly set edge color + linewidth --> DEFAULT_INTERVAL_WIDTH # Make edges more visible + linealpha --> DEFAULT_INTERVAL_EDGE_ALPHA + + xs = Float64[] + ys = Float64[] + + # build up coords: # (alternative: use @series) + for i in 1:length(xx) + (x, y) = (xx[i], yy[i]) + + # use NaNs to separate + append!(xs, [x.lb, x.ub, x.ub, x.lb, NaN]) + append!(ys, [y.lb, y.lb, y.ub, y.ub, NaN]) + end + + xs, ys +end + +### +# Plots for samples of data frames +### + +@recipe function _plot(x::Vector{Interval}) + if length(unique(x))==1 + return x[1] + else + grid --> DEFAULT_GRID + legend --> DEFAULT_LEGEND + + # xlabel --> x[1].name + ylabel --> "cdf" + N_samples = length(x) + + lows = sort(lo.(x)) + his = sort(hi.(x)) + + is = range(0, 1, length=N_samples) + + @series begin + seriestype := :steppre + color --> DEFAULT_COLOUR_LOWER + lows, is + end + + @series begin + seriestype := :steppost + color --> DEFAULT_COLOUR_UPPER + his, is + end + end +end diff --git a/test/Project.toml b/test/Project.toml index 37cc67ced..88257ec0d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,6 +5,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" diff --git a/test/inputs/imprecise/interval.jl b/test/inputs/imprecise/interval.jl index 54d202a10..6c0fe10a9 100644 --- a/test/inputs/imprecise/interval.jl +++ b/test/inputs/imprecise/interval.jl @@ -1,3 +1,5 @@ +using UncertaintyQuantification: hi, lo + @testset "Interval" begin name = :l lb = 0.14 @@ -16,6 +18,12 @@ @test !(0.17 ∈ interval) @test sprint(show, interval) == "[0.14, 0.16]" + + @test hi(interval) == interval.ub + @test lo(interval) == interval.lb + + @test hi(2.0) == 2.0 + @test lo(2.0) == 2.0 end @testset "IntervalVariable" begin diff --git a/test/plotting/plotting.jl b/test/plotting/plotting.jl new file mode 100644 index 000000000..9d456e0a6 --- /dev/null +++ b/test/plotting/plotting.jl @@ -0,0 +1,57 @@ +@testset "Plotting Recipes" begin + + @testset "RandomVariable PDF Plot" begin + rv = RandomVariable(Normal(0, 1), :X) + plt = plot(rv; cdf_on=false) + @test typeof(plt) <: Plots.Plot + @test plt.series_list[1][:seriestype] == :path + @test plt.series_list[2][:seriestype] == :path + end + + @testset "RandomVariable CDF Plot" begin + rv = RandomVariable(Normal(0, 1), :X) + plt = plot(rv; cdf_on=true) + @test typeof(plt) <: Plots.Plot + end + + @testset "IntervalVariable Plot" begin + iv = IntervalVariable(1.0, 2.0, :Y) + plt = plot(iv) + @test typeof(plt) <: Plots.Plot + end + + @testset "ProbabilityBox Plot" begin + pb = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(-1, 2), :σ => 2)), :X2) + plt = plot(pb) + @test typeof(plt) <: Plots.Plot + end + + @testset "Vector of UQInputs Plot" begin + inputs = [RandomVariable(Normal(0,1), :A), IntervalVariable(1.0, 2.0, :B)] + plt = plot(inputs) + @test typeof(plt) <: Plots.Plot + @test plt.layout isa Plots.GridLayout + end + + @testset "2D IntervalBox Plot" begin + x = Interval(1.0, 2.0) + y = Interval(3.0, 4.0) + plt = plot(x, y) + @test typeof(plt) <: Plots.Plot + @test plt.series_list[1][:seriestype] == :shape + end + + @testset "Vector of IntervalBoxes Plot" begin + xs = [Interval(1.0, 2.0), Interval(2.0, 3.0)] + ys = [Interval(3.0, 4.0), Interval(4.0, 5.0)] + plt = plot(xs, ys) + @test typeof(plt) <: Plots.Plot + @test plt.series_list[1][:seriestype] == :shape + end + + @testset "Vector of Intervals Plot" begin + intervals = [Interval(1.0, 2.0), Interval(1.5, 2.5), Interval(2.0, 3.0)] + plt = plot(intervals) + @test typeof(plt) <: Plots.Plot + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 0ba0f382d..f1f11e5cf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ using QuasiMonteCarlo using Random using StatsBase: fit, Histogram, corkendall using Test +using Plots using UncertaintyQuantification include("inputs/empiricaldistribution.jl") @@ -49,6 +50,8 @@ include("simulations/subset.jl") include("util/fourier-transform.jl") +include("plotting/plotting.jl") + if Sys.islinux() HPC = false HPC_account = "HPC_account_1"