TSP.hs

{-# OPTIONS_GHC -Wall                      #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing   #-}
{-# OPTIONS_GHC -fno-warn-type-defaults    #-}
{-# OPTIONS_GHC -fno-warn-unused-do-bind   #-}
{-# OPTIONS_GHC -fno-warn-missing-methods  #-}
{-# OPTIONS_GHC -fno-warn-orphans          #-}

module TSP ( test ) where

import Prelude hiding ( zipWith, sum, foldl, foldr, scanl, length, null, tail, zip, map, reverse, concat, take )
import qualified Prelude as P
import Data.Vector.Unboxed hiding ( replicateM )
import Data.Random.Source.PureMT
import Data.Random
import Control.Monad.State hiding ( modify )

import Debug.Trace

nNodes :: Int
nNodes = 6

xs, ys :: Vector Float
xs = generate nNodes ((*1.0) . fromIntegral)
ys = generate nNodes ((*0.0) . fromIntegral)

distance :: Int -> Int -> Float
distance i j = sqrt $ (xs!i - xs!j)^2 + (ys!i - ys!j)^2

totalDistance :: Vector Int -> Float
totalDistance v = sum $ map (uncurry distance) ms
  where
    ms = zip v (tail v)

reverseBetweenPair :: Int -> Int -> Vector Int -> Vector Int
reverseBetweenPair i j v = concat [beginning, reverse middle, end]
  where
    k = length v

    beginning = slice 0       i           v
    middle    = slice i       (j - i + 1) v
    end       = slice (j + 1) (k - j - 1) v

incDistance :: Vector Int -> Int -> Int -> Float
incDistance v i j = d + c - b - a
  where
    a = distance (v!(i - 1)) (v!i)
    b = distance (v!j)       (v!((j + 1) `mod` nNodes))
    c = distance (v!(i - 1)) (v!j)
    d = distance (v!i)       (v!((j + 1) `mod` nNodes))

expDv :: Int -> Int -> Vector Int -> Float -> Float -> Float -> Float
expDv i1 i2 v kB j t = exp(-j * (incDistance v i1 i2) / (kB * t))

randomUpdates :: Int -> Vector (Int, Int, Float)
randomUpdates m =
  fromList $
  evalState (replicateM m x)
  (pureMT 1)
  where
    x = do r <- sample (uniform (1 :: Int)    (nNodes - 2))
           c <- sample (uniform (r + 1)       (nNodes - 1))
           v <- sample (uniform (0 :: Float)           1.0)
           return (r, c, v)

kB, couplingConstant :: Float
kB               = 1.0
couplingConstant = 1.0

data McState = McState {
  mcRoute :: !(Vector Int)
  }
  deriving Show

initMcState :: McState
initMcState = McState {
  mcRoute = fromList [0,4,2,3,1,5,0]
  }

singleUpdate :: Float -> McState -> (Int, Int, Float) -> McState
singleUpdate t u (i, j, r) =
  trace (show i P.++ show j P.++ ": " P.++
         show (totalDistance v) P.++ ", " P.++
         show (totalDistance (reverseBetweenPair i j v)) P.++ ", " P.++
         show (incDistance v i j) P.++ ", " P.++
         show p P.++ ", " P.++
         show r) $ if incDistance v i j <= 0 || p > r
  then
    McState { mcRoute = reverseBetweenPair i j v }
  else
    u
  where
    v = mcRoute u
    p = expDv i j v kB couplingConstant t

test :: McState
test = foldl (singleUpdate 1.0) initMcState (randomUpdates 100)