-- | TSP attempts to find a closed loop of minimum length that visits each input
-- point exactly once.
--
-- The general workflow is:
--
--  1. Create a basic data structure with 'inputPath' or 'distances'+'tspNearestNeighbourPath'
--  2. Improve it using 'tsp2optPath', 'tsp3optPath'
--  3. Retrieve the optimized geometry with 'pathToPoints'
module Numerics.Optimization.TSP (
    -- * Setup
      distances
    , Distances
    , pathToPoints
    , TspPoint

    -- * Construction
    , inputPath
    , tspNearestNeighbourPath

    -- * Improvement

    -- ** 2-opt
    , tsp2optPath
    , tsp2optInplace
    , unsafePlan2optSwap
    , unsafeCommit2optSwap

    -- ** 3-opt
    , tsp3optPath
    , tsp3optInplace
    , Case3opt
    , unsafePlan3optSwap
    , unsafeCommit3optSwap

) where



import           Control.Monad.Primitive
import           Data.Coerce
import           Data.Vector             (Vector, (!))
import qualified Data.Vector.Extended    as V
import qualified Data.Vector.Mutable     as VM
import           Geometry.Core
import           Text.Printf



-- $setup
-- >>> import Draw
-- >>> import Geometry.Algorithms.Sampling
-- >>> import qualified Data.Vector as V
-- >>> import Graphics.Rendering.Cairo as C
-- >>> import qualified System.Random.MWC as MWC



newtype Temperature = T Double deriving (Temperature -> Temperature -> Bool
(Temperature -> Temperature -> Bool)
-> (Temperature -> Temperature -> Bool) -> Eq Temperature
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Temperature -> Temperature -> Bool
== :: Temperature -> Temperature -> Bool
$c/= :: Temperature -> Temperature -> Bool
/= :: Temperature -> Temperature -> Bool
Eq, Eq Temperature
Eq Temperature
-> (Temperature -> Temperature -> Ordering)
-> (Temperature -> Temperature -> Bool)
-> (Temperature -> Temperature -> Bool)
-> (Temperature -> Temperature -> Bool)
-> (Temperature -> Temperature -> Bool)
-> (Temperature -> Temperature -> Temperature)
-> (Temperature -> Temperature -> Temperature)
-> Ord Temperature
Temperature -> Temperature -> Bool
Temperature -> Temperature -> Ordering
Temperature -> Temperature -> Temperature
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Temperature -> Temperature -> Ordering
compare :: Temperature -> Temperature -> Ordering
$c< :: Temperature -> Temperature -> Bool
< :: Temperature -> Temperature -> Bool
$c<= :: Temperature -> Temperature -> Bool
<= :: Temperature -> Temperature -> Bool
$c> :: Temperature -> Temperature -> Bool
> :: Temperature -> Temperature -> Bool
$c>= :: Temperature -> Temperature -> Bool
>= :: Temperature -> Temperature -> Bool
$cmax :: Temperature -> Temperature -> Temperature
max :: Temperature -> Temperature -> Temperature
$cmin :: Temperature -> Temperature -> Temperature
min :: Temperature -> Temperature -> Temperature
Ord, Int -> Temperature -> ShowS
[Temperature] -> ShowS
Temperature -> String
(Int -> Temperature -> ShowS)
-> (Temperature -> String)
-> ([Temperature] -> ShowS)
-> Show Temperature
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Temperature -> ShowS
showsPrec :: Int -> Temperature -> ShowS
$cshow :: Temperature -> String
show :: Temperature -> String
$cshowList :: [Temperature] -> ShowS
showList :: [Temperature] -> ShowS
Show)
newtype Energy = E Double deriving (Energy -> Energy -> Bool
(Energy -> Energy -> Bool)
-> (Energy -> Energy -> Bool) -> Eq Energy
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Energy -> Energy -> Bool
== :: Energy -> Energy -> Bool
$c/= :: Energy -> Energy -> Bool
/= :: Energy -> Energy -> Bool
Eq, Eq Energy
Eq Energy
-> (Energy -> Energy -> Ordering)
-> (Energy -> Energy -> Bool)
-> (Energy -> Energy -> Bool)
-> (Energy -> Energy -> Bool)
-> (Energy -> Energy -> Bool)
-> (Energy -> Energy -> Energy)
-> (Energy -> Energy -> Energy)
-> Ord Energy
Energy -> Energy -> Bool
Energy -> Energy -> Ordering
Energy -> Energy -> Energy
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Energy -> Energy -> Ordering
compare :: Energy -> Energy -> Ordering
$c< :: Energy -> Energy -> Bool
< :: Energy -> Energy -> Bool
$c<= :: Energy -> Energy -> Bool
<= :: Energy -> Energy -> Bool
$c> :: Energy -> Energy -> Bool
> :: Energy -> Energy -> Bool
$c>= :: Energy -> Energy -> Bool
>= :: Energy -> Energy -> Bool
$cmax :: Energy -> Energy -> Energy
max :: Energy -> Energy -> Energy
$cmin :: Energy -> Energy -> Energy
min :: Energy -> Energy -> Energy
Ord, Int -> Energy -> ShowS
[Energy] -> ShowS
Energy -> String
(Int -> Energy -> ShowS)
-> (Energy -> String) -> ([Energy] -> ShowS) -> Show Energy
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Energy -> ShowS
showsPrec :: Int -> Energy -> ShowS
$cshow :: Energy -> String
show :: Energy -> String
$cshowList :: [Energy] -> ShowS
showList :: [Energy] -> ShowS
Show)

-- | Infix shorthand for 'mod' with a fixity of one below addition.
--
-- >>> 7+3 % 4
-- 2
(%) :: Int -> Int -> Int
% :: Int -> Int -> Int
(%) = Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod
infix 5 %

-- | Cache the distance between points i and j.
--
-- The 'Vec2's can be forgotten now, because we’re only interested in topological
-- properties, not the coordinates themselves. Once the calculation is done, the
-- result permutation can be obtained with 'pathToPoints'.
distances :: Vector Vec2 -> Distances
distances :: Vector Vec2 -> Distances
distances Vector Vec2
points = Vector (Vector Double) -> Distances
Distances (Vector (Vector Double) -> Distances)
-> Vector (Vector Double) -> Distances
forall a b. (a -> b) -> a -> b
$
    Int -> (Int -> Vector Double) -> Vector (Vector Double)
forall a. Int -> (Int -> a) -> Vector a
V.generate (Vector Vec2 -> Int
forall a. Vector a -> Int
V.length Vector Vec2
points) ((Int -> Vector Double) -> Vector (Vector Double))
-> (Int -> Vector Double) -> Vector (Vector Double)
forall a b. (a -> b) -> a -> b
$ \Int
i ->
        Int -> (Int -> Double) -> Vector Double
forall a. Int -> (Int -> a) -> Vector a
V.generate Int
i ((Int -> Double) -> Vector Double)
-> (Int -> Double) -> Vector Double
forall a b. (a -> b) -> a -> b
$ \Int
j ->
            Vec2 -> Double
norm (Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
i Vec2 -> Vec2 -> Vec2
forall v. VectorSpace v => v -> v -> v
-. Vector Vec2
pointsVector Vec2 -> Int -> Vec2
forall a. Vector a -> Int -> a
!Int
j)

-- | Lookup the distance between two points.
--
-- @
-- 'distance' ('distances' points) i j = 'norm' (points'!'i -. points'!'j)
-- @
distance :: Distances -> TspPoint -> TspPoint -> Double
distance :: Distances -> TspPoint -> TspPoint -> Double
distance (Distances Vector (Vector Double)
matrix) (TspPoint Int
i) (TspPoint Int
j) = case Int -> Int -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Int
i Int
j of
    Ordering
EQ -> Double
0
    Ordering
LT -> Vector (Vector Double)
matrixVector (Vector Double) -> Int -> Vector Double
forall a. Vector a -> Int -> a
!Int
jVector Double -> Int -> Double
forall a. Vector a -> Int -> a
!Int
i
    Ordering
GT -> Vector (Vector Double)
matrixVector (Vector Double) -> Int -> Vector Double
forall a. Vector a -> Int -> a
!Int
iVector Double -> Int -> Double
forall a. Vector a -> Int -> a
!Int
j

-- | Cache for the (symmetrical) distance between points. See 'distances'.
newtype Distances = Distances (Vector (Vector Double))

-- | Number of points
size :: Distances -> Int
size :: Distances -> Int
size (Distances Vector (Vector Double)
m) = Vector (Vector Double) -> Int
forall a. Vector a -> Int
V.length Vector (Vector Double)
m

instance Show Distances where
    show :: Distances -> String
show m :: Distances
m@(Distances Vector (Vector Double)
mxy) = do
        TspPoint
i <- [Int] -> [TspPoint]
forall a b. Coercible a b => a -> b
coerce [Int
0 .. Vector (Vector Double) -> Int
forall a. Vector a -> Int
V.length Vector (Vector Double)
mxyInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
        TspPoint
j <- [Int] -> [TspPoint]
forall a b. Coercible a b => a -> b
coerce [Int
0 .. Vector (Vector Double) -> Int
forall a. Vector a -> Int
V.length Vector (Vector Double)
mxyInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
        if TspPoint
j TspPoint -> TspPoint -> Bool
forall a. Eq a => a -> a -> Bool
== Int -> TspPoint
TspPoint Int
0
            then String -> Double -> String
forall r. PrintfType r => String -> r
printf String
"[%6.2f, " (Distances -> TspPoint -> TspPoint -> Double
distance Distances
m TspPoint
i TspPoint
j)
            else if TspPoint
j TspPoint -> TspPoint -> Bool
forall a. Ord a => a -> a -> Bool
>= Int -> TspPoint
TspPoint (Vector (Vector Double) -> Int
forall a. Vector a -> Int
V.length Vector (Vector Double)
mxy Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
                then String -> Double -> String
forall r. PrintfType r => String -> r
printf String
"%6.2f]\n" (Distances -> TspPoint -> TspPoint -> Double
distance Distances
m TspPoint
i TspPoint
j)
                else String -> Double -> String
forall r. PrintfType r => String -> r
printf String
"%6.2f, " (Distances -> TspPoint -> TspPoint -> Double
distance Distances
m TspPoint
i TspPoint
j)

-- | Permute the input vector with the result of a TSP optimization.
pathToPoints :: Vector vec2 -> Vector TspPoint -> Vector vec2
pathToPoints :: forall vec2. Vector vec2 -> Vector TspPoint -> Vector vec2
pathToPoints Vector vec2
coordinates Vector TspPoint
indices = Vector vec2 -> Vector Int -> Vector vec2
forall a. Vector a -> Vector Int -> Vector a
V.backpermute Vector vec2
coordinates (Vector TspPoint -> Vector Int
forall a b. Coercible a b => a -> b
coerce Vector TspPoint
indices)

-- | Interpret the input as a path as-is. The simplest possible TSP solver. :-)
--
-- <<docs/haddock/Numerics/Optimization/TSP/input_path.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Numerics/Optimization/TSP/input_path.svg" 300 200 $ \_ -> do
--    points <- liftIO $ do
--        gen <- MWC.initialize (V.fromList [])
--        uniform <- uniformlyDistributedPoints gen (shrinkBoundingBox 10 [zero, Vec2 300 200]) 32
--        pure uniform
--    setLineJoin LineJoinRound
--    let path = inputPath points
--        tour = Polygon (V.toList (pathToPoints points path))
--    cairoScope $ do
--        setColor (mma 0)
--        for_ points $ \p -> sketch (Circle p 3) >> fill
--    cairoScope $ do
--        setColor (mma 1)
--        sketch tour
--        stroke
-- :}
-- Generated file: size 11KB, crc32: 0xde10a233
inputPath :: Vector a -> Vector TspPoint
inputPath :: forall a. Vector a -> Vector TspPoint
inputPath = (Vector Int -> Vector TspPoint
forall a b. Coercible a b => a -> b
coerce :: Vector Int -> Vector TspPoint) (Vector Int -> Vector TspPoint)
-> (Vector a -> Vector Int) -> Vector a -> Vector TspPoint
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int -> Vector Int
forall a. Num a => a -> Int -> Vector a
V.enumFromN Int
0 (Int -> Vector Int) -> (Vector a -> Int) -> Vector a -> Vector Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector a -> Int
forall a. Vector a -> Int
V.length

-- | Starting at the \(i\)-th point, create a path by taking the nearest
-- neighbour. This algorithm is very simple and fast, but yields surprisingly short
-- paths for everyday inputs.
--
-- It is an excellent starting point for more accurate but slower algorithms, such
-- as 'tsp2optPath'.
--
-- <<docs/haddock/Numerics/Optimization/TSP/nearest_neighbour_path.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Numerics/Optimization/TSP/nearest_neighbour_path.svg" 300 200 $ \_ -> do
--    points <- liftIO $ do
--        gen <- MWC.initialize (V.fromList [])
--        uniform <- uniformlyDistributedPoints gen (shrinkBoundingBox 10 [zero, Vec2 300 200]) 32
--        pure uniform
--    setLineJoin LineJoinRound
--    let dm = distances points
--        path = tspNearestNeighbourPath dm
--        tour = Polygon (V.toList (pathToPoints points path))
--    cairoScope $ do
--        setColor (mma 0)
--        for_ points $ \p -> sketch (Circle p 3) >> fill
--    cairoScope $ do
--        setColor (mma 2)
--        sketch tour
--        stroke
-- :}
-- Generated file: size 11KB, crc32: 0x1474da44
tspNearestNeighbourPath :: Distances -> Vector TspPoint
tspNearestNeighbourPath :: Distances -> Vector TspPoint
tspNearestNeighbourPath Distances
dm = (forall s. ST s (MVector s TspPoint)) -> Vector TspPoint
forall a. (forall s. ST s (MVector s a)) -> Vector a
V.create ((forall s. ST s (MVector s TspPoint)) -> Vector TspPoint)
-> (forall s. ST s (MVector s TspPoint)) -> Vector TspPoint
forall a b. (a -> b) -> a -> b
$ do
    let n :: Int
n = Distances -> Int
size Distances
dm
    MVector s TspPoint
points <- Int
-> (Int -> TspPoint) -> ST s (MVector (PrimState (ST s)) TspPoint)
forall (m :: * -> *) a.
PrimMonad m =>
Int -> (Int -> a) -> m (MVector (PrimState m) a)
VM.generate Int
n Int -> TspPoint
TspPoint
    let loop :: MVector (PrimState m) TspPoint -> m (MVector s TspPoint)
loop MVector (PrimState m) TspPoint
vec = do
            Maybe (MVector (PrimState m) TspPoint)
result <- Distances
-> MVector (PrimState m) TspPoint
-> m (Maybe (MVector (PrimState m) TspPoint))
forall (m :: * -> *).
PrimMonad m =>
Distances
-> MVector (PrimState m) TspPoint
-> m (Maybe (MVector (PrimState m) TspPoint))
nnStep Distances
dm MVector (PrimState m) TspPoint
vec
            case Maybe (MVector (PrimState m) TspPoint)
result of
                Maybe (MVector (PrimState m) TspPoint)
Nothing -> MVector s TspPoint -> m (MVector s TspPoint)
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure MVector s TspPoint
points
                Just MVector (PrimState m) TspPoint
rest -> MVector (PrimState m) TspPoint -> m (MVector s TspPoint)
loop MVector (PrimState m) TspPoint
rest
    MVector (PrimState (ST s)) TspPoint -> ST s (MVector s TspPoint)
forall {m :: * -> *}.
PrimMonad m =>
MVector (PrimState m) TspPoint -> m (MVector s TspPoint)
loop MVector s TspPoint
MVector (PrimState (ST s)) TspPoint
points

-- | Efficient, but less obvious algorithm than the straightforward search. The
-- advantage here is that instead of maintaining a »seen« data structure that we
-- have to traverse fully on each iteration, we only have to search not-yet-seen
-- points by design in this approach.
--
-- 1. Split the input vector in (head, tail). Use the head as the current point.
-- 2. Search the tail for a closer point.
--    If one is found, swap it to the beginning of tail. The first element of
--    the tail will become the closest point to the current point.
-- 3. Return the tail; iterating this algorithm on it will sort the list.
nnStep
    :: PrimMonad m
    => Distances
    -> VM.MVector (PrimState m) TspPoint
    -> m (Maybe (VM.MVector (PrimState m) TspPoint))
nnStep :: forall (m :: * -> *).
PrimMonad m =>
Distances
-> MVector (PrimState m) TspPoint
-> m (Maybe (MVector (PrimState m) TspPoint))
nnStep Distances
_ MVector (PrimState m) TspPoint
vec | MVector (PrimState m) TspPoint -> Int
forall s a. MVector s a -> Int
VM.length MVector (PrimState m) TspPoint
vec Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1 = Maybe (MVector (PrimState m) TspPoint)
-> m (Maybe (MVector (PrimState m) TspPoint))
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure Maybe (MVector (PrimState m) TspPoint)
forall a. Maybe a
Nothing
nnStep Distances
dm MVector (PrimState m) TspPoint
vec = do
    let (MVector (PrimState m) TspPoint
hd, MVector (PrimState m) TspPoint
tl) = Int
-> MVector (PrimState m) TspPoint
-> (MVector (PrimState m) TspPoint, MVector (PrimState m) TspPoint)
forall s a. Int -> MVector s a -> (MVector s a, MVector s a)
VM.splitAt Int
1 MVector (PrimState m) TspPoint
vec
    TspPoint
current <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
hd Int
0
    let loop :: Double -> Int -> f ()
loop Double
_ Int
i | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= MVector (PrimState m) TspPoint -> Int
forall s a. MVector s a -> Int
VM.length MVector (PrimState m) TspPoint
tl = () -> f ()
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ()
        loop Double
minDist Int
i = do
            TspPoint
candidate <- MVector (PrimState f) TspPoint -> Int -> f TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
MVector (PrimState f) TspPoint
tl Int
i
            let distanceFromCurrent :: Double
distanceFromCurrent = Distances -> TspPoint -> TspPoint -> Double
distance Distances
dm TspPoint
current TspPoint
candidate
            if Double
distanceFromCurrent Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
minDist
                then MVector (PrimState f) TspPoint -> Int -> Int -> f ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> Int -> m ()
VM.swap MVector (PrimState m) TspPoint
MVector (PrimState f) TspPoint
tl Int
0 Int
i f () -> f () -> f ()
forall a b. f a -> f b -> f b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Double -> Int -> f ()
loop Double
distanceFromCurrent (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                else Double -> Int -> f ()
loop Double
minDist (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
    Double -> Int -> m ()
forall {f :: * -> *}.
(PrimState f ~ PrimState m, PrimMonad f) =>
Double -> Int -> f ()
loop (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0) Int
0
    Maybe (MVector (PrimState m) TspPoint)
-> m (Maybe (MVector (PrimState m) TspPoint))
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MVector (PrimState m) TspPoint
-> Maybe (MVector (PrimState m) TspPoint)
forall a. a -> Maybe a
Just MVector (PrimState m) TspPoint
tl)

-- | Abstract point in the TSP algorithms. Corresponds to the \(i\)-th input point.
newtype TspPoint = TspPoint Int
    deriving (TspPoint -> TspPoint -> Bool
(TspPoint -> TspPoint -> Bool)
-> (TspPoint -> TspPoint -> Bool) -> Eq TspPoint
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: TspPoint -> TspPoint -> Bool
== :: TspPoint -> TspPoint -> Bool
$c/= :: TspPoint -> TspPoint -> Bool
/= :: TspPoint -> TspPoint -> Bool
Eq, Eq TspPoint
Eq TspPoint
-> (TspPoint -> TspPoint -> Ordering)
-> (TspPoint -> TspPoint -> Bool)
-> (TspPoint -> TspPoint -> Bool)
-> (TspPoint -> TspPoint -> Bool)
-> (TspPoint -> TspPoint -> Bool)
-> (TspPoint -> TspPoint -> TspPoint)
-> (TspPoint -> TspPoint -> TspPoint)
-> Ord TspPoint
TspPoint -> TspPoint -> Bool
TspPoint -> TspPoint -> Ordering
TspPoint -> TspPoint -> TspPoint
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: TspPoint -> TspPoint -> Ordering
compare :: TspPoint -> TspPoint -> Ordering
$c< :: TspPoint -> TspPoint -> Bool
< :: TspPoint -> TspPoint -> Bool
$c<= :: TspPoint -> TspPoint -> Bool
<= :: TspPoint -> TspPoint -> Bool
$c> :: TspPoint -> TspPoint -> Bool
> :: TspPoint -> TspPoint -> Bool
$c>= :: TspPoint -> TspPoint -> Bool
>= :: TspPoint -> TspPoint -> Bool
$cmax :: TspPoint -> TspPoint -> TspPoint
max :: TspPoint -> TspPoint -> TspPoint
$cmin :: TspPoint -> TspPoint -> TspPoint
min :: TspPoint -> TspPoint -> TspPoint
Ord)

instance Show TspPoint where show :: TspPoint -> String
show (TspPoint Int
p) = Int -> String
forall a. Show a => a -> String
show Int
p

-- | How much would the 2-opt switch change the total path length?
--
-- >>> points = [zero, Vec2 100 100, Vec2 0 100, Vec2 100 0]
-- >>> dm = distances (V.fromList points)
-- >>> tsp2optFlipDelta dm (TspPoint 0) (TspPoint 1) (TspPoint 2) (TspPoint 3)
-- -82.84271247461902
--
-- >>> points = [zero, Vec2 0 100, Vec2 100 100, Vec2 100 0]
-- >>> dm = distances (V.fromList points)
-- >>> tsp2optFlipDelta dm (TspPoint 0) (TspPoint 1) (TspPoint 2) (TspPoint 3)
-- 82.84271247461902
tsp2optFlipDelta :: Distances -> TspPoint -> TspPoint -> TspPoint -> TspPoint -> Double
tsp2optFlipDelta :: Distances -> TspPoint -> TspPoint -> TspPoint -> TspPoint -> Double
tsp2optFlipDelta Distances
dm TspPoint
a TspPoint
b TspPoint
c TspPoint
d =
    let
        -- Getting these indices right cost me quite a bit of headache. Now I know
        -- I should just have written them down using pen and paper and they become
        -- quite obvious.
        --
        -- Original: a -> b -> c -> d -> a
        -- Flipped:  a -> c -> b -> d -> a
        --
        -- a <---------- d           a <---------- d
        -- |             ^            '.         .^
        -- |             |              '.     .'
        -- |             |                '. .'
        -- |   Before    |   ===>    After  X
        -- |             |                .' '.
        -- |             |              .'     '.
        -- v             |            .'         'v
        -- b ----------> c           b ----------> c
        oldLength :: Double
oldLength = Distances -> TspPoint -> TspPoint -> Double
distance Distances
dm TspPoint
a TspPoint
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Distances -> TspPoint -> TspPoint -> Double
distance Distances
dm TspPoint
c TspPoint
d
        newLength :: Double
newLength = Distances -> TspPoint -> TspPoint -> Double
distance Distances
dm TspPoint
a TspPoint
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Distances -> TspPoint -> TspPoint -> Double
distance Distances
dm TspPoint
b TspPoint
d
    in Double
newLength Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
oldLength

tsp3optFlipDelta :: Distances -> TspPoint -> TspPoint -> TspPoint -> TspPoint -> TspPoint -> TspPoint -> Vector Double
tsp3optFlipDelta :: Distances
-> TspPoint
-> TspPoint
-> TspPoint
-> TspPoint
-> TspPoint
-> TspPoint
-> Vector Double
tsp3optFlipDelta Distances
dm TspPoint
a TspPoint
b TspPoint
c TspPoint
d TspPoint
e TspPoint
f =
    let dist :: TspPoint -> TspPoint -> Double
dist = Distances -> TspPoint -> TspPoint -> Double
distance Distances
dm
        oldLength :: Double
oldLength = TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
c TspPoint
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
e TspPoint
f
        -- I don’t think I’ll paint the ASCII art for 7 diagrams here ugh.
    in [Double] -> Vector Double
forall a. [a] -> Vector a
V.fromList
        [ Double
0                                 -- Original: { a b c d e f a }
        , TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
e Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
b TspPoint
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- TspPoint -> TspPoint -> Double
dist TspPoint
e TspPoint
f   -- { a e d c b f a }
        , TspPoint -> TspPoint -> Double
dist TspPoint
c TspPoint
e Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
d TspPoint
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- TspPoint -> TspPoint -> Double
dist TspPoint
c TspPoint
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- TspPoint -> TspPoint -> Double
dist TspPoint
e TspPoint
f   -- { a b c e d f a }
        , TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
b TspPoint
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- TspPoint -> TspPoint -> Double
dist TspPoint
c TspPoint
d   -- { a c b d e f a }
        , TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
b TspPoint
e Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
d TspPoint
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
oldLength  -- { a c b e d f a }
        , TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
e Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
d TspPoint
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
c TspPoint
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
oldLength  -- { a e d b c f a }
        , TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
e TspPoint
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
b TspPoint
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
oldLength  -- { a d e c b f a }
        , TspPoint -> TspPoint -> Double
dist TspPoint
a TspPoint
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
e TspPoint
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ TspPoint -> TspPoint -> Double
dist TspPoint
c TspPoint
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
oldLength  -- { a d e b c f a }
        ]

-- | Inplace reverse part of a 'STVector'.
--
-- This is more complicated than the standard »swap around the middle« reversal
-- algorithm, because we wrap around the ends of the vector for j<i, which happens
-- when the swapping window starts at the end of the input. This is not a problem
-- for 2-gen, but in 3-gen we cannot simply reverse »the internal half«, because
-- there might be multiple overlapping swapping regions.
--
-- >>> input = V.fromList [0..10]
-- >>> V.modify (tspReverseInplace 3 6) input
-- [0,1,2,6,5,4,3,7,8,9,10]
--
-- V.modify (tspReverseInplace 8 3) input
-- [0,10,9,8,4,5,6,7,3,2,1]
--
-- >>> V.modify (tspReverseInplace 8 3) (V.init input)
-- [1,0,9,8,4,5,6,7,3,2]
--
-- >>> V.modify (tspReverseInplace 6 6) input
-- [0,1,2,3,4,5,6,7,8,9,10]
tspReverseInplace :: PrimMonad m => Int -> Int -> VM.MVector (PrimState m) a -> m ()
tspReverseInplace :: forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace Int
i Int
j MVector (PrimState m) a
vec = Int -> Int -> Int -> m ()
forall {f :: * -> *} {t}.
(PrimState f ~ PrimState m, Eq t, Num t, PrimMonad f) =>
t -> Int -> Int -> f ()
go Int
steps Int
i Int
j
  where
    n :: Int
n = MVector (PrimState m) a -> Int
forall s a. MVector s a -> Int
VM.length MVector (PrimState m) a
vec
    steps :: Int
steps = ((Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int -> Int -> Int
% Int
n) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 -- Why the +1? Beats me. I spent way too much time on this tiny alg already.
    go :: t -> Int -> Int -> f ()
go t
s Int
_ Int
_ | t
s t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== t
0 = () -> f ()
forall a. a -> f a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ()
    go t
s Int
ii Int
jj = do
        MVector (PrimState f) a -> Int -> Int -> f ()
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> Int -> m ()
VM.swap MVector (PrimState m) a
MVector (PrimState f) a
vec (Int
iiInt -> Int -> Int
%Int
n) (Int
jjInt -> Int -> Int
%Int
n)
        t -> Int -> Int -> f ()
go (t
st -> t -> t
forall a. Num a => a -> a -> a
-t
1) (Int
iiInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
jjInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)

-- | How much would the path length change if the two segments were swapped?
--
-- Does not modify the input vector, but marked as @unsafe@ because out of bounds
-- will 'error'.
--
-- For a bit more on the args’ constraints, see 'unsafeCommit3optSwap'.
unsafePlan2optSwap
    :: PrimMonad m
    => Distances -- \(n\times n\) matrix
    -> Int -- ^ \(i \in [0\ldots n-3]\)
    -> Int -- ^ \(j \in [i+2\ldots n-1]\)
    -> VM.MVector (PrimState m) TspPoint -- ^ \(n\) entries
    -> m Double -- ^ Length delta. Negative means the swap shortens the path.
unsafePlan2optSwap :: forall (m :: * -> *).
PrimMonad m =>
Distances
-> Int -> Int -> MVector (PrimState m) TspPoint -> m Double
unsafePlan2optSwap Distances
dm Int
i Int
j MVector (PrimState m) TspPoint
path = do
    let n :: Int
n = Distances -> Int
size Distances
dm
    TspPoint
a <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path Int
i
    TspPoint
b <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n)
    TspPoint
c <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path Int
j
    TspPoint
d <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n)
    Double -> m Double
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Distances -> TspPoint -> TspPoint -> TspPoint -> TspPoint -> Double
tsp2optFlipDelta Distances
dm TspPoint
a TspPoint
b TspPoint
c TspPoint
d)

-- | Perform a 2-opt swap as specified.
--
-- The requirements on \(i,j\) are not checked, hence @unsafe@. For details see
-- 'unsafeCommit3optSwap'.
unsafeCommit2optSwap
    :: PrimMonad m
    => Int -- ^ \(i \in [0\ldots n-3]\)
    -> Int -- ^ \(j \in [i+2\ldots n-1]\)
    -> VM.MVector (PrimState m) a -- ^ \(n\) entries
    -> m ()
unsafeCommit2optSwap :: forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
unsafeCommit2optSwap Int
i Int
j MVector (PrimState m) a
path = Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int
j MVector (PrimState m) a
path

-- | Inplace version of 'tsp2optPath'.
tsp2optInplace :: PrimMonad m => Distances -> VM.MVector (PrimState m) TspPoint -> m ()
tsp2optInplace :: forall (m :: * -> *).
PrimMonad m =>
Distances -> MVector (PrimState m) TspPoint -> m ()
tsp2optInplace Distances
dm MVector (PrimState m) TspPoint
path = m ()
loop0
  where
    n :: Int
n = Distances -> Int
size Distances
dm
    loop0 :: m ()
loop0 = Int -> Int -> m ()
loop Int
0 Int
2
    loop :: Int -> Int -> m ()
loop Int
i Int
j
        | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3 = () -> m ()
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ()          -- i exhausted, search done
        | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 = Int -> Int -> m ()
loop (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) -- j exhausted, increase i
    loop Int
i Int
j = do
        Double
delta <- Distances
-> Int -> Int -> MVector (PrimState m) TspPoint -> m Double
forall (m :: * -> *).
PrimMonad m =>
Distances
-> Int -> Int -> MVector (PrimState m) TspPoint -> m Double
unsafePlan2optSwap Distances
dm Int
i Int
j MVector (PrimState m) TspPoint
path
        if Double
delta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0
            then Int -> Int -> MVector (PrimState m) TspPoint -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
unsafeCommit2optSwap Int
i Int
j MVector (PrimState m) TspPoint
path m () -> m () -> m ()
forall a b. m a -> m b -> m b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> m ()
loop0
            else Int -> Int -> m ()
loop Int
i (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)

-- | Swapping cases of 3-opt. Created by 'unsafePlan3optSwap', consumed by
-- 'unsafeCommit3optSwap'.
data Case3opt
    = Case3opt1
    | Case3opt2
    | Case3opt3
    | Case3opt4
    | Case3opt5
    | Case3opt6
    | Case3opt7
    deriving (Case3opt -> Case3opt -> Bool
(Case3opt -> Case3opt -> Bool)
-> (Case3opt -> Case3opt -> Bool) -> Eq Case3opt
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Case3opt -> Case3opt -> Bool
== :: Case3opt -> Case3opt -> Bool
$c/= :: Case3opt -> Case3opt -> Bool
/= :: Case3opt -> Case3opt -> Bool
Eq, Eq Case3opt
Eq Case3opt
-> (Case3opt -> Case3opt -> Ordering)
-> (Case3opt -> Case3opt -> Bool)
-> (Case3opt -> Case3opt -> Bool)
-> (Case3opt -> Case3opt -> Bool)
-> (Case3opt -> Case3opt -> Bool)
-> (Case3opt -> Case3opt -> Case3opt)
-> (Case3opt -> Case3opt -> Case3opt)
-> Ord Case3opt
Case3opt -> Case3opt -> Bool
Case3opt -> Case3opt -> Ordering
Case3opt -> Case3opt -> Case3opt
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Case3opt -> Case3opt -> Ordering
compare :: Case3opt -> Case3opt -> Ordering
$c< :: Case3opt -> Case3opt -> Bool
< :: Case3opt -> Case3opt -> Bool
$c<= :: Case3opt -> Case3opt -> Bool
<= :: Case3opt -> Case3opt -> Bool
$c> :: Case3opt -> Case3opt -> Bool
> :: Case3opt -> Case3opt -> Bool
$c>= :: Case3opt -> Case3opt -> Bool
>= :: Case3opt -> Case3opt -> Bool
$cmax :: Case3opt -> Case3opt -> Case3opt
max :: Case3opt -> Case3opt -> Case3opt
$cmin :: Case3opt -> Case3opt -> Case3opt
min :: Case3opt -> Case3opt -> Case3opt
Ord, Int -> Case3opt -> ShowS
[Case3opt] -> ShowS
Case3opt -> String
(Int -> Case3opt -> ShowS)
-> (Case3opt -> String) -> ([Case3opt] -> ShowS) -> Show Case3opt
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Case3opt -> ShowS
showsPrec :: Int -> Case3opt -> ShowS
$cshow :: Case3opt -> String
show :: Case3opt -> String
$cshowList :: [Case3opt] -> ShowS
showList :: [Case3opt] -> ShowS
Show, Int -> Case3opt
Case3opt -> Int
Case3opt -> [Case3opt]
Case3opt -> Case3opt
Case3opt -> Case3opt -> [Case3opt]
Case3opt -> Case3opt -> Case3opt -> [Case3opt]
(Case3opt -> Case3opt)
-> (Case3opt -> Case3opt)
-> (Int -> Case3opt)
-> (Case3opt -> Int)
-> (Case3opt -> [Case3opt])
-> (Case3opt -> Case3opt -> [Case3opt])
-> (Case3opt -> Case3opt -> [Case3opt])
-> (Case3opt -> Case3opt -> Case3opt -> [Case3opt])
-> Enum Case3opt
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
$csucc :: Case3opt -> Case3opt
succ :: Case3opt -> Case3opt
$cpred :: Case3opt -> Case3opt
pred :: Case3opt -> Case3opt
$ctoEnum :: Int -> Case3opt
toEnum :: Int -> Case3opt
$cfromEnum :: Case3opt -> Int
fromEnum :: Case3opt -> Int
$cenumFrom :: Case3opt -> [Case3opt]
enumFrom :: Case3opt -> [Case3opt]
$cenumFromThen :: Case3opt -> Case3opt -> [Case3opt]
enumFromThen :: Case3opt -> Case3opt -> [Case3opt]
$cenumFromTo :: Case3opt -> Case3opt -> [Case3opt]
enumFromTo :: Case3opt -> Case3opt -> [Case3opt]
$cenumFromThenTo :: Case3opt -> Case3opt -> Case3opt -> [Case3opt]
enumFromThenTo :: Case3opt -> Case3opt -> Case3opt -> [Case3opt]
Enum)

-- | If an improvement can be made by 3 segments, which case yields the largest
-- improvement, and by what amount?
--
-- Does not modify the input vector, but marked as @unsafe@ because out of bounds
-- will 'error'.
--
-- For a bit more on the args’ constraints, see 'unsafeCommit3optSwap'.
unsafePlan3optSwap
    :: PrimMonad m
    => Distances -- ^ \(n\times n\) matrix
    -> Int -- ^ \(i \in [0\ldots n-5]\)
    -> Int -- ^ \(j \in [i+2\ldots n-3]\)
    -> Int -- ^ \(k \in [j+2\ldots n-1]\)
    -> VM.MVector (PrimState m) TspPoint -- ^ \(n\) entries
    -> m (Maybe (Case3opt, Double)) -- ^ Best case, and the (negative) length delta compared to the old path.
unsafePlan3optSwap :: forall (m :: * -> *).
PrimMonad m =>
Distances
-> Int
-> Int
-> Int
-> MVector (PrimState m) TspPoint
-> m (Maybe (Case3opt, Double))
unsafePlan3optSwap Distances
dm Int
i Int
j Int
k MVector (PrimState m) TspPoint
path = do
    let n :: Int
n = Distances -> Int
size Distances
dm
    TspPoint
a <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path Int
i
    TspPoint
b <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n)
    TspPoint
c <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path Int
j
    TspPoint
d <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n)
    TspPoint
e <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path Int
k
    TspPoint
f <- MVector (PrimState m) TspPoint -> Int -> m TspPoint
forall (m :: * -> *) a.
PrimMonad m =>
MVector (PrimState m) a -> Int -> m a
VM.read MVector (PrimState m) TspPoint
path (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n)

    let cases :: Vector Double
cases = Distances
-> TspPoint
-> TspPoint
-> TspPoint
-> TspPoint
-> TspPoint
-> TspPoint
-> Vector Double
tsp3optFlipDelta Distances
dm TspPoint
a TspPoint
b TspPoint
c TspPoint
d TspPoint
e TspPoint
f
        bestCase :: Int
bestCase = Vector Double -> Int
forall a. Ord a => Vector a -> Int
V.minIndex Vector Double
cases
        bestDelta :: Double
bestDelta = Vector Double
casesVector Double -> Int -> Double
forall a. Vector a -> Int -> a
!Int
bestCase

    Maybe (Case3opt, Double) -> m (Maybe (Case3opt, Double))
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Maybe (Case3opt, Double) -> m (Maybe (Case3opt, Double)))
-> Maybe (Case3opt, Double) -> m (Maybe (Case3opt, Double))
forall a b. (a -> b) -> a -> b
$ case Int
bestCase of
        Int
0 -> Maybe (Case3opt, Double)
forall a. Maybe a
Nothing
        Int
_ | Double
bestDelta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> -Double
1e-8 -> Maybe (Case3opt, Double)
forall a. Maybe a
Nothing -- To avoid flaky super-tiny improvements mostly introduced by numerical instabilities
        Int
_otherwise -> (Case3opt, Double) -> Maybe (Case3opt, Double)
forall a. a -> Maybe a
Just (Int -> Case3opt
forall a. Enum a => Int -> a
toEnum (Int
bestCaseInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), Double
bestDelta)

-- | Perform a 3-opt swap as specified.
--
-- The requirements on \(i,j,k\) are not checked, hence @unsafe@. Their purpose is
-- making a minimal search of the problem space possible, exploiting all symmetries
-- available by a symmetrical distance function and the 3-opt swaps:
--
--  * We can choose \(i<j<k\) WLOG, because 3-opt swaps are 3-symmetric.
--  * \(i\) needs to leave space for \(i+1,j,j+1,k,k+1\), hence the \(n-5\).
--  * Similarly, \(j\) needs to leave space for \(j+1,k,k+1\), hence the \(n-3\).
--  * \(k\) can point to the very end, leaving \(k+1\) to cycle around to the
--    beginning.
unsafeCommit3optSwap
    :: PrimMonad m
    => Case3opt
    -> Int -- ^ \(i \in [0\ldots n-5]\)
    -> Int -- ^ \(j \in [i+2\ldots n-3]\)
    -> Int -- ^ \(k \in [j+2\ldots n-1]\)
    -> Int -- ^ \(n\)
    -> VM.MVector (PrimState m) a -- ^ \(n\) points
    -> m ()
unsafeCommit3optSwap :: forall (m :: * -> *) a.
PrimMonad m =>
Case3opt
-> Int -> Int -> Int -> Int -> MVector (PrimState m) a -> m ()
unsafeCommit3optSwap Case3opt
bestCase Int
i Int
j Int
k Int
n MVector (PrimState m) a
path = case Case3opt
bestCase of
    Case3opt
Case3opt1 -> do
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n) Int
i MVector (PrimState m) a
path
    Case3opt
Case3opt2 -> do
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1    ) Int
k MVector (PrimState m) a
path
    Case3opt
Case3opt3 -> do
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1    ) Int
j MVector (PrimState m) a
path
    Case3opt
Case3opt4 -> do
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1    ) Int
k MVector (PrimState m) a
path
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1    ) Int
j MVector (PrimState m) a
path
    Case3opt
Case3opt5 -> do
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n) Int
i MVector (PrimState m) a
path
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1    ) Int
j MVector (PrimState m) a
path
    Case3opt
Case3opt6 -> do
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n) Int
i MVector (PrimState m) a
path
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1    ) Int
k MVector (PrimState m) a
path
    Case3opt
Case3opt7 -> do
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Int
% Int
n) Int
i MVector (PrimState m) a
path
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1    ) Int
j MVector (PrimState m) a
path
        Int -> Int -> MVector (PrimState m) a -> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Int -> Int -> MVector (PrimState m) a -> m ()
tspReverseInplace (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1    ) Int
k MVector (PrimState m) a
path

-- | Inplace version of 'tsp3optPath'.
tsp3optInplace :: PrimMonad m => Distances -> VM.MVector (PrimState m) TspPoint -> m ()
tsp3optInplace :: forall (m :: * -> *).
PrimMonad m =>
Distances -> MVector (PrimState m) TspPoint -> m ()
tsp3optInplace Distances
dm MVector (PrimState m) TspPoint
path = m ()
loop0
  where
    n :: Int
n = Distances -> Int
size Distances
dm
    loop0 :: m ()
loop0 = Int -> Int -> Int -> m ()
loop Int
0 Int
2 Int
4
    loop :: Int -> Int -> Int -> m ()
loop Int
i Int
j Int
k
        | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
5 = () -> m ()
forall a. a -> m a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ()                -- i exhausted, search done
        | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3 = Int -> Int -> Int -> m ()
loop (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5) -- j exhausted, increase i
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 = Int -> Int -> Int -> m ()
loop  Int
i    (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) -- k exhausted, increase j
    loop Int
i Int
j Int
k = do
        Maybe (Case3opt, Double)
result <- Distances
-> Int
-> Int
-> Int
-> MVector (PrimState m) TspPoint
-> m (Maybe (Case3opt, Double))
forall (m :: * -> *).
PrimMonad m =>
Distances
-> Int
-> Int
-> Int
-> MVector (PrimState m) TspPoint
-> m (Maybe (Case3opt, Double))
unsafePlan3optSwap Distances
dm Int
i Int
j Int
k MVector (PrimState m) TspPoint
path
        case Maybe (Case3opt, Double)
result of
            Maybe (Case3opt, Double)
Nothing -> Int -> Int -> Int -> m ()
loop Int
i Int
j (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
            Just (Case3opt
bestCase, Double
_delta) -> do
                Case3opt
-> Int
-> Int
-> Int
-> Int
-> MVector (PrimState m) TspPoint
-> m ()
forall (m :: * -> *) a.
PrimMonad m =>
Case3opt
-> Int -> Int -> Int -> Int -> MVector (PrimState m) a -> m ()
unsafeCommit3optSwap Case3opt
bestCase Int
i Int
j Int
k Int
n MVector (PrimState m) TspPoint
path
                m ()
loop0

-- | Move towards a local minimum using the 2-opt algorithm: delete 2 edges, and if
-- the reconnected version is shorter than the original. The result of this algorithm
-- is guaranteed to be at most as large as its input.
--
-- Starting this algoritm from a path optimized by 'tspNearestNeighbourPath' speeds
-- it up considerably, although this usually yields a different but similarly good
-- result than applying it directly.
--
-- <<docs/haddock/Numerics/Optimization/TSP/2-opt.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Numerics/Optimization/TSP/2-opt.svg" 300 200 $ \_ -> do
--    points <- liftIO $ do
--        gen <- MWC.initialize (V.fromList [])
--        uniform <- uniformlyDistributedPoints gen (shrinkBoundingBox 10 [zero, Vec2 300 200]) 32
--        pure uniform
--    setLineJoin LineJoinRound
--    let dm = distances points
--        path = tsp2optPath dm (tspNearestNeighbourPath dm)
--        tour = Polygon (V.toList (pathToPoints points path))
--    cairoScope $ do
--        setColor (mma 0)
--        for_ points $ \p -> sketch (Circle p 3) >> fill
--    cairoScope $ do
--        setColor (mma 3)
--        sketch tour
--        stroke
-- :}
-- Generated file: size 11KB, crc32: 0x73c4af54
tsp2optPath :: Distances -> Vector TspPoint -> Vector TspPoint
tsp2optPath :: Distances -> Vector TspPoint -> Vector TspPoint
tsp2optPath Distances
dm = (forall s. MVector s TspPoint -> ST s ())
-> Vector TspPoint -> Vector TspPoint
forall a.
(forall s. MVector s a -> ST s ()) -> Vector a -> Vector a
V.modify (Distances -> MVector (PrimState (ST s)) TspPoint -> ST s ()
forall (m :: * -> *).
PrimMonad m =>
Distances -> MVector (PrimState m) TspPoint -> m ()
tsp2optInplace Distances
dm)

-- | Move towards a local minimum using the 3-opt algorithm: delete 3 edges, and if
-- the reconnected version is shorter than the original. The result of this
-- algorithm is guaranteed to be at most as large as its input.
--
-- 3-opt yields higher quality results than 2-opt, but is considerably slower. It
-- is highly advised to only run it on already optimized paths, e.g. starting from
-- 'tspNearestNeighbourPath' or 'tsp2optPath'.
--
-- <<docs/haddock/Numerics/Optimization/TSP/3-opt.svg>>
--
-- === __(image code)__
-- >>> :{
-- haddockRender "Numerics/Optimization/TSP/3-opt.svg" 300 200 $ \_ -> do
--    points <- liftIO $ do
--        gen <- MWC.initialize (V.fromList [])
--        uniform <- uniformlyDistributedPoints gen (shrinkBoundingBox 10 [zero, Vec2 300 200]) 32
--        pure uniform
--    setLineJoin LineJoinRound
--    let dm = distances points
--        path = tsp3optPath dm (tspNearestNeighbourPath dm)
--        tour = Polygon (V.toList (pathToPoints points path))
--    cairoScope $ do
--        setColor (mma 0)
--        for_ points $ \p -> sketch (Circle p 3) >> fill
--    cairoScope $ do
--        setColor (mma 4)
--        sketch tour
--        stroke
-- :}
-- Generated file: size 11KB, crc32: 0x6c4e9493
tsp3optPath :: Distances -> Vector TspPoint -> Vector TspPoint
tsp3optPath :: Distances -> Vector TspPoint -> Vector TspPoint
tsp3optPath Distances
dm = (forall s. MVector s TspPoint -> ST s ())
-> Vector TspPoint -> Vector TspPoint
forall a.
(forall s. MVector s a -> ST s ()) -> Vector a -> Vector a
V.modify (Distances -> MVector (PrimState (ST s)) TspPoint -> ST s ()
forall (m :: * -> *).
PrimMonad m =>
Distances -> MVector (PrimState m) TspPoint -> m ()
tsp3optInplace Distances
dm)