# 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 Data.Vector (thaw, fromList, MVector, Vector, zipWith, foldl', length, (!))
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
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
let c = v ! i
let d = v ! j
swap u i j
return $x - (a*c + b*d) + (a*d + b*c)  • You do know that this problem can be solved by just sorting the vectors? Why do it in such a complicated way when it's just a reference implementation, what's wrong with Data.List.permutations? Feb 21, 2015 at 14:26 • Hjulle: it's been three years, but I think I was trying to minimize the work to compute the product for all the permutations by walking through them in transposition order. So it's O(1) work per vector product, rather than O(n), making the total algorithm O(n!) rather than O(n*n!). In retrospect, perhaps not a great savings. Feb 21, 2015 at 15:30 • @Hjulle For such a public challenge I'd suggest to update your comment not to disclose the solution as it ruins the challenge for others as well as for the OP. It'd be better to review the OP's code wrt to the idea (s)he's expressing, perhaps just mentioning something like "Note that there is a much more asymptotically efficient algorithm" and give the OP (and others) a chance to find it on their own. – Petr Jun 28, 2015 at 12:13 ## 1 Answer 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 >>=). If the product were a ST variable too, the code gets simpler (untested): adjust :: MVector s Int -> Vector Int -> STRef s Int -> (Int, Int) -> ST s Int ... msp = ... x <- 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 >=> 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.