Distribution.h

Go to the documentation of this file.

/**
 * @file Distribution.h
 * @brief Per-family polymorphic view over a continuous @c gate_rv
 *        distribution (§F.1 class hierarchy).
 *
 * Every RV evaluator dispatches through this interface (one virtual call
 * per family) and the registries below; family identity is the interned
 * @c DistributionFamily descriptor each implementation file registers --
 * there is no family enum anywhere.  A @c Distribution is a transient
 * view constructed from a parsed @c DistributionSpec via
 * @ref provsql::makeDistribution; the on-disk @c extra text encoding stores the
 * family's name token.
 *
 * Adding a family is one new self-registering implementation file under
 * @c src/distributions/ (plus its SQL constructor).  Methods a family
 * cannot answer follow the NaN-as-undecided contract so the analytic
 * paths fall back to Monte Carlo unchanged.
 */
#ifndef PROVSQL_DISTRIBUTION_H
#define PROVSQL_DISTRIBUTION_H
 
#include <cstddef>
#include <memory>
#include <optional>
#include <random>
#include <string>
#include <utility>
#include <vector>
 
#include "../RandomVariable.h"  // DistributionSpec
 
namespace provsql {
 
/** @brief A closed support interval [lo, hi] (±infinity for unbounded). */
struct DistSupport {
  double lo;
  double hi;
};
 
/**
 * @brief Abstract per-family continuous distribution.
 *
 * Concrete subclasses (Normal / Uniform / Exponential / Erlang / Gamma)
 * hold the two parameters and implement the family-specific closed forms.
 * The methods below are the family-local ones (pairwise closure /
 * comparison rules live in separate registries, not here).
 */
class Distribution {
public:
  virtual ~Distribution() = default;
 
  /** @name Identity (parameters, for the i.i.d. equality test) */
  ///@{
  /**
   * @brief The family's interned registry descriptor.
   *
   * One instance per family, so @c &family() (or comparing
   * @c family().name) is family identity.
   */
  virtual const DistributionFamily &family() const = 0;
  virtual double p1() const = 0;
  virtual double p2() const = 0;
  ///@}
 
  /** @name Closed-form moments */
  ///@{
  virtual double mean() const = 0;                 ///< E[X]
  virtual double variance() const = 0;             ///< Var(X)
  virtual double rawMoment(unsigned k) const = 0;  ///< E[X^k]
  ///@}
 
  /**
   * @brief Whether @c mean() is an affine function of the raw parameters
   *        @c (p1, p2) -- i.e. @c mean = c0 + c1·p1 + c2·p2.
   *
   * Default @c false.  A family whose mean is a linear combination of its
   * parameters (Normal @c μ, Uniform @c (a+b)/2, inverse-Gaussian @c μ)
   * overrides to @c true.  This is the authoritative signal the latent-
   * variable evaluator uses to keep the expectation of a compound (token-
   * parameterised) leaf EXACT: when the mean is affine,
   * @f$E[X] = E[\mathrm{mean}(\theta)] = \mathrm{mean}(E[\theta])@f$ by
   * linearity of expectation (no independence assumption), so the evaluator
   * recurses into the parameter wires instead of falling back to Monte
   * Carlo.  Families with a nonlinear mean (Exponential @c 1/λ, Gamma
   * @c k/λ, ...) keep the default and decline.
   */
  virtual bool meanIsAffine() const { return false; }
 
  /**
   * @brief Whether the family is discrete (integer-valued), so a point
   *        event @c X @c = @c c carries positive mass.
   *
   * Default @c false (continuous).  The discrete families (Poisson,
   * Binomial) override to @c true.  It is the authoritative, per-family
   * signal -- not inferred -- that a point observation / conditioning
   * event on a leaf of this family is @b feasible (a continuous point
   * event is measure-zero, so it is only meaningful as a likelihood
   * weight, never as a rejection interval).  @c pdf() already returns the
   * pmf for a discrete family, so the density-weight path needs no special
   * case beyond this classification.
   */
  virtual bool isDiscrete() const { return false; }
 
  /** @name Density / distribution */
  ///@{
  virtual double pdf(double x) const = 0;  ///< f(x); NaN if the family declines
  virtual double cdf(double x) const = 0;  ///< F(x); NaN if the family declines
  ///@}
 
  /** @brief Natural support interval of X. */
  virtual DistSupport support() const = 0;
 
  /**
   * @brief Finite window [lo, hi] covering essentially all of X's mass, for
   *        numerical quadrature.  Returns false (leaving @p lo / @p hi
   *        untouched) when the parameters are degenerate.
   */
  virtual bool integrationRange(double &lo, double &hi) const = 0;
 
  /**
   * @brief Plot x-window given optional truncation bounds (±infinity means
   *        "unbounded on that side"); used by the SVG curve renderer.
   */
  virtual std::pair<double, double> plotRange(double trunc_lo,
                                              double trunc_hi) const = 0;
 
  /** @brief Draw one sample using the shared MC generator. */
  virtual double sample(std::mt19937_64 &rng) const = 0;
 
  /**
   * @brief Inverse CDF @f$Q(p) = F^{-1}(p)@f$ for @f$p \in (0, 1)@f$.
   *
   * @c std::nullopt when the family has no elementary inverse CDF
   * (Erlang), so "unsupported" cannot be silently consumed -- callers
   * branch on availability explicitly.  Normal uses the
   * Beasley-Springer-Moro rational approximation (~1e-7 accuracy, ample
   * for sampling); Uniform / Exponential invert exactly.  Behaviour at
   * the endpoints is the family's (BSM diverges; clamp @p p strictly
   * inside @c (0, 1) before calling).
   */
  virtual std::optional<double> quantile(double p) const {
    return std::nullopt;
  }
 
  /**
   * @brief Closed-form raw moment @f$E[X^k \mid lo < X < hi]@f$ of the
   *        distribution truncated to @c [lo, hi] (±infinity for a
   *        semi-infinite side).
   *
   * @c std::nullopt when the family has no closed form (Erlang: needs
   * the regularised lower incomplete gamma) or the truncation is
   * degenerate (mass below a numerical floor), so the caller falls
   * through to MC rejection.  @p k is at least 1 (the caller
   * short-circuits @c k @c = @c 0).
   */
  virtual std::optional<double> truncatedRawMoment(double lo, double hi,
                                                   unsigned k) const {
    (void) lo; (void) hi; (void) k;
    return std::nullopt;
  }
 
  /**
   * @brief Draw @p n rejection-free samples of X conditioned on
   *        @c lo < X < hi (±infinity for a semi-infinite side).
   *
   * Each family uses its own exact scheme (interval intersection for
   * Uniform, memorylessness / stable inverse-CDF for Exponential,
   * inverse-CDF transform for Normal).  @c std::nullopt when the family
   * has no scheme (Erlang: needs the inverse regularised incomplete
   * gamma) or the truncation is degenerate / parameters invalid, so the
   * caller's MC-rejection fallback can emit its usual diagnostics.
   */
  virtual std::optional<std::vector<double>> sampleTruncated(
    std::mt19937_64 &rng, double lo, double hi, unsigned n) const {
    (void) rng; (void) lo; (void) hi; (void) n;
    return std::nullopt;
  }
 
  /**
   * @brief Closed-form mean of the max / min of @p n i.i.d. copies of X.
   *
   * - Uniform(a, b): @c E[max] @c = @c a + (b-a)·n/(n+1),
   *   @c E[min] @c = @c a + (b-a)/(n+1).
   * - Exponential(λ): @c E[min] @c = @c 1/(nλ), @c E[max] @c = @c H_n/λ
   *   (harmonic number), for @c λ > 0.
   *
   * @c std::nullopt when the family has no elementary order-statistic
   * mean (Normal, Erlang -- those need the 1-D layer-cake quadrature) or
   * the parameters are degenerate, so the caller falls back to
   * quadrature / Monte Carlo.
   */
  virtual std::optional<double> iidOrderStatMean(std::size_t n,
                                                 bool isMax) const {
    (void) n; (void) isMax;
    return std::nullopt;
  }
 
  /**
   * @brief Closed-form location-scale transform @c a·X @c + @c b within
   *        the family.
   *
   * Returns the transformed distribution when the family is closed under
   * these coefficients, @c nullptr when it is not: Exponential / Erlang
   * decline @c a @c <= @c 0 (the support flips) and any non-zero offset
   * (a shifted Erlang leaves the family); every family declines
   * @c a @c == @c 0 (a Dirac is not in-family -- the constant-fold path
   * is responsible for pure constants).
   */
  virtual std::unique_ptr<Distribution> affine(double a, double b) const = 0;
 
  /** @brief Closed-form @c c·X (affine with no offset). */
  std::unique_ptr<Distribution> scale(double c) const {
    return affine(c, 0.0);
  }
 
  /** @brief Closed-form @c -X (affine with coefficient -1). */
  std::unique_ptr<Distribution> negate() const {
    return affine(-1.0, 0.0);
  }
 
  /**
   * @brief The on-disk @c gate_rv @c extra text encoding
   *        (e.g. <tt>"normal:2.5,0.5"</tt>), inverse of
   *        @c parse_distribution_spec.
   */
  virtual std::string serialise() const = 0;
 
  /**
   * @brief The point-mass value when the parameters make the
   *        distribution degenerate (a Dirac), @c std::nullopt for a
   *        proper distribution.
   *
   * Only Normal reports one (σ == 0, e.g. after an underflowing scale
   * fold); the simplifier collapses such a result to a plain constant.
   * Families whose degenerate forms are still sampled as-is
   * (Uniform(a, a)) do not report.
   */
  virtual std::optional<double> asDirac() const { return std::nullopt; }
};
 
/**
 * @brief Construct the per-family @c Distribution for a parsed spec.
 *
 * Returns @c nullptr only for a spec with a null family pointer (which
 * @c parse_distribution_spec never produces).  Parameter-validity guards
 * live in the family methods (e.g. @c pdf returns NaN for a
 * non-positive σ).
 */
std::unique_ptr<Distribution> makeDistribution(const DistributionSpec &spec);
 
/**
 * @name DistributionRegistry – family descriptor table
 *
 * Family implementations self-register their descriptor at static
 * initialisation: the on-disk name token (the part before the colon in a
 * @c gate_rv @c extra), the parameter count, the display metadata, and
 * the constructor.  @c makeDistribution and @c parse_distribution_spec
 * dispatch through this table, so adding a family means one new
 * implementation file with one registrar -- no existing factory, parser,
 * or header is touched.
 */
///@{
 
/** @brief Construct a family instance from its (up to) two parameters. */
using DistributionFactory =
  std::unique_ptr<Distribution> (*)(double p1, double p2);
 
/**
 * @brief A registered family's descriptor: its complete identity.
 *
 * @c name is the on-disk token (the part before the colon in a
 * @c gate_rv @c extra); @c factory constructs an instance from the
 * parsed parameters; @c label / @c param_names are the display metadata
 * UI clients (ProvSQL Studio) read through @c provsql.rv_families() --
 * a short glance-recognisable glyph for a node circle (e.g. "N", "Exp",
 * "Γ") and the conventional parameter symbols in @c extra order
 * (e.g. {"μ", "σ"}; the second entry is @c nullptr for a 1-parameter
 * family).  Purely presentational choices beyond these (colours,
 * geometry) stay client-side.
 *
 * Each family implementation file defines exactly ONE descriptor with
 * static storage duration and registers its address, so the pointer is
 * an interned family identity (this is what @c DistributionSpec and
 * @c Distribution::family() carry; there is no family enum to extend).
 */
struct DistributionFamily {
  const char *name;
  unsigned nparams;   ///< 1 or 2 (a 1-parameter family leaves p2 = 0)
  const char *label;
  const char *param_names[2];
  DistributionFactory factory;
};
 
/** @brief Register a family; called by the registrar at static init. */
void registerDistributionFamily(const DistributionFamily &descriptor);
 
/** @brief Static-initialisation helper: one per family implementation. */
struct DistributionFamilyRegistrar {
  explicit DistributionFamilyRegistrar(const DistributionFamily &descriptor) {
    registerDistributionFamily(descriptor);
  }
};
 
/**
 * @brief Look up a family by its on-disk name token.
 *
 * Used by @c parse_distribution_spec to resolve the interned
 * descriptor and expected parameter count; @c nullptr for an unknown
 * name.
 */
const DistributionFamily *lookupDistributionFamily(const std::string &name);
 
/**
 * @brief Every registered family, sorted by name token.
 *
 * Backs the @c provsql.rv_families() catalog function.
 */
std::vector<const DistributionFamily *> listDistributionFamilies();
 
///@}
 
/**
 * @brief Numeric inverse CDF: monotone bisection of @c cdf() over the
 *        family's integration window.
 *
 * Family-agnostic fallback for families whose @c quantile() declines
 * (no elementary inverse CDF -- Erlang, Gamma): needs only @c cdf and
 * @c integrationRange, converges to ~1 ulp, and is deterministic, so
 * the quantile readout stays analytic instead of dropping to Monte
 * Carlo.  @p p outside the window's mass clamps to the window edge
 * (the window covers all but a vanishing tail).  NaN when the CDF or
 * the window is unavailable, matching the NaN-as-undecided contract of
 * the other analytic paths.
 */
double numericQuantile(const Distribution &d, double p);
 
/**
 * @name ComparatorRuleRegistry – pairwise §B.2 closed forms
 *
 * Closed-form @f$P(X < Y)@f$ for an ordered pair of independent RV
 * families.  On continuous distributions every ordered comparator reduces
 * to @f$P(X < Y)@f$ or its complement, so rules are keyed on the family
 * pair alone (no operator in the key).  Pairwise behaviour deliberately
 * stays out of the @c Distribution interface: family files self-register
 * their rules at static initialisation through
 * @ref provsql::ComparatorRuleRegistrar, so adding a family touches no existing
 * file, and a missing rule is not an error -- the driver falls back to a
 * family-agnostic quadrature.
 */
///@{
 
/**
 * @brief A pairwise closed form for @f$P(X < Y)@f$, X and Y independent.
 *
 * Returns NaN when its parameter guards fail (e.g. a non-positive rate);
 * the driver then tries the generic quadrature.
 */
using ComparatorRule = double (*)(const Distribution &X,
                                  const Distribution &Y);
 
/**
 * @brief Register the @f$P(X < Y)@f$ closed form for a family pair.
 *
 * Keyed by the families' name tokens rather than descriptor pointers so
 * a family file can register a cross-family rule without depending on
 * another file's static-initialisation order.
 */
void registerComparatorRule(const char *x, const char *y,
                            ComparatorRule rule);
 
/** @brief Static-initialisation helper: one per registered family pair. */
struct ComparatorRuleRegistrar {
  ComparatorRuleRegistrar(const char *x, const char *y,
                          ComparatorRule rule) {
    registerComparatorRule(x, y, rule);
  }
};
 
/**
 * @brief @f$P(X < Y)@f$ for two independent RVs.
 *
 * Applies the registered closed form for the family pair when there is
 * one; on a registry miss (or a rule declining with NaN) falls back to the
 * 1-D composite-Simpson quadrature
 * @f$P(X<Y) = \int (1 - F_Y(t))\, f_X(t)\, dt@f$ over X's integration
 * range.  NaN when neither decides (a density / CDF is undefined, e.g. a
 * non-integer Erlang shape), so the caller falls back to Monte Carlo.
 */
double comparatorPairLess(const Distribution &X, const Distribution &Y);
 
///@}
 
/**
 * @name ClosureRuleRegistry – pairwise family-closure folds on PLUS
 *
 * Closed-form folds of a sum of independent scalar terms into a single
 * distribution (Normal + Normal, same-rate Exponential / Erlang chains).
 * Like the comparator rules, pairwise behaviour stays out of the
 * @c Distribution interface: family files self-register at static
 * initialisation via @ref provsql::ClosureRuleRegistrar, keyed on the (ordered)
 * pair of families that may meet in the sum, and a registry miss simply
 * means the sum stays unfolded (Monte Carlo handles it).
 *
 * A rule receives the whole term list rather than reducing pairwise so
 * its accumulation arithmetic (e.g. one variance sum with a single final
 * square root) is not perturbed by intermediate re-serialisation.
 */
///@{
 
/**
 * @brief One wire of a PLUS under closure: @c a·Z @c + @c b for a base
 *        RV @c Z (@c dist non-null), or a pure additive constant @c b
 *        (@c dist null, @c a @c == @c 0).
 *
 * Structural concerns (base-RV identity, pairwise independence of the
 * @c Z's) are the caller's responsibility; rules only see distributions
 * and coefficients.
 */
struct ClosureTerm {
  const Distribution *dist;
  double a;
  double b;
};
 
/**
 * @brief A family sum-closure fold.  Returns the closed-form
 *        distribution of the summed terms, or @c nullptr when the shape
 *        is outside the closure (mixed rates, scaled / shifted terms a
 *        family cannot absorb, degenerate variance...).
 */
using ClosureRule =
  std::unique_ptr<Distribution> (*)(const std::vector<ClosureTerm> &terms);
 
/**
 * @brief Register the sum-closure rule for a family pair (name-token
 *        keyed, like the comparator rules).
 */
void registerClosureRule(const char *x, const char *y, ClosureRule rule);
 
/** @brief Static-initialisation helper: one per registered family pair. */
struct ClosureRuleRegistrar {
  ClosureRuleRegistrar(const char *x, const char *y, ClosureRule rule) {
    registerClosureRule(x, y, rule);
  }
};
 
/**
 * @brief A family product-closure fold: the distribution of the product
 *        of the (independent) @p factors, or @c nullptr when the shape
 *        is outside the closure.  Scalar factors are the caller's job
 *        (applied afterwards via @c affine).
 */
using ProductRule = std::unique_ptr<Distribution> (*)(
  const std::vector<const Distribution *> &factors);
 
/** @brief Register the product-closure rule for a family pair
 *         (name-token keyed, like the sum-closure rules). */
void registerProductRule(const char *x, const char *y, ProductRule rule);
 
/** @brief Static-initialisation helper: one per registered family pair. */
struct ProductRuleRegistrar {
  ProductRuleRegistrar(const char *x, const char *y, ProductRule rule) {
    registerProductRule(x, y, rule);
  }
};
 
/**
 * @brief Fold a product of independent factors into a single
 *        distribution when a registered closure covers every family
 *        present (same first-vs-each dispatch as @c closePlusTerms);
 *        @c nullptr on any miss.
 */
std::unique_ptr<Distribution> closeProductFactors(
  const std::vector<const Distribution *> &factors);
 
/**
 * @brief A closed-form image of a monotone transform of one family
 *        (e.g. exp of a normal is a lognormal), or @c nullptr when the
 *        family has none.
 */
using TransformRule = std::unique_ptr<Distribution> (*)(
  const Distribution &x);
 
/**
 * @brief Register the image rule for @p transform ("ln" / "exp" / ...;
 *        the evaluator maps its arith opcodes to these names, keeping
 *        this registry opcode-free) applied to @p family.
 */
void registerTransformRule(const char *transform, const char *family,
                           TransformRule rule);
 
/** @brief Static-initialisation helper: one per registered transform. */
struct TransformRuleRegistrar {
  TransformRuleRegistrar(const char *transform, const char *family,
                         TransformRule rule) {
    registerTransformRule(transform, family, rule);
  }
};
 
/**
 * @brief The image distribution of @p transform applied to @p x, when
 *        a registered rule covers @p x's family; @c nullptr otherwise.
 */
std::unique_ptr<Distribution> closeTransform(const char *transform,
                                             const Distribution &x);
 
/**
 * @brief Fold @c PLUS(terms) into a single distribution when a
 *        registered closure covers every family in the sum.
 *
 * Dispatch: the first RV term's family is looked up against itself and
 * against every other RV term's family; all lookups must resolve to the
 * same rule (this is how the Exponential / Erlang pairs share one
 * Erlang-sum rule).  Returns @c nullptr on any miss, on an
 * inconsistent pair, when no RV term is present (the constant-fold
 * path's job), or when the rule itself declines.
 */
std::unique_ptr<Distribution> closePlusTerms(
  const std::vector<ClosureTerm> &terms);
 
///@}
 
/**
 * @name ConjugateRuleRegistry – exact conjugate-prior posterior updates
 *
 * One-observation Bayesian updates for conjugate prior/likelihood pairs:
 * the exact posterior of a latent, wired into one parameter slot of an
 * observed leaf, stays in the prior's family with updated parameters.
 * The conjugate-posterior recogniser (@c ConjugatePosterior.cpp) folds
 * these rules over an @c and_agg evidence conjunction of @c gate_observe
 * atoms, replacing likelihood-weighting importance sampling with the
 * closed form where a rule chain covers every observation.
 *
 * Like the comparator / closure registries, the rules stay out of the
 * @c Distribution interface: the likelihood family's implementation file
 * self-registers its rule(s) at static initialisation, keyed by name
 * tokens (no cross-file static-init order dependence), and a registry
 * miss is a silent decline -- the caller falls back to importance
 * sampling unchanged.
 */
///@{
 
/**
 * @brief The conjugate update for one observation of a given likelihood
 *        family against the running posterior (= prior) family.
 *
 * Both callbacks receive the observed leaf's parsed template @c lik (its
 * non-wired slots carry the resolved literal parameters, e.g. the known
 * @c σ of a @c normal(θ,σ) likelihood) and the observed datum.
 */
struct ConjugateRule {
  /**
   * @brief Update the running posterior parameters @p q1 / @p q2 (a
   *        distribution of the prior's family) by one observation.
   *        Returns @c false to decline (guard failure: an invalid
   *        literal slot, a datum outside the likelihood's support...),
   *        which declines the whole recognition.
   */
  bool (*update)(double &q1, double &q2,
                 const DistributionTemplate &lik, double datum);
  /**
   * @brief Log predictive density @f$\ln m(d \mid q_1, q_2)@f$ of the
   *        datum under the running posterior BEFORE the update; summed
   *        over the fold it is the exact log marginal likelihood behind
   *        @c provsql.evidence.  NaN (or a null pointer) = unavailable,
   *        so @c evidence() keeps its Monte Carlo path.
   */
  double (*log_predictive)(double q1, double q2,
                           const DistributionTemplate &lik, double datum);
};
 
/**
 * @brief Register the conjugate update for observations of
 *        @p likelihood_family whose parameter @p wired_param (0 = p1,
 *        1 = p2) is the latent, against a prior of @p prior_family.
 */
void registerConjugateRule(const char *likelihood_family, int wired_param,
                           const char *prior_family,
                           const ConjugateRule &rule);
 
/**
 * @brief Look up the conjugate rule for (likelihood family, wired
 *        parameter position, prior family); @c nullptr on a miss.
 */
const ConjugateRule *lookupConjugateRule(const std::string &likelihood_family,
                                         int wired_param,
                                         const std::string &prior_family);
 
/** @brief Static-initialisation helper: one per registered pair. */
struct ConjugateRuleRegistrar {
  ConjugateRuleRegistrar(const char *likelihood_family, int wired_param,
                         const char *prior_family, const ConjugateRule &rule) {
    registerConjugateRule(likelihood_family, wired_param, prior_family, rule);
  }
};
 
///@}
 
}  // namespace provsql
 
#endif  // PROVSQL_DISTRIBUTION_H