shortesttrip.hs

{-
■ プログラムの実行結果(参考 data0 + 解答 data1 data2)
data0.ans (priority queue):
192.343054
35.179325 136.925993
34.681261 135.509801
34.538272 135.536512
34.695124 135.197852

data1.ans (priority queue):
2020.811047
35.640485 140.063040
35.443708 139.638026
34.741327 137.791773
35.599629 138.517272
34.778815 136.908237
35.179325 136.925993
36.146124 137.252173
34.681261 135.509801
33.893756 135.133101
34.695124 135.197852
34.538272 135.536512
36.679768 138.363255
41.083147 141.247793

data2.ans (evolutionary):
7925.326294
35.640485 140.063040
35.443708 139.638026
35.179325 136.925993
34.681261 135.509801
33.724202 135.992530
36.987514 138.577069
43.055460 141.340956
43.076125 141.363598
42.985630 141.563543
35.608932 140.124636
35.152283 140.320932
34.741327 137.791773
38.577071 140.955350
38.299248 140.145926
37.910854 139.121806
37.335416 139.610660
35.840784 140.243831
35.670650 139.771861
33.893756 135.133101
31.429200 131.005865
33.589216 130.392813
34.426380 132.743307
34.674558 132.538416
33.715295 130.446566
34.695124 135.197852
33.931408 134.511356
34.538272 135.536512
35.082676 136.902157
37.916367 139.036326
41.083147 141.247793
37.924871 139.092591
36.639753 139.065865
36.123550 137.866133
35.571384 139.373246
35.364394 139.554156

-}

{-# LANGUAGE MultiParamTypeClasses, TypeFamilies, TupleSections #-}
module Main (main) where

import Data.Function
import Data.Maybe
import Control.Monad.Trans

import Text.Parsec
import Text.Parsec.String
import Text.Printf

import qualified Data.List as L
import qualified Data.IntMap as I
import qualified Data.Map as M
import qualified Data.Heap as H
import qualified Data.Vector.Mutable as MV

import Control.Monad.Primitive
import Control.Monad.Random

import Control.Concurrent
import Control.Parallel
import Control.Parallel.Strategies
import Control.DeepSeq

-- import Debug.Trace


type Point = (Double, Double)
type Name = String


{-
  Parsing
-}
type Oder = [Name]
type Input = ([(Name, Point)], Oder)

float :: Parser Double
float = fmap read $ many (digit <|> char '.')

point :: Parser (Name, Point)
point = do
  name <- many letter
  space
  latitude <- float
  space
  longitude <- float
  newline
  return (name, (latitude, longitude))

number :: Parser Int
number = do
  ds <- many digit
  newline
  return $ read ds

points :: Parser [(Name, Point)]
points = do
  n <- number
  count n point

name :: Parser Name
name = do
    p <- many letter
    newline
    return p

oder :: Parser Oder
oder = do
  n <- number
  l <- count n name
  return l

input :: Parser Input
input = do
  p <- points
  o <- oder
  eof
  return (p, o)


{-
  Geometry
-}

type City = Int
type Path = [City]

data Geometry = Geometry {
  gCityToPoint      :: I.IntMap Point,
  gCityToName       :: I.IntMap Name,
  gNameToCities     :: M.Map Name [City],
  gcitiesToDistance :: I.IntMap (I.IntMap Double)
}

city :: [(Name, Point)] -> (I.IntMap Point, I.IntMap Name, M.Map Name [City])
city n2p = (c2p, c2n, n2cs)
  where
    cities = zipWith (\i (n,p)->(i,n,p)) [0..] n2p
    c2p  = I.fromList $ L.map (\(i, _, p)->(i,p)) cities
    c2n  = I.fromList $ L.map (\(i, n, _)->(i,n)) cities
    n2cs = M.fromListWith (++) $ L.map (\(i, n, _)->(n,[i])) cities

distance :: I.IntMap Point -> I.IntMap (I.IntMap Double)
distance c2p = I.mapWithKey dists c2f
  where
    radian x = x * pi / 180
    factor (la,ln) = (sin $ radian la, cos $ radian la, ln)
    c2f = I.map factor c2p
    dists k1 f1 =
      let dists' k2 f2
            | k1 < k2 = Just $ dist f1 f2
            | otherwise = Nothing
      in I.mapMaybeWithKey dists' c2f
    dist (sla1, cla1, ln1) (sla2, cla2, ln2) =
      acos $ sla1 * sla2 + cla1 * cla2 * cos (radian $ ln1-ln2)

geometry :: [(Name, Point)] -> Geometry
geometry n2p = Geometry c2p c2n n2c cs2d
  where
    (c2p, c2n, n2c) = city n2p
    cs2d = distance c2p

cityToPoint :: Geometry -> City -> Point
cityToPoint (Geometry c2p _ _ _) c = c2p I.! c

cityToName :: Geometry -> City -> Name
cityToName (Geometry _ c2n _ _) c = c2n I.! c

nameToCities :: Geometry -> Name -> [City]
nameToCities (Geometry _ _ n2c _) n = fromJust $ M.lookup n n2c

citiesToDistance' :: I.IntMap (I.IntMap Double) -> City -> City -> Double
citiesToDistance' cs2d x y
  | x == y = 0
  | x < y = cs2d I.! x I.! y
  | x > y = cs2d I.! y I.! x

citiesToDistance :: Geometry -> City -> City -> Double
citiesToDistance (Geometry _ _ _ cs2d) = citiesToDistance' cs2d

pathToLength :: Geometry -> Path -> Double
pathToLength (Geometry _ _ _ cs2d) path = p2l path
  where
    p2l (x:t@(y:_)) = citiesToDistance' cs2d x y + p2l t
    p2l [_] = 0
    p2l []  = 0


{-
  Priority Search
-}

data Item = Item {
    iLength :: Double,
    iPath :: Path
  } deriving Show

instance NFData Item where
  rnf (Item l p) = rnf l `pseq` rnf p


data Policy

instance H.HeapItem Policy Item where
     newtype Prio Policy Item = Pri Double deriving (Show, Eq, Ord)
     type    Val  Policy Item = Path
     split item@(Item l p) = (Pri l, p)
     merge (Pri l, p)      = Item l p


type Queue = H.Heap Policy Item

priorityStep :: Geometry -> Oder -> Item -> [Item]
priorityStep geo oder (Item len path) = parMap rdeepseq step' nexts
  where
    next :: Name
    next = head $ drop (length path) oder
    nexts :: [City]
    nexts = filter (`notElem` path) $ nameToCities geo next
    step' :: City -> Item
    step' p =
      case path of
        [] -> Item 0 [p]
        (p':_) -> Item (citiesToDistance geo p p' + len) (p:path)

priorityLoop :: Geometry -> Oder -> Queue -> Item
priorityLoop geo oder queue = loop queue
  where
    loop queue | length path == length oder = top
               | otherwise = loop $ queue' `H.union` steps
      where
        Just (top@(Item _ path), queue') = H.view queue
        steps = H.fromList $ priorityStep geo oder top

prioritySearch :: Geometry -> Oder -> Item
prioritySearch geo oder =
  let ini = H.singleton $ Item 0 []
  in priorityLoop geo oder ini



{-
  evolution Search
-}
data Machine = Machine {
    thread :: Int,
    capacityPerThread :: Int
  }

type Window = Int


getRandomCity :: (Functor m, Monad m, RandomGen g) =>
                 [City] -> RandT g m City
getRandomCity l = fmap (l L.!!) $ getRandomR (0, length l - 1)

mutation :: (Functor m, Monad m, RandomGen g) =>
            Geometry -> Path -> RandT g m Path
mutation geo path = do
    x <- getRandomCity path
    let others = filter (/=x) $ nameToCities geo $ cityToName geo x
    y <- getRandomCity others
    let swap c | x == c = y
               | y == c = x
               | otherwise = c
    return $ map swap path


getRandomItem :: (PrimMonad m, RandomGen g) =>
                 MV.MVector (PrimState m) Item
                 -> RandT g m Item
getRandomItem mv = do
  i <- getRandomR (0, MV.length mv - 1)
  lift $ MV.read mv i


evolution :: (Functor m, PrimMonad m, RandomGen g) =>
              Window -> Geometry
              -> MV.MVector (PrimState m) Item
              -> RandT g m Item
evolution win geo is = do
      is' <- choose win
      let Item _ path = L.minimumBy (compare `on` iLength) is'
      path' <- mutation geo path
      return $ Item (pathToLength geo path') path'
  where
    choose 0 = return []
    choose n = do
      h <- getRandomItem is
      t <- choose (n-1)
      return $ (h:t)


minBy :: (t -> t -> Ordering) -> t -> t -> t
minBy f x y = case f x y of
                GT -> y
                _ -> x

minByLength :: Item -> Item -> Item
minByLength = minBy (compare `on` iLength)

randomStep :: RandomGen g =>
             Machine
             -> MV.MVector RealWorld Item
             -> RandT g IO Item
             -> g -> IO Item
randomStep (Machine thr cap) result rand gen = step gen 0 thr
  where
    step g i 1 = evalRandT (step' i cap) g
    step gen i thr = do
      let (g, g') = split gen
      mvar <- newEmptyMVar
      forkIO $ do
        item <- evalRandT (step' i cap) g
        let item' = rnf item `pseq` item
        putMVar mvar item'
      item  <- step g' (i+cap) (thr-1)
      item' <- readMVar mvar
      return $ minByLength item item'

    step' i 1 = step'' i
    step' i cap = do
      item  <- step'' i
      item' <- step' (i+1) (cap-1)
      return $ minByLength item item'

    step'' i = do
      item <- rand
      lift $ MV.write result i item
      return item


progenitor :: RandomGen g =>
              Machine -> Geometry -> Oder
              -> MV.MVector RealWorld Item
              -> g -> IO Item
progenitor machine geo oder result gen =
    randomStep machine result (prog oder []) gen
  where
    prog :: (Functor m, PrimMonad m, RandomGen g) =>
            Oder -> Path -> RandT g m Item
    prog [] path = return $ Item (pathToLength geo path) path
    prog (h:t) path = do
      let next = L.filter (`notElem` path) $ nameToCities geo h
      city <- getRandomCity next
      prog t (city:path)


evolutionLoop :: RandomGen g =>
                 Machine -> Window -> Geometry
                 -> MV.MVector RealWorld Item
                 -> MV.MVector RealWorld Item
                 -> g -> Item -> Int -> IO Item
evolutionLoop machine win geo from to gen best i =
  loop from to gen best i
  where
    loop _ _ _ best 0 = return best
    loop from to gen best i = do
      let (g, g') = split gen
      item <- randomStep machine to (evolution win geo from) g
      let best' = minByLength best item
--       trace i best'
      loop to from g' best' (i-1)
--     trace i (Item len path)
--       | i `mod` 10 /= 0 = return ()
--       | otherwise = putStrLn $
--       printf "%06u" i ++ " " ++
--       (printf "%015.6f" $ 6378.137 * len) ++ " " ++
--       show path
-- 
evolutionSearch :: Machine -> Window -> Geometry -> Oder -> Int -> IO Item
evolutionSearch machine@(Machine thr cap) win geo oder loop = do
  let gen = mkStdGen $ 10^10
      (g, g') = split gen

  from <- MV.new $ thr * cap
  to <- MV.new $ thr * cap

  best <- progenitor machine geo oder from g
  evolutionLoop machine win geo from to g' best loop


{-
  main
-}

readInput :: String -> IO Input
readInput file = do
  i <- parseFromFile input file
  case i of
    Left err -> error $ show err
    Right i' -> return i'

readData :: String -> IO (String, Geometry, Oder)
readData title = do
  (table, oder) <- readInput $ title ++ ".in"
  return (title, geometry table, oder)


pretty :: Geometry -> City -> String
pretty geo city = pr $ cityToPoint geo city
  where
    pr (la, ln) = printf "%0.6f %0.6f" la ln

search :: String -> String -> Geometry -> IO Item -> IO ()
search title method geo searcher = do
  putStrLn $ title ++ ".ans (" ++ method ++ "):"
  Item len path <- searcher
  putStrLn $ printf "%0.6f" $ 6378.137 * len
  mapM_ (putStrLn.pretty geo) $ reverse $ path
  putStrLn ""

searchPriority (title, geo, oder) = do
    search title "priority queue" geo searcher
  where
    searcher = return $ prioritySearch geo oder


searchEvolution (title, geo, oder) thr cap win loop = do
    search title "evolutionary" geo searcher
  where
    searcher = evolutionSearch (Machine thr cap) win geo oder loop


main = do
  data0 <- readData "data0"
  data1 <- readData "data1"
  data2 <- readData "data2"

  searchPriority data0
  searchPriority data1
--  searchPriority data2

--  searchEvolution data0 26 100  5  5
--  searchEvolution data1 26 100  5 20
  searchEvolution data2 26 200 10 30