使用ST Monad查找最小标量产品

module MinimumScalarProduct where

import Control.Monad (replicateM, forM_, (>=>))
import Control.Monad.ST (runST, ST)
import Data.Vector (thaw, fromList, MVector, Vector, zipWith, foldl', length, (!))
import Data.Vector.Generic.Mutable (read, swap)
import Prelude hiding (read, zipWith, length)

-- sequnce of transpoitions to yield all permutations of n elts
-- http://en.wikipedia.org/wiki/Steinhaus-Johnson-Trotter_algorithm
transpositions :: Int -> [(Int, Int)]
transpositions n = runST $ do
          -- p maps index to element
          p <- thaw $ fromList [0 .. n-1]
          -- q maps element to index
          q <- thaw $ fromList [0 .. n-1]
          -- let the prefixes define themselves recursively
          foldr ((>=>) . extend p q) return [2..n] [] 

extend :: MVector s Int -> MVector s Int -> Int -> [(Int,Int)] -> ST s [(Int,Int)]
extend p q n ts = fmap ((ts++) . concat) . replicateM (n-1) $ do
  -- replace the element in the (n-1)'th position with its predecessor
  a <- read p (n-1)
  i <- read q (a-1)
  swap p (n-1) i
  swap q a (a-1)
  -- replay the earlier transpositions
  forM_ ts $ \(m,j) -> do
    b <- read p m
    c <- read p j
    swap p m j
    swap q b c
  return $ (n-1,i) : ts

-- reference implementation, takes O(n!)
msp :: Vector Int -> Vector Int -> Int
msp u v | length u /= length v = 0
msp u v = runST $ do
  let x = foldl' (+) 0 $ zipWith (*) u v
  let n = length u
  u' <- thaw u
  -- check each permutation of u'
  let steps = map (adjust u' v) $ transpositions n
  fmap minimum . sequence $ scanl (>>=) (return x) steps

-- adjust the current scalar product for the transposition of u
adjust :: MVector s Int -> Vector Int -> (Int, Int) -> Int -> ST s Int
adjust u v (i,j) x = do
  a <- read u i
  b <- read u j
  let c = v ! i
  let d = v ! j
  swap u i j
  return $ x - (a*c + b*d) + (a*d + b*c)