Merge pull request #48 from JuliaOcean/tansformvariables

Use TransformVariables + without unit conversion

Project.toml

@@ -1,7 +1,7 @@
 name = "AIBECS"
 uuid = "ace601d6-714c-11e9-04e5-89b7fad23838"
 authors = ["Benoit Pasquier <briochemc@gmail.com>"]
-version = "0.7.12"
+version = "0.7.13"
 
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
@@ -26,6 +26,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+TransformVariables = "84d833dd-6860-57f9-a1a7-6da5db126cff"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulRecipes = "42071c24-d89e-48dd-8a24-8a12d9b8861f"
@@ -48,6 +49,7 @@ OceanGrids = "0.3"
 RecipesBase = "1"
 Reexport = "0.2"
 Shapefile = "0.6"
+TransformVariables = "0.3"
 UnPack = "1"
 Unitful = "1"
 UnitfulRecipes = "0.2, 1.0"

src/AIBECS.jl

@@ -24,6 +24,7 @@ using Distances
 using NearestNeighbors
 @reexport using MetadataArrays
 using Shapefile
+using TransformVariables
 
 
 

src/Parameters.jl

@@ -399,37 +399,37 @@ If the prior is `Uniform`, then the change of variables is the logit function.
 Otherwise, it's `identity`.
 """
 function subfun(::Type{T}) where {T<:AbstractParameters}
-    return λ -> reconstruct(T, [subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, flattenable_symbols(T))])
+    ks = flattenable_symbols(T)
+    λ2p(λ) = (subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, ks))
+    return λ -> T(; zip(ks, λ2p(λ))...)
 end
 function ∇subfun(::Type{T}) where {T<:AbstractParameters}
-    return λ -> reconstruct(T, [∇subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, flattenable_symbols(T))])
+    ks = flattenable_symbols(T)
+    ∇λ2p(λ) = (∇subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, ks))
+    return λ -> T(; zip(ks, ∇λ2p(λ))...)
 end
 function ∇²subfun(::Type{T}) where {T<:AbstractParameters}
-    return λ -> reconstruct(T, [∇²subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, flattenable_symbols(T))])
+    ks = flattenable_symbols(T)
+    ∇²λ2p(λ) = (∇²subfun(T, s)(λᵢ) for (λᵢ,s) in zip(λ, ks))
+    return λ -> T(; zip(ks, ∇²λ2p(λ))...)
 end
 function invsubfun(::Type{T}) where {T<:AbstractParameters}
-    return p -> [invsubfun(T, s)(pᵢ) for (pᵢ,s) in zip(vec(p), flattenable_symbols(T))]
+    return p -> [invsubfun(T, s)(pᵢ) for (pᵢ,s) in zip(flattenable_values(p), flattenable_symbols(T))]
 end
 # substitution function (change of variables) is determined from prior distribution
 subfun(::Type{T}, s::Symbol) where {T<:AbstractParameters} = subfun(prior(T,s))
 ∇subfun(::Type{T}, s::Symbol) where {T<:AbstractParameters} = ∇subfun(prior(T,s))
 ∇²subfun(::Type{T}, s::Symbol) where {T<:AbstractParameters} = ∇²subfun(prior(T,s))
 invsubfun(::Type{T}, s::Symbol) where {T<:AbstractParameters} = invsubfun(prior(T,s))
-# Fallback rule for change of variables is identity
-subfun(::Distribution) = identity
-∇subfun(::Distribution) = x -> one(x)
-∇²subfun(::Distribution) = x -> zero(x)
-invsubfun(::Distribution) = identity
-# p = exp(λ) for LogNormal
-subfun(::LogNormal) = exp
-∇subfun(::LogNormal) = exp
-∇²subfun(::LogNormal) = exp
-invsubfun(::LogNormal) = log
-# p = logistic(λ) for Uniform
-subfun(d::Uniform) = λ -> d.a + (d.b - d.a) / (exp(-λ) + 1)
-∇subfun(d::Uniform) = λ -> (d.b - d.a) * exp(-λ) / (exp(-λ) + 1)^2
-∇²subfun(d::Uniform) = λ -> (d.a - d.b) * exp(-λ) / (exp(-λ) + 1)^2 + 2(d.b - d.a) * exp(-2λ) / (exp(-λ) + 1)^3
-invsubfun(d::Uniform) = p -> -log((d.b - d.a) / (p - d.a) - 1)
+# using TransformVariables
+lb_domain(d::ContinuousUnivariateDistribution) = (lb = support(d).lb; lb == -Inf ? TransformVariables.Infinity{false}() : lb)
+ub_domain(d::ContinuousUnivariateDistribution) = (ub = support(d).ub; ub == Inf ? TransformVariables.Infinity{true}() : ub)
+import TransformVariables: transform
+transform(d::ContinuousUnivariateDistribution) = as(Real, lb_domain(d), ub_domain(d))
+subfun(d::ContinuousUnivariateDistribution) = x -> transform(d)(x)
+∇subfun(d::ContinuousUnivariateDistribution) = x -> ForwardDiff.derivative(subfun(d), x)
+∇²subfun(d::ContinuousUnivariateDistribution) = x -> ForwardDiff.derivative(∇subfun(d), x)
+invsubfun(d::ContinuousUnivariateDistribution) = x -> inverse(transform(d))(x)
 
 export subfun, ∇subfun, ∇²subfun, invsubfun
 

test/parameters.jl

@@ -18,13 +18,12 @@ const ∞ = Inf
 end
 function prior(::Type{T}, s::Symbol) where {T<:AbstractParameters}
     if flattenable(T, s)
-        U = units(T, s)
         if limits(T, s) == (0,∞)
-            μ = log(ustrip(upreferred(initial_value(T, s) * U)))
+            μ = log(initial_value(T, s))
             return LogNormal(μ, 1.0)
         elseif limits(T, s) == (-∞,∞)
-            μ = ustrip(upreferred(initial_value(T, s) * U))
-            σ = ustrip(upreferred(10.0U)) # Assumes that a sensible unit is chosen!
+            μ = initial_value(T, s)
+            σ = 10.0 # Assumes that a sensible unit is chosen (i.e., that within 10.0 * U)
             return Normal(μ, σ)
         elseif limits(T, s) == (0,1)
             return Uniform(0,1)