From wandering to a number.

export function estimateIncrements(
  increments: Float64Array,
  dt: number,
  d: number,
  estimatorId: EstimatorId,
): Assessment<Estimate> {
  if (
    !(increments instanceof Float64Array) ||
    ![1, 2, 3].includes(d) ||
    increments.length % d ||
    increments.length > 10000 ||
    !positive(dt) ||
    ![
      "independent-increment-known-zero-drift",
      "drift-centered",
      "maximum-likelihood-centered",
    ].includes(estimatorId)
  )
    return invalid(
      "Use finite equally timed increments, 1–3 coordinates, and a registered estimator (at most 10000 coordinates).",
    );
  const M = increments.length / d,
    centered = estimatorId !== "independent-increment-known-zero-drift";
  if (M < (centered ? 2 : 1))
    return noValue(
      "underdetermined",
      centered
        ? "A second independent increment, or a known drift, is needed to separate drift from spread."
        : "No admitted increments are available.",
    );
  let scale = 0;
  for (const x of increments) {
    if (!Number.isFinite(x)) return invalid("Every displacement must be finite.");
    scale = Math.max(scale, Math.abs(x));
  }
  const means = new Float64Array(d),
    drift = new Float64Array(d);
  if (scale > 0)
    for (let i = 0; i < M; i++)
      for (let c = 0; c < d; c++) {
        const meanC = means[c] ?? 0;
        const inc = increments[i * d + c] ?? 0;
        means[c] = meanC + (inc / scale - meanC) / (i + 1);
      }
  let sum = 0,
    compensation = 0;
  if (scale > 0)
    for (let i = 0; i < M; i++)
      for (let c = 0; c < d; c++) {
        const inc = increments[i * d + c] ?? 0;
        const meanC = means[c] ?? 0;
        const x = inc / scale - (centered ? meanC : 0);
        const term = x * x - compensation,
          next = sum + term;
        compensation = next - sum - term;
        sum = next;
      }
  const observationInterval = dt,
    q = d * (centered ? M - 1 : M),
    normalization = estimatorId === "maximum-likelihood-centered" ? d * M : q;
  const sumSquares = sum * scale * scale,
    unbiasedDHat = sumSquares / (2 * observationInterval) / q,
    dHat = sumSquares / (2 * observationInterval) / normalization;
  for (let c = 0; c < d; c++) drift[c] = ((means[c] ?? 0) * scale) / dt;
  if (
    ![sumSquares, unbiasedDHat, dHat, ...drift].every(Number.isFinite) ||
    (sum > 0 && (sumSquares === 0 || dHat === 0 || unbiasedDHat === 0))
  )
    return invalid(
      "The displacement scale exceeds this binary64 calculation's representable range.",
    );
  return {
    kind: "accepted",
    data: Object.freeze({
      dHat,
      unbiasedDHat,
      M,
      d,
      q,
      sumSquares,
      normalization,
      biasFactor: q / normalization,
      drift,
      estimatorId,
    }),
  };
}

export function inverseBias(
  q: number,
  normalization = q,
): Assessment<Readonly<{ meanFactor: number; varianceFactor: number | null }>> {
  if (!positive(q) || !positive(normalization))
    return invalid("Positive degrees of freedom and normalization are required.");
  if (q <= 2)
    return noValue(
      "not-applicable",
      "The inverse estimate has no finite expectation with two or fewer degrees of freedom.",
    );
  return {
    kind: "accepted",
    data: {
      meanFactor: normalization / (q - 2),
      varianceFactor: q > 4 ? (2 * normalization ** 2) / ((q - 2) ** 2 * (q - 4)) : null,
    },
  };
}

export function invertToMolecularNumber(
  input: Conditions & { dHat: number; interval: StatisticalInterval; synthetic: boolean },
  set: ConstantSet,
): Assessment<
  Readonly<{
    estimate: number;
    interval: StatisticalInterval;
    constantSetId: string;
    semanticKind: NumberMeaning;
    consistencyRatio: number | null;
    estimatedBoltzmannConstant: number | null;
    consistencyRatioInterval: StatisticalInterval | null;
    estimatedBoltzmannConstantInterval: StatisticalInterval | null;
  }>
> {
  const conditions = checkConditions(input);
  if (conditions.kind !== "accepted") return conditions;
  const g = gas(set, input.synthetic);
  if (g.kind !== "accepted") return g;
  if (input.dHat === 0 || input.interval.lower === 0)
    return noValue(
      "underdetermined",
      "No positive diffusion scale was resolved; a finite molecular number cannot be recovered by division by zero.",
    );
  if (
    !positive(input.dHat) ||
    !positive(input.interval.lower) ||
    !positive(input.interval.upper) ||
    input.interval.lower > input.interval.upper
  )
    return invalid("A positive diffusion estimate and ordered interval are required.");
  // Live-term identifiers for show-the-code: diffusionCoefficient, avogadroNumberEstimate,
  // particleRadius, viscosity, temperature, molarGasConstant. Never a modern exact k_B here.
  const temperature = input.T,
    viscosity = input.eta,
    particleRadius = input.a,
    molarGasConstant = g.data.R,
    diffusionCoefficient = input.dHat;
  const C = (molarGasConstant * temperature) / (6 * Math.PI * viscosity * particleRadius);
  const avogadroNumberEstimate = C / diffusionCoefficient;
  const estimate = avogadroNumberEstimate;
  const lower = C / input.interval.upper,
    upper = C / input.interval.lower;
  const modern = g.data.semanticKind === "consistency-check";
  const avogadro = modern ? constantValue(set, "avogadroConstant").value : 1;
  const consistencyRatio = modern ? estimate / avogadro : null;
  const estimatedBoltzmannConstant = modern
    ? (6 * Math.PI * input.eta * input.a * input.dHat) / input.T
    : null;
  const consistencyRatioInterval = modern
    ? Object.freeze({
        ...input.interval,
        lower: lower / avogadro,
        upper: upper / avogadro,
      })
    : null;
  const estimatedBoltzmannConstantInterval = modern
    ? Object.freeze({
        ...input.interval,
        lower: (6 * Math.PI * input.eta * input.a * input.interval.lower) / input.T,
        upper: (6 * Math.PI * input.eta * input.a * input.interval.upper) / input.T,
      })
    : null;
  const modernMetrics =
    modern &&
    consistencyRatio !== null &&
    estimatedBoltzmannConstant !== null &&
    consistencyRatioInterval !== null &&
    estimatedBoltzmannConstantInterval !== null
      ? [
          consistencyRatio,
          estimatedBoltzmannConstant,
          consistencyRatioInterval.lower,
          consistencyRatioInterval.upper,
          estimatedBoltzmannConstantInterval.lower,
          estimatedBoltzmannConstantInterval.upper,
        ]
      : [];
  if (![estimate, lower, upper, ...modernMetrics].every(positive))
    return invalid("The inverse result cannot be represented at this numerical scale.");
  return {
    kind: "accepted",
    data: Object.freeze({
      estimate,
      interval: Object.freeze({ ...input.interval, lower, upper }),
      constantSetId: set.id,
      semanticKind: g.data.semanticKind,
      consistencyRatio,
      estimatedBoltzmannConstant,
      consistencyRatioInterval,
      estimatedBoltzmannConstantInterval,
    }),
  };
}

/** Conservative union-bound procedure. Marginal input intervals do NOT need
 * statistical independence. Calibration is held exact in this BM-07 slice. */
export function combinedMolecularNumberInterval(
  input: Conditions & {
    estimate: Estimate;
    alphaD: number;
    inputs: Partial<Record<"T" | "eta" | "a", DeclaredInterval>>;
    synthetic: boolean;
  },
  set: ConstantSet,
): Assessment<StatisticalInterval> {
  const conditions = checkConditions(input);
  if (conditions.kind !== "accepted") return conditions;
  const g = gas(set, input.synthetic);
  if (g.kind !== "accepted") return g;
  let alpha = input.alphaD;
  const ranges: Record<"T" | "eta" | "a", readonly [number, number]> = {
    T: [input.T, input.T],
    eta: [input.eta, input.eta],
    a: [input.a, input.a],
  };
  if (Object.keys(input.inputs).some((k) => !["T", "eta", "a"].includes(k)))
    return invalid(
      "Only temperature, viscosity and radius intervals are admitted by this procedure.",
    );
  for (const k of ["T", "eta", "a"] as const) {
    const v = input.inputs[k];
    if (v === undefined) continue;
    if (!v || !Number.isFinite(v.coverage) || v.coverage <= 0 || v.coverage > 1)
      return noValue("not-applicable", `The ${k} input interval has no valid declared coverage.`);
    if (!positive(v.lower) || !positive(v.upper) || v.lower > input[k] || v.upper < input[k])
      return invalid(`The positive ${k} interval must contain its stated point input.`);
    alpha += 1 - v.coverage;
    ranges[k] = [v.lower, v.upper];
  }
  if (!positive(input.alphaD) || input.alphaD >= 1 || !positive(1 - alpha))
    return noValue(
      "not-applicable",
      "The declared input error rates leave no positive simultaneous coverage guarantee.",
    );
  const interval = estimatorInterval(input.estimate, input.alphaD);
  if (interval.kind !== "accepted") return interval;
  if (interval.data.lower === 0)
    return noValue("underdetermined", "A positive diffusion scale is required before inversion.");
  const lower =
    (g.data.R * ranges.T[0]) / (6 * Math.PI * ranges.eta[1] * ranges.a[1] * interval.data.upper);
  const upper =
    (g.data.R * ranges.T[1]) / (6 * Math.PI * ranges.eta[0] * ranges.a[0] * interval.data.lower);
  if (![lower, upper].every(positive))
    return invalid("The combined interval cannot be represented.");
  return {
    kind: "accepted",
    data: Object.freeze({
      ...interval.data,
      lower,
      upper,
      coverage: 1 - alpha,
      coverageKind: "conservative",
      estimatorId: `${input.estimate.estimatorId}; Bonferroni with declared input intervals`,
    }),
  };
}

export function identifiabilityFamily(
  input: {
    D: number;
    T: number;
    eta: number;
    radiusRange: readonly [number, number];
    synthetic: boolean;
  },
  set: ConstantSet,
): Assessment<
  Readonly<{
    product: number;
    radii: Float64Array;
    numbers: Float64Array;
    semanticKind: NumberMeaning;
    constantSetId: string;
  }>
> {
  const g = gas(set, input.synthetic);
  if (g.kind !== "accepted") return g;
  if (input.D === 0)
    return noValue(
      "underdetermined",
      "A positive diffusion scale is needed to draw the compatible family.",
    );
  if (
    ![input.D, input.T, input.eta, ...input.radiusRange].every(positive) ||
    input.radiusRange[0] >= input.radiusRange[1]
  )
    return invalid("Use positive conditions and an increasing radius range.");
  const product = (g.data.R * input.T) / (6 * Math.PI * input.eta * input.D);
  const radii = Float64Array.from({ length: 41 }, (_, i) =>
    Math.exp(
      Math.log(input.radiusRange[0]) * (1 - i / 40) + (Math.log(input.radiusRange[1]) * i) / 40,
    ),
  );
  const numbers = Float64Array.from(radii, (a) => product / a);
  if (![product, ...numbers].every(positive))
    return invalid("The compatible family cannot be represented.");
  return {
    kind: "accepted",
    data: { product, radii, numbers, semanticKind: g.data.semanticKind, constantSetId: set.id },
  };
}