tpapp / transformvariables.jl Goto Github PK

View Code? Open in Web Editor NEW

65.0 3.0 14.0 1.21 MB

Transformations to contrained variables from ℝⁿ.

License: Other

Julia 100.00%

julia julia-package bayesian-inference bayesian-statistics

transformvariables.jl's Introduction

TransformVariables.jl

Successor of ContinuousTransformations.jl.

Features

Simple interface.
Fast implementation, unrolling when it makes sense.
Targeted to applications in statistics, mainly ML, MAP, and Bayesian inference.
Take advantage of Julia v0.7 features, especially named tuples and compiler optimizations.

Caveat

Work in progress. API will change rapidly, without warning.
Expect speed regressions until API is finalized.

transformvariables.jl's People

Contributors

Stargazers

Watchers

Forkers

simonbyrne andreasnoack ericphanson tkf sethaxen eloceanografo stjordanis cscherrer chriselrod scheidan ptiede junyuan-chen

transformvariables.jl's Issues

More tests and correctness checks

identity array transformations, especially with generalized indexing (should just copy, check that it is aligned correctly)
use the checkbounds framework for index checking with IndexInto

Add ASCII aliases for as* functions

The asℝ, asℝ₊, asℝ₊, as𝕀 are difficult to type because (as far as I know) the "Mathematical Double-Struck Capital R" character and others don't work with the backslash completion in the REPL or in IJulia. It would be nice if these were also aliased as asR, asRp, asRm, asI or similar.

Proposed transformation: Loading matrix for factor analysis

Would there be any interest in a "special array" transformation to create loading matrices for factor analysis? i.e., when n observed variables are represented by m <= n latent factors, the loading matrix L has size (n, m) and links the latent variables x to observed variables y as y ~ MvNormal(L * x, sigma). When fitting it by MLE/MCMC, all entries above the diagonal are set to zero to ensure the columns are linearly independent. The transform is thus from a vector with (n - m + 2) * m elements to an n * m matrix.

I've got a working example hacked together that I could turn into a pull request if this would be useful for others...

tag a release

TODO:

documentation (generate, polish)
examples (update DynamicHMCExamples.jl)
mention that the internals are a moving target in the manual

Feature idea: function to check if values lay in domains of transformation

As an additional feature it would be nice if we could use a transform to check if the provided inputs are legal.

For example

in_domain(t, y)

would return true if all values of y lay in the domains specified by t.

Such a function would make it convenient to check user provided arguments (i.e. model parameters).

Maybe this could be introduced when dealing with #24 or #64.

StaticArray results and input

Allow versions of transformations that return StaticArrays.

Roadmap:

include size information in types
SMatrix, SVector, and possibly MMatrix and MVector versions of UnitVector and CorrCholeskyFactor, maybe specified with the {} in the constructor, also include size there?
Array aggregator should also allow this, maybe with as(SVector{size}, ...) syntax? Possibly a unified syntax for both static and plain vanilla arrays.
Tuple and NamedTuple aggregators can probably be untouched, caching dimensions information could just lead to (minor) speedups.

Bikeshedding syntax/API is most welcome.

Compute logistic and logistic_logjac together

I noticed this bit of code:

logistic(x::Real) = inv(one(x) + exp(-x))

function logistic_logjac(x::Real)
    mx = -abs(x)
    mx - 2*log1p(exp(mx))
end

There's a neat trick that can speed this up. If we write σ=logistic (for brevity) then
σ'=σ*(1-σ)

So something like

function logistic_with_logjac(x::Real)
    y = logistic(x)
    (y, log(y*(1-y)))
end

Transformations involving `StaticArray` is not type stable

One important benefit of using StaticArray is to avoid the memory allocation required for an Array. For now, there is type instability somewhere that results in allocations for type inference, which largely undermines the advantages of using StaticArray over Array.

using TransformVariables, StaticArrays, BenchmarkTools

tr = as(SVector{3, Float64})
a = rand(3)

julia> @btime transform($tr, $a)
  39.816 ns (2 allocations: 32 bytes)
3-element SVector{3, Float64} with indices SOneTo(3):
 0.8664582577624973
 0.27172807097850415
 0.5758008759582858

From @code_warntype, I see that there is type stability issue with transform_with.

inverse(t, y) error message unhelpful for non-matching (named) tuples

Eg as in

inverse(as((a = .., b = ...)), (c = ..., d = ...))

which is currently a MethodError. Similarly for tuples of wrong type/length.

An @argcheck may be better than overspecializing the method signature, and serve a better error message.

Problems in calculating the log absolute Jacobian determinant of the CorrCholeskyFactor transformation

The way TransformVariables.jl constructs the Cholesky factor of a correlation seems similar to what's given here https://mc-stan.org/docs/reference-manual/correlation-matrix-transform.html, with the only difference being $z = tanh(y/2)$ instead of $z = tanh(y)$. The log absolute Jacobian determinant should be $$\sum_{i=1}^{K-1}\sum_{j=i+1}^{K}\frac{K-i-1}{2}\ln[1-(z_{i,j})^2]+\sum_{i=1}^{K-1}\sum_{j=i+1}^{K}\ln\frac{\partial z_{i,j}}{\partial y_{i,j}}$$ If $K=3$, it's $$\frac{1}{2}(\ln[1-(z_{1,2})^2]+\ln[1-(z_{1,3})^2])+\sum_{i=1}^{2}\sum_{j=i+1}^{3}\ln\frac{\partial z_{i,j}}{\partial y_{i,j}}$$

But the calculation in the library misses $\frac{1}{2}\ln[1-(z_{1,2})^2]$ when K=3, and this problem generalizes to larger Ks. The code below illustrates this problem for $K=3$.

t = as(Array,CorrCholeskyFactor(3))
x =  [-0.12833806457027003, 0.24721188811516634, -0.24578748441593382]
logjac_pack = TransformVariables.transform_and_logjac(t,x)[2]

#= 
Calculate the change-of-variable adjustment for LKJ transformation
=#
function den_adj(V::Vector{T}) where T
    N = length(V);
    K_X = Int(sqrt(2*N+1/4)+1/2);
    J_X = 0
    count = 1
    @inbounds for j in 2:K_X
        for i in 1:j-1
            yij =V[count]
            J_X += (K_X-i-1)/2*(log(4)+yij-2*log(exp(yij)+1))+log(2)+yij-2*log(1+exp(yij))
            count += 1
        end
    end
    return J_X
end

logjac_man = den_adj(x)

isequal(logjac_pack,den_adj(x)-1/2* (log(4)+x[1]-2*log(exp(x[1])+1)))  # this returns true

Request for documentation of the UnitVector transformation

Dear Tamas Papp,

Thanks for your efforts on TransformVariables.jl. Do you happen do know a source (paper/textbook/blog) which explains the transformation which underlies the UnitVector type?

best
Jon Eriksen

long as(::NamedTuple) fails with Enzyme

This is EnzymeAD/Enzyme.jl#1235 in the wild.

using TransformVariables, Enzyme, StaticArrays

K = 7
N_EE = N_UU = N_UE = N_EU = 20
trans = as((# common
            ω_intercept = as(SVector{K}), ω_std = as(SVector{K}, asℝ₊),
            ω_corr_factor = as(SVector{K}), ζ_std = asℝ₊, ε_std = asℝ₊,
            κ = as(Real, 0.5, 1.5),
            B1 = as(SVector{3}), B2 = as(SVector{3}), BC = as(SVector{3}),
            # EE
            α̂1_EE = as(view, N_EE), α̂2_EE = as(view, N_EE),
            β̂1_EE = as(view, N_EE), β̂2_EE = as(view, N_EE),
            M̂_EE = as(view, N_EE),
            # EU
            α̂1_EU = as(view, N_EU), α̂2_EU = as(view, N_EU),
            β̂1_EU = as(view, N_EU), β̂2_EU = as(view, N_EU),
            M̂_EU = as(view, N_EU), ŵ2_EU = as(view, N_EU),
            # UE
            α̂1_UE = as(view, N_UE), α̂2_UE = as(view, N_UE),
            β̂1_UE = as(view, N_UE), β̂2_UE = as(view, N_UE),
            M̂_UE = as(view, N_UE), ŵ1_UE = as(view, N_UE),
            # UU
            α̂1_UU = as(view, N_UU), α̂2_UU = as(view, N_UU),
            β̂1_UU = as(view, N_UU), β̂2_UU = as(view, N_UU),
            M̂_UU = as(view, N_UU),
            ŵ1_UU = as(view, N_UU), ŵ2_UU = as(view, N_UU),
        ))

_s(x::Real) = x                 # simple recursive sum, for testing
_s(x::AbstractArray) = sum(_s, x)
_s(x::NamedTuple) = sum(_s, values(x))

g(t, x) = _s(transform(t, x))
x = zeros(dimension(trans))
g(trans, x)                     # sanity check that primal call works
∂ℓ_∂x = zero(x)
_, y = Enzyme.autodiff(Enzyme.ReverseWithPrimal, g,
                       Enzyme.Active, Enzyme.Const(trans), Enzyme.Duplicated(x, ∂ℓ_∂x))

fails with

ERROR: AssertionError: Found unhandled active variable in tuple splat, jl_apply_iterate @NamedTuple{ω_intercept::TransformVariables.StaticArrayTransformation{7, Tuple{7}, TransformVariables.Identity}, ω_std::TransformVariables.StaticArrayTransformation{7, Tuple{7}, TransformVariables.ShiftedExp{true, Int64}}, ω_corr_factor::TransformVariables.StaticArrayTransformation{7, Tuple{7}, TransformVariables.Identity}, ζ_std::TransformVariables.ShiftedExp{true, Int64}, ε_std::TransformVariables.ShiftedExp{true, Int64}, κ::TransformVariables.ScaledShiftedLogistic{Float64}, B1::TransformVariables.StaticArrayTransformation{3, Tuple{3}, TransformVariables.Identity}, B2::TransformVariables.StaticArrayTransformation{3, Tuple{3}, TransformVariables.Identity}, BC::TransformVariables.StaticArrayTransformation{3, Tuple{3}, TransformVariables.Identity}, α̂1_EE::TransformVariables.ViewTransformation{1}, α̂2_EE::TransformVariables.ViewTransformation{1}, β̂1_EE::TransformVariables.ViewTransformation{1}, β̂2_EE::TransformVariables.ViewTransformation{1}, M̂_EE::TransformVariables.ViewTransformation{1}, α̂1_EU::TransformVariables.ViewTransformation{1}, α̂2_EU::TransformVariables.ViewTransformation{1}, β̂1_EU::TransformVariables.ViewTransformation{1}, β̂2_EU::TransformVariables.ViewTransformation{1}, M̂_EU::TransformVariables.ViewTransformation{1}, ŵ2_EU::TransformVariables.ViewTransformation{1}, α̂1_UE::TransformVariables.ViewTransformation{1}, α̂2_UE::TransformVariables.ViewTransformation{1}, β̂1_UE::TransformVariables.ViewTransformation{1}, β̂2_UE::TransformVariables.ViewTransformation{1}, M̂_UE::TransformVariables.ViewTransformation{1}, ŵ1_UE::TransformVariables.ViewTransformation{1}, α̂1_UU::TransformVariables.ViewTransformation{1}, α̂2_UU::TransformVariables.ViewTransformation{1}, β̂1_UU::TransformVariables.ViewTransformation{1}, β̂2_UU::TransformVariables.ViewTransformation{1}, M̂_UU::TransformVariables.ViewTransformation{1}, ŵ1_UU::TransformVariables.ViewTransformation{1}, ŵ2_UU::TransformVariables.ViewTransformation{1}}
Stacktrace:
  [1] error_if_active_iter(arg::Base.RefValue{@NamedTuple{…}})
    @ Enzyme.Compiler ~/.julia/packages/Enzyme/Dd2LU/src/rules/jitrules.jl:775
  [2] Tuple
    @ ./namedtuple.jl:200 [inlined]
  [3] values
    @ ./namedtuple.jl:379 [inlined]
  [4] transform_with
    @ ~/code/julia/TransformVariables/src/aggregation.jl:388 [inlined]
  [5] transform
    @ ~/code/julia/TransformVariables/src/generic.jl:268
  [6] g
    @ ./REPL[31]:1 [inlined]
  [7] g
    @ ./REPL[31]:0 [inlined]
  [8] augmented_julia_g_3346_inner_1wrap
    @ ./REPL[31]:0
  [9] macro expansion
    @ Enzyme.Compiler ~/.julia/packages/Enzyme/Dd2LU/src/compiler.jl:5299 [inlined]
 [10] enzyme_call(::Val{…}, ::Ptr{…}, ::Type{…}, ::Type{…}, ::Val{…}, ::Type{…}, ::Type{…}, ::Const{…}, ::Type{…}, ::Const{…}, ::Duplicated{…})
    @ Enzyme.Compiler ~/.julia/packages/Enzyme/Dd2LU/src/compiler.jl:4977
 [11] (::Enzyme.Compiler.AugmentedForwardThunk{…})(::Const{…}, ::Const{…}, ::Vararg{…})
    @ Enzyme.Compiler ~/.julia/packages/Enzyme/Dd2LU/src/compiler.jl:4930
 [12] autodiff(::ReverseMode{…}, ::Const{…}, ::Type{…}, ::Const{…}, ::Vararg{…})
    @ Enzyme ~/.julia/packages/Enzyme/Dd2LU/src/Enzyme.jl:198
 [13] autodiff(::ReverseMode{true, FFIABI}, ::typeof(g), ::Type, ::Const{TransformVariables.TransformTuple{…}}, ::Vararg{Any})    @ Enzyme ~/.julia/packages/Enzyme/Dd2LU/src/Enzyme.jl:224
 [14] top-level scope
    @ REPL[35]:1

Feature Idea: scaling factor for ShiftedExp

It would be useful, if we can optionally define a scale to ShiftedExp. For example, if we have two parameter a and b that are both positive, but have very different scales, we could write:

t = as((a = as(Real, 0, ∞, scale=1),
        b = as(Real, 0, ∞, scale=1000)))

which would map to shift + exp(x/scale).

Having similar scales helps with optimization and sampling.

(Sorry for opening so many issues today. Your package has been very useful for me!)

UnitVector only transforms to positive hemisphere

The current UnitVector transformation currently transforms an n-1 dimensional vector to an n dimensional unit vector. However, upon transformation, all resulting unit vectors lie on the "hemi-hypersphere", where the final dimension is strictly positive. This is a problem for any applications where the unit vectors across the whole sphere are needed (e.g. directional/orientational statistics)

Stan handles the unit vector transform with a simple trick. I've implemented a version of this with an explanation here and would be happy to adapt it for this package if you're willing.

migrate log density transformations into this package

tpapp/LogDensityProblems.jl#33

This package should depend on LogDensityProblems, not the other way round.

This would remove TransformVariables from the indirect dependency chain of DynamicHMC.

Domain violation in transform

E.g.

julia> t = as((σ = asℝ₊,));

julia> transform(t, [-746])
(σ = 0.0,)

This is causing problems in some optimization code where we optimize over the Vector representation and have to convert back and forth. That causes errors similar to

julia> inverse(t, transform(t, [-746]))
ERROR: DomainError with x > shift must hold. Got
x => 0.0
shift => 0.0:

simplify interface

Instead of transform_at getting a vector and an index, a view would be easier to handle for user extensions. Proposed name: transform_with(flag, t, x).

benchmark current state of the package
rewrite interface
compare benchmarks

Ordered transformations

Hi, have you looked into ordered transformations, as in Stan's approach?

I was experimenting with this, and so far this seems to work ok:

const TV = TransformVariables

struct Ordered{T <: TV.AbstractTransform} <: TV.VectorTransform
    transformation::T
    dim::Int
end

TV.dimension(t::Ordered) = t.dim

addlogjac(ℓ, Δℓ) = ℓ + Δℓ
addlogjac(::TV.NoLogJac, Δℓ) = TV.NoLogJac()


bounds(t::TV.ShiftedExp{true}) = (t.shift, TV.∞)
bounds(t::TV.ShiftedExp{false}) = (-TV.∞, t.shift)
bounds(t::TV.ScaledShiftedLogistic) = (t.shift, t.scale + t.shift)
bounds(::TV.Identity) = (-TV.∞, TV.∞)

# See https://mc-stan.org/docs/2_27/reference-manual/ordered-vector.html
function TV.transform_with(flag::TV.LogJacFlag, t::Ordered, x, index::T) where {T}
    transformation, len = (t.transformation, t.dim)
    @assert dimension(transformation) == 1
    y = similar(x, len)
        
    (lo,hi) = bounds(transformation)

    @inbounds (y[1], ℓ, _) = TV.transform_with(flag, as(Real, lo, hi), x, index)

    @inbounds for i in 2:len
        (y[i], Δℓ, _) =  TV.transform_with(flag, as(Real, y[i-1], hi), x, index)
        ℓ = addlogjac(ℓ, Δℓ)
        index += 1
    end

    return (y, ℓ, index)
end

TV.inverse_eltype(t::Ordered, y::AbstractVector) = TV.extended_eltype(y)

For example,

julia> transform(Ordered(as𝕀, 10), randn(10))
10-element Vector{Float64}:
 0.43442911666628803
 0.6801295759251248
 0.8652960112555127
 0.9378879818399162
 0.97186447061553
 0.992691488683781
 0.9954155804092922
 0.9973011716146748
 0.9984116151813937
 0.9991915044548042

julia> transform(Ordered(asℝ₊, 10), randn(10))
10-element Vector{Float64}:
 0.5238744065026507
 1.0477488130053014
 1.250352594752778
 1.6981546986067726
 2.4349325599881855
 4.314050541989461
 4.56554977859454
 5.699772340085799
 7.099209192367974
 7.7740732481190165

It does sometimes run into trouble, for example

julia> transform(Ordered(as𝕀, 100), randn(100))
ERROR: ArgumentError: the interval (1.0, 1.0) is empty
Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/ArgCheck/5xEDR/src/checks.jl:243 [inlined]
 [2] as(#unused#::Type{Real}, left::Float64, right::Float64)
   @ TransformVariables ~/.julia/packages/TransformVariables/DDsiH/src/scalar.jl:173
 [3] transform_with(flag::TransformVariables.NoLogJac, t::Ordered{TransformVariables.ScaledShiftedLogistic{Float64}}, x::Vector{Float64}, index::Int64)
   @ Main ./REPL[163]:10
 [4] transform(t::Ordered{TransformVariables.ScaledShiftedLogistic{Float64}}, x::Vector{Float64})
   @ TransformVariables ~/.julia/packages/TransformVariables/DDsiH/src/generic.jl:261
 [5] top-level scope
   @ REPL[173]:1

With some more floating point manipulation, I think this could be made to map to non-decreasing sequences, instead of requiring strict monotonicity.

Inverse transformation for empty vectors fails

I am working on an algorithm where I want to define a transformation for a vector of length p. Importantly, p may be zero. This transformation is part of a tuple of transformations. Therefore, the final full parameter vector has a non-zero dimension. What I have done in the past is to omit the zero-dimensional transform, but this leads to additional coding to "compensate" for the missing element. Also, I think it can result in type instability.

I can define a transform for a vector of dimension zero. Also, the forward transform does work:

julia> trans = as(Vector,0)
TransformVariables.ArrayTransform{TransformVariables.Identity,1}(asℝ, (0,))
julia> trans([])
0-element Array{Any,1}

However, the inverse transform fails:

julia> inverse(trans,[])
ERROR: BoundsError: attempt to access 0-element Array{Any,1} at index [1]
Stacktrace:
 [1] getindex at ./array.jl:788 [inlined]
 [2] first at ./abstractarray.jl:323 [inlined]
 [3] inverse_eltype(::TransformVariables.ArrayTransform{TransformVariables.Identity,1}, ::Array{Any,1}) at /home/student/.julia/packages/TransformVariables/gAt5r/src/aggregation.jl:78
 [4] inverse(::TransformVariables.ArrayTransform{TransformVariables.Identity,1}, ::Array{Any,1}) at /home/student/.julia/packages/TransformVariables/gAt5r/src/generic.jl:236
 [5] top-level scope at REPL[47]:1
 [6] eval(::Module, ::Any) at ./boot.jl:331
 [7] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [8] run_backend(::REPL.REPLBackend) at /home/student/.julia/packages/Revise/Pcs5V/src/Revise.jl:1073
 [9] top-level scope at none:0

It would be great if this could work.

Just as a thought - I am not sure if this makes sense in the grand scheme of things - but perhaps this and other cases could be handled by having the option to define the transform for a constant? This would also have a zero-dimensional parameter space, and zero jacobian.

inverse does not check for inputs of special arrays

For scalar transforms inverse is checking the input but not special arrays:

## scalar transform
inverse(as(Real, 5, 10), 5.5)
inverse(as(Real, 5, 10), 1.5)                  # -> error, good
inverse(as(Real, 5, 10), 15.5)                 # -> error, good

## array transform
inverse(UnitVector(3), [sqrt(1/3), sqrt(1/3), sqrt(1/3)])
inverse(UnitVector(3), [0.0, sqrt(1/3), sqrt(1/3)])      # -> no error, bad

inverse(UnitSimplex(3), [0.2, 0.2, 0.6])
inverse(UnitSimplex(3), [0.3, 0.2, 0.6])                 # -> no error, bad

Not sure if this is an omission or design decision.

Type stability and custom transformations

Using forwarddiff to calculate the log jacobian was not inferrable, #28 mitigated this somehow but a better interface would be desirable.

Feature request: Jacobian calculation for gradients

I'd like to use TransformVariables as part of an HMC sampler, but my evaluating my target distribution's probability density involves solving a partial differential equation, which is done outside Julia and thus not compatible with AutoDiff. There are a special techniques for obtaining gradients for such models which involve solving an adjoint equations, but that's a digression.

For people with use cases like mine, it would be nice if there was a function which for a transformation y = t(x) transforms grad(log(f(x)) into grad(log(f(y)) with the appropriate Jacobian correction added to each individual coordinate - similar to transform_logdensity, but for gradients of log-densities.

Nested vector transformations in tuple repeat

Sorry for the cryptic title!

I hope the example below is clearer. It seems that only the first transformations in the tuple a is used and the output repeated:

tt = as((
    a = as((as(Vector, asℝ₊, 2), as(Vector, asℝ₊, 2), as(Vector, asℝ₊, 2))),
    b = as((UnitVector(2), UnitVector(2), UnitVector(2)))
))

y = 1.0:dimension(tt)
x = tt(y)    

# this are all the same!
x.a[1]
x.a[2]
x.a[3]
x.a[1] .== x.a[2] .== x.a[3] # [true, true]

# works as expected
x.b[1]
x.b[2]
x.b[3]
x.b[1] .== x.b[2] .== x.b[3] # [false, false]

Somehow the index seems to be off.

Feature Idea: flat transform

It would be useful if transform would have an option, so that the result remains a flat vector:

transform(t, x, keep_flat=true)

A good use case is converting MCMC samples in MCMCChains.Chains objects:

import MCMCChains

samp = rand(1000, 5)            # we would get them this from a MCMC algorithm

# we get an array of named tuples, which is great to define the model but difficult to convert a `Chain`.
samp_trans1 = mapslices(s -> transform(t, s), samp, dims=2)
MCMCChains.Chains(samp_trans)              # fails

# with the new argument we would get an array
samp_trans2 = mapslice(s -> transform(t, s, keep_flat=true), samp, dims=2)
MCMCChains.Chains(samp_trans2)              # that would work

This seem related to #13

inverse of NamedTuple transformations from different ordering/superset

Currrently inverse for NamedTuple aggregators does not accept different orderings or supersets. Eg

julia> using TransformVariables
    
julia> t = as((a = asℝ, b = asℝ));

julia> inverse(t, (a = 1.0, b = 2.0))
2-element Vector{Float64}:
 1.0
 2.0

julia> inverse(t, (b = 1.0, a = 2.0))
ERROR: ArgumentError: keys(transformations) == keys(y) must hold. Got
keys(transformations) => (:a, :b)
keys(y) => (:b, :a)

julia> inverse(t, (a = 1.0, b = 2.0, c = 3.0))
ERROR: ArgumentError: keys(transformations) == keys(y) must hold. Got
keys(transformations) => (:a, :b)
keys(y) => (:a, :b, :c)

This is OK, but the user should be provided with a method to "regularize" NamedTuples (eg user input may mix up ordering).

allow generalized indexing in arguments

In principle the setup should handle it seamlessly. Test, and clarify in the docs.

TransformVariables: a subset of optics functionality?

Looks like all major functionality of TransformVariables is available as optics/lenses, as defined in Accessors. Optics are very general and composable, so maybe it's worth looking into "replacing" TransformVariables usage with optics. If some useful features aren't covered by optics yet, this most likely can be improved.

Simple example along the lines of TransformVariables docs:

julia> using AccessorsExtra, StaticArrays

# @optics may be upstreamed to Accessors proper at some point
julia> o = @optics log(_.μ) log(_.σ) log(_.τ) _.θs[∗]

julia> y = (μ=1, σ=2, τ=3, θs=SVector(4, 5))

julia> getall(y, o)
(0.0, 0.6931471805599453, 1.0986122886681098, 4, 5)

julia> setall(y, o, getall(y, o) .+ 1)
(μ = 2.718281828459045, σ = 5.436563656918091, τ = 8.154845485377137, θs = [5, 6])

Probability vectors

It would be great to be able to use this library to work with Dirichlet distributions. I have a first crack at a ProbVector transform for this purpose in this gist, but there's a lot I don't understand yet about your library. Does this look like it's headed in the right direction? I'm thinking I can use a CustomTransform to check the math in the implementation.

random transformations

something like

rand(rng::AbstractRNG, t::AbstractTransformation) = transform(t, randn(dimension(t)))

would be useful occasionally.

transform of asℝ₊ variable turns Float32 into Float64

Don't know if this is considered a bug, but I see the following behavior:

> t = as((; α = asℝ₊, β = asℝ ))
> typeof(transform(t)(Float32[1,2]))
NamedTuple{(:α, :β), Tuple{Float64, Float32}}

I was expecting both values to be of Float32 type.

Transform (flatten) as Vector{Matrix}

I have a model with observations occurring at a set of discrete time steps, where at each time one observes a matrix, and it would be great to be able to flatten from this case. More specifically, for the case of a single timestep, I can just follow one of the examples:

t = as((γ = as(Array, length(γ)),)) # identity transformation, just to get the dimension

But for multiple timesteps, it would be great to be able to do

t = as((γ = as(Vector{Array}, length(γ), length(γ[1])),))

or something like that.

(Maybe this relates to #13?)

Would you be able to give me some pointers about how to go about this?

P.S. This is a very nice set of packages!

remove workaround for copy(::TrackedArray)

from transform_with(flag::LogJacFlag, t::ArrayTransform{Identity}, x::RealVector), when FluxML/Flux.jl#416 is fixed.

`inverse` fails for transformation with nested named tuples

The inverse function fails if we have nested Named Tuples. For "normal" tuples it works well:

using TransformVariables  # v0.6.2

# transformation with nested *Named*Tuple
t1 = as((
    a = as(Real, 0, 0.1),
    b = as(
        (b1 = as(Real, 1, 5),
         b2 = as(Real, 10, 50))
    )
))

x = transform(t1, randn(dimension(t1)))
y = inverse(t1, x)              # <- fails
# ERROR: ArgumentError: length(x) == dimension(tt) must hold. Got
# length(x) => 3
# dimension(tt) => 2


# transformation with nested Tuple
t2 = as((
    a = as(Real, 0, 0.1),
    b = as(
        (as(Real, 1, 5),        # no names
         as(Real, 10, 50))
    )
))

x = transform(t2, dimension(t2))
y = inverse(t2, x)              # works :)

make transformations callable

Once JuliaLang/julia#14919 is fixed.

decouple flattening and transformations

Use separate building blocks for

splitting a vector into arrays (or other objects) of various sizes and dimensions, without transforming,
transformations.

Transformations should

have a dimension,
be able to perform log Jacobian calculations,
carry information (schemas, see below) that allows flattening and unflattening their arguments (ignoring extra payload, especially for user-defined transformations).

Transformations should be generic, eg going from StaticArray should also produce StaticArray.

Flattening/unflattening should use schemas, possibly nested.

make this work, define API.
see what's optimizable.

Vectorizing scalar transforms?

It would be nice if I can create "vectorized" versions of scalar transforms, maybe using something like as(Array, as𝕀, 10) such that

transform(as(Array, as𝕀, size(xs)...), xs) == transform.(as𝕀, xs)

Alternative syntaxes may be

as(Array, 0..1, 10)  # using IntervalSets: (..)
as(Array, 10, Real 0, 1)

Is it in the scope of this package?

explain extra "payload" in the manual

Custom transformations can return useful interim results, eg as elements of a named tuple.

[not an issue] transform for a given distribution

Given a distribution d::UnivariateContinuousDistribution from Distributions.jl, I would like to be able get a suitable change of variables and using TransformVariables seems like it provides exactly the right machinery for that. So I thought maybe I could define something like

as(Real, d) = as(Real, support(d).lb, support(d).ub)

but this does not work because Inf and ∞ are not the same. So my questions are

How do you guys go about doing that?
Is this even a good idea/use of TransformVariables, or should I look for something else?

Source links in docs point to julialang

I was browsing through the docs and noticed that if I click the "source" button associated with a docstring, I am taken to julia base instead of to the definition of the function in TransformVariables.jl

Zygote AD fails

import Zygote
using TransformVariables
f(y) = y.mu + y.sigma
t = as((mu  = as(Real, -25, 25), sigma = asℝ₊), )
Zygote.gradient(f ∘ t, zeros(2))
(f∘t)(zeros(2))

(this is the MWE for tpapp/LogDensityProblems.jl#43 (comment))

@itsdfish

dimension mismatch for inverse transformation of nested tuples

The inverse is not working if a named tuple contains (unnamed) tuples:

N = 3
tt = as((
    a = as(Tuple(as(Vector, asℝ₊, 2) for _ in 1:N)),
    b = as(Tuple(UnitVector(n) for n in 1:N)),
))

dimension(tt)  # -> 9

x = tt(randn(dimension(tt)))    # ok

y = inverse(tt, x)
# ERROR: ArgumentError: length(x) == dimension(tt) must hold. Got
# length(x) => 9
# dimension(tt) => 6

Maybe that's not how this package is intended to be used!

My motivation was that while I could use as(Vector, as(Vector, asℝ₊, 2), N) for the the field a, I could not construct a Vector containing different transformations as in field b.

rename package

Perhaps TransformDensity would be more appropriate. See tpapp/ContinuousTransformations.jl#8.

(thanks @simonbyrne)

symmetric positive definite domain

Could you add a symmetric positive definite domain transformation? It could be based on the correlation one.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Autodiff with constrained variable with Zygote

I am having trouble using DynamicHMC + Zygote on a simple example using a constrained variable. This will work with ForwardDiff, but I am interested in Zygote because of the reverse mode autodiff and its efficiency.

Of course, one solution is that Zygote will be expanded to handle this case properly. The reason for raising the problem here is that perhaps a small adjustment in the logdensity could fix this issue?

using Parameters,DynamicHMC,LogDensityProblems,TransformVariables,Random
using Zygote

struct AProblem
end
(p::AProblem)(θ) = -θ[1]^2/2
p = AProblem()
trans = as(Vector,as(Real,-1.0,1.0),1)
P = TransformedLogDensity(trans,p)
∇P = ADgradient(:Zygote,P)
results = mcmc_with_warmup(Random.GLOBAL_RNG,∇P,100)

# Output:
# ERROR: LoadError: Need an adjoint for constructor StepRange{Int64,Int64}. Gradient is of type Array{Nothing,1}
# ...
# [19] logdensity at C:\Users\sdwaele\.julia\packages\LogDensityProblems\2UHo5\src\LogDensityProblems.jl:168 [inlined]
# ...