481 lines
16 KiB
JavaScript
481 lines
16 KiB
JavaScript
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,
|
|
};
|
|
}
|