Merge pull request #3 from bonStats:dual-error

Dual-error

Project.toml

@@ -1,7 +1,7 @@
 name = "BayesScoreCal"
 uuid = "986b08a2-df14-4163-beea-1d3f6d344e0b"
 authors = ["Joshua J Bon"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"

examples/corr-ou-vi/corr-ou-test.jl

@@ -109,7 +109,7 @@ function testfun(ou::MvNormal, T::Float64, N::Int64, N_importance::Int64, N_samp
     advi = ADVI(10, 1000)
     q_approx = vi(mod, advi)
     approx_samples = vrand(q_approx, N_samples)
-    tr_approx_samples = inv(q_approx.transform).(approx_samples)
+    tr_approx_samples = inverse(q_approx.transform).(approx_samples)
 
     # calibration posterior (scale by vmultiplier on underlying transformed space)
     q_cal = transformed(
@@ -119,7 +119,7 @@ function testfun(ou::MvNormal, T::Float64, N::Int64, N_importance::Int64, N_samp
 
     # sample calibration points
     cal_points = vrand(q_cal, N_importance)
-    tr_cal_points = inv(q_cal.transform).(cal_points)
+    tr_cal_points = inverse(q_cal.transform).(cal_points)
 
     # log prior evaluated on importance distribution samples
     ℓprior = [logprior(mod, par2tuple(v)) for (i, v) in enumerate(cal_points)]

examples/lotka-volterra-sde/lotka-volterra-test.jl

@@ -143,7 +143,7 @@ select_approx_samples = getsamples(ch_approx, :pars, N_importance)
 
 # sample calibration points
 if !use_cache
-    cal_points = inv(bij).(multiplyscale(bij.(select_approx_samples), vmultiplier))
+    cal_points = inverse(bij).(multiplyscale(bij.(select_approx_samples), vmultiplier))
 
     # newx approx models: pre-allocate
     tr_approx_samples_newx = SharedArray{Float64}(length(cal_points[1]), N_energy, N_importance)
@@ -190,7 +190,7 @@ cal = Calibration(tr_cal_points, tr_approx_samples_newx)
 
 d = BayesScoreCal.dimension(cal)[1]
 tf = CholeskyAffine(d)
-M = inv(Diagonal(std(cal.μs)))
+M = inverse(Diagonal(std(cal.μs)))
 res = energyscorecalibrate!(tf, cal, is_weights, scaling = M, penalty = (0.0, 0.05))
 
 calcheck_approx = coverage(cal, checkprobs)
@@ -202,7 +202,7 @@ rmse(cal,tf)
 
 approx_samples = vec(getsamples(ch_approx, :pars))
 tr_approx_samples = bij.(approx_samples)
-tf_samples = inv(bij).(tf.(tr_approx_samples, [mean(tr_approx_samples)]))
+tf_samples = inverse(bij).(tf.(tr_approx_samples, [mean(tr_approx_samples)]))
 
 [mean(tf_samples) std(tf_samples)]
 

src/energy-score.jl

@@ -59,11 +59,10 @@ function energyscorecalibrate!(tf::T, cal::Calibration, weights::Vector{Float64}
     perm = PermuteVector(nsamples(cal))
     initialx = paramvec(tf)
 
-    f = OnceDifferentiable((x) -> negenergyscore(update!(tf,x), cal, weights, perm, β, scaling, penalty), initialx; autodiff = :forward)
+    f = OnceDifferentiable((x) -> negenergyscore(maketransform(tf,x), cal, weights, perm, β, scaling, penalty), initialx; autodiff = :forward)
 
     res = Optim.optimize(f, initialx, BFGS(), options)
 
-    # neccesary?
     update!(tf, Optim.minimizer(res))
 
     return res

src/transforms/cholesky-affine.jl

@@ -1,42 +1,40 @@
 include("pd-chol-helper.jl")
 
-mutable struct CholeskyAffine{N} <: Transform{N}
-    L::LowerTriangular
-    b::Vector{<:Real}
-    function CholeskyAffine(L::LowerTriangular, b::Vector{<:Real})
+mutable struct CholeskyAffine{N,T<:Real} <: Transform{N}
+    L::LowerTriangular{T, Matrix{T}}
+    b::Vector{T}
+    function CholeskyAffine(L::LowerTriangular{T, Matrix{T}}, b::Vector{T}) where {T<:Real}
         @assert length(b) == size(L, 1)
-        new{length(b)}(L, b)
+        new{length(b),T}(L, b)
     end
 end
 
-CholeskyAffine(vL::Vector{<:Real}, b::Vector{<:Real}) = CholeskyAffine(vec2chol(vL), b)
+CholeskyAffine(vL::Vector{T}, b::Vector{T}) where {T<:Real}  = CholeskyAffine(vec2chol(vL), b)
 CholeskyAffine(vLb::Vector{<:Real}, d::Int64) = CholeskyAffine(vLb[1:(end-d)], vLb[(end-d+1):end])
 function CholeskyAffine(d::Int64)
     par_zero = zeros(Int((d^2 + d)/2 + d))
     CholeskyAffine(par_zero, d)
 end
 
-dimension(::CholeskyAffine{N}) where {N} = (N,) # object dimension
-length(::CholeskyAffine{N}) where {N} = N # for internal use (within CholeskyAffine functions)
-nparam(::CholeskyAffine{N}) where {N} = N*(N + 3)/2 # number of params
+dimension(::CholeskyAffine{N,T}) where {N,T} = (N,) # object dimension
+length(::CholeskyAffine{N,T}) where {N,T} = N # for internal use (within CholeskyAffine functions)
+nparam(::CholeskyAffine{N,T}) where {N,T} = N*(N + 3)/2 # number of params
 
-function update!(chaf::CholeskyAffine, vL::Vector{<:Real}, b::Vector{<:Real})
-    chaf.L = vec2chol(vL)
-    chaf.b = b
-    return chaf
+# for update! CholeskyAffine may not have same type as other arguments (e.g. Float64 -> Dual)
+function update!(chaf::CholeskyAffine{N,T}, vL::Vector{T}, b::Vector{T}) where {N,T<:Real} 
+    chaf.L .= vec2chol(vL)
+    chaf.b .= b
 end
 
-function update!(chaf::CholeskyAffine, vLb::Vector{<:Real})
-    n = length(chaf)
-    chaf.L = vec2chol(vLb[1:(end-n)])
-    chaf.b = vLb[(end-n+1):end]
-    return chaf
-end
+update!(chaf::CholeskyAffine{N,T}, vLb::Vector{T}) where {N,T<:Real} = update!(chaf, vLb[1:(end-N)], vLb[(end-N+1):end])
+
+maketransform(::CholeskyAffine{N,Ti}, vLb::Vector{To}) where {N,Ti<:Real,To<:Real} = CholeskyAffine(vLb, N)
 
 scaleparamvec(tf::CholeskyAffine) = chol2vec(tf.L)
 biasparamvec(tf::CholeskyAffine) = tf.b
 paramvec(tf::CholeskyAffine) = [chol2vec(tf.L); tf.b]
 
+# for  (::CholeskyAffine) may not have same type as input (e.g. Float64 -> Dual)
 (chaf::CholeskyAffine)(x::T, μs::T) where {T} = chaf.L * (x - μs) + μs + chaf.b 
 
 (chaf::CholeskyAffine)(cal::Calibration{T}) where {T} = Calibration(cal.values, hcat([chaf.(clm, [cal.μs[i]]) for (i, clm) in enumerate(eachcol(cal.samples))]...))

src/transforms/pd-chol-helper.jl

@@ -15,9 +15,9 @@ expdiag(v::Real, i::Int64, j::Int64) = i==j ? exp(v) : v
 
 Convert a unconstrained vector v to lower-triangular cholesky decomposition of (n x n) positive definite matrix.
 """
-function vec2chol(v::Vector{<:Real}, n::Int64)
+function vec2chol(v::Vector{T}, n::Int64) where {T<:Real}
     k = 0
-    LowerTriangular([ j<=i ? (k+=1; expdiag(v[k], i , j)) : 0.0 for i=1:n, j=1:n ])
+    LowerTriangular([ j<=i ? (k+=1; expdiag(v[k], i , j)) : zero(T) for i=1:n, j=1:n ])
 end
 
 vec2chol(v::Vector{<:Real}) = vec2chol(v, cholsize(v))

src/transforms/uni-affine.jl

@@ -1,6 +1,10 @@
-mutable struct UnivariateAffine <: Transform{1}
-    σ::Real
-    b::Real
+mutable struct UnivariateAffine{T<:Real} <: Transform{1}
+    σ::T
+    b::T
+end
+
+function UnivariateAffine(v::Vector{T}) where {T<:Real}
+    UnivariateAffine(exp(v[1]),v[2])
 end
 
 function UnivariateAffine()
@@ -11,17 +15,15 @@ dimension(::UnivariateAffine) = () # object dimension
 length(::UnivariateAffine) = 1 # for internal use (within CholeskyAffine functions)
 nparam(::UnivariateAffine) = 2 # number of params
 
-function update!(uaf::UnivariateAffine, logσ::Real, b::Real)
+function update!(uaf::UnivariateAffine{T}, logσ::T, b::T) where {T<:Real}
     uaf.σ = exp(logσ)
     uaf.b = b
     return uaf
 end
 
-function update!(uaf::UnivariateAffine, v::Vector{<:Real})
-    uaf.σ = exp(v[1])
-    uaf.b = v[2]
-    return uaf
-end
+update!(uaf::UnivariateAffine{T}, v::Vector{T}) where {T<:Real} = update!(uaf, exp(v[1]), v[2])
+
+maketransform(::UnivariateAffine{Ti}, v::Vector{To}) where {Ti<:Real,To<:Real} = UnivariateAffine(v)
 
 scaleparamvec(tf::UnivariateAffine) = log(tf.σ)
 biasparamvec(tf::UnivariateAffine) = tf.b