# Example:Total circles area/Racket

## Racket

<img src="timb.net/images/RC/Total-circles-area/TCA-all-circles.png" alt="All the circles"/>
<img src="timb.net/images/RC/Total-circles-area/TCA-split-circles.png" alt="Showing the path edge from split-circles"/>

Below are an analytical solution in `TCA-analytical.rkt`, and a Monte-carlo solution in `TCA-monte-carlo.rkt`. Common code (notably types) is in `TCA-types.rkt`, and `TCA-task.rkt` (we don't need to list the circle coordinates more than once).

In order to get my head around the algorithm, I have also thrown together a `TCA-draw.rkt` which (if I can get RC to let me upload piccies) produces the attached images (otherwise you'll have to go to timb.net for them).

### `TCA-types.rkt`

This solution (with the exception of `TCA-draw.rkt`) is in typed/racket. This file defines the types, and geometrical arithmetic. In fact it encroaches on the analytical solution's magic -- but it's not so easy to separate concerns here.

`#lang typed/racket;;; TCA-types : basic geometric types, and some basic geometric operations upon them ;; Arguably even the split circles / other path operations are general enough to (provide (all-defined-out)) ; to external "drawing" module (struct Vec    ((x : Real) (y : Real)))(struct Angle2 ((a0 : Angle) (a1 : Angle))) ; Angles of the start and end points of the circle arc.(struct Circle ((c : Vec) (r : Real)))(struct Arc    ((c : Circle) (a2 : Angle2))) (define-type Angle Real)(define-type Vecs (Listof Vec))(define-type Arcs (Listof Arc))(define-type Angles (Listof Angle))(define-type Circles (Listof Circle)) (define-type Geometric (U Vec Circle Arc)) (: v-cross : Vec Vec -> Real)(: v-dot   : Vec Vec -> Real)(: v+      : Vec Vec -> Vec)(: v-      : Vec Vec -> Vec)(: v-len   : Vec -> Real)(: v-dist  : Vec Vec -> Real)(: v-scale : Vec Real -> Vec)(: v-norm  : Vec -> Vec) (define/match (v-cross v0 v1)  [((Vec a b) (Vec c d)) (- (* a d) (* b c))]) (define/match (v-dot v0 v1)  [((Vec a b) (Vec c d)) (+ (* a c) (* b d))]) (define/match (v+ v0 v1)  [((Vec a b) (Vec c d)) (Vec (+ a c) (+ b d))]) (define/match (v- v0 v1)  [((Vec a b) (Vec c d)) (Vec (- a c) (- b d))]) (define (v-len v)  (cast (sqrt (v-dot v v)) Real)) (define (v-dist a b)  (v-len (v- a b))) (define/match (v-scale v s)  [((Vec x y) s) (Vec (* x s) (* y s))]) (define/match (v-norm v)  [((and (Vec x y) (app v-len l))) (Vec (/ x l) (/ y l))]) (: v-angle : Vec -> Angle)(: a-norm  : Angle -> Angle) (define/match (v-angle v)  [((Vec x y)) (atan y x)]) (define (a-norm a)  (cond [(> a pi) (- a (* 2 pi))] [(< a (- pi)) (+ a (* 2 pi))] [else a])) (: circle-cross : Circle -> (Circle -> Angles))(define/match ((circle-cross C0) C1)  [((Circle c0 r0) (Circle c1 r1))   (define d (v-dist c0 c1))   (cond     [(>= d (+ r0 r1)) null]     [(<= d (abs (- r0 r1))) null]     [else      (define s (/ (+ r0 r1 d) 2))      (define a (sqrt (* s (- s d) (- s r0) (- s r1))))      (define h (* 2 (/ a d)))      (define dr (v- c1 c0))      (define r0^2 (sqr r0))      (define dr-ang (v-angle dr))      (define ang (+ dr-ang (if (> (+ r0^2 (sqr d)) (sqr r1)) 0 pi)))      (define da (cast (asin (/ h r0)) Real))       (map a-norm (list (- ang da) (+ ang da)))])]) (: arc-point  : Circle Angle -> Vec)(: arc-start  : Arc -> Vec)(: arc-mid    : Arc -> Vec)(: arc-end    : Arc -> Vec)(: arc-centre : Arc -> Vec)(: arc-area   : Arc -> Real) (define/match (arc-point crc arc)  [((Circle c r) a) (v+ c (Vec (* r (cos a)) (* r (sin a))))]) (define/match (arc-start arc)  [((Arc c (Angle2 a0 _)))  (arc-point c a0)])(define/match (arc-mid arc)    [((Arc c (Angle2 a0 a1))) (arc-point c (/ (+ a0 a1) 2))])(define/match (arc-end arc)    [((Arc c (Angle2 _  a1))) (arc-point c a1)])(define/match (arc-centre arc) [((Arc (Circle c _) _)) c]) (define arc-area  (match-lambda    [(Arc (Circle _ r) (Angle2 a0 a1))     (* (sqr r) (- a1 a0) 1/2)])) (: tri-area : Vec Vec Vec -> Real)(define (tri-area a b c) (/ (v-cross (v- b a) (v- c b)) 2)) ;; Boundaries -- needed both for drawing and monte-carlo(: x-min : Any -> Real)(: x-max : Any -> Real)(: y-min : Any -> Real)(: y-max : Any -> Real) ;; values is too general for the type checker... so (argmin values l) is not possible(: real-min : (Listof Real) -> Real)(: real-max : (Listof Real) -> Real) (define (real-min l)  (foldl min +Inf.0 l)) (define (real-max l)  (foldl max -Inf.0 l)) (define/match (x-min o)  [((list os ...)) (real-min (map x-min os))]  [((Vec x y)) x]  [((Circle (Vec x y) r)) (- x r)]  [((Arc c _)) (x-min c)]  [(else) +Inf.0]) (define/match (x-max o)  [((list os ...)) (real-max (map x-max os))]  [((Vec x y)) x]  [((Circle (Vec x y) r)) (+ x r)]  [((Arc c _)) (x-max c)]  [(else) -Inf.0]) (define/match (y-min o)  [((list os ...)) (real-min (map y-min os))]  [((Vec x y)) y]  [((Circle (Vec x y) r)) (- y r)]  [((Arc c _)) (y-min c)]  [(else) +Inf.0]) (define/match (y-max o)  [((list os ...)) (real-max (map y-max os))]  [((Vec x y)) y]  [((Circle (Vec x y) r)) (+ y r)]  [((Arc c _)) (y-max c)]  [(else) -Inf.0])`

### `TCA-draw.rkt`

This does not depend on any algorithms, and allows visualisation of the geometrical objects (in fact scenarios composed of lists of said objects)

`#lang racket(require "TCA-types.rkt")(require pict)(provide (all-defined-out))(define/match (draw-object dc o)  [(dc (list o ot ...)) (draw-object dc o) (draw-object dc ot)]  [(dc (list)) (void)]  [(dc (Vec x y)) (send dc draw-point x y)]  [(dc (? string? s)) (send* dc (set-pen s 0 'solid) (set-brush s 'transparent))]  [(dc (Circle (Vec cx cy) r)) (send dc draw-ellipse (- cx r) (- cy r) (* 2 r) (* 2 r))]  [(dc (and arc (Arc (Circle (Vec cx cy) r) (Angle2 start end))))   (send dc draw-arc (- cx r) (- cy r) (* 2 r) (* 2 r) (- end) (- start))]) (define (draw-situation #:fill-first? (fill-first? #f) pict-w pict-h o0 . os)  (define o* (cons o0 os))  (define x* (x-min o*))  (define y* (y-min o*))  (define w* (- (x-max o*) x*))  (define h* (- (y-max o*) y*))   (dc   (λ (dc dx dy)     (define-values (canvas-w canvas-h) (send dc get-size))     (define scl-x (/ canvas-w w*))     (define scl-y (/ canvas-h h*))     (define scl (min scl-x scl-y))     (send* dc       (set-initial-matrix (vector scl 0 0 (- scl) 0 0))       (translate (- x*) (- (+ y* h*)))       ;; outline       (set-pen "red" 0 'solid) (set-brush "red" 'transparent)       (draw-rectangle x* y* w* h*)       ;; colour of tail objects       (set-pen   "green" 0 'solid) (set-brush "green" (if fill-first? 'solid 'transparent)))          ;; drawn in reverse, so last in list is most background     (for ((o (reverse os))) (draw-object dc o))          ;; colour of first object -- drawn last, so on top     (send* dc (set-pen "cyan" 0 'solid) (set-brush "cyan" 'transparent))     (draw-object dc o0))   pict-w pict-h)) (define (show-situation #:fill-first? (f-1st? #f) #:width (w 600) #:height (h 400) o0 . os)  (show-pict (apply draw-situation w h o0 os #:fill-first? f-1st?))) ;; Saves image of situation to file(define (file-situation file-name #:fill-first? (f-1st? #f) #:width (w 600) #:height (h 400) o0 . os)  (define bmp (pict->bitmap (apply draw-situation w h o0 os #:fill-first? f-1st?)))  (send bmp save-file file-name 'png))`

### `TCA-monte-carlo.rkt`

• If `deeper` is ≤ 0, then this does a plain old Monte Carlo approximation.
• If `deeper` > 0, then the sample space is divided into 4 (2x2), and circles

that cannot possibly affect each subsample are discarded. This allows for more focussed "inside?" testing.

Nothing more sophisticated is attempted. That would be verging on analytical.

`#lang typed/racket(require "TCA-types.rkt")(provide circles-area) (: circles-area :   Circles [#:samples Nonnegative-Integer]   [#:bounds (U #f (List Real Real Real Real))]   [#:deeper Nonnegative-Integer]   -> Real) ;; Provides naive sub-sampling when deeper > 0(define (circles-area all-cs #:samples (samples 1000) #:bounds (bounds #f) #:deeper (deeper 6))  (when (>= deeper 3) (eprintf "[~a/~a]" deeper (length all-cs)) (flush-output))  (define x (or (and bounds (first bounds))  (x-min all-cs)))  (define y (or (and bounds (second bounds)) (y-min all-cs)))  (define w (or (and bounds (third bounds))  (- (x-max all-cs) x)))  (define h (or (and bounds (fourth bounds)) (- (y-max all-cs) y)))  (define w/2 (/ w 2))  (define h/2 (/ h 2))  (define x+ (+ x w/2))  (define y+ (+ y h/2))  (define bounds-c (Vec x+ y+))  (define bounds-diagonal (v-dist (Vec x y) bounds-c))   ;; we won't go into the detail of "circles that intersect the area",  ;; just circles that can't possibly overlap  ;; (their radius is further than the diagonal of the bounds)  ;; if we want more analysis on this, then we'll use the TCA-analytical solution. No?  (: circle-threatens? : Circle -> Boolean)  (define (circle-threatens? c) (< (v-dist (Circle-c c) bounds-c) (+ (Circle-r c) bounds-diagonal)))  (define cs (filter circle-threatens? all-cs))  (cond    [(null? cs) 0]    [(not (positive? deeper))     (: in-circle? : Vec -> (Circle -> Boolean))     (define ((in-circle? v) c) (< (v-dist v (Circle-c c)) (Circle-r c)))     (: in-any-circle? : Vec -> Boolean)     (define (in-any-circle? v)       (list? (memf (in-circle? v) cs)))     (: rand-vec : -> Vec)     (define (rand-vec) (Vec (+ x (* w (random))) (+ y (* h (random)))))     (define hits (for/sum : Real ((i (in-range samples)) #:when (in-any-circle? (rand-vec))) 1))     (* hits (* w h) (/ samples))]    [else     (+ (circles-area cs #:samples samples #:bounds (list x y w/2 h/2) #:deeper (sub1 deeper))        (circles-area cs #:samples samples #:bounds (list x+ y w/2 h/2) #:deeper (sub1 deeper))        (circles-area cs #:samples samples #:bounds (list x y+ w/2 h/2) #:deeper (sub1 deeper))        (circles-area cs #:samples samples #:bounds (list x+ y+ w/2 h/2) #:deeper (sub1 deeper)))]))`

### `TCA-analytical.rkt`

This provides an analytical solution. Much of the Haskell heavy lifting is rewitten in `TCA-types.rkt`.

`#lang typed/racket(require "TCA-types.rkt")(provide circles-area         split-circles ;; I want to be able to draw this         ) (: split-1-circle : (Pair Circle Angles) -> Arcs)(define/match (split-1-circle c.angs)  [((cons c angs))   (for/list : Arcs ((a1 (in-list angs)) (a2 (in-list (cdr angs)))) (Arc c (Angle2 a1 a2)))]) (: split-circles : Circles -> Arcs)(define (split-circles cs)  (: in-circle? (Arc Circle -> Boolean))  (define/match (in-circle? a c)    [((and arc (Arc (? Circle? c1) _)) (and c2 (Circle cen2 r2)))     (and (not (equal? c1 c2)) (< (v-dist (arc-mid arc) cen2) r2))])   ;; If an arc that was part of one circle is inside *another* circle,  ;; it will not be part of the zero-winding path, so reject it.  (: not-in-any-circle? : Arc -> Boolean)  (define/match (not-in-any-circle? arc)    [(arc) (not (memf (curry in-circle? arc) cs))])   (: f : Circle -> (Pair Circle Angles))  (define (f c)    (cons c (sort (list* (- pi) pi (apply append (map (circle-cross c) cs))) <)))   (define c-angs (map f cs))   (define arcs (apply append (map split-1-circle c-angs)))   (filter not-in-any-circle? arcs)) #|Given a list of arcs, build sets of closed paths from them. If one arc's end point is no more than1e-4 from another's start point, they are considered connected.  Since these start/end pointsresulted from intersecting circles earlier, they *should* be exactly the same, but floating pointprecision may cause small differences, hence the 1e-4 error margin.  When there are genuinelydifferent intersections closer than this margin, the method will backfire, badly.|# (: arc-dist-quantum Real)(define arc-dist-quantum 1e-12) (: make-paths : Arcs -> (Listof Arcs))(define (make-paths arcs)  (: join-arcs : Arcs Arcs -> (Listof Arcs))  (define join-arcs    (match-lambda**     [(a (list)) (list a)]     [((list) (cons x xs)) (join-arcs (list x) xs)]     [(a X)      #:when (< (v-dist (arc-start (first a)) (arc-end (last a))) arc-dist-quantum)      (cons a (join-arcs null X))]     [(a (cons x xs))      #:when (< (v-dist (arc-end (last a)) (arc-start x)) arc-dist-quantum)      (join-arcs (append a (list x)) xs)]     [(a (cons x xs)) (join-arcs a (append xs (list x)))]))  (join-arcs null arcs)) ;; Slice N-polygon into N-2 triangles.(: polyline-area : Vecs -> Real)(define (polyline-area v.vs)  (match-define (cons v vs) v.vs)  (for/sum : Real ((v0 : Vec (in-list vs)) (v1 : Vec (in-list (cdr vs))))    (tri-area v v0 v1))) (: path-area : Arcs -> Real)(define (path-area arcs)  (define-values (a e)    (for/fold ((a : Real 0) (e : Vecs null))      ((arc : Arc (in-list arcs)))      (values (+ a (arc-area arc)) (append e (list (arc-centre arc) (arc-end arc))))))  (+ a (polyline-area e))) (: circles-area : Circles -> Real)(define (circles-area cs)  (apply + (map path-area (make-paths (split-circles cs)))))`

### `TCA-task.rkt`

This file does the business of showing results and whatnot.

`#lang typed/racket (require "TCA-types.rkt") (require (prefix-in analytical: "TCA-analytical.rkt"))(require (prefix-in monte-carlo: "TCA-monte-carlo.rkt")) (provide (all-defined-out)) ; to external "drawing" module (: circles Circles)(define circles  (list (Circle (Vec  1.6417233788  1.6121789534) 0.0848270516)        (Circle (Vec -1.4944608174  1.2077959613) 1.1039549836)        (Circle (Vec  0.6110294452 -0.6907087527) 0.9089162485)        (Circle (Vec  0.3844862411  0.2923344616) 0.2375743054)        (Circle (Vec -0.2495892950 -0.3832854473) 1.0845181219)        (Circle (Vec  1.7813504266  1.6178237031) 0.8162655711)        (Circle (Vec -0.1985249206 -0.8343333301) 0.0538864941)        (Circle (Vec -1.7011985145 -0.1263820964) 0.4776976918)        (Circle (Vec -0.4319462812  1.4104420482) 0.7886291537)        (Circle (Vec  0.2178372997 -0.9499557344) 0.0357871187)        (Circle (Vec -0.6294854565 -1.3078893852) 0.7653357688)        (Circle (Vec  1.7952608455  0.6281269104) 0.2727652452)        (Circle (Vec  1.4168575317  1.0683357171) 1.1016025378)        (Circle (Vec  1.4637371396  0.9463877418) 1.1846214562)        (Circle (Vec -0.5263668798  1.7315156631) 1.4428514068)        (Circle (Vec -1.2197352481  0.9144146579) 1.0727263474)        (Circle (Vec -0.1389358881  0.1092805780) 0.7350208828)        (Circle (Vec  1.5293954595  0.0030278255) 1.2472867347)        (Circle (Vec -0.5258728625  1.3782633069) 1.3495508831)        (Circle (Vec -0.1403562064  0.2437382535) 1.3804956588)        (Circle (Vec  0.8055826339 -0.0482092025) 0.3327165165)        (Circle (Vec -0.6311979224  0.7184578971) 0.2491045282)        (Circle (Vec  1.4685857879 -0.8347049536) 1.3670667538)        (Circle (Vec -0.6855727502  1.6465021616) 1.0593087096)        (Circle (Vec  0.0152957411  0.0638919221) 0.9771215985))) (define exact-value 21.5650366038563989590842249288781480183977)(: run-task : (Circles -> Real) -> Void)(define (run-task f)  (define a (f circles))  (printf "circles-area (~s): ~a [error:~a]~%" f a (- a exact-value))) (run-task analytical:circles-area)(run-task monte-carlo:circles-area) (define-type file-situation-drawer-t  (->* (String (Pairof (U Geometric String) (Listof (U Geometric String))))       (#:width Nonnegative-Integer #:height Nonnegative-Integer #:fill-first? Boolean)       Void)) (require/typed "TCA-draw.rkt" [file-situation file-situation-drawer-t])(require "TCA-analytical.rkt") (file-situation "TCA-all-circles.png" (cons (car circles) (cdr circles)))(file-situation "TCA-split-circles.png"                (cons (car circles)                      (append (cdr circles) (list "black") (split-circles circles))))`
Output:

As well as this output, there are also two image files (whether or not you see them depends on my permissions at the moment!)

```circles-area (#<procedure:circles-area>): 21.5650366038564 [error:0.0]
circles-area (#<procedure:circles-area>): 21.564722220577558 [error:-0.000314383278841035]```