Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

fix xexpy(0, 1000) #80

Merged
merged 5 commits into from
Feb 16, 2024
Merged

fix xexpy(0, 1000) #80

merged 5 commits into from
Feb 16, 2024

Conversation

oscardssmith
Copy link
Contributor

now it equals 0.

@cossio
Copy link
Contributor

cossio commented Feb 14, 2024

Fixes #79

@cossio
Copy link
Contributor

cossio commented Feb 14, 2024

Maybe add xexpy(0, 1000) as a test case?

@oscardssmith
Copy link
Contributor Author

tests added.

test/basicfuns.jl Outdated Show resolved Hide resolved
@cossio
Copy link
Contributor

cossio commented Feb 14, 2024

Could we get a release with this?

@devmotion
Copy link
Member

The PR does not match what xexpy is intended to return: #79 (comment)

@oscardssmith
Copy link
Contributor Author

@devmotion I don't understand what you mean?

@devmotion
Copy link
Member

xexpy(0, Inf) is supposed to (as documented in the docstring) return 0 * exp(Inf) which is NaN. To me it seems there is no bug to fix here.

To elaborate a bit more, xexpy and xlogy are supposed to satisfy the invariant xexpy(a, b) = xlogy(exp(b), exp(a)): If b = -Inf, or equivalently exp(b) = 0, then xexpy(a, b) = xlogy(exp(b), exp(a)) = 0; otherwise they return a * exp(b) = exp(b) * log(exp(a)).

@oscardssmith
Copy link
Contributor Author

oscardssmith commented Feb 14, 2024

This PR maintains xexpy(0, Inf) == NaN. The difference is in xexpy(0, 10000) where previously exp(10000) overflowed. See the added tests:

    @test xexpy(0, 1000) == 0.0
    @test isnan(xexpy(0, Inf))
    @test isnan(xexpy(0, NaN))

@cossio
Copy link
Contributor

cossio commented Feb 14, 2024

xexpy and xlogy are supposed to satisfy the invariant xexpy(a, b) = xlogy(exp(b), exp(a))

Why? This requirement is not documented.

I would also note that xlogy(1e100, 1) == 0, while xlogy(exp(1000), exp(0)) gives NaN only because the intermediate exp(1000) overflows.

@devmotion
Copy link
Member

Why? This requirement is not documented.

The behaviour of each function is documented which implicitly documents the invariance.

I would also note that xlogy(1e100, 1) == 0, while xlogy(exp(1000), exp(0)) gives NaN only because the intermediate exp(1000) overflows.

My comment was supposed to be understood as a exact/mathematical description of the intended functions and their relation, disregarding floating-point issues such as overflow. That is, the invariance is for the mathematical functions underlying xlogy and xexpy, but not their Julia implementations.

@oscardssmith
Copy link
Contributor Author

In that case, isn't this PR correct? This PR leaves xexpy(0, Inf) as NaN.

Copy link
Member

@devmotion devmotion left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I missed the check for isfinite(y), so I think generally the PR is fine. It still disagrees with the docstring though since due to the change xexpy(x, y) is not always x * exp(y) if y != -Inf anymore. Can you adjust the docstring (e.g., by changing if `y == -Inf` to if `y == -Inf` or if `x == 0.0` and `y` is finite.) and maybe also add the xexpy(0, 1000) example as a doctest?

@oscardssmith
Copy link
Contributor Author

I interpret the docstring as describing the function mathematically, in which case it is already correct.

@devmotion
Copy link
Member

Generally, the docstrings in LogExpFunctions mention explicitly if the computation deviates from the given Julia code ("preserving numerical accuracy", "carefully evaluated", "avoiding over/underflow" ...), so IMO the docstring should either mention the additional special case or contain a comment that indicates that the computation is not just x * exp(y).

@oscardssmith
Copy link
Contributor Author

Ideally this code could avoid more overflow to deal with cases like xexpy(floatmin(Float64), 1000.0), but that is rather difficult.

src/basicfuns.jl Outdated Show resolved Hide resolved
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
@cossio
Copy link
Contributor

cossio commented Feb 15, 2024

Could this PR include a version bump? I would like to use this updated xexpy.

@devmotion devmotion merged commit 27b44d6 into JuliaStats:master Feb 16, 2024
4 checks passed
@oscardssmith oscardssmith deleted the patch-1 branch February 16, 2024 14:21
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants