My first ipelet: 1-medians

I just wrote my first ipelet, and it was an amazingly pleasant experience. Not to mention the thrill I got from marking a collection of points and whoosh! getting the answer right there on the screen.

The ipelet is for computing the discrete and geometric 1-median of a set of points. The geometric 1-median of a set of points (also called the Fermat-Weber-Steiner-Torricelli-< insert your name here > point) is the point that minimizes the sum of distances to a set of points. The discrete 1-median is the same minimization, restricted to points in the input.

Computing the discrete 1-median of a set of points in the plane is trivial in n^2 time, and I didn't do anything clever there. Computing the geometric 1-median is more nontrivial though: there is no closed form solution beyond 4 points, and it can be shown that the point cannot be computed using arithmetic and k^th roots.

There are various approximation schemes: the most pertinent one here (in the plane) is a scheme by Bose, Maheshwari and Morin that uses a cone decomposition of the plane to yield an approximation scheme. In practice, researchers often use an iterative scheme developed by Weiszfeld; it's a simple iterative scheme built around the fact that the optimal solution must be a point such that the sum of unit vectors from it to all other points is zero (formally, this is the condition on the gradient).

The cost function is convex and differentiable except at the input points themselves. Thus, the Weiszfeld iteration is guaranteed to give the optimal solution as long as it doesn't get stuck at one of the data points. This algorithm is often referred to as the "founding algorithm of location science", and there are numerous modifications, extensions, and refinements. The ipelet implements the basic Weiszfeld scheme with a fairly brain-dead termination condition.

What possessed me to do such a thing ? Well, I'm trying to find counter examples for a conjecture involving 1-medians, and needed a way of visualizing examples :).

The code is here: feel free to use it, comment on it, modify it, etc.