l2_norm <- function(x)
  sqrt(sum(x*x))

##' Angle between the segments 'ab' and 'ac'.
##'
##' @param a Central point
##' @param b First neighbor
##' @param c Second neighbor
##'
##' @return The (smaller) unsigned angle between 'ab' and 'ac'.
angle_between_points <- function(a, b, c) {
  v1 <- (a - b) / l2_norm(a - b)
  v2 <- (a - c) / l2_norm(a - c)
  angle <- acos(sum(v1 * v2))
  min(abs(angle + c(0, 2*pi, -2*pi)))
}

##' Angle feature
##'
##' @param x A TSP instance
##'
##' @return The average angle betwen the two nearest citities of every
##' city.
feature_angle <- function(x) {
  d <- as.matrix(dist(x))
  number_of_cities <- nrow(x)
  angle <- numeric(number_of_cities)
  for (city in 1:number_of_cities) {
    neighbor <- order(d[city,])[2:3] ## 2 next neighbors.
    angle[city] <- angle_between_points(x[city,],
                                        x[neighbor[1], ],
                                        x[neighbor[2], ])
  }
  list(angle_min=min(angle, na.rm=TRUE),
       angle_mean=mean(angle, na.rm=TRUE),
       angle_median=median(angle, na.rm=TRUE),
       angle_max=max(angle, na.rm=TRUE),
       angle_sd=sd(angle, na.rm=TRUE)
       )
}
