From 53efea2efc537e840c479f62e801b344dafdc80c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= <tkpapp@gmail.com>
Date: Mon, 9 Jan 2023 14:20:21 +0100
Subject: [PATCH] Add logdensity_gradient_and_hessian.

Fixes #65.
---
 docs/src/index.md         |  3 ++-
 src/LogDensityProblems.jl | 38 +++++++++++++++++++++++++++++++++-----
 src/utilities.jl          |  5 +++--
 3 files changed, 38 insertions(+), 8 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index 82371f9..babf1c1 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -102,7 +102,7 @@ LogDensityProblems.logdensity(ℓ, zeros(2))
 ```
 
 !!! note
-    Before running time-consuming algorithms like MCMC, it is advisable to test and benchmark your log density evaluations separately. The same applies to [`LogDensityProblems.logdensity_and_gradient`](@ref).
+    Before running time-consuming algorithms like MCMC, it is advisable to test and benchmark your log density evaluations separately. The same applies to [`LogDensityProblems.logdensity_and_gradient`](@ref) and [`LogDensityProblems.logdensity_gradient_and_hessian`](@ref).
 
 ## Manual unpacking and transformation
 
@@ -184,4 +184,5 @@ LogDensityProblems.LogDensityOrder
 LogDensityProblems.dimension
 LogDensityProblems.logdensity
 LogDensityProblems.logdensity_and_gradient
+LogDensityProblems.logdensity_gradient_and_hessian
 ```
diff --git a/src/LogDensityProblems.jl b/src/LogDensityProblems.jl
index 1d90208..a24ed52 100644
--- a/src/LogDensityProblems.jl
+++ b/src/LogDensityProblems.jl
@@ -45,9 +45,10 @@ $(SIGNATURES)
 
 Test if the type (or a value, for convenience) supports the log density interface.
 
-When `nothing` is returned, it doesn't support this interface.  When `LogDensityOrder{K}()`
-is returned (typically with `K == 0` or `K = 1`), derivatives up to order `K` are supported.
-*All other return values are invalid*.
+When `nothing` is returned, it doesn't support this interface. When
+`LogDensityOrder{K}()` is returned (typically with `K == 0`, `K = 1`,
+or `K == 2`), derivatives up to order `K` are supported. *All other
+return values are invalid*.
 
 # Interface description
 
@@ -59,6 +60,8 @@ The following methods need to be implemented for the interface:
 
 3. [`logdensity_and_gradient`](@ref) when `K ≥ 1`.
 
+4. [`logdensity_gradient_and_hessian`](@ref) when `K ≥ 2`.
+
 See also [`LogDensityProblems.stresstest`](@ref) for stress testing.
 """
 capabilities(T::Type) = nothing
@@ -105,8 +108,8 @@ function logdensity end
 """
     logdensity_and_gradient(ℓ, x)
 
-Evaluate the log density `ℓ` and its gradient at `x`, which has length compatible with its
-[`dimension`](@ref).
+Evaluate the log density `ℓ` and its gradient at `x`, which has length
+compatible with its [`dimension`](@ref).
 
 Return two values:
 
@@ -124,6 +127,31 @@ The first argument (the log density) can be shifted by a constant, see the note
 """
 function logdensity_and_gradient end
 
+"""
+    logdensity_gradient_and_hessian(ℓ, x)
+
+Evaluate the log density `ℓ`, its gradient, and Hessian at `x`, which
+has length compatible with its [`dimension`](@ref).
+
+Return three values:
+
+- the log density as real number, which equivalent to `logdensity(ℓ, x)`
+
+- *if* the log density is finite, the gradient, an `::AbstractVector` of real numbers,
+   otherwise this value is arbitrary and should be ignored.
+
+- *if* the log density is finite, the Hessian, an `::AbstractMatrix` of real numbers,
+   otherwise this value is arbitrary and should be ignored.
+
+!!! note
+    Caller may assume ownership of results, ie that the gradient and
+    Hessian will not be overwritten or reused for a different purpose.
+
+The first argument (the log density) can be shifted by a constant, see the note for
+[`logdensity`](@ref).
+"""
+function logdensity_gradient_and_hessian end
+
 include("utilities.jl")
 
 end # module
diff --git a/src/utilities.jl b/src/utilities.jl
index acc8864..0444169 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -37,8 +37,9 @@ Test `ℓ` with random values.
 `N` random vectors are drawn from a standard multivariate Cauchy distribution, scaled with
 `scale` (which can be a scalar or a conformable vector).
 
-Each random vector is then used as an argument in `f(ℓ, ...)`. `logdensity` and
-`logdensity_and_gradient` are  recommended for `f`.
+Each random vector is then used as an argument in `f(ℓ, ...)`. [`logdensity`](@ref),
+[`logdensity_and_gradient`](@ref), and [`logdensity_gradient_and_hessian`](@ref) are
+recommended for `f`.
 
 In case the call produces an error, the value is recorded as a failure, which are returned
 by the function.