Fix inverse_eltype for ScalarTransforms

[devmotion]

I had a use case of TransformVariables where I wanted to differentiate with respect to the parameters of a (scalar) transform. Initially this failed due to the default definition of inverse_eltype for ScalarTransforms, which is based solely on the type of the input value.

A safer but more complex - and hence in general slower - default definition would be based on inverse. AFAICT this would be the only way that would be guaranteed to be correct. Alternatively, one could require that users implement inverse_eltype and not provide a default implementation.

On the master branch

julia> using TransformVariables, Chairmarks

julia> @b (as(Real, 0, 2), 1.0) inverse_eltype(_[1], _[2])
1.137 ns

With this PR

julia> using TransformVariables, Chairmarks

julia> @b (as(Real, 0, 2), 1.0) inverse_eltype(_[1], _[2])
7.969 ns

---

Edit: I generally try to avoid Base.promote_op for the reasons listed in its docstring, but it seems to provide a fast alternative for the problematic transforms in TransformVariables. However, it means that the return type of inverse_eltype(::ScalarTransform) might change in different Julia versions, etc. The docstring of inverse_eltype only mentions that it should return a type that can be used for inverse!, so it might not be a problem.

julia> using TransformVariables, Chairmarks

julia> @b (as(Real, 0, 2), 1.0) inverse_eltype(_[1], _[2])
1.138 ns

[tpapp]

Initially this failed due to the default definition of inverse_eltype for ScalarTransforms, which is based solely on the type of the input value.

Can you please share an MWE?

BTW, the primary use case of this package is transform from homogeneous vectors to various things, for use in MCMC etc. I rarely ever use inverse, it is provided just to make things complete and does not need to be fast.

And in any case, transformations are rarely a bottleneck in real life code. Maybe we can dispose of inverse! and not even expose inverse_eltype.

[devmotion]

I added a few simple reproducers to the tests. My example was more complex but a simple cases that errors currently due to inverse_eltype is

julia> d1 = ForwardDiff.derivative(5.3) do x
               return only(inverse(as(Vector, as(Real, x, ∞), 1), [10]))
           end
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#3#4", Float64}, Float64, 1})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(::IrrationalConstants.Logπ)
   @ IrrationalConstants ~/.julia/packages/IrrationalConstants/lWTip/src/macro.jl:131
  ...

Stacktrace:
 [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{var"#3#4", Float64}, Float64, 1})
   @ Base ./number.jl:7
 [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{var"#3#4", Float64}, Float64, 1}, i::Int64)
   @ Base ./array.jl:987
 [3] inverse_at!(x::Vector{Float64}, index::Int64, t::TransformVariables.ShiftedExp{true, ForwardDiff.Dual{…}}, y::Int64)
   @ TransformVariables ~/.julia/packages/TransformVariables/9Eqxa/src/scalar.jl:25
 [4] inverse_at!
   @ ~/.julia/packages/TransformVariables/9Eqxa/src/aggregation.jl:224 [inlined]
 [5] inverse!(x::Vector{…}, transformation::TransformVariables.ArrayTransformation{…}, y::Vector{…})
   @ TransformVariables ~/.julia/packages/TransformVariables/9Eqxa/src/generic.jl:191
 [6] inverse(t::TransformVariables.ArrayTransformation{TransformVariables.ShiftedExp{…}, 1}, y::Vector{Int64})
   @ TransformVariables ~/.julia/packages/TransformVariables/9Eqxa/src/generic.jl:270
 [7] (::var"#3#4")(x::ForwardDiff.Dual{ForwardDiff.Tag{var"#3#4", Float64}, Float64, 1})
   @ Main ./REPL[5]:2
 [8] derivative(f::var"#3#4", x::Float64)
   @ ForwardDiff ~/.julia/packages/ForwardDiff/UBbGT/src/derivative.jl:14
 [9] top-level scope
   @ REPL[5]:1
Some type information was truncated. Use `show(err)` to see complete types.

[tpapp]

Given your MWE, I think that the issue is that you are differentiating the transform, and inverse_eltype(::ScalarTransform, ...) ignores this in the type calculations.

I did not expect this use case, but it should be allowed. This is an issue only with the scalar transforms, since the rest have integer parameters.

I think the quickest way to fix this may be

1. adding a type parameter to ScalarTransform, eg S,

2. taking that account with the relevant inverse_eltype method, eg

function inverse_eltype(t::ScalarTransform{S}, y::T) where {S,T <: Real}
    promote_type(S, float(T))
end

but I did not check this (not at the computer currently).

[devmotion]

Given your MWE, I think that the issue is that you are differentiating the transform, and inverse_eltype(::ScalarTransform, ...) ignores this in the type calculations.

Exactly, that's the problem I tried to address with the PR 🙂

I think the quickest way to fix this may be
adding a type parameter to ScalarTransform, eg S,

I considered this, but I was worried that in general cases beyond the ScalarTransforms in this package it might be annoying/cumbersome to specify the type parameter, in particular if there are multiple and more complex fields involved and/or one's not interested in inverse, so I tried to come up with an approach that does not require downstream changes. I also don't think promote_type is correct in general. For instance, it could lead to very wide types when unitful numbers are involved, e.g. if the parameters of the transform would be unitful numbers and the vector representation would be unitless:

julia> using Unitful

julia> typeof(1.0u"s" * exp(1.0))
Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}

julia> promote_type(typeof(1.0u"s"), Float64)
Quantity{Float64}

[tpapp]

Note that the unit types are not <: Real, so they are outside the domains of tranformation constructors and inverse (and currently fail).

Is there a case where promote_type(::Real, ::Real) fails? If yes, we can always do something like

function promote_reals(::Type{T}, ::Type{S}) where {T,S}
    typeof(exp(one(S) + zero(T)))
end

and hope that the compiler elides this.

[devmotion]

Note that the unit types are not <: Real, so they are outside the domains of tranformation constructors and inverse (and currently fail).

Yes, this was just a hypothetical example to show limitations of such a default definition of inverse_eltype(::ScalarTransform, ::Real).

My current feeling is that either the default should be based on inverse, either directly based on the values (which in general is not elided by the compiler it seems) or based on the types using eg. promote_op, or the default should be removed and definitions for each ScalarTransform in TransformVariables should be added explicitly. In the latter case, downstream users would be required to implement inverse_eltype also for ScalarTransforms, in addition to inverse. I think I'm leaning towards removing the fallback definition.

[tpapp]

First, note that ScalarTransform is not exported or public. It is an abstract type used for internal code organization, and only code within the package is supposed to subtype it. So I think that having a fallback method based on that is fine, with the explicit purpose of avoiding duplication and nothing else. It is also used already for that purpose.

Users may hook into the API of this package by defining their own methods (see the docstring of AbstractTransform), in which case they are required to define inverse_eltype. There is no "catch-all" fallback.

I am inclined to merge this as is. I am aware that Base.promote_op may be fragile in theory, but maybe we can mitigate that by peppering the tests you added with a bunch of @inferred checks.

[tpapp]

Thanks for looking into this and the thorough tests.

[tpapp]

With promote_op, on my machine currently

julia> @b (as(Real, 0, 2), 1.0) inverse_eltype(_[1], _[2])
1.530 ns

so performance-wise this is fine.

[tpapp]

LGTM, thanks!