From cdd06afc2d8f618a559c09d50ae2c6ed79a380ce Mon Sep 17 00:00:00 2001 From: Mauricio Vanzulli Date: Sun, 14 Dec 2025 11:34:12 -0300 Subject: [PATCH] Update to new julia version --- Project.toml | 12 +- .../julia/update_rules/delta_sampling.jl | 176 +++++++++--------- 2 files changed, 96 insertions(+), 92 deletions(-) diff --git a/Project.toml b/Project.toml index 6a075c5d..6f328d45 100644 --- a/Project.toml +++ b/Project.toml @@ -21,12 +21,12 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -DataStructures = "0.17, 0.18" -Documenter = "0.25, 0.26, 0.27" -ForwardDiff = "0.10" +DataStructures = "0.17, 0.18, 0.19, 0.20" +Documenter = "0.25, 0.26, 0.27, 1" +ForwardDiff = "0.10, 1" PositiveFactorizations = "0.2" Roots = "1.3, 2" -SpecialFunctions = "0.8.1, 0.9, 1.0, 2.0" -StatsBase = "0.32, 0.33" +SpecialFunctions = "0.8.1, 0.9, 1, 2" +StatsBase = "0.32, 0.33, 0.34" StatsFuns = "0.9, 1" -julia = "1.3, 1.4, 1.5, 1.6, 1.7" +julia = "1.6, 1.7, 1.8, 1.9, 1.10, 1.11" diff --git a/src/engines/julia/update_rules/delta_sampling.jl b/src/engines/julia/update_rules/delta_sampling.jl index 9e793fa8..c9eae051 100644 --- a/src/engines/julia/update_rules/delta_sampling.jl +++ b/src/engines/julia/update_rules/delta_sampling.jl @@ -1,11 +1,11 @@ export -ruleSPDeltaSOutNM, -ruleSPDeltaSIn1MN, -ruleSPDeltaSInGX, -ruleSPDeltaSOutNMX, -ruleSPDeltaSInMX, -ruleMDeltaSInMGX, -prod! + ruleSPDeltaSOutNM, + ruleSPDeltaSIn1MN, + ruleSPDeltaSInGX, + ruleSPDeltaSOutNMX, + ruleSPDeltaSInMX, + ruleMDeltaSInMGX, + prod! const default_n_samples = 1000 # Default value for the number of samples @@ -15,22 +15,22 @@ const default_n_samples = 1000 # Default value for the number of samples #---------------------- function ruleSPDeltaSOutNM(g::Function, - msg_out::Nothing, - msg_in1::Message; # Applies to any message except SampleList - dims::Any=nothing, - n_samples=default_n_samples) + msg_out::Nothing, + msg_in1::Message; # Applies to any message except SampleList + dims::Any=nothing, + n_samples=default_n_samples) samples = g.(sample(msg_in1.dist, n_samples)) - weights = ones(n_samples)/n_samples + weights = ones(n_samples) / n_samples return Message(variateType(dims), SampleList, s=samples, w=weights) end function ruleSPDeltaSOutNM(g::Function, - msg_out::Nothing, - msg_in1::Message{SampleList}; # Special case for SampleList - dims::Any=nothing, - n_samples=default_n_samples) + msg_out::Nothing, + msg_in1::Message{SampleList}; # Special case for SampleList + dims::Any=nothing, + n_samples=default_n_samples) samples = g.(msg_in1.dist.params[:s]) weights = msg_in1.dist.params[:w] @@ -39,20 +39,20 @@ function ruleSPDeltaSOutNM(g::Function, end function ruleSPDeltaSIn1MN(g::Function, - msg_out::Message, - msg_in1::Nothing; - dims::Any=nothing, - n_samples=default_n_samples) + msg_out::Message, + msg_in1::Nothing; + dims::Any=nothing, + n_samples=default_n_samples) - return Message(variateType(dims), Function, log_pdf = (z)->logPdf(msg_out.dist, g(z))) + return Message(variateType(dims), Function, log_pdf=(z) -> logPdf(msg_out.dist, g(z))) end function ruleSPDeltaSInGX(g::Function, - inx::Int64, # Index of inbound interface inx - msg_out::Message, - msgs_in::Vararg{Message{<:Gaussian}}; - dims::Any=nothing, - n_samples=default_n_samples) + inx::Int64, # Index of inbound interface inx + msg_out::Message, + msgs_in::Vararg{Message{<:Gaussian}}; + dims::Any=nothing, + n_samples=default_n_samples) # Extract joint statistics of inbound messages (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Return arrays with individual means and covariances @@ -69,7 +69,7 @@ function ruleSPDeltaSInGX(g::Function, # Marginalize joint belief on in's (m_inx, V_inx) = marginalizeGaussianMV(m_in, V_in, ds, inx) W_inx = cholinv(V_inx) # Convert to canonical statistics - xi_inx = W_inx*m_inx + xi_inx = W_inx * m_inx # Divide marginal on inx by forward message (xi_fw_inx, W_fw_inx) = unsafeWeightedMeanPrecision(msgs_in[inx].dist) @@ -80,29 +80,29 @@ function ruleSPDeltaSInGX(g::Function, end function ruleSPDeltaSOutNMX(g::Function, - msg_out::Nothing, - msgs_in::Vararg{Message}; - dims::Any=nothing, - n_samples=default_n_samples) - + msg_out::Nothing, + msgs_in::Vararg{Message}; + dims::Any=nothing, + n_samples=default_n_samples) + samples_in = [sample(msg_in.dist, n_samples) for msg_in in msgs_in] samples = g.(samples_in...) - weights = ones(n_samples)/n_samples + weights = ones(n_samples) / n_samples return Message(variateType(dims), SampleList, s=samples, w=weights) end function ruleSPDeltaSInMX(g::Function, - inx::Int64, # Index of inbound interface inx - msg_out::Message, - msgs_in::Vararg{Message}; - dims::Any=nothing, - n_samples=default_n_samples) + inx::Int64, # Index of inbound interface inx + msg_out::Message, + msgs_in::Vararg{Message}; + dims::Any=nothing, + n_samples=default_n_samples) arg_sample = (z) -> begin samples_in = [] - for i=1:length(msgs_in) - if i==inx + for i = 1:length(msgs_in) + if i == inx push!(samples_in, fill(z, n_samples)) else push!(samples_in, sample(msgs_in[i].dist, n_samples)) @@ -112,40 +112,40 @@ function ruleSPDeltaSInMX(g::Function, return samples_in end - approximate_pdf(z) = sum(exp.(logPdf.([msg_out.dist],g.(arg_sample(z)...))))/n_samples + approximate_pdf(z) = sum(exp.(logPdf.([msg_out.dist], g.(arg_sample(z)...)))) / n_samples - return Message(variateType(dims), Function, log_pdf = (z)->log(approximate_pdf(z))) + return Message(variateType(dims), Function, log_pdf=(z) -> log(approximate_pdf(z))) end # Special case for two inputs with one PointMass (no inx required) function ruleSPDeltaSInMX(g::Function, - msg_out::Message, - msg_in1::Message{PointMass}, - msg_in2::Nothing; - dims::Any=nothing, - n_samples=default_n_samples) + msg_out::Message, + msg_in1::Message{PointMass}, + msg_in2::Nothing; + dims::Any=nothing, + n_samples=default_n_samples) m_in1 = msg_in1.dist.params[:m] - return Message(variateType(dims), Function, log_pdf = (z)->logPdf(msg_out.dist, g(m_in1, z))) + return Message(variateType(dims), Function, log_pdf=(z) -> logPdf(msg_out.dist, g(m_in1, z))) end # Special case for two inputs with one PointMass (no inx required) function ruleSPDeltaSInMX(g::Function, - msg_out::Message, - msg_in1::Nothing, - msg_in2::Message{PointMass}; - dims::Any=nothing, - n_samples=default_n_samples) + msg_out::Message, + msg_in1::Nothing, + msg_in2::Message{PointMass}; + dims::Any=nothing, + n_samples=default_n_samples) m_in2 = msg_in2.dist.params[:m] - return Message(variateType(dims), Function, log_pdf = (z)->logPdf(msg_out.dist, g(z, m_in2))) + return Message(variateType(dims), Function, log_pdf=(z) -> logPdf(msg_out.dist, g(z, m_in2))) end function ruleMDeltaSInMGX(g::Function, - msg_out::Message, - msgs_in::Vararg{Message{<:Gaussian}}) + msg_out::Message, + msgs_in::Vararg{Message{<:Gaussian}}) # Extract joint statistics of inbound messages (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Return arrays with individual means and covariances @@ -173,15 +173,15 @@ function collectSumProductNodeInbounds(node::Delta{Sampling}, entry::ScheduleEnt # Push function to calling signature # This function needs to be defined in the scope of the user - push!(inbounds, Dict{Symbol, Any}(:g => node.g, - :keyword => false)) + push!(inbounds, Dict{Symbol,Any}(:g => node.g, + :keyword => false)) multi_in = isMultiIn(node) # Boolean to indicate a Delta node with multiple stochastic inbounds inx = findfirst(isequal(entry.interface), node.interfaces) - 1 # Find number of inbound interface; 0 for outbound - + if (inx > 0) && multi_in # Multi-inbound backward rule - push!(inbounds, Dict{Symbol, Any}(:inx => inx, # Push inbound identifier - :keyword => false)) + push!(inbounds, Dict{Symbol,Any}(:inx => inx, # Push inbound identifier + :keyword => false)) end interface_to_schedule_entry = current_inference_algorithm.interface_to_schedule_entry @@ -205,12 +205,12 @@ function collectSumProductNodeInbounds(node::Delta{Sampling}, entry::ScheduleEnt # Push custom arguments if defined if (node.dims !== nothing) - push!(inbounds, Dict{Symbol, Any}(:dims => node.dims[inx + 1], - :keyword => true)) + push!(inbounds, Dict{Symbol,Any}(:dims => node.dims[inx+1], + :keyword => true)) end if (node.n_samples !== nothing) - push!(inbounds, Dict{Symbol, Any}(:n_samples => node.n_samples, - :keyword => true)) + push!(inbounds, Dict{Symbol,Any}(:n_samples => node.n_samples, + :keyword => true)) end return inbounds end @@ -221,19 +221,19 @@ end #--------------------------- function prod!( - x::Distribution{V, Function}, - y::Distribution{V, Function}) where V<:VariateType # log-pdf for z cannot be predefined, because it cannot be overwritten + x::Distribution{V,Function}, + y::Distribution{V,Function}) where V<:VariateType # log-pdf for z cannot be predefined, because it cannot be overwritten - return Distribution(V, Function, log_pdf=(s)->x.params[:log_pdf](s) + y.params[:log_pdf](s)) + return Distribution(V, Function, log_pdf=(s) -> x.params[:log_pdf](s) + y.params[:log_pdf](s)) end @symmetrical function prod!( x::Distribution{Univariate}, # Includes function distributions - y::Distribution{Univariate, <:Gaussian}, - z::Distribution{Univariate, Gaussian{Precision}}=Distribution(Univariate, Gaussian{Precision}, m=0.0, w=1.0)) + y::Distribution{Univariate,<:Gaussian}, + z::Distribution{Univariate,Gaussian{Precision}}=Distribution(Univariate, Gaussian{Precision}, m=0.0, w=1.0)) # Optimize with gradient ascent - log_joint(s) = logPdf(y,s) + logPdf(x,s) + log_joint(s) = logPdf(y, s) + logPdf(x, s) d_log_joint(s) = ForwardDiff.derivative(log_joint, s) m_initial = unsafeMean(y) @@ -248,11 +248,11 @@ end @symmetrical function prod!( x::Distribution{Multivariate}, # Includes function distributions - y::Distribution{Multivariate, <:Gaussian}, - z::Distribution{Multivariate, Gaussian{Precision}}=Distribution(Multivariate, Gaussian{Precision}, m=[0.0], w=mat(1.0))) + y::Distribution{Multivariate,<:Gaussian}, + z::Distribution{Multivariate,Gaussian{Precision}}=Distribution(Multivariate, Gaussian{Precision}, m=[0.0], w=mat(1.0))) # Optimize with gradient ascent - log_joint(s) = logPdf(y,s) + logPdf(x,s) + log_joint(s) = logPdf(y, s) + logPdf(x, s) d_log_joint(s) = ForwardDiff.gradient(log_joint, s) m_initial = unsafeMean(y) @@ -278,31 +278,35 @@ function gradientOptimization(log_joint::Function, d_log_joint::Function, m_init m_old = m_initial satisfied = false step_count = 0 - m_latests = if (dim_tot == 1) Queue{Float64}() else Queue{Vector}() end + m_latests = if (dim_tot == 1) + Queue{Float64}() + else + Queue{Vector}() + end while !satisfied - m_new = m_old .+ step_size.*d_log_joint(m_old) + m_new = m_old .+ step_size .* d_log_joint(m_old) if log_joint(m_new) > log_joint(m_old) - proposal_step_size = 10*step_size - m_proposal = m_old .+ proposal_step_size.*d_log_joint(m_old) + proposal_step_size = 10 * step_size + m_proposal = m_old .+ proposal_step_size .* d_log_joint(m_old) if log_joint(m_proposal) > log_joint(m_new) m_new = m_proposal step_size = proposal_step_size end else - step_size = 0.1*step_size - m_new = m_old .+ step_size.*d_log_joint(m_old) + step_size = 0.1 * step_size + m_new = m_old .+ step_size .* d_log_joint(m_old) end step_count += 1 - enqueue!(m_latests, m_old) + push!(m_latests, m_old) if step_count > 10 - m_average = sum(x for x in m_latests)./10 - if sum(sqrt.(((m_new.-m_average)./m_average).^2)) < dim_tot*0.1 + m_average = sum(x for x in m_latests) ./ 10 + if sum(sqrt.(((m_new .- m_average) ./ m_average) .^ 2)) < dim_tot * 0.1 satisfied = true end - dequeue!(m_latests); + popfirst!(m_latests) end - if step_count > dim_tot*250 + if step_count > dim_tot * 250 satisfied = true end m_old = m_new @@ -317,7 +321,7 @@ end #-------- function logJointPdfs(m_fw_in::Vector, W_fw_in::AbstractMatrix, dist_out::Distribution, g::Function, ds::Vector) - log_joint(x) = -0.5*sum(intdim.(ds))*log(2pi) + 0.5*logdet(W_fw_in) - 0.5*(x - m_fw_in)'*W_fw_in*(x - m_fw_in) + logPdf(dist_out, g(split(x, ds)...)) + log_joint(x) = -0.5 * sum(intdim.(ds)) * log(2pi) + 0.5 * logdet(W_fw_in) - 0.5 * (x - m_fw_in)' * W_fw_in * (x - m_fw_in) + logPdf(dist_out, g(split(x, ds)...)) d_log_joint(x) = ForwardDiff.gradient(log_joint, x) return (log_joint, d_log_joint)