Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Logpdf #14

Merged
merged 3 commits into from
Mar 26, 2024
Merged

Logpdf #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b9f1eae
added macros to reduce amount of duplicate code, logpdf introduction
benjione Mar 25, 2024
70a51b7
first version of logpdf
benjione Mar 26, 2024
1b18943
fix some tests
benjione Mar 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 0 additions & 3 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,6 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
Expand Down
94 changes: 0 additions & 94 deletions src/MCMC/mcmc.jl
View file
Open in desktop

This file was deleted.

4 changes: 2 additions & 2 deletions src/Samplers/algorithms.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -282,7 +282,7 @@ end
function SelfReinforced_ML_estimation(
X::PSDDataVector{T},
model::PSDModelOrthonormal{dsub, T, S},
bridge::BridgingDensity{dsub, T},
bridge::BridgingDensity{d, T},
reference_map::ReferenceMap{d, <:Any, T};
subsample_data=false,
subsample_size=2000,
Expand Down Expand Up @@ -413,7 +413,7 @@ function Adaptive_Self_reinforced_ML_estimation(
@assert _d == d
@assert typeof(reference_map) <: ReferenceMap{d, dC, T}

bridge = BridgingDensities.DiffusionBrigdingDensity{dsub, T}()
bridge = BridgingDensities.DiffusionBrigdingDensity{d, T}()

if dsub < d
@assert subspace_reference_map !== nothing
Expand Down
16 changes: 8 additions & 8 deletions src/Samplers/mappings/MarginalMapping.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,24 +44,24 @@ function inverse_Jacobian(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
return inverse_Jacobian(m.map, u[m.subvariables])
end

function marg_pushforward(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
function marginal_pushforward(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
u::PSDdata{T}) where {d,dC,T<:Number,dsub,dCsub,Mtype<:ConditionalMapping{dsub,dCsub,T}}
u[m.subvariables[1:(dsub-dCsub)]] = marg_pushforward(m.map, u[m.subvariables[1:(dsub-dCsub)]])
u[m.subvariables[1:(dsub-dCsub)]] = marginal_pushforward(m.map, u[m.subvariables[1:(dsub-dCsub)]])
return u
end

function marg_pullback(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
function marginal_pullback(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
u::PSDdata{T}) where {d,dC,T<:Number,dsub,dCsub,Mtype<:ConditionalMapping{dsub,dCsub,T}}
u[m.subvariables[1:(dsub-dCsub)]] = marg_pullback(m.map, u[m.subvariables[1:(dsub-dCsub)]])
u[m.subvariables[1:(dsub-dCsub)]] = marginal_pullback(m.map, u[m.subvariables[1:(dsub-dCsub)]])
return u
end

function marg_Jacobian(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
function marginal_Jacobian(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
x::PSDdata{T}) where {d,dC,T<:Number,dsub,dCsub,Mtype<:ConditionalMapping{dsub,dCsub,T}}
return marg_Jacobian(m.map, x[m.subvariables[1:(dsub-dCsub)]])
return marginal_Jacobian(m.map, x[m.subvariables[1:(dsub-dCsub)]])
end

function marg_inverse_Jacobian(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
function marginal_inverse_Jacobian(m::MarginalMapping{d,dC,T,dsub,dCsub,Mtype},
u::PSDdata{T}) where {d,dC,T<:Number,dsub,dCsub,Mtype<:ConditionalMapping{dsub,dCsub,T}}
return marg_inverse_Jacobian(m.map, u[m.subvariables[1:(dsub-dCsub)]])
return marginal_inverse_Jacobian(m.map, u[m.subvariables[1:(dsub-dCsub)]])
end
101 changes: 66 additions & 35 deletions src/Samplers/mappings/ProjectionMapping.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,30 +74,47 @@ Hence, det(∂v (P * T(v) + (1-P) v)) = det(∂v T(v)) * det(∂v v)
Total expression of forward is with T being the sampler:
R^{-1}(y) = X * R_sub^{-1}(T(R_sub(X' R^{-1}(x)))) + (I - P) * R^{-1}(x)
"""
function Distributions.pdf(sampler::ProjectionMapping{<:Any, <:Any, T, <:Any, <:Any, samplerType, R, R_sub},
x::PSDdata{T}
) where {T<:Number, R, R_sub, samplerType}


function inverse_Jacobian(sampler::ProjectionMapping{<:Any, <:Any, T, R, R_sub, <:Any, <:Any, ST},
x::PSDdata{T}
) where {T<:Number, R, R_sub, ST}
xs = pullback(sampler.R_map, x) # from [0,1]^d to R^d
sub_x = sampler.P_tilde * xs # from R^d to R^d2
sub_x_maped = pushforward(sampler.R_map_sub, sub_x) # from R^d2 to [0,1]^d2
sub_u = pullback(sampler.sampler, sub_x_maped) # still in [0,1]^d2
us = sampler.X * pullback(sampler.R_map_sub, sub_u) + (I - sampler.P) * xs
# u = pushforward(sampler.R_map, us)
T_inv_jac_det = Distributions.pdf(sampler.sampler, sub_x_maped) *
T_inv_jac_det = inverse_Jacobian(sampler.sampler, sub_x_maped) *
inverse_Jacobian(sampler.R_map_sub, sub_u) *
Jacobian(sampler.R_map_sub, sub_x)

return inverse_Jacobian(sampler.R_map, x) * T_inv_jac_det * Jacobian(sampler.R_map, us)
end
Jacobian(sampler::ProjectionMapping{<:Any, <:Any, T, R, R_sub, <:Any, <:Any, ST},
x::PSDdata{T}
) where {T<:Number, R, R_sub, ST} = 1.0/inverse_Jacobian(sampler, pushforward(sampler, x))

function inverse_Jacobian(sampler::ProjectionMapping{<:Any, <:Any, T, R, R_sub, <:Any, <:Any, ST},
function inverse_log_Jacobian(sampler::ProjectionMapping{<:Any, <:Any, T, R, R_sub, <:Any, <:Any, ST},
x::PSDdata{T}
) where {T<:Number, R, R_sub, ST}
return Distributions.pdf(sampler, x)
xs = pullback(sampler.R_map, x) # from [0,1]^d to R^d
sub_x = sampler.P_tilde * xs # from R^d to R^d2
sub_x_maped = pushforward(sampler.R_map_sub, sub_x) # from R^d2 to [0,1]^d2
sub_u = pullback(sampler.sampler, sub_x_maped) # still in [0,1]^d2
us = sampler.X * pullback(sampler.R_map_sub, sub_u) + (I - sampler.P) * xs
# u = pushforward(sampler.R_map, us)
T_inv_jac_det = inverse_log_Jacobian(sampler.sampler, sub_x_maped) +
inverse_log_Jacobian(sampler.R_map_sub, sub_u) +
log_Jacobian(sampler.R_map_sub, sub_x)

return inverse_log_Jacobian(sampler.R_map, x) + T_inv_jac_det + log_Jacobian(sampler.R_map, us)
end
Jacobian(sampler::ProjectionMapping{<:Any, <:Any, T, R, R_sub, <:Any, <:Any, ST},
function log_Jacobian(sampler::ProjectionMapping{<:Any, <:Any, T, R, R_sub, <:Any, <:Any, ST},
x::PSDdata{T}
) where {T<:Number, R, R_sub, ST} = 1.0/inverse_Jacobian(sampler, pushforward(sampler, x))
) where {T<:Number, R, R_sub, ST}
return -inverse_log_Jacobian(sampler, pushforward(sampler, x))
end

###
## Helper functions
Expand Down Expand Up @@ -179,54 +196,68 @@ end
# ConditionalSampler
########################

"""
Distribution p(x) = ∫ p(x, y) d y
"""
function marg_pdf(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
x::PSDdata{T}) where {d, dC, T<:Number, dsub, dCsub, ST, R1, R2}
xs = marg_pullback(sampler.R_map, x) # from [0,1]^dx to R^dx
function marginal_inverse_Jacobian(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
x::PSDdata{T}
) where {d, dC, T<:Number, dsub, dCsub, ST, R1, R2}
xs = marginal_pullback(sampler.R_map, x) # from [0,1]^dx to R^dx
sub_x = _marg_P_tilde(sampler) * xs # from R^dx to R^d_sub_marg
sub_x_maped = marginal_pushforward(sampler.R_map_sub, sub_x) # from R^d_sub_marg to [0,1]^d_sub_marg
sub_u = marginal_pullback(sampler.sampler, sub_x_maped) # still in [0,1]^d_sub_marg
us = _marg_X(sampler) * marginal_pullback(sampler.R_map_sub, sub_u) + (I - _marg_P(sampler)) * xs
# u = pushforward(sampler.R_map, us)
T_inv_jac_det = marginal_inverse_Jacobian(sampler.sampler, sub_x_maped) *
marginal_inverse_Jacobian(sampler.R_map_sub, sub_u) *
marginal_Jacobian(sampler.R_map_sub, sub_x)

return marginal_inverse_Jacobian(sampler.R_map, x) * T_inv_jac_det * marginal_Jacobian(sampler.R_map, us)
end

function marginal_inverse_log_Jacobian(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
x::PSDdata{T}
) where {d, dC, T<:Number, dsub, dCsub, ST, R1, R2}
xs = marginal_pullback(sampler.R_map, x) # from [0,1]^dx to R^dx
sub_x = _marg_P_tilde(sampler) * xs # from R^dx to R^d_sub_marg
sub_x_maped = marg_pushforward(sampler.R_map_sub, sub_x) # from R^d_sub_marg to [0,1]^d_sub_marg
sub_u = marg_pullback(sampler.sampler, sub_x_maped) # still in [0,1]^d_sub_marg
us = _marg_X(sampler) * marg_pullback(sampler.R_map_sub, sub_u) + (I - _marg_P(sampler)) * xs
sub_x_maped = marginal_pushforward(sampler.R_map_sub, sub_x) # from R^d_sub_marg to [0,1]^d_sub_marg
sub_u = marginal_pullback(sampler.sampler, sub_x_maped) # still in [0,1]^d_sub_marg
us = _marg_X(sampler) * marginal_pullback(sampler.R_map_sub, sub_u) + (I - _marg_P(sampler)) * xs
# u = pushforward(sampler.R_map, us)
T_inv_jac_det = marg_pdf(sampler.sampler, sub_x_maped) *
marg_inverse_Jacobian(sampler.R_map_sub, sub_u) *
marg_Jacobian(sampler.R_map_sub, sub_x)
T_inv_jac_det = marginal_inverse_log_Jacobian(sampler.sampler, sub_x_maped) +
marginal_inverse_log_Jacobian(sampler.R_map_sub, sub_u) +
marginal_log_Jacobian(sampler.R_map_sub, sub_x)

return marg_inverse_Jacobian(sampler.R_map, x) * T_inv_jac_det * marg_Jacobian(sampler.R_map, us)
return marginal_inverse_log_Jacobian(sampler.R_map, x) + T_inv_jac_det + marginal_log_Jacobian(sampler.R_map, us)
end

function marg_inverse_Jacobian(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
function marginal_Jacobian(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
x::PSDdata{T}
) where {d, dC, T<:Number, dsub, dCsub, ST, R1, R2}
return marg_pdf(sampler, x)
return 1.0/marginal_inverse_Jacobian(sampler, marginal_pushforward(sampler, x))
end

function marg_Jacobian(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
function marginal_log_Jacobian(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
x::PSDdata{T}
) where {d, dC, T<:Number, dsub, dCsub, ST, R1, R2}
return 1.0/marg_inverse_Jacobian(sampler, marg_pushforward(sampler, x))
return -marginal_inverse_log_Jacobian(sampler, marginal_pushforward(sampler, x))
end

function marg_pushforward(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
function marginal_pushforward(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
u::PSDdata{T}
) where {d, dC, T<:Number, dsub, dCsub, ST, R1, R2}
us = marg_pullback(sampler.R_map, u)
us = marginal_pullback(sampler.R_map, u)
sub_u = _marg_P_tilde(sampler) * us
sub_x = marg_pushforward(sampler.sampler, marg_pushforward(sampler.R_map_sub, sub_u))
xs = _marg_X(sampler) * marg_pullback(sampler.R_map_sub, sub_x) + (I - _marg_P(sampler)) * us
x = marg_pushforward(sampler.R_map, xs)
sub_x = marginal_pushforward(sampler.sampler, marginal_pushforward(sampler.R_map_sub, sub_u))
xs = _marg_X(sampler) * marginal_pullback(sampler.R_map_sub, sub_x) + (I - _marg_P(sampler)) * us
x = marginal_pushforward(sampler.R_map, xs)
return x
end

function marg_pullback(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
function marginal_pullback(sampler::ProjectionMapping{d, dC, T, dsub, dCsub, ST, R1, R2},
x::PSDdata{T}
) where {d, dC, T<:Number, dsub, dCsub, ST, R1, R2}
xs = marg_pullback(sampler.R_map, x)
xs = marginal_pullback(sampler.R_map, x)
sub_x = _marg_P_tilde(sampler) * xs
sub_u = marg_pullback(sampler.sampler, marg_pushforward(sampler.R_map_sub, sub_x))
us = _marg_X(sampler) * marg_pullback(sampler.R_map_sub, sub_u) + (I - _marg_P(sampler)) * xs
u = marg_pushforward(sampler.R_map, us)
sub_u = marginal_pullback(sampler.sampler, marginal_pushforward(sampler.R_map_sub, sub_x))
us = _marg_X(sampler) * marginal_pullback(sampler.R_map_sub, sub_u) + (I - _marg_P(sampler)) * xs
u = marginal_pushforward(sampler.R_map, us)
return u
end
Loading
Loading