\set ECHO none
\pset format unaligned

-- Case Study 6: The City Air-Quality Sensor Network
-- Tests the random_variable surface end-to-end: gate_rv leaves for the
-- four distribution families, the planner-hook lifting of WHERE
-- comparators on RV columns (gate_cmp), Bernoulli mixtures via
-- provsql.mixture, RV aggregates (sum / avg) lowered through
-- rv_aggregate_semimod, expected() applied to a grouped RV via the
-- subquery + outer-WHERE idiom, UNION ALL on RV columns, and an
-- independent / tree-decomposition / monte-carlo cross-check.  Backs
-- doc/source/user/casestudy6.rst.

-- Pin the MC RNG and sample budget so cases that fall back to MC
-- (probabilities, conditional moments) are reproducible at a tolerance
-- loose enough to absorb sampler noise but tight enough to be
-- meaningful.
SET provsql.monte_carlo_seed = 42;
SET provsql.rv_mc_samples    = 50000;

-- ---------------------------------------------------------------------
-- Setup: schema and seed, mirroring doc/casestudy6/setup.sql.
-- ---------------------------------------------------------------------
CREATE TABLE stations (
  id        text PRIMARY KEY,
  name      text NOT NULL,
  district  text NOT NULL
);
INSERT INTO stations VALUES
  ('s1', 'City Centre',        'centre'),
  ('s2', 'Riverside Park',     'centre'),
  ('s3', 'Industrial Estate',  'east'),
  ('s4', 'Suburban Reference', 'east');
SELECT add_provenance('stations');

CREATE TABLE calibration_status (
  station_id  text PRIMARY KEY REFERENCES stations(id),
  p           double precision NOT NULL
);
INSERT INTO calibration_status VALUES
  ('s1', 0.95),
  ('s2', 0.70),
  ('s3', 0.60),
  ('s4', 1.00);
SELECT add_provenance('calibration_status');

CREATE TABLE readings (
  id          integer PRIMARY KEY,
  station_id  text NOT NULL REFERENCES stations(id),
  ts          timestamp NOT NULL,
  pm25        random_variable  -- NULL: sensor offline, no reading
);
INSERT INTO readings (id, station_id, ts, pm25) VALUES
  (1, 's1', '2026-05-12 08:00', provsql.normal(28.0, 2.0)),
  (2, 's2', '2026-05-12 08:00', provsql.uniform(10.0, 22.0)),
  (3, 's3', '2026-05-12 08:00', provsql.exponential(0.04)),
  (4, 's4', '2026-05-12 08:00', 15.0),
  (5, 's1', '2026-05-12 09:00', provsql.normal(40.0, 4.0)),
  (6, 's2', '2026-05-12 09:00', provsql.uniform(12.0, 24.0)),
  (7, 's3', '2026-05-12 09:00', provsql.erlang(3, 0.1)),
  (8, 's4', '2026-05-12 09:00', 16.5);
SELECT add_provenance('readings');

CREATE TABLE historical_readings (
  id          integer PRIMARY KEY,
  station_id  text NOT NULL REFERENCES stations(id),
  ts          timestamp NOT NULL,
  pm25        random_variable NOT NULL
);
INSERT INTO historical_readings (id, station_id, ts, pm25) VALUES
  (1, 's1', '2026-05-11 08:00', provsql.normal(34.0, 2.5)),
  (2, 's2', '2026-05-11 08:00', provsql.uniform(15.0, 28.0)),
  (3, 's3', '2026-05-11 08:00', provsql.exponential(0.03)),
  (4, 's4', '2026-05-11 08:00', 18.0),
  (5, 's1', '2026-05-11 09:00', provsql.normal(42.0, 3.0)),
  (6, 's2', '2026-05-11 09:00', provsql.uniform(20.0, 35.0)),
  (7, 's3', '2026-05-11 09:00', provsql.erlang(3, 0.08)),
  (8, 's4', '2026-05-11 09:00', 19.5);
SELECT add_provenance('historical_readings');

-- ---------------------------------------------------------------------
-- Step 1: per-row pm25 gates.  Rows 1 and 5 are Normal; row 3 is
-- Exponential; row 4 is the deterministic reference station lifted via
-- the implicit numeric -> random_variable cast (gate_value); row 7 is
-- Erlang.  One representative pick per shape.  The materialise +
-- remove_provenance dance hides the row's own auto-added provsql uuid
-- column from the output (it's non-deterministic).
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_pm25_kind AS
  SELECT id, get_gate_type(pm25::uuid) AS pm25_kind
    FROM readings
   WHERE id IN (1, 3, 4, 7);
SELECT remove_provenance('result_cs6_pm25_kind');
SELECT id, pm25_kind FROM result_cs6_pm25_kind ORDER BY id;
DROP TABLE result_cs6_pm25_kind;

-- ---------------------------------------------------------------------
-- Step 2: WHERE pm25 > 35 lifts the comparator into a gate_cmp
-- conjoined with each row's input via gate_times; the procedure body
-- of random_variable_gt is never executed, so every row survives the
-- WHERE.  The pm25 column kind on rows 5 and 7 confirms that the lift
-- did not strip the underlying distribution.
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_thresh AS
  SELECT id, station_id,
         get_gate_type(provenance())  AS prov_kind,
         get_gate_type(pm25::uuid)    AS pm25_kind
    FROM readings
   WHERE pm25 > 35;
SELECT remove_provenance('result_cs6_thresh');
SELECT id, station_id, prov_kind, pm25_kind
  FROM result_cs6_thresh ORDER BY id;
DROP TABLE result_cs6_thresh;

-- Probability that the event pm25 > 35 fires on each tracked row.  For
-- the s1 readings (Normal):
--   row 1, N(28, 2):  P(X > 35) = 1 - Phi( 3.5) ~= 2.3e-4
--   row 5, N(40, 4):  P(X > 35) = 1 - Phi(-1.25) ~= 0.8944
-- The analytic 'independent' path is exact for a single gate_cmp under
-- a deterministic input gate; MC at seed=42 must sit within 0.02 of it.
CREATE TABLE result_cs6_thresh_prob AS
  SELECT id,
         abs(probability_evaluate(provenance(), 'monte-carlo', '50000')
             - probability_evaluate(provenance(), 'independent')) < 0.02
           AS mc_matches_independent
    FROM readings
   WHERE pm25 > 35
     AND station_id = 's1';
SELECT remove_provenance('result_cs6_thresh_prob');
SELECT id, mc_matches_independent
  FROM result_cs6_thresh_prob ORDER BY id;
DROP TABLE result_cs6_thresh_prob;

-- ---------------------------------------------------------------------
-- Step 3: provsql.mixture(p, x, y) returns a gate_mixture; the
-- pm25 / 1.2 arm (the calibration correction) is a gate_arith(DIV)
-- over the per-row pm25 and a gate_value.  The constructor result
-- type is random_variable, so ::uuid grabs the gate root.
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_mixture AS
  SELECT r.id,
         get_gate_type(
           provsql.mixture(cs.p, r.pm25, r.pm25 / 1.2)::uuid
         ) AS mixture_kind,
         get_gate_type(
           (r.pm25 / 1.2)::uuid
         ) AS scaled_kind
    FROM readings r JOIN calibration_status cs USING (station_id)
   WHERE r.station_id = 's1';
SELECT remove_provenance('result_cs6_mixture');
SELECT id, mixture_kind, scaled_kind
  FROM result_cs6_mixture ORDER BY id;
DROP TABLE result_cs6_mixture;

-- ---------------------------------------------------------------------
-- Step 4: sum(pm25) and avg(pm25) over a tracked RV column lower
-- through rv_aggregate_semimod into gate_arith roots.  Aggregate
-- dispatch on the random_variable result type picks the provsql.sum /
-- provsql.avg overloads from plain sum() / avg(), so the case study's
-- SQL uses the bare names.
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_agg AS
  SELECT s.district,
         sum(r.pm25) AS total_pm25,
         avg(r.pm25) AS avg_pm25
    FROM readings r JOIN stations s ON s.id = r.station_id
   GROUP BY s.district;
SELECT remove_provenance('result_cs6_agg');

-- Structure: sum -> gate_arith(PLUS); avg -> gate_arith(DIV) over two
-- gate_arith(PLUS) subtrees (numerator / denominator).
SELECT district,
       get_gate_type(total_pm25::uuid) AS sum_root,
       get_gate_type(avg_pm25::uuid)   AS avg_root
  FROM result_cs6_agg ORDER BY district;

-- The 'centre' district has rows 1, 2, 5, 6 (N(28,2), U(10,22),
-- N(40,4), U(12,24)) with E[sum] = 28 + 16 + 40 + 18 = 102.
-- The 'east' district has rows 3, 4, 7, 8 (Exp(0.04), 15, Erl(3,0.1),
-- 16.5) with E[sum] = 25 + 15 + 30 + 16.5 = 86.5.  Linearity of
-- expectation hits exactly under the analytical evaluator.
SELECT district,
       abs(provsql.expected(total_pm25) - CASE district
             WHEN 'centre' THEN 102.0
             WHEN 'east'   THEN  86.5
                                              END) < 1e-9
         AS total_mean_exact
  FROM result_cs6_agg ORDER BY district;

DROP TABLE result_cs6_agg;

-- ---------------------------------------------------------------------
-- Step 8: UNION across the two batches dedups by (station_id, id) and
-- combines the two row provenances via gate_plus.  Both batches share
-- the same id space (1..8), so the keys line up directly.  WHERE pm25
-- > 35 lifts the comparator into each contributing row's provenance,
-- so the resulting gate_plus wraps two gate_times(input, gate_cmp)
-- subtrees.
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_union AS
SELECT station_id, id, get_gate_type(provenance()) AS prov_kind
FROM (
    (SELECT station_id, id FROM readings            WHERE pm25 > 35)
  UNION
    (SELECT station_id, id FROM historical_readings WHERE pm25 > 35)
) t;
SELECT remove_provenance('result_cs6_union');
SELECT station_id, id, prov_kind
  FROM result_cs6_union
  ORDER BY station_id, id;
DROP TABLE result_cs6_union;

-- ---------------------------------------------------------------------
-- Step 9: filter the per-district aggregates by their expected average.
-- The HAVING qual is deterministic from the planner-hook's perspective
-- (no agg_token Var, no provenance_aggregate wrapper -- avg(rv)
-- returns random_variable and expected() collapses it to a double), so
-- needs_having_lift leaves it for PostgreSQL to evaluate natively while
-- the per-group provenance still gets a delta(gate_agg) wrapper.
-- E[avg(pm25)] is 102/4 = 25.5 for centre and 86.5/4 = 21.625 for
-- east, both above 20.
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_having AS
  SELECT s.district, get_gate_type(provenance()) AS prov_kind
    FROM readings r JOIN stations s ON s.id = r.station_id
   GROUP BY s.district
  HAVING expected(avg(r.pm25)) > 20;
SELECT remove_provenance('result_cs6_having');
SELECT district, prov_kind FROM result_cs6_having ORDER BY district;
DROP TABLE result_cs6_having;

-- ---------------------------------------------------------------------
-- Step 10: independent vs monte-carlo on the threshold query.  For
-- each surviving row the closed-form 'independent' result is exact;
-- MC at n=50000 must match within 0.02.  tree-decomposition is also
-- exact and must agree with 'independent' to full precision.
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_methods AS
  SELECT id,
         abs(probability_evaluate(provenance(), 'independent')
             - probability_evaluate(provenance(), 'tree-decomposition')) < 1e-9
           AS ind_equals_td,
         abs(probability_evaluate(provenance(), 'independent')
             - probability_evaluate(provenance(), 'monte-carlo', '50000')) < 0.02
           AS ind_close_to_mc
    FROM readings WHERE pm25 > 35;
SELECT remove_provenance('result_cs6_methods');
SELECT id, ind_equals_td, ind_close_to_mc
  FROM result_cs6_methods ORDER BY id;
DROP TABLE result_cs6_methods;

-- ---------------------------------------------------------------------
-- Step 11: order statistics.  max(pm25) lowers to a gate_arith MAX root
-- over per-row mixtures (absent-row identity -Infinity); the district
-- means are MC at seed 42 (mixture children, no closed form): centre
-- ~ 40.04 (dominated by its N(40,4)), east ~ 39.63 (heavy Exp / Erlang
-- right tails).  The same-row greatest() over two independent bare
-- leaves is analytic (layer-cake quadrature), proven by running it
-- under rv_mc_samples = 0.
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_worst AS
  SELECT s.district,
         get_gate_type(max(r.pm25)::uuid) AS worst_root,
         expected(max(r.pm25))            AS worst_mean
    FROM readings r JOIN stations s ON s.id = r.station_id
   GROUP BY s.district;
SELECT remove_provenance('result_cs6_worst');
SELECT district, worst_root,
       worst_mean BETWEEN 39 AND 41 AS worst_mean_plausible
  FROM result_cs6_worst ORDER BY district;
DROP TABLE result_cs6_worst;

SET provsql.rv_mc_samples = 0;
CREATE TABLE result_cs6_greatest AS
  SELECT abs(expected(provsql.greatest(a.pm25, b.pm25)) - 28.00003) < 1e-3
           AS pairwise_max_analytic
    FROM readings a, readings b
   WHERE a.id = 1 AND b.id = 2;
SELECT remove_provenance('result_cs6_greatest');
SELECT pairwise_max_analytic FROM result_cs6_greatest;
DROP TABLE result_cs6_greatest;

-- ---------------------------------------------------------------------
-- Step 12: quantiles, all closed-form (per-family inverse CDF; Erlang
-- bisected on its regularised-incomplete-gamma CDF), still under
-- rv_mc_samples = 0.  p95 references: N(28,2) -> 28 + 2 z_95, Exp(0.04)
-- -> -ln(0.05)/0.04, N(40,4) -> 40 + 4 z_95, Erl(3,0.1) ~ 62.9579.
-- The conditional median under pm25 > 35 is the truncated-normal
-- inverse CDF at the rescaled u = F(35) + 0.5 (1 - F(35)).
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_p95 AS
  SELECT id, quantile(pm25, 0.95) AS p95
    FROM readings WHERE id IN (1, 3, 5, 7);
SELECT remove_provenance('result_cs6_p95');
SELECT id, abs(p95 - CASE id WHEN 1 THEN 31.289707253902947
                             WHEN 3 THEN 74.89330683884975
                             WHEN 5 THEN 46.57941450780589
                             WHEN 7 THEN 62.957936218719865 END) < 1e-6
         AS p95_exact
  FROM result_cs6_p95 ORDER BY id;
DROP TABLE result_cs6_p95;

CREATE TABLE result_cs6_cond_median AS
  SELECT id, quantile(pm25, 0.5, provenance()) AS median_given_high
    FROM readings
   WHERE pm25 > 35 AND station_id = 's1';
SELECT remove_provenance('result_cs6_cond_median');
SELECT id, abs(median_given_high - CASE id WHEN 1 THEN 35.36132427785059
                                           WHEN 5 THEN 40.531206716705924 END) < 1e-6
         AS cond_median_exact
  FROM result_cs6_cond_median ORDER BY id;
DROP TABLE result_cs6_cond_median;
RESET provsql.rv_mc_samples;
SET provsql.rv_mc_samples = 50000;

-- ---------------------------------------------------------------------
-- Step 13: expected excess over the threshold via CASE (gate_case,
-- evaluated per-draw by MC; guard and branch share the pm25 leaf).
-- Closed-form partial expectation for N(40,4): sigma phi(alpha) +
-- (mu-35) Phi(-alpha) ~ 5.2024; N(28,2) essentially never crosses.
-- The ELSE branch needs an explicit lift (as_random(0)): CASE branch
-- types resolve within one type category.
-- ---------------------------------------------------------------------
CREATE TABLE result_cs6_excess AS
  SELECT id,
         expected(CASE WHEN pm25 > 35 THEN pm25 - 35
                       ELSE as_random(0) END) AS expected_excess
    FROM readings
   WHERE station_id = 's1';
SELECT remove_provenance('result_cs6_excess');
SELECT id,
       CASE id WHEN 1 THEN expected_excess BETWEEN 0 AND 0.01
               WHEN 5 THEN abs(expected_excess - 5.2024) < 0.1 END
         AS excess_close_to_closed_form
  FROM result_cs6_excess ORDER BY id;
DROP TABLE result_cs6_excess;

-- ---------------------------------------------------------------------
-- Step 14: covariance / correlation / stddev.  Disjoint footprints give
-- an exact 0 covariance and exact stddev (no sampling: proven under
-- rv_mc_samples = 0); a shared plume leaf makes the pair correlated,
-- theory Cov = Var(dust) = 9 and rho = 9/sqrt(13*21) ~ 0.545, estimated
-- by MC at seed 42 (E[xy] over a shared leaf has no closed form).
-- ---------------------------------------------------------------------
SET provsql.rv_mc_samples = 0;
CREATE TABLE result_cs6_indep AS
  SELECT covariance(a.pm25, b.pm25) = 0 AS cov_indep_exact,
         stddev(a.pm25)            = 2 AS sd_exact
    FROM readings a, readings b
   WHERE a.id = 1 AND b.id = 2;
SELECT remove_provenance('result_cs6_indep');
SELECT cov_indep_exact, sd_exact FROM result_cs6_indep;
DROP TABLE result_cs6_indep;
RESET provsql.rv_mc_samples;
SET provsql.rv_mc_samples = 100000;

-- Bands are deliberately loose: the MC estimate converges to cov=9, corr~0.545
-- on every platform, but the exact draw stream differs across C++ stdlibs
-- (libstdc++ vs libc++), so a tight band is not portable. This still rejects a
-- wrong implementation (independent plume -> cov~0, corr~0).
CREATE TABLE result_cs6_plume AS
  WITH plume AS (SELECT provsql.normal(0, 3) AS dust)
  SELECT abs(covariance(a.pm25 + p.dust, b.pm25 + p.dust) - 9) < 3
           AS cov_shared_close,
         abs(correlation(a.pm25 + p.dust, b.pm25 + p.dust) - 0.545) < 0.2
           AS corr_shared_close
    FROM readings a, readings b, plume p
   WHERE a.id = 1 AND b.id = 2;
SELECT remove_provenance('result_cs6_plume');
SELECT cov_shared_close, corr_shared_close FROM result_cs6_plume;
DROP TABLE result_cs6_plume;

-- ---------------------------------------------------------------------
-- Step 15: Gamma / Weibull / Pareto, all analytic (rv_mc_samples = 0
-- throughout).  Gamma moments exact; chi-squared quantile bisected;
-- Weibull min-stability gives the two-unit first failure in closed
-- form; Pareto's alpha = 1 mean is honestly Infinity; the mixed-family
-- Pareto-vs-Normal comparison routes through the generic quadrature.
-- ---------------------------------------------------------------------
SET provsql.rv_mc_samples = 0;
SELECT expected(provsql.gamma(2.5, 0.5)) = 5   AS gamma_mean_exact,
       variance(provsql.gamma(2.5, 0.5)) = 10  AS gamma_var_exact,
       abs(quantile(provsql.chi_squared(4), 0.95) - 9.48772903678115) < 1e-6
         AS chi2_crit_exact;

SELECT abs(expected(provsql.least(provsql.weibull(1.5, 4.0),
                                  provsql.weibull(1.5, 4.0)))
           - 2.2747755945647903) < 1e-9         AS weibull_first_failure_exact,
       abs(expected(provsql.weibull(1.5, 4.0))
           - 3.6109811718037346) < 1e-9         AS weibull_lifetime_exact;

SELECT expected(provsql.pareto(30, 2.5)) = 50   AS pareto_mean_exact,
       abs(quantile(provsql.pareto(30, 2.5), 0.95) - 99.43362052019958) < 1e-6
         AS pareto_p95_exact,
       expected(provsql.pareto(30, 1.0)) = 'Infinity'::double precision
         AS pareto_divergent_mean_honest;

SELECT abs(probability(x > y) - 0.38396729360204274) < 1e-6
         AS mixed_family_quadrature_exact
  FROM (SELECT provsql.pareto(30, 2.5) AS x,
               provsql.normal(45, 5)   AS y) t;
RESET provsql.rv_mc_samples;
SET provsql.rv_mc_samples = 50000;

-- ---------------------------------------------------------------------
-- Step 16: SQL-standard statistic aggregates over RV rows under sensor
-- dropout.  Row presence is pinned to the station's calibration
-- probability; the certain reference pair (s4, values 15.0 / 16.5, p=1)
-- gives exact Dirac statistics (stddev_pop = 0.75, corr with a
-- rescaling = 1), and the district median's member set responds to
-- dropout (banded MC assertion).
-- ---------------------------------------------------------------------
DO $$ BEGIN
  PERFORM set_prob(r.provsql, cs.p)
    FROM readings r JOIN calibration_status cs USING (station_id);
END $$;

CREATE TABLE result_cs6_statagg AS
  SELECT expected(stddev_pop(pm25))     AS sd_ref,
         expected(corr(pm25, pm25 * 2)) AS corr_ref
    FROM readings
   WHERE station_id = 's4';
SELECT remove_provenance('result_cs6_statagg');
SELECT sd_ref = 0.75 AS sd_ref_exact, corr_ref = 1 AS corr_ref_exact
  FROM result_cs6_statagg;
DROP TABLE result_cs6_statagg;

CREATE TABLE result_cs6_median AS
  SELECT s.district,
         expected(percentile_cont(0.5)
                  WITHIN GROUP (ORDER BY r.pm25)) AS median_pm25
    FROM readings r JOIN stations s ON s.id = r.station_id
   GROUP BY s.district;
SELECT remove_provenance('result_cs6_median');
-- Bands are loose (cross-stdlib draw streams); they still reject a
-- median that ignores dropout or the member set.
SELECT district,
       CASE district
         WHEN 'centre' THEN median_pm25 BETWEEN 20 AND 32
         WHEN 'east'   THEN median_pm25 BETWEEN 12 AND 25
       END AS median_in_band
  FROM result_cs6_median ORDER BY district;
DROP TABLE result_cs6_median;

-- ---------------------------------------------------------------------
-- Step 17: maintenance-triage headline via CASE over aggregates
-- (agg_case).  Guards and branches are aggregates over the tracked
-- calibration log (min / sum branches: exact possible-worlds machinery;
-- an avg branch would take the Monte-Carlo path, its exact joint
-- (sum, count) arm being unconditional-only);
-- the display cell carries the actual-world CASE value, and the
-- expectation is exact (rv_mc_samples = 0 throughout).
--   One record per district is certain, so the group always exists and
--   every world is enumerable by hand:
--   centre: rows p-col {.95,.70}, probs {1.0,.8}:
--     {both}.8 -> min .70 <= .8 -> min .70; {s1}.2 -> min .95 > .8 ->
--     sum .95;  E = .8*.7 + .2*.95 = .75
--   east: rows p-col {.60,1.00}, probs {.7,1.0}:
--     {both}.7 -> min .6 -> .6; {s4}.3 -> min 1.0 > .8 -> sum 1.0
--     E = .7*.6 + .3*1 = .72
-- ---------------------------------------------------------------------
DO $$ BEGIN
  PERFORM set_prob(provenance(), 1.0) FROM stations;
  PERFORM set_prob(provenance(),
                   CASE station_id WHEN 's1' THEN 1.0
                                   WHEN 's2' THEN 0.8
                                   WHEN 's3' THEN 0.7
                                   ELSE 1.0 END)
    FROM calibration_status;
END $$;

SET provsql.rv_mc_samples = 0;
CREATE TABLE result_cs6_headline AS
  SELECT s.district,
         CASE WHEN min(cs.p) > 0.8 THEN sum(cs.p)
              ELSE min(cs.p) END AS confidence_headline
    FROM calibration_status cs JOIN stations s ON s.id = cs.station_id
   GROUP BY s.district;
SET provsql.active = off;
SELECT district,
       get_gate_type(confidence_headline::uuid) AS headline_gate,
       confidence_headline AS display
  FROM result_cs6_headline ORDER BY district;
SELECT district,
       CASE district
         WHEN 'centre' THEN abs(expected(confidence_headline) - 0.75) < 1e-9
         WHEN 'east'   THEN abs(expected(confidence_headline) - 0.72) < 1e-9
       END AS expected_headline_exact
  FROM result_cs6_headline ORDER BY district;
SET provsql.active = on;
DROP TABLE result_cs6_headline;
RESET provsql.rv_mc_samples;
SET provsql.rv_mc_samples = 50000;

-- ---------------------------------------------------------------------
-- Step 18: mutual information.  Disjoint stochastic-leaf footprints
-- give an exact 0 (no sampling: proven under rv_mc_samples = 0); the
-- shared-plume pair is estimated by the 2-D histogram plug-in over
-- coupled joint draws (Gaussian-pair theory ~ -ln(1-rho^2)/2 ~ 0.176
-- at rho ~ .545; banded loosely for cross-stdlib portability).
-- ---------------------------------------------------------------------
SET provsql.rv_mc_samples = 0;
CREATE TABLE result_cs6_mi_indep AS
  SELECT mutual_information(a.pm25, b.pm25) = 0 AS mi_indep_exact
    FROM readings a, readings b
   WHERE a.id = 1 AND b.id = 2;
SELECT remove_provenance('result_cs6_mi_indep');
SELECT mi_indep_exact FROM result_cs6_mi_indep;
DROP TABLE result_cs6_mi_indep;
RESET provsql.rv_mc_samples;
SET provsql.rv_mc_samples = 100000;

CREATE TABLE result_cs6_mi_shared AS
  WITH plume AS (SELECT provsql.normal(0, 3) AS dust)
  SELECT mutual_information(a.pm25 + p.dust, b.pm25 + p.dust) AS mi_shared
    FROM readings a, readings b, plume p
   WHERE a.id = 1 AND b.id = 2;
SELECT remove_provenance('result_cs6_mi_shared');
SELECT mi_shared BETWEEN 0.05 AND 0.5 AS mi_shared_in_band
  FROM result_cs6_mi_shared;
DROP TABLE result_cs6_mi_shared;

-- ---------------------------------------------------------------------
-- Step 19: an offline sensor reports NULL.  A comparison over the NULL
-- reading is unknown in every world (probability exactly 0, a gate_zero
-- conjunct -- never silently certain); IS NULL is an ordinary
-- deterministic predicate; and aggregates skip the NULL reading in both
-- their sum and their count.
-- ---------------------------------------------------------------------
INSERT INTO readings (id, station_id, ts, pm25)
  VALUES (9, 's2', '2026-05-12 10:00', NULL);

CREATE TABLE result_cs6_null_thresh AS
  SELECT id, ROUND(probability_evaluate(provenance())::numeric, 4) AS p
    FROM readings
   WHERE pm25 > 35 AND ts = '2026-05-12 10:00';
SELECT remove_provenance('result_cs6_null_thresh');
SELECT id, p FROM result_cs6_null_thresh;
DROP TABLE result_cs6_null_thresh;

CREATE TABLE result_cs6_null_isnull AS
  SELECT id, station_id FROM readings WHERE pm25 IS NULL;
SELECT remove_provenance('result_cs6_null_isnull');
SELECT id, station_id FROM result_cs6_null_isnull;
DROP TABLE result_cs6_null_isnull;

-- avg skips the NULL reading: same expected value with or without row 9.
CREATE TABLE result_cs6_null_avg AS
  SELECT round(expected(avg(pm25))::numeric, 4) AS with_null
    FROM readings WHERE station_id = 's2';
SELECT remove_provenance('result_cs6_null_avg');
CREATE TABLE result_cs6_null_avg2 AS
  SELECT round(expected(avg(pm25))::numeric, 4) AS without_null
    FROM readings WHERE station_id = 's2' AND id <> 9;
SELECT remove_provenance('result_cs6_null_avg2');
SELECT a.with_null = b.without_null AS avg_skips_null
  FROM result_cs6_null_avg a, result_cs6_null_avg2 b;
DROP TABLE result_cs6_null_avg;
DROP TABLE result_cs6_null_avg2;

DELETE FROM readings WHERE id = 9;

-- ---------------------------------------------------------------------
-- Cleanup.
-- ---------------------------------------------------------------------
DROP TABLE historical_readings;
DROP TABLE readings;
DROP TABLE calibration_status;
DROP TABLE stations;

RESET provsql.monte_carlo_seed;
RESET provsql.rv_mc_samples;