If we wrote that code in Haskell it would read a bit like this:

data Row = Row { shortListed :: Bool, cost :: Float, distance1 :: Int, distance2 :: Int } deriving (Show, Eq)

import Data.Ord
import Data.List
 
compareRow :: Row -> Row -> Ordering
compareRow x y = if comparing (not . shortListed) x y == EQ then 
                    if comparing cost x y == EQ then
                      if comparing distance1 x y == EQ then
                        comparing distance2 x y
                      else
                        comparing distance1 x y
                    else
                      comparing cost x y
                  else 
                    comparing (not . shortListed) x y

We can then sort a list of rows like this:

> let shortListedRow = Row {shortListed = True, cost = 10.0, distance1 = 20, distance2 = 30 }
 
> let nonShortListedRow = Row {shortListed = False, cost = 10.0, distance1 = 20, distance2 = 30 }
 
> sortBy compareRow [nonShortListedRow, shortListedRow]
[Row {shortListed = True, cost = 10.0, distance1 = 20, distance2 = 30},
 Row {shortListed = False, cost = 10.0, distance1 = 20, distance2 = 30}]

It works but it’s messy and we couldn’t see what abstraction we should be using to simplify the code.

I was continuing with my reading of Functors, Applicative Functors and Monoids yesterday and got to the section on Monoids which showed an example for simplifying this type of code.

The definition of a Monoid from the Haskell source code is:

Types with an associative binary operation that has an identity

But I prefer Dan Piponi’s definition:

In Haskell, a monoid is a type with a rule for how two elements of that type can be combined to make another element of the same type.

To be a monoid there also needs to be an element that you can think of as representing ‘nothing’ in the sense that when it’s combined with other elements it leaves the other element unchanged.

In our case we have a bunch of things of type ‘Ordering’ and we want to combine them all together and end up with a final ‘Ordering’ which takes them all into account.

For example if we were comparing the following two rows:

> let row1 = Row {shortListed = True, cost = 10.0, distance1 = 1, distance2 = 30 }
> let row2 = Row {shortListed = True, cost = 10.0, distance1 = 100, distance2 = 30 }
> compareRow row1 row2
LT

When we compare their shortListed value we get back ‘EQ’, so we compare their cost value and get back ‘EQ’ and finally we compare their distance1 value which gives back ‘LT’ which is our final value.

We can make use of the Ordering Monoid to do this rather than all the nested if statements.

Monoid is a type class defined like so:

class Monoid a where
  mempty  :: a
  mappend :: a -> a -> a
  mconcat :: [a] -> a
  mconcat = foldr mappend mempty

‘mempty’ represents the identity value for a monoid i.e. the ‘nothing’ in Dan Piponi’s definition. If we combine anything with this then we should get that thing back.

The most interesting function here is ‘mappend’ which we use to combine together two elements of a type. Each instance of Monoid needs to define this function for themselves.

The Ordering Monoid is defined like so:

instance Monoid Ordering where
  mempty         = EQ
  LT `mappend` _ = LT
  EQ `mappend` y = y
  GT `mappend` _ = GT

What makes this work for us is that we always keep the value on the left unless it’s ‘EQ’ in which case we take the value on the right.

Therefore as soon as one of our comparisons returns a non ‘EQ’ value that will be the value that eventually gets returned.

e.g.

> GT `mappend` LT `mappend` EQ
GT

Our ‘row1′/’row2′ comparison would look like this using ‘mappend’:

> EQ `mappend` EQ `mappend` LT
LT

We can then change our ‘compareRow’ function:

compareRow x y = comparing (not . shortListed) x y `mappend` comparing cost x y `mappend` comparing distance1 x y `mappend` comparing distance2 x y

We can simplify this further by making use of ‘mconcat’ which folds over a list of monoids applying ‘mappend’ each time.

For example we could replace our ‘row1′/’row2′ comparison with the following:

> mconcat [EQ, EQ, LT]
LT

And ‘compareRow’ now reads like this:

compareRow x y = mconcat [comparing (not . shortListed) x y,  comparing cost x y, comparing distance1 x y, comparing distance2 x y]

We’re still repeating the ‘comparing’ bit of code every time so I extracted that into a function:

by :: Ord a => (b -> a) -> b -> b -> Ordering
by fn x y = comparing fn x y

We then need to apply those functions to x and y to get our collection of monoids we can pass to mconcat:

compareRow x y = mconcat $ map (\fn -> fn x y) [by (not . shortListed), by cost, by distance1, by distance2]

One problem with this code is that we’re now comparing by all the parameters when we can actually stop once we’ve found a non equality.

I’m sure there’s a clever way to short circuit this but I’m not sure what it is at the moment!

As BeRwet points out in the comments what I wrote here is not actually the case, it is lazily evaluated!

Since the last time I worked in Scala I’ve been playing around with Haskell where type classes are much more common – for example if we want to compare two values we need to make sure that their type extends the ‘Eq’ type class.

Learn Me a Haskell’s chapter on type classes defines them like so:

Typeclasses are like interfaces.

A typeclass defines some behaviour (like comparing for equality, comparing for ordering, enumeration) and then types that can behave in that way are made instances of that typeclass.

The behaviour of typeclasses is achieved by defining functions or just type declarations that we then implement. So when we say that a type is an instance of a typeclass, we mean that we can use the functions that the typeclass defines with that type.

In that chapter we then go on to create a ‘YesNo’ type class which defines Boolean semantics for all different types.

We start by defining the type class like so:

class YesNo a where
  yesno :: a -> Bool

Any type which extends that type class can call this ‘yesno’ function and find out the truthyness of its value.

e.g.

instance YesNo Int where
  yesno 0 = False
  yesno _ = True

If we call that:

> yesno (0 :: Int)
False
 
> yesno (1 :: Int)
True

We get the expected result, but if we try to call it for a type which hasn’t defined an instance of ‘YesNo’:

> yesno "mark"
 
    No instance for (YesNo [Char])
      arising from a use of `yesno'
    Possible fix: add an instance declaration for (YesNo [Char])
    In the expression: yesno "mark"
    In an equation for `it': it = yesno "mark"

In Scala we can use traits and implicits to achieve the same effect. First we define the ‘YesNo’ trait:

trait YesNo[A] {
  def yesno(value:A) : Boolean
}

Then we define an implicit value in a companion object which creates an instance of the ‘YesNo’ type class for Ints:

object YesNo {
  implicit val intYesNo = new YesNo[Int] { 
    def yesno(value:Int) = 
      value match { case 0 => false;  case _ => true } }
}

We then need to call our ‘yesno’ function and the idea is that if we’ve defined a type class instance for the type we call it with it will return us a boolean value and if not then we’ll get a compilation error.

object YesNoWriter {
  def write[A](value:A)(implicit  conv: YesNo[A]) : Boolean = {
    conv.yesno(value)
  }
}

If we call that:

> YesNoWriter.write(1)
res1: Boolean = true
 
> YesNoWriter.write(0)
res2: Boolean = false

It works as expected, but if we try to call it with a type which wasn’t defined for the ‘YesNo’ type class we run into trouble:

> YesNoWriter.write("mark")
:10: error: could not find implicit value for parameter conv: YesNo[java.lang.String]
       YesNoWriter.write("mark")

We can also define YesNoWriter like this by making use of context bounds:

object YesNoWriter {
  def write[A:YesNo](value:A) : Boolean = {
    implicitly[YesNo[A]].yesno(value)
  }
}

I think this pattern is preferred when we might just be tunnelling the implicit parameter through to another method but we can still use it here and use the ‘implicitly’ method to get access to the implicit value.

I’m still not entirely sure about the use of implicits but in this case they provide another way to implement polymorphism without having to use inheritance.

Dan Rosen goes into much more detail about type classes and implicits and I wrote about how we were using them on a project I worked on last year in an earlier blog post if you want to learn more.

I first needed to install test-framework and some dependencies using cabal:

> cabal install test-framework
> cabal install test-framework-quickcheck
> cabal install test-framework-hunit

I thought it’d be interesting to try and write some tests around the windowed function that I wrote a few months ago:

Windowed.hs

module Windowed (windowed) where
 
windowed :: Int -> [a] -> [[a]]
windowed size [] = []
windowed size ls@(x:xs) = 
  if length ls >= size then (take size ls) : windowed size xs 
  else windowed size x

I wrote my first property like so:

import Windowed
 
import Test.Framework (defaultMain, testGroup)
import Test.Framework.Providers.QuickCheck (testProperty)
import Test.QuickCheck
 
main = defaultMain tests
 
tests = [
  testGroup "windowed" [
    testProperty "should not window if n is bigger than list size" prop_nothing_if_n_too_large
  ]
] 
 
prop_nothing_if_n_too_large n xs = n > length xs ==>  windowed n xs == []

And tried to compile the file:

ghc -package test-framework -package test-framework-quickcheck -threaded WindowedQC.hs -o Windowed

Which resulted in this compilation error:

WindowedQC.hs:16:17:
    No instance for (QuickCheck-1.2.0.1:Test.QuickCheck.Testable
                       (Gen Prop))
      arising from a use of `testProperty'
    Possible fix:
      add an instance declaration for
      (QuickCheck-1.2.0.1:Test.QuickCheck.Testable (Gen Prop))

According to the test-framework home page this error happens when we have more than one version of QuickCheck installed so we need to tell ‘ghc’ which one to use:

ghc -package test-framework -package test-framework-quickcheck -package QuickCheck-1.2.0.1 -threaded WindowedQC.hs -o Windowed

That solved the first compilation problem but I still had another:

WindowedQC.hs:16:17:
    Ambiguous type variable `a0' in the constraints:
      (Arbitrary a0) arising from a use of `testProperty'
                     at WindowedQC.hs:12:17-28
      (Show a0) arising from a use of `testProperty'
                at WindowedQC.hs:12:17-28
      (Eq a0) arising from a use of `prop_nothing_if_n_too_large'
              at WindowedQC.hs:12:80-106
    Probable fix: add a type signature that fixes these type variable(s)

The way to solve this problem is to explicitly state which types will be passed to the property like so:

prop_nothing_if_n_too_large n xs = n > length xs ==>  windowed n xs == []
    where types = (xs :: [Int], n :: Int)

I eventually ended up with the following:

WindowedQC.hs

import Windowed
 
import Test.Framework (defaultMain, testGroup)
import Test.Framework.Providers.QuickCheck (testProperty)
import Test.QuickCheck
 
import Test.Framework.Providers.HUnit
import Test.HUnit
 
main = defaultMain tests
 
tests = [
  testGroup "windowed" [
    testProperty "should not window if n is bigger than list size" prop_nothing_if_n_too_large,
    testProperty "should create sub arrays the size of window" prop_size_of_sub_arrays_is_n,
    testProperty "should group adjacent items into windows" prop_groups_adjacent_items,                
    testCase "should window a simple list" test_should_window_simple_list
  ]
] 
 
prop_groups_adjacent_items n xs = n < length xs ==>  
                                  not (null xs) ==>
                                  n > 0 ==> (last $  last $ windowed n xs) == last xs
  where types = (xs :: [Int], n :: Int)			
 
prop_nothing_if_n_too_large n xs = n > length xs ==>  windowed n xs == []
  where types = (xs :: [Int], n :: Int)
 
prop_size_of_sub_arrays_is_n n xs =  n > 0 ==> all (\subArray -> length subArray == n)  (windowed n xs)  
  where types = (xs :: [Int], n :: Int)
 
test_should_window_simple_list = 
  windowed 2 [1..10] @?= [[1,2], [2,3], [3,4], [4,5], [5,6], [6,7], [7,8], [8,9], [9,10]]

I initially only wrote QuickCheck properties but I wasn’t satisfied that just passing the properties I’ve written would guarantee the function is working as expected.

The first property does a simple check that the last item from the windowed function is the same as the last item in the list assuming that the list isn’t empty, that we’re using a positive n value and our window size is smaller than the list size.

The second checks that we end up with an empty list if we try to window to a size bigger than the list and the third checks that each of the sub arrays is the correct size.

I wasn’t sure how to write a QuickCheck property that would test the values of the individual sub arrays so I ended up writing a HUnit test instead.

If you know how I could achieve the same thing using QuickCheck please let me know.

We can then run all the tests like so:

> ./Windowed                                                                                                                    
windowed:
 
  should not window if n is bigger than list size: [OK, passed 100 tests]
 
  should create sub arrays the size of window: [OK, passed 100 tests]
 
  should group adjacent items into windows: [OK, passed 100 tests]
 
  should window a simple list: [OK]
 
          Properties  Test Cases  Total      
  Passed  3           1           4          
  Failed  0           0           0

Initially we tried driving the API through integration style tests but we realised that taking this approach made it quite difficult for us to imagine how an application would use it.

It also meant that we didn’t have a good way to show our business stakeholders the work we’d done – showing off pieces of JSON going back and forth probably wouldn’t have gone down too well!

We therefore decided to create our own test harness UI which we would use to show case what we’d developed so far and also use to help drive the design of the API.

Obviously it would be better if we were able to drive the API out from a real application but given that wasn’t possible the test harness approach seems to work reasonably well.

The main problem that we ran into was not knowing how much effort we should be putting into the test harness UI given that its primary purpose was to show the progress being made on the UI.

We started out not writing many tests or paying much attention to how it looked since it was assumed to be a throwaway UI.

Over time it’s become more complicated and since we use it to display progress to management stakeholders we’ve spent time making it more presentable.

Despite that I’d still say the effort we’ve put in is worthwhile because it gives us a way to show progress to non technical people.

From my experience if you can’t do that they’ll start to get nervous and wonder if you know what you’re doing which isn’t a good place to be!

I started with a simple property to check that the brute force (bf) and divide and conquer (dc) versions of the algorithm returned the same result, assuming that there were enough values in the list to have a closest pair:

prop_closest_pairs xs = length xs >= 2 ==> dcClosest xs == (fromJust $ bfClosest xs)

I could then run that as follows:

> import QuickCheck.Test
> quickCheck(prop_closest_pairs :: [(Double, Double)] -> Property)

It failed pretty quickly because the bf and dc versions of the algorithm sometimes return the pairs in a different order.

e.g. bf may say the closest pair is (2.0, 0.0), (2.1, 1.1) while dc says it’s (2.1, 1.1), (2.0, 0.0) which will lead to the quick check property failing because those values aren’t equal:

> ((2.0, 0.0), (2.1, 1.1))  == ((2.1, 1.1), (2.0, 0.0))
False

The best way I could think of to get around this problem was to create a type to represent a pair of points and then write a custom equality operator.

I initially ended up with the following:

type Point a = (a, a)
data Pair a = P (Point a) (Point a)

instance Eq (Pair a) where
	P a b == P c d = a == c && b == d || a == d && b == c

Which didn’t actually compile:

qc_test.hs:41:58:
    No instance for (Eq a)
      arising from a use of `=='
    In the second argument of `(&&)', namely `b == c'
    In the second argument of `(||)', namely `a == d && b == c'
    In the expression: a == c && b == d || a == d && b == c

The problem is that while we’ve made Pair an instance of the Equality type class there’s no guarantee that the value contained inside it is an instance of the Equality type class which means we might not be able to compare its values.

We need to add a class constraint to make sure that the value inside P is a part of Eq:

instance (Eq a) => Eq (Pair a) where
	P a b == P c d = a == c && b == d || a == d && b == c

Now we’re saying that we want to make Pair an instance of the Equality type class but only when the value that Pair contains is already an instance of the Equality type class.

In this case we’re just storing pairs of doubles inside the Pair so it will work fine.

Now if we compare the two points from above we’ll see that they’re equal:

> P (2.0, 0.0) (2.1, 1.1)  == P (2.1, 1.1) (2.0, 0.0)
True

I had to go and change the existing code to make use of this new but it didn’t take more than 5-10 minutes to do that.

This is the initial version that I wrote:

dcClosest :: (Ord a, Floating a) => [Point a] -> (Point a, Point a)
dcClosest pairs
  if length pairs <= 3 then = fromJust $ bfClosest pairs    
  else 
    foldl (\closest (p1:p2:_) -> if distance (p1, p2) < distance closest then (p1, p2) else closest) 
          closestPair 
          (windowed 2 pairsWithinMinimumDelta)
  where sortedByX = sortBy compare pairs	      
        (leftByX:rightByX:_) = chunk (length sortedByX `div` 2) sortedByX
        closestPair = if distance closestLeftPair < distance closestRightPair then closestLeftPair else closestRightPair  
          where closestLeftPair =  dcClosest leftByX
                closestRightPair = dcClosest rightByX
        pairsWithinMinimumDelta = sortBy (compare `on` snd) $ filter withinMinimumDelta sortedByX
          where withinMinimumDelta (x, _) = abs (xMidPoint - x) <= distance closestPair   
                  where (xMidPoint, _) = last leftByX

We can remove the first if statement which checks the length of the list and replace it with pattern matching code like so:

dcClosest :: (Ord a, Floating a) => [Point a] -> (Point a, Point a)
dcClosest pairs
  | length pairs <= 3 = fromJust $ bfClosest pairs    
  | otherwise = foldl (\closest (p1:p2:_) -> if distance (p1, p2) < distance closest then (p1, p2) else closest)
                      closestPair 
                      (windowed 2 pairsWithinMinimumDelta)
    ...

We can also get rid of the if statement inside the first argument passed to ‘foldl’ and replace it with a call to ‘minimumBy’:

dcClosest :: (Ord a, Floating a) => [Point a] -> (Point a, Point a)
dcClosest pairs
  | length pairs <= 3 = fromJust $ bfClosest pairs    
  | otherwise = foldl (\closest (p1:p2:_) -> minimumBy (compare `on` distance) [closest, (p1, p2)])
                      closestPair 
                      (windowed 2 pairsWithinMinimumDelta)
    ...

We can do the same to replace the if statement where we work out the closestPair which results in this final version of the code:

dcClosest :: (Ord a, Floating a) => [Point a] -> (Point a, Point a)
dcClosest pairs
  | length pairs <= 3 = fromJust $ bfClosest pairs    
  | otherwise = foldl (\closest (p1:p2:_) -> minimumBy (compare `on` distance) [closest, (p1, p2)])
                      closestPair 
                      (windowed 2 pairsWithinMinimumDelta)
  where sortedByX = sortBy compare pairs	      
        (leftByX:rightByX:_) = chunk (length sortedByX `div` 2) sortedByX    
        closestPair = minimumBy (compare `on` distance) [closestLeftPair, closestRightPair]
          where closestLeftPair =  dcClosest leftByX
                closestRightPair = dcClosest rightByX     
        pairsWithinMinimumDelta = sortBy (compare `on` snd) $ filter withinMinimumDelta sortedByX
          where withinMinimumDelta (x, _) = abs (xMidPoint - x) <= distance closestPair
                  where (xMidPoint, _) = last leftByX

It takes up marginally less space and I think the change to use pattern matching on the length of ‘pairs’ makes the biggest difference as the code is now lined up at the same level of indentation.

The other changes would have more of an impact if there were more than 2 things being compared – right now I think either of the versions of the code are equally readable.

As a result I had to update a traversal I’d written for finding the shortest path between two people in the graph.

The original traversal query I had was really simple because I had a direct connection between the people nodes:

neo = Neography::Rest.new
paths = neo.get_paths(start_node,
                      destination_node,
                      { "type" => "colleagues" },
                      depth = 3,
                      algorithm = "shortestPath")
           .map { |x| x["nodes"] }
           .uniq
paths.map { |p| p.map { |node| neo.get_node_properties(node, "name")["name"] } }

I changed the way the graph was modelled so that you needed to follow a ‘worked_on’ relationship to a project in order to go between people:

In the first version I’d written some pre processing code in Ruby to check whether or not people worked on the project at the same time before creating the relationship between the nodes.

It wasn’t possible to do that with the new structure since I was working out if there was a colleagues relationship dynamically.

I therefore added a ‘start_date’ and ‘end_date’ property to the ‘worked_on’ relationship between a person node and project node so that I’d be able to take it into account when traversing the graph.]

I initially thought it would be possible to do this using cypher and wrote the following query:

start_node_id = neo.send(:get_id, start_node)
destination_node_id = neo.send(:get_id, destination_node)
 
query =  " START a=node(#{start_node_id}), x=node(#{destination_node_id})" 
query << " MATCH p = allShortestPaths( a-[:worked_on*]-x )" 
query << " RETURN p, extract(person in nodes(p) : person.name)"
paths = neo.execute_query(query)

I wasn’t sure how to do the filtering on ‘start_date’ and ‘end_date’ and Andres pointed out that it’s not actually currently possible to take relationship properties into account when traversing a graph with cypher so we need to do the filtering on ‘start_date’ and ‘end_date’ in code.

My first attempt to do that looked like this:

paths = neo.get_paths(node1, node2, { "type" => "worked_on", "direction" => "all" }, 
                      depth = 5, algorithm = "shortestPath").uniq
matching = paths.select do |row|
  relationshipPairs = row["relationships"].each_slice(2).to_a
  relationshipPairs.all? do |pair|
    r1 = neo.get_relationship(pair[0])["data"]
    r2 = neo.get_relationship(pair[1])["data"]
    r1["start_date"] <= r2["end_date"] && r1["end_date"] >= r2["start_date"]
  end
end

The problem with this approach is that it’s really slow due to the fact that I’m pulling back every relationship back in order to check the start and end dates.

Michael Hunger suggested an alternative approach where I still used a cypher query but returned the relationships instead of just the nodes:

start_node_id = neo.send(:get_id, start_node)
destination_node_id = neo.send(:get_id, destination_node)
 
query =  " START a=node(#{start_node_id}), x=node(#{destination_node_id})" 
query << " MATCH p = allShortestPaths( a-[:worked_on*]-x )" 
query << " RETURN p, rels(p), extract(person in nodes(p) : person.name)"
paths = neo.execute_query(query)
matching = paths["data"].select do |row|
  relationship_pairs = row[1].each_slice(2).to_a
  relationship_pairs.all? do |pair|
    r1 = pair[0]["data"]
    r2 = pair[1]["data"]
    r1["start_date"] <= r2["end_date"] && r1["end_date"] >= r2["start_date"]
  end
end
 
matching.map { |x| x[2].each_with_index.select { |x,idx| idx.even? }.map(&:first) }.uniq

This approach mostly works although it runs into problems in the scenario where two people have worked on the same project but not at the same time.

In that scenario the above code will return the relationship between them but it will then be filtered out by the start date/end date logic which means we won’t see a shortest path between those nodes.

I’m not sure how to solve that problem so for the moment I’m going to take the Josh Adell’s suggestion to keep the ‘colleagues’ relationship between nodes and use that relationship for the shortest path traversal.

The ‘worked_on’ relationship is still useful for other things that I want to do but not this particular one.

For example say we want to use the read function to parse an input string into an integer.

We could do this:

> read "1" :: Int
1

But if we just evaluate the function alone and try and assign the result without casting to a type we get an exception:

> let x  = read "1"
 
<interactive>:1:10:
    Ambiguous type variable `a0' in the constraint:
      (Read a0) arising from a use of `read'
    Probable fix: add a type signature that fixes these type variable(s)
    In the expression: read "1"
    In an equation for `x': x = read "1"

sepp2k shows how we can provide a type declaration in GHCI in his Stack Over Flow answer:

> let x::Int; x = read "1"
> x
1

We can also use it when creating a list of integers to ensure they are of type ‘Int’ rather than ‘Integer’ for example:

> let y = [1,2,3]
> :t y
y :: [Integer]
 
> let y::[Int]; y = [1,2,3]
> :t y
y :: [Int]

It’s a pretty simple thing but I didn’t know it was even possible!

time ./closest_pairs 100 bf 
./closest_pairs 100 bf  0.01s user 0.00s system 87% cpu 0.016 total
 
time ./closest_pairs 1000 bf
./closest_pairs 1000 bf  3.59s user 0.01s system 99% cpu 3.597 total
 
time ./closest_pairs 5000 bf
./closest_pairs 5000  554.09s user 0.36s system 99% cpu 9:14.46 total

Luckily there’s a divide and conquer algorithm we can use which brings down the running time from O(n²

That algorithm is defined like so:

1. Sort points along the x-coordinate
2. Split the set of points into two equal-sized subsets by a vertical line x = x_mid
3. Solve the problem recursively in the left and right subsets. This will give the left-side and right-side minimal distances d_Lmin and d_Rmin respectively.
4. Find the minimal distance d_LRmin among the pair of points in which one point lies on the left of the dividing vertical and the second point lies to the right (a split pair).
5. The final answer is the minimum among d_Lmin, d_Rmin, and d_LRmin.

By step 4 in which we’ll look for a closest pair where the x and y values are in opposite subsets we don’t need to consider the whole list of points again since we already know that the closest pair of points is no further apart than dist = min(d_Lmin, d_Rmin).

We can therefore filter out a lot of the values by only keeping values which are within dist of the middle x value.

If a point is further away than that then we already know that the distance between it and any point on the other side will be greater than dist so there’s no point in considering it.

I used the C# version from Rosetta Code as a template and this is what I ended up with:

import Data.Maybe
import Data.List
import Data.List.Split
import Data.Function
 
dcClosest :: (Ord a, Floating a) => [Point a] -> (Point a, Point a)
dcClosest pairs
  if length pairs <= 3 then = fromJust $ bfClosest pairs    
  else 
    foldl (\closest (p1:p2:_) -> if distance (p1, p2) < distance closest then (p1, p2) else closest) 
          closestPair 
          (windowed 2 pairsWithinMinimumDelta)
  where sortedByX = sortBy compare pairs	      
        (leftByX:rightByX:_) = chunk (length sortedByX `div` 2) sortedByX
        closestPair = if distance closestLeftPair < distance closestRightPair then closestLeftPair else closestRightPair  
          where closestLeftPair =  dcClosest leftByX
                closestRightPair = dcClosest rightByX
        pairsWithinMinimumDelta = sortBy (compare `on` snd) $ filter withinMinimumDelta sortedByX
          where withinMinimumDelta (x, _) = abs (xMidPoint - x) <= distance closestPair   
                  where (xMidPoint, _) = last leftByX

If the number of pairs is 3 or less then we can just use the brute force algorithm since the whole splitting the list in half thing doesn’t work so well with a list of less than 4 items.

The main function is a fold which goes over all the pairs of points where we’ve determined that there might be a closest pair with one point on the left hand side of our ‘sorted by x list’ and the other on the right hand side.

We use the ‘windowed‘ function to pair up the points which in this case does something like this:

> windowed 2 [(0,0), (1,1), (2,2), (1.1, 1.2)]
[[(0.0,0.0),(1.0,1.0)],[(1.0,1.0),(2.0,2.0)],[(2.0,2.0),(1.1,1.2)]]

We pass along the closest pair that we’ve found from pairs of points on the left and right hand sides of the list as our seed value and check whether any of the split pairs are closer together than that one. By the time the fold is finished we will have the closest pair in the list.

The code in the where section is used to split the list into two halves, sort it by x and then work out which side of the list has the closest pair. That’s all done recursively and then the last bit of code works out which values we need to consider for the split pairs part of the algorithm.

I learnt about the ‘on‘ function while writing this algorithm which makes it really easy to pick a function to sort a collection with e.g. ‘compare `on` snd’ lets us sort a list in ascending order based on the second value in a tuple.

One problem with the way I’ve written this algorithm is that it places a lot of focus on the split pair bit of the code when actually a big part of why it works is that we’re sorting it into an order that makes it easier to work with and then dividing the problem by 2 each time.

My current thinking is that perhaps I should have had that code in the main body of the function rather than hiding it away in the where section.

It isn’t actually necessary to execute the foldl all the way through since we’ll often know there’s not going to be a split pair before we go through all the pairs. I thought it reads pretty nicely as it is thought and it runs pretty quickly as it is anyway!

Using the same data set as before (and getting the same answers!):

> time ./closest_pairs 1000 dc
./closest_pairs 1000 dc  0.02s user 0.00s system 91% cpu 0.024 total
 
> time ./closest_pairs 5000 dc
./closest_pairs 5000 dc  0.11s user 0.01s system 97% cpu 0.118 total

I described the code for generating random numbers in an earlier blog post but I’ll include it again to show how the algorithm is wired up:

import Control.Monad.State (State, evalState, get, put)
import System.Random (StdGen, mkStdGen, random)
import System
 
type R a = State StdGen a
rand :: R Double
rand = do
  gen <- get
  let (r, gen') = random gen
  put gen'
  return r
 
randPair :: R (Double, Double)
randPair = do
  x <- rand
  y <- rand
  return (x,y)
 
runRandom :: R a -> Int -> a
runRandom action seed = evalState action $ mkStdGen seed
 
normals :: R [(Double, Double)]
normals = mapM (\_ -> randPair) $ repeat ()
 
main = do 
	args <- getArgs
	let numberOfPairs = read (head args) :: Int
	if length args > 1 && args !! 1 == "bf" then 
		putStrLn $ show ( (bfClosest $ take numberOfPairs $ runRandom normals 42))
	else 
		putStrLn $ show ( (dcClosest $ take numberOfPairs $ runRandom normals 42))

The full code is on hpaste if anyone has any suggestions for how to improve it.

I’ve tried to work out how to create lists of random numbers in Haskell before and always ended up giving up because it seemed way more difficult than it should be but this time I came across a really good explanation of how to do it by jrockway on Stack Overflow.

I thought it’d be interesting to work up to his answer by first looking at the functions available to help us generate random numbers.

mkStdGen is a function which creates a random number generator for a given seed value and seems like a good place to start.

mkStdGen :: Int -> StdGen

The function mkStdGen provides an alternative way of producing an initial generator, by mapping an Int into a generator. Again, distinct arguments should be likely to produce distinct generators.

We can combine that with ‘randomR’ to give us a random value in a range that we choose or with ‘random’ to get a random value in a default range.

randomR :: RandomGen g => (a, a) -> g -> (a, g)

Takes a range (lo,hi) and a random number generator g, and returns a random value uniformly distributed in the closed interval [lo,hi], together with a new generator. It is unspecified what happens if lo>hi. For continuous types there is no requirement that the values lo and hi are ever produced, but they may be, depending on the implementation and the interval.

> fst $ randomR (1,10) (mkStdGen 66)
8

> let randomNumber :: Int; x = fst $ random (mkStdGen 66)
> randomNumber
-4755267682593970340

For now I’m ignoring the new random generator that is provided as the second part of the tuple that these two functions return.

If we want to create an infinite list of random numbers then something like this will do the trick:

let values :: [Int]; values = map fst $ scanl (\(r, gen) _ -> random gen) (random (mkStdGen 1)) $ repeat ()

We use the ‘repeat’ function to create an infinite sequence of nothing then run a scanl across that infinite sequence, called ‘random’ with a new instance of the random number generator each time:

> take 10 values
[7917908265643496962,-1017158127812413512,-1196564839808993555,128524678767915218,3573078269730797745,-2026762924287790163,-5402471397730582264,-8620480566299071809,5987841344909700899,1198540087579679591]
it :: [Int]

In this version we’re passing the next generator along as part of the accumulator of the ‘scanl’ function but another way to do that would be to encapsulate it inside a State monad:

Using the State monad our random number function would read like this:

1
2
3
4
5
6

myRand :: State StdGen Int
myRand = do
	gen <- get
	let (r, nextGen) = random gen
	put nextGen
	return r

The state monad contains some state (the random number generator) and creates values (the random numbers).

Line 3 gets the current generator.
Line 4 creates a new number and generator which it assigns to ‘r’ and ‘nextGen’
Line 5 updates the State monad with the new generator
Line 6 returns the value that was created.

We change ‘values’ to read like this:

let values = mapM (\_ -> myRand) $ repeat ()

‘mapM (\_ -> myRand)’ is the same as doing ‘sequence $ map (\_ -> myRand)’ and we need to do that because we want to return a State Monad with an array of Ints rather than an array of State Monads each with an Int.

In order to evaluate those values we call ‘evalState’ which returns the final result from the State monad i.e. the random value created:

> evalState values $ mkStdGen 1
[7917908265643496962,-1017158127812413512,-1196564839808993555,128524678767915218,3573078269730797745,-2026762924287790163,-5402471397730582264,-8620480566299071809,5987841344909700899,1198540087579679591,5081506216124746781,-2413363174299075397,-6804412472718217891,4559850124463334118,-362777938400029309,-23100237439495333,-2426472460098089322,...]

In our case we want to get an infinite sequence of pairs to represent the coordinates on a plane so we need to have a random function that creates those:

myRandPair :: State StdGen (Int, Int)
myRandPair = do
	x <- myRand
	y <- myRand
	return (x, y)

We then change ‘values’ accordingly:

> let values = mapM (\_ -> myRandPair) $ repeat ()
> evalState values $ mkStdGen 1
[(7917908265643496962,-1017158127812413512),(-1196564839808993555,128524678767915218),(3573078269730797745,-2026762924287790163),(-5402471397730582264,-8620480566299071809),(5987841344909700899,1198540087579679591),(5081506216124746781,-2413363174299075397),(-6804412472718217891,4559850124463334118),(-362777938400029309,-23100237439495333)...]

If we want to get back some different pairs then we just need to change the initial seed value that we pass to ‘mkStdGen’.

jrockway explains this in more detail in his post in which he also wraps the State monad inside an alias to make the code easier to understand.