Keep the data. Change what you know.

/** Inverse questions for the existing Brownian owners. No random draws occur here.
 * Camera relations: Berglund (2010), doi:10.1103/PhysRevE.82.011917.
 * These are later measurement models, not Einstein's printed derivation.
 */
import { makeRefusal } from "../../../experiments/results/refusals.ts";
import type { ConstantSet } from "../constants.ts";
import {
  type Assessment,
  combinedMolecularNumberInterval,
  type DeclaredInterval,
  type Estimate,
  estimatorInterval,
  identifiabilityFamily,
  type StatisticalInterval,
} from "../inference.ts";
import { cameraMoments } from "./observation.ts";

const invalid = (requirements: string): Assessment<never> => ({
  kind: "refused",
  refusal: makeRefusal(
    "invalid-parameter",
    { capabilityId: "diffusion.inference" },
    { details: { requirements } },
  ),
});
const positive = (value: number) => Number.isFinite(value) && value > 0;

/** The same confidence interval for D maps to an entire band, not a unique N. */
export function diffusionCompatibleBand(
  input: {
    estimate: Estimate;
    T: number;
    eta: number;
    alpha: number;
    radiusRange: readonly [number, number];
    synthetic: boolean;
  },
  constants: ConstantSet,
) {
  const interval = estimatorInterval(input.estimate, input.alpha);
  if (interval.kind !== "accepted") return interval;
  const conditions = {
    T: input.T,
    eta: input.eta,
    radiusRange: input.radiusRange,
    synthetic: input.synthetic,
  };
  const center = identifiabilityFamily({ ...conditions, D: input.estimate.dHat }, constants);
  if (center.kind !== "accepted") return center;
  const lower = identifiabilityFamily({ ...conditions, D: interval.data.upper }, constants);
  if (lower.kind !== "accepted") return lower;
  const upper = identifiabilityFamily({ ...conditions, D: interval.data.lower }, constants);
  if (upper.kind !== "accepted") return upper;
  return {
    kind: "accepted" as const,
    data: Object.freeze({
      ...center.data,
      interval: interval.data,
      lowerNumbers: lower.data.numbers,
      upperNumbers: upper.data.numbers,
    }),
  };
}

export type IndependentMeasurement = Readonly<{
  value: number;
  interval: DeclaredInterval;
  provenance: string;
}>;
/** Bonferroni allocation across D and every uncertain auxiliary input.
 * The reader declares the coverage and provenance; we cannot certify a calibration.
 * Dependence between marginal intervals does not invalidate the union bound.
 * R and timing/calibration remain exact within the declared model.
 */
export function conditionalMolecularNumberInterval(
  input: {
    estimate: Estimate;
    T: number;
    eta: number;
    alpha: number;
    synthetic: boolean;
    radius: IndependentMeasurement | null;
    radiusProvenance: "independently-declared" | "same-displacements";
    temperature?: IndependentMeasurement;
    viscosity?: IndependentMeasurement;
  },
  constants: ConstantSet,
): Assessment<StatisticalInterval> {
  if (!positive(input.alpha) || input.alpha >= 1)
    return invalid("Choose a total error probability strictly between zero and one.");
  if (!input.radius)
    return {
      kind: "no-value",
      status: "underdetermined",
      reason:
        "Diffusion fixes a radius–number product. Add an independently established radius interval before asking for a molecular-number interval.",
    };
  const measurements = [input.radius, input.temperature, input.viscosity].filter(
    (value): value is IndependentMeasurement => value !== undefined,
  );
  const alphaD = input.alpha / (measurements.length + 1);
  for (const value of measurements) {
    if (
      !value ||
      typeof value.provenance !== "string" ||
      !value.provenance.trim() ||
      value.provenance.length > 512
    )
      return invalid(
        "Every added measurement needs a provenance description of at most 512 characters.",
      );
    const interval = value.interval;
    if (
      !positive(value.value) ||
      !interval ||
      !positive(interval.lower) ||
      !positive(interval.upper) ||
      interval.lower > value.value ||
      interval.upper < value.value ||
      !Number.isFinite(interval.coverage) ||
      interval.coverage > 1 ||
      interval.coverage <= 0
    )
      return invalid(
        "Positive ordered measurement bounds must contain their point value and have a declared coverage.",
      );
    if (interval.coverage < 1 - alphaD)
      return {
        kind: "no-value",
        status: "not-applicable",
        reason: `Each auxiliary interval needs coverage of at least ${100 * (1 - alphaD)} percent for this Bonferroni allocation. A lower-coverage interval cannot be relabeled.`,
      };
  }
  return combinedMolecularNumberInterval(
    {
      estimate: input.estimate,
      alphaD,
      synthetic: input.synthetic,
      T: input.temperature?.value ?? input.T,
      eta: input.viscosity?.value ?? input.eta,
      a: input.radius.value,
      radiusProvenance: input.radiusProvenance,
      inputs: {
        a: input.radius.interval,
        ...(input.temperature ? { T: input.temperature.interval } : {}),
        ...(input.viscosity ? { eta: input.viscosity.interval } : {}),
      },
    },
    constants,
  );
}

/** Inverting one camera variance gives a segment of compatible (D, sigma²) pairs.
 * The sampled points exclude D=0, which is separately reported as a boundary.
 * Their predicted covariance is evaluated by cameraMoments, not by a view.
 */
export function cameraCompatibleFamily(input: {
  variance: number;
  dt: number;
  exposure: number;
}): Assessment<
  Readonly<{
    maximumDiffusion: number;
    maximumNoiseVariance: number;
    diffusion: Float64Array;
    noiseVariance: Float64Array;
    covariance: Float64Array;
  }>
> {
  const { variance, dt, exposure } = input;
  if (
    !positive(variance) ||
    !positive(dt) ||
    !Number.isFinite(exposure) ||
    exposure < 0 ||
    exposure > dt
  )
    return invalid(
      "Use a positive coordinate variance and frame spacing, with known uniform exposure between zero and that spacing.",
    );
  const effectiveTime = dt - exposure / 3;
  const maximumDiffusion = variance / 2 / effectiveTime;
  const maximumNoiseVariance = variance / 2;
  if (!positive(maximumDiffusion) || !positive(maximumNoiseVariance))
    return invalid("The compatible segment cannot be represented at this numerical scale.");
  const diffusion = new Float64Array(41),
    noiseVariance = new Float64Array(41),
    covariance = new Float64Array(41);
  for (let i = 0; i < diffusion.length; i++) {
    const fraction = (i + 1) / diffusion.length;
    // Parameterization avoids a negative endpoint caused by subtracting rounded products.
    const D = maximumDiffusion * fraction,
      sigma2 = maximumNoiseVariance * (1 - fraction);
    const predicted = cameraMoments({ D, dt, exposure, sigma: Math.sqrt(sigma2), drift: 0, d: 1 });
    if (predicted.kind !== "accepted") return predicted;
    diffusion[i] = D;
    noiseVariance[i] = sigma2;
    covariance[i] = predicted.data.covariance;
  }
  return {
    kind: "accepted",
    data: Object.freeze({
      maximumDiffusion,
      maximumNoiseVariance,
      diffusion,
      noiseVariance,
      covariance,
    }),
  };
}

/** Two distinct frame spacings, SAME exposure duration and localization variance.
 * These are unconstrained moment estimates, not confidence intervals. Negative
 * estimates are retained as model/sample diagnostics, never clipped to zero.
 */
export function twoIntervalCameraEstimate(input: {
  firstVariance: number;
  secondVariance: number;
  firstDt: number;
  secondDt: number;
  exposure: number;
}): Assessment<Readonly<{ D: number; sigma2: number; physical: boolean }>> {
  const { firstVariance, secondVariance, firstDt, secondDt, exposure } = input;
  if (
    ![firstVariance, secondVariance, firstDt, secondDt].every(positive) ||
    !Number.isFinite(exposure) ||
    exposure < 0 ||
    exposure > Math.min(firstDt, secondDt) ||
    secondDt <= firstDt
  )
    return invalid(
      "Use positive moments at two increasing spacings and the same known exposure duration, no longer than either spacing.",
    );
  const D = (secondVariance - firstVariance) / 2 / (secondDt - firstDt);
  const sigma2 = firstVariance / 2 - D * (firstDt - exposure / 3);
  if (![D, sigma2].every(Number.isFinite))
    return invalid("The two-interval inverse is not representable.");
  return { kind: "accepted", data: Object.freeze({ D, sigma2, physical: D >= 0 && sigma2 >= 0 }) };
}

import { makeRefusal } from "../../experiments/results/refusals.ts";
import type { ScientificResult } from "../../experiments/results/types.ts";
import type { OutputContract, Parameters } from "../../experiments/store/instanceStore.ts";
import {
  cameraCompatibleFamily,
  conditionalMolecularNumberInterval,
  diffusionCompatibleBand,
  twoIntervalCameraEstimate,
} from "../../physics/reference/inference/identifiability.ts";
import { covarianceEstimator } from "../../physics/reference/inference/observation.ts";
import { INFERENCE_CONSTANTS } from "../../physics/reference/inference/synthetic.ts";
import {
  type Assessment,
  estimateIncrements,
  estimatorInterval,
  invertToMolecularNumber,
} from "../../physics/reference/inference.ts";
import type { InferenceEvidence } from "./evidence.ts";

const contract = (
  unit: string,
  semanticKind: string,
  statuses: OutputContract["statuses"] = ["value"],
): OutputContract =>
  Object.freeze({ unit, semanticKind, ownerId: "inference.workbench", statuses });
const inferred = ["value", "underdetermined", "not-applicable", "outside-domain"] as const;
export const FAMILY_OUTPUTS = Object.freeze({
  diffusionCoefficientEstimate: contract("m2/s", "statistical-diffusion-estimate"),
  diffusionInterval: contract("m2/s", "chi-square-confidence-interval"),
  familyRadii: contract("m", "compatible-radius-grid"),
  familyNumbers: contract("1/mol", "synthetic-recovery"),
  familyLowerNumbers: contract("1/mol", "synthetic-recovery"),
  familyUpperNumbers: contract("1/mol", "synthetic-recovery"),
  radiusNumberProduct: contract("m/mol", "synthetic-recovery"),
  avogadroNumberEstimate: contract("1/mol", "synthetic-recovery", inferred),
  molecularInterval: contract("1/mol", "synthetic-recovery", inferred),
  simultaneousCoverage: contract("1", "conservative-coverage-lower-bound", inferred),
});
export const CAMERA_OUTPUTS = Object.freeze({
  sampleVariance: contract("m2", "known-drift-increment-second-moment"),
  sampleCovariance: contract("m2", "known-drift-adjacent-increment-product", inferred),
  secondVariance: contract("m2", "coarse-increment-second-moment", inferred),
  familyDiffusion: contract("m2/s", "compatible-camera-diffusivity"),
  familyNoise: contract("m2", "compatible-camera-localization-variance"),
  familyCovariance: contract("m2", "compatible-camera-covariance"),
  maximumDiffusion: contract("m2/s", "single-moment-compatible-bound"),
  maximumNoise: contract("m2", "single-moment-compatible-bound"),
  diffusionEstimate: contract("m2/s", "unconstrained-camera-moment-estimate", inferred),
  noiseEstimate: contract("m2", "unconstrained-localization-moment-estimate", inferred),
  physicalSolution: contract("1", "nonnegative-moment-solution", inferred),
  diffusionConfidenceInterval: contract("m2/s", "unavailable-camera-confidence-interval", inferred),
});
export const RADIUS_EXAMPLE = Object.freeze({
  radiusKnown: true,
  radius: 0.5e-6,
  radiusLower: 0.45e-6,
  radiusUpper: 0.55e-6,
  radiusCoverage: 0.975,
  radiusProvenance: "independently-declared",
  provenance:
    "Synthetic independent-radius interval for the worked example; not an actual calibration.",
});
export const RADIUS_DEFAULTS = Object.freeze({
  ...RADIUS_EXAMPLE,
  radiusKnown: false,
  provenance: "",
});
export const CAMERA_DEFAULTS = Object.freeze({ information: "variance" });
export type RadiusParameters = { [K in keyof typeof RADIUS_DEFAULTS]: (typeof RADIUS_DEFAULTS)[K] };
function outputsFor(contracts: Readonly<Record<string, OutputContract>>) {
  function meta(id: string) {
    const c = contracts[id];
    if (!c) throw new Error(`Unknown inference output: ${id}.`);
    return { quantityId: id, unit: c.unit, semanticKind: c.semanticKind, ownerId: c.ownerId };
  }
  return {
    value(id: string, value: number | Float64Array): ScientificResult {
      return { ...meta(id), status: "value", value };
    },
    absent(
      id: string,
      status: "underdetermined" | "not-applicable" | "outside-domain",
      reason: string,
    ): ScientificResult {
      return status === "underdetermined"
        ? { ...meta(id), status, compatibleFamily: reason, neededInformation: [reason] }
        : status === "outside-domain"
          ? {
              ...meta(id),
              status,
              reason,
              condition:
                "The inference requires positive, representable inputs in its admitted model.",
              domainKind: "model",
              boundary: { alternativeModel: "Inspect the stated measurement assumptions." },
            }
          : { ...meta(id), status, reason };
    },
  };
}
export const inferenceRefusal = (requirements: string): Assessment<never> => ({
  kind: "refused",
  refusal: makeRefusal(
    "invalid-parameter",
    { capabilityId: "diffusion.inference" },
    { details: { requirements } },
  ),
});
function matchingFields(input: unknown, expected: object): input is Parameters {
  return (
    !!input &&
    typeof input === "object" &&
    [Object.prototype, null].includes(Object.getPrototypeOf(input)) &&
    Reflect.ownKeys(input).length === Object.keys(expected).length &&
    Reflect.ownKeys(input).every(
      (key) =>
        typeof key === "string" &&
        Object.hasOwn(expected, key) &&
        Object.hasOwn(Object.getOwnPropertyDescriptor(input, key) ?? {}, "value"),
    )
  );
}
export function validateRadius(input: unknown): Assessment<Parameters> {
  if (
    !matchingFields(input, RADIUS_DEFAULTS) ||
    typeof input.radiusKnown !== "boolean" ||
    typeof input.provenance !== "string" ||
    input.provenance.length > 512 ||
    !["independently-declared", "same-displacements"].includes(String(input.radiusProvenance))
  )
    return inferenceRefusal("Declare the radius information using the named fields.");
  for (const key of ["radius", "radiusLower", "radiusUpper", "radiusCoverage"])
    if (typeof input[key] !== "number" || !Number.isFinite(input[key]) || input[key] <= 0)
      return inferenceRefusal("Use positive finite calibration inputs.");
  if (
    Number(input.radiusLower) > Number(input.radius) ||
    Number(input.radiusUpper) < Number(input.radius) ||
    Number(input.radiusCoverage) > 1
  )
    return inferenceRefusal(
      "Radius bounds must contain the point value; coverage must not exceed one.",
    );
  if (input.radiusKnown && !input.provenance.trim())
    return inferenceRefusal("State where the independent radius information came from.");
  return { kind: "accepted", data: Object.freeze({ ...input }) };
}
export function evaluateRadius(
  evidence: InferenceEvidence,
  input: unknown,
): Assessment<readonly ScientificResult[]> {
  if (evidence.kind !== "ideal")
    return inferenceRefusal("The radius case requires the ideal displacement example.");
  const valid = validateRadius(input);
  if (valid.kind !== "accepted") return valid;
  const p = valid.data;
  const estimate = estimateIncrements(
    Float64Array.from(evidence.increments),
    evidence.dt,
    evidence.d,
    "independent-increment-known-zero-drift",
  );
  if (estimate.kind !== "accepted") return estimate;
  const band = diffusionCompatibleBand(
    {
      estimate: estimate.data,
      T: evidence.T,
      eta: evidence.eta,
      alpha: 0.05,
      radiusRange: [0.1e-6, 2e-6],
      synthetic: true,
    },
    INFERENCE_CONSTANTS,
  );
  if (band.kind !== "accepted") return band;
  const { value, absent } = outputsFor(FAMILY_OUTPUTS);
  const b = band.data;
  const outputs: ScientificResult[] = [
    value("diffusionCoefficientEstimate", estimate.data.dHat),
    value("diffusionInterval", Float64Array.of(b.interval.lower, b.interval.upper)),
    value("radiusNumberProduct", b.product),
    value("familyRadii", b.radii),
    value("familyNumbers", b.numbers),
    value("familyLowerNumbers", b.lowerNumbers),
    value("familyUpperNumbers", b.upperNumbers),
  ];
  const conditional = conditionalMolecularNumberInterval(
    {
      estimate: estimate.data,
      T: evidence.T,
      eta: evidence.eta,
      alpha: 0.05,
      synthetic: true,
      radius: p.radiusKnown
        ? {
            value: Number(p.radius),
            interval: {
              lower: Number(p.radiusLower),
              upper: Number(p.radiusUpper),
              coverage: Number(p.radiusCoverage),
            },
            provenance: String(p.provenance),
          }
        : null,
      radiusProvenance: p.radiusProvenance as "independently-declared" | "same-displacements",
    },
    INFERENCE_CONSTANTS,
  );
  if (conditional.kind === "refused" || conditional.kind === "outcome") return conditional;
  if (conditional.kind === "no-value") {
    for (const id of ["avogadroNumberEstimate", "molecularInterval", "simultaneousCoverage"])
      outputs.push(absent(id, conditional.status, conditional.reason));
  } else {
    const pointInterval = estimatorInterval(estimate.data, 0.025);
    if (pointInterval.kind !== "accepted") return pointInterval;
    const point = invertToMolecularNumber(
      {
        T: evidence.T,
        eta: evidence.eta,
        a: Number(p.radius),
        radiusProvenance: "independently-declared",
        dHat: estimate.data.dHat,
        interval: pointInterval.data,
        synthetic: true,
      },
      INFERENCE_CONSTANTS,
    );
    if (point.kind !== "accepted") return point;
    outputs.push(
      value("avogadroNumberEstimate", point.data.estimate),
      value("molecularInterval", Float64Array.of(conditional.data.lower, conditional.data.upper)),
      value("simultaneousCoverage", conditional.data.coverage),
    );
  }
  return { kind: "accepted", data: outputs };
}

export function evaluateCamera(
  evidence: InferenceEvidence,
  input: unknown,
): Assessment<readonly ScientificResult[]> {
  if (
    evidence.kind !== "camera" ||
    !matchingFields(input, CAMERA_DEFAULTS) ||
    !["variance", "covariance", "second-interval"].includes(String(input.information))
  )
    return inferenceRefusal(
      "Choose variance alone, neighboring covariance, or a second interval for this fixed camera recording.",
    );
  const { value, absent } = outputsFor(CAMERA_OUTPUTS);
  const sample = covarianceEstimator(Float64Array.from(evidence.increments), evidence.dt, {
    d: evidence.d,
    exposure: evidence.exposure,
    knownDrift: evidence.knownDrift,
  });
  if (sample.kind !== "accepted") return sample;
  const family = cameraCompatibleFamily({
    variance: sample.data.variance,
    dt: evidence.dt,
    exposure: evidence.exposure,
  });
  if (family.kind !== "accepted") return family;
  const f = family.data,
    s = sample.data,
    method = String(input.information);
  const outputs: ScientificResult[] = [
    value("sampleVariance", s.variance),
    value("familyDiffusion", f.diffusion),
    value("familyNoise", f.noiseVariance),
    value("familyCovariance", f.covariance),
    value("maximumDiffusion", f.maximumDiffusion),
    value("maximumNoise", f.maximumNoiseVariance),
    absent(
      "diffusionConfidenceInterval",
      "not-applicable",
      "Correlated camera increments do not satisfy the independent-increment chi-square procedure. These moment solutions are point estimates, not confidence intervals.",
    ),
  ];
  outputs.push(
    method === "covariance"
      ? value("sampleCovariance", s.covariance)
      : absent(
          "sampleCovariance",
          "not-applicable",
          "Neighboring covariance has not been admitted to this question.",
        ),
  );
  let estimate: { D: number; sigma2: number; physical: boolean } | null = null;
  if (method === "second-interval") {
    // Pair adjacent increments from the SAME positions. No regeneration and no new random draw.
    const coarse = new Float64Array(evidence.increments.length / 2);
    for (let i = 0; i < coarse.length / evidence.d; i++)
      for (let c = 0; c < evidence.d; c++) {
        const left = evidence.increments[2 * i * evidence.d + c],
          right = evidence.increments[(2 * i + 1) * evidence.d + c];
        if (left === undefined || right === undefined)
          return inferenceRefusal("The paired observation layout is incomplete.");
        coarse[i * evidence.d + c] = left + right;
      }
    const second = covarianceEstimator(coarse, 2 * evidence.dt, {
      d: evidence.d,
      exposure: evidence.exposure,
      knownDrift: evidence.knownDrift,
    });
    if (second.kind !== "accepted") return second;
    outputs.push(value("secondVariance", second.data.variance));
    const solved = twoIntervalCameraEstimate({
      firstVariance: s.variance,
      secondVariance: second.data.variance,
      firstDt: evidence.dt,
      secondDt: 2 * evidence.dt,
      exposure: evidence.exposure,
    });
    if (solved.kind !== "accepted") return solved;
    estimate = solved.data;
  } else {
    outputs.push(
      absent(
        "secondVariance",
        "not-applicable",
        "A second observation interval has not been admitted to this question.",
      ),
    );
    if (method === "covariance")
      estimate = { D: s.D, sigma2: s.sigma2, physical: s.D >= 0 && s.sigma2 >= 0 };
  }
  if (estimate)
    outputs.push(
      value("diffusionEstimate", estimate.D),
      value("noiseEstimate", estimate.sigma2),
      value("physicalSolution", estimate.physical ? 1 : 0),
    );
  else
    for (const id of ["diffusionEstimate", "noiseEstimate", "physicalSolution"])
      outputs.push(
        absent(
          id,
          "underdetermined",
          "One variance leaves a compatible diffusion–noise family. Admit covariance or a second spacing to separate the two unknowns.",
        ),
      );
  return { kind: "accepted", data: outputs };
}