library(tidyverse)
library(magrittr)
library(scales)
library(data.table)
library(tictoc)
library(stringi)
library(arrangements)
library(ggforce)

drawDisk <- function(df) {
  # отрисуем расположение точек и действующие силы
  # если силы не заданы, создадим их по умолчанию равными 0
  for(col in c("force_x", "force_y")){
    if (!(col %in% names(df))) df[, col] <- 0
  }

  ggplot(data = df, aes(x = x, y = y)) +
    ggforce::geom_circle(aes(x0 = 0, y0 = 0, r = 1, colour = "red"), 
                         inherit.aes = FALSE) +
    geom_point(size = 2, fill = "black", shape = 21) +
    geom_text(aes(label = idx), hjust = 0.5, vjust = -1) +
    # рисуем векторное поле
    geom_segment(aes(xend = x + force_x / 10, yend = y + force_y / 10), 
                 colour = "red", 
                 arrow = arrow(length = unit(0.2,"cm")), size = 0.6) + 
    xlim(-1.5, 1.5) +
    ylim(-1.5, 1.5) +
    coord_fixed() +
    labs(x = "Ось X", y = "Ось Y") +
    theme_bw()
}

# Генерим точки-зонды внутри окружности единичного радиуса.
# Считаем, что все частицы единичного заряда, поэтому его опускаем
charges_dt <- tibble(idx = 1:13) %>%
  mutate(angle = runif(n(), min = 0, max = 2*pi), 
         r = runif(n(), min = 0, max = 1),
         x = r * cos(angle), y = r * sin(angle)) %>%
  select(idx, x, y) %>%
  setDT(key = "idx")

# для сходимости задачи генерируем также зафиксированные точки на внешней окружности
keepers_dt <- max(charges_dt$idx) %>% 
  {tibble(idx = (. + 1):(. + 40))} %>%
  mutate(angle = (idx - 1) * (2 * pi / n()),
         x = 1.3 * cos(angle), y = 1.3 * sin(angle)) %>%
  select(idx, x, y) %>%
  setDT(key = "idx")

full_dt <- rbindlist(list(charges_dt, keepers_dt))
drawDisk(full_dt)

max_force <- 10

tic("Balancing charges")
# Определяем точность моделирования!
# Будем двигать заряды пока они не застабилизируются 
# Максимальная равнодействуюущая будет близка к 0 
# и точки не стали разлетаться в бесконечность
while (max_force > 0.05 & nrow(charges_dt[x^2 + y^2 > 1.2]) == 0) {
  # общий пул координат частиц на текущую итерацию
  full_dt <- rbindlist(list(charges_dt, keepers_dt), fill = TRUE)

  ff <- function(x0, y0){
    # сила взаимодействия -- 1/r2, заряды единичные;
    # проекция силы на оси, sqrt(r2) -- гипотенуза
    # суперпозиция векторов даст результирующее воздействие на частицу

    dx = full_dt$x - x0
    dy = full_dt$y - y0
    r2 = dx^2 + dy^2
    # na.rm исключает NaN в т.ч.
    list(sum(-dx * r2^(-1.5), na.rm = TRUE),
         sum(-dy * r2^(-1.5), na.rm = TRUE))
  }  

  # проведем расчет сил, итерируем по каждой "неприбитой гвоздями" точке
  charges_dt[, c("force_x", "force_y") := ff(x0 = x, y0 = y), by = idx]
  # определим максимальную силу, действующую на частицу
  max_force <- charges_dt[, max(sqrt(force_x^2 + force_y^2), na.rm = TRUE)]
  force_scale = if_else(max_force > 1, 1 / max_force / 1e2, 1/ max_force / 5e2)

  # проводим передвижение точек
  charges_dt %>%
    .[, `:=`(x = x + force_x * force_scale, 
             y = y + force_y * force_scale)]
}
toc()

full_dt <- rbindlist(list(charges_dt, keepers_dt), fill = TRUE)

optimizePath <- function(dt) {
  # попробуем оптимизировать маршрут обхода по предоставленным точкам
  # 1. Выбираем в качестве начальной точки датчик, максимально близко расположенный к краю поля
  dt[, r0 := sqrt(x^2+y^2)] %>%
    setorder(-r0)
  n1 <- dt[1, idx]

  # теперь проводим симуляцию различных вариантов расстановки сенсоров
  # получаем последовательность номеров и убираем n1, его будем принудительно ставить первым
  points_in <- dt[idx != n1, idx]

  # для каждой точки добавим еще ближайшую точку выхода 
  # (перпендикуляр к окружности, которая единичного радиуса)
  augment_tbl <- dt %>%
    mutate_at("idx", `*`, -1) %>%
    mutate(r0 = sqrt(x^2 + y^2)) %>%
    mutate_at(vars(x, y), ~(.x/r0)) %>%
    bind_rows(dt) %>%
    select(idx, x, y)

  # однократно посчитаем матрицу расстояний между зондами
  ll_tbl <- unique(augment_tbl$idx) %>%
    tidyr::expand_grid(l = ., r = .) %>%
    filter(l != r, (l > 0 | r > 0)) %>%
    # построим уникальный идентификатор ребра
    rowwise() %>%
    # mutate(m = list(sort(c(l, r))))
    mutate(edge_id = stri_c(sort(c(l, r)), collapse = "=")) %>%
    ungroup() %>%
    distinct(edge_id, .keep_all = TRUE) %>%
    # подтягиваем координаты левой точки
    left_join(select(augment_tbl, idx, l_x = x, l_y = y), by = c("l" = "idx")) %>%
    # подтягиваем координаты левой точки
    left_join(select(augment_tbl, idx, r_x = x, r_y = y), by = c("r" = "idx")) %>%
    mutate(s = sqrt((l_x - r_x)^2 + (l_y - r_y)^2)) %>%
    arrange(l, r)

  points_seq <- arrangements::permutations(v = points_in, replace = FALSE, 
                                       layout = "column", nsample = 5000)
  # добавляем точку входа в качестве первой и соотв. точку выхода в качестве последней
  routes_lst <- points_seq %>% 
    rbind(-n1, n1, ., -tail(., 1)) %>%
    as.data.frame() %>% as.list()

  # превращаем все пути обхода в последовательности ребер
  routes_dt <- data.table(route_id = seq_along(routes_lst), route = routes_lst) %>%
    .[, .(from = unlist(route)), by = route_id] %>%
    .[, to := shift(from, type = "lead"), by = route_id] %>%
    # выкидываем все терминальные точки
    na.omit() %>%
    # строим нормализованный идентификатор ребра
    .[, edge_id := stri_c(sort(unlist(.BY)), collapse = "="), by = .(from, to)] %>%
    .[, .(route_id, edge_id)] %>%
    # подтянем информацию о длине ребра из справочника
    .[as.data.table(ll_tbl), s := i.s, on = "edge_id"]

  # считаем длину маршрутов, оставляем кратчайший
  best_routes <- routes_dt[, .(len = sum(s)), by = route_id] %>%
    setorder(len) %>%
    head(10) %T>%
    print()

  # сформируем ТП-10 лучших маршрутов
  best_routes %>%
    select(route_id) %>%
    mutate(idx = routes_lst[route_id]) %>%
    tidyr::unnest(idx) %>%
    left_join(augment_tbl) %>%
    tidyr::nest(data = -route_id) %>%
    left_join(best_routes)
}

tic("Route optimization")
best_tbl <- optimizePath(charges_dt)
toc()

best_route_tbl <- best_tbl$data[[1]]
full_dt <- rbindlist(list(best_route_tbl, keepers_dt), fill = TRUE)
gp <- drawDisk(full_dt) +
  # добавим маршрут обхода
  geom_path(arrow = arrow(type = "closed"), data = best_route_tbl)

gp