From: Frank Fischer <frank.fischer@...2001.tu-chemnitz.de>
Date: 2008-02-18T00:00:03+09:00
Subject: Re: [QUIZ][SOLUTION] The Smallest Circle (#157)

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

class Circle < Struct.new(:center, :radius)
    def to_s
        "{center:#{center}, radius:#{radius}}"
    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
    grad = 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
            grad = d*2
        end
    end

    return [y_big, grad]
end

# perform simple subgradient algorithm
# 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