RangeCheck.cpp

Go to the documentation of this file.

/**
 * @file RangeCheck.cpp
 * @brief Implementation of the support-based bound check pass.
 *        See @c RangeCheck.h for the full docstring.
 */
#include "RangeCheck.h"
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <stack>
#include <unordered_map>
#include <unordered_set>
#include <vector>
 
#include "Aggregation.h"        // ComparisonOperator + cmpOpFromOid
#include "AnalyticEvaluator.h"  // cdfAt for shape_mass under truncation
#include "CircuitFromMMap.h"    // getGenericCircuit
#include "ConjugatePosterior.h" // conjugatePosterior (observe-evidence shapes)
#include "Expectation.h"        // lift_conditioning
#include "RandomVariable.h"     // parse_distribution_spec
#include "distributions/Distribution.h"       // makeDistribution -> per-family support()
#include "provsql_utils_cpp.h"  // uuid2string
 
#include <type_traits>          // std::is_same_v in truncateShape
#include <variant>
extern "C" {
#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"            // get_call_result_type, BlessTupleDesc
#include "access/htup_details.h" // heap_form_tuple
#include "utils/uuid.h"
#include "provsql_utils.h"      // gate_type, provsql_arith_op
#include "provsql_error.h"
 
PG_FUNCTION_INFO_V1(rv_support);
}
 
namespace provsql {
 
namespace {
 
/**
 * @brief Closed interval @c [lo, hi] on the extended real line.
 *
 * @c -INFINITY / @c +INFINITY are used for unbounded ends (e.g. the
 * support of a normal RV is @c {-INF, +INF}).  Empty intervals are
 * not generated by any constructor below; comparators against an
 * empty interval would be vacuous and we consider them undecidable.
 */
struct Interval {
  double lo;
  double hi;
 
  static Interval point(double v) { return {v, v}; }
  static Interval all()           { return {-std::numeric_limits<double>::infinity(),
                                            +std::numeric_limits<double>::infinity()}; }
  bool isAll() const {
    return std::isinf(lo) && lo < 0 && std::isinf(hi) && hi > 0;
  }
};
 
Interval add(Interval a, Interval b) { return {a.lo + b.lo, a.hi + b.hi}; }
Interval sub(Interval a, Interval b) { return {a.lo - b.hi, a.hi - b.lo}; }
Interval neg(Interval a)             { return {-a.hi, -a.lo}; }
 
/* Interval product: take the min/max of the four corner products.
 * Handles signed bounds correctly (no special case for negative). */
Interval mul(Interval a, Interval b)
{
  double p1 = a.lo * b.lo, p2 = a.lo * b.hi;
  double p3 = a.hi * b.lo, p4 = a.hi * b.hi;
  return {std::min({p1, p2, p3, p4}), std::max({p1, p2, p3, p4})};
}
 
/* Interval division: if the divisor straddles zero, the result is
 * unbounded in both directions; otherwise compute via @c mul(a, 1/b).
 * The conservative all-real fallback is correct (any real value is
 * possible) but throws away precision &ndash; division by an interval
 * crossing zero is rare in our tests. */
/* exp is monotone increasing and total (exp(-inf) = 0, exp(inf) = inf). */
Interval expInt(Interval a) { return {std::exp(a.lo), std::exp(a.hi)}; }
 
/* ln is monotone increasing on [0, inf); the part of an interval below
 * 0 is a domain violation the sampler raises on, so the bound covers
 * the draws that do evaluate (lo <= 0 maps to -inf via ln 0). */
Interval lnInt(Interval a)
{
  const double lo = a.lo > 0.0 ? std::log(a.lo)
                               : -std::numeric_limits<double>::infinity();
  const double hi = a.hi > 0.0 ? std::log(a.hi)
                               : -std::numeric_limits<double>::infinity();
  return {lo, hi};
}
 
/* x^y over the interval box.  For a base interval entirely >= 0, x^y is
 * monotone in each variable separately (in x for fixed y, in y for
 * fixed x), so the extrema sit at the corners; 0^negative diverges to
 * +inf, which std::pow reports directly.  A base interval extending
 * below 0 keeps the conservative all-real bound: integer-exponent draws
 * are legitimate there, and non-integer ones raise in the sampler. */
Interval powInt(Interval b, Interval e)
{
  if (!(b.lo >= 0.0))
    return Interval::all();
  double lo = std::numeric_limits<double>::infinity();
  double hi = -std::numeric_limits<double>::infinity();
  for (double x : {b.lo, b.hi})
    for (double y : {e.lo, e.hi}) {
      const double v = std::pow(x, y);
      if (std::isnan(v))
        return Interval::all();
      lo = std::min(lo, v);
      hi = std::max(hi, v);
    }
  return {lo, hi};
}
 
Interval divInt(Interval a, Interval b)
{
  if (b.lo <= 0.0 && b.hi >= 0.0)
    return Interval::all();
  Interval inv = {1.0 / b.hi, 1.0 / b.lo};
  return mul(a, inv);
}
 
/**
 * @brief Recursively compute the interval of @p g's value across
 *        worlds.  Memoised in @p cache.
 *
 * Recognised gate types:
 * - @c gate_value: point interval on the parsed scalar.
 * - @c gate_rv:    distribution support (uniform exact, exponential
 *                  on @c [0, +∞), normal on @c (-∞, +∞)).
 * - @c gate_arith: propagated via the interval-arith helpers above.
 *
 * Anything else (e.g. an aggregate gate reached via a HAVING cmp)
 * yields the all-real interval, which downstream conservatively
 * treats as undecidable.
 */
Interval intervalOf(const GenericCircuit &gc, gate_t g,
                    std::unordered_map<gate_t, Interval> &cache)
{
  auto it = cache.find(g);
  if (it != cache.end()) return it->second;
 
  Interval result = Interval::all();
  auto type = gc.getGateType(g);
 
  switch (type) {
    case gate_value:
      /* A value RangeCheck cannot read as a double (e.g. a text constant
       * from an agg_token = text comparison) carries no numeric interval
       * constraint: leave it unconstrained rather than aborting the whole
       * load-time pass.  Mirrors the "undecidable -> all()" default. */
      try {
        result = Interval::point(parseDoubleStrict(gc.getExtra(g)));
      } catch (const CircuitException &) {
        result = Interval::all();
      }
      break;
    case gate_rv: {
      auto spec = parse_distribution_spec(gc.getExtra(g));
      if (!spec) break;
      // Natural support per family (Normal ℝ, Uniform [a,b], Exp/Erlang [0,∞)).
      const DistSupport s = makeDistribution(*spec)->support();
      result = {s.lo, s.hi};
      break;
    }
    case gate_arith: {
      auto op = static_cast<provsql_arith_op>(gc.getInfos(g).first);
      const auto &wires = gc.getWires(g);
      if (wires.empty()) break;
      Interval first = intervalOf(gc, wires[0], cache);
      switch (op) {
        case PROVSQL_ARITH_PLUS:
          result = first;
          for (std::size_t i = 1; i < wires.size(); ++i)
            result = add(result, intervalOf(gc, wires[i], cache));
          break;
        case PROVSQL_ARITH_TIMES:
          result = first;
          for (std::size_t i = 1; i < wires.size(); ++i)
            result = mul(result, intervalOf(gc, wires[i], cache));
          break;
        case PROVSQL_ARITH_MINUS:
          if (wires.size() != 2) break;
          result = sub(first, intervalOf(gc, wires[1], cache));
          break;
        case PROVSQL_ARITH_DIV:
          if (wires.size() != 2) break;
          result = divInt(first, intervalOf(gc, wires[1], cache));
          break;
        case PROVSQL_ARITH_INTDIV:
          /* Truncation toward zero is monotone: truncate the bounds. */
          if (wires.size() != 2) break;
          result = divInt(first, intervalOf(gc, wires[1], cache));
          if (std::isfinite(result.lo)) result.lo = std::trunc(result.lo);
          if (std::isfinite(result.hi)) result.hi = std::trunc(result.hi);
          break;
        case PROVSQL_ARITH_NEG:
          if (wires.size() != 1) break;
          result = neg(first);
          break;
        case PROVSQL_ARITH_MAX:
          /* max/min are monotone in each argument, so the support propagates
           * directly and soundly (independent of any correlation between the
           * children): [max(lo_i), max(hi_i)] resp. [min(lo_i), min(hi_i)]. */
          result = first;
          for (std::size_t i = 1; i < wires.size(); ++i) {
            Interval o = intervalOf(gc, wires[i], cache);
            result = { std::max(result.lo, o.lo), std::max(result.hi, o.hi) };
          }
          break;
        case PROVSQL_ARITH_MIN:
          result = first;
          for (std::size_t i = 1; i < wires.size(); ++i) {
            Interval o = intervalOf(gc, wires[i], cache);
            result = { std::min(result.lo, o.lo), std::min(result.hi, o.hi) };
          }
          break;
        case PROVSQL_ARITH_POW:
          if (wires.size() != 2) break;
          result = powInt(first, intervalOf(gc, wires[1], cache));
          break;
        case PROVSQL_ARITH_LN:
          if (wires.size() != 1) break;
          result = lnInt(first);
          break;
        case PROVSQL_ARITH_EXP:
          if (wires.size() != 1) break;
          result = expInt(first);
          break;
        case PROVSQL_ARITH_ROUND:
          /* Rounding is monotone, so the bounds carry over; the digits, if a
           * second wire gives them, are a constant we do not read here, and
           * widening by one unit in the last place of the bounds is sound for
           * any of them. */
          result = { std::floor(first.lo), std::ceil(first.hi) };
          break;
        case PROVSQL_ARITH_FLOOR:
          if (wires.size() != 1) break;
          result = { std::floor(first.lo), std::floor(first.hi) };
          break;
        case PROVSQL_ARITH_CEIL:
          if (wires.size() != 1) break;
          result = { std::ceil(first.lo), std::ceil(first.hi) };
          break;
        case PROVSQL_ARITH_ASFLOAT8:
          /* The value in its own type: a double carries it as it is, and
           * rounding to a real is monotone, so the interval follows. */
          if (wires.size() != 1) break;
          result = first;
          break;
        case PROVSQL_ARITH_ASFLOAT4:
          if (wires.size() != 1) break;
          result = { static_cast<double>(static_cast<float>(first.lo)),
                     static_cast<double>(static_cast<float>(first.hi)) };
          break;
        case PROVSQL_ARITH_ABS:
          /* |x| is not monotone: an interval straddling zero has 0 as its
           * least value, and the greatest is the farther endpoint. */
          if (wires.size() != 1) break;
          if (first.lo >= 0.0)
            result = first;
          else if (first.hi <= 0.0)
            result = { -first.hi, -first.lo };
          else
            result = { 0.0, std::max(-first.lo, first.hi) };
          break;
        case PROVSQL_ARITH_PERCENTILE:
          /* Continuous percentile over interleaved [ind, x, ...] wires:
           * every draw interpolates within the present values, so the hull
           * of the value wires' supports is a sound bound regardless of
           * which subset is present. */
          if (wires.size() < 2 || wires.size() % 2 != 0) break;
          result = intervalOf(gc, wires[1], cache);
          for (std::size_t i = 3; i < wires.size(); i += 2) {
            Interval o = intervalOf(gc, wires[i], cache);
            result = { std::min(result.lo, o.lo), std::max(result.hi, o.hi) };
          }
          break;
      }
      break;
    }
    case gate_semimod: {
      /* HAVING-style constant wrapper: semimod(gate_one, value).  The
       * semiring action of gate_one (always true) on a scalar leaves
       * the scalar unchanged in every world, so the interval of the
       * semimod equals the interval of its value child.  Other
       * semimod shapes (non-trivial k_gate) keep the conservative
       * all-real default here; @c support_intervalOf widens it to
       * a sampling-time support for rv_* callers. */
      const auto &wires = gc.getWires(g);
      if (wires.size() == 2 && gc.getGateType(wires[0]) == gate_one)
        result = intervalOf(gc, wires[1], cache);
      break;
    }
    case gate_mixture: {
      /* Support of a mixture is the union of its branch supports.
       * Two shapes:
       *  - Classic 3-wire [p_token, x_token, y_token]: the Bernoulli
       *    is a Boolean leaf and contributes nothing to the scalar
       *    interval.
       *  - Categorical N-wire [key, mul_1, ..., mul_n]: each mulinput
       *    carries its outcome value in extra; the support is the
       *    [min, max] of those values. */
      const auto &wires = gc.getWires(g);
      if (gc.isCategoricalMixture(g)) {
        double lo = std::numeric_limits<double>::infinity();
        double hi = -std::numeric_limits<double>::infinity();
        bool any = false;
        for (std::size_t i = 1; i < wires.size(); ++i) {
          double v;
          try { v = parseDoubleStrict(gc.getExtra(wires[i])); }
          catch (const CircuitException &) { any = false; break; }
          lo = std::min(lo, v);
          hi = std::max(hi, v);
          any = true;
        }
        if (any) result = {lo, hi};
      } else if (wires.size() == 3) {
        Interval ix = intervalOf(gc, wires[1], cache);
        Interval iy = intervalOf(gc, wires[2], cache);
        result = {std::min(ix.lo, iy.lo), std::max(ix.hi, iy.hi)};
      }
      break;
    }
    case gate_case: {
      /* Guarded selection [g_1, v_1, ..., g_k, v_k, default]: the result is
       * always one of the value branches, so the support is the union of the
       * values' supports (the guards only decide which one).  Values are the
       * odd-indexed wires plus the final default. */
      const auto &wires = gc.getWires(g);
      if (wires.empty()) break;
      result = intervalOf(gc, wires.back(), cache);  /* default */
      for (std::size_t i = 1; i + 1 < wires.size(); i += 2) {
        Interval o = intervalOf(gc, wires[i], cache);
        result = { std::min(result.lo, o.lo), std::max(result.hi, o.hi) };
      }
      break;
    }
    default:
      /* gate_agg is intentionally not handled here -- the empty-subset
       * NULL semantics make a flat interval misleading, so the
       * runRangeCheck loop dispatches agg-bearing cmps to a separate
       * decider that knows the asymmetry between sound FALSE and
       * unsound TRUE decisions for SUM / MIN / MAX.  All other gate
       * types fall through to the all-real default. */
      break;
  }
 
  cache[g] = result;
  return result;
}
 
/**
 * @brief Decide a @c gate_cmp from the interval of @c (lhs - rhs).
 *
 * Returns @c NaN when the comparator cannot be decided from interval
 * bounds alone (e.g. the difference straddles zero, or the comparator
 * is @c = / @c <> on overlapping continuous supports &ndash; both of
 * which need a CDF, which a downstream analytic pass can supply).
 * Otherwise returns the certain probability @c 0.0 or @c 1.0.
 */
double decideCmp(const Interval &diff, ComparisonOperator op)
{
  switch (op) {
    case ComparisonOperator::LT:
      if (diff.hi < 0.0)  return 1.0;
      if (diff.lo >= 0.0) return 0.0;
      break;
    case ComparisonOperator::LE:
      if (diff.hi <= 0.0) return 1.0;
      if (diff.lo > 0.0)  return 0.0;
      break;
    case ComparisonOperator::GT:
      if (diff.lo > 0.0)  return 1.0;
      if (diff.hi <= 0.0) return 0.0;
      break;
    case ComparisonOperator::GE:
      if (diff.lo >= 0.0) return 1.0;
      if (diff.hi < 0.0)  return 0.0;
      break;
    case ComparisonOperator::EQ:
      /* Disjoint supports ⇒ certainly false.  Overlapping supports of
       * continuous RVs would have probability zero in the measure-
       * theoretic sense, but the interval pass alone cannot tell
       * whether either side is continuous; leave that to a downstream
       * analytic-CDF pass when one is available. */
      if (diff.hi < 0.0 || diff.lo > 0.0) return 0.0;
      break;
    case ComparisonOperator::NE:
      if (diff.hi < 0.0 || diff.lo > 0.0) return 1.0;
      break;
  }
  return std::numeric_limits<double>::quiet_NaN();
}
 
/**
 * @brief Decide a @c gate_cmp where one side is a @c gate_agg, the
 *        other is a scalar constant.
 *
 * Computes a value-interval for the aggregate from its semimod
 * children's per-row values, then folds the comparator like the
 * non-agg path &ndash; but accepts only FALSE decisions, never
 * TRUE.  The reason is structural to ProvSQL's HAVING semantics:
 * the per-aggregator subset enumerators in @c subset.cpp
 * (@c count_enum, @c sum_dp, @c enumerate_exhaustive) all skip
 * the empty subset, matching SQL's "no group, no HAVING" rule.
 * So a HAVING cmp's value is the OR over the @em non-empty subsets
 * where the predicate holds.
 *
 * - When no non-empty subset satisfies the predicate (the bound is
 *   strictly disjoint from the threshold on the right side of the
 *   comparator), the cmp value is exactly @c 0 = @c gate_zero.
 *   FALSE decision: sound.
 * - When every non-empty subset satisfies the predicate, the cmp
 *   value equals "the group is non-empty" &ndash; the OR over the
 *   children's k_gates &ndash; which is a non-constant Boolean
 *   expression, @em not @c gate_one.  Returning TRUE here would
 *   replace the cmp with @c gate_one and over-count probability
 *   mass from the empty world (where the group does not exist),
 *   so TRUE decisions are blocked uniformly across all aggregators.
 *
 * Aggregators we don't bound (@c AVG, @c AND, @c OR, @c CHOOSE,
 * @c ARRAY_AGG, @c NONE) fall through to undecidable.
 *
 * @return @c 0.0 if decided to FALSE, @c NaN otherwise.
 */
double decideAggVsConstCmp(const GenericCircuit &gc, gate_t agg_gate,
                           ComparisonOperator op, double const_val,
                           bool agg_on_lhs,
                           bool *out_always_true = nullptr)
{
  AggregationOperator aop = getAggregationOperator(gc.getInfos(agg_gate).first);
 
  /* Extract per-child scalar values from the semimod children. */
  std::vector<double> values;
  for (gate_t child : gc.getWires(agg_gate)) {
    if (gc.getGateType(child) != gate_semimod)
      return std::numeric_limits<double>::quiet_NaN();
    const auto &sw = gc.getWires(child);
    if (sw.size() != 2)
      return std::numeric_limits<double>::quiet_NaN();
    gate_t value_gate = sw[1];
    if (gc.getGateType(value_gate) != gate_value)
      return std::numeric_limits<double>::quiet_NaN();
    try {
      values.push_back(parseDoubleStrict(gc.getExtra(value_gate)));
    } catch (const CircuitException &) {
      return std::numeric_limits<double>::quiet_NaN();
    }
  }
 
  Interval val_interval = Interval::all();
 
  switch (aop) {
    case AggregationOperator::COUNT:
      val_interval = {0.0, static_cast<double>(values.size())};
      break;
    case AggregationOperator::SUM: {
      double sum_neg = 0.0, sum_pos = 0.0;
      for (double v : values) {
        if (v < 0.0) sum_neg += v;
        else          sum_pos += v;
      }
      val_interval = {std::min(0.0, sum_neg), std::max(0.0, sum_pos)};
      break;
    }
    case AggregationOperator::MIN:
    case AggregationOperator::MAX:
      if (values.empty())
        return std::numeric_limits<double>::quiet_NaN();
      val_interval = {*std::min_element(values.begin(), values.end()),
                      *std::max_element(values.begin(), values.end())};
      break;
    default:
      /* AVG / AND / OR / CHOOSE / ARRAY_AGG / NONE: not decidable
       * with this pass. */
      return std::numeric_limits<double>::quiet_NaN();
  }
 
  Interval lhs = agg_on_lhs ? val_interval : Interval::point(const_val);
  Interval rhs = agg_on_lhs ? Interval::point(const_val) : val_interval;
  Interval diff = sub(lhs, rhs);
  double p = decideCmp(diff, op);
 
  /* Only FALSE decisions are universally sound (see doc comment).
   * gate_one would over-credit the empty subset, which provsql_having
   * deliberately excludes from valid worlds; the safe TRUE rewrite
   * is "the group is non-empty" = OR over the agg's K-gates, and is
   * sound only in absorptive semirings.  The TRUE signal is therefore
   * reported via @p out_always_true rather than the return value;
   * the universal load-time @c runRangeCheck caller ignores it,
   * while the probability-side @c runHavingAlwaysTrueRewriter caller
   * acts on it. */
  if (p == 0.0) return 0.0;
  if (p == 1.0 && out_always_true != nullptr) *out_always_true = true;
  return std::numeric_limits<double>::quiet_NaN();
}
 
/**
 * @brief Try to extract a scalar constant from a cmp's child.
 *
 * Recognises two shapes:
 * - bare @c gate_value: parse its @c extra as a double;
 * - HAVING-style @c gate_semimod with @c k=gate_one and
 *   @c value=gate_value: parse the value's extra.
 *
 * Returns @c NaN on any other shape.
 */
double extractScalarConst(const GenericCircuit &gc, gate_t g)
{
  auto t = gc.getGateType(g);
  if (t == gate_value) {
    try { return parseDoubleStrict(gc.getExtra(g)); }
    catch (const CircuitException &) {
      return std::numeric_limits<double>::quiet_NaN();
    }
  }
  if (t == gate_semimod) {
    const auto &w = gc.getWires(g);
    if (w.size() != 2) return std::numeric_limits<double>::quiet_NaN();
    if (gc.getGateType(w[0]) != gate_one)
      return std::numeric_limits<double>::quiet_NaN();
    if (gc.getGateType(w[1]) != gate_value)
      return std::numeric_limits<double>::quiet_NaN();
    try { return parseDoubleStrict(gc.getExtra(w[1])); }
    catch (const CircuitException &) {
      return std::numeric_limits<double>::quiet_NaN();
    }
  }
  return std::numeric_limits<double>::quiet_NaN();
}
 
/**
 * @brief Flip the sides of a comparison operator.
 *
 * @c (a op b) is equivalent to @c (b flip(op) a).  Used to normalise
 * a cmp so the random-variable side is always on the left.
 */
ComparisonOperator flipCmpOp(ComparisonOperator op)
{
  switch (op) {
    case ComparisonOperator::LT: return ComparisonOperator::GT;
    case ComparisonOperator::LE: return ComparisonOperator::GE;
    case ComparisonOperator::GT: return ComparisonOperator::LT;
    case ComparisonOperator::GE: return ComparisonOperator::LE;
    case ComparisonOperator::EQ: return ComparisonOperator::EQ;
    case ComparisonOperator::NE: return ComparisonOperator::NE;
  }
  return op;
}
 
/**
 * @brief Interpret a @c gate_cmp as a per-RV constraint @c rv op c.
 *
 * Returns @c true and fills @p rv_out, @p op_out, @p const_out when
 * exactly one side of the cmp is a @c gate_rv and the other a
 * @c gate_value with a parseable scalar; @c false otherwise (both
 * sides are RVs, both constants, an @c arith subtree appears, etc.).
 *
 * Strict-vs-non-strict inequalities are preserved as the operator;
 * the caller decides whether to treat the boundary as inclusive
 * (continuous distributions: measure-zero, irrelevant for
 * feasibility verdicts).
 */
bool asRvVsConstCmp(const GenericCircuit &gc, gate_t cmp_gate,
                    gate_t &rv_out, ComparisonOperator &op_out,
                    double &const_out)
{
  bool ok = false;
  ComparisonOperator op = cmpOpFromOid(gc.getInfos(cmp_gate).first, ok);
  if (!ok) return false;
  const auto &wires = gc.getWires(cmp_gate);
  if (wires.size() != 2) return false;
 
  /* Recognise scalar-vs-constant cmps where the scalar side is a
   * bare gate_rv (the original use case for the per-cmp resolution
   * pass) or a gate_mixture (so the conditioning walker can extract
   * intervals on mixture / categorical variables – value-vs-value
   * cmps are folded upstream by RangeCheck before they reach this
   * walker).  Dirac (gate_value) is never the scalar side of a
   * non-trivial cmp at this point; the value-vs-value pair would have
   * been resolved upstream. */
  auto isScalarRv = [](gate_type t) {
    return t == gate_rv || t == gate_mixture;
  };
  auto t0 = gc.getGateType(wires[0]);
  auto t1 = gc.getGateType(wires[1]);
  if (isScalarRv(t0) && t1 == gate_value) {
    try { const_out = parseDoubleStrict(gc.getExtra(wires[1])); }
    catch (const CircuitException &) { return false; }
    rv_out = wires[0];
    op_out = op;
    return true;
  }
  if (t0 == gate_value && isScalarRv(t1)) {
    try { const_out = parseDoubleStrict(gc.getExtra(wires[0])); }
    catch (const CircuitException &) { return false; }
    rv_out = wires[1];
    op_out = flipCmpOp(op);
    return true;
  }
  return false;
}
 
/**
 * @brief Apply a single @c rv-op-constant constraint to a running
 *        interval for the RV.
 *
 * Strict vs non-strict inequalities collapse onto the same closed
 * interval: continuous distributions assign zero mass to the
 * boundary, so the joint-feasibility verdict is unchanged whether
 * we use @c < or @c <=.  @c <> (NE) cannot be represented as a
 * single interval and is left to the per-cmp pass.
 */
Interval intersectRvConstraint(Interval current, ComparisonOperator op,
                               double c)
{
  switch (op) {
    case ComparisonOperator::LT:
    case ComparisonOperator::LE:
      current.hi = std::min(current.hi, c);
      break;
    case ComparisonOperator::GT:
    case ComparisonOperator::GE:
      current.lo = std::max(current.lo, c);
      break;
    case ComparisonOperator::EQ:
      current.lo = std::max(current.lo, c);
      current.hi = std::min(current.hi, c);
      break;
    case ComparisonOperator::NE:
      /* Cannot represent the complement of a point as a single
       * interval; leave the running interval unchanged. */
      break;
  }
  return current;
}
 
bool intervalEmpty(Interval i) { return i.lo > i.hi; }
 
/**
 * @brief Walk an AND-conjunct tree collecting per-RV interval
 *        constraints from its @c gate_cmp leaves.
 *
 * Shared between @c isAndJointlyInfeasible (which checks for an empty
 * intersection) and the public @c collectRvConstraints / conditional
 * @c compute_support paths.  Descends through @c gate_times,
 * collecting every @c gate_cmp interpretable as `rv op const` and
 * intersecting its constraint into a running interval for that RV.
 *
 * @p complete is set to @c true on entry and cleared if the walk
 * encounters any structure other than the AND-friendly set
 * (@c gate_times, @c gate_cmp, @c gate_input, @c gate_one,
 * @c gate_zero) whose footprint *might* constrain an RV
 * (i.e. excluding bare Bernoulli factors).  Callers that need a
 * tight bound (the closed-form moment shortcut) must check it; the
 * support intersection caller can use the result unconditionally
 * because dropping a disjunctive factor only loosens the interval,
 * which is sound for a superset bound on the conditional support.
 *
 * Cmps that do not interpret as `rv op const` (RV vs RV, arith on
 * either side, agg…) are silently ignored; they belong to the
 * conditioning event but don't constrain a single RV's interval.
 */
void walkAndConjunctIntervals(
    const GenericCircuit &gc, gate_t root,
    std::unordered_map<gate_t, Interval> &rv_intervals,
    std::unordered_map<gate_t, Interval> &support_cache,
    bool &complete)
{
  std::unordered_set<gate_t> seen;
  std::stack<gate_t> stk;
  stk.push(root);
  complete = true;
 
  while (!stk.empty()) {
    gate_t g = stk.top(); stk.pop();
    if (!seen.insert(g).second) continue;
 
    auto t = gc.getGateType(g);
    if (t == gate_cmp) {
      gate_t rv = static_cast<gate_t>(0);
      ComparisonOperator op = ComparisonOperator::EQ;
      double c = 0.0;
      if (!asRvVsConstCmp(gc, g, rv, op, c)) {
        /* Cmp shape we don't interpret (RV vs RV, arith involved).
         * Conservatively mark the walk incomplete: this cmp belongs
         * to the event AND could constrain an RV in a way we can't
         * fold into a single interval. */
        complete = false;
        continue;
      }
      auto it = rv_intervals.find(rv);
      Interval current = (it == rv_intervals.end())
                         ? intervalOf(gc, rv, support_cache)
                         : it->second;
      current = intersectRvConstraint(current, op, c);
      rv_intervals[rv] = current;
      continue;  /* never descend into a cmp's operands */
    }
    if (t == gate_times || t == gate_delta || g == root) {
      /* gate_delta wraps a single child as the δ-semiring identity on
       * Booleans, so the AND-conjunct walker is sound to descend
       * through it -- the wrapper carries no constraint of its own.
       * Skipping the descent would mark the walk incomplete and force
       * the moment caller to fall back to MC even when the inner
       * cmps are decidable closed-form. */
      for (gate_t c : gc.getWires(g)) stk.push(c);
      continue;
    }
    if (t == gate_input || t == gate_update || t == gate_one ||
        t == gate_zero) {
      /* Bernoulli leaf / constants: shift P(event), don't truncate
       * any continuous RV.  Skipping is sound and the walk stays
       * complete. */
      continue;
    }
    /* gate_plus (OR), gate_monus (set diff), gate_arith, gate_rv, ...:
     * could affect an RV's conditional distribution in ways that
     * don't reduce to an interval intersection.  Mark the walk
     * incomplete so a moment closed-form caller falls through to MC. */
    complete = false;
  }
}
 
/**
 * @brief Walk @p root's AND-conjunct cmps and decide whether the
 *        conjunction is jointly infeasible by per-RV interval
 *        intersection.
 *
 * For every @c gate_cmp reachable through a chain of @c gate_times
 * starting at @p root, that is interpretable as @c rv-op-constant,
 * intersect the constraint with the running interval for that RV
 * (initialised to the RV's distribution support).  As soon as any
 * RV's interval becomes empty, the AND is infeasible.
 *
 * Descends only through @c gate_times: @c gate_plus is OR (the
 * disjuncts could individually be feasible even when each is a
 * narrow constraint on the RV, so they do not contribute to the
 * conjunction's infeasibility), @c gate_monus is set difference
 * (likewise), and other gate types break the AND chain.
 *
 * Cmps that this pass cannot interpret (RV vs RV, arith on either
 * side, agg…) are simply ignored: skipping them is sound &ndash; we
 * just have fewer constraints, so we never falsely declare
 * infeasibility we cannot prove.
 */
bool isAndJointlyInfeasible(const GenericCircuit &gc, gate_t root)
{
  std::unordered_map<gate_t, Interval> rv_intervals;
  std::unordered_map<gate_t, Interval> support_cache;
  bool complete;
  walkAndConjunctIntervals(gc, root, rv_intervals, support_cache, complete);
  for (const auto &kv : rv_intervals) {
    if (intervalEmpty(kv.second)) return true;
  }
  return false;
}
 
/**
 * @brief Memoised recursive predicate: does @p g's sub-circuit
 *        produce a continuous random variable (no point-mass /
 *        Dirac component)?
 *
 * Used to widen the EQ / NE = 0 / 1 shortcut at the cmp resolution
 * site below the bare-@c gate_rv test, so multi-gate composites like
 * <tt>Exp(0.4) + Exp(0.3) = c</tt> (heterogeneous-rate exponential
 * sum, no closed-form Erlang fold) or
 * <tt>mixture(p, Normal, Uniform) = c</tt> (Bernoulli mixture over
 * two continuous arms) also resolve at load time.  Without this the
 * cmp falls through to AnalyticEvaluator (which returns NaN for
 * EQ / NE) and then to the MC marginalisation, which in finite
 * precision estimates @c P(X = Y) at 0 anyway -- but costs
 * @c provsql.rv_mc_samples iterations to do so.
 *
 * Recursion:
 *  - @c gate_rv -> true (Normal / Uniform / Exp / Erlang all have
 *    continuous densities, no point masses).
 *  - @c gate_value -> false (Dirac at the literal).
 *  - @c gate_arith -> true iff every wire has only-continuous
 *    support.  Sums, products, negations, divisions of continuous
 *    RVs stay continuous in distribution; a @c gate_value sibling
 *    poisons the result (e.g. @c X + 2 is continuous, but
 *    @c X * 0 = 0 has a Dirac at zero -- handled by the existing
 *    constant-fold pre-pass, but defensive here).
 *  - @c gate_mixture, Bernoulli 3-wire <tt>[p, X, Y]</tt> -> true
 *    iff X and Y are both continuous; the Boolean @c p only chooses
 *    an arm, so it does not affect the support type.
 *  - @c gate_mixture, categorical
 *    <tt>[key, mul_1, ..., mul_n]</tt> -> false (point masses at
 *    each mulinput's outcome value).
 *  - Any other gate type -> false (defensive: gate_plus / gate_times
 *    / gate_cmp / gate_agg are not continuous-RV containers).
 *
 * The cache is keyed on @c gate_t and may be shared across multiple
 * cmp gates inside a single @c runRangeCheck invocation.
 */
bool hasOnlyContinuousSupport(const GenericCircuit &gc, gate_t g,
                              std::unordered_map<gate_t, bool> &cache)
{
  auto it = cache.find(g);
  if (it != cache.end()) return it->second;
  /* Memoise pessimistically before recursing so a malformed cyclic
   * sub-circuit (shouldn't happen on well-formed input) returns
   * @c false rather than blowing the stack. */
  cache[g] = false;
 
  bool result = false;
  auto t = gc.getGateType(g);
  switch (t) {
    case gate_rv: {
      /* A discrete family (Poisson, Binomial) has point masses, so a point
       * event X = c carries positive mass and must NOT take the continuous
       * EQ/NE measure-zero shortcut.  isDiscrete() is the authoritative
       * per-family flag (parse_distribution_template handles a latent /
       * parametric leaf too); a malformed spec falls back to continuous. */
      auto tmpl = parse_distribution_template(gc.getExtra(g));
      result = !(tmpl && tmpl->family->factory(0.0, 0.0)->isDiscrete());
      break;
    }
    case gate_value:
      result = false;
      break;
    case gate_arith: {
      result = true;
      for (gate_t w : gc.getWires(g)) {
        if (!hasOnlyContinuousSupport(gc, w, cache)) { result = false; break; }
      }
      break;
    }
    case gate_mixture: {
      if (gc.isCategoricalMixture(g)) { result = false; break; }
      const auto &w = gc.getWires(g);
      if (w.size() != 3) { result = false; break; }
      result = hasOnlyContinuousSupport(gc, w[1], cache)
            && hasOnlyContinuousSupport(gc, w[2], cache);
      break;
    }
    default:
      result = false;
      break;
  }
 
  cache[g] = result;
  return result;
}
 
/**
 * @brief Recursive collection of the @c gate_rv and @c gate_input
 *        leaves reachable from @p g.
 *
 * The result is a sub-circuit's "random-source footprint": two
 * sub-circuits are independent iff their random-source sets are
 * disjoint.  Used to gate the exact-EQ Dirac sum-product below: the
 * factoring @c P(X = Y) = Σ_v @c P(X=v)·P(Y=v) is only valid when
 * @c X and @c Y are independent, otherwise the per-row coupling
 * (e.g. two mixtures sharing a Bernoulli @c p_token) breaks the
 * factoring and the sum-product silently produces the wrong
 * probability.
 *
 * Descent rules: @c gate_arith and @c gate_mixture descend into all
 * children (Bernoulli @c p_token, categorical key, mulinputs all
 * contribute to the random footprint).  @c gate_value is a
 * deterministic literal and contributes no random source.  Other
 * gate types (Boolean / agg / etc.) don't appear under a continuous
 * cmp side in well-formed circuits; defensively, they contribute
 * nothing.
 */
const std::unordered_set<gate_t> &
collectRandomLeaves(const GenericCircuit &gc, gate_t g,
                    std::unordered_map<gate_t, std::unordered_set<gate_t>> &cache)
{
  auto it = cache.find(g);
  if (it != cache.end()) return it->second;
  /* Insert an empty entry first so a recursive call on a cyclic
   * sub-circuit returns early.  std::unordered_map insertion does
   * not invalidate references to existing elements, but it MAY
   * rehash on growth (invalidating ALL references, including the
   * one we're about to capture).  Build the result locally, then
   * write it back in one shot at the end. */
  cache.emplace(g, std::unordered_set<gate_t>{});
 
  std::unordered_set<gate_t> out;
  auto t = gc.getGateType(g);
  if (t == gate_rv || t == gate_input) {
    out.insert(g);
  } else if (t == gate_arith || t == gate_mixture) {
    for (gate_t w : gc.getWires(g)) {
      const auto &child = collectRandomLeaves(gc, w, cache);
      out.insert(child.begin(), child.end());
    }
  }
 
  /* Overwrite the placeholder; locate by find() to avoid a fresh
   * insertion that could rehash and invalidate other iterators in
   * upstream frames. */
  auto fit = cache.find(g);
  fit->second = std::move(out);
  return fit->second;
}
 
using DiracMap    = std::unordered_map<double, double>;
using DiracMapOpt = std::optional<DiracMap>;
 
/**
 * @brief Recursive extraction of @p g's Dirac mass map (value -> mass).
 *
 * Returns @c std::nullopt when the sub-circuit's discrete component
 * is not statically extractable (e.g. an opaque @c gate_arith over
 * mixtures, a Bernoulli mixture whose @c p_token is a compound
 * Boolean, etc.).  When the sub-circuit is purely continuous the
 * map is well-defined but empty (no Diracs, no masses).
 *
 * Used by the exact EQ shortcut below: for independent @c X, @c Y
 * with extractable mass maps @c M_X, @c M_Y:
 * <tt>P(X = Y) = Σ_{v ∈ M_X ∩ M_Y} M_X[v] · M_Y[v]</tt>.  Continuous
 * components contribute zero by measure-zero arguments (Dirac vs
 * continuous and continuous vs continuous), so they need not appear
 * in the sum.
 *
 * Shape rules:
 * - @c gate_value:v: a Dirac at the literal with mass @c 1.
 * - @c gate_rv: continuous in every supported family, empty map.
 * - categorical @c gate_mixture <tt>[key, mul_1, ..., mul_n]</tt>:
 *   sum @c getProb(mul_i) into @c map[parseDouble(extra(mul_i))].
 *   Multiple mulinputs at the same outcome (which the constructor
 *   doesn't produce but is sound to handle) merge masses.
 * - Bernoulli @c gate_mixture <tt>[p_token, X, Y]</tt> with
 *   @c p_token a bare @c gate_input: pull @c π = @c getProb(p_token)
 *   and recurse into X, Y to get @c M_X, @c M_Y; result is
 *   <tt>π·M_X[v] + (1-π)·M_Y[v]</tt> per outcome value.  Compound
 *   Boolean @c p_tokens (whose probability would have to come from
 *   a recursive @c probability_evaluate call) bail.
 * - Anything else: @c std::nullopt.
 */
DiracMapOpt
collectDiracMassMap(const GenericCircuit &gc, gate_t g,
                    std::unordered_map<gate_t, DiracMapOpt> &cache)
{
  auto it = cache.find(g);
  if (it != cache.end()) return it->second;
  /* Pessimistic cycle guard, same reasoning as @c collectRandomLeaves. */
  cache.emplace(g, std::nullopt);
 
  DiracMapOpt result;
  auto t = gc.getGateType(g);
  switch (t) {
    case gate_value: {
      try {
        DiracMap m;
        m[parseDoubleStrict(gc.getExtra(g))] = 1.0;
        result = std::move(m);
      } catch (const CircuitException &) {
        /* unparseable extra: bail */
      }
      break;
    }
    case gate_rv: {
      /* A CONTINUOUS leaf has no point masses (an empty map is exact: the
       * Dirac sum-product then contributes zero, correct by measure zero).
       * A DISCRETE leaf (Poisson, Binomial) DOES have point masses, but
       * they are not statically enumerable (infinite / large support), so
       * decline (nullopt) -- claiming an empty map here would make the
       * sum-product read "no overlap" and wrongly fold X = c to false.
       * isDiscrete() is the authoritative per-family flag. */
      auto tmpl = parse_distribution_template(gc.getExtra(g));
      if (tmpl && tmpl->family->factory(0.0, 0.0)->isDiscrete())
        result = std::nullopt;
      else
        result = DiracMap{};
      break;
    }
    case gate_mixture: {
      const auto &w = gc.getWires(g);
      if (gc.isCategoricalMixture(g)) {
        DiracMap m;
        bool ok = true;
        for (std::size_t i = 1; i < w.size(); ++i) {
          double v;
          try { v = parseDoubleStrict(gc.getExtra(w[i])); }
          catch (const CircuitException &) { ok = false; break; }
          const double p = gc.getProb(w[i]);
          if (!std::isfinite(p) || p < 0.0 || p > 1.0) { ok = false; break; }
          m[v] += p;
        }
        if (ok) result = std::move(m);
      } else if (w.size() == 3
                 && gc.getGateType(w[0]) == gate_input) {
        const double pi = gc.getProb(w[0]);
        if (std::isfinite(pi) && pi >= 0.0 && pi <= 1.0) {
          auto mx = collectDiracMassMap(gc, w[1], cache);
          auto my = collectDiracMassMap(gc, w[2], cache);
          if (mx && my) {
            DiracMap m;
            for (const auto &[v, mass] : *mx) m[v] += pi * mass;
            for (const auto &[v, mass] : *my) m[v] += (1.0 - pi) * mass;
            result = std::move(m);
          }
        }
      }
      break;
    }
    default:
      break;
  }
 
  auto fit = cache.find(g);
  fit->second = result;
  return result;
}
 
}  // namespace
 
unsigned runRangeCheck(GenericCircuit &gc)
{
  std::unordered_map<gate_t, Interval> cache;
  /* Shared across all cmp gates in this @c runRangeCheck invocation.
   * Keyed on gate_t and immutable across cmp iterations because
   * resolving one cmp only changes the cmp's own gate type, not
   * the sub-circuit underneath @c wires[0..1] of other cmps. */
  std::unordered_map<gate_t, bool> continuous_support_cache;
  std::unordered_map<gate_t, DiracMapOpt> dirac_cache;
  std::unordered_map<gate_t, std::unordered_set<gate_t>> leaf_cache;
  unsigned resolved = 0;
 
  /* Snapshot the cmp gate ids before we start mutating: in-place
   * resolution turns a @c gate_cmp into a @c gate_input, but
   * @c getNbGates only grows, never shrinks, so iterating by index
   * over the original count is safe.  We re-check the type at each
   * step to skip already-resolved slots. */
  const auto nb = gc.getNbGates();
  std::vector<gate_t> cmps;
  for (std::size_t i = 0; i < nb; ++i) {
    auto g = static_cast<gate_t>(i);
    if (gc.getGateType(g) == gate_cmp)
      cmps.push_back(g);
  }
 
  for (gate_t c : cmps) {
    if (gc.getGateType(c) != gate_cmp) continue;  /* defensive */
 
    bool ok = false;
    ComparisonOperator op = cmpOpFromOid(gc.getInfos(c).first, ok);
    if (!ok) continue;
 
    const auto &wires = gc.getWires(c);
    if (wires.size() != 2) continue;
 
    /* Identity shortcut: when both sides of the cmp are the same
     * gate (same UUID), the sampler's per-iteration memoisation
     * guarantees both reads return identical values, so the
     * comparator collapses to a constant.  Universal across gate
     * types and semirings; runs first so neither the continuous
     * EQ/NE shortcut nor the interval-based path needs an explicit
     * @c lhs != rhs guard. */
    if (wires[0] == wires[1]) {
      double p = std::numeric_limits<double>::quiet_NaN();
      switch (op) {
        case ComparisonOperator::EQ:
        case ComparisonOperator::LE:
        case ComparisonOperator::GE:
          p = 1.0; break;
        case ComparisonOperator::NE:
        case ComparisonOperator::LT:
        case ComparisonOperator::GT:
          p = 0.0; break;
      }
      gc.resolveCmpToBernoulli(c, p);
      ++resolved;
      continue;
    }
 
    /* Continuous EQ / NE shortcut: P(X = c) = 0 and P(X != c) = 1
     * exactly when at least one side has a continuous distribution
     * (point equality has measure zero under any continuous
     * distribution).  Universal across semirings: the gate_zero /
     * gate_one rewrite is meaningful in every semiring (not just
     * probability), so the resolution belongs here rather than in
     * AnalyticEvaluator.
     *
     * @c hasOnlyContinuousSupport widens the test beyond a bare
     * @c gate_rv leaf: heterogeneous-rate exponential sums, products
     * of independent continuous RVs, and Bernoulli mixtures over
     * two continuous arms all qualify because their distribution
     * has no point-mass component.  Categorical mixtures (point
     * masses at each outcome value) and pure-deterministic
     * @c gate_value sub-circuits do NOT qualify and fall through to
     * the agg / interval / AnalyticEvaluator paths.
     *
     * The @c wires[0] == @c wires[1] case is already handled by the
     * identity shortcut above. */
    if (op == ComparisonOperator::EQ ||
        op == ComparisonOperator::NE) {
      bool lhs_continuous = hasOnlyContinuousSupport(gc, wires[0],
                                                     continuous_support_cache);
      bool rhs_continuous = hasOnlyContinuousSupport(gc, wires[1],
                                                     continuous_support_cache);
      if (lhs_continuous || rhs_continuous) {
        double p = (op == ComparisonOperator::EQ) ? 0.0 : 1.0;
        gc.resolveCmpToBernoulli(c, p);
        ++resolved;
        continue;
      }
 
      /* Exact Dirac sum-product.  When both sides have extractable
       * @c (value -> mass) maps AND the two sub-circuits are
       * independent (random-leaf footprints disjoint), the
       * convolution at zero of @c (X - Y) has support exactly on
       * @c Dirac(X) ∩ Dirac(Y) with mass
       * <tt>M_X(v) · M_Y(v)</tt> per overlapping value; the
       * continuous and continuous-vs-Dirac contributions vanish by
       * measure zero.  This generalises the bare-disjoint case to
       * any pair of statically-known discrete distributions:
       * <tt>P(categorical(a) = categorical(b))</tt> with overlapping
       * outcomes, mixtures with @c as_random branches, etc.
       *
       * The independence test is essential: two mixtures sharing a
       * Bernoulli @c p_token are correlated and the sum-product
       * factoring breaks (the actual @c P(X=Y) cannot be recovered
       * from the marginals alone).  @c collectRandomLeaves'
       * footprint-disjoint check is the gate.
       *
       * When both maps are empty (purely continuous on both sides)
       * the existing branch above already fired, so the sum-product
       * path here only runs for at-least-one-discrete shapes. */
      auto m_l = collectDiracMassMap(gc, wires[0], dirac_cache);
      auto m_r = collectDiracMassMap(gc, wires[1], dirac_cache);
      if (m_l && m_r) {
        const auto &leaves_l = collectRandomLeaves(gc, wires[0], leaf_cache);
        const auto &leaves_r = collectRandomLeaves(gc, wires[1], leaf_cache);
        bool independent = true;
        for (gate_t leaf : leaves_l) {
          if (leaves_r.count(leaf)) { independent = false; break; }
        }
        if (independent) {
          double p_eq = 0.0;
          /* Iterate over the smaller map to keep the sum at
           * O(min(|M_l|, |M_r|)) lookups. */
          const DiracMap *small = (m_l->size() <= m_r->size()) ? &*m_l : &*m_r;
          const DiracMap *large = (m_l->size() <= m_r->size()) ? &*m_r : &*m_l;
          for (const auto &[v, mass] : *small) {
            auto fit = large->find(v);
            if (fit != large->end()) p_eq += mass * fit->second;
          }
          /* Clamp into @c [0, 1] defensively: floating-point summation
           * of masses (each in [0, 1]) might overshoot by an ULP, and
           * @c resolveCmpToBernoulli requires a strict probability. */
          if (p_eq < 0.0) p_eq = 0.0;
          if (p_eq > 1.0) p_eq = 1.0;
          double p = (op == ComparisonOperator::EQ) ? p_eq : 1.0 - p_eq;
          gc.resolveCmpToBernoulli(c, p);
          ++resolved;
          continue;
        }
      }
    }
 
    /* HAVING-style cmp: agg on one side, scalar constant on the
     * other.  Decide via the agg-aware path which is cheaper than
     * intervalOf + decideCmp and which knows the empty-subset NULL
     * semantics for SUM / MIN / MAX (see decideAggVsConstCmp). */
    bool lhs_is_agg = gc.getGateType(wires[0]) == gate_agg;
    bool rhs_is_agg = gc.getGateType(wires[1]) == gate_agg;
    if (lhs_is_agg != rhs_is_agg) {
      gate_t agg_side   = lhs_is_agg ? wires[0] : wires[1];
      gate_t const_side = lhs_is_agg ? wires[1] : wires[0];
      double const_val = extractScalarConst(gc, const_side);
      if (!std::isnan(const_val)) {
        double p = decideAggVsConstCmp(gc, agg_side, op, const_val,
                                       lhs_is_agg);
        if (!std::isnan(p)) {
          gc.resolveCmpToBernoulli(c, p);
          ++resolved;
          continue;
        }
      }
    }
 
    /* Interval-based path for non-agg cmps (RV, gate_arith, value). */
    Interval lhs = intervalOf(gc, wires[0], cache);
    Interval rhs = intervalOf(gc, wires[1], cache);
    /* Skip if both sides are unbounded; @c decideCmp would never
     * return a decision and the work is wasted. */
    if (lhs.isAll() && rhs.isAll()) continue;
 
    Interval diff = sub(lhs, rhs);
    double p = decideCmp(diff, op);
    if (!std::isnan(p)) {
      gc.resolveCmpToBernoulli(c, p);
      ++resolved;
    }
  }
 
  /* Joint-conjunction pass: walk every @c gate_times and check
   * whether its AND-conjunct cmps, viewed together, constrain some
   * shared RV to an empty interval.  Catches the joint-infeasibility
   * case the per-cmp pass above cannot see (each cmp individually
   * leaves a non-empty range, but their intersection is empty).
   *
   * Snapshot the gate_times indices first: @c resolveGateToZero
   * mutates the type, so iterating the live vector while resolving
   * would skip slots.  The post-snapshot type re-check guards against
   * a @c gate_times that the per-cmp pass somehow already collapsed
   * (currently not possible, but cheap insurance for future passes). */
  const auto nb_after = gc.getNbGates();
  std::vector<gate_t> times_gates;
  for (std::size_t i = 0; i < nb_after; ++i) {
    auto g = static_cast<gate_t>(i);
    if (gc.getGateType(g) == gate_times)
      times_gates.push_back(g);
  }
  for (gate_t t : times_gates) {
    if (gc.getGateType(t) != gate_times) continue;  /* defensive */
    if (isAndJointlyInfeasible(gc, t)) {
      gc.resolveGateToZero(t);
      ++resolved;
    }
  }
 
  return resolved;
}
 
/**
 * @brief Probability-side pre-pass: rewrite HAVING-style @c gate_cmp
 *        gates that are provably TRUE on the agg's value-interval
 *        into an OR over the agg's per-row K-gates.
 *
 * Companion to @c runCountCmpEvaluator's Poisson-binomial pre-pass:
 * where that one resolves @c COUNT op C to a closed-form Bernoulli,
 * this one catches the always-true sub-case (e.g. @c COUNT <= K with
 * @c K >= N inputs, or any aggregator whose value-interval entirely
 * satisfies the predicate) and replaces the cmp with @c gate_plus
 * over the agg's K-gates -- the "group is non-empty" indicator.
 *
 * Why a separate pass: @c runRangeCheck deliberately blocks TRUE
 * decisions because @c gate_one is universally unsound for HAVING
 * (it would credit the empty world).  The safe TRUE rewrite
 * "OR of K-gates" requires absorptive @c gate_plus semantics
 * (probability, Boolean, formula, why, which, max-min, max-max), so
 * the pass is restricted to the probability-evaluate path where
 * absorption is guaranteed by the downstream BoolExpr translation.
 *
 * Fires regardless of @c provsql.cmp_probability_evaluation: when
 * the Poisson-binomial path is disabled (developer A/B testing),
 * this lighter shortcut still catches the always-true case and
 * spares the d-DNNF compiler the 2^N-clause DNF that
 * @c provsql_having's @c enumerate_valid_worlds would otherwise emit.
 *
 * Same matching contract as @c decideAggVsConstCmp for the agg side:
 * cmp wires must be {gate_agg, scalar-const-encoded-as-semimod}, the
 * agg's children must all be @c gate_semimod, and the agg kind must
 * be one of COUNT / SUM / MIN / MAX (the only kinds with an
 * interval).  Mismatches leave the cmp untouched.
 *
 * @param gc  Circuit to mutate in place.
 * @return    Number of comparators rewritten to gate_plus.
 */
/**
 * @brief Decide every comparison whose two sides are the very same gate.
 *
 * Two aggregates that read the same rows and contribute the same value on
 * each are ONE gate, the circuit being hash-consed: @c count(v) and
 * @c count(w) over a group both count its rows.  The comparison is then
 * settled by the operator alone, whatever the data -- but settling it TRUE
 * is not @c gate_one.  A tautology over a group still says that the group is
 * there, which is the OR over the k-gates; @c gate_one would credit the world
 * where the group is empty and SQL returns no row at all.  That is the
 * over-credit the doc comment on @c decideAggVsConstCmp describes, reached by
 * another road.
 *
 * Runs at the FRONT of the resolution pipeline: this is a structural fact
 * about the circuit, not a probabilistic one, and the value simplifier that
 * runs early would otherwise fold the comparison away as a constant and lose
 * the group with it.
 */
unsigned runReflexiveCmpRewriter(GenericCircuit &gc)
{
  unsigned resolved = 0;
  const auto nb = gc.getNbGates();
  std::vector<gate_t> cmps;
  for (std::size_t i = 0; i < nb; ++i) {
    auto g = static_cast<gate_t>(i);
    if (gc.getGateType(g) == gate_cmp)
      cmps.push_back(g);
  }
 
  for (gate_t c : cmps) {
    if (gc.getGateType(c) != gate_cmp) continue;  /* defensive */
    bool ok = false;
    ComparisonOperator op = cmpOpFromOid(gc.getInfos(c).first, ok);
    if (!ok) continue;
    const auto &wires = gc.getWires(c);
    if (wires.size() != 2) continue;
    const bool lhs_is_agg = gc.getGateType(wires[0]) == gate_agg;
    /* The very same gate on both sides.  Two aggregates that read the same
   * rows and contribute the same value on each are ONE gate, the circuit
   * being hash-consed: count(v) and count(w) over a group both count its
   * rows, so the comparison is decided by the operator alone, whatever the
   * data.  Deciding it true must not give gate_one, though -- a tautology
   * over a group still says that the group is there, which is the OR over
   * the k-gates.  gate_one would credit the world where the group is empty,
   * where SQL returns no row at all; that is the same over-credit the doc
   * comment on decideAggVsConstCmp describes, reached by another road. */
  if (!(wires[0] == wires[1] && lhs_is_agg)) continue;
  {
    const bool reflexive_true = op == ComparisonOperator::EQ ||
                                op == ComparisonOperator::LE ||
                                op == ComparisonOperator::GE;
    gate_t agg = wires[0];
    std::vector<gate_t> ks;
    bool shape_ok = true;
 
    /* Every contribution a value: the aggregate then has a value on any
     * non-empty group, so the comparison is not the UNKNOWN of a NULL
     * aggregate, which SQL filters out. */
    for (gate_t ch : gc.getWires(agg)) {
      if (gc.getGateType(ch) != gate_semimod) { shape_ok = false; break; }
      const auto &sw = gc.getWires(ch);
      if (sw.size() != 2 || gc.getGateType(sw[1]) != gate_value) {
        shape_ok = false;
        break;
      }
      ks.push_back(sw[0]);
    }
    if (!shape_ok) continue;
 
    if (!reflexive_true) {
      gc.resolveGateToZero(c);  /* x < x, x > x, x <> x: in no world */
      ++resolved;
      continue;
    }
    /* A scalar COUNT has a value over no row too (it is 0), so its
     * tautology holds in every world, the empty one included. */
    if ((gc.getInfos(agg).second & PROVSQL_AGG_SCALAR_FLAG) != 0 &&
        getAggregationOperator(gc.getInfos(agg).first) ==
          AggregationOperator::COUNT) {
      gc.resolveCmpToBernoulli(c, 1.0);
      ++resolved;
      continue;
    }
    if (ks.empty()) continue;
    gc.resolveCmpToPlusOfKGates(c, ks);
    ++resolved;
    continue;
  }
 
  }
  return resolved;
}
 
unsigned runHavingAlwaysTrueRewriter(GenericCircuit &gc)
{
  unsigned resolved = 0;
  const auto nb = gc.getNbGates();
 
  std::vector<gate_t> cmps;
  cmps.reserve(nb / 8);  /* rough guess */
  for (std::size_t i = 0; i < nb; ++i) {
    auto g = static_cast<gate_t>(i);
    if (gc.getGateType(g) == gate_cmp)
      cmps.push_back(g);
  }
 
  for (gate_t c : cmps) {
    if (gc.getGateType(c) != gate_cmp) continue;  /* defensive */
 
    bool ok = false;
    ComparisonOperator op = cmpOpFromOid(gc.getInfos(c).first, ok);
    if (!ok) continue;
 
    const auto &wires = gc.getWires(c);
    if (wires.size() != 2) continue;
 
    bool lhs_is_agg = gc.getGateType(wires[0]) == gate_agg;
    bool rhs_is_agg = gc.getGateType(wires[1]) == gate_agg;
 
    if (lhs_is_agg == rhs_is_agg) continue;  /* both agg or neither */
 
    gate_t agg_side   = lhs_is_agg ? wires[0] : wires[1];
    gate_t const_side = lhs_is_agg ? wires[1] : wires[0];
 
    double const_val = extractScalarConst(gc, const_side);
    if (std::isnan(const_val)) continue;
 
    bool always_true = false;
    double p = decideAggVsConstCmp(gc, agg_side, op, const_val,
                                   lhs_is_agg, &always_true);
    if (!always_true) {
      /* p might be 0.0 (already handled by runRangeCheck at load time
       * if simplify_on_load is on); skip either way. */
      (void)p;
      continue;
    }
 
    /* Scalar aggregation (no GROUP BY): the single result row always exists, so
     * a tautological predicate is gate_one -- probability 1, including the
     * empty-input world -- but only where the aggregate has a value there.
     * count(*) does (it is 0, and count >= 0 holds); sum, min, max and avg are
     * NULL over no row, the comparison is then unknown and the row is filtered
     * out, so the empty world must be excluded exactly as it is for a group.
     * The "group is non-empty" rewrite below does that, and is what the doc
     * comment on decideAggVsConstCmp calls the empty-world over-credit only
     * where the aggregate is defined on the empty input. */
    if ((gc.getInfos(agg_side).second & PROVSQL_AGG_SCALAR_FLAG) != 0 &&
        getAggregationOperator(gc.getInfos(agg_side).first) ==
          AggregationOperator::COUNT) {
      gc.resolveCmpToBernoulli(c, 1.0);
      ++resolved;
      continue;
    }
 
    /* Gather the per-row K-gates from the agg's semimod children. */
    std::vector<gate_t> ks;
    bool shape_ok = true;
    ks.reserve(gc.getWires(agg_side).size());
    for (gate_t ch : gc.getWires(agg_side)) {
      if (gc.getGateType(ch) != gate_semimod) { shape_ok = false; break; }
      const auto &sw = gc.getWires(ch);
      if (sw.size() != 2) { shape_ok = false; break; }
      ks.push_back(sw[0]);  /* K side; M side is sw[1] = gate_value */
    }
    if (!shape_ok || ks.empty()) continue;
 
    gc.resolveCmpToPlusOfKGates(c, ks);
    ++resolved;
  }
 
  return resolved;
}
 
namespace {
 
/* Sampling-time support for aggregation gates.  Unlike @c intervalOf
 * (shared with the @c runRangeCheck cmp dispatcher, which must stay
 * conservative on agg to respect SQL NULL semantics), this widens
 * gate_agg and non-trivial gate_semimod to the actual range of
 * scalar MC samples @c MonteCarloSampler::evalScalar can produce.
 * Used only by @c compute_support, called from rv_support /
 * rv_histogram / rv_moment fallbacks, so the cmp decider remains
 * untouched.
 *
 * Empty-group convention (see test/sql/continuous_aggregation §5):
 * COUNT and SUM yield 0; MIN / MAX / AVG yield NaN.  NaN sits
 * outside any real interval, so callers binning samples drop those
 * worlds automatically; the moment averagers in
 * Expectation::mc_raw_moment also skip them.
 */
Interval aggSupportOf(const GenericCircuit &gc, gate_t root,
                      std::unordered_map<gate_t, Interval> &cache)
{
  const auto type = gc.getGateType(root);
  if (type == gate_semimod) {
    const auto &wires = gc.getWires(root);
    if (wires.size() != 2) return Interval::all();
    Interval vi = intervalOf(gc, wires[1], cache);
    if (gc.getGateType(wires[0]) == gate_one) return vi;
    /* Boolean k child: per-iteration scalar is value · 1_{k fires},
     * so the support is the union of {0} and the value's range. */
    return Interval{std::min(0.0, vi.lo), std::max(0.0, vi.hi)};
  }
  if (type != gate_agg) return intervalOf(gc, root, cache);
 
  const auto &wires = gc.getWires(root);
  AggregationOperator op =
    getAggregationOperator(gc.getInfos(root).first);
 
  std::vector<gate_t> sm_children;
  sm_children.reserve(wires.size());
  for (gate_t c : wires)
    if (gc.getGateType(c) == gate_semimod) sm_children.push_back(c);
 
  auto child_value_iv = [&](gate_t sm) -> Interval {
    const auto &sw = gc.getWires(sm);
    if (sw.size() != 2) return Interval::all();
    return intervalOf(gc, sw[1], cache);
  };
  auto child_always_fires = [&](gate_t sm) -> bool {
    const auto &sw = gc.getWires(sm);
    return sw.size() == 2 && gc.getGateType(sw[0]) == gate_one;
  };
 
  const auto inf = std::numeric_limits<double>::infinity();
  switch (op) {
    case AggregationOperator::COUNT:
      /* [0, n_rows].  Each semimod contributes 0 or 1 to the count. */
      return Interval{0.0, static_cast<double>(sm_children.size())};
    case AggregationOperator::SUM: {
      /* Per row, contribution is value if k fires, else 0; sum the
       * per-row support intervals.  Always-firing rows contribute
       * their value interval verbatim; possibly-firing rows
       * contribute [min(0, v.lo), max(0, v.hi)]. */
      double lo = 0.0, hi = 0.0;
      for (gate_t sm : sm_children) {
        Interval vi = child_value_iv(sm);
        if (vi.isAll()) return Interval::all();
        if (child_always_fires(sm)) {
          lo += vi.lo;
          hi += vi.hi;
        } else {
          lo += std::min(0.0, vi.lo);
          hi += std::max(0.0, vi.hi);
        }
      }
      return Interval{lo, hi};
    }
    case AggregationOperator::MIN:
    case AggregationOperator::MAX: {
      /* MIN / MAX of values where k_i fires.  The actual MIN (or MAX)
       * is some value from one of the firing rows, so the support is
       * the union of the children's value intervals.  Empty-group
       * worlds finalise to NaN, which sits outside the real
       * interval. */
      if (sm_children.empty()) return Interval::all();
      double lo =  inf;
      double hi = -inf;
      for (gate_t sm : sm_children) {
        Interval vi = child_value_iv(sm);
        if (vi.isAll()) return Interval::all();
        lo = std::min(lo, vi.lo);
        hi = std::max(hi, vi.hi);
      }
      if (lo > hi) return Interval::all();
      return Interval{lo, hi};
    }
    default:
      /* AVG: ratio depends on the world's row count.  AND / OR /
       * CHOOSE / ARRAY_AGG: not numeric carriers rv_* surfaces.
       * Keep the conservative all-real default. */
      return Interval::all();
  }
}
 
}  // namespace
 
std::pair<double, double>
compute_support(const GenericCircuit &gc, gate_t root,
                std::optional<gate_t> event_root)
{
  std::unordered_map<gate_t, Interval> cache;
  Interval iv = aggSupportOf(gc, root, cache);
 
  /* Conditional path: intersect with the event's AND-conjunct
   * constraints on @p root.  Walks event_root collecting `rv op c`
   * cmps; non-target constraints are ignored (they affect P(event)
   * but not the truncation of root's distribution).  Even if the
   * walk is "incomplete" (gate_plus / gate_monus / arith encountered)
   * the result is sound: we're computing a SUPERSET bound on the
   * conditional support, and the unconditional support is already a
   * superset, so the intersection of the collected constraints with
   * the unconditional is also a superset. */
  if (event_root.has_value()) {
    std::unordered_map<gate_t, Interval> rv_intervals;
    bool complete;
    walkAndConjunctIntervals(gc, *event_root, rv_intervals, cache, complete);
    auto it = rv_intervals.find(root);
    if (it != rv_intervals.end()) {
      iv.lo = std::max(iv.lo, it->second.lo);
      iv.hi = std::min(iv.hi, it->second.hi);
      /* Defensively clamp to avoid an inverted interval if a buggy
       * walker produced one; should not happen but cheap. */
      if (iv.lo > iv.hi) iv.lo = iv.hi;
    }
  }
 
  return {iv.lo, iv.hi};
}
 
std::optional<std::pair<double, double>>
collectRvConstraints(const GenericCircuit &gc, gate_t event_root,
                     gate_t target_rv)
{
  std::unordered_map<gate_t, Interval> rv_intervals;
  std::unordered_map<gate_t, Interval> support_cache;
  bool complete;
  walkAndConjunctIntervals(gc, event_root, rv_intervals, support_cache,
                           complete);
  if (!complete) return std::nullopt;
  /* If the walk found no cmp constraining target_rv, the conditional
   * support is the unconditional support (the event is independent
   * of target_rv along the recognised structure).  Returning the
   * unconditional interval lets the moment closed-form path
   * short-circuit to the unconditional moment, matching the
   * mathematical truth. */
  auto it = rv_intervals.find(target_rv);
  Interval iv;
  if (it != rv_intervals.end()) {
    iv = it->second;
    /* Intersect with the RV's own support to be safe (event may
     * over-constrain past the support, e.g. `Exp(λ) < -1`). */
    Interval base = intervalOf(gc, target_rv, support_cache);
    iv.lo = std::max(iv.lo, base.lo);
    iv.hi = std::min(iv.hi, base.hi);
    if (iv.lo > iv.hi) iv.lo = iv.hi;
  } else {
    iv = intervalOf(gc, target_rv, support_cache);
  }
  return std::make_pair(iv.lo, iv.hi);
}
 
/**
 * @brief Parse a @c gate_value's @c extra as a finite @c float8.
 *
 * Sibling of @c extract_constant_string in @c having_semantics.cpp but
 * parsing a double, with a const @c GenericCircuit ref (used in the
 * closed-form shape detector path).  Bails on @c NaN / @c ±Infinity so a downstream
 * stem renderer never sees a non-finite x coordinate.
 */
static bool extract_finite_double(const GenericCircuit &gc, gate_t x,
                                  double &out)
{
  if (gc.getGateType(x) != gate_value) return false;
  const std::string &s = gc.getExtra(x);
  if (s.empty()) return false;
  try {
    size_t idx = 0;
    double v = std::stod(s, &idx);
    if (idx != s.size() || !std::isfinite(v)) return false;
    out = v;
    return true;
  } catch (...) {
    return false;
  }
}
 
/** @brief Same parsing applied to a mulinput's outcome label (categorical). */
static bool extract_mulinput_value(const GenericCircuit &gc, gate_t mul,
                                   double &out)
{
  if (gc.getGateType(mul) != gate_mulinput) return false;
  const std::string &s = gc.getExtra(mul);
  if (s.empty()) return false;
  try {
    size_t idx = 0;
    double v = std::stod(s, &idx);
    if (idx != s.size() || !std::isfinite(v)) return false;
    out = v;
    return true;
  } catch (...) {
    return false;
  }
}
 
std::optional<TruncatedSingleRv>
matchTruncatedSingleRv(const GenericCircuit &gc, gate_t root,
                       std::optional<gate_t> event_root)
{
  if (gc.getGateType(root) != gate_rv) return std::nullopt;
  auto spec = parse_distribution_spec(gc.getExtra(root));
  if (!spec) return std::nullopt;
 
  /* Natural support per family.  Normal is unbounded both sides;
   * Uniform sits exactly on its parameters; Exp / Erlang on
   * [0, +inf).  Used both as the unconditional case and as the
   * intersection seed for collectRvConstraints (which already
   * intersects internally, but the bare-natural case still needs
   * a baseline). */
  const DistSupport nat_support = makeDistribution(*spec)->support();
  double nat_lo = nat_support.lo;
  double nat_hi = nat_support.hi;
 
  /* Unconditional path: return natural support, mark untruncated. */
  if (!event_root.has_value()
      || gc.getGateType(*event_root) == gate_one) {
    return TruncatedSingleRv{*spec, nat_lo, nat_hi, /*truncated=*/false};
  }
 
  /* Infeasible event resolved upstream by RangeCheck: the cmp was
   * folded to gate_zero, the conditional distribution is undefined.
   * @c collectRvConstraints would silently fall back to the natural
   * support here (its walker skips gate_zero like gate_one), so we
   * have to detect this explicitly. */
  if (gc.getGateType(*event_root) == gate_zero) return std::nullopt;
 
  auto iv = collectRvConstraints(gc, *event_root, root);
  if (!iv.has_value()) return std::nullopt;
  if (!(iv->first < iv->second)) return std::nullopt;
 
  return TruncatedSingleRv{*spec, iv->first, iv->second, /*truncated=*/true};
}
 
bool eventIsProvablyInfeasible(const GenericCircuit &gc, gate_t root,
                               std::optional<gate_t> event_root)
{
  if (!event_root.has_value()) return false;
  const auto et = gc.getGateType(*event_root);
  if (et == gate_one) return false;
  /* RangeCheck folded the event to false upstream – universal
   * signal, independent of root gate type (a constant scalar
   * value paired with an impossible cmp lands here too). */
  if (et == gate_zero) return true;
  /* Walk the event's AND-conjuncts; an empty intersection with the
   * RV's natural support is the second infeasibility signal that
   * @c matchTruncatedSingleRv collapses into @c std::nullopt.  Only
   * applicable when the root is itself a bare gate_rv that the
   * walker recognises. */
  if (gc.getGateType(root) != gate_rv) return false;
  auto iv = collectRvConstraints(gc, *event_root, root);
  if (!iv.has_value()) return false;
  return !(iv->first < iv->second);
}
 
/**
 * @brief Unconditional probability mass of a shape over the
 *        interval @c [lo, hi].
 *
 * @c TruncatedSingleRv arms supplied here must carry
 * @c truncated == @c false (the unconditional shape); the helper
 * uses the natural support to compute the CDF endpoints, so calling
 * with an already-truncated input would double-truncate.
 *
 * Recursive: a Bernoulli mixture's mass is the Bernoulli-weighted
 * combination of its arms' masses.  Categorical mass is the sum of
 * outcome masses falling in the interval.  Dirac mass is 1 iff the
 * Dirac value sits in the interval, else 0.  Returns @c std::nullopt
 * when a leaf's spec defeats the closed-form CDF (e.g. non-integer
 * Erlang shape – @c cdfAt returns NaN there).
 */
static std::optional<double>
shape_mass(const ClosedFormShape &s, double lo, double hi)
{
  return std::visit([&](const auto &v) -> std::optional<double> {
    using T = std::decay_t<decltype(v)>;
    if constexpr (std::is_same_v<T, TruncatedSingleRv>) {
      const double a = std::max(lo, v.lo);
      const double b = std::min(hi, v.hi);
      if (!(a < b)) return 0.0;
      const double cl = std::isfinite(a) ? cdfAt(v.spec, a) : 0.0;
      const double ch = std::isfinite(b) ? cdfAt(v.spec, b) : 1.0;
      if (std::isnan(cl) || std::isnan(ch)) return std::nullopt;
      return ch - cl;
    } else if constexpr (std::is_same_v<T, DiracShape>) {
      return (v.value >= lo && v.value <= hi) ? 1.0 : 0.0;
    } else if constexpr (std::is_same_v<T, CategoricalShape>) {
      double m = 0.0;
      for (const auto &pr : v.outcomes)
        if (pr.first >= lo && pr.first <= hi) m += pr.second;
      return m;
    } else if constexpr (std::is_same_v<T, BernoulliMixtureShape>) {
      auto L = shape_mass(*v.left, lo, hi);
      auto R = shape_mass(*v.right, lo, hi);
      if (!L || !R) return std::nullopt;
      return v.p * (*L) + (1.0 - v.p) * (*R);
    }
    return std::nullopt;
  }, s);
}
 
/**
 * @brief Conditional shape after truncating the underlying variable
 *        to @c [lo, hi].
 *
 * Bare-RV arm: intersects its natural / current truncation with
 * @c [lo, hi] and marks the result truncated so downstream
 * @c shape_pdf renormalises by the truncated CDF.  Dirac: keep iff
 * value ∈ interval, otherwise nullopt (infeasible).  Categorical:
 * keep outcomes in interval, renormalise masses.  Bernoulli mixture:
 * recursively truncate each arm and reweight the Bernoulli by the
 * ratio of arm masses (the standard
 * @f$ \pi' = \pi Z_L / (\pi Z_L + (1-\pi) Z_R) @f$ update); a
 * fully-eliminated arm degenerates to the surviving one.  Returns
 * @c nullopt when the truncated shape has zero mass (caller can
 * raise infeasibility).
 */
static std::optional<ClosedFormShape>
truncateShape(const ClosedFormShape &s, double lo, double hi)
{
  return std::visit([&](const auto &v) -> std::optional<ClosedFormShape> {
    using T = std::decay_t<decltype(v)>;
    if constexpr (std::is_same_v<T, TruncatedSingleRv>) {
      const double a = std::max(lo, v.lo);
      const double b = std::min(hi, v.hi);
      if (!(a < b)) return std::nullopt;
      return ClosedFormShape{TruncatedSingleRv{v.spec, a, b, /*trunc=*/true}};
    } else if constexpr (std::is_same_v<T, DiracShape>) {
      if (v.value < lo || v.value > hi) return std::nullopt;
      return ClosedFormShape{v};
    } else if constexpr (std::is_same_v<T, CategoricalShape>) {
      CategoricalShape out;
      double total = 0.0;
      for (const auto &pr : v.outcomes) {
        if (pr.first >= lo && pr.first <= hi) {
          out.outcomes.emplace_back(pr.first, pr.second);
          total += pr.second;
        }
      }
      if (out.outcomes.empty() || !(total > 0.0)) return std::nullopt;
      for (auto &pr : out.outcomes) pr.second /= total;
      return ClosedFormShape{std::move(out)};
    } else if constexpr (std::is_same_v<T, BernoulliMixtureShape>) {
      auto mL = shape_mass(*v.left, lo, hi);
      auto mR = shape_mass(*v.right, lo, hi);
      if (!mL || !mR) return std::nullopt;
      const double pL = v.p         * (*mL);
      const double pR = (1.0 - v.p) * (*mR);
      const double Z = pL + pR;
      if (!(Z > 0.0)) return std::nullopt;
      auto Lt = truncateShape(*v.left,  lo, hi);
      auto Rt = truncateShape(*v.right, lo, hi);
      /* Either arm eliminated by the truncation collapses to the
       * surviving arm (its mass was already 0 in shape_mass, so the
       * reweighted p_arm is 1). */
      if (!Lt && !Rt) return std::nullopt;
      if (!Lt) return Rt;
      if (!Rt) return Lt;
      BernoulliMixtureShape m;
      m.p     = pL / Z;
      m.left  = std::make_shared<ClosedFormShape>(std::move(*Lt));
      m.right = std::make_shared<ClosedFormShape>(std::move(*Rt));
      return ClosedFormShape{std::move(m)};
    }
    return std::nullopt;
  }, s);
}
 
std::optional<ClosedFormShape>
matchClosedFormDistribution(const GenericCircuit &gc, gate_t root,
                            std::optional<gate_t> event_root)
{
  /* Test "event is trivial true": either absent, or resolved to
   * gate_one by load-time simplification. */
  const bool event_trivial = !event_root.has_value()
                          || gc.getGateType(*event_root) == gate_one;
 
  /* Bare gate_rv root: delegate to the existing single-RV matcher
   * so the truncation logic (collectRvConstraints) is the single
   * source of truth across the closed-form-shape surface. */
  if (gc.getGateType(root) == gate_rv) {
    /* Conjugate observe-evidence: the posterior is itself a bare
     * distribution of the prior's family, so the shape is the
     * (untruncated) posterior -- exact pdf/CDF for the histogram and
     * curve renderers.  Declines to the truncation matcher on any
     * mismatch (whose AND-conjunct walker treats a gate_observe as an
     * uninterpretable factor and declines in turn). */
    if (!event_trivial)
      if (auto post = conjugatePosterior(gc, root, *event_root)) {
        const DistSupport sup = makeDistribution(*post)->support();
        return ClosedFormShape{
          TruncatedSingleRv{*post, sup.lo, sup.hi, /*truncated=*/false}};
      }
    auto m = matchTruncatedSingleRv(gc, root, event_root);
    if (!m) return std::nullopt;
    return ClosedFormShape{*m};
  }
 
  /* Helper: match the shape unconditionally first, then if the event
   * is non-trivial extract an interval via collectRvConstraints and
   * apply truncateShape.  Used by the Dirac / categorical / mixture
   * branches below so all three honour conditioning through the same
   * pipeline. */
  auto with_optional_truncation =
    [&](std::optional<ClosedFormShape> unc)
      -> std::optional<ClosedFormShape> {
    if (!unc) return std::nullopt;
    if (event_trivial) return unc;
    auto iv = collectRvConstraints(gc, *event_root, root);
    if (!iv.has_value()) return std::nullopt;
    if (!(iv->first < iv->second)) return std::nullopt;
    return truncateShape(*unc, iv->first, iv->second);
  };
 
  /* Dirac point: a gate_value with extra parseable as a finite
   * float8 (the underlying form of as_random(c)).  Conditioning on
   * a constant is normally folded upstream by RangeCheck to
   * gate_one / gate_zero, but a probabilistic event whose footprint
   * doesn't constrain the constant lands here untouched (the cmp
   * walker returns the unconditional support); truncateShape then
   * keeps the Dirac iff its value falls in the recognised interval. */
  if (gc.getGateType(root) == gate_value) {
    double v;
    if (!extract_finite_double(gc, root, v)) return std::nullopt;
    return with_optional_truncation(ClosedFormShape{DiracShape{v}});
  }
 
  /* gate_mixture: either the explicit categorical form
   * (isCategoricalMixture) or the classic Bernoulli triple
   * [p_token, x_token, y_token]. */
  if (gc.getGateType(root) == gate_mixture) {
    const auto &w = gc.getWires(root);
 
    if (gc.isCategoricalMixture(root)) {
      CategoricalShape cs;
      cs.outcomes.reserve(w.size() - 1);
      for (std::size_t i = 1; i < w.size(); ++i) {
        double v;
        if (!extract_mulinput_value(gc, w[i], v)) return std::nullopt;
        double p = gc.getProb(w[i]);
        if (!std::isfinite(p) || p < 0.0 || p > 1.0) return std::nullopt;
        cs.outcomes.emplace_back(v, p);
      }
      if (cs.outcomes.empty()) return std::nullopt;
      return with_optional_truncation(ClosedFormShape{std::move(cs)});
    }
 
    /* Classic Bernoulli mixture: 3 wires, [p_token, x_token, y_token]
     * with p_token a bare gate_input; compound Boolean p bails (the
     * generic path would need a probability-over-Boolean-circuit
     * pre-pass we deliberately do not run here). */
    if (w.size() != 3) return std::nullopt;
    if (gc.getGateType(w[0]) != gate_input) return std::nullopt;
    double p = gc.getProb(w[0]);
    if (!std::isfinite(p) || p < 0.0 || p > 1.0) return std::nullopt;
 
    auto left  = matchClosedFormDistribution(gc, w[1], std::nullopt);
    auto right = matchClosedFormDistribution(gc, w[2], std::nullopt);
    if (!left || !right) return std::nullopt;
 
    BernoulliMixtureShape m;
    m.p     = p;
    m.left  = std::make_shared<ClosedFormShape>(std::move(*left));
    m.right = std::make_shared<ClosedFormShape>(std::move(*right));
    return with_optional_truncation(ClosedFormShape{std::move(m)});
  }
 
  return std::nullopt;
}
 
}  // namespace provsql
 
extern "C" {
 
/**
 * @brief SQL: rv_support(token uuid, prov uuid, OUT lo float8, OUT hi float8)
 *
 * Loads the persisted circuit rooted at @p token, intersects with the
 * AND-conjunct cmps in @p prov constraining @p token, and returns the
 * resulting @c [lo, hi] support interval.  When @p prov resolves to
 * @c gate_one (the unconditional default after load-time
 * simplification), the conditional path is skipped and the bare
 * unconditional support of @p token is returned.
 *
 * @c -Infinity / @c +Infinity float8 represent unbounded ends (e.g.
 * the support of a normal RV is @c [-Infinity, +Infinity]).
 */
Datum rv_support(PG_FUNCTION_ARGS)
{
  try {
    pg_uuid_t *token = PG_GETARG_UUID_P(0);
    pg_uuid_t *prov  = PG_GETARG_UUID_P(1);
 
    gate_t root_gate, event_gate;
    auto gc = getJointCircuit(*token, *prov, root_gate, event_gate);
 
    /* gate_one as event-side means the conditioning is the trivial
     * "always true" event (either the user passed gate_one() directly
     * or load-time simplification collapsed the event to it).  Take
     * the unconditional path. */
    std::optional<gate_t> event_opt;
    if (gc.getGateType(event_gate) != gate_one)
      event_opt = event_gate;
 
    /* A stored "X | C" arrives here as a conditioned root: peel it to the
     * bare target and fold the condition into the event, so the support is
     * the conditional (truncated) one rather than the unconditional. */
    root_gate = provsql::lift_conditioning(gc, root_gate, event_opt);
 
    auto iv = provsql::compute_support(gc, root_gate, event_opt);
 
    TupleDesc tupdesc;
    Datum values[2];
    bool nulls[2] = {false, false};
 
    if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
      provsql_error("rv_support: expected composite return type");
    tupdesc = BlessTupleDesc(tupdesc);
 
    values[0] = Float8GetDatum(iv.first);
    values[1] = Float8GetDatum(iv.second);
 
    PG_RETURN_DATUM(HeapTupleGetDatum(heap_form_tuple(tupdesc, values, nulls)));
  } catch (const std::exception &e) {
    provsql_error("rv_support: %s", e.what());
  } catch (...) {
    provsql_error("rv_support: unknown exception");
  }
  PG_RETURN_NULL();
}
 
}  // extern "C"