«

On Feb 20, 7:19=A0pm, Lionel Bouton <lionel-subscript...@bouton.name>
wrote:
> John Joyce 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 =A0not 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(, :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(, :y)
def self.random(type=0)
case type
when 1
ro = rand / 2
theta = Math:I * 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:oint.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:I * 2 * rand
[ ro * Math::cos(theta), ro * Math::sin(theta)]
}
},
"Random on circle" =>
(1..num_runs).map {(1..num_samples).map {
theta = Math:I * 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:I * 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 Kelly, Thomas ML and Lionel Bouton; 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:

> 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....
>>
>> 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 Seifert did provide
such a solution.

The solution I am going to examine here is that of Frank Fischer,
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/~rmuham.../centercli.htm
[2] http://www.cgal.org/Manual/3.2/doc_h...spheres_d.html
[3] http://www.inf.ethz.ch/personal/gaertner/miniball.html
[4] http://www.cs.mcgill.ca/~cs507/proje...b/problem.html
[5]
http://citeseer.ist.psu.edu/cache/pa...circles-an.pdf
[6] http://citeseer.ist.psu.edu/235065.html

Apologies for the typos in the summary's code bits. Google Groups has
decided it wants to protect a few email addresses from spam, despite
there being no email addresses. (It's seeing the @ character of the
instance member variables.)

Where you see "..." immediately preceding "@", realize those dots
should not be there. In one of the three cases (the to_s method),
those dots should be replaced with a left curly brace "{".

For the future, I'll send the email directly to the list rather than
going through Google.