Planes In Space : Random Sampling

;; My recursion tells me that n planes can divide space into
;; 2, 4, 8, 15, 26, and finally 42 regions.

;; But I'd like some sort of empirical confirmation of this result.

;; Forms in 3-space can be represented as ordered sets of three numbers
(def a-form [1 2 3])

;; As can vectors
(def a-point [2 3 5])

;; And we evaluate the form on the vector
(defn contract [form vector]
  (reduce + (map * form vector)))
;; (or if you prefer, take the dot product, or contract the tensors)

(contract a-form a-point) ;-> 23

;; Any plane can be defined in terms of a 1-form and a number.  
(def a-plane {:form [1 2 3] :val 4})

;; Now if we have a plane or a vector, we can evaluate the form on the
;; vector, and compare the result with the value. This tells us which
;; side of the plane the vector is on.

(defn side [plane point]
  (- (contract (plane :form) point) (plane :val)))


(side a-plane [2 3  5]) ;-> 19  (this point is on the positive side)
(side a-plane [2 3 -5]) ;-> -11 (on the negative side)
(side a-plane [2 3 -4/3]) ;-> 0N (in the plane itself)


;; Ok, now we need a way of taking vectors and forms at random. 
;; The cauchy distribution is easy to sample from and has nice fat tails

(defn cauchy[]
  (Math/tan (* Math/PI (- (rand) 0.5))))

(repeatedly 20 cauchy) ;-> (-0.43989542100474244 -0.6517139433588255 1.58518947555566 0.001268073580101198 3.6164981498788262 0.44928717656825584 0.3365831420089349 0.4646894211443393 0.8802485518044282 1.8146747880005754 0.1608864471756546 -0.23538854021056904 8.836583912257565 3.8174659241864703 0.5387819323291936 -0.18830386363467239 -1.0430272980416788 0.3310448308016341 -0.10735190850334911 0.3426157380908667)

(defn make-point []
  (repeatedly 3 cauchy))

(defn make-plane []
  {:form (repeatedly 3 cauchy) :val (cauchy)})

(make-point) ;-> (33.032354006369815 -29.428219536044043 -37.796430533111334)
(make-plane) ;-> {:form (-45.36910184399889 -1.6741101969009575 9.952054197916382), :val 0.9505471630252558}

(def points (repeatedly #(make-point)))
(def planes (repeatedly #(make-plane)))

;; And we'll need a function to tell us the sign of a number
(defn sign[x]
  (if (< x 0) '- '+))

(map sign [ -1 -2 -3 0 -0.5 1.3]) ;-> (- - - + - +)

;; Now if we take a set of planes and a point, 
(defn sig [point planes]
  (for [p planes] (sign (side p point))))

;; We can check which side of each plane the point is on
(sig (first points) (take 3 planes)) ;-> (+ - +)

;; Every different region gives a different signature. 
;; The more planes, the more signatures.

(count (frequencies (take 10 (map #(sig % (take 1 planes)) points)))) ;-> 2
(count (frequencies (take 10 (map #(sig % (take 2 planes)) points)))) ;-> 4
(count (frequencies (take 10 (map #(sig % (take 3 planes)) points)))) ;-> 6
(count (frequencies (take 10 (map #(sig % (take 4 planes)) points)))) ;-> 7
(count (frequencies (take 10 (map #(sig % (take 5 planes)) points)))) ;-> 7
(count (frequencies (take 10 (map #(sig % (take 6 planes)) points)))) ;-> 7

;; But the more planes we have, the smaller the smallest regions are
;; and thus the chance of a point falling in every one goes down.


;; The more points we take, the more likely we are to get one in every region

(count (frequencies (take 100 (map #(sig % (take 1 planes)) points)))) ;-> 2
(count (frequencies (take 100 (map #(sig % (take 2 planes)) points)))) ;-> 4
(count (frequencies (take 100 (map #(sig % (take 3 planes)) points)))) ;-> 7
(count (frequencies (take 100 (map #(sig % (take 4 planes)) points)))) ;-> 11
(count (frequencies (take 100 (map #(sig % (take 5 planes)) points)))) ;-> 18
(count (frequencies (take 100 (map #(sig % (take 6 planes)) points)))) ;-> 21

(count (frequencies (take 1000 (map #(sig % (take 1 planes)) points)))) ;-> 2
(count (frequencies (take 1000 (map #(sig % (take 2 planes)) points)))) ;-> 4
(count (frequencies (take 1000 (map #(sig % (take 3 planes)) points)))) ;-> 8
(count (frequencies (take 1000 (map #(sig % (take 4 planes)) points)))) ;-> 15
(count (frequencies (take 1000 (map #(sig % (take 5 planes)) points)))) ;-> 26
(count (frequencies (take 1000 (map #(sig % (take 6 planes)) points)))) ;-> 38

(count (frequencies (take 10000 (map #(sig % (take 1 planes)) points)))) ; 2 
(count (frequencies (take 10000 (map #(sig % (take 2 planes)) points)))) ; 4 
(count (frequencies (take 10000 (map #(sig % (take 3 planes)) points)))) ; 8 
(count (frequencies (take 10000 (map #(sig % (take 4 planes)) points)))) ; 15 
(count (frequencies (take 10000 (map #(sig % (take 5 planes)) points)))) ; 26 
(count (frequencies (take 10000 (map #(sig % (take 6 planes)) points)))) ; 40 

(count (frequencies (take 100000 (map #(sig % (take 1 planes)) points)))) ; 2 
(count (frequencies (take 100000 (map #(sig % (take 2 planes)) points)))) ; 4 
(count (frequencies (take 100000 (map #(sig % (take 3 planes)) points)))) ; 8 
(count (frequencies (take 100000 (map #(sig % (take 4 planes)) points)))) ; 15
(count (frequencies (take 100000 (map #(sig % (take 5 planes)) points)))) ; 26
(count (frequencies (take 100000 (map #(sig % (take 6 planes)) points)))) ; 41

(count (frequencies (take 1000000 (map #(sig % (take 1 planes)) points)))) ; 2
(count (frequencies (take 1000000 (map #(sig % (take 2 planes)) points)))) ; 4
(count (frequencies (take 1000000 (map #(sig % (take 3 planes)) points)))) ; 8
(count (frequencies (take 1000000 (map #(sig % (take 4 planes)) points)))) ; 15
(count (frequencies (take 1000000 (map #(sig % (take 5 planes)) points)))) ; 26
(count (frequencies (take 1000000 (map #(sig % (take 6 planes)) points)))) ; 42 

;; I'm painfully conscious of having stopped the experiment at the
;; exact point where I got the answer I expected. But my poor little
;; computer is not going to be up to running this for 10000000 points.

;; But this can't just be coincidence, surely?

Fizz Buzz : An Interview Question

;; Fizz Buzz

;; My sources
;; http://blog.codinghorror.com/why-cant-programmers-program/
;; inform me that:

;; The majority of computer science graduates can't solve this problem:

;; Write a program that prints the numbers from 1 to 100.  But for
;; multiples of three print "Fizz" instead of the number and for the
;; multiples of five print "Buzz". For numbers which are multiples of
;; both three and five print "FizzBuzz".

;; And Brother Downing of this parish, who actually hires people to
;; program in Java and Clojure learns me that he does indeed use this
;; to screen job applicants, and that most of them can't do it.

;; It is hard to read a thing like that without thinking: 'hang on, is that harder than it looks?'

;; So I did it, just to check:

;; I decided to use pull it out your ass driven development, where
;; you just pull the answer out of your ass.

;; First bit, print out the numbers from 1 to 100
(range 100) ;-> (0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 ...)

;; Bugger
(range 1 101) ;-> (1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 ...)

;; paranoid now
(last (range 1 101)) ;-> 100

;; Although I guess really (print (range 1 101)) is what I want
;; here. In PIOYADD, we defer this important user interface question for later.

;; Next, print 'Fizz' instead of all the multiples of three
(map (= 0 #(% quot 3)) (range 1 101))
;; ClassCastException java.lang.Boolean cannot be cast to clojure.lang.IFn  clojure.core/map/fn--4207 (core.clj:2485)

;; bugger
(map #(= 0 (quot % 3)) (range 1 101)) ;-> (true true false false false false false false false false false false false false false false false false false false false false false false false false false ...)

;; ok
(map #(if (= 0 (quot % 3)) "Fizz" %) (range 1 101)) ;-> ("Fizz" "Fizz" 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 ...)

;; oh for fuck's sake
(map #(if (= 0 (mod % 3)) "Fizz" %) (range 1 101)) ;-> (1 2 "Fizz" 4 5 "Fizz" 7 8 "Fizz" 10 11 "Fizz" 13 14 "Fizz" 16 17 "Fizz" 19 20 "Fizz" 22 23 "Fizz" 25 26 "Fizz" ...)

;; payday
(map #(case (= 0 (mod % 3)) "Fizz" (= 0 (mod % 5)) "Buzz" %) (range 1 101))
;; IllegalArgumentException No matching clause: false  user/eval1234/fn--1235 (NO_SOURCE_FILE:1)

;; ok, I always screw that up
(map #(cond (= 0 (mod % 3)) "Fizz" (= 0 (mod % 5)) "Buzz" %) (range 1 101))
;; CompilerException java.lang.IllegalArgumentException: cond requires an even number of forms, compiling:(NO_SOURCE_PATH:1:7) 

;; twice usually
(map #(cond (= 0 (mod % 3)) "Fizz" (= 0 (mod % 5)) "Buzz" :else %) (range 1 101)) 
;-> (1 2 "Fizz" 4 "Buzz" "Fizz" 7 8 "Fizz" "Buzz" 11 "Fizz" 13 14 "Fizz" 16 17 "Fizz" 19 "Buzz" "Fizz" 22 23 "Fizz" "Buzz" 26 "Fizz" ...)

;; Bwahhahhhahhh! No, Mr Bond, I expect you to die.

;; And now, for the tricky bit of the problem, we unsheathe the
;; superweapon, copy-and-paste-driven development
(map #(cond (or (= 0 (mod % 3)) (= 0 (mod % 3))) "FizzBuzz" (= 0 (mod % 5)) "Buzz" :else %) (range 1 101)) 
;-> (1 2 "FizzBuzz" 4 "Buzz" "FizzBuzz" 7 8 "FizzBuzz" "Buzz" 11 "FizzBuzz" 13 14 "FizzBuzz" 16 17 "FizzBuzz" 19 "Buzz" "FizzBuzz" 22 23 "FizzBuzz" "Buzz" 26 "FizzBuzz" ...)

;; hmmm
(map #(cond (and (= 0 (mod % 3)) (= 0 (mod % 3))) "FizzBuzz" (= 0 (mod % 5)) "Buzz" :else %) (range 1 101)) 
;-> (1 2 "FizzBuzz" 4 "Buzz" "FizzBuzz" 7 8 "FizzBuzz" "Buzz" 11 "FizzBuzz" 13 14 "FizzBuzz" 16 17 "FizzBuzz" 19 "Buzz" "FizzBuzz" 22 23 "FizzBuzz" "Buzz" 26 "FizzBuzz" ...)

;; still wrong
(map #(cond (and (= 0 (mod % 3)) (= 0 (mod % 5))) "FizzBuzz" (= 0 (mod % 5)) "Buzz" :else %) (range 1 101)) 
;-> (1 2 3 4 "Buzz" 6 7 8 9 "Buzz" 11 12 13 14 "FizzBuzz" 16 17 18 19 "Buzz" 21 22 23 24 "Buzz" 26 27 ...)

;; aargh, where did the fizzes go?
(map #(cond (and (= 0 (mod % 3)) (= 0 (mod % 5))) "FizzBuzz" (= 0 (mod % 3)) "Fizz" (= 0 (mod % 5)) "Buzz" :else %) (range 1 101)) 
;-> (1 2 "Fizz" 4 "Buzz" "Fizz" 7 8 "Fizz" "Buzz" 11 "Fizz" 13 14 "FizzBuzz" 16 17 "Fizz" 19 "Buzz" "Fizz" 22 23 "Fizz" "Buzz" 26 "Fizz" ...)

;; Wahey! That looks done. Three minutes.

;; And now, pretending you're not a filthy hacker driven development:
(map 
 #(cond (and (= 0 (mod % 3)) (= 0 (mod % 5))) "FizzBuzz" 
        (= 0 (mod % 3)) "Fizz" 
        (= 0 (mod % 5)) "Buzz" 
        :else %) (range 1 101)) 
;-> (1 2 "Fizz" 4 "Buzz" "Fizz" 7 8 "Fizz" "Buzz" 11 "Fizz" 13 14 "FizzBuzz" 16 17 "Fizz" 19 "Buzz" "Fizz" 22 23 "Fizz" "Buzz" 26 "Fizz" ...)

(defn divides? [n m] (= 0 (mod n m)))
(divides? 3 15) ;-> false

;; sigh
(defn divides? [n m] (= 0 (mod m n)))
(divides? 3 15) ;-> true
(divides? 15 3) ;-> false
;; rah!

;; and so, behold: beauty is truth, and truth, beauty
(map 
 #(cond (divides? 15 %) "FizzBuzz" 
        (divides? 3 %)  "Fizz" 
        (divides? 5 %)  "Buzz" 
        :else %) (range 1 101)) 
;-> (1 2 "Fizz" 4 "Buzz" "Fizz" 7 8 "Fizz" "Buzz" 11 "Fizz" 13 14 "FizzBuzz" 16 17 "Fizz" 19 "Buzz" "Fizz" 22 23 "Fizz" "Buzz" 26 "Fizz" ...)

;; finally, smugness driven development:
(def fizzbuzz
  (map 
   #(cond (divides? 15 %) "FizzBuzz" 
          (divides? 3 %)  "Fizz" 
          (divides? 5 %)  "Buzz" 
          :else %) (map inc (range))))

;; There are those who would call this 'premature abstraction', but they deserve not the names of men.
(print (take 100 fizzbuzz))
;; (1 2 Fizz 4 Buzz Fizz 7 8 Fizz Buzz 11 Fizz 13 14 FizzBuzz 16 17 Fizz 19 Buzz Fizz 22 23 Fizz Buzz 26 Fizz ...)

;; And I suppose a regression test would be nice, if I'm trying to
;; give some sort of professional impression:
(= 
 (take 27 fizzbuzz)
 (list 1 2 "Fizz" 4 "Buzz" "Fizz" 7 8 "Fizz" "Buzz" 
       11 "Fizz" 13 14 "FizzBuzz" 16 17 "Fizz" 19 "Buzz"
       "Fizz" 22 23 "Fizz" "Buzz" 26 "Fizz"))

;; And some paranoid checks, if I'm going to put this travesty on my blog:
(count (filter #(= "Fizz" %) (take 1000 fizzbuzz))) ; -> 267
(count (filter #(= "FizzBuzz" %) (take 1000 fizzbuzz))) ;-> 66
(count (filter #(= "Buzz" %) (take 1000 fizzbuzz))) ;-> 134

(* (+ 134 66) 5) ;-> 1000
(* (+ 267 66) 3) ;-> 999

;; bah, that's close enough for government work. I declare myself
;; done. Three minutes of flail and two minutes of tidying up and
;; checking it works.

;; So my question to the wider community is: Does my three minutes of
;; flailing trying to remember the semantics of my favourite language
;; (which I can perfectly imagine looking dreadful at an interview)
;; count as a fail or a pass?

;; In a previously unknown assembly language, but given a library to
;; print numbers on the screen, I can imagine taking half an hour to
;; get this working. Without the library, it's a research project.

;; And obviously, if I tried to do it in Haskell, I'd have to spend
;; three weeks remembering what a monad was in order to print the
;; damned thing out, but at least the type system would magically
;; guarantee the correctness of the final program.

;; In other words, is fizzbuzz really actually quite hard, or is
;; everyone out there a complete idiot?

Planes in Space : Checking by Hand

;; Ok, so I reckon:

(defn regions [n m]
  (cond (= n 1) (list m (inc m))
        (= m 0) (concat (repeat n 0) '(1))
        :else (map + 
                   (regions n (dec m))
                   (concat '(0) (regions (dec n) (dec m))) 
                   (concat (regions (dec n) (dec m)) '(0)))))

;; In fact, I reckon that I can prove it. 

;; It's quite easy to prove something that isn't true. Usually when
;; you run across a counterexample, that shows you where the
;; unsuspected false step in your reasoning was.

;; So I'd like to verify the formula on lots of examples.

;; But how?

;; These I can count in my head.

(regions 3 0) ;-> (0 0 0 1) ; just space, no planes
(regions 3 1) ;-> (0 0 1 2) ; a plane cuts space in half
(regions 3 2) ;-> (0 1 4 4) ; two planes intersect in one line, cutting space into four
(regions 3 3) ;-> (1 6 12 8) ; 1 point where three planes meet, 
                             ; six coordinate half-axes, 
                             ; three coordinate planes divided into four quadrants each,
                             ; and eight octants

;; I am sort of confident that four planes make fifteen regions and intersect in four points.
(regions 3 4) ;-> (4 18 28 15)
;; But every attempt to count the 18 lines and 28 regions I make ends up relying on 
;; the sort of arguments I made to make the recursion in the first place.

;; And at this point my intuition breaks.
(regions 3 5) ;-> (10 40 55 26)

;; I mean, five planes, any three intersect in a point,
;; 10 ways to choose 3 planes from five, so 10 points,
;; but after that I'm dead.

;; And this? Six choose 3 is 20, I can see that....
(regions 3 6) ;-> (20 75 96 42)

;; And I defy anyone to even picture
(regions 5 8) ;-> (56 350 896 1176 792 219)

;; 8 choose 5 is (/ (* 8 7 6) (* 1 2 3)) = 56, which is fair enough.

;; 8 choose 4 is (/ (* 8 7 6 5) (* 1 2 3 4)) is 70

;; If we take any 4 hyperplanes from our 8, they'll define a line.
;; The four remaining hyperplanes then divide each line like:
(regions 1 4) ;-> (4 5)
;; So if each of those lines is sliced into 5 pieces then that's our 350.
;; But I'm just using the same recursive argument again.

;; So I don't even know if that's evidence or not.

;; One nice thing about it, the formula has an alternating sum property, 
;; like the Euler index.

(defn signature [lst]
  (reduce + 
          (map * 
               (apply concat (repeat '(+1 -1))) 
               (reverse lst))))

(signature (regions 5 8)) ;-> 1

(def regions (memoize regions))
(regions 7 23) ;-> (245157 1817046 5787628 10271800 10973116 7057688 2531288 390656)
(signature (regions 7 23)) ;-> 1

(regions 23 20) ;-> (0 0 0 1 40 760 9120 77520 496128 2480640 9922560 32248320 85995520 189190144 343982080 515973120 635043840 635043840 508035072 317521920 149422080 49807360 10485760 1048576)
(signature (regions 23 20)) ;-> 1

;; But annoyingly, you can just read that property straight from the recursion!