lat-long.clj

(require '[clojure.math.combinatorics :as combo]
         '[clojure.core.async :as async])

(def zeus-route
  ["Kadena"
   "Pohang"
   "Pusan"
   "Kimhae"
   "Sachon"
   "Taegu"
   "Kumi"
   "R601"
   "Yechon"
   "Yongju"
   "Choongwon"
   "R605"
   "Hoengsong"
   "R419"
   "Kangnung"
   "Sokcho"
   "Chunchon"
   "R217"
   "R222"
   "R218"
   "Mandumi"
   "R107"
   "R113"
   "Kimpo"
   "R103"
   "Seoul"
   "Singal"
   "Suwon"
   "Osan"
   "P'Yong'Taeg"
   "Songwhan"
   "Chongju"
   "Chongwon"
   "R505"
   "R110"
   "Seosan"
   "Kwangju"
   "Kunsan"])

(def best-route
  ["Kadena"
   "Pohang"
   "Pusan"
   "Kimhae"
   "Sachon"
   "Taegu"
   "Kumi"
   "R601"
   "Yechon"
   "Yongju"
   "R605"
   "Choongwon"
   "Hoengsong"
   "R419"
   "Kangnung"
   "Sokcho"
   "Chunchon"
   "R217"
   "R218"
   "Mandumi"
   "R222"
   "R113"
   "R107"
   "Kimpo"
   "R103"
   "Seoul"
   "Singal"
   "Suwon"
   "Osan"
   "P'Yong'Taeg"
   "Songwhan"
   "Chongju"
   "Chongwon"
   "R505"
   "R110"
   "Seosan"
   "Kwangju"
   "Kunsan"])

(def airbases
  {"Chongju"     [[36 48.100] [128 36.314]]
   "Choongwon"   [[36 59.565] [129 08.769]]
   "Hoengsong"   [[37 27.477] [129 13.341]]
   "Kangnung"    [[37 47.315] [130 28.979]]
   "Kimhae"      [[35 12.856] [130 11.734]]
   "Kimpo"       [[37 36.868] [127 42.801]]
   "Kunsan"      [[35 57.554] [127 24.492]]
   "Kwangju"     [[35 12.906] [127 40.157]]
   "Osan"        [[37 04.141] [128 00.881]]
   "Pohang"      [[36 02.779] [130 48.826]]
   "Pusan"       [[35 14.353] [130 29.144]]
   "P'Yong'Taeg" [[36 57.931] [128 03.884]]
   "Sachon"      [[35 07.972] [129 10.731]]
   "Seosan"      [[36 41.807] [127 19.884]]
   "Seoul"       [[37 26.719] [128 07.510]]
   "Sokcho"      [[38 06.185] [130 07.622]]
   "Suwon"       [[37 16.341] [127 57.800]]
   "Taegu"       [[35 57.999] [129 55.962]]
   "Yechon"      [[36 41.576] [129 40.163]]
   "Chongwon"    [[36 37.265] [128 35.318]]
   "Chunchon"    [[37 53.781] [128 56.751]]
   "Kumi"        [[36 12.633] [129 28.087]]
   "Mandumi"     [[38 04.259] [127 57.586]]
   "R103"        [[37 28.610] [127 39.140]]
   "R107"        [[37 46.237] [127 28.889]]
   "R110"        [[36 49.119] [127 39.642]]
   "R113"        [[37 43.159] [127 49.588]]
   "R217"        [[37 58.092] [128 18.822]]
   "R218"        [[38 03.480] [128 04.842]]
   "R222"        [[37 52.703] [128 06.146]]
   "R419"        [[37 41.926] [129 12.823]]
   "R505"        [[36 37.265] [128 19.205]]
   "R601"        [[36 41.731] [129 25.817]]
   "R605"        [[37 10.673] [129 31.878]]
   "Singal"      [[37 19.295] [128 07.928]]
   "Songwhan"    [[36 54.508] [128 09.622]]
   "Yongju"      [[36 51.275] [129 55.131]]
   "Kadena"      [[34 39.350] [132 55.367]]})

;; 25
(def dprk-airbases
  {"Haeju"          [[38 02.223] [126 36.807]]
   "Hwangju"        [[38 42.959] [126 32.227]]
   "Hwangsuwon"     [[40 43.756] [129 50.507]]
   "Hyon-Ni"        [[38 41.300] [128 37.565]]
   "Iwon"           [[40 24.718] [130 26.856]]
   "Kaech'on"       [[39 47.400] [126 42.427]]
   "Koksan"         [[38 48.157] [127 32.883]]
   "Kuum-ni"        [[38 54.190] [129 17.199]]
   "Kwail"          [[38 26.721] [125 29.149]]
   "Manp'o"         [[41 11.801] [127 27.140]]
   "Mirim"          [[39 03.544] [126 39.615]]
   "Onch'on"        [[38 57.412] [125 45.757]]
   "Ongjin"         [[37 59.158] [126 01.761]]
   "Orang"          [[41 27.626] [131 48.613]]
   "Panghyon"       [[39 58.840] [125 51.027]]
   "Pukch'ang-up"   [[39 37.849] [126 44.436]]
   "Samjiyon"       [[41 55.243] [130 13.612]]
   "Sondok"         [[39 47.577] [128 47.130]]
   "Sunan"          [[39 15.720] [126 21.840]]
   "Sunch'on"       [[39 26.670] [126 42.318]]
   "T'aech'on"      [[39 57.415] [126 10.724]]
   "Taetan"         [[38 08.607] [125 51.628]]
   "Toksan"         [[40 03.690] [128 59.915]]
   "Uiju"           [[40 15.017] [125 01.967]]
   "Wonsan"         [[39 13.084] [128 45.577]]})

;; 19
(def dprk-airstrips
  "Ayang-Ni"           [[38 18.184] [126 43.238]]
  "Hoeyang SE"         [[38 47.282] [128 54.712]]
  "Ich'on"             [[38 36.505] [127 51.762]]
  "Kaech'on SW"        [[39 51.249] [126 37.207]]
  "Kilchu"             [[40 59.144] [131 18.570]]
  "Kojo"               [[38 52.132] [129 11.725]]
  "Koksan"             [[38 47.665] [127 35.473]]
  "Kwaksan"            [[39 43.861] [125 44.733]]
  "Kyongsong Chuul"    [[41 35.786] [131 48.478]]
  "Nuch'on-Ni"         [[38 16.568] [127 07.865]]
  "Okpyong-Ni"         [[39 19.613] [128 30.958]]
  "Panghyon"           [[39 56.254] [125 46.633]]
  "Pongsan"            [[38 32.194] [126 32.937]]
  "Sangwon"            [[38 53.209] [126 50.598]]
  "Sonch'on"           [[39 47.633] [125 19.635]]
  "Sugam-Ni"           [[41 46.179] [131 56.882]]
  "Taebukpo-Ri"        [[38 21.417] [127 46.622]]
  "Tangch'on"          [[40 25.352] [130 36.617]]
  "Yonghung"           [[39 33.084] [128 33.424]])

(def russia
  "Nachodka" [[42 59.916] [134 00.718]])

;; Where's the missing one?

(def aviator2-stops
  (-> (merge dprk-airbases dprk-airstrips russia)
      (disj "Nachodka")
      (disj "Wonsan")))

(def aviator1-stops (-> airbases keys set (disj "Kadena") (disj "Kunsan")))

(defn deg->rad
  [d]
  (/ (* d Math/PI) 180.0))

(defn rad->deg
  [r]
  (/ (* r 180.0) Math/PI))

(defn dm->deg [[d m]]
  (+ d (/ m 60.0)))

(defn dms->deg [[d m s]]
  (+ d (/ m 60.0) (/ s 3600.0)))

(defn dms->rad [dms]
  (-> dms dms->deg deg->rad))

(defn dm->rad [dm]
  (-> dm dm->deg deg->rad))

(defn distance [[lata longa] [latb longb]]
  (let [theta (- (dm->rad longa) (dm->rad longb))]
    (-> (+ (* (Math/sin (dm->rad lata))
              (Math/sin (dm->rad latb)))
           (* (Math/cos (dm->rad lata))
              (Math/cos (dm->rad latb))
              (Math/cos theta)))
        Math/acos
        rad->deg
        (* 60 1.1515))))

(defn ab-distance
  [node1 node2]
  (distance (airbases node1) (airbases node2)))

(def ab-distance*
  (memoize ab-distance))

(defn routes
  "Return all the permutations of nodes that begin at `start` and
  terminate at `end`"
  [nodes start end]
  (->> nodes
       combo/permutations
       (map #(conj (into [start] %) end))))

(defn route-distance
  [route]
  (->> (map (fn [a b]
              (ab-distance* a b))
            route
            (drop 1 route))
       (reduce +)))

(defn all-routes
  [nodes start end]
  (for [route (routes nodes start end)]
    {:distance (route-distance route)
     :route route}))

(defn first-route
  [nodes start end]
  (loop [remaining (set nodes)
         route [start]
         last-node start]
    (if (empty? remaining)
      (conj route end)
      (let [closest (->> remaining
                         (map (fn [node]
                                {:distance (ab-distance* node last-node)
                                 :node node}))
                         (apply min-key :distance)
                         :node)]
        (recur (disj remaining closest)
               (conj route closest)
               closest)))))

(defn optimal-route
  "Shortest route among the passed nodes"
  [nodes start end]
  (apply min-key :distance (all-routes nodes start end)))

(defn improve
  [generator scorer score model index]
  (->> (repeatedly #(generator model index))
       (map (fn [m] [(scorer m) m]))
       (drop-while (fn [[s m]] (<= s score)))
       first))

(defn- model-stream*
  "Internal implementation of model-stream"
  [model score scorer generator index]
  (let [[next-score next-model] (improve generator scorer score model index)]
    (lazy-seq
     (cons model
           (when next-model
             (model-stream* next-model next-score scorer generator (inc index)))))))

(defn model-stream
  "Returns an infinite sequence of models starting with `initial`. Each
  successive model will have a higher fitness (as determined by
  calling the `fitness-fn` function on it). Creates candidate models
  with `generator`, a function of the previous model and the iteration
  number.

  `(generator model index)` will be called repeatedly until it either
  returns a more-fit model or nil."
  [initial fitness-fn generator]
  (model-stream* initial (fitness-fn initial) fitness-fn generator 0))

(defn rand-between
  "Return a random integer between a and b, inclusive"
  [a b]
  (-> (- b a) inc rand-int (+ a)))

(defn pairwise-exchange
  "Returns a new route with a random subset reversed."
  [route _]
  (let [n (count route)
        a (rand-between 1 (- n 2))
        b (rand-between (inc a) (- n 1))]
    (reduce into
            []
            [(subvec route 0 a)
             (reverse (subvec route a b))
             (subvec route b n)])))

(defn route-fitness
  [route]
  (- (route-distance route)))

(defn scramble-fn
  [rounds p-skip]
  (fn [route _]
    (->> (iterate (fn [r]
                    (if (< (rand) p-skip)
                      route
                      (pairwise-exchange r nil)))
                  route)
         (drop rounds)
         first)))

(comment
  (all-routes ["Seoul" "Taegu"] "Kimpo" "Kunsan")

  (let [nodes (-> (keys airbases)
                  set
                  (disj "Kadena")
                  (disj "Kunsan"))]
    (->> (all-routes nodes "Kadena" "Kunsan")
         (take 100000)
         dorun
         time))

  (let [nodes (-> (keys airbases)
                  set
                  (disj "Kadena")
                  (disj "Kunsan"))]
    (first-route nodes "Kadena" "Kunsan"))

  (let [nodes (-> (keys airbases)
                  set
                  (disj "Kadena")
                  (disj "Kunsan"))]
    (optimal-route nodes "Kadena" "Kunsan"))

  (map (fn [a b]
         [[a b] (ab-distance a b)])
       zeus-route
       (drop 1 zeus-route))

  (route-distance zeus-route)           ; => 1210.0821092501515
  (route-distance best-route)           ; => 1200.3717719438573

  (route-distance (first-route aviator1-stops "Kadena" "Kunsan")) ; => 1371.3848846665865

  ;; Aviator 1 computation
  (let [n 8
        l (Object.)]
    (doseq [[rounds p-skip] (for [rounds [1 2 4 8]
                                  p-skip [0.1 0.5]]
                              [rounds p-skip])]
      (future (doseq [route (model-stream best-route #_zeus-route #_(first-route stops "Kadena" "Kunsan") route-fitness (scramble-fn rounds p-skip))]
                (locking l
                  (println (route-distance route) rounds p-skip route))))))

  ;; Aviator 2 computation
  (let [n 8
        l (Object.)]
    (doseq [[rounds p-skip] (for [rounds [1 2 4 8]
                                  p-skip [0.1 0.5]]
                              [rounds p-skip])]
      (future (doseq [route (model-stream #_best-route #_zeus-route (first-route aviator2-stops "Nachodka" "Wonsan") route-fitness (scramble-fn rounds p-skip))]
                (locking l
                  (println (route-distance route) rounds p-skip route))))))

  )