const FLOAT64_BYTES = Float64Array.BYTES_PER_ELEMENT;

export const DEFAULT_MAX_COVARIANCE_WORK_BYTES = 32 * 1024 * 1024;
export const DEFAULT_MAX_COVARIANCE_COMPLEX_ACCUMULATIONS = 250_000_000;
export const DEFAULT_COVARIANCE_EIGENVALUE_ROUNDOFF_FACTOR = 64;

function assertPositiveInteger(value, label) {
  if (!Number.isSafeInteger(value) || value < 1) {
    throw new TypeError(`${label} must be a positive safe integer.`);
  }
}

function checkedProduct(factors, label) {
  let product = 1;
  for (const factor of factors) {
    if (factor !== 0 && product > Number.MAX_SAFE_INTEGER / factor) {
      throw new RangeError(`${label} exceeds the safe integer range.`);
    }
    product *= factor;
  }
  return product;
}

function roundoffGamma(operationCount) {
  const product = operationCount * (Number.EPSILON / 2);
  if (!(product < 1)) {
    throw new RangeError("roundoff bound is undefined for this workload.");
  }
  return product / (1 - product);
}

/** Estimate covariance output allocation and upper-triangle accumulation work. */
export function estimateComplexCovarianceWorkload({ sensorCount, snapshotCount }) {
  assertPositiveInteger(sensorCount, "sensorCount");
  assertPositiveInteger(snapshotCount, "snapshotCount");
  const inputElementCount = checkedProduct(
    [sensorCount, snapshotCount],
    "covariance input element count",
  );
  const matrixElementCount = checkedProduct(
    [sensorCount, sensorCount],
    "covariance matrix element count",
  );
  const upperTriangleElementCount =
    checkedProduct([sensorCount, sensorCount + 1], "covariance upper triangle") / 2;
  const complexAccumulations = checkedProduct(
    [upperTriangleElementCount, snapshotCount],
    "covariance accumulation count",
  );
  const allocationBytes = checkedProduct(
    [2, matrixElementCount, FLOAT64_BYTES],
    "covariance allocation",
  );
  return Object.freeze({
    sensorCount,
    snapshotCount,
    inputElementCount,
    matrixElementCount,
    upperTriangleElementCount,
    complexAccumulations,
    allocationBytes,
  });
}

function compensatedAdd(sum, compensation, value) {
  const next = sum + value;
  const correction = Math.abs(sum) >= Math.abs(value)
    ? (sum - next) + value
    : (value - next) + sum;
  return [next, compensation + correction];
}

/**
 * Compute R = (1/T) sum_t x[t] x[t]^H in Float64.
 *
 * `snapshot-major` stores x(t,m) at `t * sensorCount + m`; `sensor-major`
 * stores it at `m * snapshotCount + t`. Only the upper triangle is
 * accumulated, with Neumaier compensation, and the lower triangle is filled
 * as its exact complex conjugate. No centering or diagonal loading is applied.
 */
export function computeComplexSampleCovariance({
  real,
  imaginary,
  sensorCount,
  snapshotCount,
  layout = "snapshot-major",
  maximumWorkBytes = DEFAULT_MAX_COVARIANCE_WORK_BYTES,
  maximumComplexAccumulations = DEFAULT_MAX_COVARIANCE_COMPLEX_ACCUMULATIONS,
}) {
  const workload = estimateComplexCovarianceWorkload({ sensorCount, snapshotCount });
  if (layout !== "snapshot-major" && layout !== "sensor-major") {
    throw new RangeError("layout must be 'snapshot-major' or 'sensor-major'.");
  }
  if (real == null || imaginary == null) {
    throw new TypeError("real and imaginary snapshots are required.");
  }
  if (
    real.length !== workload.inputElementCount ||
    imaginary.length !== workload.inputElementCount
  ) {
    throw new RangeError(
      `snapshot components must each contain ${workload.inputElementCount} values.`,
    );
  }
  if (!(maximumWorkBytes >= 0) || !Number.isFinite(maximumWorkBytes)) {
    throw new TypeError("maximumWorkBytes must be finite and nonnegative.");
  }
  if (
    !(maximumComplexAccumulations >= 0) ||
    !Number.isFinite(maximumComplexAccumulations)
  ) {
    throw new TypeError("maximumComplexAccumulations must be finite and nonnegative.");
  }
  if (workload.allocationBytes > maximumWorkBytes) {
    throw new RangeError(
      `covariance requires ${workload.allocationBytes} typed-array bytes; limit is ${maximumWorkBytes}.`,
    );
  }
  if (workload.complexAccumulations > maximumComplexAccumulations) {
    throw new RangeError(
      `covariance requires ${workload.complexAccumulations} complex accumulations; limit is ${maximumComplexAccumulations}.`,
    );
  }

  for (let index = 0; index < workload.inputElementCount; index += 1) {
    if (!Number.isFinite(real[index]) || !Number.isFinite(imaginary[index])) {
      throw new TypeError("snapshot components must be finite.");
    }
  }

  const matrixReal = new Float64Array(workload.matrixElementCount);
  const matrixImaginary = new Float64Array(workload.matrixElementCount);
  const sampleIndex = layout === "snapshot-major"
    ? (sensor, snapshot) => snapshot * sensorCount + sensor
    : (sensor, snapshot) => sensor * snapshotCount + snapshot;
  let matrixScale = 0;
  let traceSum = 0;
  let traceCompensation = 0;

  for (let row = 0; row < sensorCount; row += 1) {
    for (let column = row; column < sensorCount; column += 1) {
      let sumReal = 0;
      let compensationReal = 0;
      let sumImaginary = 0;
      let compensationImaginary = 0;
      for (let snapshot = 0; snapshot < snapshotCount; snapshot += 1) {
        const rowIndex = sampleIndex(row, snapshot);
        const columnIndex = sampleIndex(column, snapshot);
        const rowReal = real[rowIndex];
        const rowImaginary = imaginary[rowIndex];
        const columnReal = real[columnIndex];
        const columnImaginary = imaginary[columnIndex];
        const productReal =
          rowReal * columnReal + rowImaginary * columnImaginary;
        const productImaginary =
          rowImaginary * columnReal - rowReal * columnImaginary;
        if (!Number.isFinite(productReal) || !Number.isFinite(productImaginary)) {
          throw new RangeError("covariance product exceeds the finite Float64 range.");
        }
        [sumReal, compensationReal] = compensatedAdd(
          sumReal,
          compensationReal,
          productReal,
        );
        [sumImaginary, compensationImaginary] = compensatedAdd(
          sumImaginary,
          compensationImaginary,
          productImaginary,
        );
      }
      const valueReal = (sumReal + compensationReal) / snapshotCount;
      const valueImaginary = row === column
        ? 0
        : (sumImaginary + compensationImaginary) / snapshotCount;
      if (!Number.isFinite(valueReal) || !Number.isFinite(valueImaginary)) {
        throw new RangeError("covariance value exceeds the finite Float64 range.");
      }
      const upper = row * sensorCount + column;
      const lower = column * sensorCount + row;
      matrixReal[upper] = valueReal;
      matrixImaginary[upper] = valueImaginary;
      matrixReal[lower] = valueReal;
      matrixImaginary[lower] = row === column ? 0 : -valueImaginary;
      matrixScale = Math.max(
        matrixScale,
        Math.abs(valueReal),
        Math.abs(valueImaginary),
      );
      if (row === column) {
        [traceSum, traceCompensation] = compensatedAdd(
          traceSum,
          traceCompensation,
          valueReal,
        );
      }
    }
  }

  const trace = traceSum + traceCompensation;
  const covarianceRoundoffBound =
    roundoffGamma(4 * snapshotCount + 32) * trace;

  return {
    sensorCount,
    snapshotCount,
    inputLayout: layout,
    normalization: "1/T",
    centered: false,
    matrix: {
      real: matrixReal,
      imaginary: matrixImaginary,
      layout: "row-major",
    },
    matrixScale,
    trace,
    covarianceRoundoffBound,
    workload,
  };
}

/**
 * Audit covariance eigenvalues against a scale-aware backward-roundoff bound.
 * Strictly positive values are never floored. Only negative values within the
 * bound are represented as zero in `roundoffClampedEigenvalues`; raw values
 * are always retained separately.
 */
export function auditCovarianceEigenvalues({
  eigenvalues,
  sensorCount = eigenvalues?.length,
  matrixScale,
  roundoffTolerance,
  covarianceRoundoffBound,
  eigensolverResidualFrobeniusNorm = 0,
  roundoffFactor = DEFAULT_COVARIANCE_EIGENVALUE_ROUNDOFF_FACTOR,
  rejectIndefinite = true,
}) {
  if (eigenvalues == null) throw new TypeError("eigenvalues are required.");
  assertPositiveInteger(sensorCount, "sensorCount");
  if (eigenvalues.length !== sensorCount) {
    throw new RangeError("eigenvalue count must equal sensorCount.");
  }
  if (!(roundoffFactor > 0) || !Number.isFinite(roundoffFactor)) {
    throw new TypeError("roundoffFactor must be finite and positive.");
  }
  if (matrixScale !== undefined && (!(matrixScale >= 0) || !Number.isFinite(matrixScale))) {
    throw new TypeError("matrixScale must be finite and nonnegative.");
  }
  if (
    covarianceRoundoffBound !== undefined &&
    (!(covarianceRoundoffBound >= 0) || !Number.isFinite(covarianceRoundoffBound))
  ) {
    throw new TypeError("covarianceRoundoffBound must be finite and nonnegative.");
  }
  if (
    !(eigensolverResidualFrobeniusNorm >= 0) ||
    !Number.isFinite(eigensolverResidualFrobeniusNorm)
  ) {
    throw new TypeError(
      "eigensolverResidualFrobeniusNorm must be finite and nonnegative.",
    );
  }

  let eigenvalueScale = matrixScale ?? 0;
  let minimumEigenvalue = Infinity;
  for (const value of eigenvalues) {
    if (!Number.isFinite(value)) throw new TypeError("eigenvalues must be finite.");
    eigenvalueScale = Math.max(eigenvalueScale, Math.abs(value));
    minimumEigenvalue = Math.min(minimumEigenvalue, value);
  }
  const scaleDerivedRoundoffTolerance = eigenvalueScale === 0
    ? 0
    : roundoffFactor * Number.EPSILON * Math.max(1, sensorCount) * eigenvalueScale;
  const effectiveRoundoffTolerance = roundoffTolerance ?? (
    covarianceRoundoffBound === undefined
      ? scaleDerivedRoundoffTolerance
      : covarianceRoundoffBound + eigensolverResidualFrobeniusNorm
  );
  if (
    !(effectiveRoundoffTolerance >= 0) ||
    !Number.isFinite(effectiveRoundoffTolerance)
  ) {
    throw new TypeError("roundoffTolerance must be finite and nonnegative.");
  }

  const rawEigenvalues = Float64Array.from(eigenvalues);
  const roundoffClampedEigenvalues = Float64Array.from(eigenvalues);
  let clampedNegativeCount = 0;
  let firstIndefiniteIndex = -1;
  for (let index = 0; index < sensorCount; index += 1) {
    const value = rawEigenvalues[index];
    if (value < -effectiveRoundoffTolerance) {
      if (firstIndefiniteIndex < 0) firstIndefiniteIndex = index;
    } else if (value < 0 || Object.is(value, -0)) {
      roundoffClampedEigenvalues[index] = 0;
      clampedNegativeCount += 1;
    }
  }
  const isPositiveSemidefinite = firstIndefiniteIndex < 0;
  if (!isPositiveSemidefinite && rejectIndefinite) {
    const error = new RangeError(
      `covariance is genuinely indefinite: eigenvalue ${rawEigenvalues[firstIndefiniteIndex]} is below -${effectiveRoundoffTolerance}.`,
    );
    error.eigenvalueIndex = firstIndefiniteIndex;
    error.eigenvalue = rawEigenvalues[firstIndefiniteIndex];
    error.roundoffTolerance = effectiveRoundoffTolerance;
    throw error;
  }

  return {
    rawEigenvalues,
    roundoffClampedEigenvalues,
    roundoffTolerance: effectiveRoundoffTolerance,
    roundoffFactor,
    roundoffDerivation: roundoffTolerance !== undefined
      ? "explicit"
      : covarianceRoundoffBound === undefined
        ? "eigenvalue-scale"
        : "covariance-accumulation-plus-eigensolver-residual",
    covarianceRoundoffBound: covarianceRoundoffBound ?? null,
    eigensolverResidualFrobeniusNorm,
    eigenvalueScale,
    minimumEigenvalue,
    clampedNegativeCount,
    firstIndefiniteIndex,
    isPositiveSemidefinite,
  };
}

function unavailableSourceCount(criterion, reason, details = {}) {
  return {
    available: false,
    criterion,
    sourceCount: null,
    score: null,
    scores: null,
    reason,
    ...details,
  };
}

/**
 * Select source count with the Wax-Kailath AIC or MDL criterion.
 *
 * This pure function deliberately declines an answer unless T > M, snapshots
 * are declared independent, and every eigenvalue used by the candidate set is
 * strictly above the supplied roundoff bound. It performs no loading, rank
 * repair, or cross-frequency source-count fusion.
 */
export function selectWaxKailathSourceCount({
  eigenvalues,
  snapshotCount,
  independentSnapshots,
  criterion = "mdl",
  roundoffTolerance,
}) {
  if (criterion !== "aic" && criterion !== "mdl") {
    throw new RangeError("criterion must be 'aic' or 'mdl'.");
  }
  if (eigenvalues == null || eigenvalues.length < 1) {
    throw new TypeError("at least one eigenvalue is required.");
  }
  assertPositiveInteger(snapshotCount, "snapshotCount");
  const sensorCount = eigenvalues.length;
  assertPositiveInteger(sensorCount, "sensorCount");
  let scale = 0;
  for (const value of eigenvalues) {
    if (!Number.isFinite(value)) throw new TypeError("eigenvalues must be finite.");
    scale = Math.max(scale, Math.abs(value));
  }
  const effectiveRoundoffTolerance = roundoffTolerance ?? (
    scale === 0
      ? 0
      : DEFAULT_COVARIANCE_EIGENVALUE_ROUNDOFF_FACTOR *
        Number.EPSILON * Math.max(1, sensorCount) * scale
  );
  if (
    !(effectiveRoundoffTolerance >= 0) ||
    !Number.isFinite(effectiveRoundoffTolerance)
  ) {
    throw new TypeError("roundoffTolerance must be finite and nonnegative.");
  }

  if (independentSnapshots !== true) {
    return unavailableSourceCount(criterion, "snapshots-not-independent", {
      sensorCount,
      snapshotCount,
      roundoffTolerance: effectiveRoundoffTolerance,
    });
  }
  if (snapshotCount <= sensorCount) {
    return unavailableSourceCount(criterion, "snapshot-count-not-greater-than-sensors", {
      sensorCount,
      snapshotCount,
      roundoffTolerance: effectiveRoundoffTolerance,
    });
  }

  for (let index = 0; index < sensorCount; index += 1) {
    const value = eigenvalues[index];
    if (index > 0 && value > eigenvalues[index - 1] + effectiveRoundoffTolerance) {
      throw new RangeError("eigenvalues must be sorted in descending order.");
    }
    if (!(value > effectiveRoundoffTolerance)) {
      return unavailableSourceCount(
        criterion,
        "noise-eigenvalue-not-strictly-above-roundoff",
        {
          sensorCount,
          snapshotCount,
          roundoffTolerance: effectiveRoundoffTolerance,
          failingEigenvalueIndex: index,
          failingEigenvalue: value,
        },
      );
    }
  }

  const scores = new Float64Array(sensorCount);
  let selectedSourceCount = 0;
  let selectedScore = Infinity;
  let logRatioRoundoffClampCount = 0;
  for (let sourceCount = 0; sourceCount < sensorCount; sourceCount += 1) {
    const noiseCount = sensorCount - sourceCount;
    let arithmeticSum = 0;
    let arithmeticCompensation = 0;
    let logarithmSum = 0;
    let logarithmCompensation = 0;
    for (let index = sourceCount; index < sensorCount; index += 1) {
      [arithmeticSum, arithmeticCompensation] = compensatedAdd(
        arithmeticSum,
        arithmeticCompensation,
        eigenvalues[index],
      );
      [logarithmSum, logarithmCompensation] = compensatedAdd(
        logarithmSum,
        logarithmCompensation,
        Math.log(eigenvalues[index]),
      );
    }
    const arithmeticMean =
      (arithmeticSum + arithmeticCompensation) / noiseCount;
    const logarithmicGeometricMean =
      (logarithmSum + logarithmCompensation) / noiseCount;
    let logarithmicRatio = Math.log(arithmeticMean) - logarithmicGeometricMean;
    const ratioRoundoffTolerance =
      64 * Number.EPSILON * Math.max(1, noiseCount) *
      Math.max(1, Math.abs(Math.log(arithmeticMean)), Math.abs(logarithmicGeometricMean));
    if (logarithmicRatio < 0 && logarithmicRatio >= -ratioRoundoffTolerance) {
      logarithmicRatio = 0;
      logRatioRoundoffClampCount += 1;
    }
    if (!(logarithmicRatio >= 0) || !Number.isFinite(logarithmicRatio)) {
      throw new RangeError("Wax-Kailath arithmetic/geometric mean ratio is invalid.");
    }

    const likelihoodTerm = snapshotCount * noiseCount * logarithmicRatio;
    const parameterCount = sourceCount * (2 * sensorCount - sourceCount);
    const score = criterion === "aic"
      ? 2 * likelihoodTerm + 2 * parameterCount
      : likelihoodTerm + 0.5 * parameterCount * Math.log(snapshotCount);
    scores[sourceCount] = score;
    if (score < selectedScore) {
      selectedScore = score;
      selectedSourceCount = sourceCount;
    }
  }

  return {
    available: true,
    criterion,
    sourceCount: selectedSourceCount,
    score: selectedScore,
    scores,
    reason: null,
    sensorCount,
    snapshotCount,
    roundoffTolerance: effectiveRoundoffTolerance,
    logRatioRoundoffClampCount,
  };
}