The Smallest Circle (#157)

Bill K. wrote:

in fact, pass the same point data to each solution per run.
I had hoped by averaging multiple runs (10) - each run having
different point data - that the rankings would have been pretty
stable.

However, JUSTIN and LIONEL still tend to trade places on
different invocations:

Here’s a 157_benchmark_2.rb which benchs on several point distributions
(square, disk, circle and gaussian) and my code adapted to the benchmark
(the convex hull was done outside of encircle so it wasn’t active for
the benchmark).

Lionel

Array#combination is ruby19 method. But it’s probably a good thing you
didn’t get my “solution” to run. If I remove irrelevant points in
advance, I get somewhat close to Frank’s solution with respect to
speed but it still lags behind.

Anyway, I can offer some accuracy tests with random point sets (size
1…100).

 s_tml2.rb | s_billk.rb | s_billk2.rb | s_lionel.rb | s_frank.rb |

s_justin.rb
Min: 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0004 |
0.0000
Max: 0.0225 | 0.0795 | 0.0422* | 0.0140 | 0.0063* |
0.0817
Avg: 0.0027 | 0.0297 | ------- | 0.0017 | 0.0058 |
0.0115

[*] Slightly invalid results, maybe the center isn’t quite where it
should be?

x/y-range is -0.5…0.5

Since the number of test runs is rather small, you should take into
account a margin of 30% or so.

s_lionel is Lionel’s very first submission because I have difficulties
to integrate the later solutions.

s_tml2 is my second solution that starts at the center of the longest
cord. The first solution yields results mostly identical to Justin’s
submission.

Regards,
Thomas.

From: “ThoML” [email protected]

should be?
I don’t quite understand how to interpret the numbers?

But I have to ask… which is billk vs. billk2 ? If billk2 is less
accurate, would I be correct in assuming it is the naive bounding
box version?

Thanks,

Bill

Bill K. wrote:

From: “Lionel B.” [email protected]

Here’s a 157_benchmark_2.rb which benchs on several point
distributions (square, disk, circle and gaussian) and my code adapted
to the benchmark (the convex hull was done outside of encircle so it
wasn’t active for the benchmark).

Neat! :slight_smile:

Oups : the 2D Gaussian was the same as the circle one (stupid me)… The
gradient walking couldn’t possibly be so slow on such a distribution :slight_smile:

I updated my submission with :

  • a new heuristic for the initial step,
  • an optimized step ratio for the gradient algorithm,
  • a bug fix (it didn’t run out of the benchmark).

Lionel

From: “Lionel B.” [email protected]

Here’s a 157_benchmark_2.rb which benchs on several point distributions
(square, disk, circle and gaussian) and my code adapted to the benchmark
(the convex hull was done outside of encircle so it wasn’t active for
the benchmark).

Neat! :slight_smile:

– Random on disk –
user system total real
FRANK 23.407000 0.313000 23.720000 ( 23.734000)
JUSTIN 4.359000 0.000000 4.359000 ( 4.360000)
DOUG 3.719000 0.031000 3.750000 ( 3.750000)
PHILIPP 2.703000 0.000000 2.703000 ( 2.704000)
LIONEL 0.734000 0.031000 0.765000 ( 0.765000)
BILL 0.282000 0.000000 0.282000 ( 0.281000)
– Random on square –
user system total real
FRANK 27.734000 0.282000 28.016000 ( 28.016000)
JUSTIN 7.422000 0.000000 7.422000 ( 7.422000)
DOUG 3.703000 0.000000 3.703000 ( 3.703000)
PHILIPP 1.422000 0.000000 1.422000 ( 1.422000)
LIONEL 0.750000 0.000000 0.750000 ( 0.750000)
BILL 0.344000 0.000000 0.344000 ( 0.343000)
– 2D Gaussian –
user system total real
FRANK 27.297000 0.093000 27.390000 ( 27.406000)
LIONEL 16.141000 0.000000 16.141000 ( 16.219000)
JUSTIN 8.047000 0.000000 8.047000 ( 8.094000)
PHILIPP 0.422000 0.000000 0.422000 ( 0.422000)
BILL 0.172000 0.000000 0.172000 ( 0.172000)
DOUG 0.140000 0.000000 0.140000 ( 0.140000)
– Random on circle –
user system total real
FRANK 28.187000 0.079000 28.266000 ( 28.344000)
LIONEL 16.547000 0.000000 16.547000 ( 16.657000)
JUSTIN 7.297000 0.000000 7.297000 ( 7.328000)
PHILIPP 0.468000 0.000000 0.468000 ( 0.469000)
DOUG 0.172000 0.000000 0.172000 ( 0.171000)
BILL 0.141000 0.000000 0.141000 ( 0.141000)

Regards,

Bill

Lionel B. wrote:

Oups : the 2D Gaussian was the same as the circle one (stupid me)…
The gradient walking couldn’t possibly be so slow on such a
distribution :slight_smile:

I updated my submission with :

  • a new heuristic for the initial step,
  • an optimized step ratio for the gradient algorithm,
  • a bug fix (it didn’t run out of the benchmark).

Lionel

Note that if you filter the points with the convex hull, Justin’s
algorithm becomes competitive with the gradient walking with roughly the
same performance profile:

– Random on disk –
user system total real
FRANK 36.730000 0.040000 36.770000 ( 37.330030)
DOUG 6.570000 0.010000 6.580000 ( 6.680385)
PHILIPP 5.290000 0.000000 5.290000 ( 5.377509)
JUSTIN 1.090000 0.010000 1.100000 ( 1.107219)
LIONEL 0.900000 0.010000 0.910000 ( 0.915325)
BILL 0.610000 0.000000 0.610000 ( 0.622287)
– Random on square –
user system total real
FRANK 35.420000 0.030000 35.450000 ( 36.056951)
DOUG 6.090000 0.010000 6.100000 ( 6.147446)
PHILIPP 2.720000 0.000000 2.720000 ( 2.766508)
JUSTIN 1.060000 0.000000 1.060000 ( 1.081301)
LIONEL 1.810000 0.010000 1.820000 ( 1.857297)
BILL 0.540000 0.000000 0.540000 ( 0.536318)
– 2D Gaussian –
user system total real
FRANK 36.370000 0.050000 36.420000 ( 37.079492)
DOUG 4.730000 0.000000 4.730000 ( 4.816418)
PHILIPP 3.000000 0.010000 3.010000 ( 3.036339)
JUSTIN 1.270000 0.000000 1.270000 ( 1.279535)
LIONEL 0.810000 0.010000 0.820000 ( 0.828373)
BILL 0.540000 0.000000 0.540000 ( 0.544507)
– Random on circle –
user system total real
FRANK 36.540000 0.030000 36.570000 ( 37.410887)
DOUG 0.260000 0.000000 0.260000 ( 0.255157)
PHILIPP 0.390000 0.000000 0.390000 ( 0.395307)
JUSTIN 11.350000 0.010000 11.360000 ( 11.568574)
LIONEL 12.730000 0.010000 12.740000 ( 13.040802)
BILL 0.170000 0.000000 0.170000 ( 0.166975)

Lionel

ThoML wrote:

I didn’t include Douglas’s and Philip’s solutions since they are
precise.

I’m not sure that you can draw a clear line between Douglas, Philip,
Justin and my solution versus precision. I’ve not studied the results of
each precisely and I’ve only read through the algorithms quickly, but my
guess is that:

  • Douglas and Philip are theoretically correct and Justin and mine are
    theoretically incorrect.
  • that said you must keep in mind that floating point computations are
    inherently incorrect most of the time.

So solutions relying on defining circles by 3 points (which uses line
intersections with a good deal of add and multiply on floating point
values) might and most probably often do give circles that don’t include
one of the points it is supposed to pass through.

In fact Justin and me have solutions that often give suboptimal
solutions but:

  • within a predictable range of the ideal solution (a small multiple of
    the floating point error).
  • as they rely on the max_distance from a point to the point cloud the
    circle is guaranteed to include all points.

Lionel

I think the benchmarking code is getting more interesting than the quiz
solutions :slight_smile:

-Doug

I don’t quite understand how to interpret the numbers?

But I have to ask… which is billk vs. billk2 ? If billk2 is less
accurate, would I be correct in assuming it is the naive bounding
box version?

billk2 is the revised version. The maximum distance is lower, thus
it’s more accurate. Basically, you’d want a 0.0 in all cells.

min … minimal distance to the winner (smallest circle)
max … maximal distance to the winner
avg … average (arithmetic mean) distance to the winner

The table doesn’t contain info about the diameter. I suppose a
relative measure would be more useful. Hm.

Here are two runs with relative numbers. The relevant measure is the
circle radius. The numbers are the distance in per cent of the correct
solution’s radius.

 s_tml.rb | s_billk.rb | s_billk2.rb | s_lionel.rb | s_frank.rb |

s_justin.rb
Min: 0.00 | 0.00 | 0.00 | 0.00 | 0.15 |
0.00
Max: 15.77 | 17.72 | 7.81 | 4.98 | 2.13 |
9.28
Avg: 0.54 | 4.79 | 1.18 | 0.43 | 1.03 |
1.38

  s_tml.rb | s_billk.rb | s_billk2.rb | s_lionel.rb | s_frank.rb |

s_justin.rb
Min: 0.00 | 0.00 | 0.00 | 0.00 | 0.11 |
0.00
Max: 5.37 | 15.99 | 7.39 | 5.80 | 2.01 |
12.10
Avg: 0.36 | 5.07 | 1.07 | 0.35 | 1.01 |
1.38

Okay, now with Lionel’s current version:

 s_tml.rb | s_billk.rb | s_billk2.rb | s_lionel.rb | s_lionel4.rb

| s_frank.rb | s_justin.rb
Min: 0.00 | 0.00 | 0.00 | 0.00 | 0.00
| 0.26 | 0.00
Max: 4.79 | 18.35 | 9.15 | 3.33 | 4.66
| 1.95 | 17.57
Avg: 0.54 | 5.44 | 1.20 | 0.22 | 0.36
| 0.99 | 1.76

As you can see, may solution fails noticably on certain point sets but
does more or less well in average. At least for random point sets. I
think Justin’s solution suffers from the same problem since his
approach is quite similar, only that he moves the center gradually bit
by bit while my “solution” … well, let’s not talk about that.

The performance of Frank’s solution is quite constant albeit slightly
imprecise.

I didn’t include Douglas’s and Philip’s solutions since they are
precise.

Regards,
Thomas.

From: “Douglas A. Seifert” [email protected]

I think the benchmarking code is getting more interesting than the quiz
solutions :slight_smile:

I dunno if there’s a better way to do it - but I’ll admit I had a lot of
fun
writing this line:

namespaces.each {|mod| solutions[mod.name] = Class.new { include mod
}.new}

:slight_smile:

Ruby is just so amazingly dynamic…

Regards,

Bill

On 2/20/08, Rick DeNatale [email protected] wrote:

First point
that it must lie on the circumference of the new circle which will
contain all of the points examined thus far. We now find the set of
points which are furthest from the new point (again within some
epsilon). There will of course be one or more of these points.

Now the center of the new circle is the average of those furthest
points plus the new points, and the radius is the distance between the
center and the new point.

Oh, well, it was a good try I guess, but I ran afoul of the case
exposed by Alex Shugin. The points farthest from the new point don’t
necessarily (all) have to be on the circumference of the circle.


Rick DeNatale

My blog on Ruby
http://talklikeaduck.denhaven2.com/

Okay, I know that this probably isn’t the best first impression I
could make, but is there someone who would mind doing a summary for
this week’s quiz?

I developed a fever of about 101F last night and have it today, will
probably carry it into tomorrow. I can still do a summary, though it
will likely show up Friday and be shorter than perhaps is usual,
because I can barely stay focused on something for more than 5 minutes
at the moment.

If someone is up for it, just reply here and have at it. I will
appreciate it enormously.

Here’s my solution.

I also eschewed doing any google research and started from scratch. I
think I’ve taken a slightly different approach than any already
described.

The basic approach is to start with a circle containing only the first
point, and then add the points one by one changing the circle as
necessary. So

First point
Set center to the point and radius to zero.

Subsequent points.
If the point is already in the circle then the point is simply
added to the collection of points and the circle is unchanged. The
point is in the circle if the distance from the center of the circle
is less than the radius (within some epsilon to deal with floating
point comparisons).

If the point is further from the center than the radius, the we know
that it must lie on the circumference of the new circle which will
contain all of the points examined thus far. We now find the set of
points which are furthest from the new point (again within some
epsilon). There will of course be one or more of these points.

Now the center of the new circle is the average of those furthest
points plus the new points, and the radius is the distance between the
center and the new point.

The code is in several files:

==== lib/circle.rb =====
class Circle < Struct.new(:center, :radius)
def to_s
“{center:#{center}, radius:#{radius}}”
end
end

==== lib/point.rb =====
class Point < Struct.new(:x, :y)

def self.random
new(rand, rand)
end

def to_s
“Point(#{x}, #{y})”
end

def inspect
to_s
end

def -(aPoint)
Point.new(x - aPoint.x, y - aPoint.y)
end

def dist_squared(other_point)
diff = self - other_point
diff.x * diff.x + diff.y * diff.y
end

def dist(other_point)
Math.sqrt(dist_squared(other_point))
end

def ==(other_point)
x = other_point.x && y == other_point.y
end

def center(other_point)
Point.new((x + other_point.x)/2, (y + other_point.y)/2)
end

end

==== lib/point_collection.rb ====
require File.dirname(FILE) + ‘/point’
require File.dirname(FILE) + ‘/circle’

class PointCollection
attr_reader :radius, :center

def initialize
@points = []
@radius = 0
end

def circle
Circle.new(center, radius)
end

def add_point(a_point)
@points << a_point
recalc_after_adding(a_point)
end

def empty?
@points.empty?
end

private

Recalculate the center and radius if a point has been added.

If the new point is within the current circle then no change

is needed.

If it isn’t then, it must be on the circumference of the new circle

In this case find the point or points furthest from the new point,

The new center will be the center of the new point plus these far

points

The new radius will be the distance of the new point from the new

center.
def recalc_after_adding(a_point)
if @center
if a_point.dist_squared(center) > @radius * @radius
circ_points = (farthest_from(a_point) + [a_point]).uniq
@center = center_of(circ_points)
@radius = @center.dist(a_point)
end
else
@center = a_point
end
end

Given an array of points find the center.

def center_of(an_array_of_points)
x_sum = y_sum = 0.0
an_array_of_points.each do |point|
x_sum += point.x
y_sum += point.y
end
Point.new(x_sum / an_array_of_points.length, y_sum /
an_array_of_points.length)
end

def farthest_from(a_point)
max_dist_squared = 0.0
delta = 0.000001
farthest = []
p @points
@points.each do |point|
dist_squared = point.dist_squared(a_point)
diff = dist_squared - max_dist_squared
if diff.abs <= delta
max_dist_squared = dist_squared if diff > 0
farthest << point
elsif diff > delta
max_dist_squared = dist_squared
farthest = [ point ]
end
end
farthest
end
end

==== spec/point_spec.rb
require File.dirname(FILE) + ‘/…/lib/point’

def reflect(point)
Point.new(-point.x, -point.y)
end

describe Point, “(4, 5)” do

before(:each) do
@point = Point.new(4, 5)
end

it “should have a dist_squared of 25.0 from point(1, 1)” do
@point.dist_squared(Point.new(1,1)).should == 25.0
end

it “should have a dist of 5 from point(1,1)” do
@point.dist(Point.new(1,1)).should be_close(5.0, 0.0000001)

end

it “should equal a new point with the same coordinates” do
@point.should == Point.new(4, 5)
end

it “should be its own center when paired with self” do
@point.center(@point).should == @point
end

it “should have a center of zero when paired with its reflection” do
@point.center(reflect(@point)).should == Point.new(0,0)
end
end

==== spec/point_collection_spec.rb =====
require File.dirname(FILE) + ‘/…/lib/point’
require File.dirname(FILE) + ‘/…/lib/point_collection’

describe PointCollection do
before(:each) do
@pc = PointCollection.new
end

it “should be empty initially” do
@pc.should be_empty
end

it “should should have a radius of zero” do
@pc.radius.should == 0
end

it “should not have a center” do
@pc.center.should == nil
end

describe PointCollection, “with one point” do
before(:each) do
@first_point = Point.new(1,1)
@pc.add_point(@first_point)
end

it "should not be empty" do
  @pc.should_not be_empty
end

it "should have the first point as the center" do
  @pc.center.should == @first_point
end

end

describe PointCollection, “with two points” do
before(:each) do
@first_point = Point.new(1,1)
@pc.add_point(@first_point)
@second_point = Point.new(5, 5)
@pc.add_point(@second_point)
end

it "should have the center of the two points as its center" do
  @pc.center.should == @first_point.center(@second_point)
end

it "should have a radius half the distance between the two points" 

do
@pc.radius.should be_close(@first_point.dist(@second_point) / 2.0,
0.001)
end

describe PointCollection, "adding a point inside" do

  before(:each) do
    @old_center, @old_radius = @pc.center, @pc.radius
    @pc.add_point(Point.new(4,4))
  end

  it "should have the same radius as before" do
    @pc.radius.should == @old_radius
  end

  it "should have the same center as before" do
    @pc.center.should == @old_center
  end
end

describe PointCollection, "adding a point ouside" do
  before(:each) do
    @old_center, @old_radius = @pc.center, @pc.radius
    @third_point = Point.new(6,6)
    @pc.add_point(@third_point)
  end

  it "should calculate the radius correctly" do
    @pc.radius.should

be_close(@first_point.dist(@third_point)/2.0, 0.0001)
end
end
end
end

This much was done about 5 days ago, I’ve spent time sporadically
since then getting rmagick working with X11 on OS X Leopard so that I
could visually convince myself that this was working.

Here’s the final piece, with the visualization code adapted from
someone elses submission (sorry but I can’t seem to find it).

==== lib/main.rb =====
#! /usr/bin/env ruby

require ‘rubygems’
require File.dirname(FILE) + ‘/point’
require File.dirname(FILE) + ‘/circle’
require File.dirname(FILE) + ‘/point_collection’

def generate_samples(n)
(1…n).map { Point.random }
end

def encircle(points)
point_collection = PointCollection.new
points.each do |point|
point_collection.add_point(point)
end
point_collection.circle
end

def scale(value)
$size / 6 + value * $size / 1.5
end

def draw_circle(circle, color=false)
center = circle.center
radius = circle.radius
$gc.stroke_color(color.to_s) if color
$gc.circle(scale(center.x), scale(center.y), scale(center.x),
scale(center.y+radius))
end

def draw_point(center, color=false)
$gc.stroke_color(color.to_s) if color
$gc.circle(scale(center.x), scale(center.y), scale(center.x),
scale(center.y)+2)
end

def draw_line(p0, p1, color=false)
$gc.stroke_color(color.to_s) if color
$gc.line(scale(p0.x), scale(p0.y), scale(p1.x), scale(p1.y))
end

require ‘RMagick’
$size = 600
$gc = Magick::Draw.new
$gc.fill_opacity(0)
$gc.stroke_width(2)
$gc.stroke_opacity(1)

samples = generate_samples(ARGV[0].to_i)
circle = encircle(samples)
puts circle

samples.each { |p| draw_point(p, :black) }
draw_circle(circle)
canvas = Magick::ImageList.new
canvas.new_image($size, $size)
$gc.draw(canvas)
canvas.display


Rick DeNatale

My blog on Ruby
http://talklikeaduck.denhaven2.com/

On Wed, Feb 20, 2008 at 2:49 PM, Matthew M. [email protected]
wrote:

If someone is up for it, just reply here and have at it. I will
appreciate it enormously.

Hi Matthew,

I think it might be better if you wait until you are better and do it.
Just
delay the next quiz. Or have the person who submitted the quiz (if it
wasn’t you) write the summary. I don’t think a free-for-all on the
summary
is a good idea.

Get well soon!

Eric

On Feb 20, 2008, at 2:49 PM, Matthew M. wrote:

If someone is up for it, just reply here and have at it. I will
appreciate it enormously.

Matthew, get well soon, I hope you didn’t catch the flu that i’m just
now recovering from!
I suggest just delaying the summary. It’s not the end of the world.
JJ

John J. wrote:

Matthew, get well soon, I hope you didn’t catch the flu that i’m just
now recovering from!
I suggest just delaying the summary. It’s not the end of the world.
JJ

Same here, get well soon and don’t stress yourself over the summary. I
think we can wait: the summary is more for people who didn’t follow the
quiz closely or who will find it searching the web.

Lionel

On Wed, Feb 20, 2008 at 3:05 PM, Eric M. [email protected]
wrote:

On Wed, Feb 20, 2008 at 2:49 PM, Matthew M. [email protected]
wrote:

Okay, I know that this probably isn’t the best first impression I
could make, but is there someone who would mind doing a summary for
this week’s quiz?

Hi Matthew,

Another thing… To ease the burden on you for writing the summaries,
you
might have the quiz submitter write the summary by default. Since that
person wrote the quiz, hopefully they should be a good person to write
the
summary. You could be a backup if they are unable for some reason. To
make
this work well, you’d want quiz ideas to come mainly from others instead
of
yourself.

Thanks for taking on the new role of quizmaster!

Eric

On Feb 20, 7:19 pm, Lionel B. [email protected]
wrote:

John J. wrote:

Matthew, get well soon, I hope you didn’t catch the flu that i’m just
now recovering from!
I suggest just delaying the summary. It’s not the end of the world.
JJ

Same here, get well soon and don’t stress yourself over the summary. I
think we can wait: the summary is more for people who didn’t follow the
quiz closely or who will find it searching the web.

Thanks for the support… I decided, after much lack of cowbell, that
I will delay the summary one week.

Hi,

Since Matthew has a cold, I decided to take another look at my
“solution”. Accuracy hasn’t changed because the idea is exactly the
same as before, but speed is at least measurable now. Well, overall
performance is still abominable but since I started it (i.e. posted
the previous version) I thought I could also post this slightly more
complex version. I also included some hopefully correct code for ruby
1.8 compatibility. The use of complex could be easily avoided of
course.

Maybe of more interest to you: I also include a revised version of my
accuracy checker that makes use of Lionel’s more sophisticated sample
generator. My “solution” performs especially bad on “random on disk”
and “2d gaussian” it seems. Oh well.

Regards,
Thomas.

========================================================================

require ‘complex’

class Point < Struct.new(:x, :y)
def self.random
Point.new(rand, rand)
end

def to_s
    "(#{x}, #{y})"
end

def self.from_complex(c)
    self.new(c.real, c.image)
end

end

class Circle < Struct.new(:center, :radius)
def to_s
“{center:#{center}, radius:#{radius}}”
end
end

def encircle(points) # takes array of Point objects
return if points.nil? or points.empty?
return Circle.new(points[0], 0) if points.size == 1

points = prepare(points)
a, ar  = points_max_dist_center(points)
return Circle.new(Point.from_complex(a), ar) unless points.size >

2

points = points.sort_by {|p| -(a - p).abs}
f      = points[0]
points.delete(f)
zz = nil
zr = 0
points.each do |p|
    aang = (a - f).angle
    pang = (p - f).angle
    phi  = pang - aang
    rd   = (f - p).abs / 2.0
    rr   = rd / Math.cos(phi)
    if rr > zr
        ra = f.real  + rr * Math.cos(aang)
        rb = f.image + rr * Math.sin(aang)
        zz = Complex(ra, rb)
        zr = rr
    end
end
return Circle.new(Point.from_complex(zz), zr)

end

def points_max_dist_center(points)
m, n = points.combination(2).sort_by {|a, b| -(a - b).abs}[0]
a = (m + n) / 2.0
d = (a - m).abs
return [a, d]
end

def prepare(points)
points = points.uniq
points.map! {|p| Complex(*p.to_a)}
i = 3
while i < points.size
p1, p2, p3 = points[i-3, i]
points.delete_if {|p| in_triangle?(p1, p2, p3, p)}
i += 1
end
return points
end

def in_triangle?(a, b, c, p)
return false if a.nil? or b.nil? or c.nil? or p.nil? or a == p or
b == p or c == p
u = (a - b).angle
l = (a - c).angle
x = (a - p).angle
return false unless x.between?(l, u)
u = (c - a).angle
l = (c - b).angle
x = (c - p).angle
return x.between?(l, u)
end

if RUBY_VERSION < ‘1.9.0’
class ::Array

    def combination(n)
        out = []
        size.times do |i|
            head  = self[i]
            if n > 1
                rest = self[i + 1, size]
                if rest.size >= n - 1
                    tails = rest.combination(n - 1)
                    out += tails.map {|t| t.unshift(head)}
                end
            else
                out << [head]
            end
        end
        out
    end

end

end

if FILE == $0
points = []
ARGV.each {|p| points << Point.new(*p.split(/,/).map{|c| c.to_f})}
c = encircle(points)
p c
end

========================================================================

#!/usr/bin/env ruby19

require ‘mathn’

E = 1e-10

num_runs = 10
num_samples = 100

num_runs = 10

num_samples = 1000

module PHILIPP ; eval(File.read(‘s_philip.rb’)) ; end
module PHILIPP2; eval(File.read(‘s_philip2.rb’)) ; end
module LIONEL ; eval(File.read(‘s_lionel.rb’)) ; end
module LIONEL4; eval(File.read(‘s_lionel4.rb’)) ; end
module DOUG ; eval(File.read(‘s_douglas.rb’)) ; end
module FRANK ; eval(File.read(‘s_frank.rb’)) ; end
module JUSTIN ; eval(File.read(‘s_justin.rb’)) ; end
module BILL ; eval(File.read(‘s_billk.rb’)) ; end
module BILL2 ; eval(File.read(‘s_billk2.rb’)) ; end
module THOMAS ; eval(File.read(‘quiz157l.rb’)) ; end
module THOMAS2 ; eval(File.read(‘quiz157s.rb’)) ; end

namespaces = [
THOMAS,
THOMAS2,
FRANK,
JUSTIN,
LIONEL,
LIONEL4,
DOUG,
PHILIPP,
PHILIPP2,
BILL,
BILL2,
]

class Point < Struct.new(:x, :y)
def self.random(type=0)
case type
when 1
ro = rand / 2
theta = Math::PI * 2 * rand
Point.new(ro * Math::cos(theta) + 0.5, ro *
Math::sin(theta) + 0.5)
else
Point.new(rand, rand)
end
end

def to_s
    "(#{x}, #{y})"
end

end

def generate_samples(mod, coords)
coords.map {|xy| mod::Point.new(*xy)}
end

def distance(p1, p2)
Math.sqrt((p1.x - p2.x) ** 2 + (p1.y - p2.y) ** 2)
end

distributions = {
“Random on square” =>
(1…num_runs).map {(1…num_samples).map {[rand, rand]}},
“Random on disk” =>
(1…num_runs).map {(1…num_samples).map {
ro = rand(0.5)
theta = Math::PI * 2 * rand
[ ro * Math::cos(theta), ro * Math::sin(theta)]
}
},
“Random on circle” =>
(1…num_runs).map {(1…num_samples).map {
theta = Math::PI * 2 * rand
[ Math::cos(theta), Math::sin(theta)]
}
},
“2D Gaussian” =>
(1…num_runs).map {(1…num_samples).map {
ro = Math::sqrt(-2 * Math::log(rand))
theta = Math::PI * 2 * rand
[ ro * Math::cos(theta), ro * Math::sin(theta)]
}
},
}

def prepare_pointsets(namespaces, raw_coords, num_runs)
pointsets = {}
namespaces.each do |mod|
pointsets[mod.name] = (0…num_runs).map do |i|
generate_samples(mod, raw_coords[i])
end
end
pointsets
end

puts
distributions.each { |name,raw_coords|
results = Hash.new {|h, k| h[k] = {}}
winners = Hash.new {|h,k| h[k] = 0}

puts "-- #{name} --"
pointsets = prepare_pointsets(namespaces, raw_coords, num_runs)
solutions = {}
namespaces.each {|mod| solutions[mod.name] = Class.new { include

mod }.new}
namespaces.each do |mod|
name = mod.name
solution = solutions[name]
num_runs.times do |i|
points = pointsets[name][i]
val = solution.encircle(points)
results[i][mod] = val
end
end

f_dist = Hash.new {|h,k| h[k] = []}

results.each do |i, s|
    cwinner = namespaces.sort_by {|f| s[f] ? s[f].radius :

100000000}
best_radius = s[cwinner[0]].radius
winners[cwinner[0]] += 1
cwinner.each_cons(2) do |a,b|
if (s[a].radius - s[b].radius).abs <= E
winners[b] += 1
else
break
end
end

    namespaces.each do |fav|
        f_d = ((s[fav].radius - best_radius) / best_radius).abs *

100
f_dist[fav] << f_d
end
end

puts "Winners (best circle):"
winners.each do |f, n|
    puts '%10s %d' % [f, n]
end
puts

puts "Fav: %s" % namespaces.join(', ')
puts 'Min: %s' % f_dist.map {|fav, d| '%6.2f' % d.min}.join(', ')
puts 'Max: %s' % f_dist.map {|fav, d| '%6.2f' % d.max}.join(', ')
puts 'Avg: %s' % f_dist.map {|fav, d| '%6.2f' % (d.inject(0.0){|

a,d| a + d} / d.size)}.join(’, ')
puts

}

First, I want to thank everyone participating and watching this week
for being patient with me while I was down with the flu. I am almost
completely recovered at this point, and I hope that doesn’t come back
for a good, long time.

Second, thanks to all quizzers who made an attempt at the quiz, even
after the scary warning I made about the difficulty.

Third, I had hoped to provide a set of tests to help quizzers check
their solutions, but the flu hit me pretty hard and I was unable. A
special thanks go to Bill K., Thomas ML and Lionel B.; they
took time to benchmark the solutions for both speed and accuracy, when
I didn’t have the ability to do it.

Now, about the smallest enclosing circle problem… I warned from the
start that the problem might prove more difficult than it appears.
It’s a problem that has been attacked and solved numerous times. [1][2]
[3][4][5]

One easy to understand algorithm was originally described by Welzl[6],
and Rick DeNatale almost came up with the same answer:

If the point is further from the center than the radius, the we know that
it must lie on the circumference of the new circle which will contain all
of the points examined thus far.

The point where Rick’s propsed algorithm fails is in determining what
other points will be part of the circle’s boundary.

I don’t intend to get into the finer details of these algorithms here;
I’ve provided a number of references for those who are interested in
crafting an exact solution. In particular, Douglas S. did provide
such a solution.

The solution I am going to examine here is that of Frank F.,
whose solution is actually the slowest according to the benchmarks.
However, his solution was consistently close, simple to understand,
and provides a technique that can be used for many kinds of problems,
perhaps not always for speed but sometimes for sanity.

Frank’s Circle class is identical to that provided by the quiz, so
I’ll skip that. Let’s look at his modified Point class:

class Point
    def initialize( *coords )
        @coords = coords.map{|x| x.to_f}
    end

The initializer takes an array of numbers, rather than just the two
that my original Point class (derived from Struct) provided. Frank’s
goal is to support more than just two dimensions; in fact, his
solution will work with any dimensions. The asterisk placed in front
of the parameter name means that the caller doesn’t need to
specifically provide an array… just the content. So the call:

p = Point.new(1, 2, 3, 4, 5)

Will set coords to the array [1, 2, 3, 4, 5] inside the initializer.
The initializer stores these numbers away after ensuring they all look
like floats.

    def size
        @coords.size
    end

    def []( index )
        @coords[index]
    end

    def []= ( index, x )
        @coords[index] = x
    end

    def self.random( n = 2 )
        Point.new( *(1..n).map{ 0.5 - rand } )
    end

    def +(p)
        return Point.new( *(0...size).map{|i| @coords[i] +

p[i]} )
end

    def -(p)
        return Point.new( *(0...size).map{|i| @coords[i] -

p[i]} )
end

    def *(a)
        return Point.new( *@coords.map{|x| x * a} )
    end

    def /(a)
        return Point.new( *@coords.map{|x| x / a} )
    end

    # calculate the inner product if this point with p
    def ip(p)
        (0...size).inject(0) { |sum, i| sum + @coords[i] * p[i] }
    end

    def to_s
        "(#{@coords.join(', ')})"
    end
end

Frank adds a number of arithmetic support functions for his Point
class. Looking at this reminds me of the standard library class
Vector… and it’s not surprising, as the intent of the Point class
overlaps significantly with the Vector class. I hoped someone might
make use of Vector; alas, I think the starting code I provided was
probably a bit too suggestive. Lesson learned.

Next we break down Frank’s evaluation function.

def evaluate( m, points )
    y_big = nil
    grad = nil

Here, m is a point to evaluate as the center of a circle, and points
is the list of points that need to be enclosed. This function aims to
calculate two things. First, y_big is the minimum radius of a circle
centered at m that would enclose all the points. (Actually, it’s the
square of the radius, which is a simpler calculation that avoids the
costly Math.sqrt function.) Second, grad is the gradient, which I’ll
explain later.

y_big and grad both start out invalid, but provided we have at least
one point to enclose, they will both have values at the end. Now we
loop over all the points.

    points.each do |p|
        d = (m - p)

Here, d is the delta between the circle’s center and the current
point. It’s a vector, which means it has direction (pointing at m from
p) and length (the distance from p to m).

        y = d.ip( d )

Function ip is documented by Frank as the inner product (also known as
the dot product). The inner product has a number of uses. In this
particular case, taking the inner product of d against itself gives us
the squared length of d. If we come back to the idea that we’re trying
to find a circle, if m is the center and p is on the circle, then the
length of d would be the radius. So y is the radius squared.

        if y_big.nil? or y > y_big
            y_big = y
            grad = d*2
        end
    end

Now we check to see if y is the largest squared radius we’ve found so
far; if so, remember it. Otherwise, keep the existing largest squared
radius y_big. Hit the end of the loop and keep checking points.

In order to enclose all the points, we remember the largest squared
radius y_big. One way to shrink that radius in a future iteration is
to move the center towards the point farthest away from our current
center m. Above, we mentioned that d was a vector that pointed away
from p towards m. That’s exactly the direction we’re looking for…
well, exactly the opposite direction, but that’s easily remedied by
using subtraction later rather than addition. So d is scaled (by two,
maintaining direction but with length equal to the diameter) and
assigned to grad.

    return [y_big, grad]
end

Once all points have been examined, we return the largest squared
radius and the corresponding gradient. Now let’s look at Frank’s final
function that makes use of this evaluation.

def encircle( points,
              x_start = Point.new( *([0]*points[0].size) ),
              max_iter = 100 )

Frank’s encircle function takes extra arguments, but provides sensible
defaults so that it still matches the quiz’s requirements. x_start is
a point at the origin (i.e. all components are zero) which Frank uses
as the initial guess for the circle’s center. max_iter is the number
of iterations performed in an attempt to narrow down the best answer.

    x = x_start
    y, g = nil, nil

    for k in 1..max_iter do
        y, g = evaluate( x, points )     # step 1
        x = x - g/k                      # step 2
    end

We run in a loop max_iter times, repeating two simple steps.

First, we evaluate the “fitness” of the current circle center x given
the list of points that we need to enclose. This function we described
above, and gives us back y as the squared radius of the circle and g
as the gradient used for improving upon our initial guess.

Second, we compute a new circle center x by moving it along the
gradient. As mentioned earlier, the gradient points toward the
circle’s center and away from the furthest point; by subtracting it
here, that moves the new circle’s center towards the furthest point in
the hope that will reduce the largest radius.

If we kept moving the circle’s center about by the full length of the
gradient, we’d keep jumping around without ever narrowing in on a
solution. The answer to this problem is to scale down the gradient by
k which increases each iteration. So initial jumps are large, but each
additional jump is less and less.

    return Circle.new(x, Math.sqrt(y))
end

After all iterations are done, we now have a circle. Since y has
always been the square of the radius, we finally take its square root
to get the true radius.

Note that this is not an exact solution to the quiz problem, but it
has shown to consistently come reasonably close to the exact solution.
And sometimes reasonably close is sufficient. The reason I picked
Frank’s solution is not for its exactness nor its speed, but because
it is an elegant demonstration of a general technique. If you can
provide a fitness function (in our case, the circle radius) and some
sort of gradient, you can consider this technique for other problems.

Thanks to everyone who participated in this first quiz of mine.
Tomorrow we’ll have something very simple that’ll get everyone saying,
“Hello!”.

[1]
http://www.personal.kent.edu/~rmuhamma/Compgeometry/MyCG/CG-Applets/Center/centercli.htm
[2]
http://www.cgal.org/Manual/3.2/doc_html/cgal_manual/Optimisation_ref/Class_Min_sphere_of_spheres_d.html
[3] Prof. Dr. Bernd Gärtner
[4] Minimal Enclosing Circle Problem
[5]
http://citeseer.ist.psu.edu/cache/papers/cs/30066/http:zSzzSzpage.inf.fu-berlin.dezSz~svenzSzown_workzSzmin-circle_impl_tr-b-98-04.pdf/smallest-enclosing-circles-an.pdf
[6] http://citeseer.ist.psu.edu/235065.html