A few days ago I came across this really interesting problem on the ‘dailyprogrammer‘. Here’s the question.

*Imagine you’ve got a square in 2D space, with axis values between 0 and 1, like this diagram. The challenge today is conceptually simple: can you place a circle within the square such that exactly half of the points in the square lie within the circle and half lie outside the circle, like here? You’re going to write a program which does this – but you also need to find the smallest circle which solves the challenge, ie. has the minimum area of any circle containing exactly half the points in the square.*

*This is a hard challenge so we have a few constraints:*

*Your circle must lie entirely within the square (the circle may touch the edge of the square, but no point within the circle may lie outside of the square).**Points on the edge of the circle count as being inside it.**There will always be an even number of points.*

*There are some inputs which cannot be solved. If there is no solution to this challenge then your solver must indicate this – for example, in this scenaro, there’s no “dividing sphere” which lies entirely within the square.*

I mulled this over for quite a while and the first ‘*light going on*‘ was that the minimum circle would have to have some of the points on its circumference – because if it didn’t then it wouldn’t be the minimum circle and so could be made smaller. Combining this idea with the fact that a circle can be drawn through any three points (or the degenerate case of two points being coincident and diametrically opposed to the other point) gave me the idea of creating all possible circles from three and then from two points and counting how many points are are in these circles.

So let’s start with some types…

1 2 3 4 |
-- x,y type Point = (Float, Float) -- center (x,y) and radius. type Circle = (Point, Float) |

and now a function to create a circle from three points. Notice that we haven’t implemented it. We’re just exploring a solution based on types and function signatures. Where an implementation is simple or obvious we’ll do it, otherwise it will be done later on a second iteration.

1 2 |
makeCircle :: Point -> Point -> Point -> Circle makeCircle = undefined |

As we can make one circle let’s use it to create all possible circles from a list of points using the *makeCircle* function.

1 2 |
allCircles :: [Point] -> [Circle] allCircles ps = [makeCircle p1 p2 p3 | p1 <- ps, p2 <- ps, p3 <- ps, p1 /= p2 && p1 /= p3 && p2 /= p3] |

All we’re doing is taking the points three at a time, making sure that they’re all different, then calling *makeCircle* with them. Even though *makeCircle* isn’t implemented this all compiles and type checks.

We also need the ‘two point circles’ – these can be obtained by calling *makeCircle* with the the second and third point being the same. So *allCircles* now looks like…

1 2 3 |
allCircles :: [Point] -> [Circle] allCircles ps = [makeCircle p1 p2 p3 | p1 <- ps, p2 <- ps, p3 <- ps, p1 /= p2 && p1 /= p3 && p2 /= p3] ++ [makeCircle p1 p2 p2 | p1 <- ps, p2 <- ps, p1 /= p2 && p2 /= p1] |

At some point we will need to test if a given point is in a given circle. This is simple enough and is based on the equation of a circle and the implementation is a couple of lines.

1 2 |
isInCircle :: Point -> Circle -> Bool isInCircle (x, y) ((cx, cy), r) = (x - cx)^2 + (y - cy)^2 - r^2 < 0.0001 |

where 0.0001 is fairly arbitrary constant that we could perhaps ‘tweak’.

We can now write a function that takes each point in turn and checks to see if it is in a given circle.

1 2 3 4 5 |
countPointsIn :: [Point] -> Circle -> Int countPointsIn [] _ = 0 countPointsIn (p:ps) c | isInCircle p c = 1 + countPointsIn ps c | otherwise = countPointsIn ps c |

Line 2 is the base case, if we have no points then a total of zero points are in any circle.

The general case checks for the point at the head of the list being in the circle and acts accordingly.

The next function is

1 2 |
isMinCircle :: [Point] -> Circle -> Bool isMinCircle ps c = countPointsIn ps c == min where min = div (length ps) 2 |

Part of condition for the circle to be the minimum circle is that it has exactly half of the points inside of it. Which is what *isMinCircle* succinctly expresses.

The next function is really the key to it all. The *getMinCirlce* first gets all the circles and, crucially, **sorts them on the ascending radius of the circle**. Then it takes each circle in turn until it finds one that has the desired number of points inside of it then stops. This will be the smallest circle as they are sorted by smallest radius first.

1 2 3 4 5 6 7 8 9 10 |
getMinCircle :: [Point] -> Maybe Circle getMinCircle ps = result where -- smallest radius first circs = sortBy (comparing snd) . allCircles $ ps -- keep dropping circle while they don't have correct number of points inside maybeList = dropWhile (\x -> not $ isMinCircle ps x ) circs -- pull result out result = f maybeList where f [] = Nothing f ls = Just (head ls) |

The *main* function is below along with a couple of helpers, one to show the result and the other to create the points from an input file.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
toPoint :: [String] -> Point toPoint (x:y:[]) = (read x :: Float, read y :: Float) showResult :: Maybe Circle -> String showResult c = case c of Nothing -> "Nothing" Just ((x, y), d ) -> show x ++ " " ++ show y ++ " " ++ show d main :: IO () main = do input <- readFile "points.txt" let points = map toPoint . map words . lines $ input print $ showResult . getMinCircle $ points |

and the input text file is like this…

1 2 3 4 5 |
0.17741 0.3643 0.70953 0.64191 0.8507 0.40874 0.56568 0.35134 ... etc. |

which is slightly different from that on the puzzle link in that it does not have an integer at the start giving the total number of points. (We get that from the length of the list of points).

Finally we need to revisit the function
makeCircle :: Point -> Point -> Point -> Circle which is really an exercise in coordinate geometry. Take the general equation of a circle, plug in the points and solve the resulting equations. The page show the details which I have implemented directly below.

1 2 3 4 5 6 |
makeCircle :: Point -> Point -> Point -> Circle makeCircle (x1, y1) (x2, y2) (x3, y3) = ((x, y), r) where k = 2 * (x1 * (y2 - y3) - y1 * (x2 - x3) + x2 * y3 - x3 * y2) x = ((x1^2 + y1^2) * (y2 - y3) + (x2^2 + y2^2) * (y3 - y1) + (x3^2 + y3^2) * (y1 - y2)) / k y = ((x1^2 + y1^2) * (x3 - x2) + (x2^2 + y2^2) * (x1 - x3) + (x3^2 + y3^2) * (x2 - x1)) / k r = sqrt ((x - x1)^2 + (y - y1)^2) |

Here is the complete code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 |
import Data.Ord import Data.List (sortBy) -- x,y type Point = (Float, Float) -- center (x,y) and radius. type Circle = (Point, Float) makeCircle :: Point -> Point -> Point -> Circle makeCircle (x1, y1) (x2, y2) (x3, y3) = ((x, y), r) where k = 2 * (x1 * (y2 - y3) - y1 * (x2 - x3) + x2 * y3 - x3 * y2) x = ((x1^2 + y1^2) * (y2 - y3) + (x2^2 + y2^2) * (y3 - y1) + (x3^2 + y3^2) * (y1 - y2)) / k y = ((x1^2 + y1^2) * (x3 - x2) + (x2^2 + y2^2) * (x1 - x3) + (x3^2 + y3^2) * (x2 - x1)) / k r = sqrt ((x - x1)^2 + (y - y1)^2) isInCircle :: Point -> Circle -> Bool isInCircle (x, y) ((cx, cy), r) = (x - cx)^2 + (y - cy)^2 - r^2 < 0.0001 allCircles :: [Point] -> [Circle] allCircles ps = [makeCircle p1 p2 p3 | p1 <- ps, p2 <- ps, p3 <- ps, p1 /= p2 && p1 /= p3 && p2 /= p3] ++ [makeCircle p1 p2 p2| p1 <- ps, p2 <- ps, p1 /= p2 && p2 /= p1] countPointsIn :: [Point] -> Circle -> Int countPointsIn [] _ = 0 countPointsIn (p:ps) c | isInCircle p c = 1 + countPointsIn ps c | otherwise = countPointsIn ps c isMinCircle :: [Point] -> Circle -> Bool isMinCircle ps c = countPointsIn ps c == min where min = div (length ps) 2 getMinCircle :: [Point] -> Maybe Circle getMinCircle ps = result where circs = sortBy (comparing snd) . allCircles $ ps maybeList = dropWhile (\x -> not $ isMinCircle ps x ) circs result = f maybeList where f [] = Nothing f ls = Just (head ls) toPoint :: [String] -> Point toPoint (x:y:[]) = (read x :: Float, read y :: Float) showResult :: Maybe Circle -> String showResult c = case c of Nothing -> "Nothing" Just ((x, y), d ) -> show x ++ " " ++ show y ++ " " ++ show d main :: IO () main = do input <- readFile "points.txt" let points = map toPoint . map words . lines $ input print $ showResult . getMinCircle $ points |

Test input of 50 points.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 |
0.17741 0.3643 0.70953 0.64191 0.8507 0.40874 0.56568 0.35134 0.56568 0.76805 0.96336 0.43929 0.91685 0.67614 0.61353 0.74929 0.32178 0.57213 0.90087 0.04679 0.47266 0.80022 0.35282 0.5678 0.12193 0.78317 0.27018 0.70962 0.60974 0.67577 0.00205 0.78244 0.30361 0.38056 0.88588 0.01516 0.85813 0.46199 0.16135 0.02435 0.09417 0.19733 0.96364 0.07113 0.97072 0.12663 0.70131 0.12721 0.63128 0.28139 0.08434 0.72938 0.73983 0.58714 0.08425 0.0076 0.50742 0.65794 0.9657 0.36221 0.98576 0.90219 0.94785 0.12917 0.10915 0.31807 0.61073 0.6223 0.17409 0.72275 0.53869 0.77369 0.10818 0.73822 0.57002 0.0511 0.18551 0.40974 0.95258 0.76844 0.77214 0.06548 0.93926 0.62406 0.8117 0.78065 0.61618 0.70183 0.51838 0.48385 0.20152 0.23971 0.68713 0.89812 0.65714 0.91779 0.95237 0.93368 0.45216 0.16395 |

and the output of running the code.

1 2 |
"0.65400547 0.5596951 0.37332553" [Finished in 3.1s] |

I checked these results using the link from the puzzle and it passed! 🙂

This is certainly not an optimised solution but what I like about it is that the problem can be ‘seen’ in the solution. It also shows how Haskell lets you build up a solution by writing small, very focussed functions, then composing them to solve a non-trivial problem. Overall I really enjoyed doing this and may well explore extending it to three or more dimensions – perhaps exploiting the potential abstraction that may be gained by taking a vector based approach. As usual the code is on Github.

Thanks for reading!