Random steps to diffusion

export function kernelDiffusivity(
  kernel: StepKernel,
  tau: number,
  symmetryTolerance = 1e-12,
): Readonly<{
  diffusion: ScientificResult;
  drift: ScientificResult;
  centeredDiffusion: ScientificResult;
}> {
  const m = kernelMoments(kernel);
  const bad = (id: string, unit: string, reason: string) =>
    outside(id, unit, "kernelDiffusivity", "symmetric-finite-variance", reason);
  if (
    !Number.isFinite(tau) ||
    tau <= 0 ||
    m.mean.status !== "value" ||
    m.variance.status !== "value" ||
    typeof m.mean.value !== "number" ||
    typeof m.variance.value !== "number"
  ) {
    const reason =
      "A positive step interval and finite step moments are required; the Cauchy variance is not finite.";
    return {
      diffusion: bad("diffusionCoefficient", "m2/s", reason),
      drift: bad("driftVelocity", "m/s", reason),
      centeredDiffusion: bad("centeredDiffusionCoefficient", "m2/s", reason),
    };
  }
  const centered = value(
    "centeredDiffusionCoefficient",
    "m2/s",
    "kernelDiffusivity",
    m.variance.value / (2 * tau),
    m.variance.value > 0,
  );
  return {
    diffusion:
      Math.abs(m.mean.value) <= symmetryTolerance
        ? value(
            "diffusionCoefficient",
            "m2/s",
            "kernelDiffusivity",
            m.variance.value / (2 * tau),
            m.variance.value > 0,
          )
        : bad(
            "diffusionCoefficient",
            "m2/s",
            "A nonzero mean step produces drift. The pure-diffusion argument omits that first-order term.",
          ),
    drift: value("driftVelocity", "m/s", "kernelDiffusivity", m.mean.value / tau),
    centeredDiffusion: centered,
  };
}

export function randomWalkMoments(
  stepRms: number,
  tau: number,
  n: number,
): Computation<{
  mean: number;
  meanSquare: number;
  rms: number;
  elapsedTime: number;
  diffusion: number;
}> {
  if (!isCount(n) || !Number.isFinite(stepRms) || stepRms <= 0 || !Number.isFinite(tau) || tau <= 0)
    return invalid("Use positive finite step RMS and interval, and 0–10000 steps.");
  const meanSquare = n * stepRms * stepRms,
    diffusion = (stepRms * stepRms) / (2 * tau),
    elapsedTime = n * tau;
  if (
    ![meanSquare, diffusion, elapsedTime].every(Number.isFinite) ||
    diffusion <= 0 ||
    (n > 0 && (meanSquare === 0 || elapsedTime === 0))
  )
    return numerical("The moment calculation exceeded binary64 range.");
  return {
    kind: "accepted",
    data: { mean: 0, meanSquare, rms: Math.sqrt(meanSquare), elapsedTime, diffusion },
  };
}

/** Exact integer coefficients through n=64; positive center-out recurrence otherwise.
 * Central normalization avoids cancellation and factorial overflow. Extreme tails
 * beyond binary64 at large n round to zero; no such tails occur in the n<=400 fixtures.
 */
export function coinWalkDistribution(n: number, ell: number): Computation<CoinDistribution> {
  if (!isCount(n) || !Number.isFinite(ell) || ell <= 0 || !Number.isFinite(n * ell))
    return invalid("Use 0–10000 steps and a positive finite step scale.");
  const positions = Float64Array.from({ length: n + 1 }, (_, k) => (2 * k - n) * ell),
    probabilities = new Float64Array(n + 1);
  if (n <= 64) {
    const coefficients: string[] = [];
    const denominator = 1n << BigInt(n);
    let coefficient = 1n;
    for (let k = 0; k <= n; k++) {
      coefficients.push(coefficient.toString());
      probabilities[k] = Number(coefficient) / Number(denominator);
      if (k < n) coefficient = (coefficient * BigInt(n - k)) / BigInt(k + 1);
    }
    return {
      kind: "accepted",
      data: { positions, probabilities, coefficients, denominator: denominator.toString() },
    };
  }
  const mid = Math.floor(n / 2);
  probabilities[mid] = 1;
  for (let k = mid; k < n; k++) {
    const pk = probabilities[k] ?? 0;
    probabilities[k + 1] = (pk * (n - k)) / (k + 1);
  }
  for (let k = mid; k > 0; k--) {
    const pk = probabilities[k] ?? 0;
    probabilities[k - 1] = (pk * k) / (n - k + 1);
  }
  let total = 0,
    correction = 0;
  for (const p of probabilities) {
    const y = p - correction,
      t = total + y;
    correction = t - total - y;
    total = t;
  }
  for (let k = 0; k <= n; k++) {
    const pk = probabilities[k] ?? 0;
    probabilities[k] = pk / total;
  }
  return { kind: "accepted", data: { positions, probabilities } };
}

export async function recordWalks(
  input: WalkSetup,
  selectedStep: number,
  options: WalkExecutionOptions = {},
): Promise<Computation<WalkRecording>> {
  const valid = validate(input);
  if (valid.kind !== "accepted") return valid;
  const p = valid.data;
  if (!Number.isInteger(selectedStep) || selectedStep < 0 || selectedStep > p.runSteps)
    return invalid("Observe an integer step within the declared recording.");
  let config: ReturnType<typeof execution>;
  try {
    config = execution(options);
  } catch {
    return invalid("Use a bounded deterministic chunk size.");
  }
  if (options.cancelled?.()) return cancelled();
  const protectedSteps = new Set([
    0,
    ...WALK_MILESTONES.filter((n) => n <= p.runSteps),
    p.runSteps,
  ]);
  const checkpoints = new Map(
    [...new Set([...protectedSteps, selectedStep])]
      .sort((a, b) => a - b)
      .map((n) => [n, new Float64Array(p.walkers)]),
  );
  const traceValues = new Float64Array(Math.min(WALK_TRACE_COUNT, p.walkers) * (p.runSteps + 1));
  let work = 0;
  for (let walker = 0; walker < p.walkers; walker++) {
    const next = sampler(p, walker);
    let x = 0;
    for (let n = 1; n <= p.runSteps; n++) {
      if (work % config.chunk === 0) {
        if (options.cancelled?.()) return cancelled();
        await config.yieldControl();
        if (options.cancelled?.()) return cancelled();
      }
      x += next();
      work++;
      if (!Number.isFinite(x))
        return failed(
          "A walk coordinate exceeded the representable range; no partial trial was accepted.",
        );
      const checkpoint = checkpoints.get(n);
      if (checkpoint) checkpoint[walker] = x;
      if (walker < WALK_TRACE_COUNT) traceValues[walker * (p.runSteps + 1) + n] = x;
    }
  }
  if (options.cancelled?.()) return cancelled();
  return {
    kind: "accepted",
    data: Object.freeze({
      setup: p,
      checkpoints,
      traceValues,
      protectedSteps,
      allocationId: BM05_ALLOCATION.allocationId,
      normalVersion: HOST_NORMAL_VERSION,
      draws: work * WALK_KERNELS[p.kernel].drawsPerStep,
      retainedBytes:
        traceValues.byteLength +
        [...checkpoints.values()].reduce((sum, a) => sum + a.byteLength, 0),
    }),
  };
}