Ever Decreasing Circles. 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.

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.

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

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…

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.

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.

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

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.

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.

and the input text file is like this…

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.

Here is the complete code:

Test input of 50 points.

and the output of running the code.

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.

• matrixmike says:

Excellent work – I am using this as a self-propelled learning exercise. Just gingerly stepping into Haskell and revisiting Maths etc – I started documenting the (your) original code on my fork then discovered this piece of beautiful work. Turns out I have been tinkering with your code for a while…
I could add my website later when I feel more confident …

• BanditPig says:

Thanks again! Haskell does become slightly addictive when exploring interesting maths problems like this.
What are you doing on your website?

• matrixmike says:

In one of your posts, I have started to read a few, but want to save some, as a treat when I need a boost, you ask about code improvements. Did you ever consider using hlint? I use it when I explore the code from other people because it produces a standard format or should I say style. I can sometimes see, before I use it, that the code was written by someone used to other styles and/or languages. If the original format is useful I add the old way as a comment.

• BanditPig says:

Thanks for the compliment. Yes, I have started using hlint. I recently added it as a plugin to my VScode setup. Sometimes the hlint suggestion seems (to me) overly complex but this isn’t too often and usually the change is an improvement.