{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.4.0"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":9110664,"sourceType":"datasetVersion","datasetId":5498833},{"sourceId":9177563,"sourceType":"datasetVersion","datasetId":5546655},{"sourceId":13553229,"sourceType":"datasetVersion","datasetId":8608145}],"dockerImageVersionId":30751,"isInternetEnabled":true,"language":"r","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# NeuralIPS Ariel Data Challenge","metadata":{}},{"cell_type":"markdown","source":"#### installing necessary packages","metadata":{}},{"cell_type":"code","source":"install.packages(c( \"tensorflow\", \"reticulate\", \"ggplot2\", \"abind\"))","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true,"execution":{"iopub.status.busy":"2025-10-30T13:48:33.806685Z","iopub.execute_input":"2025-10-30T13:48:33.808521Z","iopub.status.idle":"2025-10-30T13:50:11.155653Z","shell.execute_reply":"2025-10-30T13:50:11.15371Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"install.packages(\"keras3\")  # Not \"keras\"\nkeras3::install_keras()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-30T13:50:11.158303Z","iopub.execute_input":"2025-10-30T13:50:11.197576Z","iopub.status.idle":"2025-10-30T13:52:06.049024Z","shell.execute_reply":"2025-10-30T13:52:06.047368Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#### Import libraries","metadata":{}},{"cell_type":"code","source":"py_require(\"tensorflow\")\npy_require(\"keras3\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-10-30T13:52:06.051436Z","iopub.execute_input":"2025-10-30T13:52:06.053186Z","iopub.status.idle":"2025-10-30T13:52:06.095776Z","shell.execute_reply":"2025-10-30T13:52:06.064831Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"library(keras3)\nlibrary(tensorflow)\nlibrary(reticulate)\nlibrary(ggplot2)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#### Set random Seed","metadata":{}},{"cell_type":"code","source":"set.seed(42)\ntensorflow::tf$random$set_seed(42L)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Data Loading","metadata":{}},{"cell_type":"code","source":"data_folder <- '/kaggle/input/binned-dataset-v3/'\nauxiliary_folder <- '/kaggle/input/my-data/'\n\n# Load numpy arrays using reticulate\nnp <- import(\"numpy\")\ndata_train <- np$load(paste0(data_folder, 'data_train.npy'))\ndata_train_FGS <- np$load(paste0(data_folder, 'data_train_FGS.npy'))\n\n# Load train labels\ntrain_solution <- read.csv(paste0(auxiliary_folder, 'train_labels.csv'))\ntargets <- as.matrix(train_solution[, -1])  # Remove first column (IDs)\ntargets_mean <- rowMeans(targets[, -1])  # Exclude FGS column\nN <- nrow(targets)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Data Preprocessing","metadata":{}},{"cell_type":"code","source":"# Process the data (simplified - adjust based on actual dimensions)\nsignal_AIRS <- data_train\nsignal_FGS <- data_train_FGS\n\n# Sum along appropriate dimensions\nFGS_column <- apply(signal_FGS, c(1, 2), sum)\ndataset <- apply(signal_AIRS, c(1, 2, 3), sum)\n\n# Normalize by star spectrum\nnorm_star_spectrum <- function(signal) {\n  n_samples <- dim(signal)[1]\n  n_time <- dim(signal)[2]\n  n_wave <- dim(signal)[3]\n  \n  img_star <- apply(signal[, 1:50, , drop = FALSE], c(1, 3), mean) + \n              apply(signal[, (n_time-49):n_time, , drop = FALSE], c(1, 3), mean)\n  \n  result <- array(0, dim = dim(signal))\n  for (i in 1:n_samples) {\n    for (j in 1:n_time) {\n      result[i, j, ] <- signal[i, j, ] / img_star[i, ]\n    }\n  }\n  return(result)\n}\n\ndataset_norm <- norm_star_spectrum(dataset)\ndataset_norm <- aperm(dataset_norm, c(1, 3, 2))  # Transpose\ncat(\"Normalized dataset shape:\", dim(dataset_norm), \"\\n\\n\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### TRAIN/VALIDATION SPLIT","metadata":{}},{"cell_type":"code","source":"split_data <- function(data, targets, N_train) {\n  n_total <- dim(data)[1]\n  \n  cat(\"Total samples:\", n_total, \"\\n\")\n  cat(\"Requested N_train:\", N_train, \"\\n\")\n  \n  if (N_train > n_total) {\n    warning(\"N_train exceeds available samples. Using \", floor(0.8 * n_total))\n    N_train <- floor(0.8 * n_total)\n  }\n  \n  train_indices <- sample(1:n_total, N_train)\n  \n  list(\n    train_data = data[train_indices, , , drop = FALSE],\n    valid_data = data[-train_indices, , , drop = FALSE],\n    train_targets = targets[train_indices, , drop = FALSE],\n    valid_targets = targets[-train_indices, , drop = FALSE],\n    train_idx = train_indices\n  )\n}\n\nN_train <- floor(8 * dim(dataset_norm)[1] / 10)\ncat(\"Calculated N_train:\", N_train, \"\\n\")\n\nsplit_result <- split_data(dataset_norm, targets, N_train)\n\ntrain_obs <- split_result$train_data\nvalid_obs <- split_result$valid_data\ntrain_targets <- split_result$train_targets\nvalid_targets <- split_result$valid_targets\ntrain_idx <- split_result$train_idx\n\ncat(\"\\n\")\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### WHITE CURVE PROCESSING","metadata":{}},{"cell_type":"code","source":"cat(\"Processing white light curves...\\n\")\n\nsignal_AIRS_binned <- apply(signal_AIRS, c(1, 2, 3), sum)\n\n# Calculate white curve properly\nn_samples <- dim(signal_AIRS_binned)[1]\nn_time <- dim(signal_AIRS_binned)[2]\nn_wave <- dim(signal_AIRS_binned)[3]\n\nwc_mean <- numeric(n_samples)\nwhite_curve <- matrix(0, nrow = n_samples, ncol = n_time)\n\nfor (i in 1:n_samples) {\n  wc_mean[i] <- mean(signal_AIRS_binned[i, , ])\n  white_curve[i, ] <- rowSums(signal_AIRS_binned[i, , ])\n}\nwhite_curve <- white_curve / wc_mean\n\ncat(\"White curve shape:\", dim(white_curve), \"\\n\")\n\n# Split white curve data\ntrain_wc <- white_curve[train_idx, , drop = FALSE]\nvalid_wc <- white_curve[-train_idx, , drop = FALSE]\ntrain_targets_wc <- targets_mean[train_idx]\nvalid_targets_wc <- targets_mean[-train_idx]\n\ncat(\"train_wc shape:\", dim(train_wc), \"\\n\")\ncat(\"train_targets_wc length:\", length(train_targets_wc), \"\\n\")\n\n# Normalize white curve\nwlc_min <- min(train_wc)\nwlc_max <- max(train_wc)\ntrain_wc <- (train_wc - wlc_min) / (wlc_max - wlc_min)\nvalid_wc <- (valid_wc - wlc_min) / (wlc_max - wlc_min)\n\n# Normalize targets\nmin_train_valid_wc <- min(train_targets_wc)\nmax_train_valid_wc <- max(train_targets_wc)\ntrain_targets_wc_norm <- (train_targets_wc - min_train_valid_wc) / (max_train_valid_wc - min_train_valid_wc)\nvalid_targets_wc_norm <- (valid_targets_wc - min_train_valid_wc) / (max_train_valid_wc - min_train_valid_wc)\n\n# Reshape for Keras\ntrain_wc <- array(train_wc, dim = c(nrow(train_wc), ncol(train_wc), 1))\nvalid_wc <- array(valid_wc, dim = c(nrow(valid_wc), ncol(valid_wc), 1))\n\ncat(\"Reshaped train_wc:\", dim(train_wc), \"\\n\")\ncat(\"train_targets_wc_norm length:\", length(train_targets_wc_norm), \"\\n\\n\")\n\n# Normalize targets\nmin_train_valid_wc <- min(train_targets_wc)\nmax_train_valid_wc <- max(train_targets_wc)\ntrain_targets_wc_norm <- (train_targets_wc - min_train_valid_wc) / (max_train_valid_wc - min_train_valid_wc)\nvalid_targets_wc_norm <- (valid_targets_wc - min_train_valid_wc) / (max_train_valid_wc - min_train_valid_wc)\n\n# ADD THESE TWO LINES:\ntrain_targets_wc_norm <- as.array(train_targets_wc_norm)\nvalid_targets_wc_norm <- as.array(valid_targets_wc_norm)\n\n# Reshape for Keras\ntrain_wc <- array(train_wc, dim = c(nrow(train_wc), ncol(train_wc), 1))\nvalid_wc <- array(valid_wc, dim = c(nrow(valid_wc), ncol(valid_wc), 1))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 1D CNN MODEL FOR WHITE CURVE","metadata":{}},{"cell_type":"code","source":"py_require(\"keras\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cat(\"Building 1D CNN model...\\n\")\n\n# FIXED: Use functional API properly for keras3\ninput_wc <- layer_input(shape = c(ncol(white_curve), 1))\n\noutput_wc <- input_wc %>%\n  layer_conv_1d(filters = 32, kernel_size = 3, activation = 'relu') %>%\n  layer_max_pooling_1d(pool_size = 2) %>%\n  layer_batch_normalization() %>%\n  layer_conv_1d(filters = 64, kernel_size = 3, activation = 'relu') %>%\n  layer_max_pooling_1d(pool_size = 2) %>%\n  layer_conv_1d(filters = 128, kernel_size = 3, activation = 'relu') %>%\n  layer_max_pooling_1d(pool_size = 2) %>%\n  layer_conv_1d(filters = 256, kernel_size = 3, activation = 'relu') %>%\n  layer_max_pooling_1d(pool_size = 2) %>%\n  layer_flatten() %>%\n  layer_dense(units = 500, activation = 'relu') %>%\n  layer_dropout(rate = 0.2) %>%\n  layer_dense(units = 100, activation = 'relu') %>%\n  layer_dropout(rate = 0.1) %>%\n  layer_dense(units = 1, activation = 'linear')\n\nmodel_wc <- keras_model(inputs = input_wc, outputs = output_wc)\n\n# FIXED: Proper compile syntax for keras3\nmodel_wc$compile(\n  optimizer = optimizer_sgd(learning_rate = 0.001),\n  loss = 'mse',\n  metrics = list('mae')\n)\n\nmodel_wc$summary()\n\n# Learning rate scheduler\nlr_schedule_callback <- callback_learning_rate_scheduler(\n  schedule = function(epoch, lr) {\n    decay_rate <- 0.2\n    decay_step <- 200\n    if (epoch %% decay_step == 0 && epoch > 0) {\n      return(lr * decay_rate)\n    }\n    return(lr)\n  }\n)\n\n# Checkpoint\ncheckpoint_callback <- callback_model_checkpoint(\n  filepath = 'output/model_1dcnn.keras',\n  monitor = 'val_loss',\n  save_best_only = TRUE,\n  mode = 'min',\n  verbose = 1\n)\n\n# Create output directory\nif (!dir.exists('output')) {\n  dir.create('output')\n}\n\n# RIGHT BEFORE model_wc$fit(), add these:\ncat(\"DEBUG INFO:\\n\")\ncat(\"train_wc dimensions:\", dim(train_wc), \"\\n\")\ncat(\"train_targets_wc_norm length:\", length(train_targets_wc_norm), \"\\n\")\ncat(\"valid_wc dimensions:\", dim(valid_wc), \"\\n\")\ncat(\"valid_targets_wc_norm length:\", length(valid_targets_wc_norm), \"\\n\")\ncat(\"train_idx length:\", length(train_idx), \"\\n\")\ncat(\"white_curve dimensions:\", dim(white_curve), \"\\n\")\n\n# Train model\ncat('\\nTraining 1D CNN model...\\n')\nhistory_wc <- model_wc$fit(\n  x = train_wc,\n  y = train_targets_wc_norm,\n  validation_data = list(valid_wc, valid_targets_wc_norm),\n  batch_size = 16L,\n  epochs = 1200L,\n  verbose = 1L,\n  callbacks = list(checkpoint_callback, lr_schedule_callback)\n)\n\n# Load best model\nmodel_wc <- load_model('output/model_1dcnn.keras')\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### MC Dropout for uncertainty","metadata":{}},{"cell_type":"code","source":"MC_dropout_WC <- function(model, data, nb_dropout = 1000) {\n  predictions <- matrix(0, nrow = nb_dropout, ncol = nrow(data))\n  \n  cat('Running MC Dropout with', nb_dropout, 'iterations...\\n')\n  for (i in 1:nb_dropout) {\n    if (i %% 100 == 0) cat('Iteration', i, '/', nb_dropout, '\\n')\n    # In Keras 3, need to call model directly with training=TRUE for dropout\n    pred <- model(data, training = TRUE)\n    predictions[i, ] <- as.vector(pred$numpy())\n  }\n  \n  return(predictions)\n}\n\nunstandardize <- function(data, min_val, max_val) {\n  return(data * (max_val - min_val) + min_val)\n}\n\n# Run MC Dropout\ndo_mcdropout_wc <- TRUE\nnb_dropout_wc <- 1000\n\nif (do_mcdropout_wc) {\n  prediction_valid_wc <- MC_dropout_WC(model_wc, valid_wc, nb_dropout_wc)\n  spectre_valid_wc_all <- unstandardize(prediction_valid_wc, min_train_valid_wc, max_train_valid_wc)\n  spectre_valid_wc <- colMeans(spectre_valid_wc_all)\n  spectre_valid_std_wc <- apply(spectre_valid_wc_all, 2, sd)\n} else {\n  pred <- model_wc(valid_wc, training = FALSE)\n  spectre_valid_wc <- as.vector(pred$numpy())\n  spectre_valid_wc <- unstandardize(spectre_valid_wc, min_train_valid_wc, max_train_valid_wc)\n  spectre_valid_std_wc <- 0.1 * abs(spectre_valid_wc)\n}\n\n# Calculate residuals\nresiduals <- spectre_valid_wc - valid_targets_wc\nmse_ppm <- sqrt(mean(residuals^2)) * 1e6\ncat('White Curve MSE:', mse_ppm, 'ppm\\n')\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 2D CNN MODEL preprocessing","metadata":{}},{"cell_type":"code","source":"cat(\"Preparing data for 2D CNN...\\n\")\n\nsuppress_mean <- function(targets, mean_vals) {\n  targets - matrix(mean_vals, nrow = length(mean_vals), ncol = ncol(targets))\n}\n\ntrain_targets_shift <- suppress_mean(train_targets, targets_mean[train_idx])\nvalid_targets_shift <- suppress_mean(valid_targets, targets_mean[-train_idx])\n\n# Normalize\ndata_min <- min(train_targets_shift)\ndata_max <- max(train_targets_shift)\ntargets_abs_max <- max(abs(c(data_min, data_max)))\ntrain_targets_norm <- train_targets_shift / targets_abs_max\nvalid_targets_norm <- valid_targets_shift / targets_abs_max\n\n# Extract in-transit\ningress <- 75\negress <- 115\n\n# Adjust if needed\nif (egress > dim(train_obs)[2]) {\n  ingress <- max(1, floor(dim(train_obs)[2] * 0.4))\n  egress <- min(dim(train_obs)[2], floor(dim(train_obs)[2] * 0.6))\n  cat(\"Adjusted ingress/egress to:\", ingress, \"/\", egress, \"\\n\")\n}\n\ntrain_obs_in <- train_obs[, ingress:egress, , drop = FALSE]\nvalid_obs_in <- valid_obs[, ingress:egress, , drop = FALSE]\n\n# Subtract mean\nsubstract_data_mean <- function(data) {\n  result <- array(0, dim = dim(data))\n  for (i in 1:dim(data)[1]) {\n    result[i, , ] <- data[i, , ] - mean(data[i, , ])\n  }\n  return(result)\n}\n\ntrain_obs_2d_mean <- substract_data_mean(train_obs_in)\nvalid_obs_2d_mean <- substract_data_mean(valid_obs_in)\n\n# Normalize\ndata_min <- min(train_obs_2d_mean)\ndata_max <- max(train_obs_2d_mean)\ndata_abs_max <- max(abs(c(data_min, data_max)))\ntrain_obs_norm <- train_obs_2d_mean / data_abs_max\nvalid_obs_norm <- valid_obs_2d_mean / data_abs_max\n\n# Add channel dimension\ntrain_obs_norm <- array(train_obs_norm, dim = c(dim(train_obs_norm), 1))\nvalid_obs_norm <- array(valid_obs_norm, dim = c(dim(valid_obs_norm), 1))\n\ncat(\"2D CNN input shape:\", dim(train_obs_norm), \"\\n\\n\")\n\n# Normalize\ndata_min <- min(train_targets_shift)\ndata_max <- max(train_targets_shift)\ntargets_abs_max <- max(abs(c(data_min, data_max)))\ntrain_targets_norm <- train_targets_shift / targets_abs_max\nvalid_targets_norm <- valid_targets_shift / targets_abs_max\n\n# ADD THESE TWO LINES:\ntrain_targets_norm <- as.array(train_targets_norm)\nvalid_targets_norm <- as.array(valid_targets_norm)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"2d Model","metadata":{}},{"cell_type":"code","source":"cat(\"Building 2D CNN model...\\n\")\n\ninput_obs <- layer_input(shape = dim(train_obs_norm)[-1])\n\noutput_2d <- input_obs %>%\n  layer_conv_2d(filters = 32, kernel_size = c(3, 1), activation = 'relu', padding = 'same') %>%\n  layer_max_pooling_2d(pool_size = c(2, 1)) %>%\n  layer_batch_normalization() %>%\n  layer_conv_2d(filters = 64, kernel_size = c(3, 1), activation = 'relu', padding = 'same') %>%\n  layer_max_pooling_2d(pool_size = c(2, 1)) %>%\n  layer_conv_2d(filters = 128, kernel_size = c(3, 1), activation = 'relu', padding = 'same') %>%\n  layer_max_pooling_2d(pool_size = c(2, 1)) %>%\n  layer_conv_2d(filters = 256, kernel_size = c(3, 1), activation = 'relu', padding = 'same') %>%\n  layer_conv_2d(filters = 32, kernel_size = c(1, 3), activation = 'relu', padding = 'same') %>%\n  layer_max_pooling_2d(pool_size = c(1, 2)) %>%\n  layer_batch_normalization() %>%\n  layer_conv_2d(filters = 64, kernel_size = c(1, 3), activation = 'relu', padding = 'same') %>%\n  layer_max_pooling_2d(pool_size = c(1, 2)) %>%\n  layer_conv_2d(filters = 128, kernel_size = c(1, 3), activation = 'relu', padding = 'same') %>%\n  layer_max_pooling_2d(pool_size = c(1, 2)) %>%\n  layer_conv_2d(filters = 256, kernel_size = c(1, 3), activation = 'relu', padding = 'same') %>%\n  layer_max_pooling_2d(pool_size = c(1, 2)) %>%\n  layer_flatten() %>%\n  layer_dense(units = 700, activation = 'relu') %>%\n  layer_dropout(rate = 0.2) %>%\n  layer_dense(units = ncol(train_targets), activation = 'linear')\n\nmodel_2d <- keras_model(inputs = input_obs, outputs = output_2d)\n\n# FIXED: Use $compile for keras3\nmodel_2d$compile(\n  optimizer = optimizer_adam(learning_rate = 0.001),\n  loss = 'mse',\n  metrics = list('mae')\n)\n\nmodel_2d$summary()\n\n# Checkpoint\ncheckpoint_callback_2d <- callback_model_checkpoint(\n  filepath = 'output/model_2dcnn.keras',\n  monitor = 'val_loss',\n  save_best_only = TRUE,\n  mode = 'min',\n  verbose = 1\n)\n\n# Train\ncat('\\nTraining 2D CNN model...\\n')\nhistory_2d <- model_2d$fit(\n  x = train_obs_norm,\n  y = train_targets_norm,\n  validation_data = list(valid_obs_norm, valid_targets_norm),\n  batch_size = 32L,\n  epochs = 200L,\n  verbose = 1L,\n  callbacks = list(checkpoint_callback_2d)\n)\n\nmodel_2d <- load_model('output/model_2dcnn.keras')\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### predictions","metadata":{}},{"cell_type":"code","source":"cat(\"\\nRunning MC Dropout for 2D CNN...\\n\")\n\nNN_uncertainty <- function(model, x_test, targets_abs_max, T = 5) {\n  n_samples <- dim(x_test)[1]\n  n_features <- ncol(train_targets)\n  predictions <- array(0, dim = c(T, n_samples, n_features))\n  \n  for (i in 1:T) {\n    cat('Iteration', i, '/', T, '\\n')\n    pred_norm <- model(x_test, training = TRUE)\n    pred <- as.array(pred_norm) * targets_abs_max\n    predictions[i, , ] <- pred\n  }\n  \n  list(mean = apply(predictions, c(2, 3), mean),\n       std = apply(predictions, c(2, 3), sd))\n}\n\ndo_mcdropout <- TRUE\nnb_dropout <- 5\n\nif (do_mcdropout) {\n  result_2d <- NN_uncertainty(model_2d, valid_obs_norm, targets_abs_max, T = nb_dropout)\n  spectre_valid_shift <- result_2d$mean\n  spectre_valid_shift_std <- result_2d$std\n} else {\n  pred_valid_norm <- model_2d(valid_obs_norm, training = FALSE)\n  spectre_valid_shift <- as.array(pred_valid_norm) * targets_abs_max\n  spectre_valid_shift_std <- abs(spectre_valid_shift) * 0.1\n}\n\n\npredictions_valid <- spectre_valid_shift + matrix(spectre_valid_wc, \n                                                   nrow = length(spectre_valid_wc), \n                                                   ncol = ncol(spectre_valid_shift))\n\npredictions_std_valid <- sqrt(matrix(spectre_valid_std_wc^2, \n                                     nrow = length(spectre_valid_std_wc), \n                                     ncol = ncol(spectre_valid_shift_std)) + \n                              spectre_valid_shift_std^2)\n\nresiduals_final <- valid_targets - predictions_valid\nmse_final_ppm <- sqrt(mean(residuals_final^2)) * 1e6\ncat('Final MSE:', mse_final_ppm, 'ppm\\n\\n')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### PLOTTING FUNCTIONS","metadata":{}},{"cell_type":"code","source":"wavelengths <- read.csv(paste0(auxiliary_folder, 'wavelengths.csv'))\nwl <- wavelengths[, 2]\n\nplot_sample <- function(sample_idx, predictions, targets, std, wavelengths) {\n  df <- data.frame(\n    wavelength = wavelengths,\n    prediction = predictions[sample_idx, ],\n    target = targets[sample_idx, ],\n    lower = predictions[sample_idx, ] - std[sample_idx, ],\n    upper = predictions[sample_idx, ] + std[sample_idx, ]\n  )\n  \n  p <- ggplot(df, aes(x = wavelength)) +\n    geom_ribbon(aes(ymin = lower, ymax = upper), fill = 'gray', alpha = 0.5) +\n    geom_point(aes(y = prediction, color = 'Prediction'), size = 1) +\n    geom_line(aes(y = target, color = 'Target'), linewidth = 0.8) +\n    scale_color_manual(values = c('Prediction' = 'black', 'Target' = 'tomato')) +\n    labs(title = paste('Sample', sample_idx),\n         x = 'Wavelength (μm)', y = '(Rp/Rs)²', color = '') +\n    theme_minimal() +\n    theme(legend.position = 'bottom')\n  \n  print(p)\n}\n\n# Plot samples\nfor (i in 1:min(3, nrow(predictions_valid))) {\n  plot_sample(i, predictions_valid, valid_targets, predictions_std_valid, wl)\n}\n\ncat(\"\\n=== COMPLETE ===\\n\")\ncat(\"Final MSE:\", round(mse_final_ppm, 2), \"ppm\\n\")","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}