#' Matrix Metrics
#'
#' Measure the performance of a covariance or correlation
#' matrix estimate with various metrics.
#'
#' @param S_true The true covariance matrix.
#' @param R_hat The estimate correlation matrix.
#' @param strong Eigenvalues strength threshold (default 0.75).
#'
#' @return A data.frame with 5 metrics.
#'
#' @examples
#' X_obj <- simulate_data(n = 10, p = 10)
#' X <- X_obj$X
#' S_true <- X_obj$S
#' R_hat <- cor(X)
#' results <- matrix_metrics(S_true, R_hat)
#'
#' @export
matrix_metrics <- function(S_true, R_hat, strong = 0.75){
	p <- nrow(S_true)
	D <- S_true - R_hat

	rmse <- sqrt(mean(D^2))
	# NOTE: Serra's text is ambiguous on squared-vs-not. This uses ||.||_F / sqrt(p),
	# the reading under which the identity matrix normalises to exactly 1. If you
	# want the squared version, use sum(D^2)/p instead. Verify against RSC source.
	L_F <- sqrt(sum(D^2)) / sqrt(p)

	ev_t <- sort(eigen(S_true, symmetric = TRUE, only.values = TRUE)$values)
	ev_h <- sort(eigen(R_hat,  symmetric = TRUE, only.values = TRUE)$values)
	L_S  <- sum(abs(ev_t - ev_h))

	lt <- lower.tri(S_true)          # unique off-diagonal pairs (v < l)
	t0 <- S_true[lt]
	th <- R_hat[lt]

	E1 <- sum( (t0 == 0) & (abs(th) > strong) )
	E2 <- sum( (th == 0) & (abs(t0) > strong) )
	SS <- sum( (t0 > 0 & th < 0) | (t0 < 0 & th > 0) )

	data.frame(
		rmse = rmse,
		L_F = L_F,
		L_S = L_S,
		E1 = E1,
		E2 = E2,
		SS = SS
	)
}