diff --git a/Project.toml b/Project.toml
index a3a2fae..9c4e093 100644
--- a/Project.toml
+++ b/Project.toml
@@ -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"
diff --git a/src/MCMC/mcmc.jl b/src/MCMC/mcmc.jl
deleted file mode 100644
index 6d2ef92..0000000
--- a/src/MCMC/mcmc.jl
+++ /dev/null
@@ -1,94 +0,0 @@
-module MCMC
-"""
-Implements MCMC based on the `AbstractMCMC` interface using the
-sampler from `Sampler` in this package.
-"""
-
-using ..SequentialMeasureTransport
-using ..SequentialMeasureTransport: AbstractSampler
-import AbstractMCMC
-import MCMCChains
-using Distributions
-using Random
-using LinearAlgebra
-
-mutable struct MCMCSampler{T, S<:AbstractSampler{<:Any, <:Any, T}} <: AbstractMCMC.AbstractSampler
-    current_state::Vector{T}
-    sampler::S
-end
-
-struct MCMCModel{F} <: AbstractMCMC.AbstractModel
-    model::F
-    function MCMCModel(model::F) where {F}
-        return new{F}(model)
-    end
-end
-function (a::MCMCModel{F})(x::Vector{T}) where {F, T<:Number}
-    return a.model(x)
-end
-function (a::MCMCModel{F})(x::Vector{T}, 
-                sampler::AbstractSampler{<:Any, <:Any, T}
-            ) where {F, T<:Number}
-    return SequentialMeasureTransport.pullback(sampler, x->a.model(x))(x)
-end
-
-function proposal(sampler::MCMCSampler{T}, x::Vector{T}) where {T}
-    return x + rand(MvNormal(zeros(length(sampler.current_state)), diagm(ones(length(sampler.current_state)))))
-end
-
-q(x::Vector{T}, y::Vector{T}) where {T<:Number} = pdf(MvNormal(zeros(length(x)), diagm(ones(length(x)))), x-y)
-
-function AbstractMCMC.step(rng::Random.AbstractRNG,
-        model::MCMCModel,
-        samp::MCMCSampler{T};
-        kwargs...) where {T <: Number}
-    next_θ = SequentialMeasureTransport.sample(samp.sampler)
-    return next_θ, next_θ
-end
-function AbstractMCMC.step(rng::Random.AbstractRNG,
-        model::MCMCModel,
-        samp::MCMCSampler{T},
-        θ_prev::Vector{T};
-        kwargs...) where {T <: Number}
-    # θ_prev = samp.current_state
-    # Generate a new proposal.
-    # θ = propose(spl, model, θ_prev)
-    θ = proposal(samp, θ_prev)
-
-    # Calculate the acceptance probability.
-    α = model(θ, samp.sampler) *  q(Θ, θ_prev) / (model(θ_prev, samp.sampler) * q(θ_prev, θ))
-
-    # Decide whether to return the previous θ or the new one.
-    next_Θ = if rand() < min(α, 1.0)
-        θ
-    else
-        θ_prev
-    end
-    samp.current_state = next_Θ
-    return next_Θ, next_Θ
-end
-
-# A basic chains constructor that works with the Transition struct we defined.
-function AbstractMCMC.bundle_samples(
-            ts::Vector{Vector{T}}, 
-            ℓ::MCMCModel, 
-            s::MCMCSampler{T}, 
-            state, 
-            chain_type::Type{Any};
-            param_names=missing,
-            kwargs...
-        ) where {T<:Number}
-    # Turn all the transitions into a vector-of-vectors.
-    vals = reduce(hcat,ts)'
-
-    # Check if we received any parameter names.
-    if ismissing(param_names)
-        param_names = String["Parameter $i" for i in 1:size(vals,2)]
-    end
-
-    # Bundle everything up and return a Chains struct.
-    return MCMCChains.Chains(vals, param_names)
-end
-
-
-end # module MCMC
\ No newline at end of file
diff --git a/src/Samplers/algorithms.jl b/src/Samplers/algorithms.jl
index 69485bf..bfa64e8 100644
--- a/src/Samplers/algorithms.jl
+++ b/src/Samplers/algorithms.jl
@@ -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,
@@ -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
diff --git a/src/Samplers/mappings/MarginalMapping.jl b/src/Samplers/mappings/MarginalMapping.jl
index a42cd27..08fabbb 100644
--- a/src/Samplers/mappings/MarginalMapping.jl
+++ b/src/Samplers/mappings/MarginalMapping.jl
@@ -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
\ No newline at end of file
diff --git a/src/Samplers/mappings/ProjectionMapping.jl b/src/Samplers/mappings/ProjectionMapping.jl
index 0c2f55e..754c722 100644
--- a/src/Samplers/mappings/ProjectionMapping.jl
+++ b/src/Samplers/mappings/ProjectionMapping.jl
@@ -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
@@ -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
diff --git a/src/Samplers/reference_maps/algebraic.jl b/src/Samplers/reference_maps/algebraic.jl
index ae8635d..26f61ac 100644
--- a/src/Samplers/reference_maps/algebraic.jl
+++ b/src/Samplers/reference_maps/algebraic.jl
@@ -11,13 +11,38 @@ struct AlgebraicReference{d, dC, T} <: ReferenceMap{d, dC, T}
     end
 end
 
+function _pushforward(m::AlgebraicReference{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number}
+    return ((x./sqrt.(1 .+ x.^2)).+1.0)/2.0
+end
+
+function _pullback(m::AlgebraicReference{<:Any, <:Any, T}, u::PSDdata{T}) where {T<:Number}
+    ξ = 2.0*(u .- 0.5)
+    return ξ./sqrt.(1.0 .- ξ.^2)
+end
+
+function _Jacobian(m::AlgebraicReference{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number}
+    return mapreduce(xi->0.5/(1+(xi)^2)^(3/2), *, x)
+end
+
+function _inverse_Jacobian(m::AlgebraicReference{<:Any, <:Any, T}, u::PSDdata{T}) where {T<:Number}
+    return mapreduce(ui->2.0/(1-((2*ui - 1.0)^2))^(3/2), *, u)
+end
+
+function _log_Jacobian(m::AlgebraicReference{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number}
+    return mapreduce(xi->log(0.5) - (3/2) * log(1+xi^2), +, x)
+end
+
+function _inverse_log_Jacobian(m::AlgebraicReference{<:Any, <:Any, T}, u::PSDdata{T}) where {T<:Number}
+    return mapreduce(ui->log(2.0) - (3/2) * log(1.0-((2*ui - 1.0)^2)), +, u)
+end
+
 
 function SMT.pushforward(
         m::AlgebraicReference{d, <:Any, T}, 
         x::PSDdata{T}
     ) where {d, T<:Number}
     @assert length(x) == d
-    return ((x./sqrt.(1 .+ x.^2)).+1.0)/2.0
+    return _pushforward(m, x)
 end
 
 
@@ -26,58 +51,88 @@ function SMT.pullback(
         u::PSDdata{T}
     ) where {d, T<:Number}
     @assert length(u) == d
-    ξ = 2.0*(u .- 0.5)
-    return ξ./sqrt.(1.0 .- ξ.^2)
+    return _pullback(m, u)
 end
 
-
 function SMT.Jacobian(
         m::AlgebraicReference{d, <:Any, T}, 
         x::PSDdata{T}
     ) where {d, T<:Number}
     @assert length(x) == d
-    return mapreduce(xi->0.5/(1+(xi)^2)^(3/2), *, x)
+    return _Jacobian(m, x)
 end
 
 
 function SMT.inverse_Jacobian(
+        m::AlgebraicReference{d, <:Any, T}, 
+        u::PSDdata{T}
+    ) where {d, T<:Number}
+    @assert length(u) == d
+    return _inverse_Jacobian(m, u)
+end
+
+function SMT.log_Jacobian(
+        m::AlgebraicReference{d, <:Any, T}, 
+        x::PSDdata{T}
+    ) where {d, T<:Number}
+    @assert length(x) == d
+    return _log_Jacobian(m, x)
+end
+
+function SMT.inverse_log_Jacobian(
         mapping::AlgebraicReference{d, <:Any, T}, 
         u::PSDdata{T}
     ) where {d, T<:Number}
     @assert length(u) == d
-    return mapreduce(ui->2.0/(1-((2*ui - 1.0)^2))^(3/2), *, u)
+    return _inverse_log_Jacobian(mapping, u)
 end
 
 
-function SMT.marg_pushforward(
+function SMT.marginal_pushforward(
         m::AlgebraicReference{d, dC, T}, 
         x::PSDdata{T}
     ) where {d, dC, T<:Number}
     @assert length(x) == d-dC
-    return ((x./sqrt.(1 .+ x.^2)).+1.0)/2.0
+    return _pushforward(m, x)
 end
 
-function SMT.marg_pullback(
+function SMT.marginal_pullback(
         m::AlgebraicReference{d, dC, T}, 
         u::PSDdata{T}
     ) where {d, dC, T<:Number}
     @assert length(u) == d-dC
     ξ = 2.0*(u .- 0.5)
-    return ξ./sqrt.(1.0 .- ξ.^2)
+    return _pullback(m, u)
 end
 
-function SMT.marg_Jacobian(
+function SMT.marginal_Jacobian(
         m::AlgebraicReference{d, dC, T}, 
         x::PSDdata{T}
     ) where {d, dC, T<:Number}
     @assert length(x) == d-dC
-    return mapreduce(xi->0.5/(1+(xi)^2)^(3/2), *, x)
+    return _Jacobian(m, x)
 end
 
-function SMT.marg_inverse_Jacobian(
+function SMT.marginal_inverse_Jacobian(
         mapping::AlgebraicReference{d, dC, T}, 
         u::PSDdata{T}
     ) where {d, dC, T<:Number}
     @assert length(u) == d-dC
-    return mapreduce(ui->2.0/(1-((2*ui - 1.0)^2))^(3/2), *, u)
+    return _inverse_Jacobian(mapping, u)
+end
+
+function SMT.marginal_log_Jacobian(
+        m::AlgebraicReference{d, dC, T}, 
+        x::PSDdata{T}
+    ) where {d, dC, T<:Number}
+    @assert length(x) == d-dC
+    return _log_Jacobian(m, x)
+end
+
+function SMT.marginal_inverse_log_Jacobian(
+        mapping::AlgebraicReference{d, dC, T}, 
+        u::PSDdata{T}
+    ) where {d, dC, T<:Number}
+    @assert length(u) == d-dC
+    return _inverse_log_Jacobian(mapping, u)
 end
\ No newline at end of file
diff --git a/src/Samplers/reference_maps/gaussian.jl b/src/Samplers/reference_maps/gaussian.jl
index f90f30c..af0d101 100644
--- a/src/Samplers/reference_maps/gaussian.jl
+++ b/src/Samplers/reference_maps/gaussian.jl
@@ -67,7 +67,23 @@ function SMT.inverse_Jacobian(
     return 1/Jacobian(mapping, pullback(mapping, u))
 end
 
-function SMT.marg_pushforward(
+function SMT.log_Jacobian(
+        m::GaussianReference{d, <:Any, T}, 
+        x::PSDdata{T}
+    ) where {d, T<:Number}
+    @assert length(x) == d
+    mapreduce(xi->Distributions.logpdf(Distributions.Normal(0, m.σ), xi), +, x)
+end
+
+function SMT.inverse_log_Jacobian(
+        mapping::GaussianReference{d, <:Any, T}, 
+        u::PSDdata{T}
+    ) where {d, T<:Number}
+    # inverse function theorem
+    return -SMT.log_Jacobian(mapping, SMT.pullback(mapping, u))
+end
+
+function SMT.marginal_pushforward(
         m::GaussianReference{d, dC, T}, 
         x::PSDdata{T}
     ) where {d, dC, T<:Number}
@@ -75,7 +91,7 @@ function SMT.marg_pushforward(
     return 0.5 * (1 .+ erf.(x ./ (m.σ * sqrt(2))))
 end
 
-function SMT.marg_pullback(
+function SMT.marginal_pullback(
         m::GaussianReference{d, dC, T}, 
         u::PSDdata{T}
     ) where {d, dC, T<:Number}
@@ -83,7 +99,7 @@ function SMT.marg_pullback(
     return sqrt(2) * m.σ * erfcinv.(2.0 .- 2*u)
 end
 
-function SMT.marg_Jacobian(
+function SMT.marginal_Jacobian(
         m::GaussianReference{d, dC, T}, 
         x::PSDdata{T}
     ) where {d, dC, T<:Number}
@@ -91,10 +107,26 @@ function SMT.marg_Jacobian(
     mapreduce(xi->Distributions.pdf(Distributions.Normal(0, m.σ), xi), *, x)
 end
 
-function SMT.marg_inverse_Jacobian(
+function SMT.marginal_inverse_Jacobian(
+        mapping::GaussianReference{<:Any, <:Any, T}, 
+        u::PSDdata{T}
+    ) where {T<:Number}
+    # inverse function theorem
+    return 1/SMT.marginal_Jacobian(mapping, SMT.marginal_pullback(mapping, u))
+end
+
+function SMT.marginal_log_Jacobian(
+        m::GaussianReference{d, dC, T}, 
+        x::PSDdata{T}
+    ) where {d, dC, T<:Number}
+    @assert length(x) == d-dC
+    mapreduce(xi->Distributions.logpdf(Distributions.Normal(0, m.σ), xi), +, x)
+end
+
+function SMT.marginal_inverse_log_Jacobian(
         mapping::GaussianReference{<:Any, <:Any, T}, 
         u::PSDdata{T}
     ) where {T<:Number}
     # inverse function theorem
-    return 1/SMT.marg_Jacobian(mapping, SMT.marg_pullback(mapping, u))
+    return -SMT.marginal_log_Jacobian(mapping, SMT.marginal_pullback(mapping, u))
 end
\ No newline at end of file
diff --git a/src/Samplers/reference_maps/reference_maps.jl b/src/Samplers/reference_maps/reference_maps.jl
index 8cab8b7..f52c562 100644
--- a/src/Samplers/reference_maps/reference_maps.jl
+++ b/src/Samplers/reference_maps/reference_maps.jl
@@ -87,28 +87,28 @@ function SMT.inverse_Jacobian(
     throw(error("Not implemented"))
 end
 
-function SMT.marg_pushforward(
+function SMT.marginal_pushforward(
         mapping::ReferenceMap{d, dC, T}, 
         x::PSDdata{T}
     ) where{d, dC, T<:Number}
     throw(error("Not implemented"))
 end
 
-function SMT.marg_pullback(
+function SMT.marginal_pullback(
         mapping::ReferenceMap{d, dC, T}, 
         u::PSDdata{T}
     ) where {d, dC, T<:Number}
     throw(error("Not implemented"))
 end
 
-function SMT.marg_Jacobian(
+function SMT.marginal_Jacobian(
         mapping::ReferenceMap{d, dC, T}, 
         x::PSDdata{T}
     ) where {d, dC, T<:Number}
     throw(error("Not implemented"))
 end
 
-function SMT.marg_inverse_Jacobian(
+function SMT.marginal_inverse_Jacobian(
         mapping::ReferenceMap{d, dC, T}, 
         u::PSDdata{T}
     ) where {d, dC, T<:Number}
diff --git a/src/Samplers/reference_maps/scaling.jl b/src/Samplers/reference_maps/scaling.jl
index 85222e9..3eb67db 100644
--- a/src/Samplers/reference_maps/scaling.jl
+++ b/src/Samplers/reference_maps/scaling.jl
@@ -70,7 +70,7 @@ function SMT.inverse_Jacobian(
     return mapping.V
 end
 
-function SMT.marg_pushforward(
+function SMT.marginal_pushforward(
         m::ScalingReference{d, dC, T}, 
         x::PSDdata{T}
     ) where {d, dC, T<:Number}
@@ -78,7 +78,7 @@ function SMT.marg_pushforward(
     return (x .- m.L[1:dC]) ./ (m.R[1:dC] .- m.L[1:dC])
 end
 
-function SMT.marg_pullback(
+function SMT.marginal_pullback(
         m::ScalingReference{d, dC, T}, 
         u::PSDdata{T}
     ) where {d, dC, T<:Number}
@@ -86,7 +86,7 @@ function SMT.marg_pullback(
     return u .* (m.R[1:dC] .- m.L[1:dC]) .+ m.L[1:dC]
 end
 
-function SMT.marg_Jacobian(
+function SMT.marginal_Jacobian(
         mapping::ScalingReference{d, dC, T}, 
         x::PSDdata{T}
     ) where {d, dC, T<:Number}
@@ -94,7 +94,7 @@ function SMT.marg_Jacobian(
     1/mapping.V_margin
 end
 
-function SMT.marg_inverse_Jacobian(
+function SMT.marginal_inverse_Jacobian(
         mapping::ScalingReference{d, dC, T}, 
         u::PSDdata{T}
     ) where {d, dC, T<:Number}
diff --git a/src/Samplers/sampler.jl b/src/Samplers/sampler.jl
index 3acd152..0bcb8ce 100644
--- a/src/Samplers/sampler.jl
+++ b/src/Samplers/sampler.jl
@@ -34,6 +34,8 @@ pushforward(sampler::ConditionalMapping, u::PSDdata) = throw("Not Implemented")
 pullback(sampler::ConditionalMapping, x::PSDdata) = throw("Not Implemented")
 Jacobian(sampler::ConditionalMapping, x::PSDdata) = throw("Not Implemented")
 inverse_Jacobian(sampler::ConditionalMapping, u::PSDdata) = throw("Not Implemented")
+log_Jacobian(sampler::ConditionalMapping, x::PSDdata) = log(Jacobian(sampler, x))
+inverse_log_Jacobian(sampler::ConditionalMapping, u::PSDdata) = log(inverse_Jacobian(sampler, u))
 function pushforward(sampler::ConditionalMapping{d,<:Any,T}, π::Function) where {d,T<:Number}
     π_pushed = let sampler = sampler, π = π
         (x) -> π(pullback(sampler, x)) * inverse_Jacobian(sampler, x)
@@ -46,30 +48,52 @@ function pullback(sampler::ConditionalMapping{d,<:Any,T}, π::Function) where {d
     end
     return π_pulled
 end
+## log pushforward functions
+function log_pushforward(sampler::ConditionalMapping{d,<:Any,T}, log_π::Function) where {d,T<:Number}
+    π_pushed = let sampler = sampler, log_π = log_π
+        (x) -> log_π(pullback(sampler, x)) + inverse_log_Jacobian(sampler, x)
+    end
+    return π_pushed
+end
+function log_pullback(sampler::ConditionalMapping{d,<:Any,T}, log_π::Function) where {d,T<:Number}
+    π_pulled = let sampler = sampler, log_π = log_π
+        (u) -> log_π(pushforward(sampler, u)) + log_Jacobian(sampler, u)
+    end
+    return π_pulled
+end
 
 ## Methods interface for ConditionalMapping
 
-marg_pushforward(sampler::ConditionalMapping, u::PSDdata) = throw("Not Implemented")
-marg_pullback(sampler::ConditionalMapping, x::PSDdata) = throw("Not Implemented")
-marg_Jacobian(mapping::ConditionalMapping, u::PSDdata) = throw("Not Implemented")
-marg_inverse_Jacobian(mapping::ConditionalMapping, x::PSDdata) = throw("Not Implemented")
-function cond_Jacobian(sampler::ConditionalMapping{<:Any, <:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {T}
-    x = marg_pullback(sampler, x)
-    return Jacobian(sampler, [x; y]) / marg_Jacobian(sampler, x)
+marginal_pushforward(sampler::ConditionalMapping, u::PSDdata) = throw("Not Implemented")
+marginal_pullback(sampler::ConditionalMapping, x::PSDdata) = throw("Not Implemented")
+marginal_Jacobian(mapping::ConditionalMapping, u::PSDdata) = throw("Not Implemented")
+marginal_inverse_Jacobian(mapping::ConditionalMapping, x::PSDdata) = throw("Not Implemented")
+marginal_log_Jacobian(mapping::ConditionalMapping, u::PSDdata) = log(marginal_Jacobian(mapping, u))
+marginal_inverse_log_Jacobian(mapping::ConditionalMapping, x::PSDdata) = log(marginal_inverse_Jacobian(mapping, x))
+function conditional_Jacobian(sampler::ConditionalMapping{<:Any, <:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {T}
+    x = marginal_pullback(sampler, x)
+    return Jacobian(sampler, [x; y]) / marginal_Jacobian(sampler, x)
+end
+function conditional_inverse_Jacobian(sampler::ConditionalMapping{<:Any, <:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {T}
+    return inverse_Jacobian(sampler, [x; y]) / marginal_inverse_Jacobian(sampler, x)
 end
-function cond_inverse_Jacobian(sampler::ConditionalMapping{<:Any, <:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {T}
-    return inverse_Jacobian(sampler, [x; y]) / marg_inverse_Jacobian(sampler, x)
+function conditional_log_Jacobian(sampler::ConditionalMapping{<:Any, <:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {T}
+    x = marginal_pullback(sampler, x)
+    return log_Jacobian(sampler, [x; y]) - marginal_log_Jacobian(sampler, x)
+end
+function conditional_inverse_log_Jacobian(sampler::ConditionalMapping{<:Any, <:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {T}
+    return log_inverse_Jacobian(sampler, [x; y]) - marginal_inverse_log_Jacobian(sampler, x)
 end
 
 """
 Pullback of a conditional mapping.
 x is from x ∼ π_x
 """
-function cond_pushforward(sampler::ConditionalMapping{d,dC,T},
+function conditional_pushforward(sampler::ConditionalMapping{d,dC,T},
     u::PSDdata{T},
     x::PSDdata{T}
 ) where {d,T,dC}
-    x = marg_pullback(sampler, x)
+    x = marginal_pullback(sampler, x)
     xu = pushforward(sampler, [x; u])
     return xu[_d_marg(sampler)+1:d]
 end
@@ -77,53 +101,96 @@ end
 Pullback of a conditional mapping.
 x is from x ∼ π_x
 """
-function cond_pullback(sra::ConditionalMapping{d,dC,T},
+function conditional_pullback(sra::ConditionalMapping{d,dC,T},
     y::PSDdata{T},
     x::PSDdata{T}
 ) where {d,T<:Number,dC}
     yx = pullback(sra, [x; y])
     return yx[_d_marg(sra)+1:end]
 end
-function marg_pullback(sampler::ConditionalMapping{d,dC,T}, π::Function) where {d,dC,T}
+
+##
+## pushforward and pullback of functions
+##
+
+function marginal_pullback(sampler::ConditionalMapping{d,dC,T}, π::Function) where {d,dC,T}
     π_pb = let sampler = sampler, π = π
         (x) -> begin
-            π(marg_pushforward(sampler, x)) * marg_Jacobian(sampler, x)
+            π(marginal_pushforward(sampler, x)) * marginal_Jacobian(sampler, x)
         end
     end
     return π_pb
 end
-function marg_pushforward(sampler::ConditionalMapping{d,dC,T}, π::Function) where {d,dC,T}
+function marginal_pushforward(sampler::ConditionalMapping{d,dC,T}, π::Function) where {d,dC,T}
     π_pf = let sampler = sampler, π = π
         (u) -> begin
-            π(marg_pullback(sampler, u)) * marg_inverse_Jacobian(sampler, u)
+            π(marginal_pullback(sampler, u)) * marginal_inverse_Jacobian(sampler, u)
         end
     end
     return π_pf
 end
-function cond_pushforward(sampler::ConditionalMapping{d,dC,T},
+
+function marginal_log_pushforward(sampler::ConditionalMapping{d,dC,T}, log_π::Function) where {d,dC,T}
+    log_π_pf = let sampler = sampler, log_π = log_π
+        (u) -> begin
+            log_π(marginal_pullback(sampler, u)) + marginal_inverse_log_Jacobian(sampler, u)
+        end
+    end
+    return log_π_pf
+end
+function marginal_log_pullback(sampler::ConditionalMapping{d,dC,T}, log_π::Function) where {d,dC,T}
+    log_π_pb = let sampler = sampler, log_π = log_π
+        (x) -> begin
+            log_π(marginal_pushforward(sampler, x)) + marginal_log_Jacobian(sampler, x)
+        end
+    end
+    return log_π_pb
+end
+
+function conditional_pushforward(sampler::ConditionalMapping{d,dC,T},
     π::Function,
     x::PSDdata{T}
 ) where {d,T<:Number,dC}
     π_pf = let sampler = sampler, π = π, x = x
         (y::PSDdata{T}) -> begin
-            π(cond_pullback(sampler, y, x)) * cond_inverse_Jacobian(sampler, y, x)
+            π(conditional_pullback(sampler, y, x)) * conditional_inverse_Jacobian(sampler, y, x)
         end
     end
     return π_pf
 end
-
-function cond_pullback(sampler::ConditionalMapping{d,dC,T},
+function conditional_pullback(sampler::ConditionalMapping{d,dC,T},
     π::Function,
     x::PSDdata{T}
 ) where {d,T<:Number,dC}
     π_pb = let sampler = sampler, π = π, x = x
         (y::PSDdata{T}) -> begin
-            π(cond_pushforward(sampler, y, x)) * cond_Jacobian(sampler, y, x)
+            π(conditional_pushforward(sampler, y, x)) * conditional_Jacobian(sampler, y, x)
         end
     end
     return π_pb
 end
-
+function conditional_log_pushforward(sampler::ConditionalMapping{d,dC,T},
+    log_π::Function,
+    x::PSDdata{T}
+) where {d,T<:Number,dC}
+    log_π_pf = let sampler = sampler, log_π = log_π, x = x
+        (y::PSDdata{T}) -> begin
+            log_π(conditional_pullback(sampler, y, x)) + conditional_inverse_log_Jacobian(sampler, y, x)
+        end
+    end
+    return log_π_pf
+end
+function conditional_log_pullback(sampler::ConditionalMapping{d,dC,T},
+    log_π::Function,
+    x::PSDdata{T}
+) where {d,T<:Number,dC}
+    log_π_pb = let sampler = sampler, log_π = log_π, x = x
+        (y::PSDdata{T}) -> begin
+            log_π(conditional_pushforward(sampler, y, x)) + conditional_log_Jacobian(sampler, y, x)
+        end
+    end
+    return log_π_pb
+end
 
 # include reference maps
 include("reference_maps/reference_maps.jl")
@@ -154,6 +221,7 @@ end
 
 ## Interface for Sampler
 Distributions.pdf(sampler::AbstractCondSampler, x::PSDdata) = throw(NotImplementedError())
+Distributions.logpdf(sampler::AbstractCondSampler, x::PSDdata) = throw(NotImplementedError())
 
 ## methods not necessarily implemented by concrete implementations of Sampler:
 Base.rand(sampler::AbstractCondSampler) = sample(sampler)
@@ -182,14 +250,15 @@ end
 """
 Distribution p(x) = ∫ p(x, y) d y
 """
-marg_pdf(sampler::AbstractCondSampler, x::PSDdata) = throw(NotImplementedError())
-function marg_sample(sampler::AbstractCondSampler{d, dC, T}) where {d, dC, T<:Number}
-    return marg_pushforward(sampler, sample_reference(sampler)[1:_d_marg(sampler)])
+marginal_pdf(sampler::AbstractCondSampler, x::PSDdata) = throw(NotImplementedError())
+marginal_logpdf(sampler::AbstractCondSampler, x::PSDdata) = throw(NotImplementedError())
+function marginal_sample(sampler::AbstractCondSampler{d, dC, T}) where {d, dC, T<:Number}
+    return marginal_pushforward(sampler, sample_reference(sampler)[1:_d_marg(sampler)])
 end
-function marg_sample(sampler::AbstractCondSampler{<:Any,<:Any,T}, amount::Int; threading=true) where {T}
+function marginal_sample(sampler::AbstractCondSampler{<:Any,<:Any,T}, amount::Int; threading=true) where {T}
     res = Vector{Vector{T}}(undef, amount)
     @_condusethreads threading for i = 1:amount
-        res[i] = marg_sample(sampler)
+        res[i] = marginal_sample(sampler)
     end
     return res
 end
@@ -197,27 +266,32 @@ end
 """
 PDF p(y|x) = p(x, y) / p(x)
 """
-function cond_pdf(sampler::AbstractCondSampler{d,<:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {d,T}
-    return Distributions.pdf(sampler, [x; y]) / marg_pdf(sampler, x)
+function conditional_pdf(sampler::AbstractCondSampler{d,<:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {d,T}
+    # almost always better to use logpdf and exp
+    return exp(conditional_logpdf(sampler, y, x))
+    # return Distributions.pdf(sampler, [x; y]) / marginal_pdf(sampler, x)
+end
+function conditional_logpdf(sampler::AbstractCondSampler{d,<:Any,T}, y::PSDdata{T}, x::PSDdata{T}) where {d,T}
+    return Distributions.logpdf(sampler, [x; y]) - marginal_logpdf(sampler, x)
 end
 
-function cond_sample(sampler::AbstractCondSampler{d,<:Any,T},
+function conditional_sample(sampler::AbstractCondSampler{d,<:Any,T},
     X::PSDDataVector{T};
     threading=true
 ) where {d,T<:Number}
     if threading == false
-        return PSDdata{T}[cond_sample(sampler, x) for x in X]
+        return PSDdata{T}[conditional_sample(sampler, x) for x in X]
     else
         res = Vector{PSDdata{T}}(undef, length(X))
         Threads.@threads for i = 1:length(X)
-            res[i] = cond_sample(sampler, X[i])
+            res[i] = conditional_sample(sampler, X[i])
         end
         return res
     end
 end
-function cond_sample(sampler::AbstractCondSampler{d,dC,T,R}, x::PSDdata{T}) where {d,dC,T<:Number,R}
+function conditional_sample(sampler::AbstractCondSampler{d,dC,T,R}, x::PSDdata{T}) where {d,dC,T<:Number,R}
     dx = _d_marg(sampler)
-    return cond_pushforward(sampler, sample_reference(sampler)[dx+1:d], x)
+    return conditional_pushforward(sampler, sample_reference(sampler)[dx+1:d], x)
 end
 
 ## methods for Reference distribution
@@ -234,8 +308,8 @@ end
 @inline reference_pdf(sampler::AbstractCondSampler{d,<:Any,T,R}, x) where {d,T,R<:ReferenceMap} = _ref_Jacobian(sampler, x)
 @inline reference_pdf(_::AbstractCondSampler{d,<:Any,T,Nothing}, x) where {d,T} = all(1.0 .> x .> 0) ? 1.0 : 0.0
 
-@inline marg_reference_pdf(sampler::AbstractCondSampler{d,<:Any,T,R}, x) where {d,T,R<:ReferenceMap} = marg_Jacobian(sampler.R1_map, x)
-@inline marg_reference_pdf(_::AbstractCondSampler{d,<:Any,T,Nothing}, x) where {d,T} = all(1.0 .> x .> 0) ? 1.0 : 0.0
+@inline marginal_reference_pdf(sampler::AbstractCondSampler{d,<:Any,T,R}, x) where {d,T,R<:ReferenceMap} = marginal_Jacobian(sampler.R1_map, x)
+@inline marginal_reference_pdf(_::AbstractCondSampler{d,<:Any,T,Nothing}, x) where {d,T} = all(1.0 .> x .> 0) ? 1.0 : 0.0
 
 include("mappings/ProjectionMapping.jl")
 include("mappings/MarginalMapping.jl")
diff --git a/src/Samplers/samplers/PSDModelSampler.jl b/src/Samplers/samplers/PSDModelSampler.jl
index d2213ad..cfd0b8e 100644
--- a/src/Samplers/samplers/PSDModelSampler.jl
+++ b/src/Samplers/samplers/PSDModelSampler.jl
@@ -115,6 +115,20 @@ function Distributions.pdf(
     return sar.model(x)::T
 end
 
+function Distributions.logpdf(
+        sar::PSDModelSampler{d, <:Any, T},
+        x::PSDdata{T}
+    ) where {d, T<:Number}
+    return log(sar.model(x))::T
+end
+
+
+function marginal_pdf(sampler::PSDModelSampler{d, dC, T, S}, x::PSDdata{T}) where {d, dC, T<:Number, S}
+    return sampler.margins[d-dC](x)
+end
+
+marginal_logpdf(sampler::PSDModelSampler{d, dC, T, S}, x::PSDdata{T}) where {d, dC, T<:Number, S} = log(marginal_pdf(sampler, x))
+
 @inline pushforward(sampler::PSDModelSampler{d, <:Any, T, S}, 
                     u::PSDdata{T}) where {d, T<:Number, S} = return _pushforward_first_n(sampler, u, d)
 
@@ -131,26 +145,22 @@ end
 
 ## Methods for satisfying ConditionalSampler interface
 
-function marg_pdf(sampler::PSDModelSampler{d, dC, T, S}, x::PSDdata{T}) where {d, dC, T<:Number, S}
-    return sampler.margins[d-dC](x)
+function marginal_Jacobian(sampler::PSDModelSampler{d, dC, T, S}, x::PSDdata{T}) where {d, dC, T<:Number, S}
+    return 1/marginal_inverse_Jacobian(sampler, marginal_pushforward(sampler, x))
 end
 
-function marg_Jacobian(sampler::PSDModelSampler{d, dC, T, S}, x::PSDdata{T}) where {d, dC, T<:Number, S}
-    return 1/marg_inverse_Jacobian(sampler, marg_pushforward(sampler, x))
+function marginal_inverse_Jacobian(sampler::PSDModelSampler{d, dC, T, S}, x::PSDdata{T}) where {d, dC, T<:Number, S}
+    return marginal_pdf(sampler, x)
 end
 
-function marg_inverse_Jacobian(sampler::PSDModelSampler{d, dC, T, S}, x::PSDdata{T}) where {d, dC, T<:Number, S}
-    return marg_pdf(sampler, x)
-end
-
-@inline marg_pushforward(sampler::PSDModelSampler{d, dC, T, S}, 
+@inline marginal_pushforward(sampler::PSDModelSampler{d, dC, T, S}, 
                          u::PSDdata{T}) where {d, T<:Number, S, dC} = 
                             return _pushforward_first_n(sampler, u, d-dC)
-@inline marg_pullback(sampler::PSDModelSampler{d, dC, T, S}, 
+@inline marginal_pullback(sampler::PSDModelSampler{d, dC, T, S}, 
                       x::PSDdata{T}) where {d, T<:Number, S, dC} = 
                             return _pullback_first_n(sampler, x, d-dC)
 
-function cond_pushforward(sampler::PSDModelSampler{d, dC, T, S}, u::PSDdata{T}, x::PSDdata{T}) where {d, T<:Number, S, dC}
+function conditional_pushforward(sampler::PSDModelSampler{d, dC, T, S}, u::PSDdata{T}, x::PSDdata{T}) where {d, T<:Number, S, dC}
     dx = d-dC
     # @assert length(u) == dC
     # @assert length(x) == dx
@@ -161,7 +171,7 @@ function cond_pushforward(sampler::PSDModelSampler{d, dC, T, S}, u::PSDdata{T},
                 sampler.margins[k+dx-1]([x; y[1:k-1]])) - u[k]
     end
     if S<:OMF
-        for k=1:d
+        for k=1:dC
             y[k] = find_zero(f(k), zero(T))
         end
     else
@@ -170,11 +180,11 @@ function cond_pushforward(sampler::PSDModelSampler{d, dC, T, S}, u::PSDdata{T},
             y[k] = find_zero(f(k), (left, right))
         end
     end
-    return invpermute!(y, sampler.variable_ordering[dx+1:end].-dx) # might be wrong, subtract dx from variable_ordering
+    return invpermute!(y, sampler.variable_ordering[dx+1:end].-dx)
 end
 
 
-function cond_pullback(sampler::PSDModelSampler{d, dC, T, S}, 
+function conditional_pullback(sampler::PSDModelSampler{d, dC, T, S}, 
                         y::PSDdata{T},
                         x::PSDdata{T}) where {d, T<:Number, S, dC}
     dx = d-dC
diff --git a/src/Samplers/samplers/Sampler.jl b/src/Samplers/samplers/Sampler.jl
index ec0c839..9e26e42 100644
--- a/src/Samplers/samplers/Sampler.jl
+++ b/src/Samplers/samplers/Sampler.jl
@@ -55,31 +55,6 @@ function Base.show(io::IO, sra::CondSampler{d, <:Any, T}) where {d, T<:Number}
     println(io, "  reference map: $(sra.R1_map), $(sra.R2_map)")
 end
 
-## Overwrite pdf function from Distributions
-function Distributions.pdf(
-        sar::CondSampler{d, <:Any, T}, 
-        x::PSDdata{T}
-    ) where {d, T<:Number}
-    pdf_func = pushforward(sar, x->reference_pdf(sar, x))
-    return pdf_func(x)
-end
-
-function pullback(
-        sra::CondSampler{d, <:Any, T},
-        pdf_tar::Function; # defined on [a, b]^d
-        layers=nothing
-    ) where {d, T<:Number}
-    _layers = layers === nothing ? (1:length(sra.samplers)) : layers
-    # apply map R (domain transformation)
-    pdf_func = pushforward(sra.R2_map, pdf_tar)
-    for sampler in sra.samplers[_layers]
-        # apply map T_i
-        pdf_func = pullback(sampler, pdf_func)
-    end
-    return pullback(sra.R1_map, pdf_func)
-end
-
-
 function _broadcasted_pullback_pdf_function(
         sra::CondSampler{d, <:Any, T},
         broad_pdf_tar::Function; # defined on [a, b]^d
@@ -96,54 +71,102 @@ function _broadcasted_pullback_pdf_function(
     return pdf_func
 end
 
-function pushforward(
-        sra::CondSampler{d, <:Any, T},
-        pdf_ref::Function;
-        layers=nothing
-    ) where {d, T<:Number}
-    # from last to first
-    _layers = layers === nothing ? (1:length(sra.samplers)) : layers
-    pdf_func = pushforward(sra.R2_map, pdf_ref)
-    for sampler in reverse(sra.samplers[_layers])
-        # apply map T_i
-        pdf_func = pushforward(sampler, pdf_func)
+#### Macros for pushforward and pullback functions
+
+macro pushforward_function(prefix)
+    quote
+        function $(esc(Symbol(prefix, "pushforward")))(
+            sra::CondSampler{d, <:Any, T},
+            pdf_ref::Function;
+            layers=nothing
+        ) where {d, T<:Number}
+            _layers = layers === nothing ? (1:length(sra.samplers)) : layers
+            pdf_func = $(esc(Symbol(prefix, "pushforward")))(sra.R2_map, pdf_ref)
+            for j=reverse(_layers) # reverse order
+                pdf_func = $(esc(Symbol(prefix, "pushforward")))(sra.samplers[j], pdf_func)
+            end
+            return $(esc(Symbol(prefix, "pullback")))(sra.R1_map, pdf_func)
+        end
     end
-    return pullback(sra.R1_map, pdf_func)
 end
 
+macro pullback_function(prefix)
+    quote
+        function $(esc(Symbol(prefix, "pullback")))(
+            sra::CondSampler{d, <:Any, T},
+            π::Function;
+            layers=nothing
+        ) where {d, T<:Number}
+            _layers = layers === nothing ? (1:length(sra.samplers)) : layers
+            pdf_func = $(esc(Symbol(prefix, "pushforward")))(sra.R1_map, π)
+            for j=_layers
+                pdf_func = $(esc(Symbol(prefix, "pullback")))(sra.samplers[j], pdf_func)
+            end
+            return $(esc(Symbol(prefix, "pullback")))(sra.R2_map, pdf_func)
+        end
+    end
+end
 
-function marg_pushforward(
-        sra::CondSampler{d, <:Any, T},
-        pdf_ref::Function;
-        layers=nothing
-    ) where {d, T<:Number}
-    # from last to first
-    _layers = layers === nothing ? (1:length(sra.samplers)) : layers
-    pdf_func = marg_pushforward(sra.R2_map, pdf_ref)
-    for sampler in reverse(sra.samplers[_layers])
-        # apply map T_i
-        pdf_func = marg_pushforward(sampler, pdf_func)
+macro pushforward_sample(prefix)
+    quote
+        function $(esc(Symbol(prefix, "pushforward")))(
+            sra::CondSampler{d, <:Any, T}, 
+            u::PSDdata{T};
+            layers=nothing
+        ) where {d, T<:Number}
+            _layers = layers === nothing ? (1:length(sra.samplers)) : layers
+            u = $(esc(Symbol(prefix, "pushforward")))(sra.R2_map, u)
+            for j=reverse(_layers) # reverse order
+                u = $(esc(Symbol(prefix, "pushforward")))(sra.samplers[j], u)
+            end
+            u = $(esc(Symbol(prefix, "pullback")))(sra.R1_map, u)
+            return u::Vector{T}
+        end     
     end
-    pdf_func = marg_pullback(sra.R1_map, pdf_func)
-    return pdf_func
 end
 
-function marg_pullback(
-        sra::CondSampler{d, <:Any, T},
-        pdf_ref::Function;
-        layers=nothing
-    ) where {d, T<:Number}
-    # from last to first
-    _layers = layers === nothing ? (1:length(sra.samplers)) : layers
-    pdf_func = marg_pushforward(sra.R1_map, pdf_ref)
-    for sampler in sra.samplers[_layers]
-        # apply map T_i
-        pdf_func = marg_pullback(sampler, pdf_func)
+macro pullback_sample(prefix)
+    quote
+        function $(esc(Symbol(prefix, "pullback")))(
+            sra::CondSampler{d, <:Any, T}, 
+            u::PSDdata{T};
+            layers=nothing
+        ) where {d, T<:Number}
+            _layers = layers === nothing ? (1:length(sra.samplers)) : layers
+            u = $(esc(Symbol(prefix, "pushforward")))(sra.R1_map, u)
+            for j=_layers
+                u = $(esc(Symbol(prefix, "pullback")))(sra.samplers[j], u)
+            end
+            u = $(esc(Symbol(prefix, "pullback")))(sra.R2_map, u)
+            return u::Vector{T}
+        end     
     end
-    pdf_func = marg_pullback(sra.R2_map, pdf_func)
-    return pdf_func
 end
 
+###########################################
+# Function pushforward and pullback
+
+@pushforward_function ""
+@pullback_function ""
+
+@pushforward_function "log_"
+@pullback_function "log_"
+
+@pushforward_function "marginal_"
+@pullback_function "marginal_"
+
+@pushforward_function "marginal_log_"
+@pullback_function "marginal_log_"
+
+###########################################
+# Sample pushforward and pullback
+
+@pushforward_sample ""
+@pullback_sample ""
+
+@pushforward_sample "marginal_"
+@pullback_sample "marginal_"
+
 
 function add_layer!(
         sra::CondSampler{d, dC, T},
@@ -154,81 +177,54 @@ function add_layer!(
 end
 
 
-## Interface of Sampler
+"""
+Jacobian defintions
+"""
 
-function pushforward(
-        sra::CondSampler{d, <:Any, T}, 
-        u::PSDdata{T};
-        layers=nothing
-    ) where {d, T<:Number}
-    _layers = layers === nothing ? (1:length(sra.samplers)) : layers
-    u = pushforward(sra.R2_map, u)
-    for j=reverse(_layers) # reverse order
-        u = pushforward(sra.samplers[j], u)
-    end
-    u = pullback(sra.R1_map, u)
-    return u::Vector{T}
-end
+@inline inverse_Jacobian(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number} = pushforward(sra, x->one(T))(x)
+@inline Jacobian(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number} = pullback(sra, x->one(T))(x)
 
-function pullback(
-        sra::CondSampler{d, <:Any, T}, 
-        x::PSDdata{T};
-        layers=nothing
-    ) where {d, T<:Number}
-    _layers = layers === nothing ? (1:length(sra.samplers)) : layers
-    x = pushforward(sra.R1_map, x)
-    for j=_layers
-        x = pullback(sra.samplers[j], x)
-    end
-    x = pullback(sra.R2_map, x)
-    return x::Vector{T}
-end
+@inline marginal_inverse_Jacobian(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number} = marginal_pushforward(sra, x->one(T))(x)
+@inline marginal_Jacobian(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number} = marginal_pullback(sra, x->one(T))(x)
+
+@inline log_Jacobian(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number} = log_pullback(sra, x->zero(T))(x)
+@inline inverse_log_Jacobian(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number} = log_pushforward(sra, x->zero(T))(x)
+
+@inline marginal_log_Jacobian(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number} = marginal_log_pullback(sra, x->zero(T))(x)
+@inline marginal_inverse_log_Jacobian(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number} = marginal_log_pushforward(sra, x->zero(T))(x)
 
 """
-TODO: can be done prettier
+PDF definitions
 """
-function inverse_Jacobian(sra::CondSampler{d, <:Any, T},
-            x::PSDdata{T}) where {d, T<:Number}
-    pushforward(sra, x->one(T))(x)
+## Overwrite pdf function from Distributions
+function Distributions.pdf(
+        sar::CondSampler{d, <:Any, T}, 
+        x::PSDdata{T}
+    ) where {d, T<:Number}
+    log_pdf_func = log_pushforward(sar, x->log(reference_pdf(sar, x)))
+    return exp(log_pdf_func(x))
+    # pdf_func = pushforward(sar, x->reference_pdf(sar, x))
+    # return pdf_func(x)
 end
-@inline Jacobian(sra::CondSampler{d, <:Any, T},
-            x::PSDdata{T}) where {d, T<:Number} = pullback(sra, x->one(T))(x)
 
-## Interface of ConditionalSampler
+function Distributions.logpdf(
+        sar::CondSampler{<:Any, <:Any, T}, 
+        x::PSDdata{T}
+    ) where {T<:Number}
+    log_pdf_func = log_pushforward(sar, x->log(reference_pdf(sar, x)))
+    return log_pdf_func(x)
+end
 
-function marg_pdf(sra::CondSampler{d, dC, T, R1, R2}, x::PSDdata{T}) where {d, dC, T<:Number, R1, R2}
+function marginal_pdf(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number}
     @assert length(x) == _d_marg(sra)
     # reference Jacobian flexible for dimension
-    pdf_func = marg_pushforward(sra, x->marg_reference_pdf(sra, x))
+    pdf_func = marginal_pushforward(sra, x->marginal_reference_pdf(sra, x))
     return pdf_func(x)
 end
 
-function marg_inverse_Jacobian(sra::CondSampler{d, dC, T, R1, R2}, x::PSDdata{T}) where {d, dC, T<:Number, R1, R2}
-    return marg_pushforward(sra, x->one(T))(x)
-end
-
-function marg_Jacobian(sra::CondSampler{d, dC, T, R1, R2}, x::PSDdata{T}) where {d, dC, T<:Number, R1, R2}
-    return marg_pullback(sra, x->one(T))(x)
-end
-
-function marg_pushforward(sra::CondSampler{d, <:Any, T}, u::PSDdata{T};
-                        layers=nothing) where {d, T<:Number}
-    _layers = layers === nothing ? (1:length(sra.samplers)) : layers
-    u = marg_pushforward(sra.R2_map, u)
-    for j=reverse(_layers) # reverse order
-        u = marg_pushforward(sra.samplers[j], u)
-    end
-    u = marg_pullback(sra.R1_map, u)
-    return u
-end
-
-function marg_pullback(sra::CondSampler{d, <:Any, T}, x::PSDdata{T};
-                        layers=nothing) where {d, T<:Number}
-    _layers = layers === nothing ? (1:length(sra.samplers)) : layers
-    x = marg_pushforward(sra.R1_map, x)
-    for j=_layers
-        x = marg_pullback(sra.samplers[j], x)
-    end
-    x = marg_pullback(sra.R2_map, x)
-    return x
-end
+function marginal_logpdf(sra::CondSampler{<:Any, <:Any, T}, x::PSDdata{T}) where {T<:Number}
+    @assert length(x) == _d_marg(sra)
+    # reference Jacobian flexible for dimension
+    log_pdf_func = marginal_log_pushforward(sra, x->log(marginal_reference_pdf(sra, x)))
+    return log_pdf_func(x)
+end
\ No newline at end of file
diff --git a/src/SequentialMeasureTransport.jl b/src/SequentialMeasureTransport.jl
index d44d938..221febc 100644
--- a/src/SequentialMeasureTransport.jl
+++ b/src/SequentialMeasureTransport.jl
@@ -54,8 +54,4 @@ include("Samplers/sampler.jl")
 include("statistics.jl")
 using .Statistics
 
-# for AbstractMCMC interface
-include("MCMC/mcmc.jl")
-# using .MCMC
-
 end # module PositiveSemidefiniteModels
diff --git a/src/optimization/manopt_optimizer.jl b/src/optimization/manopt_optimizer.jl
deleted file mode 100644
index 1eea391..0000000
--- a/src/optimization/manopt_optimizer.jl
+++ /dev/null
@@ -1,57 +0,0 @@
-import Manopt
-import Manifolds
-
-struct ManOptProp{T} <: OptProp{T}
-    initial::AbstractMatrix{T}
-    loss::Function
-    grad_loss!::Function
-    trace::Bool
-    function ManOptProp(
-            initial::AbstractMatrix{T}, 
-            loss::Function, 
-            grad_loss!::Function;
-            trace=false
-        ) where {T<:Number}
-        new{T}(initial, loss, grad_loss!, trace)
-    end
-end
-function ManOptProp(
-            initial::AbstractMatrix{T}, 
-            loss::Function;
-            kwargs...
-        ) where {T<:Number}
-    grad_loss! = let loss = loss 
-        (M, G, x) -> begin
-            G .= Manifolds.project(M, x, FD.gradient(loss, x))
-            return G
-        end
-    end
-    loss_manopt = let loss=loss
-        (M, x) -> loss(x)
-    end
-    ManOptProp(initial, loss_manopt, grad_loss!; kwargs...)
-end
-
-function optimize(
-        prob::ManOptProp{T}
-    ) where {T<:Number}
-    n = size(prob.initial, 1)
-    M = Manifolds.SymmetricPositiveDefinite(n)
-    x0 = Matrix(prob.initial)
-    # x0 = Manifolds.rand(M)
-    Manopt.gradient_descent!(M, prob.loss,
-            prob.grad_loss!, x0;
-            evaluation=Manopt.InplaceEvaluation(),
-            stepsize=Manopt.ArmijoLinesearch(
-                M;
-                initial_stepsize=1.0,
-                contraction_factor=0.9,
-                sufficient_decrease=0.05,
-                stop_when_stepsize_less=1e-9,
-            ),
-            debug=[:Iteration,(:Change, "|Δp|: %1.10f |"),(:GradientChange, "|Δg|: %1.10f |"),
-                    (:Cost, " F(x): %1.11f | "), "\n", :Stop],
-            stopping_criterion = Manopt.StopWhenGradientChangeLess(1e-10) | Manopt.StopAfterIteration(200),
-        )
-    return Hermitian(x0)
-end
\ No newline at end of file
diff --git a/src/optimization/optimization.jl b/src/optimization/optimization.jl
index 412182a..a186197 100644
--- a/src/optimization/optimization.jl
+++ b/src/optimization/optimization.jl
@@ -2,7 +2,6 @@ abstract type OptProp{T} end
 
 include("SDP_optimizer.jl")
 include("JuMP_optimizer.jl")
-include("manopt_optimizer.jl")
 
 const _optimize_PSD_kwargs = 
     (:convex, :trace, :maxit, :tol, 
@@ -30,9 +29,7 @@ function create_SoS_opt_problem(
             loss::Function;
             kwargs...
         ) where {T<:Number}
-    if method == :manopt
-        return create_SoS_manopt_problem(initial, loss; kwargs...)
-    elseif method == :SDP
+    if method == :SDP
         return create_SoS_SDP_problem(initial, loss; kwargs...)
     else
         throw(error("Optimization method $method not implemented."))
@@ -40,26 +37,6 @@ function create_SoS_opt_problem(
     end
 end
 
-function create_SoS_manopt_problem(
-            initial::AbstractMatrix{T}, 
-            loss::Function;
-            grad_loss=nothing,
-            convex::Bool = true,
-            trace::Bool=false,
-            maxit::Int=5000,
-            tol::Real=1e-6,
-            optimizer=nothing,
-            vectorize_matrix::Bool=true,
-            normalization_constraint::Bool=false,
-            fixed_variables=nothing,
-        ) where {T<:Number}
-    if grad_loss === nothing
-        return ManOptProp(initial, loss)
-    else
-        return ManOptProp(initial, loss, grad_loss)
-    end
-end
-
 function create_SoS_SDP_problem(
             initial::AbstractMatrix{T}, 
             loss::Function;
diff --git a/test/Sampler/Sampler_test.jl b/test/Sampler/Sampler_test.jl
index 407e264..4144138 100644
--- a/test/Sampler/Sampler_test.jl
+++ b/test/Sampler/Sampler_test.jl
@@ -54,7 +54,7 @@ end
     end
 
     @testset "indefinite domain" begin
-        f(x) = pdf(MvNormal([0.0, 0.0], 4.0*diagm(ones(2))), x)
+        f(x) = pdf(MvNormal([0.0, 0.0], 4.0*I), x)
         model = PSDModel(Legendre(0.0..1.0)^2, 
                     :downward_closed, 3)
         sra = SelfReinforcedSampler(
diff --git a/test/Sampler/conditional_sampler_test.jl b/test/Sampler/conditional_sampler_test.jl
index f9d860e..aba1a6a 100644
--- a/test/Sampler/conditional_sampler_test.jl
+++ b/test/Sampler/conditional_sampler_test.jl
@@ -9,7 +9,7 @@ using SequentialMeasureTransport.BridgingDensities
     fit!(model, X, Y, trace=false)
     SMT.normalize!(model)
     cond_sampler = ConditionalSampler(model, 1)
-    Y_samp = SMT.cond_sample(cond_sampler, [[rand() * 2 - 1] for _=1:1000]; threading=false)
+    Y_samp = SMT.conditional_sample(cond_sampler, [[rand() * 2 - 1] for _=1:1000]; threading=false)
     @test all([-1≤xi[1]≤1 for xi in Y_samp])
 
     ## test pdf and marginal pdf
@@ -18,7 +18,7 @@ using SequentialMeasureTransport.BridgingDensities
     f_pdf(x) = 3/8 * f(x)
     f_marg_pdf(x) = 3/4 * x.^2 .+ 1/4
     @test norm(pdf.(Ref(cond_sampler), Y) - vcat(f_pdf.(Y)...), 2)/norm(vcat(f_pdf.(Y)...), 2) < 0.01
-    @test norm(SMT.marg_pdf.(Ref(cond_sampler), Y_marg) - vcat(f_marg_pdf.(Y_marg)...), 2)/norm(vcat(f_pdf.(Y)...), 2) < 0.01
+    @test norm(SMT.marginal_pdf.(Ref(cond_sampler), Y_marg) - vcat(f_marg_pdf.(Y_marg)...), 2)/norm(vcat(f_pdf.(Y)...), 2) < 0.01
 
     ## test conditional pdf
     Y = [[rand() * 2 - 1] for _=1:100]
@@ -26,7 +26,7 @@ using SequentialMeasureTransport.BridgingDensities
     for _=1:10
         x = [rand() * 2 - 1]
         vec1 = Y .|> y -> f_cond_pdf(y, x)
-        vec2 = Y .|> y -> SMT.cond_pdf(cond_sampler, y, x)
+        vec2 = Y .|> y -> SMT.conditional_pdf(cond_sampler, y, x)
         @test norm(vec1 - vec2, 2)/norm(vec1, 2) < 0.01
     end
 end
@@ -46,7 +46,7 @@ end
     f_pdf(x) = 3/8 * f(x)
     f_marg_pdf(x) = 3/4 * x.^2 .+ 1/4
     @test norm(pdf.(Ref(sra), Y) - vcat(f_pdf.(Y)...), 2)/norm(vcat(f_pdf.(Y)...), 2) < 0.1
-    @test norm(SMT.marg_pdf.(Ref(sra), Y_marg) - vcat(f_marg_pdf.(Y_marg)...), 2)/norm(vcat(f_marg_pdf.(Y)...), 2) < 0.1
+    @test norm(SMT.marginal_pdf.(Ref(sra), Y_marg) - vcat(f_marg_pdf.(Y_marg)...), 2)/norm(vcat(f_marg_pdf.(Y)...), 2) < 0.1
 
     ## test conditional pdf
     Y = [[rand() * 2 - 1] for _=1:100]
@@ -54,15 +54,21 @@ end
     for _=1:10
         x = [rand() * 2 - 1]
         vec1 = Y .|> y -> f_cond_pdf(y, x)
-        vec2 = Y .|> y -> SMT.cond_pdf(sra, y, x)
+        vec2 = Y .|> y -> SMT.conditional_pdf(sra, y, x)
         @test norm(vec1 - vec2, 2)/norm(vec1, 2) < 0.1
     end
+
+    ## test marginal pushforward normalized
+    marg_pf = SMT.marginal_pushforward(sra, x->0.5)
+    rng = range(-1.0, 1.0, length=1000)
+    @test isapprox((2/1000)*sum(marg_pf([x]) for x in rng), 1.0, atol=1e-2)
+    @test isapprox((2/1000)*sum(SMT.marginal_pdf(sra, [x]) for x in rng), 1.0, atol=1e-2)
 end
 
 
 @testset "2 layer conditional sampler indefinite domain" begin
     @testset "Gaussian reference" begin
-        f(x) = pdf(MvNormal(zeros(2), ones(2)), x)
+        f(x) = pdf(MvNormal(zeros(2), I), x)
         model = PSDModel(Legendre(0.0..1.0)^2, :downward_closed, 3)
         
         sra = SelfReinforcedSampler(f, model, 2, 
@@ -76,52 +82,65 @@ end
         f_pdf(x) = f(x)
         f_marg_pdf(x) = pdf(Normal(0, 1), x)
         @test norm(pdf.(Ref(sra), Y) - vcat(f_pdf.(Y)...), 2)/norm(vcat(f_pdf.(Y)...), 2) < 0.1
-        @test norm(SMT.marg_pdf.(Ref(sra), Y_marg) - vcat(f_marg_pdf.(Y_marg)...), 2)/norm(vcat(f_marg_pdf.(Y)...), 2) < 0.1
+        @test norm(SMT.marginal_pdf.(Ref(sra), Y_marg) - vcat(f_marg_pdf.(Y_marg)...), 2)/norm(vcat(f_marg_pdf.(Y)...), 2) < 0.1
 
         ## test conditional pdf
         Y = [[rand() * 2 - 1] for _=1:100]
         f_cond_pdf(y, x) = f_pdf([x; y]) / f_marg_pdf(x)[1]
         for _=1:10
             x = [rand() * 2 - 1]
-            cond_pdf2 = SMT.cond_pushforward(sra, y->SMT.marg_Jacobian(SMT.GaussianReference{2, 1, Float64}(2.0), y), x)
+            cond_pdf2 = SMT.conditional_pushforward(sra, y->SMT.marginal_Jacobian(SMT.GaussianReference{2, 1, Float64}(2.0), y), x)
             vec1 = Y .|> y -> f_cond_pdf(y, x)
-            vec2 = Y .|> y -> SMT.cond_pdf(sra, y, x)
+            vec2 = Y .|> y -> SMT.conditional_pdf(sra, y, x)
             vec3 = Y .|> y -> cond_pdf2(y)
             @test norm(vec1 - vec2, 2)/norm(vec1, 2) < 0.1
             @test norm(vec1 - vec3, 2)/norm(vec1, 2) < 0.1
         end
+
+        ## test marginal pushforward normalized
+        marg_pf = SMT.marginal_pushforward(sra, x->pdf(Normal(0, 1), x)[1])
+        rng = range(-4.0, 4.0, length=1000)
+        @test isapprox((8/1000)*sum(marg_pf([x]) for x in rng), 1.0, atol=1e-2)
+        @test isapprox((8/1000)*sum(SMT.marginal_pdf(sra, [x]) for x in rng), 1.0, atol=1e-2)
     end
 
     @testset "Algebraic reference" begin
-        f(x) = pdf(MvNormal(zeros(2), 0.5*ones(2)), x)
-        model = PSDModel(Legendre(0.0..1.0)^2, :downward_closed, 5)
+        variance = 1.0
+        f(x) = pdf(MvNormal(zeros(2), variance*I), x)
+        model = PSDModel(Legendre(0.0..1.0)^2, :downward_closed, 6)
         
-        bridge = AlgebraicBridgingDensity{2}(f, [0.5, 1.0])
-        sra = SelfReinforcedSampler(bridge, model, 2, 
+        bridge = AlgebraicBridgingDensity{2}(f, [1.0, 1.0])
+        sra = SelfReinforcedSampler(bridge, model, length(bridge.β_list), 
                     :Chi2, SMT.AlgebraicReference{2, 1, Float64}(); 
-                    N_sample=1000,
+                    N_sample=2000,
                     dC=1)
 
         ## test pdf and marginal pdf
-        Y_marg = [[rand() * 2 - 1] for _=1:100]
-        Y = [rand(2) * 2 .- 1 for _=1:100]
+        Y_marg = [[rand() * 3 - 1] for _=1:100]
+        Y = [rand(2) * 3 .- 1 for _=1:100]
         f_pdf(x) = f(x)
-        f_marg_pdf(x) = pdf(Normal(0, 0.5), x)
+        f_marg_pdf(x) = pdf(Normal(0, variance), x)
         @test norm(pdf.(Ref(sra), Y) - vcat(f_pdf.(Y)...), 2)/norm(vcat(f_pdf.(Y)...), 2) < 0.1
-        @test norm(SMT.marg_pdf.(Ref(sra), Y_marg) - vcat(f_marg_pdf.(Y_marg)...), 2)/norm(vcat(f_marg_pdf.(Y)...), 2) < 0.1
+        @test norm(SMT.marginal_pdf.(Ref(sra), Y_marg) - vcat(f_marg_pdf.(Y_marg)...), 2)/norm(vcat(f_marg_pdf.(Y_marg)...), 2) < 0.1
 
         ## test conditional pdf
         Y = [[rand() * 2 - 1] for _=1:100]
         f_cond_pdf(y, x) = f_pdf([x; y]) / f_marg_pdf(x)[1]
         for _=1:10
             x = [rand() * 2 - 1]
-            cond_pdf2 = SMT.cond_pushforward(sra, y->SMT.marg_Jacobian(SMT.AlgebraicReference{2, 1, Float64}(), y), x)
+            cond_pdf2 = SMT.conditional_pushforward(sra, y->SMT.marginal_Jacobian(SMT.AlgebraicReference{2, 1, Float64}(), y), x)
             vec1 = Y .|> y -> f_cond_pdf(y, x)
-            vec2 = Y .|> y -> SMT.cond_pdf(sra, y, x)
+            vec2 = Y .|> y -> SMT.conditional_pdf(sra, y, x)
             vec3 = Y .|> y -> cond_pdf2(y)
             @test norm(vec1 - vec2, 2)/norm(vec1, 2) < 0.1
             @test norm(vec1 - vec3, 2)/norm(vec1, 2) < 0.1
         end
+
+        ## test marginal pushforward normalized
+        marg_pf = SMT.marginal_pushforward(sra, x->pdf(Normal(0, 1), x)[1])
+        rng = range(-4.0, 4.0, length=1000)
+        @test isapprox((8/1000)*sum(marg_pf([x]) for x in rng), 1.0, atol=1e-2)
+        @test isapprox((8/1000)*sum(SMT.marginal_pdf(sra, [x]) for x in rng), 1.0, atol=1e-2)
     end
 end
 
@@ -138,8 +157,8 @@ end
     samples = [x for x in eachcol(_samples)]
     model = PSDModel(Legendre(0.0..1.0)^2, :downward_closed, 3)
     bridge = DiffusionBrigdingDensity{2}(x->1.0, [1.0, 0.5, 0.25, 0.1, 0.05, 0.02, 0.01, 0.0], 2.5)
-    ref_map = SequentialMeasureTransport.ReferenceMaps.GaussianReference{2, 1, Float64}(2.5)
-    sra = SequentialMeasureTransport.SelfReinforced_ML_estimation(
+    ref_map = SMT.ReferenceMaps.GaussianReference{2, 1, Float64}(2.5)
+    sra = SMT.SelfReinforced_ML_estimation(
         samples,
         model,
         bridge,
@@ -155,6 +174,6 @@ end
         end
     end
     for x in range(-4.0, 4.0, length=100)
-        @test abs(SequentialMeasureTransport.marg_pdf(sra, [x]) - f_marg(x)) < 0.2
+        @test abs(SMT.marginal_pdf(sra, [x]) - f_marg(x)) < 0.2
     end
 end
diff --git a/test/Sampler/reference_maps_test.jl b/test/Sampler/reference_maps_test.jl
index de83567..4b099ee 100644
--- a/test/Sampler/reference_maps_test.jl
+++ b/test/Sampler/reference_maps_test.jl
@@ -27,6 +27,47 @@ import SequentialMeasureTransport as SMT
             end
         end
     end
+    @testset "log(T^♯ T_♯ f) = log(f)" begin
+        for d=1:5
+            L = rand(d)
+            R = L + 2.0 * rand(d) .+ 0.1
+            ref_map = SMT.ScalingReference{d}(L, R)
+            f_app = SMT.log_pullback(ref_map, x->log(1.0))
+            f_app_2 = SMT.log_pushforward(ref_map, f_app)
+
+            for i=1:100
+                x = rand(d)
+                @test exp(f_app_2(x)) ≈ 1.0
+            end
+        end
+    end
+    @testset "T^♯ f normalized" begin
+        for d=1:2
+            L = rand(d)
+            R = L + 2.0 * rand(d) .+ 0.1
+            V = prod(R - L)
+            ref_map = SMT.ScalingReference{d}(L, R)
+            f_app = SMT.pushforward(ref_map, x->1/V)
+            rng = 0.001:0.001:0.999
+            iter = Iterators.product(fill(rng, d)...)
+            @test isapprox(0.001^d*sum(f_app([x...]) for x in iter), 1.0, atol=1e-2)
+        end
+    end
+    @testset "exp( log(T^♯ f)) normalized" begin
+        for d=1:2
+            L = rand(d)
+            R = L + 2.0 * rand(d) .+ 0.1
+            V = prod(R - L)
+            ref_map = SMT.ScalingReference{d}(L, R)
+            f_app = SMT.log_pushforward(ref_map, x->log(1/V))
+            f_app = let f_app=f_app 
+                x->exp(f_app(x))
+            end
+            rng = 0.001:0.001:0.999
+            iter = Iterators.product(fill(rng, d)...)
+            @test isapprox(0.001^d*sum(f_app([x...]) for x in iter), 1.0, atol=1e-2)
+        end
+    end
 end
 
 @testset "Gaussian Reference" begin
@@ -55,6 +96,42 @@ end
             end
         end
     end
+    @testset "log(T^♯ T_♯ f) = lof(f)" begin
+        for d=1:5
+            for _=1:5
+                s = rand() * 10.0
+                ref_map = SMT.GaussianReference{d, Float64}(s)
+                f_app = SMT.log_pullback(ref_map, x->0.0)
+                f_app_2 = SMT.log_pushforward(ref_map, f_app)
+        
+                for i=1:100
+                    x = rand(d)
+                    @test f_app_2(x) ≈ 0.0
+                end
+            end
+        end
+    end
+    @testset "T^♯ f normalized" begin
+        for d=1:2
+            ref_map = SMT.AlgebraicReference{d, Float64}()
+            f_app = SMT.pushforward(ref_map, x->pdf(MvNormal(zeros(d), I), x))
+            rng = 0.001:0.001:0.999
+            iter = Iterators.product(fill(rng, d)...)
+            @test isapprox(0.001^d*sum(f_app([x...]) for x in iter), 1.0, atol=1e-2)
+        end
+    end
+    @testset "exp(log(T^♯ f)) normalized" begin
+        for d=1:2
+            ref_map = SMT.AlgebraicReference{d, Float64}()
+            f_app = SMT.log_pushforward(ref_map, x->logpdf(MvNormal(zeros(d), I), x))
+            f_app = let f_app=f_app 
+                x->exp(f_app(x))
+            end
+            rng = 0.001:0.001:0.999
+            iter = Iterators.product(fill(rng, d)...)
+            @test isapprox(0.001^d*sum(f_app([x...]) for x in iter), 1.0, atol=1e-2)
+        end
+    end
 end
 
 @testset "Algebraic Reference" begin
@@ -80,4 +157,37 @@ end
             end
         end
     end
+    @testset "log(T^♯ T_♯ f) = log(f)" begin
+        for d=1:5
+            ref_map = SMT.AlgebraicReference{d, Float64}()
+            f_app = SMT.log_pullback(ref_map, x->0.0)
+            f_app_2 = SMT.log_pushforward(ref_map, f_app)
+    
+            for i=1:100
+                x = rand(d)
+                @test isapprox(f_app_2(x), 0.0, atol=1e-13)
+            end
+        end
+    end
+    @testset "T^♯ f normalized" begin
+        for d=1:2
+            ref_map = SMT.AlgebraicReference{d, Float64}()
+            f_app = SMT.pushforward(ref_map, x->pdf(MvNormal(zeros(d), I), x))
+            rng = 0.001:0.001:0.999
+            iter = Iterators.product(fill(rng, d)...)
+            @test isapprox(0.001^d*sum(f_app([x...]) for x in iter), 1.0, atol=1e-2)
+        end
+    end
+    @testset "exp(log(T^♯ f)) normalized" begin
+        for d=1:2
+            ref_map = SMT.AlgebraicReference{d, Float64}()
+            f_app = SMT.log_pushforward(ref_map, x->logpdf(MvNormal(zeros(d), I), x))
+            f_app = let f_app=f_app 
+                x->exp(f_app(x))
+            end
+            rng = 0.001:0.001:0.999
+            iter = Iterators.product(fill(rng, d)...)
+            @test isapprox(0.001^d*sum(f_app([x...]) for x in iter), 1.0, atol=1e-3)
+        end
+    end
 end
\ No newline at end of file
diff --git a/test/Sampler/subset_sampler_test.jl b/test/Sampler/subset_sampler_test.jl
index b6ece10..b1bbb65 100644
--- a/test/Sampler/subset_sampler_test.jl
+++ b/test/Sampler/subset_sampler_test.jl
@@ -2,7 +2,7 @@
 
 
 @testset "Single Projection" begin
-    f(x) = pdf(MvNormal(0.5.+zeros(2), 0.2*diagm(ones(2))), x)
+    f(x) = pdf(MvNormal(0.5.+zeros(2), 0.2*I), x)
     f_marg(x) = pdf(Normal(0.5, 0.2), x)[1]
     model = PSDModel(Legendre(0.0..1.0), :downward_closed, 3)
     X = rand(1, 1000)
@@ -29,7 +29,8 @@
     #compare pdf with marginal
     rng = [[x, y] for x in rand(20), y in rand(20)]
     rng = reshape(rng, length(rng))
-    vec1 = pdf.(Ref(proj_map), rng)
+    pdf_func = SMT.pushforward(proj_map, x->1.0)
+    vec1 = pdf_func.(rng)
     @test norm(vec1 - f_marg.(rng), 2)/norm(f_marg.(rng), 2) < 0.4
 end
 
@@ -37,9 +38,9 @@ end
     variance = 0.5
     scaling_1d = cdf(Normal(0.5, variance), 1.0) - cdf(Normal(0.5, variance), 0.0)
     rng = range(0, 1, length=200)
-    scaling_2d = 1/length(rng)^2 * sum(pdf(MvNormal(0.5.+zeros(2), variance*diagm(ones(2))), [x, y]) for x in rng, y in rng)
-    f(x) = pdf(MvNormal(0.5.+zeros(3), variance*diagm(ones(3))), x)
-    f_marg(x) = (1/scaling_2d) * pdf(MvNormal(0.5.+zeros(2), variance*diagm(ones(2))), x)
+    scaling_2d = 1/length(rng)^2 * sum(pdf(MvNormal(0.5.+zeros(2), variance*I), [x, y]) for x in rng, y in rng)
+    f(x) = pdf(MvNormal(0.5.+zeros(3), variance*I), x)
+    f_marg(x) = (1/scaling_2d) * pdf(MvNormal(0.5.+zeros(2), variance*I), x)
     f_margin_single(x) = (1/scaling_1d) * pdf(Normal(0.5, variance), x[1])
     model = PSDModel(Legendre(0.0..1.0)^2, :downward_closed, 5)
     X = rand(2, 1000)
@@ -68,11 +69,11 @@ end
     ## for the marginal
     for i=1:100
         x = rand(2)
-        @test isapprox(x, SMT.marg_pullback(proj_map, SMT.marg_pushforward(proj_map, x)), atol=1e-9)
+        @test isapprox(x, SMT.marginal_pullback(proj_map, SMT.marginal_pushforward(proj_map, x)), atol=1e-9)
     end
 
     ## for the marginal
-    pb_pf_func = SMT.marg_pullback(proj_map, SMT.marg_pushforward(proj_map, x->1.0))
+    pb_pf_func = SMT.marginal_pullback(proj_map, SMT.marginal_pushforward(proj_map, x->1.0))
     for i=1:100
         x = rand(2)
         @test isapprox(pb_pf_func(x), 1.0, atol=1e-9)
@@ -83,13 +84,15 @@ end
     #compare pdf with marginal
     rng = [[x, y, z] for x in rand(20), y in rand(20), z in rand(20)]
     rng = reshape(rng, length(rng))
-    vec1 = pdf.(Ref(proj_map), rng)
+    pdf_func = SMT.pushforward(proj_map, x->1.0)
+    vec1 = pdf_func.(rng)
     @test norm(vec1 - f_marg2.(rng), 2)/norm(f_marg2.(rng), 2) < 0.2
 
     # compare marginal with second marginal
     rng = [[x, y] for x in rand(20), y in rand(20)]
     rng = reshape(rng, length(rng))
-    vec1 = SMT.marg_pdf.(Ref(proj_map), rng)
+    marginal_pdf_func = SMT.marginal_pushforward(proj_map, x->1.0)
+    vec1 = marginal_pdf_func.(rng)
     @test norm(vec1 - f_margin_single.(rng), 2)/norm(f_margin_single.(rng), 2) < 0.2
 end
 
@@ -106,7 +109,7 @@ end
     X2 = rand(distr2, N2)
     X = hcat(X1, X2)
     T = Float64
-    bridge = DiffusionBrigdingDensity{1}(f, T[1.5, 1.3, 1.0, 0.8, 0.75, 0.5, 0.3, 0.25, 0.18, 0.13, 0.1, 0.07, 0.02, 0.01, 0.005, 0.0], T(2.0))
+    bridge = DiffusionBrigdingDensity{2}(f, T[1.5, 1.3, 1.0, 0.8, 0.75, 0.5, 0.3, 0.25, 0.18, 0.13, 0.1, 0.07, 0.02, 0.01, 0.005, 0.0], T(2.0))
     ref_map = SMT.ReferenceMaps.GaussianReference{2, T}(T(2.0))
     to_subspace_ref_map = SMT.ReferenceMaps.GaussianReference{2, T}(T(2.0))
     subspace_ref_map = SMT.ReferenceMaps.GaussianReference{1, T}(T(2.0))
@@ -160,10 +163,10 @@ end
     X2 = rand(distr2, N2)
     X = hcat(X1, X2)
     T = Float64
-    bridge = DiffusionBrigdingDensity{2}(f, T[1.8, 1.5, 1.3, 1.2, 1.1, 1.0, 0.8, 0.75, 
+    bridge = DiffusionBrigdingDensity{3}(f, T[1.8, 1.5, 1.3, 1.1, 1.0, 0.8, 0.75, 
                                             0.6, 0.5, 0.25, 0.18, 0.13,
                                             0.1, 0.07, 0.05, 0.02, 0.001, 0.007, 
-                                            0.005, 0.003, 0.001, 0.0], T(1.0))
+                                            0.005, 0.001, 0.0], T(1.0))
     # ref_map = SequentialMeasureTransport.ReferenceMaps.GaussianReference{3, T}(T(2.5))
     # to_subspace_ref_map = SequentialMeasureTransport.ReferenceMaps.GaussianReference{3, T}(T(3.0))
     # subspace_ref_map = SequentialMeasureTransport.ReferenceMaps.GaussianReference{2, T}(T(3.0))
@@ -198,8 +201,8 @@ end
         ]
     rng = reshape(rng, length(rng))
     rng_marg = reshape(rng_marg, length(rng_marg))
-    @test norm(pdf.(Ref(sra_sub), rng) - f.(rng), 2)/norm(f.(rng), 2) < 0.4
-    @test norm(SMT.marg_pdf.(Ref(sra_sub), rng_marg) - f_marg.(rng_marg), 2)/norm(f_marg.(rng_marg), 2) < 0.4
+    @test norm(pdf.(Ref(sra_sub), rng) - f.(rng), 2)/norm(f.(rng), 2) < 0.45
+    @test norm(SMT.marginal_pdf.(Ref(sra_sub), rng_marg) - f_marg.(rng_marg), 2)/norm(f_marg.(rng_marg), 2) < 0.45
     
 
     ## check that conditional is normalized
@@ -211,37 +214,44 @@ end
         else
             rand(marg_distr2)
         end
-        cond_pdf_func = SMT.cond_pushforward(sra_sub, y->SMT.Jacobian(SMT.AlgebraicReference{1, 0, Float64}(), y), x)
+        cond_pdf_func = SMT.conditional_pushforward(sra_sub, y->SMT.Jacobian(SMT.AlgebraicReference{1, 0, Float64}(), y), x)
         int_cond = Δrng*sum(cond_pdf_func([y]) for y in rng)
-        int_cond2 = Δrng*sum(SMT.cond_pdf(sra_sub, [y], x) for y in rng)
+        int_cond2 = Δrng*sum(SMT.conditional_pdf(sra_sub, [y], x) for y in rng)
         @test isapprox(int_cond, 1.0, atol=0.05)
         @test isapprox(int_cond2, 1.0, atol=0.05)
     end
 
+    rng = [[x...] for x in Iterators.product(
+                range(-5, 5, length=N),
+                range(-5, 5, length=N),
+                range(-5, 5, length=N)
+            )
+        ]
     # comparison of conditional difficult, use conditional negative log likelihood
-    model_c_vec = rng .|> (x)->SMT.cond_pdf(sra_sub, x[3:3], x[1:2])
-    c_vec = rng .|> (x)->f_cond(x[1:2], x[3:3])
-    @test norm(model_c_vec - c_vec, 2)/norm(c_vec, 2) < 0.1
+    # model_c_vec = rng .|> (x)->SMT.conditional_pdf(sra_sub, x[3:3], x[1:2])
+    # c_vec = rng .|> (x)->f_cond(x[1:2], x[3:3])
+    # @test norm(model_c_vec - c_vec, 2)/norm(c_vec, 2) < 0.1
+
     X1 = rand(distr1, N1)
     X2 = rand(distr2, N2)
     X = hcat(X1, X2)
     cond_KL = (1/N) * sum(log.(f_cond(x[1:2], x[3:3])) for x in eachcol(X))
 
-    cond_neg_log_likelihood = (1/N) * sum(log.(SMT.cond_pdf(sra_sub, x[3:3], x[1:2])) for x in eachcol(X))
+    cond_neg_log_likelihood = (1/N) * sum(SMT.conditional_logpdf(sra_sub, x[3:3], x[1:2]) for x in eachcol(X))
 
     KL_cond = cond_KL - cond_neg_log_likelihood
 
-    @assert KL_cond < 8.0
+    @assert KL_cond < 10.0
 
 
     ## test that pushforward is pullback
     for i=1:100
         x = rand(2)
-        @test isapprox(x, SMT.marg_pullback(sra_sub, SMT.marg_pushforward(sra_sub, x)), atol=1e-9)
+        @test isapprox(x, SMT.marginal_pullback(sra_sub, SMT.marginal_pushforward(sra_sub, x)), atol=1e-9)
     end
 
     ## test that marginal pushforward is pullback
-    pb_pf_func = SMT.marg_pullback(sra_sub, SMT.marg_pushforward(sra_sub, x->1.0))
+    pb_pf_func = SMT.marginal_pullback(sra_sub, SMT.marginal_pushforward(sra_sub, x->1.0))
     for i=1:100
         x = rand(2)
         @test isapprox(pb_pf_func(x), 1.0, atol=1e-9)
@@ -252,11 +262,11 @@ end
         x = rand(2)
         for i=1:100
             y = rand(1)
-            y_pf = SMT.cond_pushforward(sra_sub, y, x)
-            @test isapprox(y, SMT.cond_pullback(sra_sub, y_pf, x), atol=1e-9)
+            y_pf = SMT.conditional_pushforward(sra_sub, y, x)
+            @test isapprox(y, SMT.conditional_pullback(sra_sub, y_pf, x), atol=1e-9)
         end
 
-        pb_pf_func = SMT.cond_pullback(sra_sub, SMT.cond_pushforward(sra_sub, x->1.0, x), x)
+        pb_pf_func = SMT.conditional_pullback(sra_sub, SMT.conditional_pushforward(sra_sub, x->1.0, x), x)
         for i=1:100
             y = rand()
             @test isapprox(pb_pf_func([y]), 1.0, atol=1e-9)
@@ -271,15 +281,15 @@ end
         model = PSDModel(Legendre(-5.0..5.0), :downward_closed, 4)
         X = rand(1, 1000)
         Y = map(x->f(x)[1], eachcol(X))
-        SequentialMeasureTransport.Chi2_fit!(model, eachcol(X), Y, trace=false)
+        SMT.Chi2_fit!(model, eachcol(X), Y, trace=false)
         normalize!(model)
         smp = Sampler(model)
 
-        smp_marg = SequentialMeasureTransport.MarginalMapping{2, 0}(smp, [1])
+        smp_marg = SMT.MarginalMapping{2, 0}(smp, [1])
 
         # compare pdfs
-        f_app = SequentialMeasureTransport.pushforward(smp, x->1.0)
-        f_app_marg = SequentialMeasureTransport.pushforward(smp_marg, x->1.0)
+        f_app = SMT.pushforward(smp, x->1.0)
+        f_app_marg = SMT.pushforward(smp_marg, x->1.0)
 
         rng = range(-5, 5, length=100)
         for x in rng