From ded4523b9a5208721e059af8dc7804ac7f74d400 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 9 Sep 2022 11:17:13 -0400 Subject: [PATCH 01/11] Rename Julia interface, fix exports --- julia/Project.toml | 2 +- julia/src/{Bridgestan.jl => BridgeStan.jl} | 17 +++++--- julia/test/runtests.jl | 50 +++++++++++----------- 3 files changed, 36 insertions(+), 33 deletions(-) rename julia/src/{Bridgestan.jl => BridgeStan.jl} (96%) diff --git a/julia/Project.toml b/julia/Project.toml index d528ef73..cc7d7ce7 100644 --- a/julia/Project.toml +++ b/julia/Project.toml @@ -1,4 +1,4 @@ -name = "Bridgestan" +name = "BridgeStan" uuid = "c88b6f0a-829e-4b0b-94b7-f06ab5908f5a" authors = ["Brian Ward ", "Bob Carpenter "] version = "0.1.0" diff --git a/julia/src/Bridgestan.jl b/julia/src/BridgeStan.jl similarity index 96% rename from julia/src/Bridgestan.jl rename to julia/src/BridgeStan.jl index c6598b83..6a2704cf 100644 --- a/julia/src/Bridgestan.jl +++ b/julia/src/BridgeStan.jl @@ -1,15 +1,18 @@ -module Bridgestan +module BridgeStan export StanModel, - log_density_gradient!, - K, + name, param_num, - param_constrain!, - dims, param_unc_num, - param_unconstrain!, - destroy + param_names, + param_unc_names, + param_constrain, + param_unconstrain, + param_unconstrain_json, + log_density, + log_density_gradient, + log_density_hessian mutable struct StanModelStruct end diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 145b0e3c..7d156208 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -1,4 +1,4 @@ -using Bridgestan +using BridgeStan using Test @testset "bernoulli" begin @@ -12,30 +12,30 @@ using Test lib = joinpath(@__DIR__, "../../stan/bernoulli/bernoulli_model.so") data = joinpath(@__DIR__, "../../stan/bernoulli/bernoulli.data.json") - model = Bridgestan.StanModel(lib, data) + model = BridgeStan.StanModel(lib, data) - @test Bridgestan.name(model) == "bernoulli_model" + @test BridgeStan.name(model) == "bernoulli_model" y = [0, 1, 0, 0, 0, 0, 0, 0, 0, 1] R = 1000 for _ in 1:R - x = rand(Bridgestan.param_num(model)) + x = rand(BridgeStan.param_num(model)) q = @. log(x / (1 - x)) # unconstrained scale - (log_density, gradient) = Bridgestan.log_density_gradient(model, q, jacobian = 0) + (log_density, gradient) = BridgeStan.log_density_gradient(model, q, jacobian = 0) p = x[1] @test isapprox(log_density, bernoulli(y, p)) - constrained_parameters = Bridgestan.param_constrain(model, q) + constrained_parameters = BridgeStan.param_constrain(model, q) @test isapprox(constrained_parameters, x) - unconstrained_parameters= Bridgestan.param_unconstrain(model, constrained_parameters) + unconstrained_parameters= BridgeStan.param_unconstrain(model, constrained_parameters) @test isapprox(unconstrained_parameters, q) end - @test isapprox(Bridgestan.param_num(model), 1) - @test isapprox(Bridgestan.param_unc_num(model), 1) + @test isapprox(BridgeStan.param_num(model), 1) + @test isapprox(BridgeStan.param_unc_num(model), 1) end @@ -55,7 +55,7 @@ end data = joinpath(@__DIR__, "../../stan/multi/multi.data.json") nt = Threads.nthreads() - model = Bridgestan.StanModel(lib, data) + model = BridgeStan.StanModel(lib, data) R = 1000 ld = Vector{Bool}(undef, R) @@ -63,8 +63,8 @@ end @sync for it in 1:nt Threads.@spawn for r in it:nt:R - x = randn(Bridgestan.param_num(model)) - (lp, grad) = Bridgestan.log_density_gradient(model, x) + x = randn(BridgeStan.param_num(model)) + (lp, grad) = BridgeStan.log_density_gradient(model, x) ld[r] = isapprox(lp, gaussian(x)) g[r] = isapprox(grad, grad_gaussian(x)) @@ -83,20 +83,20 @@ end lib = joinpath(@__DIR__, "../../stan/gaussian/gaussian_model.so") data = joinpath(@__DIR__, "../../stan/gaussian/gaussian.data.json") - model = Bridgestan.StanModel(lib, data) + model = BridgeStan.StanModel(lib, data) theta = [0.2, 1.9] theta_unc = [0.2, log(1.9)] - theta_test = Bridgestan.param_constrain(model, theta_unc) + theta_test = BridgeStan.param_constrain(model, theta_unc) @test isapprox(theta, theta_test) - theta_unc_test = Bridgestan.param_unconstrain(model, theta) + theta_unc_test = BridgeStan.param_unconstrain(model, theta) @test isapprox(theta_unc, theta_unc_test) theta_json = "{\"mu\": 0.2, \"sigma\": 1.9}" - theta_unc_j_test = Bridgestan.param_unconstrain_json(model, theta_json) + theta_unc_j_test = BridgeStan.param_unconstrain_json(model, theta_json) @test isapprox(theta_unc, theta_unc_j_test) end @@ -117,25 +117,25 @@ end lib = joinpath(@__DIR__, "../../stan/fr_gaussian/fr_gaussian_model.so") data = joinpath(@__DIR__, "../../stan/fr_gaussian/fr_gaussian.data.json") - model = Bridgestan.StanModel(lib, data) + model = BridgeStan.StanModel(lib, data) size = 16 unc_size = 10 - @test isapprox(size, Bridgestan.param_num(model, include_tp=true, include_gq=true)) - @test isapprox(unc_size, Bridgestan.param_unc_num(model)) + @test isapprox(size, BridgeStan.param_num(model, include_tp=true, include_gq=true)) + @test isapprox(unc_size, BridgeStan.param_unc_num(model)) D = 4 a = randn(unc_size) - b = Bridgestan.param_constrain(model, a) + b = BridgeStan.param_constrain(model, a) B = reshape(b, (D,D)) B_expected = _covariance_constrain_transform(a, D) @test isapprox(B_expected, B) - c = Bridgestan.param_unconstrain(model, b) + c = BridgeStan.param_unconstrain(model, b) @test isapprox(a, c) - names = Bridgestan.param_names(model, include_tp=true, include_gq=true) + names = BridgeStan.param_names(model, include_tp=true, include_gq=true) name_eq = Vector{Bool}(undef, size) pos = 1 for j = 1:4 @@ -146,7 +146,7 @@ end end @test all(name_eq) - unc_names = Bridgestan.param_unc_names(model) + unc_names = BridgeStan.param_unc_names(model) name_unc_eq = Vector{Bool}(undef, unc_size) for n = 1:10 name_unc_eq[n] = unc_names[n] == ("Omega." * string(n)) @@ -159,11 +159,11 @@ end lib = joinpath(@__DIR__, "../../stan/simple/simple_model.so") data = joinpath(@__DIR__, "../../stan/simple/simple.data.json") - model = Bridgestan.StanModel(lib, data) + model = BridgeStan.StanModel(lib, data) D = 5 y = rand(D) - lp, grad, hess = Bridgestan.log_density_hessian(model, y) + lp, grad, hess = BridgeStan.log_density_hessian(model, y) @test isapprox(-y, grad) using LinearAlgebra From b940fc9522254d508b0fcaa93e532c15b7d4ace9 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 9 Sep 2022 11:18:46 -0400 Subject: [PATCH 02/11] Remove unnecessary types --- julia/src/BridgeStan.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index 6a2704cf..dd6f1696 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -90,7 +90,7 @@ end function param_constrain(sm::StanModel, theta_unc; include_tp=false, include_gq=false) - out::Vector{Float64} = zeros(param_num(sm, include_tp=include_tp, include_gq=include_gq)) + out = zeros(param_num(sm, include_tp=include_tp, include_gq=include_gq)) rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_constrain"), Cint, (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}), @@ -104,7 +104,7 @@ end function param_unconstrain(sm::StanModel, theta) - out::Vector{Float64} = zeros(param_unc_num(sm)) + out = zeros(param_unc_num(sm)) rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_unconstrain"), Cint, (Ptr{StanModelStruct}, Ref{Cdouble}, Ref{Cdouble}), @@ -117,7 +117,7 @@ function param_unconstrain(sm::StanModel, theta) end function param_unconstrain_json(sm::StanModel, theta::String) - out::Vector{Float64} = zeros(param_unc_num(sm)) + out = zeros(param_unc_num(sm)) rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_unconstrain_json"), Cint, (Ptr{StanModelStruct}, Cstring, Ref{Cdouble}), @@ -144,7 +144,7 @@ end function log_density_gradient(sm::StanModel, q; propto = true, jacobian = true) lp = Ref{Float64}(0.0) - grad::Vector{Float64} = zeros(param_unc_num(sm)) + grad = zeros(param_unc_num(sm)) rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density_gradient"), Cint, @@ -162,8 +162,8 @@ end function log_density_hessian(sm::StanModel, q; propto = true, jacobian = true) lp = Ref{Float64}(0.0) dims = param_unc_num(sm) - grad::Vector{Float64} = zeros(dims) - hess::Vector{Float64} = zeros(dims * dims) + grad = zeros(dims) + hess = zeros(dims * dims) rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density_hessian"), Cint, From fb51bf30f7e0ed0c91a2d275403e9fa0a5ad63df Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 9 Sep 2022 11:58:23 -0400 Subject: [PATCH 03/11] Add out! versions --- julia/src/BridgeStan.jl | 78 ++++++++++++++++++++++++++++++++--------- 1 file changed, 61 insertions(+), 17 deletions(-) diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index dd6f1696..ff06ffd9 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -7,6 +7,11 @@ export param_unc_num, param_names, param_unc_names, + param_constrain!, + param_unconstrain!, + param_unconstrain_json!, + log_density_gradient!, + log_density_hessian!, param_constrain, param_unconstrain, param_unconstrain_json, @@ -79,7 +84,6 @@ function param_names(sm::StanModel; include_tp = false, include_gq = false) [string(s) for s in split(unsafe_string(cstr), ',')] end - function param_unc_names(sm::StanModel) cstr = ccall(Libc.Libdl.dlsym(sm.lib, "param_unc_names"), Cstring, @@ -88,9 +92,11 @@ function param_unc_names(sm::StanModel) [string(s) for s in split(unsafe_string(cstr), ',')] end - -function param_constrain(sm::StanModel, theta_unc; include_tp=false, include_gq=false) - out = zeros(param_num(sm, include_tp=include_tp, include_gq=include_gq)) +function param_constrain!(sm::StanModel, theta_unc, out::Vector{Float64}; include_tp=false, include_gq=false) + dims = param_num(sm; include_tp=include_tp, include_gq=include_gq) + if length(out) != dims + error("out must be same size as number of constrained parameters") + end rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_constrain"), Cint, (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}), @@ -102,9 +108,18 @@ function param_constrain(sm::StanModel, theta_unc; include_tp=false, include_gq= end end +function param_constrain(sm::StanModel, theta_unc; include_tp=false, include_gq=false) + out = zeros(param_num(sm, include_tp=include_tp, include_gq=include_gq)) + param_constrain!(sm, theta_unc, out; include_tp=include_tp, include_gq=include_gq) +end + + +function param_unconstrain!(sm::StanModel, theta, out::Vector{Float64}) + dims = param_unc_num(sm) + if length(out) != dims + error("out must be same size as number of unconstrained parameters") + end -function param_unconstrain(sm::StanModel, theta) - out = zeros(param_unc_num(sm)) rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_unconstrain"), Cint, (Ptr{StanModelStruct}, Ref{Cdouble}, Ref{Cdouble}), @@ -116,8 +131,17 @@ function param_unconstrain(sm::StanModel, theta) end end -function param_unconstrain_json(sm::StanModel, theta::String) +function param_unconstrain(sm::StanModel, theta) out = zeros(param_unc_num(sm)) + param_unconstrain!(sm, theta, out) +end + +function param_unconstrain_json!(sm::StanModel, theta::String, out::Vector{Float64}) + dims = param_unc_num(sm) + if length(out) != dims + error("out must be same size as number of unconstrained parameters") + end + rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_unconstrain_json"), Cint, (Ptr{StanModelStruct}, Cstring, Ref{Cdouble}), @@ -129,6 +153,11 @@ function param_unconstrain_json(sm::StanModel, theta::String) end end +function param_unconstrain_json(sm::StanModel, theta::String) + out = zeros(param_unc_num(sm)) + param_unconstrain_json!(sm, theta, out) +end + function log_density(sm::StanModel, q; propto = true, jacobian = true) lp = Ref{Float64}(0.0) rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density"), @@ -142,39 +171,54 @@ function log_density(sm::StanModel, q; propto = true, jacobian = true) end end -function log_density_gradient(sm::StanModel, q; propto = true, jacobian = true) +function log_density_gradient!(sm::StanModel, q, out::Vector{Float64}; propto = true, jacobian = true) lp = Ref{Float64}(0.0) - grad = zeros(param_unc_num(sm)) + dims = param_unc_num(sm) + if length(out) != dims + error("out must be same size as number of unconstrained parameters") + end rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density_gradient"), Cint, (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), - sm.stanmodel, propto, jacobian, q, lp, grad) + sm.stanmodel, propto, jacobian, q, lp, out) if rc != 0 error("log_density_gradient failed on C++ side; see stderr for messages") else - (lp[], grad) + (lp[], out) end end +function log_density_gradient(sm::StanModel, q; propto = true, jacobian = true) + grad = zeros(param_unc_num(sm)) + log_density_gradient!(sm, q, grad; propto=propto, jacobian=jacobian) +end - -function log_density_hessian(sm::StanModel, q; propto = true, jacobian = true) +function log_density_hessian!(sm::StanModel, q, out_grad::Vector{Float64}, out_hess::Vector{Float64}; propto = true, jacobian = true) lp = Ref{Float64}(0.0) dims = param_unc_num(sm) - grad = zeros(dims) - hess = zeros(dims * dims) + if length(out_grad) != dims + error("out_grad must be same size as number of unconstrained parameters") + elseif length(out_hess) != dims*dims + error("out_hess must be same size as (number of unconstrained parameters)^2") + end rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density_hessian"), Cint, (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), - sm.stanmodel, propto, jacobian, q, lp, grad, hess) + sm.stanmodel, propto, jacobian, q, lp, out_grad, out_hess) if rc != 0 error("log_density_hessian failed on C++ side; see stderr for messages") else - (lp[], grad, reshape(hess, (dims, dims))) + (lp[], out_grad, reshape(out_hess, (dims, dims))) end end +function log_density_hessian(sm::StanModel, q; propto = true, jacobian = true) + dims = param_unc_num(sm) + grad = zeros(dims) + hess = zeros(dims * dims) + log_density_hessian!(sm,q, grad, hess; propto=propto, jacobian=jacobian) +end end From a7ac69d48061b21f94610b6b747655a819a64127 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 9 Sep 2022 15:06:21 -0400 Subject: [PATCH 04/11] Add extensive testing --- julia/src/BridgeStan.jl | 28 ++- julia/test/Manifest.toml | 9 +- julia/test/Project.toml | 1 + julia/test/runtests.jl | 349 ++++++++++++++++++++++++++++++--- python/test/test_bridgestan.py | 12 +- 5 files changed, 354 insertions(+), 45 deletions(-) diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index ff06ffd9..fb09a87d 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -29,16 +29,27 @@ mutable struct StanModel const seed::UInt32 const chain_id::UInt32 - function StanModel(stanlib_::String, datafile_::String, seed_ = 204, chain_id_ = 0) + function StanModel(stanlib_::String, datafile_::String = "", seed_ = 204, chain_id_ = 0) seed = convert(UInt32, seed_) chain_id = convert(UInt32, chain_id_) + + if !isfile(stanlib_) + throw(SystemError("Dynamic library file not found")) + end + + if datafile_ != "" && !isfile(datafile_) + throw(SystemError("Data file not found")) + end + lib = Libc.Libdl.dlopen(stanlib_) stanmodel = ccall(Libc.Libdl.dlsym(lib, "construct"), Ptr{StanModelStruct}, (Cstring, UInt32, UInt32), datafile_, seed, chain_id) - + if stanmodel == C_NULL + error("could not construct model RNG") + end sm = new(lib, stanmodel, datafile_, seed, chain_id) @@ -47,7 +58,6 @@ mutable struct StanModel UInt32, (Ptr{StanModelStruct},), sm.stanmodel) - @async Libc.Libdl.dlclose(sm.lib) end finalizer(f, sm) @@ -95,7 +105,7 @@ end function param_constrain!(sm::StanModel, theta_unc, out::Vector{Float64}; include_tp=false, include_gq=false) dims = param_num(sm; include_tp=include_tp, include_gq=include_gq) if length(out) != dims - error("out must be same size as number of constrained parameters") + throw(DimensionMismatch("out must be same size as number of constrained parameters")) end rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_constrain"), Cint, @@ -117,7 +127,7 @@ end function param_unconstrain!(sm::StanModel, theta, out::Vector{Float64}) dims = param_unc_num(sm) if length(out) != dims - error("out must be same size as number of unconstrained parameters") + throw(DimensionMismatch("out must be same size as number of unconstrained parameters")) end rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_unconstrain"), @@ -139,7 +149,7 @@ end function param_unconstrain_json!(sm::StanModel, theta::String, out::Vector{Float64}) dims = param_unc_num(sm) if length(out) != dims - error("out must be same size as number of unconstrained parameters") + throw(DimensionMismatch("out must be same size as number of unconstrained parameters")) end rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_unconstrain_json"), @@ -175,7 +185,7 @@ function log_density_gradient!(sm::StanModel, q, out::Vector{Float64}; propto = lp = Ref{Float64}(0.0) dims = param_unc_num(sm) if length(out) != dims - error("out must be same size as number of unconstrained parameters") + throw(DimensionMismatch("out must be same size as number of unconstrained parameters")) end rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density_gradient"), @@ -198,9 +208,9 @@ function log_density_hessian!(sm::StanModel, q, out_grad::Vector{Float64}, out_h lp = Ref{Float64}(0.0) dims = param_unc_num(sm) if length(out_grad) != dims - error("out_grad must be same size as number of unconstrained parameters") + throw(DimensionMismatch("out_grad must be same size as number of unconstrained parameters")) elseif length(out_hess) != dims*dims - error("out_hess must be same size as (number of unconstrained parameters)^2") + throw(DimensionMismatch("out_hess must be same size as (number of unconstrained parameters)^2")) end rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density_hessian"), diff --git a/julia/test/Manifest.toml b/julia/test/Manifest.toml index c5ac9232..53ac1e72 100644 --- a/julia/test/Manifest.toml +++ b/julia/test/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.8.0" manifest_format = "2.0" -project_hash = "70c6548fc0267b7c924ca6e56c4af9fd2abca604" +project_hash = "fc8a4bda90ee849287f19231f95018cefc4f54fb" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -38,6 +38,10 @@ deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" version = "0.3.20+0" +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + [[deps.Random]] deps = ["SHA", "Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -53,6 +57,9 @@ uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" diff --git a/julia/test/Project.toml b/julia/test/Project.toml index a2cd2f8d..00cebca5 100644 --- a/julia/test/Project.toml +++ b/julia/test/Project.toml @@ -1,3 +1,4 @@ [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 7d156208..0cc6b5c0 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -1,18 +1,328 @@ using BridgeStan using Test +using Printf + +function load_test_model(name::String, with_data=true) + lib = joinpath(@__DIR__, @sprintf("../../stan/%s/%s_model.so", name, name)) + if with_data + data = joinpath(@__DIR__, @sprintf("../../stan/%s/%s.data.json", name, name)) + else + data = "" + end + return BridgeStan.StanModel(lib, data) +end -@testset "bernoulli" begin - # Bernoulli - # CMDSTAN=/path/to/cmdstan/ make stan/bernoulli/bernoulli_model.so +@testset "constructor" begin + # no data + load_test_model("stdnormal", false) + # missing DSO + @test_throws SystemError BridgeStan.StanModel("Not going to find it") + # missing data + @test_throws SystemError load_test_model("stdnormal") + # exception in constructor + @test_throws ErrorException load_test_model("throw_data", false) +end + +@testset "name" begin + b = load_test_model("stdnormal", false) + @test BridgeStan.name(b) == "stdnormal_model" +end + +@testset "param_num" begin + b = load_test_model("full", false) + @test BridgeStan.param_num(b) == 1 + @test BridgeStan.param_num(b; include_tp=false) == 1 + @test BridgeStan.param_num(b; include_gq=false) == 1 + @test BridgeStan.param_num(b; include_tp=false, include_gq=false) == 1 + @test BridgeStan.param_num(b; include_gq=true) == 3 + @test BridgeStan.param_num(b; include_tp=false, include_gq=true) == 3 + @test BridgeStan.param_num(b; include_tp=true) == 2 + @test BridgeStan.param_num(b; include_tp=true, include_gq=false) == 2 + @test BridgeStan.param_num(b; include_tp=true, include_gq=true) == 4 +end + +@testset "param_unc_num" begin + b = load_test_model("simplex", false) + @test BridgeStan.param_num(b) == 5 + @test BridgeStan.param_unc_num(b) == 4 +end - function bernoulli(y, p) - sum(yn -> yn * log(p) + (1 - yn) * log(1 - p), y) +@testset "param_names" begin + b = load_test_model("matrix", false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_names(b) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_names(b; include_tp=false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_names(b; include_gq=false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_names(b; include_tp=false, include_gq=false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "B.1.1", "B.2.1", "B.3.1", "B.1.2", "B.2.2", "B.3.2"] == BridgeStan.param_names(b; include_tp=true) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "B.1.1", "B.2.1", "B.3.1", "B.1.2", "B.2.2", "B.3.2",] == BridgeStan.param_names(b; include_tp=true, include_gq=false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "c"] == BridgeStan.param_names(b; include_gq=true) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "c"] == BridgeStan.param_names(b; include_tp=false, include_gq=true) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "B.1.1", "B.2.1", "B.3.1", "B.1.2", "B.2.2", "B.3.2", "c",] == BridgeStan.param_names(b; include_tp=true, include_gq=true) +end + +@testset "param_unc_names" begin + b1 = load_test_model("matrix", false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_unc_names(b1) + + b2 = load_test_model("simplex", false) + @test ["theta.1", "theta.2", "theta.3", "theta.4"] == BridgeStan.param_unc_names(b2) + +end + +function _covariance_constrain_transform(v, D) + k = 0 + L = [j >= i ? (k += 1; v[k]) : 0 for i in 1:D, j in 1:D]' + for d in 1:D + L[d, d] = exp(L[d, d]) end + return L * L' +end - lib = joinpath(@__DIR__, "../../stan/bernoulli/bernoulli_model.so") - data = joinpath(@__DIR__, "../../stan/bernoulli/bernoulli.data.json") +@testset "param_constrain" begin + model = load_test_model("fr_gaussian") + D = 4 + size = 16 + unc_size = 10 + a = randn(unc_size) + B_expected = _covariance_constrain_transform(a, D) - model = BridgeStan.StanModel(lib, data) + b = BridgeStan.param_constrain(model, a) + B = reshape(b, (D, D)) + @test isapprox(B, B_expected) + + b = BridgeStan.param_constrain(model, a; include_tp=false) + B = reshape(b, (D, D)) + @test isapprox(B, B_expected) + + b = BridgeStan.param_constrain(model, a; include_gq=false) + B = reshape(b, (D, D)) + @test isapprox(B, B_expected) + + b = BridgeStan.param_constrain(model, a; include_tp=false, include_gq=false) + B = reshape(b, (D, D)) + @test isapprox(B, B_expected) + + # out test + scratch = zeros(16) + b = BridgeStan.param_constrain!(model, a, scratch) + @test b === scratch + B = reshape(b, (D, D)) + @test isapprox(B, B_expected) + scratch_wrong = zeros(10) + @test_throws DimensionMismatch BridgeStan.param_constrain!(model, a, scratch_wrong) + + + model2 = load_test_model("full", false) + @test 1 == length(BridgeStan.param_constrain(model2, a)) + @test 2 == length(BridgeStan.param_constrain(model2, a; include_tp=true)) + @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq=true)) + @test 4 == length(BridgeStan.param_constrain(model2, a; include_tp=true, include_gq=true)) + + # exception handling + model3 = load_test_model("throw_tp", false) + y = rand(1) + BridgeStan.param_constrain(model3, y) + @test_throws ErrorException BridgeStan.param_constrain(model3, y; include_tp=true) + + model4 = load_test_model("throw_gq", false) + BridgeStan.param_constrain(model4, y) + @test_throws ErrorException BridgeStan.param_constrain(model4, y; include_gq=true) +end + +@testset "param_unconstrain" begin + model = load_test_model("fr_gaussian") + unc_size = 10 + a = randn(unc_size) + b = BridgeStan.param_constrain(model, a) + c = BridgeStan.param_unconstrain(model, b) + @test isapprox(a, c) + + scratch = zeros(unc_size) + c2 = BridgeStan.param_unconstrain!(model, b, scratch) + @test c2 === scratch + @test isapprox(a, c2) + + scratch_wrong = zeros(16) + @test_throws DimensionMismatch BridgeStan.param_unconstrain!(model, b, scratch_wrong) +end + +@testset "param_unconstrain_json" begin + model = load_test_model("gaussian") + theta_unc = [0.2, log(1.9)] + theta_json = "{\"mu\": 0.2, \"sigma\": 1.9}" + theta_unc_test = BridgeStan.param_unconstrain_json(model, theta_json) + @test isapprox(theta_unc, theta_unc_test) + + scratch = zeros(2) + theta_unc_test2 = BridgeStan.param_unconstrain_json!(model, theta_json, scratch) + @test theta_unc_test2 === scratch + @test isapprox(theta_unc_test2, theta_unc) + + scratch_wrong = zeros(10) + @test_throws DimensionMismatch BridgeStan.param_unconstrain_json!(model, theta_json, scratch_wrong) + +end + +function _log_jacobian(p) + log.(p .* (1 .- p)) +end + +function _bernoulli(y, p) + sum(yn -> yn .* log.(p) + (1 - yn) .* log.(1 .- p), y) +end + +function _bernoulli_jacobian(y, p) + _bernoulli(y, p) .+ _log_jacobian(p) +end + +@testset "log_density" begin + model = load_test_model("bernoulli") + y = [0, 1, 0, 0, 0, 0, 0, 0, 0, 1] + x = rand(BridgeStan.param_unc_num(model)) + x_unc = log(x / (1 .- x)) + lp = BridgeStan.log_density(model, x_unc; propto=false, jacobian=false) + @test isapprox([lp], _bernoulli(y, x)) + lp2 = BridgeStan.log_density(model, x_unc; propto=false, jacobian=true) + @test isapprox([lp2], _bernoulli_jacobian(y, x)) + lp3 = BridgeStan.log_density(model, x_unc; propto=true, jacobian=true) + @test isapprox([lp2], _bernoulli_jacobian(y, x)) + lp4 = BridgeStan.log_density(model, x_unc; propto=true, jacobian=false) + @test isapprox([lp4], _bernoulli(y, x)) + + model2 = load_test_model("throw_lp", false) + y2 = rand(1) + @test_throws ErrorException BridgeStan.log_density(model2, y2) +end + + +@testset "jacobian model tests" begin + + function _logp(y_unc) + y = exp.(y_unc) + -0.5 .* y .^ 2 + end + + function _propto_false(y_unc) + -0.5 .* log(2 * pi) + end + + function _jacobian_true(y_unc) + y_unc + end + + function _grad_logp(y_unc) + y = exp.(y_unc) + return -1 .* (y .^ 2) + end + + function _grad_propto_false(y_unc) + 0 + end + + function _grad_jacobian_true(y_unc) + 1 + end + + function _hess_logp(y_unc) + y = exp.(y_unc) + -2.0 .* y .^ 2 + end + + function _hess_propto_false(y_unc) + 0 + end + + function _hess_jacobian_true(y_unc) + 0 + end + + @testset "log_density_gradient" begin + model = load_test_model("jacobian", false) + # test value, gradient, all combos +/- propto, +/- jacobian + y = abs.(randn(1)) + y_unc = log.(y) + ld, grad = BridgeStan.log_density_gradient(model, y_unc; propto=true, jacobian=true) + @test isapprox(_logp(y_unc) .+ _jacobian_true(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc) .+ _grad_jacobian_true(y_unc), grad) + + ld, grad = BridgeStan.log_density_gradient(model, y_unc; propto=true, jacobian=false) + @test isapprox(_logp(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc), grad) + + ld, grad = BridgeStan.log_density_gradient(model, y_unc; propto=false, jacobian=true) + @test isapprox(_logp(y_unc) .+ _propto_false(y_unc) .+ _jacobian_true(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc) .+ _grad_propto_false(y_unc) .+ _grad_jacobian_true(y_unc), grad) + + ld, grad = BridgeStan.log_density_gradient(model, y_unc; propto=false, jacobian=false) + @test isapprox(_logp(y_unc) .+ _propto_false(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc) .+ _grad_propto_false(y_unc), grad) + + # out parameters + scratch = zeros(BridgeStan.param_unc_num(model)) + ld, grad = BridgeStan.log_density_gradient!(model, y_unc, scratch; propto=true, jacobian=true) + @test grad === scratch + @test isapprox(_logp(y_unc) .+ _jacobian_true(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc) .+ _grad_jacobian_true(y_unc), grad) + + scratch_bad = zeros(BridgeStan.param_unc_num(model) + 10) + @test_throws DimensionMismatch BridgeStan.log_density_gradient!(model, y_unc, scratch_bad; propto=true, jacobian=true) + + end + + @testset "log_density_hessian" begin + + model = load_test_model("jacobian", false) + # test value, gradient, hessian, all combos +/- propto, +/- jacobian + y = abs.(randn(1)) + y_unc = log.(y) + ld, grad, hess = BridgeStan.log_density_hessian(model, y_unc; propto=true, jacobian=true) + @test isapprox(_logp(y_unc) .+ _jacobian_true(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc) .+ _grad_jacobian_true(y_unc), grad) + @test isapprox(_hess_logp(y_unc) .+ _hess_jacobian_true(y_unc), hess) + + ld, grad, hes = BridgeStan.log_density_hessian(model, y_unc; propto=true, jacobian=false) + @test isapprox(_logp(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc), grad) + @test isapprox(_hess_logp(y_unc), hess) + + ld, grad, hess = BridgeStan.log_density_hessian(model, y_unc; propto=false, jacobian=true) + @test isapprox(_logp(y_unc) .+ _propto_false(y_unc) .+ _jacobian_true(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc) .+ _grad_propto_false(y_unc) .+ _grad_jacobian_true(y_unc), grad) + @test isapprox(_hess_logp(y_unc) .+ _hess_propto_false(y_unc) .+ _hess_jacobian_true(y_unc), hess) + + ld, grad, hess = BridgeStan.log_density_hessian(model, y_unc; propto=false, jacobian=false) + @test isapprox(_logp(y_unc) .+ _propto_false(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc) .+ _grad_propto_false(y_unc), grad) + @test isapprox(_hess_logp(y_unc) .+ _hess_propto_false(y_unc), hess) + + # out parameters + dims = BridgeStan.param_unc_num(model) + grad_scratch = zeros(dims) + hess_scratch = zeros(dims ^ 2) + ld, grad, hess = BridgeStan.log_density_hessian!(model, y_unc, grad_scratch, hess_scratch; propto=true, jacobian=true) + @test grad === grad_scratch + @test isapprox(_logp(y_unc) .+ _jacobian_true(y_unc), [ld]) + @test isapprox(_grad_logp(y_unc) .+ _grad_jacobian_true(y_unc), grad) + @test isapprox(_hess_logp(y_unc) .+ _hess_jacobian_true(y_unc), hess) + @test isapprox(_hess_logp(y_unc) .+ _hess_jacobian_true(y_unc), reshape(hess_scratch, (dims, dims))) + + # force changes to prove underlying memory is same + hess_scratch[1] = 3 + @test reshape(hess, (dims ^ 2,)) == hess_scratch + + scratch_bad = zeros(BridgeStan.param_unc_num(model) + 10) + @test_throws DimensionMismatch BridgeStan.log_density_hessian!(model, y_unc, scratch_bad, hess_scratch; propto=true, jacobian=true) + @test_throws DimensionMismatch BridgeStan.log_density_hessian!(model, y_unc, grad_scratch, scratch_bad; propto=true, jacobian=true) + + end + +end + +@testset "bernoulli" begin + # Bernoulli + # CMDSTAN=/path/to/cmdstan/ make stan/bernoulli/bernoulli_model.so + + model = load_test_model("bernoulli") @test BridgeStan.name(model) == "bernoulli_model" @@ -22,15 +332,15 @@ using Test for _ in 1:R x = rand(BridgeStan.param_num(model)) q = @. log(x / (1 - x)) # unconstrained scale - (log_density, gradient) = BridgeStan.log_density_gradient(model, q, jacobian = 0) + (log_density, gradient) = BridgeStan.log_density_gradient(model, q, jacobian=0) p = x[1] - @test isapprox(log_density, bernoulli(y, p)) + @test isapprox(log_density, _bernoulli(y, p)) constrained_parameters = BridgeStan.param_constrain(model, q) @test isapprox(constrained_parameters, x) - unconstrained_parameters= BridgeStan.param_unconstrain(model, constrained_parameters) + unconstrained_parameters = BridgeStan.param_unconstrain(model, constrained_parameters) @test isapprox(unconstrained_parameters, q) end @@ -39,7 +349,7 @@ using Test end -@testset "multi" begin +@testset "threaded model: multi" begin # Multivariate Gaussian # CMDSTAN=/path/to/cmdstan/ make stan/multi/multi_model.so @@ -105,15 +415,6 @@ end # Full rank Gaussian # CMDSTAN=/path/to/cmdstan/ make stan/fr_gaussian/fr_gaussian_model.so - function _covariance_constrain_transform(v, D) - k = 0 - L = [j >= i ? (k += 1; v[k]) : 0 for i in 1:D, j in 1:D]' - for d in 1:D - L[d, d] = exp(L[d, d]) - end - return L * L' - end - lib = joinpath(@__DIR__, "../../stan/fr_gaussian/fr_gaussian_model.so") data = joinpath(@__DIR__, "../../stan/fr_gaussian/fr_gaussian.data.json") @@ -128,7 +429,7 @@ end D = 4 a = randn(unc_size) b = BridgeStan.param_constrain(model, a) - B = reshape(b, (D,D)) + B = reshape(b, (D, D)) B_expected = _covariance_constrain_transform(a, D) @test isapprox(B_expected, B) @@ -140,8 +441,8 @@ end pos = 1 for j = 1:4 for i = 1:4 - name_eq[pos] = names[pos] == ("Omega." * string(i) * "." * string(j)) - pos = pos + 1 + name_eq[pos] = names[pos] == ("Omega." * string(i) * "." * string(j)) + pos = pos + 1 end end @test all(name_eq) diff --git a/python/test/test_bridgestan.py b/python/test/test_bridgestan.py index 16713845..aa124160 100644 --- a/python/test/test_bridgestan.py +++ b/python/test/test_bridgestan.py @@ -294,6 +294,7 @@ def test_log_density(): lp2 = bridge.log_density(np.array([x_unc]), propto=False, jacobian=True) np.testing.assert_allclose(lp2, _bernoulli_jacobian(y, x)) lp3 = bridge.log_density(np.array([x_unc]), propto=True, jacobian=True) + print(_bernoulli_jacobian(y,x)) np.testing.assert_allclose(lp3, _bernoulli_jacobian(y, x)) lp4 = bridge.log_density(np.array([x_unc]), propto=True, jacobian=False) np.testing.assert_allclose(lp4, _bernoulli(y, x)) @@ -356,17 +357,6 @@ def _grad_jacobian_true(y_unc): grad[0], ) # - logdensity, grad = bridge.log_density_gradient( - y_unc_arr, propto=False, jacobian=True - ) - np.testing.assert_allclose( - _logp(y_unc) + _propto_false(y_unc) + _jacobian_true(y_unc), logdensity - ) - np.testing.assert_allclose( - _grad_logp(y_unc) + _grad_propto_false(y_unc) + _grad_jacobian_true(y_unc), - grad[0], - ) - # logdensity, grad = bridge.log_density_gradient( y_unc_arr, propto=False, jacobian=False ) From 86748f25d323725de07bab3f112430230229ba79 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 9 Sep 2022 15:09:39 -0400 Subject: [PATCH 05/11] Update example, apply JuliaFormatter.jl --- julia/MCMC.jl | 17 +-- julia/example.jl | 22 ++-- julia/src/BridgeStan.jl | 244 ++++++++++++++++++++++++++++------------ julia/test/runtests.jl | 229 +++++++++++++++++++++++++++---------- 4 files changed, 362 insertions(+), 150 deletions(-) diff --git a/julia/MCMC.jl b/julia/MCMC.jl index adcab878..84d40e79 100644 --- a/julia/MCMC.jl +++ b/julia/MCMC.jl @@ -1,5 +1,5 @@ struct HMCDiag - model::Bridgestan.StanModel + model::BridgeStan.StanModel stepsize::Float64 steps::Int64 metric::Vector{Float64} @@ -11,22 +11,23 @@ function HMCDiag(model, stepsize, steps) model, stepsize, steps, - ones(Bridgestan.param_unc_num(model)), - randn(Bridgestan.param_unc_num(model))) + ones(BridgeStan.param_unc_num(model)), + randn(BridgeStan.param_unc_num(model)), + ) end function joint_logp(hmc::HMCDiag, theta, rho) - logp, _ = Bridgestan.log_density_gradient(hmc.model, theta) + logp, _ = BridgeStan.log_density_gradient(hmc.model, theta) return logp - 0.5 * rho' * (hmc.metric .* rho) end function leapfrog(hmc::HMCDiag, theta, rho) e = hmc.stepsize .* hmc.metric - lp, grad = Bridgestan.log_density_gradient(hmc.model, theta) + lp, grad = BridgeStan.log_density_gradient(hmc.model, theta) rho_p = rho + 0.5 * hmc.stepsize .* grad - for n in 1:hmc.steps + for n = 1:hmc.steps theta .+= e .* rho_p - lp, grad = Bridgestan.log_density_gradient(hmc.model, theta) + lp, grad = BridgeStan.log_density_gradient(hmc.model, theta) if n != hmc.steps rho_p .+= e .* grad end @@ -36,7 +37,7 @@ function leapfrog(hmc::HMCDiag, theta, rho) end function sample(hmc::HMCDiag) - rho = randn(Bridgestan.param_unc_num(model)) + rho = randn(BridgeStan.param_unc_num(model)) logp = joint_logp(hmc, hmc.theta, rho) theta_prop, rho_prop = leapfrog(hmc, hmc.theta, rho) logp_prop = joint_logp(hmc, theta_prop, rho_prop) diff --git a/julia/example.jl b/julia/example.jl index f1faf0ef..7a5bd60c 100644 --- a/julia/example.jl +++ b/julia/example.jl @@ -1,4 +1,6 @@ -using Bridgestan +using BridgeStan + +const BS = BridgeStan # Bernoulli # CMDSTAN=/path/to/cmdstan/ make stan/bernoulli/bernoulli @@ -6,11 +8,11 @@ using Bridgestan bernoulli_lib = joinpath(@__DIR__, "../stan/bernoulli/bernoulli_model.so") bernoulli_data = joinpath(@__DIR__, "../stan/bernoulli/bernoulli.data.json") -smb = Bridgestan.StanModel(bernoulli_lib, bernoulli_data); -x = rand(Bridgestan.param_unc_num(smb)); +smb = BS.StanModel(bernoulli_lib, bernoulli_data); +x = rand(BS.param_unc_num(smb)); q = @. log(x / (1 - x)); # unconstrained scale -lp, grad = Bridgestan.log_density_gradient(smb, q, jacobian = 0) +lp, grad = BS.log_density_gradient(smb, q, jacobian = 0) println() println("log_density and gradient of Bernoulli model:") @@ -25,10 +27,10 @@ println() multi_lib = joinpath(@__DIR__, "../stan/multi/multi_model.so") multi_data = joinpath(@__DIR__, "../stan/multi/multi.data.json") -smm = Bridgestan.StanModel(multi_lib, multi_data) -x = randn(Bridgestan.param_unc_num(smm)); +smm = BS.StanModel(multi_lib, multi_data) +x = randn(BS.param_unc_num(smm)); -lp, grad = Bridgestan.log_density_gradient(smm, x) +lp, grad = BS.log_density_gradient(smm, x) println("log_density and gradient of Multivariate Gaussian model:") println((lp, grad)) @@ -39,15 +41,15 @@ println() include("./MCMC.jl") using Statistics -model = Bridgestan.StanModel(multi_lib, multi_data); +model = BS.StanModel(multi_lib, multi_data); stepsize = 0.25 steps = 10 hmcd = HMCDiag(model, stepsize, steps); M = 10_000 -theta = zeros(M, Bridgestan.param_unc_num(model)) -for m in 1:M +theta = zeros(M, BS.param_unc_num(model)) +for m = 1:M theta[m, :] .= sample(hmcd) end diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index fb09a87d..2930d8c0 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -1,7 +1,6 @@ module BridgeStan -export - StanModel, +export StanModel, name, param_num, param_unc_num, @@ -19,8 +18,7 @@ export log_density_gradient, log_density_hessian -mutable struct StanModelStruct -end +mutable struct StanModelStruct end mutable struct StanModel lib::Ptr{Nothing} @@ -37,16 +35,20 @@ mutable struct StanModel throw(SystemError("Dynamic library file not found")) end - if datafile_ != "" && !isfile(datafile_) + if datafile_ != "" && !isfile(datafile_) throw(SystemError("Data file not found")) end lib = Libc.Libdl.dlopen(stanlib_) - stanmodel = ccall(Libc.Libdl.dlsym(lib, "construct"), - Ptr{StanModelStruct}, - (Cstring, UInt32, UInt32), - datafile_, seed, chain_id) + stanmodel = ccall( + Libc.Libdl.dlsym(lib, "construct"), + Ptr{StanModelStruct}, + (Cstring, UInt32, UInt32), + datafile_, + seed, + chain_id, + ) if stanmodel == C_NULL error("could not construct model RNG") end @@ -54,10 +56,12 @@ mutable struct StanModel sm = new(lib, stanmodel, datafile_, seed, chain_id) function f(sm) - ccall(Libc.Libdl.dlsym(sm.lib, "destruct"), + ccall( + Libc.Libdl.dlsym(sm.lib, "destruct"), UInt32, (Ptr{StanModelStruct},), - sm.stanmodel) + sm.stanmodel, + ) end finalizer(f, sm) @@ -65,52 +69,80 @@ mutable struct StanModel end function name(sm::StanModel) - cstr = ccall(Libc.Libdl.dlsym(sm.lib, "name"), - Cstring, - (Ptr{StanModelStruct},), - sm.stanmodel) + cstr = ccall( + Libc.Libdl.dlsym(sm.lib, "name"), + Cstring, + (Ptr{StanModelStruct},), + sm.stanmodel, + ) unsafe_string(cstr) end function param_num(sm::StanModel; include_tp = false, include_gq = false) - ccall(Libc.Libdl.dlsym(sm.lib, "param_num"), - Cint, - (Ptr{StanModelStruct}, Cint, Cint), - sm.stanmodel, include_tp, include_gq) + ccall( + Libc.Libdl.dlsym(sm.lib, "param_num"), + Cint, + (Ptr{StanModelStruct}, Cint, Cint), + sm.stanmodel, + include_tp, + include_gq, + ) end function param_unc_num(sm::StanModel) - ccall(Libc.Libdl.dlsym(sm.lib, "param_unc_num"), - Cint, - (Ptr{StanModelStruct},), - sm.stanmodel) + ccall( + Libc.Libdl.dlsym(sm.lib, "param_unc_num"), + Cint, + (Ptr{StanModelStruct},), + sm.stanmodel, + ) end function param_names(sm::StanModel; include_tp = false, include_gq = false) - cstr = ccall(Libc.Libdl.dlsym(sm.lib, "param_names"), - Cstring, - (Ptr{StanModelStruct}, Cint, Cint), - sm.stanmodel, include_tp, include_gq) + cstr = ccall( + Libc.Libdl.dlsym(sm.lib, "param_names"), + Cstring, + (Ptr{StanModelStruct}, Cint, Cint), + sm.stanmodel, + include_tp, + include_gq, + ) [string(s) for s in split(unsafe_string(cstr), ',')] end function param_unc_names(sm::StanModel) - cstr = ccall(Libc.Libdl.dlsym(sm.lib, "param_unc_names"), - Cstring, - (Ptr{StanModelStruct},), - sm.stanmodel) + cstr = ccall( + Libc.Libdl.dlsym(sm.lib, "param_unc_names"), + Cstring, + (Ptr{StanModelStruct},), + sm.stanmodel, + ) [string(s) for s in split(unsafe_string(cstr), ',')] end -function param_constrain!(sm::StanModel, theta_unc, out::Vector{Float64}; include_tp=false, include_gq=false) - dims = param_num(sm; include_tp=include_tp, include_gq=include_gq) +function param_constrain!( + sm::StanModel, + theta_unc, + out::Vector{Float64}; + include_tp = false, + include_gq = false, +) + dims = param_num(sm; include_tp = include_tp, include_gq = include_gq) if length(out) != dims - throw(DimensionMismatch("out must be same size as number of constrained parameters")) + throw( + DimensionMismatch("out must be same size as number of constrained parameters"), + ) end - rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_constrain"), - Cint, - (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}), - sm.stanmodel, include_tp, include_gq, theta_unc, out) + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "param_constrain"), + Cint, + (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}), + sm.stanmodel, + include_tp, + include_gq, + theta_unc, + out, + ) if rc != 0 error("param_constrain failed on C++ side; see stderr for messages") else @@ -118,22 +150,30 @@ function param_constrain!(sm::StanModel, theta_unc, out::Vector{Float64}; includ end end -function param_constrain(sm::StanModel, theta_unc; include_tp=false, include_gq=false) - out = zeros(param_num(sm, include_tp=include_tp, include_gq=include_gq)) - param_constrain!(sm, theta_unc, out; include_tp=include_tp, include_gq=include_gq) +function param_constrain(sm::StanModel, theta_unc; include_tp = false, include_gq = false) + out = zeros(param_num(sm, include_tp = include_tp, include_gq = include_gq)) + param_constrain!(sm, theta_unc, out; include_tp = include_tp, include_gq = include_gq) end function param_unconstrain!(sm::StanModel, theta, out::Vector{Float64}) dims = param_unc_num(sm) if length(out) != dims - throw(DimensionMismatch("out must be same size as number of unconstrained parameters")) + throw( + DimensionMismatch( + "out must be same size as number of unconstrained parameters", + ), + ) end - rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_unconstrain"), - Cint, - (Ptr{StanModelStruct}, Ref{Cdouble}, Ref{Cdouble}), - sm.stanmodel, theta, out) + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "param_unconstrain"), + Cint, + (Ptr{StanModelStruct}, Ref{Cdouble}, Ref{Cdouble}), + sm.stanmodel, + theta, + out, + ) if rc != 0 error("param_unconstrain failed on C++ side; see stderr for messages") else @@ -149,13 +189,21 @@ end function param_unconstrain_json!(sm::StanModel, theta::String, out::Vector{Float64}) dims = param_unc_num(sm) if length(out) != dims - throw(DimensionMismatch("out must be same size as number of unconstrained parameters")) + throw( + DimensionMismatch( + "out must be same size as number of unconstrained parameters", + ), + ) end - rc = ccall(Libc.Libdl.dlsym(sm.lib, "param_unconstrain_json"), - Cint, - (Ptr{StanModelStruct}, Cstring, Ref{Cdouble}), - sm.stanmodel, theta, out) + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "param_unconstrain_json"), + Cint, + (Ptr{StanModelStruct}, Cstring, Ref{Cdouble}), + sm.stanmodel, + theta, + out, + ) if rc != 0 error("param_unconstrain_json failed on C++ side; see stderr for messages") else @@ -170,10 +218,16 @@ end function log_density(sm::StanModel, q; propto = true, jacobian = true) lp = Ref{Float64}(0.0) - rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density"), - Cint, - (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}), - sm.stanmodel, propto, jacobian, q, lp) + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "log_density"), + Cint, + (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}), + sm.stanmodel, + propto, + jacobian, + q, + lp, + ) if rc != 0 error("log_density failed on C++ side; see stderr for messages") else @@ -181,17 +235,34 @@ function log_density(sm::StanModel, q; propto = true, jacobian = true) end end -function log_density_gradient!(sm::StanModel, q, out::Vector{Float64}; propto = true, jacobian = true) +function log_density_gradient!( + sm::StanModel, + q, + out::Vector{Float64}; + propto = true, + jacobian = true, +) lp = Ref{Float64}(0.0) dims = param_unc_num(sm) if length(out) != dims - throw(DimensionMismatch("out must be same size as number of unconstrained parameters")) + throw( + DimensionMismatch( + "out must be same size as number of unconstrained parameters", + ), + ) end - rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density_gradient"), - Cint, - (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), - sm.stanmodel, propto, jacobian, q, lp, out) + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "log_density_gradient"), + Cint, + (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), + sm.stanmodel, + propto, + jacobian, + q, + lp, + out, + ) if rc != 0 error("log_density_gradient failed on C++ side; see stderr for messages") else @@ -201,22 +272,53 @@ end function log_density_gradient(sm::StanModel, q; propto = true, jacobian = true) grad = zeros(param_unc_num(sm)) - log_density_gradient!(sm, q, grad; propto=propto, jacobian=jacobian) + log_density_gradient!(sm, q, grad; propto = propto, jacobian = jacobian) end -function log_density_hessian!(sm::StanModel, q, out_grad::Vector{Float64}, out_hess::Vector{Float64}; propto = true, jacobian = true) +function log_density_hessian!( + sm::StanModel, + q, + out_grad::Vector{Float64}, + out_hess::Vector{Float64}; + propto = true, + jacobian = true, +) lp = Ref{Float64}(0.0) dims = param_unc_num(sm) if length(out_grad) != dims - throw(DimensionMismatch("out_grad must be same size as number of unconstrained parameters")) - elseif length(out_hess) != dims*dims - throw(DimensionMismatch("out_hess must be same size as (number of unconstrained parameters)^2")) + throw( + DimensionMismatch( + "out_grad must be same size as number of unconstrained parameters", + ), + ) + elseif length(out_hess) != dims * dims + throw( + DimensionMismatch( + "out_hess must be same size as (number of unconstrained parameters)^2", + ), + ) end - rc = ccall(Libc.Libdl.dlsym(sm.lib, "log_density_hessian"), - Cint, - (Ptr{StanModelStruct}, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), - sm.stanmodel, propto, jacobian, q, lp, out_grad, out_hess) + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "log_density_hessian"), + Cint, + ( + Ptr{StanModelStruct}, + Cint, + Cint, + Ref{Cdouble}, + Ref{Cdouble}, + Ref{Cdouble}, + Ref{Cdouble}, + ), + sm.stanmodel, + propto, + jacobian, + q, + lp, + out_grad, + out_hess, + ) if rc != 0 error("log_density_hessian failed on C++ side; see stderr for messages") else @@ -228,7 +330,7 @@ function log_density_hessian(sm::StanModel, q; propto = true, jacobian = true) dims = param_unc_num(sm) grad = zeros(dims) hess = zeros(dims * dims) - log_density_hessian!(sm,q, grad, hess; propto=propto, jacobian=jacobian) + log_density_hessian!(sm, q, grad, hess; propto = propto, jacobian = jacobian) end end diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 0cc6b5c0..65a9b33f 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -2,7 +2,7 @@ using BridgeStan using Test using Printf -function load_test_model(name::String, with_data=true) +function load_test_model(name::String, with_data = true) lib = joinpath(@__DIR__, @sprintf("../../stan/%s/%s_model.so", name, name)) if with_data data = joinpath(@__DIR__, @sprintf("../../stan/%s/%s.data.json", name, name)) @@ -31,14 +31,14 @@ end @testset "param_num" begin b = load_test_model("full", false) @test BridgeStan.param_num(b) == 1 - @test BridgeStan.param_num(b; include_tp=false) == 1 - @test BridgeStan.param_num(b; include_gq=false) == 1 - @test BridgeStan.param_num(b; include_tp=false, include_gq=false) == 1 - @test BridgeStan.param_num(b; include_gq=true) == 3 - @test BridgeStan.param_num(b; include_tp=false, include_gq=true) == 3 - @test BridgeStan.param_num(b; include_tp=true) == 2 - @test BridgeStan.param_num(b; include_tp=true, include_gq=false) == 2 - @test BridgeStan.param_num(b; include_tp=true, include_gq=true) == 4 + @test BridgeStan.param_num(b; include_tp = false) == 1 + @test BridgeStan.param_num(b; include_gq = false) == 1 + @test BridgeStan.param_num(b; include_tp = false, include_gq = false) == 1 + @test BridgeStan.param_num(b; include_gq = true) == 3 + @test BridgeStan.param_num(b; include_tp = false, include_gq = true) == 3 + @test BridgeStan.param_num(b; include_tp = true) == 2 + @test BridgeStan.param_num(b; include_tp = true, include_gq = false) == 2 + @test BridgeStan.param_num(b; include_tp = true, include_gq = true) == 4 end @testset "param_unc_num" begin @@ -49,20 +49,67 @@ end @testset "param_names" begin b = load_test_model("matrix", false) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_names(b) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_names(b; include_tp=false) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_names(b; include_gq=false) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_names(b; include_tp=false, include_gq=false) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "B.1.1", "B.2.1", "B.3.1", "B.1.2", "B.2.2", "B.3.2"] == BridgeStan.param_names(b; include_tp=true) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "B.1.1", "B.2.1", "B.3.1", "B.1.2", "B.2.2", "B.3.2",] == BridgeStan.param_names(b; include_tp=true, include_gq=false) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "c"] == BridgeStan.param_names(b; include_gq=true) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "c"] == BridgeStan.param_names(b; include_tp=false, include_gq=true) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "B.1.1", "B.2.1", "B.3.1", "B.1.2", "B.2.2", "B.3.2", "c",] == BridgeStan.param_names(b; include_tp=true, include_gq=true) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == + BridgeStan.param_names(b) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == + BridgeStan.param_names(b; include_tp = false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == + BridgeStan.param_names(b; include_gq = false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == + BridgeStan.param_names(b; include_tp = false, include_gq = false) + @test [ + "A.1.1", + "A.2.1", + "A.3.1", + "A.1.2", + "A.2.2", + "A.3.2", + "B.1.1", + "B.2.1", + "B.3.1", + "B.1.2", + "B.2.2", + "B.3.2", + ] == BridgeStan.param_names(b; include_tp = true) + @test [ + "A.1.1", + "A.2.1", + "A.3.1", + "A.1.2", + "A.2.2", + "A.3.2", + "B.1.1", + "B.2.1", + "B.3.1", + "B.1.2", + "B.2.2", + "B.3.2", + ] == BridgeStan.param_names(b; include_tp = true, include_gq = false) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "c"] == + BridgeStan.param_names(b; include_gq = true) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2", "c"] == + BridgeStan.param_names(b; include_tp = false, include_gq = true) + @test [ + "A.1.1", + "A.2.1", + "A.3.1", + "A.1.2", + "A.2.2", + "A.3.2", + "B.1.1", + "B.2.1", + "B.3.1", + "B.1.2", + "B.2.2", + "B.3.2", + "c", + ] == BridgeStan.param_names(b; include_tp = true, include_gq = true) end @testset "param_unc_names" begin b1 = load_test_model("matrix", false) - @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == BridgeStan.param_unc_names(b1) + @test ["A.1.1", "A.2.1", "A.3.1", "A.1.2", "A.2.2", "A.3.2"] == + BridgeStan.param_unc_names(b1) b2 = load_test_model("simplex", false) @test ["theta.1", "theta.2", "theta.3", "theta.4"] == BridgeStan.param_unc_names(b2) @@ -71,8 +118,8 @@ end function _covariance_constrain_transform(v, D) k = 0 - L = [j >= i ? (k += 1; v[k]) : 0 for i in 1:D, j in 1:D]' - for d in 1:D + L = [j >= i ? (k += 1; v[k]) : 0 for i = 1:D, j = 1:D]' + for d = 1:D L[d, d] = exp(L[d, d]) end return L * L' @@ -90,15 +137,15 @@ end B = reshape(b, (D, D)) @test isapprox(B, B_expected) - b = BridgeStan.param_constrain(model, a; include_tp=false) + b = BridgeStan.param_constrain(model, a; include_tp = false) B = reshape(b, (D, D)) @test isapprox(B, B_expected) - b = BridgeStan.param_constrain(model, a; include_gq=false) + b = BridgeStan.param_constrain(model, a; include_gq = false) B = reshape(b, (D, D)) @test isapprox(B, B_expected) - b = BridgeStan.param_constrain(model, a; include_tp=false, include_gq=false) + b = BridgeStan.param_constrain(model, a; include_tp = false, include_gq = false) B = reshape(b, (D, D)) @test isapprox(B, B_expected) @@ -114,19 +161,21 @@ end model2 = load_test_model("full", false) @test 1 == length(BridgeStan.param_constrain(model2, a)) - @test 2 == length(BridgeStan.param_constrain(model2, a; include_tp=true)) - @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq=true)) - @test 4 == length(BridgeStan.param_constrain(model2, a; include_tp=true, include_gq=true)) + @test 2 == length(BridgeStan.param_constrain(model2, a; include_tp = true)) + @test 3 == length(BridgeStan.param_constrain(model2, a; include_gq = true)) + @test 4 == length( + BridgeStan.param_constrain(model2, a; include_tp = true, include_gq = true), + ) # exception handling model3 = load_test_model("throw_tp", false) y = rand(1) BridgeStan.param_constrain(model3, y) - @test_throws ErrorException BridgeStan.param_constrain(model3, y; include_tp=true) + @test_throws ErrorException BridgeStan.param_constrain(model3, y; include_tp = true) model4 = load_test_model("throw_gq", false) BridgeStan.param_constrain(model4, y) - @test_throws ErrorException BridgeStan.param_constrain(model4, y; include_gq=true) + @test_throws ErrorException BridgeStan.param_constrain(model4, y; include_gq = true) end @testset "param_unconstrain" begin @@ -159,7 +208,11 @@ end @test isapprox(theta_unc_test2, theta_unc) scratch_wrong = zeros(10) - @test_throws DimensionMismatch BridgeStan.param_unconstrain_json!(model, theta_json, scratch_wrong) + @test_throws DimensionMismatch BridgeStan.param_unconstrain_json!( + model, + theta_json, + scratch_wrong, + ) end @@ -180,13 +233,13 @@ end y = [0, 1, 0, 0, 0, 0, 0, 0, 0, 1] x = rand(BridgeStan.param_unc_num(model)) x_unc = log(x / (1 .- x)) - lp = BridgeStan.log_density(model, x_unc; propto=false, jacobian=false) + lp = BridgeStan.log_density(model, x_unc; propto = false, jacobian = false) @test isapprox([lp], _bernoulli(y, x)) - lp2 = BridgeStan.log_density(model, x_unc; propto=false, jacobian=true) + lp2 = BridgeStan.log_density(model, x_unc; propto = false, jacobian = true) @test isapprox([lp2], _bernoulli_jacobian(y, x)) - lp3 = BridgeStan.log_density(model, x_unc; propto=true, jacobian=true) + lp3 = BridgeStan.log_density(model, x_unc; propto = true, jacobian = true) @test isapprox([lp2], _bernoulli_jacobian(y, x)) - lp4 = BridgeStan.log_density(model, x_unc; propto=true, jacobian=false) + lp4 = BridgeStan.log_density(model, x_unc; propto = true, jacobian = false) @test isapprox([lp4], _bernoulli(y, x)) model2 = load_test_model("throw_lp", false) @@ -241,31 +294,50 @@ end # test value, gradient, all combos +/- propto, +/- jacobian y = abs.(randn(1)) y_unc = log.(y) - ld, grad = BridgeStan.log_density_gradient(model, y_unc; propto=true, jacobian=true) + ld, grad = + BridgeStan.log_density_gradient(model, y_unc; propto = true, jacobian = true) @test isapprox(_logp(y_unc) .+ _jacobian_true(y_unc), [ld]) @test isapprox(_grad_logp(y_unc) .+ _grad_jacobian_true(y_unc), grad) - ld, grad = BridgeStan.log_density_gradient(model, y_unc; propto=true, jacobian=false) + ld, grad = + BridgeStan.log_density_gradient(model, y_unc; propto = true, jacobian = false) @test isapprox(_logp(y_unc), [ld]) @test isapprox(_grad_logp(y_unc), grad) - ld, grad = BridgeStan.log_density_gradient(model, y_unc; propto=false, jacobian=true) + ld, grad = + BridgeStan.log_density_gradient(model, y_unc; propto = false, jacobian = true) @test isapprox(_logp(y_unc) .+ _propto_false(y_unc) .+ _jacobian_true(y_unc), [ld]) - @test isapprox(_grad_logp(y_unc) .+ _grad_propto_false(y_unc) .+ _grad_jacobian_true(y_unc), grad) + @test isapprox( + _grad_logp(y_unc) .+ _grad_propto_false(y_unc) .+ _grad_jacobian_true(y_unc), + grad, + ) - ld, grad = BridgeStan.log_density_gradient(model, y_unc; propto=false, jacobian=false) + ld, grad = + BridgeStan.log_density_gradient(model, y_unc; propto = false, jacobian = false) @test isapprox(_logp(y_unc) .+ _propto_false(y_unc), [ld]) @test isapprox(_grad_logp(y_unc) .+ _grad_propto_false(y_unc), grad) # out parameters scratch = zeros(BridgeStan.param_unc_num(model)) - ld, grad = BridgeStan.log_density_gradient!(model, y_unc, scratch; propto=true, jacobian=true) + ld, grad = BridgeStan.log_density_gradient!( + model, + y_unc, + scratch; + propto = true, + jacobian = true, + ) @test grad === scratch @test isapprox(_logp(y_unc) .+ _jacobian_true(y_unc), [ld]) @test isapprox(_grad_logp(y_unc) .+ _grad_jacobian_true(y_unc), grad) scratch_bad = zeros(BridgeStan.param_unc_num(model) + 10) - @test_throws DimensionMismatch BridgeStan.log_density_gradient!(model, y_unc, scratch_bad; propto=true, jacobian=true) + @test_throws DimensionMismatch BridgeStan.log_density_gradient!( + model, + y_unc, + scratch_bad; + propto = true, + jacobian = true, + ) end @@ -275,22 +347,32 @@ end # test value, gradient, hessian, all combos +/- propto, +/- jacobian y = abs.(randn(1)) y_unc = log.(y) - ld, grad, hess = BridgeStan.log_density_hessian(model, y_unc; propto=true, jacobian=true) + ld, grad, hess = + BridgeStan.log_density_hessian(model, y_unc; propto = true, jacobian = true) @test isapprox(_logp(y_unc) .+ _jacobian_true(y_unc), [ld]) @test isapprox(_grad_logp(y_unc) .+ _grad_jacobian_true(y_unc), grad) @test isapprox(_hess_logp(y_unc) .+ _hess_jacobian_true(y_unc), hess) - ld, grad, hes = BridgeStan.log_density_hessian(model, y_unc; propto=true, jacobian=false) + ld, grad, hes = + BridgeStan.log_density_hessian(model, y_unc; propto = true, jacobian = false) @test isapprox(_logp(y_unc), [ld]) @test isapprox(_grad_logp(y_unc), grad) @test isapprox(_hess_logp(y_unc), hess) - ld, grad, hess = BridgeStan.log_density_hessian(model, y_unc; propto=false, jacobian=true) + ld, grad, hess = + BridgeStan.log_density_hessian(model, y_unc; propto = false, jacobian = true) @test isapprox(_logp(y_unc) .+ _propto_false(y_unc) .+ _jacobian_true(y_unc), [ld]) - @test isapprox(_grad_logp(y_unc) .+ _grad_propto_false(y_unc) .+ _grad_jacobian_true(y_unc), grad) - @test isapprox(_hess_logp(y_unc) .+ _hess_propto_false(y_unc) .+ _hess_jacobian_true(y_unc), hess) - - ld, grad, hess = BridgeStan.log_density_hessian(model, y_unc; propto=false, jacobian=false) + @test isapprox( + _grad_logp(y_unc) .+ _grad_propto_false(y_unc) .+ _grad_jacobian_true(y_unc), + grad, + ) + @test isapprox( + _hess_logp(y_unc) .+ _hess_propto_false(y_unc) .+ _hess_jacobian_true(y_unc), + hess, + ) + + ld, grad, hess = + BridgeStan.log_density_hessian(model, y_unc; propto = false, jacobian = false) @test isapprox(_logp(y_unc) .+ _propto_false(y_unc), [ld]) @test isapprox(_grad_logp(y_unc) .+ _grad_propto_false(y_unc), grad) @test isapprox(_hess_logp(y_unc) .+ _hess_propto_false(y_unc), hess) @@ -298,21 +380,45 @@ end # out parameters dims = BridgeStan.param_unc_num(model) grad_scratch = zeros(dims) - hess_scratch = zeros(dims ^ 2) - ld, grad, hess = BridgeStan.log_density_hessian!(model, y_unc, grad_scratch, hess_scratch; propto=true, jacobian=true) + hess_scratch = zeros(dims^2) + ld, grad, hess = BridgeStan.log_density_hessian!( + model, + y_unc, + grad_scratch, + hess_scratch; + propto = true, + jacobian = true, + ) @test grad === grad_scratch @test isapprox(_logp(y_unc) .+ _jacobian_true(y_unc), [ld]) @test isapprox(_grad_logp(y_unc) .+ _grad_jacobian_true(y_unc), grad) @test isapprox(_hess_logp(y_unc) .+ _hess_jacobian_true(y_unc), hess) - @test isapprox(_hess_logp(y_unc) .+ _hess_jacobian_true(y_unc), reshape(hess_scratch, (dims, dims))) + @test isapprox( + _hess_logp(y_unc) .+ _hess_jacobian_true(y_unc), + reshape(hess_scratch, (dims, dims)), + ) # force changes to prove underlying memory is same hess_scratch[1] = 3 - @test reshape(hess, (dims ^ 2,)) == hess_scratch + @test reshape(hess, (dims^2,)) == hess_scratch scratch_bad = zeros(BridgeStan.param_unc_num(model) + 10) - @test_throws DimensionMismatch BridgeStan.log_density_hessian!(model, y_unc, scratch_bad, hess_scratch; propto=true, jacobian=true) - @test_throws DimensionMismatch BridgeStan.log_density_hessian!(model, y_unc, grad_scratch, scratch_bad; propto=true, jacobian=true) + @test_throws DimensionMismatch BridgeStan.log_density_hessian!( + model, + y_unc, + scratch_bad, + hess_scratch; + propto = true, + jacobian = true, + ) + @test_throws DimensionMismatch BridgeStan.log_density_hessian!( + model, + y_unc, + grad_scratch, + scratch_bad; + propto = true, + jacobian = true, + ) end @@ -329,10 +435,10 @@ end y = [0, 1, 0, 0, 0, 0, 0, 0, 0, 1] R = 1000 - for _ in 1:R + for _ = 1:R x = rand(BridgeStan.param_num(model)) q = @. log(x / (1 - x)) # unconstrained scale - (log_density, gradient) = BridgeStan.log_density_gradient(model, q, jacobian=0) + (log_density, gradient) = BridgeStan.log_density_gradient(model, q, jacobian = 0) p = x[1] @test isapprox(log_density, _bernoulli(y, p)) @@ -340,7 +446,8 @@ end constrained_parameters = BridgeStan.param_constrain(model, q) @test isapprox(constrained_parameters, x) - unconstrained_parameters = BridgeStan.param_unconstrain(model, constrained_parameters) + unconstrained_parameters = + BridgeStan.param_unconstrain(model, constrained_parameters) @test isapprox(unconstrained_parameters, q) end @@ -371,8 +478,8 @@ end ld = Vector{Bool}(undef, R) g = Vector{Bool}(undef, R) - @sync for it in 1:nt - Threads.@spawn for r in it:nt:R + @sync for it = 1:nt + Threads.@spawn for r = it:nt:R x = randn(BridgeStan.param_num(model)) (lp, grad) = BridgeStan.log_density_gradient(model, x) @@ -423,7 +530,7 @@ end size = 16 unc_size = 10 - @test isapprox(size, BridgeStan.param_num(model, include_tp=true, include_gq=true)) + @test isapprox(size, BridgeStan.param_num(model, include_tp = true, include_gq = true)) @test isapprox(unc_size, BridgeStan.param_unc_num(model)) D = 4 @@ -436,7 +543,7 @@ end c = BridgeStan.param_unconstrain(model, b) @test isapprox(a, c) - names = BridgeStan.param_names(model, include_tp=true, include_gq=true) + names = BridgeStan.param_names(model, include_tp = true, include_gq = true) name_eq = Vector{Bool}(undef, size) pos = 1 for j = 1:4 From 15fbc6741cf5f3e287452e392e73ff393ff68f75 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 9 Sep 2022 15:11:16 -0400 Subject: [PATCH 06/11] Remove debug print --- python/test/test_bridgestan.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/test/test_bridgestan.py b/python/test/test_bridgestan.py index aa124160..da893789 100644 --- a/python/test/test_bridgestan.py +++ b/python/test/test_bridgestan.py @@ -294,7 +294,6 @@ def test_log_density(): lp2 = bridge.log_density(np.array([x_unc]), propto=False, jacobian=True) np.testing.assert_allclose(lp2, _bernoulli_jacobian(y, x)) lp3 = bridge.log_density(np.array([x_unc]), propto=True, jacobian=True) - print(_bernoulli_jacobian(y,x)) np.testing.assert_allclose(lp3, _bernoulli_jacobian(y, x)) lp4 = bridge.log_density(np.array([x_unc]), propto=True, jacobian=False) np.testing.assert_allclose(lp4, _bernoulli(y, x)) From beca54e839ce901b8da0b28cae0fb74e0357b8bb Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 9 Sep 2022 15:13:29 -0400 Subject: [PATCH 07/11] Clean up older tests --- julia/test/runtests.jl | 20 ++++---------------- 1 file changed, 4 insertions(+), 16 deletions(-) diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index 65a9b33f..ac347200 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -468,11 +468,8 @@ end return -x end - lib = joinpath(@__DIR__, "../../stan/multi/multi_model.so") - data = joinpath(@__DIR__, "../../stan/multi/multi.data.json") - + model = load_test_model("multi") nt = Threads.nthreads() - model = BridgeStan.StanModel(lib, data) R = 1000 ld = Vector{Bool}(undef, R) @@ -497,10 +494,7 @@ end # Guassian with positive constrained standard deviation # CMDSTAN=/path/to/cmdstan/ make stan/gaussian/gaussian_model.so - lib = joinpath(@__DIR__, "../../stan/gaussian/gaussian_model.so") - data = joinpath(@__DIR__, "../../stan/gaussian/gaussian.data.json") - - model = BridgeStan.StanModel(lib, data) + model = load_test_model("gaussian") theta = [0.2, 1.9] theta_unc = [0.2, log(1.9)] @@ -522,10 +516,7 @@ end # Full rank Gaussian # CMDSTAN=/path/to/cmdstan/ make stan/fr_gaussian/fr_gaussian_model.so - lib = joinpath(@__DIR__, "../../stan/fr_gaussian/fr_gaussian_model.so") - data = joinpath(@__DIR__, "../../stan/fr_gaussian/fr_gaussian.data.json") - - model = BridgeStan.StanModel(lib, data) + model = load_test_model("fr_gaussian") size = 16 unc_size = 10 @@ -564,10 +555,7 @@ end @testset "simple" begin - lib = joinpath(@__DIR__, "../../stan/simple/simple_model.so") - data = joinpath(@__DIR__, "../../stan/simple/simple.data.json") - - model = BridgeStan.StanModel(lib, data) + model = load_test_model("simple") D = 5 y = rand(D) From 9f4d7fa18cddc5dd172152136b386d653898dc7b Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 9 Sep 2022 15:41:28 -0400 Subject: [PATCH 08/11] Typo fix --- julia/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index ac347200..c43ae900 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -238,7 +238,7 @@ end lp2 = BridgeStan.log_density(model, x_unc; propto = false, jacobian = true) @test isapprox([lp2], _bernoulli_jacobian(y, x)) lp3 = BridgeStan.log_density(model, x_unc; propto = true, jacobian = true) - @test isapprox([lp2], _bernoulli_jacobian(y, x)) + @test isapprox([lp3], _bernoulli_jacobian(y, x)) lp4 = BridgeStan.log_density(model, x_unc; propto = true, jacobian = false) @test isapprox([lp4], _bernoulli(y, x)) From 57351e7d9a2c0009066b9ef19f807e7e055bf733 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 13 Sep 2022 09:34:15 -0400 Subject: [PATCH 09/11] Changes per review --- julia/src/BridgeStan.jl | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index 2930d8c0..b074eb46 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -107,7 +107,7 @@ function param_names(sm::StanModel; include_tp = false, include_gq = false) include_tp, include_gq, ) - [string(s) for s in split(unsafe_string(cstr), ',')] + string.(split(unsafe_string(cstr), ',')) end function param_unc_names(sm::StanModel) @@ -117,7 +117,7 @@ function param_unc_names(sm::StanModel) (Ptr{StanModelStruct},), sm.stanmodel, ) - [string(s) for s in split(unsafe_string(cstr), ',')] + string.(split(unsafe_string(cstr), ',')) end function param_constrain!( @@ -145,9 +145,8 @@ function param_constrain!( ) if rc != 0 error("param_constrain failed on C++ side; see stderr for messages") - else - out end + out end function param_constrain(sm::StanModel, theta_unc; include_tp = false, include_gq = false) @@ -176,9 +175,8 @@ function param_unconstrain!(sm::StanModel, theta, out::Vector{Float64}) ) if rc != 0 error("param_unconstrain failed on C++ side; see stderr for messages") - else - out end + out end function param_unconstrain(sm::StanModel, theta) @@ -206,9 +204,8 @@ function param_unconstrain_json!(sm::StanModel, theta::String, out::Vector{Float ) if rc != 0 error("param_unconstrain_json failed on C++ side; see stderr for messages") - else - out end + out end function param_unconstrain_json(sm::StanModel, theta::String) @@ -217,7 +214,7 @@ function param_unconstrain_json(sm::StanModel, theta::String) end function log_density(sm::StanModel, q; propto = true, jacobian = true) - lp = Ref{Float64}(0.0) + lp = Ref(0.0) rc = ccall( Libc.Libdl.dlsym(sm.lib, "log_density"), Cint, @@ -230,9 +227,8 @@ function log_density(sm::StanModel, q; propto = true, jacobian = true) ) if rc != 0 error("log_density failed on C++ side; see stderr for messages") - else - lp[] end + lp[] end function log_density_gradient!( @@ -242,7 +238,7 @@ function log_density_gradient!( propto = true, jacobian = true, ) - lp = Ref{Float64}(0.0) + lp = Ref(0.0) dims = param_unc_num(sm) if length(out) != dims throw( @@ -265,9 +261,8 @@ function log_density_gradient!( ) if rc != 0 error("log_density_gradient failed on C++ side; see stderr for messages") - else - (lp[], out) end + (lp[], out) end function log_density_gradient(sm::StanModel, q; propto = true, jacobian = true) @@ -283,7 +278,7 @@ function log_density_hessian!( propto = true, jacobian = true, ) - lp = Ref{Float64}(0.0) + lp = Ref(0.0) dims = param_unc_num(sm) if length(out_grad) != dims throw( @@ -321,9 +316,9 @@ function log_density_hessian!( ) if rc != 0 error("log_density_hessian failed on C++ side; see stderr for messages") - else - (lp[], out_grad, reshape(out_hess, (dims, dims))) + end + (lp[], out_grad, reshape(out_hess, (dims, dims))) end function log_density_hessian(sm::StanModel, q; propto = true, jacobian = true) From b3a01ebb2a06b9347323c11de390d252dec65b9c Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 13 Sep 2022 12:37:26 -0400 Subject: [PATCH 10/11] Mark all arguments passed to C as Vector{Float64} --- julia/src/BridgeStan.jl | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index b074eb46..58f53d21 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -122,7 +122,7 @@ end function param_constrain!( sm::StanModel, - theta_unc, + theta_unc::Vector{Float64}, out::Vector{Float64}; include_tp = false, include_gq = false, @@ -149,13 +149,18 @@ function param_constrain!( out end -function param_constrain(sm::StanModel, theta_unc; include_tp = false, include_gq = false) +function param_constrain( + sm::StanModel, + theta_unc::Vector{Float64}; + include_tp = false, + include_gq = false, +) out = zeros(param_num(sm, include_tp = include_tp, include_gq = include_gq)) param_constrain!(sm, theta_unc, out; include_tp = include_tp, include_gq = include_gq) end -function param_unconstrain!(sm::StanModel, theta, out::Vector{Float64}) +function param_unconstrain!(sm::StanModel, theta::Vector{Float64}, out::Vector{Float64}) dims = param_unc_num(sm) if length(out) != dims throw( @@ -179,7 +184,7 @@ function param_unconstrain!(sm::StanModel, theta, out::Vector{Float64}) out end -function param_unconstrain(sm::StanModel, theta) +function param_unconstrain(sm::StanModel, theta::Vector{Float64}) out = zeros(param_unc_num(sm)) param_unconstrain!(sm, theta, out) end @@ -213,7 +218,7 @@ function param_unconstrain_json(sm::StanModel, theta::String) param_unconstrain_json!(sm, theta, out) end -function log_density(sm::StanModel, q; propto = true, jacobian = true) +function log_density(sm::StanModel, q::Vector{Float64}; propto = true, jacobian = true) lp = Ref(0.0) rc = ccall( Libc.Libdl.dlsym(sm.lib, "log_density"), @@ -233,7 +238,7 @@ end function log_density_gradient!( sm::StanModel, - q, + q::Vector{Float64}, out::Vector{Float64}; propto = true, jacobian = true, @@ -265,14 +270,19 @@ function log_density_gradient!( (lp[], out) end -function log_density_gradient(sm::StanModel, q; propto = true, jacobian = true) +function log_density_gradient( + sm::StanModel, + q::Vector{Float64}; + propto = true, + jacobian = true, +) grad = zeros(param_unc_num(sm)) log_density_gradient!(sm, q, grad; propto = propto, jacobian = jacobian) end function log_density_hessian!( sm::StanModel, - q, + q::Vector{Float64}, out_grad::Vector{Float64}, out_hess::Vector{Float64}; propto = true, @@ -321,7 +331,12 @@ function log_density_hessian!( (lp[], out_grad, reshape(out_hess, (dims, dims))) end -function log_density_hessian(sm::StanModel, q; propto = true, jacobian = true) +function log_density_hessian( + sm::StanModel, + q::Vector{Float64}; + propto = true, + jacobian = true, +) dims = param_unc_num(sm) grad = zeros(dims) hess = zeros(dims * dims) From 5cac08d69ade12f2e487b811110587830c3f7cf0 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 13 Sep 2022 12:42:32 -0400 Subject: [PATCH 11/11] Fix test --- julia/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/julia/test/runtests.jl b/julia/test/runtests.jl index c43ae900..27cd2eba 100644 --- a/julia/test/runtests.jl +++ b/julia/test/runtests.jl @@ -232,7 +232,7 @@ end model = load_test_model("bernoulli") y = [0, 1, 0, 0, 0, 0, 0, 0, 0, 1] x = rand(BridgeStan.param_unc_num(model)) - x_unc = log(x / (1 .- x)) + x_unc = @. log(x / (1 - x)) lp = BridgeStan.log_density(model, x_unc; propto = false, jacobian = false) @test isapprox([lp], _bernoulli(y, x)) lp2 = BridgeStan.log_density(model, x_unc; propto = false, jacobian = true)