feat: Add methods to Psi to calculate D-optimality

[mhovd]

Resolves #238

Note: if we convert psi to contain log-likelihoods, then the calculation of D-optimality must be updated!

[Copilot]

Pull request overview

This PR implements D-optimality calculations for the Psi struct to measure convergence in the NPML/NPOD algorithm, resolving issue #238. The D-optimality criterion should approach 0 at convergence, indicating that no support point can further improve the likelihood.

Key changes:

• Added two methods for computing D-optimality: a public d_optimality() that returns the maximum value across all support points, and an internal d_optimality_spp() that returns per-support-point values
• Comprehensive test suite with 8 test cases covering edge cases including uniform/non-uniform weights, convergence scenarios, dimension mismatches, and zero probability errors

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

Direct floating-point comparison with zero may be unreliable due to floating-point arithmetic precision issues. Consider using a small epsilon threshold instead, such as checking if the absolute value is less than a small tolerance (e.g., 1e-12 or f64::EPSILON).

Suggested change

	// Check for zero probabilities
	for (i, &p) in pyl.iter().enumerate() {
	if p == 0.0 {
	// Check for (effectively) zero probabilities using a small tolerance
	let eps = 1e-12_f64;
	for (i, &p) in pyl.iter().enumerate() {
	if p.abs() < eps {

[Copilot]

The mathematical notation using dollar signs (e.g., $j$, $$D(\theta_j) = ...$$) is LaTeX/MathJax syntax that does not render correctly in rustdoc. Rustdoc does not natively support LaTeX math expressions. Consider using plain text mathematical notation or code formatting instead, such as "D(theta_j) = sum(psi_ij / p_lambda(y_i)) - n" in a code block or plain description.

Suggested change

	/// The D-optimality value for support point $j$ is:
	/// $$D(\theta_j) = \sum_{i=1}^{n} \frac{\psi_{ij}}{p_\lambda(y_i)} - n$$
	/// The D-optimality value for support point j is:
	/// ```text
	/// D(theta_j) = sum_{i=1}^n psi_ij / p_lambda(y_i) - n
	/// ```

[Copilot]

The nested loop pattern for computing weighted probabilities could be more efficient using matrix-vector multiplication. The faer library used for the matrix likely provides optimized operations for this. Consider using matrix-vector multiplication: pyl = psi_mat * weights.weights() instead of manual iteration, which would be both more performant and more readable.

Suggested change

	let mut pyl = vec![0.0; nsub];
	for i in 0..nsub {
	for (j, w_j) in weights.iter().enumerate() {
	pyl[i] += psi_mat.get(i, j) * w_j;
	}
	// Build a column vector (nspp x 1) of weights in column-major order.
	let weights_vec: Vec<f64> = weights.iter().copied().collect();
	// Treat the weights as a (nspp x 1) matrix so we can use faer matrix multiplication.
	let weights_mat = Mat::from_column_major_slice(nspp, 1, &weights_vec);
	// psi_mat has shape (nsub x nspp), so the product has shape (nsub x 1).
	let pyl_mat = &psi_mat * &weights_mat;
	let mut pyl = Vec::with_capacity(nsub);
	for i in 0..nsub {
	// Extract the single column of the result matrix.
	pyl.push(*pyl_mat.get(i, 0));