From 44e3079696b2816ca10d1c8203b10d8afa8c4bd2 Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Mon, 23 Mar 2026 15:36:52 +1030 Subject: [PATCH 1/4] debugging --- DESCRIPTION | 2 +- R/design_utils.R | 20 +++++-- R/speed.R | 5 +- vignettes/met.qmd | 146 +++++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 162 insertions(+), 11 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2959766..542fbb3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: speed Title: Generate Spatially Efficient Experimental Designs -Version: 0.0.4 +Version: 0.0.5 Authors@R: c(person(given = "Sam", family = "Rogers", diff --git a/R/design_utils.R b/R/design_utils.R index fe2e494..d2bfe11 100644 --- a/R/design_utils.R +++ b/R/design_utils.R @@ -52,7 +52,12 @@ generate_single_swap_neighbour <- function(design, swap, swap_within, swap_count # Perform swaps in selected blocks for (block in blocks_to_swap) { # Get indices of plots in this block - block_indices <- which(all_blocks == block & !is.na(design[[swap]])) + block_indices <- which(all_blocks == block & !is.na(all_blocks) & !is.na(design[[swap]])) + # if (block == 'a') { + # + # print(block) + # print(block_indices) + # } if (length(block_indices) >= 2) { # Need at least 2 plots to swap @@ -75,6 +80,10 @@ generate_single_swap_neighbour <- function(design, swap, swap_within, swap_count } } + # if (block == 'a' && '9' %in% to_be_swapped) { + # print(paste(block, paste(swap_pair, collapse = ','), paste(to_be_swapped, collapse = ','))) + # } + # Perform the swap only if we have valid treatments to swap if (!is.null(to_be_swapped)) { new_design[[swap]][rev(swap_pair)] <- to_be_swapped @@ -111,7 +120,7 @@ generate_multi_swap_neighbour <- function(design, swap, swap_within, swap_count, # Perform swaps in selected groups for (group in groups_to_swap) { # Get unique treatments within this group - group_filter <- new_design[[swap_within]] == group + group_filter <- new_design[[swap_within]] == group & !is.na(new_design[[swap_within]]) group_data <- new_design[group_filter & !is.na(new_design[[swap]]), ] group_treatments <- unique(group_data[[swap]]) @@ -339,9 +348,10 @@ shuffle_items <- function(design, swap, swap_within, seed = NULL) { set.seed(seed) } - for (i in unique(design[[swap_within]])) { - items <- design[design[[swap_within]] == i, ][[swap]] - design[design[[swap_within]] == i, ][[swap]] <- sample(items) + for (i in levels(design[[swap_within]])) { + swap_within_filter <- design[[swap_within]] == i & !is.na(design[[swap_within]]) + items <- design[swap_within_filter, ][[swap]] + design[swap_within_filter, ][[swap]] <- sample(items) } return(design) diff --git a/R/speed.R b/R/speed.R index e9afe4b..b4f81b7 100644 --- a/R/speed.R +++ b/R/speed.R @@ -188,6 +188,7 @@ speed <- function(data, if (inferred$inferred) { # Sort the data frame to start with to ensure consistency in calculating the adjacency later data <- data[do.call(order, data[c(row_column, col_column)]), ] + rownames(data) <- seq_len(nrow(data)) } # dummy group for swapping within whole design @@ -281,8 +282,8 @@ speed_hierarchical <- function(data, optimise, quiet, seed, ...) { } # Generate new design by swapping treatments at this level - new_design <- generate_neighbour(current_design,opt$swap, opt$swap_within, current_swap_count, - current_swap_all_blocks,opt$swap_all) + new_design <- generate_neighbour(current_design, opt$swap, opt$swap_within, current_swap_count, + current_swap_all_blocks, opt$swap_all) # Calculate new score new_score_obj <- opt$obj_function(new_design$design,opt$swap, spatial_cols, adj_weight = adj_weight, diff --git a/vignettes/met.qmd b/vignettes/met.qmd index c18a93a..d86d95b 100644 --- a/vignettes/met.qmd +++ b/vignettes/met.qmd @@ -113,12 +113,13 @@ assess the quality of optimisation at each hierarchy level. str(met_result) ``` -No duplicated treatments along any row or column. +No duplicated treatments along any row, column, or block. ```{r met-no-dupes} df <- met_result$design_df max(table(df$treatment, df$site_col)) max(table(df$treatment, paste0(df$site, df$row))) +max(table(df$treatment, df$site_block)) ``` ## Visualise the Output @@ -163,7 +164,7 @@ met_design <- initialise_design_df( met_design$site_col <- paste(met_design$site, met_design$col, sep = "_") met_design$site_block <- paste(met_design$site, met_design$block, sep = "_") met_design$allocation <- "free" -met_design$allocation[met_design$site == "a"] <- "fixed" +met_design$allocation[met_design$site == "a"] <- NA head(met_design) ``` @@ -215,12 +216,13 @@ assess the quality of optimisation at each hierarchy level. str(met_result) ``` -No duplicated treatments along any row or column. +No duplicated treatments along any row, column, or block. ```{r met2-no-dupes} df <- met_result$design_df max(table(df$treatment, df$site_col)) max(table(df$treatment, paste0(df$site, df$row))) +max(table(df$treatment, df$site_block)) ``` Site "a" maintains 3 replicates and no missing treatments in any other sites. @@ -249,6 +251,144 @@ ggplot2::ggplot(met_result$design_df, ggplot2::aes(col, row, fill = treatment)) ``` This design has now been optimised at both the connectivity between sites and the balance within each site. + +# Optimising Partial Allocation across Sites + +## Setting Up MET Design with speed + +Now we can create a data frame representing a MET design. Note that we can specify different dimensions for each +site in the `designs` argument. Also, some treatments are pre-allocated to each site: + +- Site 'a': Treatments 1-5 +- Site 'b': Treatments 1-4 +- Site 'c': Treatments 1-7 +- Site 'd': Treatments 1-3, 8-9 +- Site 'e': Treatments 1-5 + +```{r met3-df} +fixed_treatments <- list( + a = 1:5, + b = 1:4, + c = 1:7, + d = c(1:3, 8:9), + e = 1:5 +) +non_fixed_treatments <- c(rep(10:18, 8), rep(19:60, 7)) +met_design <- initialise_design_df( + items = 1, + designs = list( + a = list(nrows = 27, ncols = 6, block_nrows = 9, block_ncols = 6), + b = list(nrows = 20, ncols = 3, block_nrows = 10, block_ncols = 3), + c = list(nrows = 20, ncols = 4, block_nrows = 10, block_ncols = 4), + d = list(nrows = 15, ncols = 4, block_nrows = 5, block_ncols = 4), + e = list(nrows = 22, ncols = 3, block_nrows = 11, block_ncols = 3) + ) +) +met_design$allocation <- "free" +for (site in unique(met_design$site)) { + df_site <- met_design[met_design$site == site, ] + n_blocks <- max(unique(df_site$block)) + n_fixed <- length(fixed_treatments[[site]]) * n_blocks + non_fixed_indices <- 1:(nrow(df_site) - n_fixed) + + treatments <- c(rep(fixed_treatments[[site]], n_blocks), non_fixed_treatments[non_fixed_indices]) + met_design[met_design$site == site, ]$treatment <- treatments + met_design[met_design$site == site & met_design$treatment %in% fixed_treatments[[site]], ]$allocation <- NA + + non_fixed_treatments <- non_fixed_treatments[-non_fixed_indices] +} + +met_design$site_col <- paste(met_design$site, met_design$col, sep = "_") +met_design$site_block <- paste(met_design$site, met_design$block, sep = "_") + +head(met_design) +``` + +```{r met3-plot1} +#| echo: false +met_design$block <- factor(met_design$block) + +ggplot2::ggplot(met_design, ggplot2::aes(col, row, fill = block)) + + ggplot2::geom_tile(color = "black") + + ggplot2::scale_fill_viridis_d() + + ggplot2::facet_wrap(~site, scales = "free") + + ggplot2::scale_x_continuous(expand = c(0, 0), breaks = 1:max(met_design$col)) + + ggplot2::scale_y_continuous(expand = c(0, 0), breaks = 1:max(met_design$row), trans = scales::reverse_trans()) + + ggplot2::theme_bw() + + ggplot2::theme(panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank()) +``` + +## Performing the Optimisation + +For MET designs, we use lists of named arguments to specify the hierarchical structure. The `optimise` parameter +defines what to optimise and constraints at each level. + +```{r met3-example} +optimise <- list( + connectivity = list(swap_within = "allocation", spatial_factors = ~site), + balance = list(swap_within = "site", spatial_factors = ~ site_col + site_block, iterations = 50) +) + +met_result <- speed( + data = met_design, + swap = "treatment", + early_stop_iterations = 10000, + iterations = 30000, + optimise = optimise, + optimise_params = optim_params(random_initialisation = TRUE, adj_weight = 0), + seed = 112 +) + +met_result +``` + +## Output of the Optimisation + +The output shows optimisation results for the design. The score and iterations are combined for the entire +design, while the treatments, and stopping criteria are reported separately for each level, allowing you to +assess the quality of optimisation at each hierarchy level. + +```{r met3-output} +str(met_result) +``` + +No duplicated treatments along any row, column, or block. + +```{r met3-no-dupes} +df <- met_result$design_df +max(table(df$treatment, df$site_col)) +max(table(df$treatment, paste0(df$site, df$row))) +max(table(df$treatment, df$site_block)) +``` + + + + + + + + + +## Visualise the Output + +```{r met3-plot2} +met_result$design_df$treatment = as.numeric(met_result$design_df$treatment) +ggplot2::ggplot(met_result$design_df, ggplot2::aes(col, row, fill = treatment)) + + ggplot2::geom_tile(color = "black") + + ggplot2::scale_fill_viridis_c() + + ggplot2::facet_wrap(~site, scales = "free") + + ggplot2::scale_x_continuous(expand = c(0, 0), breaks = 1:max(met_design$col)) + + ggplot2::scale_y_continuous(expand = c(0, 0), breaks = 1:max(met_design$row), trans = scales::reverse_trans()) + + ggplot2::theme_bw() + + ggplot2::theme( + legend.position = "none", + panel.grid.major = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank() + ) +``` + +This design has now been optimised at both the connectivity between sites and the balance within each site. + # Spatial Design Considerations ## Field Shape and Orientation From 76cbdf00ebbfa7f790a741ddaae82e2e88cf9991 Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Tue, 24 Mar 2026 18:24:22 +1030 Subject: [PATCH 2/4] fix shuffle bug for complex designs --- R/design_utils.R | 29 +++++++++++++++++++--------- vignettes/met.qmd | 48 ++++++++++++++++++++++++----------------------- 2 files changed, 45 insertions(+), 32 deletions(-) diff --git a/R/design_utils.R b/R/design_utils.R index d2bfe11..8306e54 100644 --- a/R/design_utils.R +++ b/R/design_utils.R @@ -53,11 +53,6 @@ generate_single_swap_neighbour <- function(design, swap, swap_within, swap_count for (block in blocks_to_swap) { # Get indices of plots in this block block_indices <- which(all_blocks == block & !is.na(all_blocks) & !is.na(design[[swap]])) - # if (block == 'a') { - # - # print(block) - # print(block_indices) - # } if (length(block_indices) >= 2) { # Need at least 2 plots to swap @@ -80,10 +75,6 @@ generate_single_swap_neighbour <- function(design, swap, swap_within, swap_count } } - # if (block == 'a' && '9' %in% to_be_swapped) { - # print(paste(block, paste(swap_pair, collapse = ','), paste(to_be_swapped, collapse = ','))) - # } - # Perform the swap only if we have valid treatments to swap if (!is.null(to_be_swapped)) { new_design[[swap]][rev(swap_pair)] <- to_be_swapped @@ -373,6 +364,22 @@ random_initialise <- function(design, optimise, seed = NULL, ...) { return(design) } + if (length(optimise) > 1) { + groups <- c() + for (i in seq_along(optimise)) { + groups <- c(groups, optimise[[i]]$swap_within) + if (i == 1) { + next + } + + now <- as.numeric(Sys.time()) + dummy_col <- paste0(paste(groups, collapse = "_"), "_", now) + optimise[[i]]$swap_within <- dummy_col + design[[dummy_col]] <- apply(design[, groups], 1, paste, collapse = "-") |> + factor() + } + } + best_score <- Inf best_design <- design for (i in seq_len(random_initialisation)) { @@ -407,6 +414,10 @@ random_initialise <- function(design, optimise, seed = NULL, ...) { } } + for (opt in optimise[-1]) { + design[[opt$swap_within]] <- NULL + } + return(best_design) } diff --git a/vignettes/met.qmd b/vignettes/met.qmd index d86d95b..879a9c0 100644 --- a/vignettes/met.qmd +++ b/vignettes/met.qmd @@ -116,10 +116,14 @@ str(met_result) No duplicated treatments along any row, column, or block. ```{r met-no-dupes} +check_no_dupes <- function(df) { + cat(max(table(df$treatment, df$site_col)), sep = "", ", ") + cat(max(table(df$treatment, paste0(df$site, df$row))), sep = "", ", ") + cat(max(table(df$treatment, df$site_block))) +} + df <- met_result$design_df -max(table(df$treatment, df$site_col)) -max(table(df$treatment, paste0(df$site, df$row))) -max(table(df$treatment, df$site_block)) +check_no_dupes(df) ``` ## Visualise the Output @@ -197,9 +201,10 @@ optimise <- list( met_result <- speed( data = met_design, swap = "treatment", - early_stop_iterations = 5000, + early_stop_iterations = 10000, + iterations = 50000, optimise = optimise, - optimise_params = optim_params(random_initialisation = TRUE, adj_weight = 0), + optimise_params = optim_params(random_initialisation = 10, adj_weight = 0), seed = 112 ) @@ -220,9 +225,7 @@ No duplicated treatments along any row, column, or block. ```{r met2-no-dupes} df <- met_result$design_df -max(table(df$treatment, df$site_col)) -max(table(df$treatment, paste0(df$site, df$row))) -max(table(df$treatment, df$site_block)) +check_no_dupes(df) ``` Site "a" maintains 3 replicates and no missing treatments in any other sites. @@ -326,16 +329,16 @@ defines what to optimise and constraints at each level. ```{r met3-example} optimise <- list( connectivity = list(swap_within = "allocation", spatial_factors = ~site), - balance = list(swap_within = "site", spatial_factors = ~ site_col + site_block, iterations = 50) + balance = list(swap_within = "site", spatial_factors = ~ site_col + site_block) ) met_result <- speed( data = met_design, swap = "treatment", - early_stop_iterations = 10000, - iterations = 30000, + early_stop_iterations = 5000, + iterations = 20000, optimise = optimise, - optimise_params = optim_params(random_initialisation = TRUE, adj_weight = 0), + optimise_params = optim_params(random_initialisation = 30, adj_weight = 0), seed = 112 ) @@ -356,23 +359,22 @@ No duplicated treatments along any row, column, or block. ```{r met3-no-dupes} df <- met_result$design_df -max(table(df$treatment, df$site_col)) -max(table(df$treatment, paste0(df$site, df$row))) -max(table(df$treatment, df$site_block)) +check_no_dupes(df) ``` - - - - - - - +All sites maintain pre-allocated treatments. + +```{r met3-reps} +treatment_count <- table(df$site, df$treatment) +for (site in names(fixed_treatments)) { + print(treatment_count[site, as.character(fixed_treatments[[site]]), drop = FALSE]) +} +``` ## Visualise the Output ```{r met3-plot2} -met_result$design_df$treatment = as.numeric(met_result$design_df$treatment) +met_result$design_df$treatment <- as.numeric(met_result$design_df$treatment) ggplot2::ggplot(met_result$design_df, ggplot2::aes(col, row, fill = treatment)) + ggplot2::geom_tile(color = "black") + ggplot2::scale_fill_viridis_c() + From 992925ba87c03b38a73c030ed57023ad97281bbe Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Tue, 24 Mar 2026 18:26:42 +1030 Subject: [PATCH 3/4] minor change --- vignettes/met.qmd | 2 ++ 1 file changed, 2 insertions(+) diff --git a/vignettes/met.qmd b/vignettes/met.qmd index 879a9c0..025f465 100644 --- a/vignettes/met.qmd +++ b/vignettes/met.qmd @@ -268,6 +268,8 @@ site in the `designs` argument. Also, some treatments are pre-allocated to each - Site 'd': Treatments 1-3, 8-9 - Site 'e': Treatments 1-5 +Note that version 0.0.5 is required. Otherwise, there would be a bug caused by random initialisation later on. + ```{r met3-df} fixed_treatments <- list( a = 1:5, From 4d047a419784591e7edbf43903b98c209bc02c2a Mon Sep 17 00:00:00 2001 From: Wasin Pipattungsakul Date: Tue, 24 Mar 2026 18:29:51 +1030 Subject: [PATCH 4/4] fix minor bug --- R/design_utils.R | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/R/design_utils.R b/R/design_utils.R index 8306e54..9e4df47 100644 --- a/R/design_utils.R +++ b/R/design_utils.R @@ -368,9 +368,7 @@ random_initialise <- function(design, optimise, seed = NULL, ...) { groups <- c() for (i in seq_along(optimise)) { groups <- c(groups, optimise[[i]]$swap_within) - if (i == 1) { - next - } + if (i == 1) next now <- as.numeric(Sys.time()) dummy_col <- paste0(paste(groups, collapse = "_"), "_", now) @@ -415,7 +413,7 @@ random_initialise <- function(design, optimise, seed = NULL, ...) { } for (opt in optimise[-1]) { - design[[opt$swap_within]] <- NULL + best_design[[opt$swap_within]] <- NULL } return(best_design)