const math = require('mathjs');

/**
 * Compute PCA (eigenvectors and eigenvalues) for deviation data to assess geometric repeatability.
 * @param {Array<Array<number>>} deviations - Array of [point, line, i_deviation, j_deviation]
 * @returns {Object} - { eigenvalues, eigenvectors, rms, anisotropy, primaryDirection }
 */
function computePCA(deviations) {
	// Extract i_deviation and j_deviation
	const deviationMatrix = deviations.map(row => [row[2], row[3]]);

	// Convert to mathjs matrix
	const D = math.matrix(deviationMatrix);

	// Compute mean for centering (1 x 2 matrix)
	const mean = math.mean(D, 0);

	// Explicitly repeat mean to match D's shape (n x 2)
	const n = deviationMatrix.length;
	const meanRepeated = math.repmat(mean, n, 1);

	// Center the data
	const centered = math.subtract(D, meanRepeated);

	// Compute covariance matrix: (1/(n-1)) * (D_centered^T * D_centered)
	const covMatrix = math.multiply(
		math.multiply(1 / (n - 1), math.transpose(centered)),
		centered
	);

	// Perform eigen decomposition
	const result = math.eigs(covMatrix);
	let eigenvalues = result.values;
	const evObjs = result.eigenvectors;

	// Convert eigenvalues to array if it's a matrix
	eigenvalues = Array.isArray(eigenvalues) ? eigenvalues : eigenvalues.toArray();

	// Create pairs and convert vector to array if necessary
	const pairs = eigenvalues.map((val, i) => {
		let vec = evObjs[i].vector;
		if (vec.toArray) vec = vec.toArray();
		return { val, vec };
	});

	// Sort by descending eigenvalues
	pairs.sort((a, b) => b.val - a.val);

	// Sorted eigenvalues
	const sortedEigenvalues = pairs.map(p => p.val);

	// Build eigenvector matrix: rows as components, columns as eigenvectors
	const dimension = pairs[0].vec.length; // e.g., 2
	const evecRows = [];
	for (let comp = 0; comp < dimension; comp++) {
		evecRows.push(pairs.map(p => p.vec[comp]));
	}
	const sortedEigenvectors = math.matrix(evecRows);

	// Compute RMS errors along principal axes
	const rms = sortedEigenvalues.map(val => Math.sqrt(Math.max(val, 0)));

	// Compute anisotropy (ratio of major to minor axis variance)
	const anisotropy = sortedEigenvalues[0] / (sortedEigenvalues[1] || 1); // Avoid division by zero

	// Primary direction (angle in degrees of major eigenvector)
	const primaryVector = sortedEigenvectors.subset(math.index([0, 1], 0)).toArray();
	const primaryDirection = Math.atan2(primaryVector[1], primaryVector[0]) * 180 / Math.PI;

	return {
		eigenvalues: sortedEigenvalues,
		eigenvectors: sortedEigenvectors.toArray(),
		rms: rms, // RMS errors along major/minor axes
		anisotropy: anisotropy, // Ratio of variances
		primaryDirection: primaryDirection // Angle of major axis (degrees)
	};
}


module.exports = { computePCA };
