library(dplyr) basis <- expand.grid(c3 = 0:3, c5 = 0:5) basis$id <- 1:nrow(basis) fill_3 <- function (c3, c5) { list(c3 = 3, c5 = c5) } fill_5 <- function (c3, c5) { list(c3 = c3, c5 = 5) } dump_3 <- function (c3, c5) { list(c3 = 0, c5 = c5) } dump_5 <- function (c3, c5) { list(c3 = c3, c5 = 0) } from_3_to_5 <- function (c3, c5) { transfer <- min(5 - c5, c3) list(c3 = c3 - transfer, c5 = c5 + transfer) } from_5_to_3 <- function (c3, c5) { transfer <- min(3 - c3, c5) list(c3 = c3 + transfer, c5 = c5 - transfer) } functions <- list(fill_3, fill_5, dump_3, dump_5, from_3_to_5, from_5_to_3) function_to_matrix <- function (f) { m <- matrix(0, nrow = nrow(basis), ncol = nrow(basis)) for (i in 1:nrow(basis)) { elem_in <- basis[i, ] res <- f(elem_in$c3, elem_in$c5) j <- filter(basis, c3 == res$c3, c5 == res$c5)$id m[j, i] <- 1 } return (m) } d3 <- function_to_matrix(dump_3) d5 <- function_to_matrix(dump_5) f3 <- function_to_matrix(fill_3) f5 <- function_to_matrix(fill_5) t35 <- function_to_matrix(from_5_to_3) t53 <- function_to_matrix(from_3_to_5) stopifnot(all(d3 %*% d3 == d3)) stopifnot(all(d5 %*% d5 == d5)) stopifnot(all(d5 %*% d3 == d3 %*% d5)) vs <- numeric(nrow(basis)) vs[1] <- 1 v1 <- (basis$c3 == 1 | basis$c5 == 1) * 1 v2 <- (basis$c3 == 2 | basis$c5 == 2) * 1 v3 <- (basis$c3 == 3 | basis$c5 == 3) * 1 v4 <- (basis$c3 == 4 | basis$c5 == 4) * 1 operators <- list(d3 = d3, d5 = d5, f3 = f3, f5 = f5, t53 = t35, t35 = t53) orbits <- list(list(ops = c(), vec = vs)) for (iter in 1:100) { for (op_id in 1:length(operators)) { op_name <- names(operators)[op_id] op_matrix <- operators[[op_id]] for (or_id in 1:length(orbits)) { or_ops <- orbits[[or_id]]$ops or_vec <- orbits[[or_id]]$vec new_vec <- op_matrix %*% or_vec found <- FALSE for (or2 in orbits) { if (all(new_vec == or2$vec)) { found <- TRUE break } } if (!found) { orbits <- c(orbits, list(list(ops = c(or_ops, op_name), vec = new_vec))) break } } } } endings <- lapply(orbits, function (orbit) { b <- basis[c(orbit$vec == 1), ] stopifnot(nrow(b) == 1) list(c3 = b$c3, c5 = b$c5, operations = paste(orbit$ops, collapse = ' ')) }) endings <- do.call(rbind, lapply(endings, as.data.frame)) %>% mutate(Sum = c3 + c5) arranged <- endings %>% arrange(c3, c5) print(arranged) u <- unique(c(endings$c3, endings$c5)) us <- unique(endings$Sum) print(u[order(u)]) print(us[order(us)]) edges <- list() for (i in 1:nrow(arranged)) { r <- arranged[i, ] cat(sprintf('"%d, %d" [fillcolor = "#ffaa7f"]\n', r$c3, r$c5)) } for (i in 1:nrow(basis)) { before <- basis[i, ] for (f in functions) { after <- f(before$c3, before$c5) edges <- c(edges, list(list(sprintf('%d, %d', before$c3, before$c5), sprintf('%d, %d', after$c3, after$c5)))) } } for (edge in edges) { if (edge[[1]] != edge[[2]]) { cat(sprintf('"%s" -> "%s"\n', edge[[1]], edge[[2]])) } } matrix_sum = (Reduce(`+`, operators, matrix(0, nrow = 24, ncol = 24)) != 0) * 1 print(matrix_sum) pm <- paste(apply(matrix_sum, 1, paste, collapse = ' '), collapse = '\n') cat(pm, '\n') library(igraph) g <- graph_from_adjacency_matrix(matrix_sum) reachable <- sapply(1:nrow(basis), function (to) { length(all_shortest_paths(g, 1, to, mode = 'in')$res) })