Velocity Reviews > Ruby > [QUIZ] The Smallest Circle (#157)

# [QUIZ] The Smallest Circle (#157)

Frank Fischer
Guest
Posts: n/a

 02-17-2008
Hi,

here's my solution of this nice quiz. It's heavily inspired by the last
week's quiz, i.e. solving a non-linear equation. This time the program
solves the non-smooth convex optimization problem

min{max{|p_i-x|^2: p_i the given points}: x in R^k},

where k=1,2,3,... (dimension), using a simple subgradient algorithm.
Therefore, it works for an arbitrary number of points
(up to 10000 in about 40s on my 5 year old notebook) and an
arbitrary dimension (number of coordinates) of the points. The solutions
are usually a good approximation of the optimal ones.

================================================== ==============
# extends point-class with an arbitrary dimension
# and simple vector-operations
class Point
def initialize( *coords )
@coords = coords.map{|x| x.to_f}
end

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

def to_s
end
end

def generate_samples(n, dim = 2)
(1..n).map { Point.random(dim) }
end

# calculate value and subgradient of
# f(x,y) = max{(x-x_i)^2 + (y-y_i)^2: (x_i,y_i) in points}
def evaluate( m, points )
y_big = nil
points.each do |p|
d = (m - p)
y = d.ip( d ) # y = sum (m_i-p_i)^2
if y_big.nil? or y > y_big
y_big = y
end
end

end

# with given starting-point and number of iterations
def encircle( points,
x_start = Point.new( *([0]*points[0].size) ),
max_iter = 100 )
x = x_start
y, g = nil, nil

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

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

puts encircle( generate_samples(10000, 3) )
================================================== ==============

Bye,
Frank

Douglas A. Seifert
Guest
Posts: n/a

 02-17-2008
Mine is here:

http://www.pastie.org/153343

Lifted mostly from http://www.cs.mcgill.ca/~cs507/projects/1998/jacob/.
Classic Convex Hull based solution that seems to scale O(n).

-Doug Seifert

Justin Ethier
Guest
Posts: n/a

 02-17-2008
[Note: parts of this message were removed to make it a legal post.]

Hello,

Had a lot of fun with this one My solution works by taking the average of
all points as the center, and the farthest outlying point as the radius. An
optimal solution is then found by moving the circle center towards the
outlier as the radius is reduced, while all points still are within the
radius. The same principles could be applied for higher dimensions as well.

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

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

def to_s
end
end

# Calculate distance between points
def distance(pt_a, pt_b)
Math.sqrt((pt_a.x - pt_b.x) * (pt_a.x - pt_b.x) +
(pt_a.y - pt_b.y) * (pt_a.y - pt_b.y))
end

# Determine if given points are all within a circle
def inside_circle?(points, circle)
for point in points
dist = distance(point, circle.center)
return false if dist > circle.radius
end
true
end

def encircle(points) # takes array of Point objects
# find the average midpoint of all points
mid = Point.new(
points.inject(0){|sum, pt| sum += pt.x} / (points.size * 1.0),
points.inject(0){|sum, pt| sum += pt.y} / (points.size * 1.0))

# sort points by longest distance from midpoint
# longest point (index 0) is the initial radius
points.sort!{|a,b| distance(mid, a) <=> distance(mid, b) }.reverse!

# Taking the average midpoint does not necessarily work because the points
may
# be weighted more heavily to one side. We correct for this by sliding the
circle
# along the line from its original center to the outmost point, decreasing
the
# radius as we go. Keep doing this until the circle can be made no
smaller.
point = points[0]
slope = (mid.y - point.y) / (mid.x - point.x)
new_pt, delta, sign = mid, 0.0, 1.0
sign = -1.0 if mid.x > point.x
while inside_circle?(points, Circle.new(new_pt, distance(new_pt, point)))
mid = new_pt
delta += 0.000001 * sign
new_pt = Point.new(mid.x + delta, mid.y + (slope * (delta)))
end

return Circle.new(mid, distance(mid, point))
end

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

Here are Pasties of the program and a RMagick script to generate images:

Algorithm: http://pastie.caboo.se/153387
Graphics code: http://pastie.caboo.se/153388
Images:

Thanks,

Justin

On Feb 15, 2008 8:19 AM, Matthew D Moss <(E-Mail Removed)> wrote:

> The three rules of Ruby Quiz 2:
>
> 1. Please do not post any solutions or spoiler discussion for this
> quiz until 48 hours have passed from the time on this message.
>
> 2. Support Ruby Quiz 2 by submitting ideas as often as you can! (A
> permanent, new website is in the works for Ruby Quiz 2. Until then,
> please visit the temporary website at <http://
>
> 3. Enjoy!
>
> Suggestion: A [QUIZ] in the subject of emails about the problem
> helps everyone
> quiz message,
> if you can.
>
> -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
> =-=-=-
>
> The Smallest Circle
> by Matthew Moss
>
> Your task this week sounds simple enough, but may be more difficult
> than it first appears. Given a set of points on a plane, your goal is
> to find the smallest circle that encloses those points.
>
> You are to provide a function, *encircle*, that takes an array of
> points and returns the smallest circle surrounding those points.
> Start with the following base code and extend as needed to solve the
> problem:
>
> class Point < Struct.new(, :y)
> def self.random
> Point.new(rand, rand)
> end
>
> def to_s
> "(#{x}, #{y})"
> end
> end
>
> class Circle < Struct.new(:center, :radius)
> def to_s
> end
> end
>
> def encircle(points) # takes array of Point objects
> # returns a Circle object
> end
>
>
> I will be running several tests on the submitted solutions, with
> various point sets, to see how well they perform at this task. I
> recommend you you test your algorithm with a variety of sample sets,
> from small sets consisting of just 1-5 points, up to medium and
> larger sets, containing a few thousand points.
>
> To generate an array of random points, start with the above code and
>
> def generate_samples(n)
> (1..n).map { Point.random }
> end
>
>
> And then you may test your implementation like this:
>
> # encircle 10 random points
> puts encircle( generate_samples(10) )
>
>
> As mentioned, this one may be more difficult than it seems. Feel free
> to estimate the smallest circle, if you get stuck on getting the
> exact solution.
>
>

Frank Fischer
Guest
Posts: n/a

 02-17-2008
On 2008-02-17, Justin Ethier <(E-Mail Removed)> wrote:
> Hello,
>
> Had a lot of fun with this one My solution works by taking the average of
> all points as the center, and the farthest outlying point as the radius. An
> optimal solution is then found by moving the circle center towards the
> outlier as the radius is reduced, while all points still are within the
> radius. The same principles could be applied for higher dimensions as well.

Nice solution, this was also one of my first ideas. The problem with
this approach is if many points are close to each other and only few
points are far away. In this case the center is near that mass of points
and the line-search may fail to find an optimal solution. For example:

puts encircle( ([[0.0,0.0]] * 100 + [[0.0,1.0], [1.0,0.0]]).map{|p|
Point.new(*p)})

(100 points at (0,0), one at (1,0) and one at (0,1)) gives a solution of
but optimal would be (1/2, 1/2) and radius 1/sqrt(2).

Frank

ThoML
Guest
Posts: n/a

 02-17-2008
> Had a lot of fun with this one My solution works by taking the average of
> all points as the center, and the farthest outlying point as the radius. An
> optimal solution is then found by moving the circle center towards the
> outlier as the radius is reduced

Well, I didn't have the time to check out how to do this the right way
(as Douglas did it). I nevertheless tried several approaches. This one
has some imprecision built-in but it turns out (in comparison with
Douglas' solution) that for random point sets, it performs more or
less well (for an incorrect semi-solution).

I also attach some code to run several solution against each other,
which is fun if you exclude the precise solutions. This requires the
script lib from http://redshift.sourceforge.net/script/.

Cool quiz BTW that made me refresh my knowledge on triangles etc. I
algorithms.

Regards,
Thomas.

My semi-solution that uses a similar approach as Justin:

def encircle(points) # takes array of Point objects
points = points.uniq
return if points.nil? or points.empty?

x = points.inject(0.0) {|a, p| a += p.x} / points.size
y = points.inject(0.0) {|a, p| a += p.y} / points.size
a = Point.new(x, y)
return Circle.new(a, 0) unless points.size > 1

f = points.sort_by {|p| distance(a, p)}[-1]
df = distance(f, a)
b = med(f, a)
e = 1e-12
1000.times do
db = distance(f, b)
if points.all? {|p| distance(b, p) <= db + e}
da = distance(f, a)
if (da - db).abs < e
return Circle.new(b, db)
else
a, b = b, med(f, b)
end
else
b = med(a, b)
end
end
raise RuntimeError
end

def med(a, b)
Point.new((b.x - a.x) / 2 + a.x, (b.y - a.y) / 2 + a.y)
end

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

The code to run the tests:

#!/usr/bin/env ruby19

require 'script'

files = [
's_tml.rb',
's_frank.rb',
's_justin.rb',
's_douglas.rb',
]

def generate_samples(n)
(1..n).map { [rand - 0.5, rand - 0.5] }
end

# Numer of samples
# N = 10
N = 50
# N = 100
# N = 200

# Tolerance
E = 1e-10

# Your favourite solution, you might want to change this
F = 's_tml.rb'

sets = Array.new(N) {|i| generate_samples(i + 2)}
solutions = Hash.new {|h, k| h[k] = {}}

# Best choice per sample
winners = Hash.new {|h,k| h[k] = 0}

# Number of invalid solutions
losers = Hash.new {|h,k| h[k] = 0}

puts
files.each do |f|
print f
sets.each_with_index do |ps, i|
if i % 10 == 0
print '.'; STDOUT.flush
end
ps = ps.map {|p| script:oint.new(*p)}
v = script.encircle(ps)
if ps.any? {|p| distance(v.center, p) > v.radius + E}
losers[f] += 1
print 'x'; STDOUT.flush
end
solutions[i][f] = v
end
puts
end
puts

f_dist = 0

puts " %s" % files.map {|f| '%10s' % f}.join(' ')
solutions.each do |i, s|
winner = files.sort_by {|f| s[f] ? s[f].radius : 100000000}[0]
winners[winner] += 1
f_dist = f_d if f_d > f_dist
puts '%5d %10s %s (%6.4f)' % [
i,
'(%s)' % winner,
files.map {|f| s[f] && '%10.4f' % s[f].radius}.join(' '),
f_d,
]
end
puts

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

puts "Losers (invalid solutions):"
losers.each do |f, n|
puts '%10s %d' % [f, n]
end
puts

puts "Favourite vs Winner:"
puts 'Max distance: %6.4f' % f_dist
puts

Bill Kelly
Guest
Posts: n/a

 02-17-2008
Greetings,

Not sure if my solution is worth posting.... I took a naive
axis-aligned bounding box approach.

http://tastyspleen.net/~billk/ruby/q..._naive_aabb.rb

I knew it would be inaccurate, but I was curious how bad it
would be, so I wrote an OpenGL/GLUT based visualizer:

http://tastyspleen.net/~billk/ruby/q...sualization.rb

(The visualizer will work with any quiz solution that
provides the encircle() and generate_samples(n) methods,
and Point#x and Point#y.)

Occasionally the naive approach lucks out with a dataset
that allows it to look good:

http://tastyspleen.net/~billk/ruby/q...le/157_vis.png

But more often it looks something like this:

http://tastyspleen.net/~billk/ruby/q..._vis_naive.png

I had planned to improve my solution by taking an approach
similar to Justin Ethier's. . . . By iterating and "moving the
circle center towards the outlier as the radius is reduced,
while all points still are within the radius."

I may still code that if I have time... but then again, it's

Regards,

Bill

ThoML
Guest
Posts: n/a

 02-17-2008
And again I reply to my own post. The test code requires a distance()
method:

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

In case somebody cares.

Lionel Bouton
Guest
Posts: n/a

 02-17-2008
Alex Shulgin wrote:
>> http://tastyspleen.net/~billk/ruby/q...rcle/157_small...
>>

>
> I'm trying to utilize gnuplot for visualization. With any luck I can
> post the visual part without spoiling the quiz.
>

Glad to see that I'm not the sole old-school minded

Lionel

Lionel Bouton
Guest
Posts: n/a

 02-17-2008
Hi,

here's my solution. I initially thought that the best circle would be
either defined by a diameter (2 points) or one combinations of 3 points
through which it passes. The idea was that given three points defining a
circle a fourth point is either in the circle or defines a circle with 2
from the 3 originals that includes the remaining point.

This would be neat as I wasn't so happy about a gradient walking
solution. In this case it can be proven that you can deterministically
maximums) but it's not the general case. Unfortunately it can be proven
that such a circle can encircle all points, but not that it's the
smallest (and my experiments proved that it isn't). Finally this is even
slower than my proposed solution below, ...

So here's a more brute force method : minimizing the maximum distance
with some optimizations (heuristic for initial values and a not so naive
2D-space walker). There's probably some more optimizations for walking
the gradient but I'm happy with the current speed (~10s for 10000 points
on my 1.06GHz Core Duo).

This seems similar in spirit to Frank's solution, but less clean from
the Math PoV and more optimized for speed.

It outputs the solution and creates a gnuplot plotting file for you to
verify (you have the path used by the algorithm to find the center). I'm
no Gnuplot expert so unfortunately fancy colors are out and crappy
legends in...

Lionel

=============================================
include Math

# Uses 10 points by default
POINT_COUNT = (ARGV[0].to_i == 0) ? 10 : ARGV[0].to_i

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

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

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

class PotentialCenter < Struct.new(, :y, oints)
# Square distance with a point
def square_distance(point)
((point.x - x) ** 2) + ((point.y - y) ** 2)
end

# Maximum square distance with my points
# It's cached because it's called several times and can be very costly
def max_square_distance
@sq_dist ||= (points.map { |point| square_distance(point) }.max)
end

# Compute the best move with given distance (ro)
# and angles (thetas)
# returns nil if none is better than the current position
def best_move(ro, thetas)
potential_moves = thetas.map { |theta|
new_x = x + (ro * cos(theta))
new_y = y + (ro * sin(theta))
PotentialCenter.new(new_x,new_y,points)
}
choose_move_1(potential_moves)
end

# Dumber, faster
def choose_move_1(potential_moves)
potential_moves.detect { |move|
max_square_distance > move.max_square_distance
}
end

# Look for the shortest path, slower (~1.4x on average with many points)
def choose_move_2(potential_moves)
potential_moves.reject! { |move|
max_square_distance <= move.max_square_distance
}
potential_moves.min { |c1,c2|
c1.max_square_distance <=> c2.max_square_distance
}
end

# Check that the given offset is big enough to move
# the current position
# (floating point rounding errors prevent infinite precision...)
def can_be_moved_by(offset)
((x + offset) != x) || ((y + offset) != y)
end

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

def to_s
end
end

def chose_original_center_and_move_amount(points)
# Start with a good potential (around the middle of the point cloud)
max_x = points.map{|p| p.x}.max
min_x = points.map{|p| p.x}.min
x = (max_x + min_x) / 2
max_y = points.map{|p| p.y}.max
min_y = points.map{|p| p.y}.min
y = (max_y + min_y) / 2
center = PotentialCenter.new(x, y, points)
# The original displacement value is based on the fact that a truly
# random distribution gets a mean value with a precision of this order.
# The actual center is probably reachable with a few of these steps
move_amount = ((max_y - min_y) + (max_x - min_x)) / points.size
return [center, move_amount]
end

# Look for the best center by minimizing
# its maximum distance with the given points
# This is the same as minimizing its square
# (the maximum of the square distances)
# This can be seen as a dichotomy on a 2-dimensional space
def encircle(points)
center, ro = chose_original_center_and_move_amount(points)
# We'll move with polar coordinates (distance + angle : ro/theta)
# How many directions do we follow?
# each new angle creates a new (costly) max_square_distance call
# so we take the minimum needed to move around and climb the
angle_count = 3
thetas = (1..angle_count).map { |c| (2 * c * Math:I) / angle_count}
iter = 0
# Trace our moves
center_list = []

loop do
center_list << center
iter += 1
puts "#{iter}: #{center.max_square_distance}"
# Is there a best move available ?
if new_center = center.best_move(ro, thetas)
center = new_center
else
ro /= 2
# We stop when ro becomes too small to move the center
break unless center.can_be_moved_by(ro)
end
end
puts "Circle found in #{iter} iterations"
return [ Circle.new(center, sqrt(center.max_square_distance)),
center_list ]
end

points = Point.generate_samples(POINT_COUNT)
t = Time.now
circle, center_list = encircle(points)
puts circle
puts "Found in: #{Time.now - t}s"

# Generate a GnuPlot command file to see the result
f = File.open('plotcommand.cmd', 'w')
f.print <<EOF
set parametric
set size square
set xrange [-0.2:1.2]
set yrange [-0.2:1.2]
set multiplot
plot "-"
#{points.map { |p| "#{p.x} #{p.y}" }.join("\n")}
end
plot "-" with lines
#{center_list.map { |p| "#{p.x} #{p.y}" }.join("\n")}
end
plot [0:2*pi]
set nomultiplot
EOF
f.close

puts <<EOF
use 'gnuplot -persist plotcommand.cmd' to see
the circle, points and path to the center
EOF

Bill Kelly
Guest
Posts: n/a

 02-17-2008

From: "Bill Kelly" <(E-Mail Removed)>
>
> I had planned to improve my solution by taking an approach
> similar to Justin Ethier's. . . . By iterating and "moving the
> circle center towards the outlier as the radius is reduced,
> while all points still are within the radius."
>
> I may still code that if I have time... but then again, it's

Ah, what the heck... Here 'tis..

http://tastyspleen.net/~billk/ruby/q...lest_circle.rb

# Algorithm used: First enclose the points in a simple
# axis-aligned bounding box, and choose the center of the
# box as the initial center of the circle. Then, binary
# search to determine how far we can move the center point
# toward the furthest outlying point, while keeping the
# rest of the points still inside the circle.

Regards,

Bill