auto conversion

Project.toml

@@ -30,9 +30,10 @@ julia = "1.0"
 [extras]
 ControlSystems = "a6e380b2-a6ca-5380-bf3e-84a91bcd477e"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SLEEFPirates = "476501e8-09a2-5ece-8869-fb82de89a1fa"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ControlSystems", "ForwardDiff", "SLEEFPirates", "Test", "Plots"]
+test = ["ControlSystems", "ForwardDiff", "Measurements", "SLEEFPirates", "Test", "Plots"]

docs/src/index.md

@@ -10,7 +10,7 @@ This package facilitates working with probability distributions by means of Mont
 
 Although several interesting use cases for doing calculations with probability distributions have popped up (see [Examples](https://github.com/baggepinnen/MonteCarloMeasurements.jl#examples-1)), the original goal of the package is similar to that of [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl), to propagate the uncertainty from input of a function to the output. The difference compared to a `Measurement` is that `Particles` represent the distribution using a vector of unweighted particles, and can thus represent arbitrary distributions and handle nonlinear uncertainty propagation well. Functions like `f(x) = x²`, `f(x) = sign(x)` at `x=0` and long-time integration, are examples that are not handled well using linear uncertainty propagation ala [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl). MonteCarloMeasurements also support correlations between quantities.
 
-A number of type `Particles` behaves just as any other `Number` while partaking in calculations. After a calculation, an approximation to the **complete distribution** of the output is captured and represented by the output particles. `mean`, `std` etc. can be extracted from the particles using the corresponding functions. `Particles` also interact with [Distributions.jl](https://github.com/JuliaStats/Distributions.jl), so that you can call, e.g., `Normal(p)` and get back a `Normal` type from distributions or `fit(Gamma, p)` to get a `Gamma`distribution. Particles can also be iterated, asked for `maximum/minimum`, `quantile` etc. If particles are plotted with `plot(p)`, a histogram is displayed. This requires Plots.jl. A kernel-density estimate can be obtained by `density(p)` is StatsPlots.jl is loaded.
+A number of type `Particles` behaves just as any other `Number` while partaking in calculations. Particles also behave like a distribution, so after a calculation, an approximation to the **complete distribution** of the output is captured and represented by the output particles. `mean`, `std` etc. can be extracted from the particles using the corresponding functions. `Particles` also interact with [Distributions.jl](https://github.com/JuliaStats/Distributions.jl), so that you can call, e.g., `Normal(p)` and get back a `Normal` type from distributions or `fit(Gamma, p)` to get a `Gamma`distribution. Particles can also be iterated, asked for `maximum/minimum`, `quantile` etc. If particles are plotted with `plot(p)`, a histogram is displayed. This requires Plots.jl. A kernel-density estimate can be obtained by `density(p)` is StatsPlots.jl is loaded. A `Measurements.Measurements` can be converted to particles by calling the `Particles` constructor.
 
 Below, we show an example where an input uncertainty is propagated through `σ(x)`
 
@@ -100,7 +100,7 @@ The most basic constructor of [`Particles`](@ref) acts more or less like `randn(
 One can also call (`Particles/StaticParticles`)
 - `Particles(v::Vector)` pre-sampled particles
 - `Particles(N = 500, d::Distribution = Normal(0,1))` samples `N` particles from the distribution `d`.
-- The [`±`](@ref) operator (`\pm`) (similar to [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl)). We have `μ ± σ = μ + σ*Particles(DEFAUL_NUM_PARTICLES)`, where the global constant `DEFAUL_NUM_PARTICLES = 500`. You can change this if you would like, or simply define your own `±` operator like `±(μ,σ) = μ + σ*Particles(my_default_number, my_default_distribution)`. The upside-down operator [`∓`](@ref) (`\mp`) instead creates a `StaticParticles(100)`.
+- The [`±`](@ref) operator (`\pm`) (similar to [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl)). We have `μ ± σ = μ + σ*Particles(DEFAULT_NUM_PARTICLES)`, where the global constant `DEFAULT_NUM_PARTICLES = 500`. You can change this if you would like, or simply define your own `±` operator like `±(μ,σ) = μ + σ*Particles(my_default_number, my_default_distribution)`. The upside-down operator [`∓`](@ref) (`\mp`) instead creates a `StaticParticles(100)`.
 - The `..` binary infix operator creates uniformly sampled particles, e.g., `2..3 = Particles(Uniform(2,3))`
 
 **Common univariate distributions are sampled systematically**, meaning that a single random number is drawn and used to seed the sample. This will reduce the variance of the sample. If this is not desired, call `Particles(N, [d]; systematic=false)` The systematic sample can maintain its originally sorted order by calling `Particles(N, permute=false)`, but the default is to permute the sample so as to not have different `Particles` correlate strongly with each other.

src/MonteCarloMeasurements.jl

@@ -6,8 +6,8 @@ import Base: add_sum
 using Distributions, StatsBase, Requires
 
 
-const DEFAUL_NUM_PARTICLES = 500
-const DEFAUL_STATIC_NUM_PARTICLES = 100
+const DEFAULT_NUM_PARTICLES = 500
+const DEFAULT_STATIC_NUM_PARTICLES = 100
 
 const COMPARISON_FUNCTION = Ref{Function}(mean)
 const USE_UNSAFE_COMPARIONS = Ref(false)
@@ -93,6 +93,7 @@ Base.:\(x::AbstractVecOrMat{<:AbstractParticles}, y::AbstractVecOrMat{<:Abstract
 function __init__()
     @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")
     @require SLEEFPirates="476501e8-09a2-5ece-8869-fb82de89a1fa" include("sleefpirates.jl")
+    @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" include("measurements.jl")
 end
 
 end

src/measurements.jl

@@ -0,0 +1,16 @@
+
+for PT in (:Particles, :StaticParticles)
+    @eval begin
+        function $PT(N, m::Measurements.Measurement{T})::$PT{T,N} where T
+            $PT(N, Normal(Measurements.value(m),Measurements.uncertainty(m)))
+        end
+    end
+end
+
+function Particles(m::Measurements.Measurement{T}) where T
+    Particles(DEFAULT_NUM_PARTICLES, m)
+end
+
+function StaticParticles(m::Measurements.Measurement{T}) where T
+    StaticParticles(DEFAULT_STATIC_NUM_PARTICLES, m)
+end

src/particles.jl

@@ -1,7 +1,7 @@
 """
     μ ± σ
 
-Creates $DEFAUL_NUM_PARTICLES `Particles` with mean `μ` and std `σ`.
+Creates $DEFAULT_NUM_PARTICLES `Particles` with mean `μ` and std `σ`.
 If `μ` is a vector, the constructor `MvNormal` is used, and `σ` is thus treated as std if it's a scalar, and variances if it's a matrix or vector.
 See also [`∓`](@ref), [`..`](@ref)
 """
@@ -10,25 +10,25 @@ See also [`∓`](@ref), [`..`](@ref)
 """
     μ ∓ σ
 
-Creates $DEFAUL_STATIC_NUM_PARTICLES `StaticParticles` with mean `μ` and std `σ`.
+Creates $DEFAULT_STATIC_NUM_PARTICLES `StaticParticles` with mean `μ` and std `σ`.
 If `μ` is a vector, the constructor `MvNormal` is used, and `σ` is thus treated as std if it's a scalar, and variances if it's a matrix or vector.
 See also [`±`](@ref), [`⊗`](@ref)
 """
 ∓
 
 
-±(μ::Real,σ) = Particles{promote_type(float(typeof(μ)),float(typeof(σ))),DEFAUL_NUM_PARTICLES}(systematic_sample(DEFAUL_NUM_PARTICLES,Normal(μ,σ); permute=true))
-±(μ::AbstractVector,σ) = Particles(DEFAUL_NUM_PARTICLES, MvNormal(μ, σ))
-∓(μ::Real,σ) = StaticParticles{promote_type(float(typeof(μ)),float(typeof(σ))),DEFAUL_STATIC_NUM_PARTICLES}(systematic_sample(DEFAUL_STATIC_NUM_PARTICLES,Normal(μ,σ); permute=true))
-∓(μ::AbstractVector,σ) = StaticParticles(DEFAUL_STATIC_NUM_PARTICLES, MvNormal(μ, σ))
+±(μ::Real,σ) = Particles{promote_type(float(typeof(μ)),float(typeof(σ))),DEFAULT_NUM_PARTICLES}(systematic_sample(DEFAULT_NUM_PARTICLES,Normal(μ,σ); permute=true))
+±(μ::AbstractVector,σ) = Particles(DEFAULT_NUM_PARTICLES, MvNormal(μ, σ))
+∓(μ::Real,σ) = StaticParticles{promote_type(float(typeof(μ)),float(typeof(σ))),DEFAULT_STATIC_NUM_PARTICLES}(systematic_sample(DEFAULT_STATIC_NUM_PARTICLES,Normal(μ,σ); permute=true))
+∓(μ::AbstractVector,σ) = StaticParticles(DEFAULT_STATIC_NUM_PARTICLES, MvNormal(μ, σ))
 
 """
     a .. b
 
-Creates $DEFAUL_NUM_PARTICLES `Particles` with a `Uniform` distribution between `a` and `b`.
+Creates $DEFAULT_NUM_PARTICLES `Particles` with a `Uniform` distribution between `a` and `b`.
 See also [`±`](@ref), [`⊗`](@ref)
 """
-(..)(a,b) = Particles(DEFAUL_NUM_PARTICLES, Uniform(a,b))
+(..)(a,b) = Particles(DEFAULT_NUM_PARTICLES, Uniform(a,b))
 
 """
     ⊗(μ,σ) = outer_product(Normal.(μ,σ))
@@ -57,7 +57,7 @@ function outer_product(dists::AbstractVector{<:Distribution}, N=100_000)
     end
 end
 
-# StaticParticles(N::Integer = DEFAUL_NUM_PARTICLES; permute=true) = StaticParticles{Float64,N}(SVector{N,Float64}(systematic_sample(N, permute=permute)))
+# StaticParticles(N::Integer = DEFAULT_NUM_PARTICLES; permute=true) = StaticParticles{Float64,N}(SVector{N,Float64}(systematic_sample(N, permute=permute)))
 
 
 function print_functions_to_extend()

src/types.jl

@@ -56,7 +56,7 @@ for PT in (:Particles, :StaticParticles)
 
         $PT{T,N}(p::$PT{T,N}) where {T,N} = p
 
-        function $PT(N::Integer=DEFAUL_NUM_PARTICLES, d::Distribution{<:Any,VS}=Normal(0,1); permute=true, systematic=VS==Continuous) where VS
+        function $PT(N::Integer=DEFAULT_NUM_PARTICLES, d::Distribution{<:Any,VS}=Normal(0,1); permute=true, systematic=VS==Continuous) where VS
             if systematic
                 v = systematic_sample(N,d; permute=permute)
             else
@@ -77,11 +77,11 @@ for PT in (:Particles, :StaticParticles)
 end
 
 function Particles(d::Distribution;kwargs...)
-    Particles(DEFAUL_NUM_PARTICLES, d; kwargs...)
+    Particles(DEFAULT_NUM_PARTICLES, d; kwargs...)
 end
 
 function StaticParticles(d::Distribution;kwargs...)
-    StaticParticles(DEFAUL_STATIC_NUM_PARTICLES, d; kwargs...)
+    StaticParticles(DEFAULT_STATIC_NUM_PARTICLES, d; kwargs...)
 end
 
 

test/runtests.jl

@@ -9,6 +9,8 @@ import Plots
 Random.seed!(0)
 
 @testset "MonteCarloMeasurements.jl" begin
+    @info "Testing MonteCarloMeasurements"
+
     # σ/√N = σm
     @time @testset "sampling" begin
         @info "Testing sampling"
@@ -526,6 +528,7 @@ Random.seed!(0)
     include("test_forwarddiff.jl")
     include("test_deconstruct.jl")
     include("test_sleefpirates.jl")
+    include("test_measurements.jl")
 
 end
 

test/test_deconstruct.jl

@@ -7,6 +7,7 @@ ControlSystems.TransferFunction(matrix::Array{<:ControlSystems.SisoRational,2},
 
 
 @testset "deconstruct" begin
+    @info "Testing deconstruct"
     unsafe_comparisons()
     N = 50
     P = tf(1 +0.1StaticParticles(N), [1, 1+0.1StaticParticles(N)])

test/test_forwarddiff.jl

@@ -2,6 +2,7 @@ using MonteCarloMeasurements, ForwardDiff, Test
 const FD = ForwardDiff
 
 @testset "forwarddiff" begin
+    @info "Testing ForwardDiff"
     c = 1 ± 0.1 # These are the uncertain parameters
     d = 1 ± 0.1 # These are the uncertain parameters
     # In the cost function below, we ensure that $cx+dy > 10 \; ∀ \; c,d ∈ P$ by looking at the worst case

test/test_measurements.jl

@@ -0,0 +1,23 @@
+
+@testset "Measurements" begin
+    @info "Testing Measurements"
+    import Measurements
+
+    m = Measurements.measurement(1,2)
+    @test Particles(m) isa Particles{Float64, MonteCarloMeasurements.DEFAULT_NUM_PARTICLES}
+    @test StaticParticles(m) isa StaticParticles{Float64, MonteCarloMeasurements.DEFAULT_STATIC_NUM_PARTICLES}
+
+
+    @test Particles(10, m) isa Particles{Float64, 10}
+    @test StaticParticles(10, m) isa StaticParticles{Float64, 10}
+
+
+
+    @test Particles.([m, m]) isa Vector{Particles{Float64, MonteCarloMeasurements.DEFAULT_NUM_PARTICLES}}
+    @test StaticParticles.([m, m]) isa Vector{StaticParticles{Float64, MonteCarloMeasurements.DEFAULT_STATIC_NUM_PARTICLES}}
+
+    @test Particles(m) ≈ 1 ± 2
+    @test std(Particles(m)) ≈ 2 atol=1e-3
+    @test mean(Particles(m)) ≈ 1 atol=1e-3
+
+end

test/test_sleefpirates.jl

@@ -1,11 +1,16 @@
-using SLEEFPirates
-a = rand(500)
-@test map(exp, a) ≈ exp(Particles(a)).particles
-@test map(log, a) ≈ log(Particles(a)).particles
-@test map(cos, a) ≈ cos(Particles(a)).particles
+@testset "SLEEFPirates" begin
+    @info "Testing SLEEFPirates"
 
+    using SLEEFPirates
+    a = rand(500)
+    @test map(exp, a) ≈ exp(Particles(a)).particles
+    @test map(log, a) ≈ log(Particles(a)).particles
+    @test map(cos, a) ≈ cos(Particles(a)).particles
 
-a = rand(Float32, 500)
-@test map(exp, a) ≈ exp(Particles(a)).particles
-@test map(log, a) ≈ log(Particles(a)).particles
-@test map(cos, a) ≈ cos(Particles(a)).particles
+
+    a = rand(Float32, 500)
+    @test map(exp, a) ≈ exp(Particles(a)).particles
+    @test map(log, a) ≈ log(Particles(a)).particles
+    @test map(cos, a) ≈ cos(Particles(a)).particles
+
+end