#' Comprehensive analysis pipeline for embedding method comparison
#'
#' @param data_gen_func Function to generate the dataset
#' @param dataset_name Name for the dataset (for labeling)
#' @param n_objects Number of objects in the dataset
#' @param n_runs Number of runs for stochastic methods
run_comprehensive_analysis <- function(data_gen_func, dataset_name, n_objects = 50, n_runs = 50) {
cat("=== ANALYSIS FOR", toupper(dataset_name), "===\n")
# Step 1: Data Generation
cat("\n1. Generating", dataset_name, "dataset...\n")
data_gen_output <- data_gen_func(n_objects = n_objects, missing_fraction = 0.3)
complete_distances_for_evaluation <- data_gen_output$complete_distances
incomplete_distances_for_embedding <- data_gen_output$incomplete_distances
object_names <- data_gen_output$object_names
actual_missing_percentage <- sum(is.na(incomplete_distances_for_embedding)) /
(n_objects * (n_objects-1)) * 100
data_quality <- calculate_data_quality_metrics(incomplete_distances_for_embedding)
cat("Generated dataset with", n_objects, "objects and",
round(actual_missing_percentage, 1), "% missing values\n")
# Step 2: Assess Non-Euclidean Character
cat("\n2. Assessing non-Euclidean character...\n")
D_squared <- complete_distances_for_evaluation^2
n <- nrow(D_squared)
J <- diag(n) - (1/n) * matrix(1, n, n)
B <- -0.5 * J %*% D_squared %*% J
eigenvals <- eigen(B, symmetric = TRUE)$values
numerical_tolerance <- 1e-12
positive_eigenvals <- eigenvals[eigenvals > numerical_tolerance]
negative_eigenvals <- eigenvals[eigenvals < -numerical_tolerance]
total_variance <- sum(abs(eigenvals))
negative_variance <- sum(abs(negative_eigenvals))
positive_variance <- sum(positive_eigenvals)
if(positive_variance > 0) {
deviation_score <- negative_variance / positive_variance
} else {
deviation_score <- 1.0
}
negative_variance_fraction <- negative_variance / total_variance
cumulative_positive_variance <- cumsum(positive_eigenvals) / positive_variance
dims_for_90_percent <- which(cumulative_positive_variance >= 0.90)[1]
if(is.na(dims_for_90_percent)) dims_for_90_percent <- length(positive_eigenvals)
cat("Non-Euclidean deviation score:", round(deviation_score, 4), "\n")
cat("Dimensions for 90% variance:", dims_for_90_percent, "\n")
# Step 3: Parameter Optimization using Euclidify
cat("\n3. Running automated parameter optimization...\n")
# Create temporary directory for optimization
temp_dir <- tempfile()
dir.create(temp_dir, recursive = TRUE)
euclidify_results <- tryCatch({
Euclidify(
dissimilarity_matrix = incomplete_distances_for_embedding,
output_dir = temp_dir,
ndim_range = c(2, min(10, dims_for_90_percent + 3)),
n_initial_samples = 20, # Reduced for vignette speed
n_adaptive_samples = 40, # Reduced for vignette speed
folds = 20,
mapping_max_iter = 300, # Reduced for speed
clean_intermediate = TRUE,
verbose = "off",
fallback_to_defaults = TRUE,
save_results = FALSE
)
}, error = function(e) {
cat("Euclidify failed, using default parameters\n")
NULL
})
# Clean up temp directory
unlink(temp_dir, recursive = TRUE)
if(!is.null(euclidify_results)) {
optimal_params <- euclidify_results$optimal_params
euclidify_positions <- euclidify_results$positions
target_dims <- optimal_params$ndim
cat("Optimal parameters found - dimensions:", target_dims, "\n")
} else {
# Fallback parameters
target_dims <- max(2, min(5, dims_for_90_percent))
optimal_params <- list(
ndim = target_dims,
k0 = 5.0,
cooling_rate = 0.01,
c_repulsion = 0.02
)
euclidify_positions <- NULL
cat("Using fallback parameters - dimensions:", target_dims, "\n")
}
# Step 4: Three-Method Comparison
cat("\n4. Running three-method comparison...\n")
# Initialize storage
topolow_results <- vector("list", n_runs)
iterative_mds_results <- vector("list", n_runs)
classical_mds_result <- NULL
topolow_errors <- numeric(n_runs)
iterative_mds_errors <- numeric(n_runs)
classical_mds_error <- NA
# Topolow runs
cat("Running Topolow embeddings...\n")
for(i in 1:n_runs) {
if(i == 1 && !is.null(euclidify_positions)) {
# Use Euclidify result for first run
topolow_coords <- euclidify_positions
} else {
# Run new embedding
topolow_result <- tryCatch({
euclidean_embedding(
dissimilarity_matrix = incomplete_distances_for_embedding,
ndim = optimal_params$ndim,
mapping_max_iter = 300,
k0 = optimal_params$k0,
cooling_rate = optimal_params$cooling_rate,
c_repulsion = optimal_params$c_repulsion,
relative_epsilon = 1e-6,
convergence_counter = 3,
verbose = FALSE
)
}, error = function(e) NULL)
if(!is.null(topolow_result)) {
topolow_coords <- topolow_result$positions
} else {
topolow_coords <- NULL
}
}
if(!is.null(topolow_coords)) {
topolow_coords <- topolow_coords[order(row.names(topolow_coords)), ]
if(validate_coordinates(topolow_coords, "Topolow", n_objects, target_dims)) {
topolow_coords <- scale(topolow_coords, center = TRUE, scale = FALSE)
topolow_results[[i]] <- topolow_coords
embedded_distances <- as.matrix(dist(topolow_coords))
rownames(embedded_distances) <- rownames(topolow_coords)
colnames(embedded_distances) <- rownames(topolow_coords)
valid_mask <- !is.na(complete_distances_for_evaluation)
distance_errors <- abs(complete_distances_for_evaluation[valid_mask] -
embedded_distances[valid_mask])
topolow_errors[i] <- sum(distance_errors)
} else {
topolow_results[[i]] <- NULL
topolow_errors[i] <- NA
}
} else {
topolow_results[[i]] <- NULL
topolow_errors[i] <- NA
}
}
# Classical MDS
cat("Running Classical MDS...\n")
imputation_result <- improved_missing_data_imputation(incomplete_distances_for_embedding)
classical_mds_matrix <- imputation_result$matrix
tryCatch({
classical_mds_result_raw <- cmdscale(classical_mds_matrix, k = target_dims, eig = TRUE)
classical_mds_coords <- classical_mds_result_raw$points
rownames(classical_mds_coords) <- object_names
classical_mds_coords <- classical_mds_coords[order(row.names(classical_mds_coords)), ]
if(validate_coordinates(classical_mds_coords, "Classical MDS", n_objects, target_dims)) {
classical_mds_coords <- scale(classical_mds_coords, center = TRUE, scale = FALSE)
classical_mds_result <- classical_mds_coords
embedded_distances <- as.matrix(dist(classical_mds_coords))
rownames(embedded_distances) <- rownames(classical_mds_coords)
colnames(embedded_distances) <- rownames(classical_mds_coords)
valid_mask <- !is.na(complete_distances_for_evaluation)
distance_errors <- abs(complete_distances_for_evaluation[valid_mask] -
embedded_distances[valid_mask])
classical_mds_error <- sum(distance_errors)
}
}, error = function(e) {
cat("Classical MDS failed:", e$message, "\n")
})
# Iterative MDS
cat("Running Iterative MDS...\n")
for(i in 1:n_runs) {
tryCatch({
if(requireNamespace("smacof", quietly = TRUE)) {
max_dist <- max(classical_mds_matrix, na.rm = TRUE)
scaled_matrix <- classical_mds_matrix / max_dist
iterative_mds_result_raw <- smacof::smacofSym(
delta = scaled_matrix,
ndim = target_dims,
type = "interval",
init = "torgerson",
verbose = FALSE,
itmax = 1000,
eps = 1e-5
)
iterative_mds_coords <- iterative_mds_result_raw$conf
current_max_dist <- max(dist(iterative_mds_coords))
if(current_max_dist > 0) {
scale_factor <- max_dist / current_max_dist
iterative_mds_coords <- iterative_mds_coords * scale_factor
}
} else {
# Fallback to isoMDS
iterative_mds_result_raw <- MASS::isoMDS(classical_mds_matrix, k = target_dims, trace = FALSE)
iterative_mds_coords <- iterative_mds_result_raw$points
}
rownames(iterative_mds_coords) <- object_names
iterative_mds_coords <- iterative_mds_coords[order(row.names(iterative_mds_coords)), ]
if(validate_coordinates(iterative_mds_coords, "Iterative MDS", n_objects, target_dims)) {
iterative_mds_coords <- scale(iterative_mds_coords, center = TRUE, scale = FALSE)
iterative_mds_results[[i]] <- iterative_mds_coords
embedded_distances <- as.matrix(dist(iterative_mds_coords))
rownames(embedded_distances) <- rownames(iterative_mds_coords)
colnames(embedded_distances) <- rownames(iterative_mds_coords)
valid_mask <- !is.na(complete_distances_for_evaluation)
distance_errors <- abs(complete_distances_for_evaluation[valid_mask] -
embedded_distances[valid_mask])
iterative_mds_errors[i] <- sum(distance_errors)
} else {
iterative_mds_results[[i]] <- NULL
iterative_mds_errors[i] <- NA
}
}, error = function(e) {
iterative_mds_results[[i]] <- NULL
iterative_mds_errors[i] <- NA
})
}
# Step 5: Compile Results
valid_topolow_results <- !is.na(topolow_errors)
valid_iterative_results <- !is.na(iterative_mds_errors)
results <- list(
dataset_name = dataset_name,
data_characteristics = list(
n_objects = n_objects,
missing_percentage = actual_missing_percentage,
deviation_score = deviation_score,
dims_90_percent = dims_for_90_percent,
negative_variance_fraction = negative_variance_fraction
),
optimal_params = optimal_params,
topolow_errors = topolow_errors,
topolow_results = topolow_results,
valid_topolow_results = valid_topolow_results,
classical_mds_error = classical_mds_error,
classical_mds_result = classical_mds_result,
iterative_mds_errors = iterative_mds_errors,
iterative_mds_results = iterative_mds_results,
valid_iterative_results = valid_iterative_results,
complete_distances = complete_distances_for_evaluation,
incomplete_distances = incomplete_distances_for_embedding,
target_dims = target_dims
)
return(results)
}
