Add NoLogAbsDetJacobian

[oschulz]

Closes #6, same approach as using in InverseFunctions.

@devmotion would you have time to give this a quick look (I have an immediate use case for it)?

[codecov]

Codecov Report

Merging #7 (360ad86) into master (4ecf0a9) will not change coverage.
The diff coverage is 100.00%.

❗ Current head 360ad86 differs from pull request most recent head acd2a0c. Consider uploading reports for the commit acd2a0c to get more accurate results

@@            Coverage Diff            @@
##            master        #7   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files            2         2           
  Lines           53        53           
=========================================
  Hits            53        53           

Impacted Files	Coverage Δ
src/with_ladj.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4ecf0a9...acd2a0c. Read the comment docs.

[oschulz]

Fixed the tests.

[devmotion]

In principle it seems reasonable to add something that allows to check if the interface is implemented. I'm not sure about the design in this PR though, it seems a bit complicated. Maybe a singleton type would be sufficient, without any custom error messages for some hand-picked methods?

[devmotion]

It seems it would be sufficient to use a singleton type here and remove all fields and type parameters? In InverseFunction, f is included in NoInverse to be able to recover the inverse of NoInverse but here it seems that it's only used for the error messages and hence does not seem strictly necessary.

[oschulz]

I think it would be very valuable to the user to include the function and the argument type. The NoLogAbsDetJacobian may be generated deep in a package used by a package used by the user's application - if the error is just "no ladj available", it will be very hard for the user (especially the non-expert user) to track down the cause. The function may be one the user supplied to a chain of algorithms, or a function derived from it, the user need to know what function was missing volume element information and for which argument type.

[devmotion]

I see. In this case the names, as strings or symbols, would be sufficient though since only they are used in the error message? Hence the type could be non-parameteric which would also remove the need for @nospecialize.

[oschulz]

In this case the names, as strings or symbols, would be sufficient though since only they are used in the error message?

That would actually be really bad for my use case. I need to check if ladj is available and then go down one path of computation or another - but it can be time critical, so I really don't want string generation and mem allocation there.

I would like string handling and message generation to only occur if the exception is thrown, not if the application handles NoLogAbsDetJacobian (that's what it's for, after all).

[devmotion]

Are symbols equally bad?

[oschulz]

Yes, I think Symbol also comes with a cost. And I think NoLogAbsDetJacobian{F,T} will be simple, cheap and effective.

[oschulz]

Also, in applications with autodiff, any code run during the generation of the NoLogAbsDetJacobian instance would lead to additional costs because Zygote would process it (or we'd need to a a custom rrule to ignore it).

[devmotion]

This selection of methods seems a bit arbitrary, why are these specifically defined?

In any case I think it would maybe be a bit simpler to use @eval or to use a function that throws the error. Creating the arguments and splatting them seems a bit complicated.

[oschulz]

This selection of methods seems a bit arbitrary, why are these specifically defined?

I thought it would be good if the typical ways (defined as test cases in "test_with_ladj.jl") most users would use the result of r = with_logabsdet_jacobian(f,x) will throw the exception, so

• r[1], r[2] -> getindex
• first(r), y, ladj = r -> iterate
• last(r) -> lastindex

firstindex I added because it seemed weird to only specialize lastindex.

maybe be a bit simpler to use @eval

I used @eval in my first draft, but getindex takes additional arguments and the other's don't, it seemed easiest and most readable like this in the end.

use a function that throws the error. Creating the arguments and splatting them seems a bit complicated.

I thought it would be most user-friendly to have the exception originating directly from the point of call.

[devmotion]

This selection of methods seems a bit arbitrary, why are these specifically defined?

I thought it would be good if the typical ways (defined as test cases in "test_with_ladj.jl") most users would use the result of r = with_logabsdet_jacobian(f,x) will throw the exception, so

• r[1], r[2] -> getindex
• first(r), y, ladj = r -> iterate
• last(r) -> lastindex

firstindex I added because it seemed weird to only specialize lastindex.

I can see the motivation but I'm still not completely sure about it. It reminds me of the tangent types in CRC which allow some operations but not others depending on the type and throw errors otherwise - but my impression has always been that it's a bit of a mess.

I think one problem is that the suggestion here mixes two different use cases:

• Allowing developers to dispatch on whether with_logabsdet_jacobian is defined
• Throwing reasonably informative errors if with_logabsdet_jacobian is not defined for a user provided function

It seems for the first use case the type information is not needed (since no errors are thrown) and in the second case it might be clearer if the errors are thrown right away instead of hoping that users call one of the functions for which error messages are implemented (in other cases less informative MethodImplemented errors will be thrown or some undesirable fallback is used).

maybe be a bit simpler to use @eval
I used @eval in my first draft, but getindex takes additional arguments and the other's don't, it seemed easiest and most readable like this in the end.

I see. I think it would be possible though to handle this case, e.g. by something like

for f in (...)
    args = f === Base.getindex ? :(x::Vararg) : :()
    @eval begin
        function $f(x::NoAbsLogDetJacobian, $args)
            error(...)
...
end

[devmotion]

Maybe we could

• add a non-parameteriv singleton type NoLogAbsDetJacobian and
• add a function

  function check_with_logabsdet jacobian(f, x)
      y_ladj = with_logabsdet_jacobian(f, x)
      if y_ladj isa NoLogAbsDetJacobian
          error(...)
      end
      return y_ladj
  end

to cover both use cases.

[oschulz]

I can see the motivation but I'm still not completely sure about it. It reminds me of the tangent types in CRC which allow some operations but not others

I don't think it's the same here - NoAbsLogDetJacobian doesn't support any operations, it's just throws a more user friendly exception for those operations that users are likely to use on it. I would argue that in this specific cases that's a very finite list.

I think one problem is that the suggestion here mixes two different use cases [...] It seems for the first use case the type information is not needed (since no errors are thrown)

We did take the same approach for InverseFunctions as in this PR - return NoInverse, with information about the function. Allowing inversion of that to recover the original function was actually a secondary consideration if I remember correctly.

I think the current approach is quite simple and effective

I see. I think it would be possible though to handle this case, e.g. by something like for f in ...

I'm not strictly against it, but will it really be easier to read/maintain than the current four lines of code in this PR? In JuliaMath/InverseFunctions.jl#16 we decided to go with explicit instead of generated code because that wouldn't have saved much.

[devmotion]

I think the situation here is different from the one in InverseFunctions since there users will always peform a specific operations (namely calling the NoInverse) and hence it is clear for which method an error should be thrown and that all possible cases are covered.

[oschulz]

Ok, then how about now throwing any exceptions (I've updated the PR accordingly)? The user will now see this:

julia> _, _ = with_logabsdet_jacobian(sum, 5)
ERROR: MethodError: no method matching iterate(::NoLogAbsDetJacobian{typeof(sum), Int64})

I think that's pretty clear in itself, without a fancy error string.

[oschulz]

Forgot to remove _no_ladj_errmsg, fixed now.

[oschulz]

I've removed the content of NoLogAbsDetJacobian, turning it into a singleton, but with type parameters. I think you're right, it's different from inverse in that the return value of with_logabsdet_jacobian is typically split up directly and won't propagate through a call stack unused. This way we can keep things simple and the user will still get an informative error:

julia> function foo(f, x)
           y, ladj = with_logabsdet_jacobian(f, x)
           cos(y  + ladj)
       end
foo (generic function with 1 method)

julia> foo(sum, (2, 3))
ERROR: MethodError: no method matching iterate(::NoLogAbsDetJacobian{typeof(sum), Tuple{Int64, Int64}})

I think providing an additional function check_with_logabsdet would be unnecessarily complicated and also problematic - package code might often just use with_logabsdet_jacobian, expecting it to be defined, and then the user won't get an informative error if it's not (and they won't have control over the code at the point of call of with_logabsdet_jacobian).

I don't see the harm in NoLogAbsDetJacobian having type parameters. I don't think there's a cost and in any case - as the return value of with_logabsdet_jacobian will always almost be used directly (either via getindex or iterate, or if the code checks it, via isa NoLogAbsDetJacobian) - the compiler should elide generation of theNoLogAbsDetJacobian object. Here's a dummy example:

julia> function bar(f, x)
           r = with_logabsdet_jacobian(f, x)
           if r isa NoLogAbsDetJacobian
               sin(f(x))
           else
               y, ladj = r
               cos(y  + ladj)
           end
       end
bar (generic function with 1 method)

julia> @code_typed bar(sum, (2, 3))
CodeInfo(
1 ─ %1 = Core.getfield(x, 1)::Int64
│   %2 = Core.getfield(x, 2)::Int64
│   %3 = Base.add_int(%1, %2)::Int64
│   %4 = Base.sitofp(Float64, %3)::Float64
│   %5 = invoke Base.Math.sin(%4::Float64)::Float64
└──      return %5
) => Float64

julia> @code_typed bar(log, 2)
CodeInfo(
1 ─ %1 = Base.sitofp(Float64, x)::Float64
│   %2 = invoke Base.Math._log(%1::Float64, $(QuoteNode(Val{:ℯ}()))::Val{:ℯ}, :log::Symbol)::Float64
│   %3 = Base.neg_float(%2)::Float64
└──      goto #3 if not false
2 ─      nothing::Nothing
3 ┄ %6 = Base.add_float(%2, %3)::Float64
│   %7 = invoke Main.cos(%6::Float64)::Float64
└──      return %7
) => Float64

@devmotion is this approach Ok with you?

[devmotion]

I guess probably it's easier to start without the (somewhat arbitrary) selection of exceptions and improve error handling later if needed (possibly with a different function?).

[oschulz]

I guess probably it's easier to start without the (somewhat arbitrary) selection of exceptions

@devmotion, dhould I re-add the exceptions, or do we stay with the

julia> _, _ = with_logabsdet_jacobian(sum, 5)
ERROR: MethodError: no method matching iterate(::NoLogAbsDetJacobian{typeof(sum), Int64})

error response for now?

[devmotion]

LGTM 👍

[oschulz]

Ok, then I'll merge #7 and #8 as-is and release them together?