Rerum Cognoscere Causas

I don't actually know the answer to this question, or even whether it can be answered. Yet.

;; Rerum Cognoscere Causas

;; A question of epistemology, with reference to the detection of causality, incorporating a slightly smug trick for generating the same random things over and over again.

;; In Al'no'oth, people are rated by their strength (STR) and intelligence (INT)

;; Both are easily determined with simple tests. For instance, a
;; person of STR 11 will win 45% of arm-wrestles against a person of
;; STR 12.

;; The philosophers of Al'no'oth have long suspected that the Gods
;; determine STR and INT by rolling six-sided dice:

(def ^:dynamic *randomizer* (java.util.Random. 0))

(defn rand-int [n]
  (.nextInt *randomizer* n))

(defn D6 [] (inc (rand-int 6)))

;; In fact, the traditional belief is that the Gods use three dice:

(defn three-D6 []
  (reduce + (for [i (range 3)] (D6))))

;; But although all unenhanced humans have scores between 3 and 18 in
;; both characteristics, affirming the suspicions of the philosophers,
;; the higher scores are slightly more common than they should be were
;; 3 dice simply added.

;; Another method the Gods might use would be to roll four dice, and
;; to discard the lowest of them:

(defn three-from-four-D6 []
  (reduce + (drop 1 (sort (for [i (range 4)] (D6))))))


;; At this point, the philosophers disagree. 

;; Those of the first e'dition maintain that all scores are generated using 3 dice.

(defn first-edition [] {:str (three-D6) :int (three-D6)})

;; Those of the second e'dition maintain that, occasionally, by caprice, the Gods use the modified procedure.

(defn mixed []
  (if (zero? (rand-int 10))
    (three-from-four-D6)
    (three-D6)))

(defn second-edition [] {:str (mixed) :int (mixed)})

;; Those of the third e'dition believe that one person in ten is a
;; "player-character".  These characters are awe'some.

;; Awe'someness is not directly detectable, but a player-character has
;; all his characteristics generated using the modified procedure.

(defn third-edition []
  (if (zero? (rand-int 10))
    {:str (three-from-four-D6) :int (three-from-four-D6)}
    {:str (three-D6) :int (three-D6)}))


;; With huge effort, and a very large number of arm-wrestles and
;; riddling contests, a great sage has tabulated the characteristics
;; of an entire village:

(def village
  (binding [*randomizer* (java.util.Random. 0)]
    (doall (repeatedly 100 (case (rand-int 3)
                    0 first-edition 
                    1 second-edition
                    2 third-edition)))))


village ;-> ({:str 13, :int 18} {:str 11, :int 18} {:str 14, :int 15} {:str 6, :int 12} {:str 14, :int 13} {:str 18, :int 10} {:str 15, :int 11} {:str 12, :int 15} {:str 7, :int 8} {:str 16, :int 12} {:str 8, :int 7} {:str 9, :int 14} {:str 10, :int 9} {:str 11, :int 10} {:str 5, :int 10} {:str 7, :int 9} {:str 9, :int 13} {:str 12, :int 9} {:str 13, :int 9} {:str 5, :int 9} {:str 8, :int 13} {:str 9, :int 11} {:str 13, :int 14} {:str 12, :int 14} {:str 12, :int 17} {:str 14, :int 9} {:str 10, :int 11} {:str 18, :int 17} {:str 11, :int 9} {:str 8, :int 9} {:str 15, :int 13} {:str 8, :int 5} {:str 11, :int 9} {:str 10, :int 8} {:str 9, :int 12} {:str 5, :int 11} {:str 10, :int 7} {:str 9, :int 14} {:str 11, :int 9} {:str 11, :int 12} {:str 12, :int 13} {:str 15, :int 9} {:str 12, :int 12} {:str 6, :int 13} {:str 5, :int 4} {:str 12, :int 13} {:str 15, :int 10} {:str 14, :int 14} {:str 11, :int 4} {:str 12, :int 9} {:str 10, :int 12} {:str 7, :int 12} {:str 8, :int 11} {:str 10, :int 10} {:str 9, :int 8} {:str 8, :int 12} {:str 7, :int 9} {:str 13, :int 3} {:str 14, :int 9} {:str 8, :int 9} {:str 10, :int 11} {:str 15, :int 4} {:str 10, :int 11} {:str 8, :int 10} {:str 15, :int 10} {:str 8, :int 13} {:str 12, :int 5} {:str 8, :int 16} {:str 4, :int 8} {:str 10, :int 18} {:str 12, :int 12} {:str 11, :int 10} {:str 12, :int 8} {:str 12, :int 13} {:str 8, :int 12} {:str 9, :int 12} {:str 12, :int 10} {:str 15, :int 10} {:str 8, :int 11} {:str 7, :int 11} {:str 4, :int 8} {:str 12, :int 11} {:str 13, :int 9} {:str 14, :int 13} {:str 5, :int 9} {:str 17, :int 10} {:str 8, :int 13} {:str 9, :int 10} {:str 5, :int 14} {:str 15, :int 12} {:str 13, :int 13} {:str 11, :int 8} {:str 8, :int 6} {:str 12, :int 8} {:str 10, :int 3} {:str 14, :int 9} {:str 15, :int 12} {:str 15, :int 14} {:str 6, :int 10} {:str 16, :int 13})

;; How should the philosophers reason amongst themselves, should their object be to determine the System of their World?

Clojure's Reader is Unsafe

;; Clojure's Reader is Unsafe

;; dependencies:
;; [org.clojure/clojure "1.4.0"]
;; [org.clojure/tools.reader "0.7.0"]


;; By a weird coincidence this week I wanted to read a file of data
;; coming from a web app, and I was about to use Clojure's reader to
;; do it. I remembered that there was some variable that I needed to
;; bind to stop it executing arbitrary code, so I went googling.

;; And found that there's been a fair bit of debate on just this topic recently.

;; Essentially, Clojure's reader is the thing that turns strings into
;; data structures when you're reading in a program or typing at the
;; REPL. And you have it available programmatically.

(read-string "(+ 1 2)") ; (+ 1 2)

;; Not so impressive, you say?

;; Consider:

(first "(+ 1 2)") ; \(
(second "(+ 1 2)") ; \+
(count "(+ 1 2)") ; 7
(type "(+ 1 2)") ; java.lang.String
(type (first "(+ 1 2)")) ; java.lang.Character

;; "(+ 1 2)" is a String, a flat data structure of seven Characters

(first (read-string "(+ 1 2)")) ; + 
(second (read-string "(+ 1 2)")) ; 1
(count (read-string "(+ 1 2)")) ; 3 
(type  (read-string "(+ 1 2)")) ; clojure.lang.PersistentList
(type (first (read-string "(+ 1 2)"))) ; clojure.lang.Symbol

;; Clearly some sort of sea-change has occurred.

;; Even better, the value of a string is a similar string
(eval "(+ 1 2)") ; "(+ 1 2)"
(identical? (eval "(+ 1 2)") "(+ 1 2)") ; true
;; But the thing that read-string makes is a program which can run
(eval (read-string "(+ 1 2)")) ; 3
;; Another way to make such a program is:
(eval (list (quote +) (quote 1) (quote 2))) ; 3

;; In my innocence, I was going to use this magic to read some data out of a file that had been produced
;; by a web server.

(def file "{:username leethaxor :score 5 :alignment :chaotic-evil}")

(:username (read-string file)) ;-> leethaxor


;; IT TURNS OUT THAT YOU SHOULD NEVER DO THIS, EVEN THOUGH THE ABILITY
;; TO DO THIS IS ONE OF THE THINGS THAT IS NICE ABOUT LISPS.

;; There are various reasons why:

;; Firstly, note that programs running on a computer can affect that computer

;; If you happen to have a file called precious.txt in the directory
;; where clojure is running, then you are strongly advised not to
;; execute any of the following code. Especially if it contains your
;; first novel and you haven't got round to backing up yet this year.

(clojure.java.shell/sh "touch" "precious.txt") ;-> {:exit 0, :out "", :err ""}
(clojure.java.shell/sh "cat" "precious.txt") ;-> {:exit 0, :out "", :err ""}


;; Secondly, note that there is special syntax to cause things to happen while
;; reading: 

(read-string "    #=(clojure.java.shell/sh \"rm\" \"precious.txt\")    ") ;-> {:exit 0, :out "", :err ""}

;; That did not get turned into a data structure. That got executed. precious.txt is gone.

(clojure.java.shell/sh "cat" "precious.txt") ;-> {:exit 1, :out "", :err "cat: precious.txt: No such file or directory\n"}

;; Luckily we can rebuild it:
(clojure.java.shell/sh "touch" "precious.txt") ;-> {:exit 0, :out "", :err ""}

(clojure.java.shell/sh "cat" "precious.txt") ;-> {:exit 0, :out "", :err ""}

;; And thirdly note that not everyone on the internet is on your side:
(def file-of-evil " {:username #=(clojure.java.shell/sh \"rm\" \"precious.txt\") :score 5 :alignment :chaotic-evil}")

(read-string file-of-evil) ;-> {:username {:exit 0, :out "", :err ""}, :score 5, :alignment :chaotic-evil}

;; This is a thing that you're supposed to guard against when reading data structures.

(clojure.java.shell/sh "touch" "precious.txt") ;-> {:exit 0, :out "", :err ""}

;; It is traditional to do that like this:

(binding [*read-eval* false] (read-string "(+ 1 2)")) ;-> (+ 1 2)
(binding [*read-eval* false] (read-string file-of-evil))
; -> RuntimeException EvalReader not allowed when *read-eval* is false.  clojure.lang.Util.runtimeException (Util.java:170)

;; And in my innocence, I was about to actually do this, which would have been a mistake:

(def string-of-evil "#java.io.FileWriter[\"precious.txt\"]")

(spit "precious.txt" "precious content") ;-> nil
(slurp "precious.txt") ;-> "precious content"
(binding [*read-eval* false] (clojure.core/read-string string-of-evil)) ;-> #<FileWriter java.io.FileWriter@17b2712>
(slurp "precious.txt") ;-> ""

;; bugger

;; It turns out that even with *read-eval* bound to false, the Clojure
;; reader can be persuaded to execute arbitrary java constructors, and
;; some of those have side-effects.

;; I believe that that particular loophole is going to be fixed in the
;; next version of clojure, but the general thinking seems to be that
;; the reader is just intrinsically an unsafe thing and that you
;; shouldn't use it to read strings which you aren't sure are friendly.

;; Luckily, there is a library function that does what I thought the
;; reader would do, and that is apparently Absolutely Guaranteed by
;; Act of Parliament not to do bad things:

(require 'clojure.tools.reader.edn)

(spit "precious.txt" "precious content") ;-> nil
(slurp "precious.txt") ;-> "precious content"
(clojure.tools.reader.edn/read-string string-of-evil)
;; -> ExceptionInfo No reader function for tag java.io.FileWriter  clojure.core/ex-info (core.clj:4227)
(clojure.tools.reader.edn/read-string file-of-evil)
;; -> ExceptionInfo No reader function for tag =  clojure.core/ex-info (core.clj:4227)

;; But which is still usable for the purposes of good:
(clojure.tools.reader.edn/read-string file) ;-> {:username leethaxor, :score 5, :alignment :chaotic-evil}
(clojure.tools.reader.edn/read-string "(+ 1 2)") ;-> (+ 1 2)
(eval (clojure.tools.reader.edn/read-string "(+ 1 2)")) ; 3

(slurp "precious.txt") ;-> "precious content"

;; So that's Super News!

;; We've got a function, read-string, which can turn strings into data structures while doing no harm whatsoever,
;; and another function, read-string, which can turn strings into data structures while potentially causing arbitrary side effects depending on the content of the strings.

;; Clearly nothing could possibly go wrong with this arrangement,
;; especially since read-string will be fixed in clojure 1.5 so that
;; binding *read-eval* to false around it will make it safe, as far as
;; anyone knows, even though the advice is still never to use it to
;; read untrusted strings, for some reason.

Linear Congruential Random Number Generators

;; Linear Congruential Random Number Generators

;; I was wondering how pseudo-random numbers are actually generated.

;; Apparently it's still done as I remember it done on the ZX80, using multiply, add and chop.

;; There was a nice article on Rosetta Code :
;; http://rosettacode.org/wiki/Linear_congruential_generator

;; But it didn't have a clojure example, so I made an account and added one:

;; This code took me a couple of minutes to write, and replicates a couple of classic generators

(defn iterator [a b c]
  (fn[x] (mod (+ (* a x) b) (bit-shift-left 1 c))))

(def bsd (drop 1 (iterate (iterator 1103515245 12345 31) 0)))

(def ms (drop 1 (for [x (iterate  (iterator 214013 2531011 31) 0)] (bit-shift-right x 16))))

(take 10 bsd) ;-> (12345 1406932606 654583775 1449466924 229283573 1109335178 1051550459 1293799192 794471793 551188310)
(take 10 ms) ;-> (38 7719 21238 2437 8855 11797 8365 32285 10450 30612)

;; What's really nice is that the natural expression of this idea in
;; Clojure is (I think) one of the neatest pieces of code on there.

;; Java uses a similar generator, but it takes some proving since it won't let you see the internal state directly:

(defn java [seed] 
  (drop 1 (iterate (iterator 0x5DEECE66D 0xB 48) seed )))

(def java-random-ints (map #(unchecked-int (bit-shift-right (long %) 16)) (java 0x5DEECE66DN)))

(take 10 java-random-ints) ;-> (-1155484576 -723955400 1033096058 -1690734402 -1557280266 1327362106 -1930858313 502539523 -1728529858 -938301587)

;; We can also get them direct from java:

(def java-util-random
  (let [r (java.util.Random. 0)]
    (for [i (range)] (.nextInt r))))

(take 10 java-util-random) ;-> (-1155484576 -723955400 1033096058 -1690734402 -1557280266 1327362106 -1930858313 502539523 -1728529858 -938301587)


;; And compare a few:

(drop-while zero? (take 100000 (map - java-util-random java-random-ints)))  ;-> ()


;; Incidentally while messing around here I discovered how to call a protected method of a java class from clojure:

(.nextInt (java.util.Random. 0) ) ;-> -1155484576
(.next (java.util.Random. 0) 32)  ; java.lang.IllegalArgumentException: No matching method found: next for class java.util.Random

;; It's lying.
;; But I know it's in there somewhere:

(defn myrandom [seed] (proxy [java.util.Random][seed]
                    (next [a] (proxy-super next a))))

(def r (myrandom 0)) ;-> #'user/r
(.nextInt (myrandom 0)) ;-> -1155484576
(.next (myrandom 0) 32) ;-> -1155484576

;; Stupid language.

Packing and Unpacking Numbers : Deterministic Randomness

;; Packing and Unpacking Numbers

;; This post is kind of pointless, since while writing it I discovered a better way to do the same thing.

;; But I also think it's kind of cool, so I'm putting it up in case anyone else finds a use for it.

;; Suppose we've got a function which gives random results (I mean, deliberately...)

(defn randomfunction []
  (list 
   (case (rand-int 2) 0 '+ 1 '-)
   (rand-int 100)
   (list 
    (case (rand-int 2) 0 '+ 1 '-)
    (rand-int 100)
    (rand-int 100))))

(randomfunction) ;-> (+ 43 (- 71 75))
(randomfunction) ;-> (- 18 (- 13 22))
(randomfunction) ;-> (- 36 (+ 60 3))
(randomfunction) ;-> (+ 56 (+ 34 4))
(randomfunction) ;-> (- 87 (- 61 76))

;; And it comes to testing it:

(= (randomfunction) '(+ 26 (- 69 0))) ;-> false

;; We'd prefer it to be deterministic. 

;; So we can separate the concerns thus:

(defn deterministicfunction [a b c d e]
  (list 
   (case a 0 '+ 1 '-)
   b
   (list 
    (case c 0 '+ 1 '-)
    d
    e)))

(defn randomfunction []
  (deterministicfunction (rand-int 2) (rand-int 100) (rand-int 2) (rand-int 100) (rand-int 100)))

;; tadaa!

(= (deterministicfunction 0 26 1 69 0) '(+ 26 (- 69 0))) ;-> true

(randomfunction) ;-> (+ 30 (+ 52 95))
(randomfunction) ;-> (- 97 (- 57 61))
(randomfunction) ;-> (+ 88 (+ 26 15))

;; Although useful, this is a right royal pain in the neck, and horribly error-prone

;; It occurred to me that it might be made less painful by setting the
;; seed of the random number generator before running tests, but I
;; also had a need to make 'deterministic random' values in the code
;; itself, and having everything hanging off a global variable that
;; any old function occasionally resets struck me as a recipe for
;; disaster.

;; So I thought, well why not pack all the random parameters into one big number:

(def << bit-shift-left)

;; If all the random ranges are powers of two then this is easy:

(let [[a b c d] [(rand-int 8) (rand-int 4) (rand-int 8) (rand-int 2)]]
 [[a b c d] (+ (<< a (+ 1 3 2)) (<< b (+ 1 3)) (<< c 1) (<< d 0))]) 
; [[2 2 4 0] 168] 
; [[3 1 4 0] 216] 
; [[5 3 6 0] 380]

;; Looking at the last one
;; 5/3/6/0 is 
;; 101/11/110/0 in binary

;; and binary 101111100 is:
2r101111100 ; 380

;; in decimal

;; To extract the numbers, we can do:

(mod 380 2) ; 0 
(/ (- 380 (mod 380 2)) 2) ; 190 
(mod 190 8) ; 6 
(/ (- 190 (mod 190 8)) 8) ; 23
(mod 23 4) ; 3
(/ (- 23 (mod 23 4)) 4) ; 5
(mod 5 8) ; 5

;; And this motivates the following recursion:

(defn ^:dynamic unpack [rl n]
  (if (empty? rl) '() 
      (let [[rh & rr] rl]
        (cons (mod n rh) (unpack rr (/ (- n (mod n rh)) rh))))))

;; Which behaves as I was hoping

(unpack '() 0) ;-> ()
(unpack '() 1) ;-> ()
(unpack '(8) 1) ;-> (1)
(unpack '(8) 9) ;-> (1)
(unpack '(8 8) 9) ;-> (1 1)

(unpack '(2 8 4 8) 380) ;-> (0 6 3 5)
(unpack '(2 8 4 8) 216) ;-> (0 4 1 3)
(unpack '(2 8 4 8) 168) ;-> (0 4 2 2)

;; And that then motivates the following packing function

(defn ^:dynamic pack [rlist vlist]
  (if (empty? rlist) 0
      (let [[rh & rr] rlist
            [vh & vr] vlist] 
        (+ (* rh (pack rr vr)) vh))))

;; Which also works
(pack '()()) ;-> 0
(pack '(8) '(1)) ;-> 1
(pack '(8 8) '(1 1)) ;-> 9
(pack '(2 8 4 8) '(0 6 3 5)) ;-> 380
(pack '(2 8 4 8) '(0 4 1 3)) ;-> 216
(pack '(2 8 4 8) '(0 4 2 2)) ;-> 168

;; But the cool thing is that it works whether or not the ranges are powers of two:

(use 'clojure.tools.trace)

(let [a (for [i (range (inc (rand-int 7)))] (rand-int 30))
      b (for [i a] (rand-int i))]
  (clojure.tools.trace/dotrace [unpack pack] (unpack a (pack a b))))

;; TRACE t3191: (pack (12 23 5 5 8 3) (4 9 3 4 0 0))
;; TRACE t3192: | (pack (23 5 5 8 3) (9 3 4 0 0))
;; TRACE t3193: | | (pack (5 5 8 3) (3 4 0 0))
;; TRACE t3194: | | | (pack (5 8 3) (4 0 0))
;; TRACE t3195: | | | | (pack (8 3) (0 0))
;; TRACE t3196: | | | | | (pack (3) (0))
;; TRACE t3197: | | | | | | (pack nil nil)
;; TRACE t3197: | | | | | | => 0
;; TRACE t3196: | | | | | => 0
;; TRACE t3195: | | | | => 0
;; TRACE t3194: | | | => 4
;; TRACE t3193: | | => 23
;; TRACE t3192: | => 538
;; TRACE t3191: => 6460
;; TRACE t3198: (unpack (12 23 5 5 8 3) 6460)
;; TRACE t3199: | (unpack (23 5 5 8 3) 538)
;; TRACE t3200: | | (unpack (5 5 8 3) 23)
;; TRACE t3201: | | | (unpack (5 8 3) 4)
;; TRACE t3202: | | | | (unpack (8 3) 0)
;; TRACE t3203: | | | | | (unpack (3) 0)
;; TRACE t3204: | | | | | | (unpack nil 0)
;; TRACE t3204: | | | | | | => ()
;; TRACE t3203: | | | | | => (0)
;; TRACE t3202: | | | | => (0 0)
;; TRACE t3201: | | | => (4 0 0)
;; TRACE t3200: | | => (3 4 0 0)
;; TRACE t3199: | => (9 3 4 0 0)
;; TRACE t3198: => (4 9 3 4 0 0)
;; (4 9 3 4 0 0)

;; Which I found slightly smug-making

;; So now we can say:
(defn deterministicfunction [big-rand]
  (let [[a b c d e] (unpack [2 100 1 100 100] big-rand)]
    (list 
     (case a 0 '+ 1 '-)
     b
     (list 
      (case c 0 '+ 1 '-)
      d
      e))))

(defn randomfunction []
  (deterministicfunction (rand-int (* 2 100 2 100 100))))

;; Behold the many benefits of this scheme:

;; We can enumerate all possible varieties
(deterministicfunction 1) ;-> (- 0 (+ 0 0))
(deterministicfunction 2) ;-> (+ 1 (+ 0 0))
(deterministicfunction 3) ;-> (- 1 (+ 0 0))

;; We can generate examples using a single number
(deterministicfunction 10345778979) ;-> (- 89 (+ 94 88))

;; Or we can get random examples
(randomfunction) ;-> (+ 99 (+ 72 8))

;; If we are not desperately fussed about the statistical properties of things
;; we do not need to worry about the exact range of the random number at all:
;; I think that this will produce only tiny biases:

(defn randomfunction []
  (deterministicfunction (rand-int (<< 1 30))))

;; And it still looks pretty random to me
(randomfunction) ;-> (+ 76 (+ 10 87))
(randomfunction) ;-> (+ 76 (+ 81 3))
(randomfunction) ;-> (- 27 (+ 18 63))

Clojure Emacs Eval Paste Keyboard Macros / Generating your Regression Tests Automatically while Writing your Functions

;; Clojure Emacs Eval Paste Keyboard Macros / 
;; Generating your Regression Tests Automatically while Writing your Functions

;; There's a phenomenally useful feature of emacs when using clojure
;; and nrepl where you can evaluate the result of an
;; expression and paste the result directly into a buffer.

;; However my attempts to use it in keyboard macros have always led to frustration.

;; I asked about it on stack overflow:
;; http://stackoverflow.com/questions/14959155/how-can-i-make-emacs-keyboard-macros-work-properly-when-pasting-the-results-of-c#14960321

;; and a nice chap called sds helped me come up with the following work-around:

;; Say I've defined the function:

(defn << [a b]
  (if (zero? b) a
      (<< (* 2 a) (dec b))))


;; And now I want to check that it works:

;; I evaluate the following two forms in emacs lisp (in the *scratch* buffer)

(defun clojure-eval-paste-test ()
  (interactive)
  (next-line)
  (move-beginning-of-line nil)
  (insert "(is= ")
  (forward-sexp)
  (insert " ")
  (nrepl-eval-last-expression 't)
  (sleep-for 0.1)
  (insert ")"))

(global-set-key [f5] #'clojure-eval-paste-test)


;; Now I type:

(<< 1 0)
(<< 1 1)
(<< 1 2)
(<< 2 0)
(<< 2 2)
(<< 2 3)

;; And then I move point to just above the first line, and press f5 six times, and I get:


(is= (<< 1 0) 1)
(is= (<< 1 1) 2)
(is= (<< 1 2) 4)
(is= (<< 2 0) 2)
(is= (<< 2 2) 8)
(is= (<< 2 3) 16)

;; Which both reassures me that it does work, and provides the skeleton of a regression test for it!



;; There are more elegant ways of doing this, presumably, one being to
;; do as sds suggests, and code up the function using lower level
;; calls rather than the (sleep-for 0.1) thing that will obviously go
;; wrong if the eval takes too long.

;; And another way would be to fix it so that nrepl-eval-last-expression actually works in keyboard macros. 

;; But this will do for now, and gets rid of a papercut bug that's been annoying me for ages.

Mathematics as an Experimental Science

;; Mathematics as an experimental science

;; Beta-Bernoulli distribution / Machine Learning

;; I'm reading Kevin Murphy's Machine Learning: A Probabilistic Perspective.

;; As one of the exercises in chapter three, I've just proved that:

;; "For the Beta-Bernoulli model, the probability of the data is the
;; ratio of the normalizing constants of the prior distribution and
;; the posterior distribution". 

;; What I think this actually means is:

;; If you start with Coins which are biased, and the bias is chosen
;; uniformly (so that for instance a fair coin is just as likely as a
;; coin which comes up heads 1/3 of the time, or a coin which always
;; comes up heads).

;; And you pick out the fair ones by tossing them all 20 times, and
;; throw away all the ones that don't come up heads exactly 10 times
;; and tails 10 times.

;; Then, because your test has not quite guaranteed you fair coins,
;; but only coins which are quite a bit fairer than you started with:

;; You are ever so slightly more likely to see three heads (or three
;; tails) in three tosses than you would be if the coins were truly
;; fair.

;; In fact the bias is now the same as if the coins had been chosen
;; from the Beta(11,11) distribution.

;; And so the chances of getting three heads in a row is not, as you
;; might have naively expected, one in eight,

;; But rather B(14,11)/B(11,11), where B is the Beta function

;; And that works out to be G(14)G(11)G(22)/G(25)G(11)G(11), where G is the Gamma function.

;; And that is 13!10!21!/24!10!10!, where ! is the factorial function

;; And that is (11*12*13)/(22*23*24), where / and * need no introduction if you have got this far in the post.

;; And that is:

(/ (* 11 12 13) 22 23 24) ;-> 13/92

;; or 

(float (/ 13 92)) ; 0.14130434

;; Which is to say, slightly more than one in eight.

;; And I am inordinately pleased with myself, not for having proved
;; this result, which was easy, but for having worked out what the
;; mysterious squiggles in the book might actually mean in practice.

;; And it occurs to me that Mathematics is in fact an experimental
;; science, which makes definite predictions about physical things,
;; such as the movements of electrons in my computer when I type the
;; following things:

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Let us make coins with biases distributed uniformly between 0 and 1:

(defn make-coin []
  (let [theta (rand)]
    (fn [] (if (< (rand) theta) :h :t))))

;; Here is one:
(def a (make-coin))

;; Let us test it
(frequencies (for [i (range 20)] (a))) ;-> {:h 11, :t 9}

;; It is no good. To the scrapheap with it:
(ns-unmap *ns* 'a)

;; Let us instead make a large number of coins:
(def coins (for [i (range)] (make-coin)))

;; And throw away all the ones that do not satisfy our criterion
(def fair-coins
  (filter (fn[coin] (= {:t 10 :h 10 } (frequencies (for [i (range 20)] (coin))))) coins))

;; Let us toss each of our fair coins three times:
(def results (for [coin fair-coins] (frequencies (for [i (range 3)] (coin)))))

;; And add up all the results
(def collected-results (drop 1 (reductions (fn[m f] (assoc m f (inc (get m f 0)))) {} results)))

;; And calculate the empirical distribution of those results:
(def empirical-distribution (map 
                             (fn[m] [(float (/ (m {:h 3} 0) (reduce + (vals m))))
                                     (float (/ (m {:t 3} 0) (reduce + (vals m))))])
                             collected-results))

;; It does seem to me that after a while, these numbers settle down to something near 14%
(doseq [e empirical-distribution] (println e))

A Naive Bayesian Classifier

Anyway this very evening Hermann Hauser of this parish is telling me that Machine Learning is the coming next big thing and those skills are highly in demand.

So if anybody wants to give me a fuckload of money for this old rope, they must feel free

;; A Naive Bayesian Classifier

;; Ed Jackson ( http://boss-level.com ) and I are currently working
;; our way through Kevin Murphy's book:
;; Machine Learning: A Probabilistic Perspective.

;; Ed tells me that there is big money in this Machine Learning game.

;; We have just read chapter 3, which involves a hair-raising amount
;; of Greek, but when all is said and done boils down to this:

;; We got Martians, who tend to be thin, tall and green
(defn make-martian []
  #{:martian
    (if (< (rand) 0.2) :fat :thin)
    (if (< (rand) 0.7) :tall :short)
    (if (< (rand) 0.8) :green :blue)})

; And we got Venusians, who tend to be fat and blue
(defn make-venusian []
  #{:venusian
    (if (< (rand) 0.8) :fat :thin)
    (if (< (rand) 0.5) :tall :short)
    (if (< (rand) 0.3) :green :blue)})

; They visit Harry's 
(defn make-space-bar-patron []
  (if (< (rand) 0.3) (make-martian) (make-venusian)))

; Harry takes notes
(def client-database
  (for [i (range 10000)] (make-space-bar-patron)))

; Like this:
(take 5 client-database)
(#{:tall :martian :blue :thin} 
 #{:blue :venusian :short :fat} 
 #{:blue :venusian :short :thin} 
 #{:tall :martian :green :fat} 
 #{:martian :blue :short :thin})

; And tabulates them so:
(def martians (filter :martian client-database))
(def venusians (filter :venusian client-database))

(count martians) ; 3054
(count (filter :thin martians)) ; 2454
(count (filter :green martians)) ; 2438
(count (filter :tall martians)) ; 2177


(count venusians) ; 6946
(count (filter :fat venusians)) ; 5628
(count (filter :blue venusians)) ; 4919
(count (filter :tall venusians)) ; 3535


; A guy walks into Harry's
(def new-guy (make-space-bar-patron))

; But he is not wearing his uniform
(reduce disj new-guy [:martian :venusian])  ; #{:blue :short :fat}

; And so a sweepstake comes into being:
(defn probability [characteristic given-class prior-for prior-against]
  (let [class (if given-class (filter given-class client-database) client-database)
        for (count (filter characteristic class))
        against (- (count class) for)
        total-for (+ prior-for for)
        total-against (+ prior-against against)]
    (/ total-for (+ total-for total-against))))



(def m (* (probability :martian nil 1 1)      ; most people come in here are venusians
          (probability :blue  :martian 1 1)   ; and martians tend to be green
          (probability :short :martian 1 1)   ; also he's short. Martians are tall, often
          (probability :fat   :martian 1 1))) ; and fat. All the fat guys is from Venusburg

(def v (* 
        (probability :venusian nil 1 1)       ; most likely he's from venus
        (probability :blue  :venusian 1 1)    ; the blueness bears this out
        (probability :short :venusian 1 1)    ; sometimes they are short, sometimes they are tall
        (probability :fat   :venusian 1 1)))  ; but they are mostly fat

(float (/ m (+ m v))) ; 0.017495114
(float (/ v (+ m v))) ; 0.9825049

;; 45:1 he's from Venus !

new-guy ; #{:blue :short :venusian :fat}

;; As indeed he is.

Clojure is Fast / Is Clojure Still Fast? / Clojure is as Fast as Java

;; Is Clojure Still Fast ?

;; Once upon a time I wrote a blog post saying that clojure was fast.  It still is, and optimizing
;; it is now much easier than it used to be, but it doesn't seem to be *quite* as fast as it once
;; was.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; I optimized a program to solve a differential equation in a
;; simple-minded way:

;; The equation was: dy/dt = f(t,y) where f(t, y) = t-y and y=0 when t=0
;; ( the exact solution is: y=e^(-t)+t-1 )

;; Here's a program to solve it using Euler's method

(defn f [t y] (- t y))

(defn solveit [t0 y0 h its]
  (if (> its 0) 
    (let [t1 (+ t0 h)
          y1 (+ y0 (* h (f t0 y0)))]
      (recur t1 y1 h (dec its)))
    [t0 y0 h its]))

;; For every go round the loop we have to:
;; compare its with 0,
;; branch depending on the result,
;; add t0 to h,
;; call f with t0 and y0,
;; multiply h and the result,
;; add that to y0,
;; jump.

;; So if this was an assembly language program that worked the way
;; you'd expect, each loop would take 7 cycles.

;; According to /proc/cpuinfo, my netbook runs at 1.662 GHz.

(def *cpuspeed* 1.662)

;; We care about how many cycles each iteration of the solver takes:

(defmacro cyclesperit [expr its]
  `(let [start# (. System (nanoTime))
         ret# ( ~@expr (/ 1.0 ~its) ~its )
         finish# (. System (nanoTime))]
     (int (/ (* *cpuspeed* (- finish# start#)) ~its))))


;; With the program as written, this estimate turns out to have been a little optimistic.

;; The figures in the original post were on a desktop machine that I no longer have, which was more
;; powerful than my current netbook.

;; To two significant figures, the results of the timing expression
(cyclesperit (solveit 0.0 1.0) 1000000) 
;; are:

;; On my old desktop with clojure 1.2 : 2300 cycles

;; On my netbook with clojure 1.2 : 2800 cycles
;; On my netbook with clojure 1.3 : 2500 cycles
;; On my netbook with clojure 1.4 : 2400 cycles

;; So it looks like my netbook is not only slower in clock speed than my desktop was, but also in
;; terms of cycles/iteration. That's not surprising as the netbook has an Atom processor, optimized
;; for low power rather than for speed.

;; But it also looks as though clojure has been speeding up slightly, which has almost made up for that.

;; I'm also assuming that the JVM itself hasn't changed much since the original blog post.  The
;; netbook timings were all done today on the same JVM, but the desktop timings are from a while
;; back, so that might account for some differences.

;; So we're looking at a slowdown of about 300 times over what we could probably
;; achieve coding in assembler or in C with a good optimizing compiler (and of
;; course I'm assuming that floating point operations take one cycle each)

;; This is about the sort of speed that you'd expect from a dynamic language
;; without any optimization or type hinting.

;; In the original blog post I went through a number of faster versions:

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; In solveit-2 I explicitly typed the loop variables, inlined the function call, and made an
;; internal target for recur. That speeded things up considerably:

(defn solveit-2 [t0 y0 h its]
  (loop [t0 (double t0), y0 (double y0), h (double h), its (int its)]
    (if (> its 0) 
      (let [t1 (+ t0 h)
            y1 (+ y0 (* h (- t0 y0)))]
        (recur t1 y1 h (dec its)))
      [t0 y0 h its])))

;; Let's time that and see how it goes:
(cyclesperit (solveit-2 0.0 1.0) 10000000) 

;; On my old desktop with clojure 1.2 : 490 cycles

;; On my netbook with clojure 1.2 : 600 cycles
;; On my netbook with clojure 1.3 : 44 cycles
;; On my netbook with clojure 1.4 : 44 cycles


;; Wow!

;; Again we see that my netbook is around 20% slower, but the performance change between clojure 1.2
;; and clojure 1.3 was incredible.

;; There's now about a factor of 6 between the clojure version and the imaginary assembler program
;; on the imaginary cpu in my head. And all I've done is to declare the types and inline the functions.

;; This is actually so impressive that I want to examine the three things in detail, to see how the
;; three things combine to cause the speedup. As I remember from clojure 1.2, you needed all three
;; changes to see any great difference. In 1.4 you can type hint separately from making an internal loop target.
;; All timings with clojure 1.4.0.

(defn f [t y] (- t y))

(defn solveit [t0 y0 h its]
  (if (> its 0) 
    (let [t1 (+ t0 h)
          y1 (+ y0 (* h (f t0 y0)))]
      (recur t1 y1 h (dec its)))
    [t0 y0 h its]))

(cyclesperit (solveit 0.0 1.0) 10000000) ; 2481

;; type hints: 

(defn ^double tf [^double t ^double y] (- t y))

(defn solveit-1-1 [^double t0 ^double y0 ^double h ^long its]
  (if (> its 0) 
    (let [t1 (+ t0 h)
          y1 (+ y0 (* h (tf t0 y0)))]
      (recur t1 y1 h (dec its)))
    [t0 y0 h its]))

(cyclesperit (solveit-1-0 0.0 1.0) 10000000) ; 289

;; Inline f
(defn solveit-1-2 [ ^double t0 ^double y0 ^double h ^long its]
  (if (> its 0) 
    (let [t1 (+ t0 h)
          y1 (+ y0 (* h (- t0 y0)))]
      (recur t1 y1 h (dec its)))
    [t0 y0 h its]))

(cyclesperit (solveit-1-2 0.0 1.0) 10000000) ; 44

;; Internal loop target
(defn solveit-1-3 [ ^double t0 ^double y0 ^double h ^long its]
  (loop [t0 t0, y0 y0, h h, its its]
    (if (> its 0) 
      (let [t1 (+ t0 h)
            y1 (+ y0 (* h (- t0 y0)))]
        (recur t1 y1 h (dec its)))
      [t0 y0 h its])))

(cyclesperit (solveit-1-3 0.0 1.0) 10000000) ; 44

;; Original version
(defn solveit-2 [t0 y0 h its]
  (loop [t0 (double t0), y0 (double y0), h (double h), its (int its)]
    (if (> its 0) 
      (let [t1 (+ t0 h)
            y1 (+ y0 (* h (- t0 y0)))]
        (recur t1 y1 h (dec its)))
      [t0 y0 h its])))

;; Let's time that and see how it goes:
(cyclesperit (solveit-2 0.0 1.0) 10000000) ; 44


;; This is pretty awesome. Type hints alone are giving us a factor of eight speedup, and inlining
;; the function then gives us another factor of 6. The internal loop target, which was originally a
;; bit of a hack to allow clojure to deduce types doesn't make any difference to the hinted version,
;; but in fact that still works, and gives the compiler enough information to get the same speed.

;; So far I really could not be more impressed.

;; In the original post, I then used a profiler to find out where the loop was running slow, and 
;; did some strange things to make it fast:

(defn solveit-3 [t0 y0 h its]
  (loop [t0 (double t0), y0 (double y0), h (double h), its (int its)]
    (if (> its (int 0)) 
      (let [t1 (+ t0 h)
            y1 (+ y0 (* h (- t0 y0)))]
        (recur t1 y1 h (dec its)))
      [t0 y0 h its])))

(cyclesperit (solveit-3 0.0 1.0) 10000000) 

;; On my old desktop with clojure 1.2 : 70 cycles

;; On my netbook with clojure 1.2 : 90 cycles
;; On my netbook with clojure 1.3 : 100 cycles
;; On my netbook with clojure 1.4 : 100 cycles

;; Originally, that (int 0) instead of 0 in the comparison made all the difference, causing a huge
;; and unexpected speedup.

;; The odd thing here is that the same program actually runs slightly slower in clojure 1.3/1.4 than
;; it does in clojure 1.2. Not much, and they're all considerably slower than clojure 1.4's best
;; shot so far.

;; What I think is happening here is that that int cast is taking up unnecessary time in the later clojures,
;; whereas in 1.2 it was the final piece of the puzzle as far as the compiler was concerned.

;; To test this:
(defn solveit-3-a [ t0 y0 h its]
  (loop [t0 (double t0), y0 (double y0), h (double h), its (long its)]
    (if (> its (long 0 )) 
      (let [t1 (+ t0 h)
            y1 (+ y0 (* h (- t0 y0)))]
        (recur t1 y1 h (dec its)))
      [t0 y0 h its])))

(cyclesperit (solveit-3-a 0.0 1.0) 10000000) ; 44

;; solveit-3-a runs at the same 44 cycles/iteration as we've been seeing before. I get the
;; impression that modern clojure prefers longs to ints.


;; The final craziness was to take the (int 0) out of the loop entirely:

(defn solveit-4 [t0 y0 h its]
  (let [zero (int 0)]
    (loop [t0 (double t0) y0 (double y0) h (double h) its (int its)]
      (if (> its zero) 
        (let [t1 (+ t0 h)
              y1 (+ y0 (* h (- t0 y0)))]
          (recur t1 y1 h (dec its)))
      [t0 y0 h its]))))

(cyclesperit (solveit-4 0.0  1.0) 100000000)

;; And this is where it gets weird:

;; On my old desktop with clojure 1.2 : 23 cycles

;; On my netbook with clojure 1.2 : 32 cycles
;; On my netbook with clojure 1.3 : 45 cycles
;; On my netbook with clojure 1.4 : 45 cycles

;; Clojure 1.3/1.4 is now running at its usual (fast) speed, despite the ints rather than longs in solveit-4

;; But Clojure 1.2 on this program is able to run about 33% faster than the modern versions.

;; So it actually looks as though whatever awesomeness has caused the phenomenal speedup between
;; clojure 1.2 and 1.3 has actually slightly slowed down the optimized version.

;; I'm imagining that this is something to do with using longs rather than ints. 

;; But I haven't done any kind of close profiling to see if I can make 1.4 run any faster than 44 cycles/loop.

;; Someone in a comment on my earlier post said that they were seeing 8 cycles/second on their
;; Macintosh, so it might just be that there's something weird about the Atom processor.

;; At any rate, it looks as though optimizing numeric code in Clojure is now dead easy. You just
;; tell it the types, inline function calls, and that's as good as it gets.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; For those in real need of cycles and willing to take risks to save them, there's:

(set! *unchecked-math* true)

;; Which I don't think has an equivalent in 1.2, which cuts a few cycles off:

(defn solveit-4 [t0 y0 h its]
  (let [zero (long 0)]
    (loop [t0 (double t0) y0 (double y0) h (double h) its (long its)]
      (if (> its zero) 
        (let [t1 (+ t0 h)
              y1 (+ y0 (* h (- t0 y0)))]
          (recur t1 y1 h (dec its)))
      [t0 y0 h its]))))

(cyclesperit (solveit-4 0.0  1.0) 100000000) ; 37

;; I do wonder what is going on here. I'm not sure what checking is here to be turned off.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Incidentally, the following two java programs also produce the numbers that are seen above:

;; public class Euler{
;;     public static void main (String[] args){
;;         double cpuspeed = 1.662;
;;         int its = 10000000;
;;         double t=0;
;;         double y=0;
;;         double h=1.0/its;

;;         Long start = System.nanoTime();
;;         for(int i = 0; i < its; i++){
;;             y = y+h*(t-y);
;;             t = t+h;
;;         }
;;         Long finish = System.nanoTime();

;;         System.out.println("y=" + y + " t=" +t); 
;;         System.out.println("cycles/iteration: " +  ((int) ((cpuspeed * (finish - start)) / its)));
;;     }
;; }


;; y=0.3678794227282174 t=0.99999999975017
;; cycles/iteration: 32


;; public class Euler{
;;     public static void main (String[] args){
;;         double cpuspeed = 1.662;
;;         long its = 10000000;
;;         double t=0;
;;         double y=0;
;;         double h=1.0/its;

;;         long start = System.nanoTime();
;;         for(long i = 0; i < its; i++){
;;             y = y+h*(t-y);
;;             t = t+h;
;;         }
;;         long finish = System.nanoTime();

;;         System.out.println("y=" + y + " t=" +t); 
;;         System.out.println("cycles/iteration: " +  ((int) ((cpuspeed * (finish - start)) / its)));
;;     }
;; }

;; y=0.3678794227282174 t=0.99999999975017
;; cycles/iteration: 37