Edit page
Note: The code on this page may fail. (Unknown error.)

Church version:

(define (make-plate pos size table)
  (list pos size table))
(define plate-position first)
(define plate-size second)
(define plate-table third)

(define (make-table size num-plates)
  (list size num-plates))
(define table-size first)
(define table-num-plates second)

(define (plate-area plate)
  (* my-pi (* (plate-size plate) (plate-size plate))))

(define (total-area plates)
  (apply + (map plate-area plates)))

(define (table-area table)
  (* my-pi (* (table-size table) (table-size table))))

(define (vec-dist xy1 xy2)
  (let* ([x (list-ref xy1 0)]
         [y (list-ref xy1 1)]
         [z (list-ref xy2 0)]
         [w (list-ref xy2 1)])
    (expt (+ (expt (- x z) 2) (expt (- y w) 2)) 0.5)))

(define (dist-from-origin xy)
  (vec-dist xy (list 0 0)))

(define (circular-area r)
  (* my-pi (expt r 2)))

(define (area-outside-table plate)
  (let* ([d (dist-from-origin (plate-position plate))]
         [table-limit (table-size (plate-table plate))]
         [difference (max 0 (- d table-limit))])
    (circular-area difference)))

(define (overlap-area p1 p2)
  (let* ([dist-between (vec-dist (plate-position p1) (plate-position p2))]
         [overlap-amt (max 0 (- (+ (plate-size p1) (plate-size p2)) dist-between))])
    (circular-area overlap-amt)))

(define (select-k-subsets n l)
  (letrec ([loop (lambda (l ln n prev-els accum)
                   (cond
                     ((<= n 0) (cons prev-els accum))
                     ((< ln n) accum)
                     ((= ln n) (cons (append l prev-els) accum))
                     ((= ln (+ 1 n)) 
                      (letrec ([fold (lambda (l seen accum)
                        (if (null? l) accum
                          (fold (cdr l) (cons (car l) seen)
                                (cons
                                  (append (cdr l) seen)
                                  accum))))])
                        (fold l prev-els accum)))
                     ((= n 1)
                      (letrec ([fold (lambda (l accum)
                        (if (null? l) accum
                          (fold (cdr l) (cons (cons (car l) prev-els) accum))))])
                        (fold l accum)))
                     (else
                       (loop (cdr l) (- ln 1) n prev-els
                             (loop (cdr l) (- ln 1) (- n 1) (cons (car l) prev-els) accum)))))])
    (loop l (length l) n '() '())))

(define (pairs xs) (select-k-subsets 2 xs))

(define (randint low high)
  (+ low (sample-integer (- high low))))

(define (sample-table)
  (let* ([num-plates (randint 1 10)]
         [table-size (uniform 15 25)])
    (make-table table-size num-plates)))

(define (sample-plate table)
  (let* ([posx (uniform -10 10)]
         [posy (uniform -10 10)]
         [plate-size (uniform 0.1 2.0)])
    (make-plate (list posx posy) plate-size table)))

(define (tables-and-plates)
  (let* ([table (sample-table)]
         [plates (map (lambda (i) (sample-plate table)) (iota (table-num-plates table)))]
         [occupy-area (softeq 0.7 (/ (total-area plates) (table-area table)))]
         [inside (map (lambda (p) (softeq 0.0 (area-outside-table p))) plates)]
         [non-overlap (map (lambda (pq) (softeq 0.0 (overlap-area (first pq) (second pq)))) (pairs plates))])
    (list (apply + (list occupy-area (apply + inside) (apply + non-overlap))) table plates)))

(for-each display
          (mh-query 100 100
                    (define asn (tables-and-plates))
                    (list 'score (car asn) 'num-plates (table-num-plates (cadr asn)))
                    true))

Terra version:

terralib.require("prob")
local mem = terralib.require("mem")
local Vec = terralib.require("linalg").Vec
local Vector = terralib.require("vector")
local inheritance = terralib.require("inheritance")
local rand = terralib.require("prob.random")
local cmath = terralib.includec("math.h")

local C = terralib.includecstring [[
#include <stdio.h>
]]

local Vec2 = Vec(double, 2)

local struct Circle { pos: Vec2, rad: double }
terra Circle:area() return [math.pi]*self.rad*self.rad end
terra Circle:intersectArea(other: &Circle)
	var r = self.rad
	var R = other.rad
	var d = self.pos:dist(other.pos)
	-- No intersection
	if d > r+R then
		return 0.0
	end
	if R < r then
		r = other.rad
		R = self.rad
	end
	-- Complete containment
	if d < R-r then
		return [math.pi]*r*r
	end
	var d2 = d*d
	var r2 = r*r
	var R2 = R*R
	var x1 = r2*cmath.acos((d2 + r2 - R2)/(2*d*r))
	var x2 = R2*cmath.acos((d2 + R2 - r2)/(2*d*R))
	var x3 = 0.5*cmath.sqrt((-d+r+R)*(d+r-R)*(d-r+R)*(d+r+R))
	return x1 + x2 - x3
end

local struct Table { plates: Vector(Circle) }
inheritance.staticExtend(Circle, Table)
terra Table:__construct() mem.init(self.plates) end
terra Table:__destruct() mem.destruct(self.plates) end
terra Table:__copy(t: &Table) self.plates = mem.copy(t.plates) end
mem.addConstructors(Table)

local function layoutModel()

	-- Soft equality factor function
	local softEq = macro(function(x, target, softness)
		return `[rand.gaussian_logprob(real)](x, target, softness)
	end)

	-- Generate random tables (with plates)
	local plateNums = Vector.fromNums(0, 1, 2, 3, 4)
	local makeTables = pfn()
	makeTables:define(terra(numTables: int, tables: &Vector(Table)) : {}
		if numTables > 0 then
			-- Generate table with random posiiton and size
			tables:resize(tables.size + 1)
			var newTable = tables:backPointer()
			newTable.pos = Vec2.stackAlloc(uniform(0.0, 50.0), uniform(0.0, 50.0))
			newTable.rad = uniform(5.0, 15.0)

			-- Make random number of plates with random 
			--    position and size
			var numPlates = int(uniformDraw(&plateNums))
			newTable.plates:resize(numPlates)
			for i=0,numPlates do
				var posvariance = newTable.rad / 3.0
				var perturb = Vec2.stackAlloc(gaussian(0.0, posvariance), gaussian(0.0, posvariance))
				newTable.plates(i).pos = newTable.pos + perturb
				newTable.plates(i).rad = uniform(1.0, 2.0)

				-- Encourage plate to be fully on the table
				var area = newTable.plates(i):area()
				var isectArea = newTable.plates(i):intersectArea(newTable)
				factor(softEq(isectArea, area, 2.0))
			end

			-- Encourage plates not to overlap
			if numPlates > 0 then
				for i=0,numPlates-1 do
					for j=i+1, numPlates do
						factor(softEq(newTable.plates(i):intersectArea(newTable.plates:getPointer(j)), 0.0, 0.1))
					end
				end
			end

			-- Recursively generate remaining tables
			makeTables(numTables-1, tables)
		end
	end)

	-- Program to perform inference on
	return terra()
		-- Generate random number of tables
		var tables = [Vector(Table)].stackAlloc()
		var numTables = poisson(6)
		makeTables(numTables, &tables)

		-- Encourage tables not to overlap
		for i=0,tables.size-1 do
			for j=i+1, tables.size do
				factor(softEq(tables(i):intersectArea(tables:getPointer(j)), 0.0, 0.1))
			end
		end

		return tables
	end
end

References: