"""Log-likelihood of the multivariate exponential-kernel Hawkes process.
This module is a typed, documented face over the compiled core. It does the
shape checks the compiled core cannot express conveniently and then hands the
arrays over unchanged; it performs no arithmetic of its own, on purpose, so that
the number a user gets is the number the C++ tests verified.
"""
from __future__ import annotations
from typing import TypedDict
import numpy as np
import numpy.typing as npt
from hawkes_engine import _core
from hawkes_engine._arrays import as_real_1d, as_real_2d, as_uint8_1d, flat
__all__ = ["LoglikResult", "loglik"]
[docs]
class LoglikResult(TypedDict, total=False):
"""Return value of :func:`loglik`.
``value`` is always present. The gradient blocks are present only when
``grad=True``: ``dmu`` has shape ``(M,)``, ``dalpha`` and the pair ``dbeta``
have shape ``(M, M)``, and the shared-decay ``dbeta`` has shape ``(M,)``.
When the intensity is not strictly positive at some event, ``value`` is
``-inf`` and every gradient entry is ``nan``. Nothing is clipped: a fit that
reaches that region has found something wrong, and it should be able to see
it.
"""
value: float
dmu: npt.NDArray[np.float64]
dalpha: npt.NDArray[np.float64]
dbeta: npt.NDArray[np.float64]
def _prepare(
times: npt.ArrayLike,
types: npt.ArrayLike,
mu: npt.ArrayLike,
alpha: npt.ArrayLike,
beta: npt.ArrayLike,
shared_decay: bool,
) -> tuple[
npt.NDArray[np.float64],
npt.NDArray[np.uint8],
npt.NDArray[np.float64],
npt.NDArray[np.float64],
npt.NDArray[np.float64],
int,
]:
"""Convert every argument to the form the core reads, and check the shapes.
The conversion is deliberate and happens here, once, where the user can see
it in a traceback. The core itself refuses anything it would have to copy,
so a silent conversion is impossible from this entry point.
"""
t = as_real_1d(times, "times")
k = as_uint8_1d(types, "types")
m = as_real_1d(mu, "mu")
a = as_real_2d(alpha, "alpha")
dim = int(m.shape[0])
if a.shape != (dim, dim):
raise ValueError(f"alpha must have shape ({dim}, {dim}), got {a.shape}")
expected_beta = (dim,) if shared_decay else (dim, dim)
b = as_real_2d(beta, "beta") if not shared_decay else as_real_1d(beta, "beta")
if b.shape != expected_beta:
raise ValueError(f"beta must have shape {expected_beta}, got {b.shape}")
if t.shape[0] != k.shape[0]:
raise ValueError(
f"times and types must have the same length, got {t.shape[0]} and {k.shape[0]}"
)
return t, k, m, flat(a), flat(b), dim
[docs]
def loglik(
times: npt.ArrayLike,
types: npt.ArrayLike,
mu: npt.ArrayLike,
alpha: npt.ArrayLike,
beta: npt.ArrayLike,
T: float,
*,
grad: bool = False,
shared_decay: bool = False,
) -> LoglikResult:
"""Evaluate the log-likelihood of a multivariate Hawkes process.
Implements ``docs/spec/mathematics.md`` section 3:
.. math::
\\ell(\\theta) = \\sum_i \\left[ \\sum_k \\log \\lambda_i(t^i_k)
- \\int_0^T \\lambda_i(t)\\,dt \\right]
Parameters
----------
times
Event times, shape ``(N,)``, non-decreasing, seconds from the start of
the session. Convert epoch timestamps to a session-relative origin
before calling (ADR-003).
types
Type index of each event, shape ``(N,)``, in ``0..M-1``.
mu
Baseline intensities, shape ``(M,)``, non-negative.
alpha
Excitation block, shape ``(M, M)``, non-negative. Row is the target
type, column the source (ADR-002).
beta
Decay rates, shape ``(M, M)`` for a distinct decay per pair, or ``(M,)``
when ``shared_decay`` is true, in which case ``beta_ij = beta_j``. All
entries strictly positive.
T
End of the observation window. Must be at least the last event time.
grad
Also compute the gradient with respect to the natural-space parameters.
shared_decay
Use the source-only decay parameterization, O(N*M) instead of O(N*M*M).
Returns
-------
LoglikResult
``value`` and, when requested, ``dmu``, ``dalpha``, ``dbeta``.
Raises
------
ValueError
If the events are not sorted or lie outside ``[0, T]``, a parameter is
outside its domain, or a shape disagrees with ``M``.
Notes
-----
The gradient is in natural space. The optimizer's log-space gradient is
``natural * parameter`` entrywise; that conversion is applied once, in the
objective wrapper, not here (ADR-002).
Two events at the same timestamp do not excite each other: the intensity
sums over strictly earlier events.
"""
t, k, m, a, b, dim = _prepare(times, types, mu, alpha, beta, shared_decay)
out = _core.loglik(t, k, m, a, b, float(T), dim, shared_decay, grad)
result: LoglikResult = {"value": float(out["value"])}
if grad:
result["dmu"] = out["dmu"]
result["dalpha"] = out["dalpha"]
result["dbeta"] = out["dbeta"]
return result
def loglik_grad_log_space(
times: npt.ArrayLike,
types: npt.ArrayLike,
mu: npt.ArrayLike,
alpha: npt.ArrayLike,
beta: npt.ArrayLike,
T: float,
*,
shared_decay: bool = False,
) -> LoglikResult:
"""Log-likelihood with the gradient in log-space coordinates.
Same arguments as :func:`loglik` with ``grad=True``; the gradient blocks are
``d/d(log x)``, not ``d/dx``. This is the objective the optimizer steps in
(ADR-002); it is not part of the public surface.
"""
t, k, m, a, b, dim = _prepare(times, types, mu, alpha, beta, shared_decay)
out = _core.loglik_grad_log_space(t, k, m, a, b, float(T), dim, shared_decay)
return {
"value": float(out["value"]),
"dmu": out["dmu"],
"dalpha": out["dalpha"],
"dbeta": out["dbeta"],
}