Finding minimum scalar product using ST Monad Inspired by a Stack Overflow question, I decided to practice my Haskell by taking a crack at the Google Code Jam's Minimum Scalar Product problem:


  Given two vectors \$\mathbf{v_1}=(x_1,x_2,\ldots,x_n)\$ and \$\mathbf{v_2}=(y_1,y_2,\ldots,y_n)\$. If you could permute the coordinates of each vector, what is the minimum scalar product \$\mathbf{v_1} \cdot \mathbf{v_2}\$?
  
  Constraints:
  \$100 \le n \le 800\$
  \$-100000 \le x_i, y_i \le 100000\$


I'm not claiming any algorithmic awesomeness (this is just a reference implementation for checking correctness later).

This is my first time using Vectors and the ST monad, so what I really want is a sanity check that I'm using both correctly, and that I'm using the correct tools for the job.

module MinimumScalarProduct where

import Control.Monad (replicateM, forM_, (&gt;=&gt;))
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 -&gt; [(Int, Int)]
transpositions n = runST $ do
          -- p maps index to element
          p &lt;- thaw $ fromList [0 .. n-1]
          -- q maps element to index
          q &lt;- thaw $ fromList [0 .. n-1]
          -- let the prefixes define themselves recursively
          foldr ((&gt;=&gt;) . extend p q) return [2..n] [] 

extend :: MVector s Int -&gt; MVector s Int -&gt; Int -&gt; [(Int,Int)] -&gt; 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 &lt;- read p (n-1)
  i &lt;- read q (a-1)
  swap p (n-1) i
  swap q a (a-1)
  -- replay the earlier transpositions
  forM_ ts $ \(m,j) -&gt; do
    b &lt;- read p m
    c &lt;- read p j
    swap p m j
    swap q b c
  return $ (n-1,i) : ts

-- reference implementation, takes O(n!)
msp :: Vector Int -&gt; Vector Int -&gt; 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' &lt;- thaw u
  -- check each permutation of u'
  let steps = map (adjust u' v) $ transpositions n
  fmap minimum . sequence $ scanl (&gt;&gt;=) (return x) steps

-- adjust the current scalar product for the transposition of u
adjust :: MVector s Int -&gt; Vector Int -&gt; (Int, Int) -&gt; Int -&gt; ST s Int
adjust u v (i,j) x = do
  a &lt;- read u i
  b &lt;- 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)
Somewhat old question, but still it's a pity it hasn't been answered :). Putting aside the asymptotically better solution suggested in the comments, and focusing just on the code:

Pluses: Top-level types for all funtcions, that's very good, as well as comments. The code shows good understanding of various Haskell functions and idioms.

My most general comment is that there is disproportion between code verbosity and complexity. Usually more complex code should be more verbose and more descriptive. Visually, most of the code is just simple read/write/swap ST operations, but the important part is squeezed into not-so-easy-to-understand and uncommon combinations of folds/scans/sequence/>=>/>>= etc. I'd probably separate the simple parts (read/swap) into helper functions and expand/simplify the complex parts a bit.

Function adjust somewhat mixes pure and mutable approaches. While it mutates its first argument, it keeps the scalar product pure and passes it as an argument in and out. This is somewhat confusing and complicates the computation of the minimum (scanl with &gt;&gt;=). If the product were a ST variable too, the code gets simpler (untested):

adjust :: MVector s Int -&gt; Vector Int -&gt; STRef s Int -&gt; (Int, Int) -&gt; ST s Int
...
msp = ...
  x &lt;- newSTRef $ foldl' (+) 0 $ zipWith (*) u v
  fmap minimum . traverse (adjust u' v x) . transpositions $ n


Alternatively you could even keep the result in an ST variable.

Also foldr with &gt;=&gt; can be replaced with foldM (untested):

foldM (flip $ extend p q) [] [2..n]


(Then it'd make sense to refactor the type of extend to get rid of the flip.)

Some questions:


Is replaying all ts in extend necessary? As it's executed repeatedly, it'd make sense to compute the combined permutation of ts separately or once and then just apply it.
Which variant of the algorithm is actually used? Using two vectors doesn't seem to match anything in the linked Wikipedia article. While the general idea of the algorithm is clear, it'd be easier for someone not familiar with it to get more precise guidelines.


Nit: The order of functions is somewhat unnatural: transpositions uses extend, msp uses transpositions and adjust. There is no visible order in this. One good option is that functions are always defined before being used (with the obvious exception of mutually-recursive functions). Or the other way around - the main/exported functions are first and helper functions follow. Or, separate functions into (possibly nested) sections, then ordering of functions doesn't really matter, ordering of sections becomes important.