Feature/add met example

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",

R/design_utils.R

@@ -52,7 +52,7 @@ 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 (length(block_indices) >= 2) {
       # Need at least 2 plots to swap
@@ -111,7 +111,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 +339,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)
@@ -363,6 +364,20 @@ 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)) {
@@ -397,6 +412,10 @@ random_initialise <- function(design, optimise, seed = NULL, ...) {
     }
   }
 
+  for (opt in optimise[-1]) {
+    best_design[[opt$swap_within]] <- NULL
+  }
+
   return(best_design)
 }
 

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,

vignettes/met.qmd

@@ -113,12 +113,17 @@ 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}
+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)))
+check_no_dupes(df)
 ```
 
 ## Visualise the Output
@@ -163,7 +168,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)
 ```
@@ -196,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
 )
 
@@ -215,12 +221,11 @@ 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)))
+check_no_dupes(df)
 ```
 
 Site "a" maintains 3 replicates and no missing treatments in any other sites.
@@ -249,6 +254,145 @@ 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
+
+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,
+  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)
+)
+
+met_result <- speed(
+  data = met_design,
+  swap = "treatment",
+  early_stop_iterations = 5000,
+  iterations = 20000,
+  optimise = optimise,
+  optimise_params = optim_params(random_initialisation = 30, 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
+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)
+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