Learning Clojure

Clojure Setup Tutorial with EMACS: clojure, clojure.contrib, swank, slime, maven and maven-clojure-plugin in a couple of seconds

2011-09-22T13:28:00.000+01:00

As of 22nd September 2011, this is still working fine, this time on Fedora 14 (substituting yum for apt-get, but that's the only difference)

------------------------------------------------------------------------------------------------------------------------------------

I've just (15th February 2011) had to set this up again for someone on an Ubuntu 10.10 box.

And to my speechless amazement (16th February 2011), it worked absolutely the same way on Radek Ostrowski (a complete stranger) 's mac at dev8d, thus condemning him to years of suffering at the hands of Aquamacs.

This method is still my favourite, still the one I actually use in practice, and still the easiest, and it's been stable for about a year now.

Since I've just confirmed that it still works, I've moved it back to the top.

Install maven:

$sudo apt-get install maven2

Then create a pom.xml file to tell maven which repositories to use. There's an example below to copy and paste.

Once you've got maven installed, and made a pom.xml file, then:

$mvn clojure:repl

Will start a REPL with clojure 1.2 and clojure-contrib 1.2 on the classpath. The REPL should use jLine to give it editing and history on all platforms.

If you like to use EMACS instead, then

$mvn clojure:swank

Will start a swank server, which you can connect to with M-x slime-connect.

That's it from the clojure side.

To set up emacs, and equip it to talk to the clojure swank server:

$sudo apt-get install emacs

Then install the emacs lisp package archive (see http://tromey.com/elpa/) by evaluating this code:

(cut and paste it into the scratch buffer, put the cursor in the middle, and use M-C-x):

--emacs lisp to install elpa---------------

(let ((buffer (url-retrieve-synchronously
        "http://tromey.com/elpa/package-install.el")))
  (save-excursion
    (set-buffer buffer)
    (goto-char (point-min))
    (re-search-forward "^$" nil 'move)
    (eval-region (point) (point-max))
    (kill-buffer (current-buffer))))

---end of emacs lisp to install elpa-------------

Now use M-x package-list-packages to bring up the list of packages, and use i and then x to mark and then install slime, slime-repl, and clojure-mode.

Then connect emacs to the already running clojure image with M-x slime-connect

That should be it. You should now be at a running clojure 1.2 repl inside emacs.

Here's an example of a pom.xml file that pulls in clojure and clojure-contrib 1.2 . Just cut and paste it.

As well as the essentials, I've also added: the clojars.org repository, where many useful clojure packages live; the maven versions plugin, which helps with keeping everything cutting edge; and jline so that command line repls work better (mvn clojure:repl)

I like to have my startup repls conditioned a little, so if you have a startup script that you always want to run, add this snippet

<replScript>startup.clj</replScript>

</configuration>

to the clojure-maven-plugin section so that when maven starts a repl, the code in startup.clj is loaded as the first action. This is a good place to set print-length and print-level, so that they will be set before the swank server starts, which means that you won't hang emacs by evaluating an infinite sequence.

I also like to require everything on the classpath, so that I can use things like find-doc to find out about everything.

pom.xml files look terrifying, but they're really not. It's just that xml is such a godawful verbose way to write things out.

It contains: the addresses of three repositories which hold vital code; the names and versions of four vital jar files that you need; and the names and versions of two helpful maven plugins.

pom.xml

<project>

  <modelVersion>4.0.0</modelVersion>
  <groupId>com.example</groupId>
  <artifactId>hello-maven-clojure-swank</artifactId>

  <version>1.0-SNAPSHOT</version>
  <name>hello-maven</name>
  <description>maven, clojure, emacs: together at last</description>

  <repositories>
    <repository>
      <id>clojars</id>
      <url>http://clojars.org/repo/</url>
    </repository>
    <repository>
      <id>clojure</id>
      <url>http://build.clojure.org/releases</url>
    </repository>
    <repository>
      <id>central</id>
      <url>http://repo1.maven.org/maven2</url>
    </repository>
  </repositories>

  <dependencies>
    <dependency>
      <groupId>org.clojure</groupId>
      <artifactId>clojure</artifactId>
      <version>1.2.0</version>
    </dependency>
    <dependency>
      <groupId>org.clojure</groupId>
      <artifactId>clojure-contrib</artifactId>
      <version>1.2.0</version>
    </dependency>
    <dependency>
      <groupId>jline</groupId>
      <artifactId>jline</artifactId>
      <version>0.9.94</version>
    </dependency>
    <dependency>
      <groupId>swank-clojure</groupId>
      <artifactId>swank-clojure</artifactId>
      <version>1.2.1</version>
    </dependency>
  </dependencies>

  <build>
    <plugins>
      <plugin>
         <groupId>com.theoryinpractise</groupId>
            <artifactId>clojure-maven-plugin</artifactId>
         <version>1.3.3</version>
      </plugin>
      <plugin>
        <groupId>org.codehaus.mojo</groupId>
          <artifactId>versions-maven-plugin</artifactId>
        <version>1.2</version>
      </plugin>
    </plugins>

  </build>

</project>

Numerical Integration: Better Refinements?

2011-05-29T00:28:00.002+01:00


;; Numerical Integration: Better Refinements?

;; Here are some very simple functions which we might want to test integration
;; methods on:
(defn square  [x] (* x x))
(defn sine    [x] (Math/sin x))
(defn step    [x] (if (< x 1/2) 0.0 1.0))
(defn inverse [x] (/ x))

;; Here are some Newton-Cotes formulae for approximate integration:

(defn trapezium-rule [f a b]
  (* 1/2 (- b a) (+ (f a) (f b))))

(defn simpson-rule [f a b]
  (let [midpoint (+ a (/ (- b a) 2))]
    (* 1/6 (- b a) (+ (f a) (* 4 (f midpoint)) (f b)))))

(defn simpson38-rule [f a b]
  (let [midpoint1 (/ (+ a a b) 3)
        midpoint2 (/ (+ a b b) 3)]
    (* 1/8 (- b a) (+ (f a) (* 3 (f midpoint1)) (* 3 (f midpoint2)) (f b)))))

(defn booles-rule [f a b]
  (let [midpoint1 (/ (+ a a a b) 4)
        midpoint2 (/ (+ a a b b) 4)
        midpoint3 (/ (+ a b b b) 4)]
    (* 1/90 (- b a) (+ (* 7 (f a)) (* 32 (f midpoint1)) (* 12 (f midpoint2)) (* 32 (f midpoint3)) (* 7 (f b))))))

;; We can use any of these rules to get estimate of the integral of a function over an interval:
(simpson-rule inverse 1 3) ; 10/9

;; If we halve the interval and use the rule over both halves, then we can use the rule to
;; get a better estimate by adding the estimates for the half-intervals
(+
 (simpson-rule inverse 1 2)
 (simpson-rule inverse 2 3)) ; 11/10

;; We can guess at the error involved in the estimate by taking the difference
;; between these two estimates, on the basis that splitting the interval usually
;; makes most of the error go away.

(- (simpson-rule inverse 1 3)
   (+
    (simpson-rule inverse 1 2)
    (simpson-rule inverse 2 3))) ; 1/90

;; So we'd expect that the first estimate is out by a bit more than 1/90, and
;; that the second is out by rather less than 1/90

;; For the inverse function, which can be integrated symbolically, we know the
;; true answer:
(- (Math/log 3) (Math/log 1)) ; 1.0986122886681098
(/ 10.0 9) ; 1.1111111111111112
(/ 11.0 10) ; 1.1

;; So the errors are really:
(- 1.0986122 10/9)  ; -0.0124989111111109  ; which is ~ 1/90
(- 1.0986122 11/10) ; -0.00138780000000005 ; which is ~ 1/900

;; This method of guessing the error is deeply suspect, and can go wrong, but I
;; won't go into details.

;; I think it's good enough for our purposes as long as the functions we want to
;; integrate are reasonably well behaved and we take small enough intervals.

;; So we can easily make a function which gives us the more refined of the two
;; estimates, together with a guess as to how close it is to the truth.
(defn approx-with-error[rule f a b]
  (let [guess (rule f a b)
        midpoint (/ (+ a b) 2)
        better-guess (+ (rule f a midpoint) (rule f midpoint b))
        error-estimate (- guess better-guess)
        abs-error-estimate (if (> error-estimate 0) error-estimate (- error-estimate))]
    [better-guess abs-error-estimate]))


;; Let's try it out on a few cases, on the particularly nasty integral of 1/x over [0.01,100]

;; This is the true answer
(- (Math/log 100) (Math/log 0.01)) ; 9.210340371976184

(approx-with-error trapezium-rule inverse 0.01 100) ; [2500.999775019998 2499.000174980002]
;; We guess 2500, and we think we're out by at most 2499, which is just true
(approx-with-error simpson-rule inverse 0.01 100) ; [835.4437770204454 832.5559396728856]
(approx-with-error simpson38-rule inverse 0.01 100) ; [627.4427811912234 624.2442845054817]
(approx-with-error booles-rule inverse 0.01 100) ; [391.7824297523125 388.1576179566068]

;; When we split the interval into two halves [0.01, 50.05] [50.05,100]
(approx-with-error trapezium-rule inverse 0.01 50.05) ; [1252.2495505293746 1250.2503495705255]
(approx-with-error trapezium-rule inverse 50.05 100) ; [0.7072645364881702 0.041486462512828726]

;; Our guess tells us that the great majority of the error is in the first sub interval
;; We might want to refine that first, before bothering with the other one:
;; We'll now split [0.01, 25.025][25.025,50.05][50.05,100]
(approx-with-error trapezium-rule inverse 0.01 25.025) ; [626.6241012183343 624.6256989814656]
(approx-with-error trapezium-rule inverse 25.025 50.05) ; [0.7083333333333333 0.04166666666666663]
(approx-with-error trapezium-rule inverse 50.05 100) ; [0.7072645364881702 0.041486462512828726]

;; Again, one subinterval seems to be responsible for the majority of our errors.

;; We could keep a list of intervals, sorted by the estimated error, and always refine the one
;; with the largest guessed error.

(defn interval->errorstruct [rule f [a b]]
  (let [[guess error-guess] (approx-with-error rule f a b)]
    [error-guess, guess, [a,b]]))

(def errorstructs (map (partial interval->errorstruct trapezium-rule inverse)
                       [[0.01,25.025][25.025 50.05][50.05 100]]))

errorstructs
;; ([624.6256989814656 626.6241012183343     [0.01   25.025]]
;;  [0.04166666666666663 0.7083333333333333  [25.025 50.05]]
;;  [0.041486462512828726 0.7072645364881702 [50.05  100]])

;; And now we need a function to refine the interval with the largest error

(defn refine[rule f errorstructs]
  (let [sortedstructs (reverse (sort errorstructs))
        [_ _ [a b]] (first sortedstructs)
        remains (rest sortedstructs)
        midpoint (/ (+ a b) 2)
        subinterval1 (interval->errorstruct rule f [a midpoint])
        subinterval2 (interval->errorstruct rule f [midpoint b])] 
        (cons subinterval1 (cons subinterval2 remains))))

;; Now with every call to refine, we refine the interval with the largest error estimate
(refine trapezium-rule inverse errorstructs)
;; ([311.93889676733556 313.93570379188156 [0.01 12.5175]]
;;  [0.04159457274964806 0.7079060863675477 [12.5175 25.025]]
;;  [0.04166666666666663 0.7083333333333333 [25.025 50.05]]
;;  [0.041486462512828726 0.7072645364881702 [50.05 100]])


(def successive-trapezium-refinements (iterate (partial refine trapezium-rule inverse) errorstructs))


;; Here's what it looks like after a few iterations
(nth successive-trapezium-refinements 5)
;; ([18.81475746721543 20.764864658796014 [0.01 0.7917187499999999]]
;;  [0.040533241059784286 0.7015625070966475 [0.7917187499999999 1.5734374999999998]]
;;  [0.04109463731966401 0.7049306354620377 [1.5734374999999998 3.136875]]
;;  [0.041379306898238544 0.7066276281534741 [3.136875 6.26375]]
;;  [0.04152264848816445 0.7074793872568832 [6.26375 12.5175]]
;;  [0.04159457274964806 0.7079060863675477 [12.5175 25.025]]
;;  [0.04166666666666663 0.7083333333333333 [25.025 50.05]]
;;  [0.041486462512828726 0.7072645364881702 [50.05 100]])

        
;; We can get our best guess for the whole thing and our total error estimate
;; by reducing this list

(reduce + (map first  (nth successive-trapezium-refinements 5))) ; 19.104035002910425
(reduce + (map second (nth successive-trapezium-refinements 5))) ; 25.708968772954105


;; After a hundred refinements..
(reduce + (map first  (nth successive-trapezium-refinements 100))) ; 0.010431101535137086
(reduce + (map second (nth successive-trapezium-refinements 100))) ; 9.213824736866899

;; After a thousand refinements..
(reduce + (map first  (nth successive-trapezium-refinements 1000))) ; 1.0913238861095381E-4
(reduce + (map second (nth successive-trapezium-refinements 1000))) ; 9.210376750235199

;; That's not bad, (the real answer is 9.210340371976184), but it's running very slowly.

;; We could try with a higher order rule

(def successive-boole-refinements (iterate (partial refine booles-rule inverse) errorstructs))
(reduce + (map first   (nth successive-boole-refinements 1000))) ; 4.420942778526893E-15
(reduce + (map second  (nth successive-boole-refinements 1000))) ; 9.210340371976176

;; In this case, that seems to work very well, but the run time is appalling.

;; The problem is that we have a longer and longer list of intervals at every
;; step, and every step, we have to sort this list. That's an n^2 algorithm,
;; which won't scale well.

;; What we should do here is use a priority queue. Clojure doesn't have an
;; immutable version, although it's possible to fake one with a sorted map.

;; But rather than do that, I'm going to drop out of the functional paradigm
;; altogther, and use the heap implementation from Java in a mutable fashion,
;; looping and popping and adding.

(defn improve-loop [rule f a b count]
  (let [pq (java.util.PriorityQueue. count (comparator (fn[a b](> (first a)(first b)))))]
    (.add pq (interval->errorstruct rule f [a b]))
    (loop [pq pq count count]
      (if (zero? count) pq
          (let [[err val [a b]] (.poll pq)
                midpoint (/ (+ a b) 2)
                aa (interval->errorstruct rule f [a midpoint])
                bb (interval->errorstruct rule f [midpoint b])]
            (doto pq
              (.add aa)
              (.add bb))
            (recur pq (dec count)))))))

;; Now we can do our calculation much faster
(defn integrate [rule f a b count]
  (let [pq (improve-loop rule f a b count)]
    [(reduce + (map first pq))
     (reduce + (map second pq))]))

;; We'll ask for a thousand refinements, and get back the error estimate, and the answer.
(integrate booles-rule inverse 0.01 100 1000) ; [4.455637248046429E-15 9.21034037197618]


;; Let's try the same integral over the very nasty range [0.0000001, 10000000] which caused serious
;; problems for our previous methods.
;; The real answer is
(- (Math/log 10000000) (Math/log 0.00000001)) ; 34.538776394910684

;; And our approximations are:
(integrate booles-rule inverse 0.00000001 10000000 10) ; [3.797743055486256E10 3.797743056542089E10]
(integrate booles-rule inverse 0.00000001 10000000 100) ; [3.3430724324184924E-5 34.53877704296225]
(integrate booles-rule inverse 0.00000001 10000000 1000) ; [4.549938203979309E-11 34.53877639491147]
(integrate booles-rule inverse 0.00000001 10000000 10000) ; [9.361001557239845E-16 34.53877639491065]

;; For the non-stiff integrals that we started playing with, Boole's rule is great:
;; It's exact for quadratics, and several higher powers
(integrate booles-rule square 0 2 10) ; [0 8/3]
(integrate booles-rule (fn[x] (* x x x x)) 0 2 10) ; [0 32/5]
(integrate booles-rule (fn[x] (* x x x x x)) 0 2 10) ; [0 32/3]

;; and very good for higher powers, even with very few refinements
(integrate booles-rule (fn[x] (* x x x x x x)) 0 2 10) ; [969/8589934592 471219269093/25769803776]
(integrate booles-rule (fn[x] (* x x x x x x)) 0 2 20) ; [2127/1099511627776 60316066438099/3298534883328]

;; convergence is great for sine
(integrate booles-rule sine 0 Math/PI 10) ; [1.7383848804897184E-9 1.9999999999725113]
(integrate booles-rule sine 0 Math/PI 100) ; [3.1931922384043077E-15 1.9999999999999991]
(integrate booles-rule sine 0 Math/PI 1000) ; [2.526233413538588E-17 1.999999999999999]
(integrate booles-rule sine 0 Math/PI 10000) ; [6.32722651455846E-18 2.0]


;; But I'm still quite worried about the error estimate that we made. It's only
;; a guess, and it can be a bad guess.  Here are some functions that are
;; deliberately designed to screw things up.

;; This function is extremely vibratey near the origin.
(defn sineinverse[x] (Math/sin (/ x)))

;; The error estimates are clearly wrong here, but the answers seem to settle down to something that looks plausible.

(integrate booles-rule sineinverse 0.001 10 1) ; [0.09752245288534744 3.189170812427795]
(integrate booles-rule sineinverse 0.001 10 10) ; [0.014802407142066881 2.725700351059874]
(integrate booles-rule sineinverse 0.001 10 100) ; [2.666579898821515E-4 2.7262259059929814]
(integrate booles-rule sineinverse 0.001 10 1000) ; [7.84363268117651E-9 2.7262019887881457]
(integrate booles-rule sineinverse 0.001 10 10000) ; [8.311387750656713E-15 2.726201989096135]

;; I'm slightly reassured that if we use the trapezium rule, which should be
;; slower converging but less sensitive to high derivatives, we seem to settle down to the same thing:

(integrate trapezium-rule sineinverse 0.001 10 1) ; [0.3493754261290617 2.9603246206316127]
(integrate trapezium-rule sineinverse 0.001 10 10) ; [0.12037643221528535 2.759621850911819]
(integrate trapezium-rule sineinverse 0.001 10 100) ; [0.011584111090323689 2.728470051290524]
(integrate trapezium-rule sineinverse 0.001 10 1000) ; [7.174438961790802E-4 2.7265014884997414]
(integrate trapezium-rule sineinverse 0.001 10 10000) ; [1.0854830311799172E-5 2.7262023437746654]

;; Since I don't actually know what the integral of sin(1/x) is, I've no idea
;; whether this answer is correct. Since both rules seem to settle down to the
;; same answer, I tend to believe that it is.

;; Here's a weird function, which looks like it should be even worse that sin(1/x) on its own
(defn strange[x] (- (Math/sin (/ x)) (/ (Math/cos (/ x)) x)))

;; In fact it's the derivative of x sin(1/x), so we can calculate the real answer over [0.001, 10]
;; which should be:
(- (* 10 (Math/sin 1/10)) (* 0.001 (Math/sin 1000))) ; 0.9975072869277496

;; Interestingly, the error estimates look sound for this one:
(integrate booles-rule strange 0.001 10 1) ; [109.91706304856582 -108.12753277035351]
(integrate booles-rule strange 0.001 10 10) ; [0.07641821305025362 -1.0276123964492345]
(integrate booles-rule strange 0.001 10 100) ; [0.0798435700032961 1.0469088424961843]
(integrate booles-rule strange 0.001 10 1000) ; [2.0359056110968434E-6 0.9975072871949854]
(integrate booles-rule strange 0.001 10 10000) ; [1.9224976990340685E-12 0.997507286927752]

;; Since we seem to have dealt well with some nasty functions, we might be
;; getting confident in our rule. I know I was!

;; But this innocuous looking function
(defn sine80squared[x] (square (Math/sin (* x 80))))

;; Is a complete nightmare. It looks as though the method is converging well:
(integrate booles-rule sine80squared 0 Math/PI 1) ; [1.1091279850485843E-28 3.7074689566598855E-28]
(integrate booles-rule sine80squared 0 Math/PI 10) ; [0.013089969389960716 0.7853981633974437]
(integrate booles-rule sine80squared 0 Math/PI 100) ; [1.7991469747360833E-12 0.7853981633974478]
(integrate booles-rule sine80squared 0 Math/PI 1000) ; [3.207733520089899E-13 0.7853981633974484]
(integrate booles-rule sine80squared 0 Math/PI 10000) ; [3.1095566849572033E-15 0.7853981633974481]

;; But if we use a different rule, it also seems to converge, but to a completely different answer
(integrate trapezium-rule sine80squared 0 Math/PI 1) ; [3.226319244612108E-28 4.262878793991289E-28]
(integrate trapezium-rule sine80squared 0 Math/PI 10) ; [0.01573134053904405 0.19774924859401588]
(integrate trapezium-rule sine80squared 0 Math/PI 100) ; [5.4528883422580115E-5 0.19634904414812832]
(integrate trapezium-rule sine80squared 0 Math/PI 1000) ; [4.6327740574637914E-7 0.19634954102557595]
(integrate trapezium-rule sine80squared 0 Math/PI 10000) ; [5.200068066397881E-9 0.19634954084967793]

;; In fact both answers are wrong. We can calculate the real integral of this
;; function over the interval [0,pi] which should be:

(/ (Math/PI) 2) ; 1.5707963267948966

;; In fact if we use our much earlier 'divide every interval evenly' algorithm:
(defn iterated-rule [rule f a b N]
  (if (= N 0)
    (rule f a b)
    (let [midpoint (+ a (/ (- b a) 2))]
      (+ (iterated-rule rule f a midpoint (dec N))
         (iterated-rule rule f midpoint b (dec N))))))

;; We very quickly get surprisingly good answers:
(iterated-rule trapezium-rule sine80squared 0 Math/PI 1) ; 1.13079223568638E-28
(iterated-rule trapezium-rule sine80squared 0 Math/PI 2) ; 4.262878793991289E-28
(iterated-rule trapezium-rule sine80squared 0 Math/PI 3) ; 3.5626606693542546E-28
(iterated-rule trapezium-rule sine80squared 0 Math/PI 4) ; 3.6535380881685736E-28
(iterated-rule trapezium-rule sine80squared 0 Math/PI 5) ; 1.5707963267948966
(iterated-rule trapezium-rule sine80squared 0 Math/PI 6) ; 1.570796326794897
(iterated-rule trapezium-rule sine80squared 0 Math/PI 7) ; 1.5707963267948972
(iterated-rule trapezium-rule sine80squared 0 Math/PI 8) ; 1.5707963267948966
(iterated-rule trapezium-rule sine80squared 0 Math/PI 9) ; 1.5707963267948957
(iterated-rule trapezium-rule sine80squared 0 Math/PI 10) ; 1.5707963267948963

;; With Booles' rule the story is the same:
(iterated-rule booles-rule sine80squared 0 Math/PI 1) ; 3.1974110932118413E-28
(iterated-rule booles-rule sine80squared 0 Math/PI 2) ; 3.7074689566598855E-28
(iterated-rule booles-rule sine80squared 0 Math/PI 3) ; 2.2340214425527414
(iterated-rule booles-rule sine80squared 0 Math/PI 4) ; 1.5358897417550046
(iterated-rule booles-rule sine80squared 0 Math/PI 5) ; 1.570796326794895
(iterated-rule booles-rule sine80squared 0 Math/PI 6) ; 1.5707963267948961
(iterated-rule booles-rule sine80squared 0 Math/PI 7) ; 1.5707963267948966
(iterated-rule booles-rule sine80squared 0 Math/PI 8) ; 1.5707963267948968
(iterated-rule booles-rule sine80squared 0 Math/PI 9) ; 1.5707963267948963


;; So the nice rule that we'd come up with, which worked so well for the stiff problem
;; of integrating 1/x near the origin, is completely broken on something that the obvious
;; recursion integrates (suspiciously) well.

;; The problem is that we used an error estimate that we can't trust to control
;; the refinement process.

;; If we guess that a certain interval contains hardly any error, then it will
;; never get refined at all, so we'll never find out that our guess is wrong.

Numerical Integration: Better Refinements?

2011-05-26T01:51:00.002+01:00


;; Numerical Integration: Better Refinements?

;; Here are some Newton-Cotes formulae for approximate integration:

(defn trapezium-rule [f a b]
  (* 1/2 (- b a) (+ (f a) (f b))))

(defn simpson-rule [f a b]
  (let [midpoint (+ a (/ (- b a) 2))]
    (* 1/6 (- b a) (+ (f a) (* 4 (f midpoint)) (f b)))))

(defn simpson38-rule [f a b]
  (let [midpoint1 (/ (+ a a b) 3)
        midpoint2 (/ (+ a b b) 3)]
    (* 1/8 (- b a) (+ (f a) (* 3 (f midpoint1)) (* 3 (f midpoint2)) (f b)))))

(defn booles-rule [f a b]
  (let [midpoint1 (/ (+ a a a b) 4)
        midpoint2 (/ (+ a a b b) 4)
        midpoint3 (/ (+ a b b b) 4)]
    (* 1/90 (- b a) (+ (* 7 (f a)) (* 32 (f midpoint1)) (* 12 (f midpoint2)) (* 32 (f midpoint3)) (* 7 (f b))))))


;; And here is a way to apply them to (power 2 N) subintervals
(defn iterated-rule [rule f a b N]
  (if (= N 0)
    (rule f a b)
    (let [midpoint (+ a (/ (- b a) 2))]
      (+ (iterated-rule rule f a midpoint (dec N))
         (iterated-rule rule f midpoint b (dec N))))))

;; Here are some very simple functions which we might want to test integration
;; methods on:
(defn square  [x] (* x x))
(defn sine    [x] (Math/sin x))
(defn step    [x] (if (< x 1/2) 0.0 1.0))
(defn inverse [x] (/ x))

;; We might notice that with our Newton-Cotes formulae, the change in the
;; estimate when we make a refinement is often around the same size as the
;; actual error.

;; Let's try using that to produce a method where we can say what error we'd
;; like, and if the change from refining the guess is too big, we should refine
;; more. And if we do refine more, we'll ask for half the desired error on each
;; of the two subintervals, which should mean that the total error is at least
;; comparable to the error we asked for!

(defn adaptive-rule-recurse [rule f a b desired-error]
  (let [guess (rule f a b)
        midpoint (/ (+ a b) 2)
        better-guess (+ (rule f a midpoint) (rule f midpoint b))
        error-estimate (- guess better-guess)
        abs-error-estimate (if (> error-estimate 0) error-estimate (- error-estimate))]
    (if (< abs-error-estimate desired-error) better-guess
        (let [half-desired-error (/ desired-error 2)]
          (+ (adaptive-rule-recurse rule f a midpoint half-desired-error)
             (adaptive-rule-recurse rule f midpoint b half-desired-error))))))

(adaptive-rule-recurse trapezium-rule square 0. 2 0.1) ; 2.6875
(adaptive-rule-recurse trapezium-rule square 0. 2 0.01) ; 2.66796875
(adaptive-rule-recurse trapezium-rule square 0. 2 0.001) ; 2.6669921875

(adaptive-rule-recurse trapezium-rule step 0.0001 1 0.1) ; 0.5
(adaptive-rule-recurse trapezium-rule step 0.0001 1 0.01) ; 0.5
(adaptive-rule-recurse trapezium-rule step 0.0001 1 0.001) ; 0.5
(adaptive-rule-recurse trapezium-rule step 0.0001 1 0.0001) ; 0.5

;; Here we're using the trapezium rule on the integral that we were previously unable to get
;; a good answer for: remember that the correct answer is 9.210340371976182

(adaptive-rule-recurse trapezium-rule inverse 0.0001 1 0.1)    ; 9.234002964342716
(adaptive-rule-recurse trapezium-rule inverse 0.0001 1 0.01)   ; 9.211881820961814
(adaptive-rule-recurse trapezium-rule inverse 0.0001 1 0.001)  ; 9.210518109406467
(adaptive-rule-recurse trapezium-rule inverse 0.0001 1 0.0001) ; 9.210358164670637

;; At this point, we're getting much better answers than we ever got before, but they've started
;; taking noticeable time to compute.

;; However, we still retain the option of using the higher order formulae:

(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.1) ; 9.210347324857782
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.01) ; 9.210345994304586
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.001) ; 9.210344014376869
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.0001) ; 9.21034116413936
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.00001) ; 9.210340404907965
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.000001) ; 9.210340376142819
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.0000001) ; 9.210340372376441
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.00000001) ; 9.210340372016345
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.000000001) ; 9.21034037198001
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.0000000001) ; 9.210340371976589
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.00000000001) ; 9.210340371976214
(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.000000000001) ; 9.210340371976185

;; These answers come back instantaneously, even when we're calculating the
;; previously impossible answer to the question to 12 decimal places.

;; But beware! If we ask for too much accuracy, then the effect of floating
;; point noise means that our recursion may never terminate:

(adaptive-rule-recurse booles-rule inverse 0.0001 1 0.0000000000001) ; freezes REPL!


;; At this point, I'm very tempted to say:

(defn integrate [f a b]
  (adaptive-rule-recurse booles-rule f a b 0.00000001))

;; Famously, the integral of 1/x diverges (very slowly) as x gets large
;; which makes it quite difficult to calculate for large intervals without using special tricks.

;; The integral from 1 to 1000000 of 1/x is log 1000000 - log 1

(- (Math/log 1000000) (Math/log 1)) ; 13.815510557964274

;; Let's try calculating that with our new integrate function:

(time
 (integrate (fn[x] (/ 1.0 x)) 1 1000000))

;; "Elapsed time: 293.219062 msecs"
;; 13.81551055800455

;; I think that's pretty good.

;; Unfortunately, we can still lock up the computer by asking the wrong question

(- (Math/log 100000000) (Math/log 1)) ; 18.420680743952367
;; but:
(time
 (integrate (fn[x] (/ 1.0 x)) 1 100000000)) ; Freezes REPL


;; Can you work out what's going on and what we could do about it?

Numerical Integration: Better Rules?

2011-05-26T00:27:00.001+01:00


;; Numerical Integration: Better Rules?

;; So far we've found a couple of approximations to 'the area under a graph'.

(defn trapezium-rule [f a b]
  (* 1/2 (- b a) (+ (f a) (f b))))

(defn simpson-rule [f a b]
  (let [midpoint (+ a (/ (- b a) 2))]
    (* 1/6 (- b a) (+ (f a) (* 4 (f midpoint)) (f b)))))

;; And a way of repeatedly splitting intervals and applying the rules to the sub-intervals
;; to produce better and better approximations.
(defn iterated-rule [rule f a b N]
  (if (= N 0)
    (rule f a b)
    (let [midpoint (+ a (/ (- b a) 2))]
      (+ (iterated-rule rule f a midpoint (dec N))
         (iterated-rule rule f midpoint b (dec N))))))


;; For these two functions, which are nice and smooth, with derivatives that
;; don't get too large, these rules produce good approximations
(defn square[x] (* x x))
(defn sine [x] (Math/sin x))

;; For this one, which isn't smooth, the approximations converge slowly:
(defn step [x]  (if (< x 1/2) 0.0 1.0))

;; For this, which is smooth, but which blows up near 0, the approximations are bad
(defn inverse [x] (/ x))


;; One approach which is often taken is to use 'more accurate' rules.

;; The trapezium rule and Simpson's rule are members of a family called 'Newton-Cotes' formulas.
;; The more complicated a Newton-Cotes formula is, the more polynomials it can integrate exactly.

;; But in fact, for nice smooth functions, they tend to be more accurate the
;; more polynomials they can integrate exactly.

;; Here are two more examples, and their estimates for the integral of sine over
;; a half circle (which should be exactly 2):

(defn simpson38-rule [f a b]
  (let [midpoint1 (/ (+ a a b) 3)
        midpoint2 (/ (+ a b b) 3)]
    (* 1/8 (- b a) (+ (f a) (* 3 (f midpoint1)) (* 3 (f midpoint2)) (f b)))))

(simpson38-rule sine 0 Math/PI) ; 2.040524284763495

(defn booles-rule [f a b]
  (let [midpoint1 (/ (+ a a a b) 4)
        midpoint2 (/ (+ a a b b) 4)
        midpoint3 (/ (+ a b b b) 4)]
    (* 1/90 (- b a) (+ (* 7 (f a)) (* 32 (f midpoint1)) (* 12 (f midpoint2)) (* 32 (f midpoint3)) (* 7 (f b))))))

(booles-rule sine 0 Math/PI) ; 1.9985707318238355

;; We can use the same approach to getting better estimates by subdividing: For
;; sine, as well as getting better estimates to start with, these rules have
;; high rates of convergence. It doesn't take many subdivisions before we start
;; to run into the limits of floating point accuracy.

(iterated-rule booles-rule sine 0 Math/PI 0) ; 1.9985707318238355
(iterated-rule booles-rule sine 0 Math/PI 1) ; 1.9999831309459855
(iterated-rule booles-rule sine 0 Math/PI 2) ; 1.9999997524545716
(iterated-rule booles-rule sine 0 Math/PI 3) ; 1.9999999961908446
(iterated-rule booles-rule sine 0 Math/PI 4) ; 1.9999999999407074
(iterated-rule booles-rule sine 0 Math/PI 5) ; 1.999999999999074
(iterated-rule booles-rule sine 0 Math/PI 6) ; 1.9999999999999853
(iterated-rule booles-rule sine 0 Math/PI 7) ; 1.9999999999999993
(iterated-rule booles-rule sine 0 Math/PI 8) ; 1.9999999999999998
(iterated-rule booles-rule sine 0 Math/PI 9) ; 1.9999999999999998


(iterated-rule simpson38-rule sine 0 Math/PI 0) ; 2.040524284763495
(iterated-rule simpson38-rule sine 0 Math/PI 1) ; 2.002009846628558
(iterated-rule simpson38-rule sine 0 Math/PI 2) ; 2.000119386415226
(iterated-rule simpson38-rule sine 0 Math/PI 3) ; 2.000007370036249
(iterated-rule simpson38-rule sine 0 Math/PI 4) ; 2.000000459216732
(iterated-rule simpson38-rule sine 0 Math/PI 5) ; 2.0000000286790867
(iterated-rule simpson38-rule sine 0 Math/PI 6) ; 2.0000000017921
(iterated-rule simpson38-rule sine 0 Math/PI 7) ; 2.000000000112001
(iterated-rule simpson38-rule sine 0 Math/PI 8) ; 2.0000000000069997
(iterated-rule simpson38-rule sine 0 Math/PI 9) ; 2.000000000000438


;; For the step function, however, boole's rule isn't really that much better
;; than the trapezium rule

(iterated-rule booles-rule step 0 1 0) ; 0.5666666666666667
(iterated-rule booles-rule step 0 1 1) ; 0.5388888888888889
(iterated-rule booles-rule step 0 1 2) ; 0.5194444444444445
(iterated-rule booles-rule step 0 1 3) ; 0.5097222222222222
(iterated-rule booles-rule step 0 1 4) ; 0.5048611111111111
(iterated-rule booles-rule step 0 1 5) ; 0.5024305555555555
(iterated-rule booles-rule step 0 1 6) ; 0.5012152777777777
(iterated-rule booles-rule step 0 1 7) ; 0.5006076388888889
(iterated-rule booles-rule step 0 1 8) ; 0.5003038194444445
(iterated-rule booles-rule step 0 1 9) ; 0.5001519097222222

;; Both seem to need double the number of points to halve their error.

(iterated-rule trapezium-rule step 0 1 0) ; 0.5
(iterated-rule trapezium-rule step 0 1 1) ; 0.75
(iterated-rule trapezium-rule step 0 1 2) ; 0.625
(iterated-rule trapezium-rule step 0 1 3) ; 0.5625
(iterated-rule trapezium-rule step 0 1 4) ; 0.53125
(iterated-rule trapezium-rule step 0 1 5) ; 0.515625
(iterated-rule trapezium-rule step 0 1 6) ; 0.5078125
(iterated-rule trapezium-rule step 0 1 7) ; 0.50390625
(iterated-rule trapezium-rule step 0 1 8) ; 0.501953125
(iterated-rule trapezium-rule step 0 1 9) ; 0.5009765625

;; Performance is similarly bad for boole's rule on the inverse function

(iterated-rule booles-rule inverse 0.0001 1 0) ; 779.9400477089192
(iterated-rule booles-rule inverse 0.0001 1 1) ; 391.7824297523124
(iterated-rule booles-rule inverse 0.0001 1 2) ; 198.04899407170583
(iterated-rule booles-rule inverse 0.0001 1 3) ; 101.52648013049377
(iterated-rule booles-rule inverse 0.0001 1 4) ; 53.607079025964396
(iterated-rule booles-rule inverse 0.0001 1 5) ; 29.984599583713347
(iterated-rule booles-rule inverse 0.0001 1 6) ; 18.501553032592827
(iterated-rule booles-rule inverse 0.0001 1 7) ; 13.071085113384171
(iterated-rule booles-rule inverse 0.0001 1 8) ; 10.635953348374398
(iterated-rule booles-rule inverse 0.0001 1 9) ; 9.647557264415854


(iterated-rule trapezium-rule inverse 0.0001 1 0) ; 4999.99995
(iterated-rule trapezium-rule inverse 0.0001 1 1) ; 2500.9997750199977
(iterated-rule trapezium-rule inverse 0.0001 1 2) ; 1251.8327765203333
(iterated-rule trapezium-rule inverse 0.0001 1 3) ; 627.5916420975113
(iterated-rule trapezium-rule inverse 0.0001 1 4) ; 315.8156999790102
(iterated-rule trapezium-rule inverse 0.0001 1 5) ; 160.27209317742054
(iterated-rule trapezium-rule inverse 0.0001 1 6) ; 82.84288624590312
(iterated-rule trapezium-rule inverse 0.0001 1 7) ; 44.46707207824598
(iterated-rule trapezium-rule inverse 0.0001 1 8) ; 25.61044440290146
(iterated-rule trapezium-rule inverse 0.0001 1 9) ; 16.499072595032356


;; So although it seems these higher order Newton-Cotes formulae are much more
;; accurate and faster converging on well behaved functions, they don't seem to
;; help much in integrating anything tricky.

Numerical Integration: Harder Functions

2011-05-25T23:26:00.001+01:00


;; Numerical Integration: Harder Functions

;; So far we've found a couple of approximations to 'the area under a graph'.

(defn trapezium-rule [f a b]
  (* 1/2 (- b a) (+ (f a) (f b))))

(defn simpson-rule [f a b]
  (let [midpoint (+ a (/ (- b a) 2))]
    (* 1/6 (- b a) (+ (f a) (* 4 (f midpoint)) (f b)))))

;; And a way of repeatedly splitting intervals and applying the rules to the sub-intervals
;; to produce better and better approximations.
(defn iterated-rule [rule f a b N]
  (if (= N 0)
    (rule f a b)
    (let [midpoint (+ a (/ (- b a) 2))]
      (+ (iterated-rule rule f a midpoint (dec N))
         (iterated-rule rule f midpoint b (dec N))))))

;; We know that Simpson's rule is exact for our original function

(defn square[x] (* x x))

(simpson-rule square 0 2) ; 8/3
(iterated-rule simpson-rule square 0 2 1) ; 8/3
(iterated-rule simpson-rule square 0 2 2) ; 8/3

;; How do we do with some other functions?

;; Another function where we can calculate the answer directly is sine
(defn sine [x] (Math/sin x))

;; The integral of sine over a half-turn (i.e. over the interval [0,pi]) is exactly 2.

;; here are the first few approximations with the trapezium rule:
(take 10 (map (partial iterated-rule trapezium-rule sine 0 Math/PI) (iterate inc 0)))
;; (1.9236706937217898E-16 1.5707963267948968 1.8961188979370398
;; 1.9742316019455508 1.993570343772339 1.9983933609701445 1.9995983886400377
;; 1.9998996001842024 1.9999749002350526 1.9999937250705755)

;; After an inauspicious start, the answer is clearly settling down to 2.

;; Simpson's rule also settles down to 2, but it's much quicker.
(take 10 (map (partial iterated-rule simpson-rule sine 0 Math/PI) (iterate inc 0)))
;; (2.094395102393196 2.0045597549844216 2.000269169948388 2.0000165910479364
;; 2.000001033369414 2.0000000645300027 2.000000004032258 2.000000000252003
;; 2.000000000015751 2.000000000000985)

;; So it looks like we're onto a winner so far.

;; Another type of function that doesn't do so well is:

(defn step [x]  (if (< x 1/2) 0 1))

;; It should be fairly obvious that the integral of this over [0,1] is 0.5, but
;; in fact the convergence of the methods is very slow compared to what we had
;; for sine or square.

(take 10 (map (partial iterated-rule trapezium-rule step 0. 1) (iterate inc 0)))
;; (0.5 0.75 0.625 0.5625 0.53125 0.515625 0.5078125 0.50390625 0.501953125 ...

(take 10 (map (partial iterated-rule simpson-rule step 0. 1) (iterate inc 0)))
;; (0.8333333333333336 0.5833333333333335 0.5416666666666667 0.5208333333333335
;; 0.5104166666666667 0.5052083333333335 0.5026041666666667 0.5013020833333335
;; 0.5006510416666667 0.5003255208333335 )

;; Notice how if we make the tiniest possible change to the function, which doesn't change the integral at all:
(defn evil-step [x]  (if (<= x 1/2) 0 1))

;; the approximations are all changed quite a lot, although they do seem to
;; still converge to the correct answer.

(take 10 (map (partial iterated-rule trapezium-rule evil-step 0. 1) (iterate inc 0)))
;; (0.5 0.25 0.375 0.4375 0.46875 0.484375 0.4921875 0.49609375 0.498046875 0.4990234375)
(take 10 (map (partial iterated-rule simpson-rule evil-step 0. 1) (iterate inc 0)))
;; (0.1666666666666667 0.4166666666666668 0.4583333333333335 0.4791666666666668
;; 0.4895833333333335 0.4947916666666668 0.4973958333333335 0.4986979166666668
;; 0.4993489583333335 0.4996744791666668)



;; What about if we're integrating a (moderately) badly behaved function?:

(defn inverse [x] (/ x))

;; the real answer, by devious mathematical trickery is (log a) - (log b)

;; So if we integrate over the region [0.0001, 1], where the values of 1/x are
;; sometimes very large, we should get:

(- (Math/log 1) (Math/log 0.0001)) ; 9.210340371976182

(take 10 (map (partial iterated-rule trapezium-rule inverse 0.0001 1) (iterate inc 0)))
;; (4999.99995 2500.9997750199977 1251.8327765203333 627.5916420975113
;; 315.8156999790102 160.27209317742054 82.84288624590312 44.46707207824598
;; 25.61044440290146 16.499072595032356)

(take 10 (map (partial iterated-rule simpson-rule inverse 0.0001 1) (iterate inc 0)))
;; (1667.9997166933313 835.4437770204453 419.5112639565708 211.89038593950997
;; 108.42422424355732 57.03315060206397 31.67513402236028 19.324901844453297
;; 13.461948659075995 10.812578055293251)

;; It looks like both the rules, after starting off with appalling first
;; guesses, are converging to something, and maybe even to the right answer, but
;; the convergence is very slow.  Remember that to get the tenth element of this
;; sequence we're splitting the interval up into 1024 pieces, and we haven't
;; even got the answer right to one significant figure.

;; Numerical Integration can give very misleading answers if it's used naively.
;; Quite often, the results of a careless numerical simulation will just 'look
;; wrong' to a trained eye, and in those cases it's usually the eye that's in
;; the right.

;; One way to be reassured that your answer is roughly correct is to alter
;; expected degree of accuracy of the method, and to check that the answer
;; doesn't change by much!

Numerical Integration: Better Approximations

2011-05-25T22:37:00.002+01:00


;; Numerical Integration: Better Approximations

;; Remember our example function:

(defn square[x] (* x x))

;; We're trying to find ways to calculate the 'area under its graph between 0
;; and 2', or the 'integral over the interval [0,2]'

;; A better estimate of this area under the graph is to imagine a trapezium
;; which has its base corners at (0,0) and (2,0), and the top corners at
;; (0,(square 0)) (2, (square 2)), and calculate the area of that.

;; More generally, if the interval we're interested is [a,b], and the function's
;; values there are fa and fb, then the area of the trapezium will just be:

(defn trapezium [a fa b fb]
  (* 1/2 (- b a) (+ fa fb)))

;; Another way to think about that is that it's the length of the interval
;; multiplied by the average of the values of the function at the ends.

;; So another approximation to the integral of f over [a,b] is:

(defn trapezium-rule [f a b]
  (trapezium a (f a) b (f b)))

(trapezium-rule square 0 2) ; 4

;; We can make another approximation by using the trapezium rule on the two
;; subintervals [0,1] and [1,2] and adding the results

(+ (trapezium-rule square 0 1)
   (trapezium-rule square 1 2)) ; 3

;; And an even better one by splitting those subintervals in half

(+ (trapezium-rule square 0 1/2)
   (trapezium-rule square 1/2 1)
   (trapezium-rule square 1 3/2)
   (trapezium-rule square 3/2 2)) ; 11/4

;; And so on...
(defn iterated-rule [rule f a b N]
  (if (= N 0)
    (rule f a b)
    (let [midpoint (+ a (/ (- b a) 2))]
      (+ (iterated-rule rule f a midpoint (dec N))
         (iterated-rule rule f midpoint b (dec N))))))

;; This converges fairly nicely:
(map (partial iterated-rule trapezium-rule square 0 2) (range 10))
;; (4 3 11/4 43/16 171/64 683/256 2731/1024 10923/4096 43691/16384 174763/65536)

(map (partial - 8/3)
     (map (partial iterated-rule trapezium-rule square 0 2) (range 10)))
;; (-4/3 -1/3 -1/12 -1/48 -1/192 -1/768 -1/3072 -1/12288 -1/49152 -1/196608)

;; We now only need a thousand samples of the function to get the answer
;; accurate to one part in 100000.

;; But an even nicer approximation is Simpson's rule, which involves fitting a
;; parabola to the two end points and the midpoint, and calculating the area
;; under that.

;; That's equivalent to taking a weighted sum of the values of f at the
;; beginning, midpoint, and end of the interval, with weights 1:4:1

(defn simpson-rule [f a b]
  (let [midpoint (+ a (/ (- b a) 2))]
    (* 1/6 (- b a) (+ (f a) (* 4 (f midpoint)) (f b)))))

;; For the square function, which is itself a parabola, this rule actually
;; calculates the area exactly!
(simpson-rule square 0 2) ; 8/3

;; That about wraps it up for numerical approximations to the integrals of
;; quadratics, but they are easy to calculate exactly anyway.

Numerical Integration: What is an Integral?

2011-05-25T22:07:00.001+01:00

I've found myself needing to take a few difficult integrals recently, as part of some machine learning algorithms, and I ended up calculating them with numerical approximations.

Frustration at the amount of time they were taking lead me to progressively improve my methods, which got very much faster as I played with them.

Weirdly, my initial instincts had been to translate my programs into either C or Java, which might well have resulted in 100x speed ups over the obvious implementations in Clojure, but in fact the clarity of expression afforded by writing in a language where recursion is a natural thing to do lead me, by making fairly obvious changes, to some very fast algorithms.

Of course, I could still translate them to C or Java, or optimized Clojure, and they'd still get a lot faster, but it hardly seems worth the bother now.

;; Numerical Integration I: What is an integral, anyway?

;; What does it mean to 'evaluate the integral' of a function over an interval?

;; Let's take an example function

(defn square[x] (* x x))

;; We could evaluate it at a few input values:
(map square [0 1 2]) ; (0 1 4)

;; The integral of a function over an interval is a sort of
;; 'average sum over all the possible values in the interval'.

;; A first approximation to the integral of square over the interval [0,2]
;; would be to say:

;; Let's pretend that the value of square between 0 and 1 is like (square 0)
;; and the value between 1 and 2 is like (square 1)

;; So to find our approximation we add (square 0) to (square 1):

(reduce + (map square [0 1])) ; 1

;; Another approximation, equally valid, would be to say that the value between
;; 0 and 1 is like (square 1) so we could equally well add the value at 1 to the
;; value at 2

(reduce + (map square [1 2])) ; 5

;; These answers are quite different! One is too low, because we paid too much
;; attention to the numbers at the start of the sub-ranges. One is too high,
;; because we added up numbers from the ends.

;; We can make both answers more accurate by sampling the function at more
;; points, and dividing the sum by an appropriate factor to make up for having
;; more points.

(/ (reduce + (map square [0 1/2 1 3/2]  )) 2) ; 7/4
(/ (reduce + (map square   [1/2 1 3/2 2])) 2) ; 15/4

;; We can continue the procedure, sampling the function at twice as many points again:

(/ (reduce + (map square [0 1/4 1/2 3/4 1 5/4 6/4 7/4]  )) 4) ; 35/16
(/ (reduce + (map square   [1/4 1/2 3/4 1 5/4 6/4 7/4 2])) 4) ; 51/16

;; It's clear that these values will get closer to one another the more points
;; we sample at

;; If we continue the procedure indefinitely:
(defn riemann-lower-sum [f a b N]
  (let [points (range a b (/ N))]
    (/ (reduce + (map f points)) N)))

(map (partial riemann-lower-sum square 0 2) (iterate inc 1)) ; (1 7/4 55/27 35/16 57/25 253/108 ...)
;; We find that the answers settle down to a particular number, and in this case
;; it looks like they're settling down to 8/3

;; However we do need to take quite a lot of samples to make the approximation get close!
(take 15 (map (partial - 8/3)
              (map (partial riemann-lower-sum square 0 2)
                   (iterate (partial * 2) 1))))
;; (5/3 11/12 23/48 47/192 95/768 191/3072 383/12288 767/49152 1535/196608
;; 3071/786432 6143/3145728 12287/12582912 24575/50331648 49151/201326592
;; 98303/805306368)

;; The last number here is what we get if we keep splitting our sub-intervals in
;; half fifteen times, and that's quite a lot of sub-intervals. (1, 2, 4, 8, 16 ...)

;; Strictly speaking, the riemann lower sum is a sum of the lowest values of
;; the function on the subintervals, but because the square function is
;; monotonically increasing over our range, the function as we defined it has
;; the same value.

;; If this procedure has a limit, then the integral is defined as this limit.
;; For well behaved functions, this limit will exist.

;; The limit can be thought of as the area under the graph of f from a to b.

;; At school we learn to calculate this limit by finding an antiderivative to
;; the original function, and calculating how much it changes over the interval.

;; An antiderivative of (fn[x] (* x x)) is (fn[x] (* 1/3 x x x))

(let [anti (fn[x] (* 1/3 x x x))]
  (- (anti 2) (anti 0))) ; 8/3

;; Although this is clearly a much better way to calculate integrals than by
;; taking millions of function samples and adding them up, it's often quite hard
;; to find antiderivatives even for quite simple functions, and for some quite
;; simple functions, no easily expressible antiderivative exists.

;; So we quite often find ourselves needing to calculate approximations by computer,
;; and this is called numerical integration.

Clojure Inference (reduce, cl-format, Poisson distribution, Bayes)

2011-05-04T01:06:00.000+01:00

Apologies for the mathematical nature of this post. It's just what I happen to be thinking about at the moment.

Although actually apart from the formula for a Poisson distribution, it's all just philosophy and counting and adding up and taking averages.

I published it here because:
(a) It's a nice example of how lovely a language clojure is for exploring mathematics.
(b) It's got reduce in it. a lot.
(c) I use cl-format to print lists of lists of floats, and I had to look through the common lisp docs to remember how to do that. I couldn't find any nice clojure examples to crib from.

;; A Poisson Process

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

;; Imagine a very lot of things, each of which has a very small chance of
;; happening, so that on average, only a few things happen every second.

;; A classic case is radioactivity, where we've got a very large number of
;; atoms, each of which has a very small chance of exploding in any given
;; second.

;; We might know the average number of decays per second without having the
;; slightest idea what the large number or the small chance are.

;; Alternatively, we might just know what the sample did over a one second
;; period, and be trying to guess what the constant is.

;; In probability, we model this situation with a Poisson process.

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

;; This cryptic function gives the probabilities, for a given decay rate l
;; of each of the numbers we might read on our Geiger counter.

;; I'm ashamed to say that although I wrote this function and can prove these
;; are the probabilities, I can't come up with an easy way to explain why it
;; gives the right answers. You have to imagine what the probabilities are for
;; larger and larger numbers but smaller and smaller chances, keeping the
;; average the same, and notice that they settle down as the numbers get
;; bigger. These are the limits:

(defn poisson [l]
  (map first
       (iterate
        (fn [[a c]] [(/(* a l) (inc c)) (inc c)])
        [(Math/exp (- l)) 0])))

;; So let's just take this as the definition of a Poisson process, and have a
;; look at what it means.

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

;; Clojure has inherited common lisp's excellent format function, so by casting
;; one cryptic spell we can print out the probabilities for the first few geiger
;; counter readings, if the average decay rate is one click/second.

(clojure.pprint/cl-format nil "~{~$~^, ~}" (take 10 (poisson 1)))
;; "0.37, 0.37, 0.18, 0.06, 0.02, 0.00, 0.00, 0.00, 0.00, 0.00"

;; One way of looking at this is to say that if we ran the experiment one
;; hundred times, then we'd expect that 37 times we'd hear no clicks, 37 times
;; we'd get one click, 18 times we'd get two, 6 times we'd get three, and twice
;; we'd get four. Sort of. Although actually we'd be quite surprised if that
;; happened.

;; As a sanity check, let's work out what our expected number of clicks is here.
;; We know it should be one.

(reductions +
            (map *
                 (poisson 1)
                 (range)))

;; (0.0 0.36787944117144233 0.7357588823428847 0.9196986029286058 0.9810118431238463 0.9963401531726563 0.9994058151824183 0.999916758850712 0.9999897508033253 0.999998874797402 0.9999998885745216 0.9999999899522336 0.9999999991683892 0.9999999999364022 0.9999999999954802 0.9999999999997 0.9999999999999813 0.9999999999999989 0.9999999999999999 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ...)

;; If we include enough terms we do indeed get one.

;; I enjoyed that. Let's try it again for a process with constant 7.5

(reductions +
            (map *
                 (poisson 7.5)
                 (range)))

;;(0.0 0.0041481327761087525 0.0352591285969244 0.15192536292498304 0.44359094874512967 0.9904639221579046 1.810773382277067 2.83616020742602 3.934788948657041 4.964753393561123 5.823057097647858 6.46678487571291 6.905690178939081 7.180005993455438 7.338265117214875 7.4230467906574304 7.465437627378708 7.485308332091806 7.494074819465232 7.497727522537493 7.499169379013385 7.499710075191845 7.499903180969866 7.499969012485101 7.4999904792835475 7.499997187658062 7.499999200170416 7.499999780702826 7.499999941961828 7.4999999851562045 7.499999996327164 7.499999999119903 7.499999999795566 7.499999999953925 7.499999999989916 7.499999999997855 7.499999999999556 7.49999999999991 7.499999999999982 7.4999999999999964 7.499999999999999 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 ...)

;; Neat huh?

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

;; The philosophical foundations of probability theory are completely incoherent
;; as far as I can tell, so I've always felt free to believe whatever I like as
;; long as it means my intuition leads me to the right answers.

;; I've found that for me, intuition is best served by imagining that we have a
;; vast number of parallel universes, and in them the results of the experiments
;; are spread out in the proportions given by the formula.

;; It's fair to say that this position is controversial.

;; If we take the numbers to three significant figures:

(clojure.pprint/cl-format nil "~{~3$~^, ~}" (take 10 (poisson 1)))
;; "0.368, 0.368, 0.184, 0.061, 0.015, 0.003, 0.001, 0.000, 0.000, 0.000"

;; Then we'll notice that there are some very rare universes where we get five
;; clicks, or even six from a process where the average number of clicks/second
;; is one.

;; In fact none of the numbers is really zero. There will be some universes
;; where the whole radioactive sample decays at once.  It's probably best not to
;; be near the sample in those universes, but fortunately they're very rare.

(clojure.pprint/cl-format nil "~{~6$~^, ~}" (take 0 (poisson 1)))
;; "0.367879, 0.367879, 0.183940, 0.061313, 0.015328, 0.003066, 0.000511, 0.000073, 0.000009, 0.000001"

;; It turns out we need to sift through a million universes to have a good
;; chance of finding one where nine decays happen in one second. It would be a
;; cold day in hell before you found one where there were twenty.

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

;; Imagine that we're pointing a photon detector at a distant star, and we want
;; to make a guess at how bright it is.

;; Photons happening to hit your detector is very much the same sort of problem
;; as radioactive decays hitting your geiger counter.

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

;; Let's imagine that we've got an absolute stack of parallel universes.  In
;; some of them, the constant for the star is 1, in some it's 2, etc, right up
;; to 20, and those constants are all equally likely

;; We can calculate the proportions of those universes where we see all the
;; different numbers of photons.

;; Photons across the top, 0-9. Star types down the side, 1-20.

(clojure.pprint/cl-format nil "~{~{~$~^,~}~% ~}" (map #(take 10 (poisson %)) (range 1 21)))
;; 0.37,0.37,0.18,0.06,0.02,0.00,0.00,0.00,0.00,0.00
;; 0.14,0.27,0.27,0.18,0.09,0.04,0.01,0.00,0.00,0.00
;; 0.05,0.15,0.22,0.22,0.17,0.10,0.05,0.02,0.01,0.00
;; 0.02,0.07,0.15,0.20,0.20,0.16,0.10,0.06,0.03,0.01
;; 0.01,0.03,0.08,0.14,0.18,0.18,0.15,0.10,0.07,0.04
;; 0.00,0.01,0.04,0.09,0.13,0.16,0.16,0.14,0.10,0.07
;; 0.00,0.01,0.02,0.05,0.09,0.13,0.15,0.15,0.13,0.10
;; 0.00,0.00,0.01,0.03,0.06,0.09,0.12,0.14,0.14,0.12
;; 0.00,0.00,0.00,0.01,0.03,0.06,0.09,0.12,0.13,0.13
;; 0.00,0.00,0.00,0.01,0.02,0.04,0.06,0.09,0.11,0.13
;; 0.00,0.00,0.00,0.00,0.01,0.02,0.04,0.06,0.09,0.11
;; 0.00,0.00,0.00,0.00,0.01,0.01,0.03,0.04,0.07,0.09
;; 0.00,0.00,0.00,0.00,0.00,0.01,0.02,0.03,0.05,0.07
;; 0.00,0.00,0.00,0.00,0.00,0.00,0.01,0.02,0.03,0.05
;; 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.01,0.02,0.03
;; 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.01,0.01,0.02
;; 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.01,0.01
;; 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.01
;; 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00


;; The way I read this table is:

;; We had a look in several hundred parallel universes where we did the experiment.

;; In thirty seven of them, the star constant was one, and we saw no photons in our detector.

;; In another thirty seven, the star constant was one, but we saw one photon.

;; There were four universes where the star constant was five, and we saw nine photons.

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

;; So what should we believe if we point a photodetector at a star and detect 9
;; photons in one second?

;; Well, what we're trying to do is work out what the star constant is.

;; And what we know is that we're in one of the universes where we counted nine photons.

;; I.e. we know that we're in the last row of the table.

;; Let's have a look at that last row:

(clojure.pprint/cl-format nil "~{~$~^,~}" (map #(nth (poisson %) 9) (range 1 21)))

;; "0.00,0.00,0.00,0.01,0.04,0.07,0.10,0.12,0.13,0.13,0.11,0.09,0.07,0.05,0.03,0.02,0.01,0.01,0.00,0.00"

;; What is the total number of the universes that we could be in?
(reduce + (map #(nth (poisson %) 9) (range 1 21)))

;; 0.9963258860986628

;; Call it 100

;; Counting along, we see that there's about a 13% chance that we're in a
;; universe where the constant is 9

;; And about a 13% chance that the constant is 10

;; And a one percent chance that the constant is 4 and that was a particularly lucky second to measure.

;; And a one percent chance that the constant is 18 and it was a particularly bad second to measure in.

;; We might even want to say that there's an 88% chance that the constant is between 6 and 14.

;; Or a half chance that it's 8, 9, 10, or 11

;; As you might have guessed in advance, nine, the number of photons that we
;; actually saw in a second, is probably our best guess at the average number of
;; photons we'll see per second from now on. In classical statistics we'd call
;; that the maximum likelihood estimator, or an unbiased estimator, but in fact
;; there's only a 13% chance that it's the right answer.

;; Probably better just to keep the distribution above in mind.

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

;; What's the expected number of photons per second we'll see if we now run our
;; counter for a while?

;; Remember that we're assuming that stars come in discrete varieties 1-20
;; and that we saw 9 photons.

(defn expected-rate [ stars photons ]
  (/ 
   (reduce + (map *
                  stars
                  (map #(nth (poisson %) photons) stars)))
   (reduce + (map #(nth (poisson %) photons) stars))))

(expected-rate (map inc (range 20)) 9) ; 9.955123984139409

;; This strikes me as very weird.

;; We've seen 9 photons. Surely, given no more information than that, our best
;; guess should be that we'll see 9 photons/second?

;; Consider some edge cases.

;; What about if the stars are only ever of brightness one?

(expected-rate '(1) 9) ; 1.0

;; That's definitely right. We must have just been really lucky to catch nine
;; photons at once. Since we're *sure* that the star is brightness one, the long
;; term estimate has to be one.

;; Here we're in complete agreement with the maximum likelihood estimator, and
;; disagree violently with the unbiased estimator.

;; What about if stars come in brightnesses 1 or 10?

(expected-rate '(1 10) 9) ; 9.999927072835586

;; Again, fair enough. Chances are it was a brightness 10 star behaving
;; normally, but there's a very small chance that it was a brightness 1 star
;; having a very good second.

;; If we run our experiment often enough, the average value will be pulled down
;; a tiny bit by those cases.

;; What if we saw 5 photons?

(expected-rate '(1 10) 5) ; 9.325386911283335

;; Pretty much the same. But now we have to take more seriously the chance now
;; that it was a weak star having a field day rather than a bright star having a
;; half-hearted go.

;; If you're having to budget for photons for some reason, then 9.32 is an
;; unbelievably bad answer. Your budget is going to need to be 10 or 1 exactly.

;; But if you're placing a bet at a bookmaker's, or buying photon insurance from
;; a company with a presence in lots of parallel universes, then 9.33 might well
;; be a fair price. Everyone's got to make a profit.

(expected-rate '(1 10) 2) ; 1.1097148350481423

;; Same sort of thing. About one time in a hundred, two photons will have been
;; emitted by a brightness 10 star rather than a brightness 1 star. However that
;; occasional much brighter star will skew the results noticeably.

;; These answers all look about right to me:
(map (partial expected-rate '(1 10)) (range 20)) ; (1.001110551183876 1.0110931922809403 1.1097148350481423 1.9886759335192858 5.97152862794323 9.325386911283335 9.927658434429095 9.992713129077847 9.999270781535703 9.999927072835586 9.999992707230374 9.999999270722506 9.999999927072245 9.999999992707224 9.999999999270722 9.999999999927073 9.999999999992706 9.999999999999272 9.999999999999927 9.999999999999993)

;; The dodgy case is when we see 4 photons.  Brightness one stars have about a
;; 2% chance of producing 4 photons, and brightness ten stars have about a 2%
;; chance as well.

;; So half the time, we've got a 1, and half the time, we've got a 10.

;; So we're expecting to get a long term average of 10 half the time, and 1 the
;; other half.  And therefore we expect to see around 5.5 photons/second if we
;; add up all our experiments.

;; But I still find this weird:
(expected-rate (map inc (range 20)) 9) ; 9.955123984139409

;; However the whole range of photons that we can receive behaves correctly:

(map (partial expected-rate (map inc (range 20))) (range 40) ) ; (1.5819766656462542 2.163952866521654 3.0148662995404725 3.9961819317191947 4.9990309464433125 5.9993253251689005 6.997434541967013 7.992370962274987 8.980384871785288 9.955123984139409 10.90739733739242 11.825629164415346 12.697213993063654 13.510414954047969 14.256168445529575 14.92921796885468 15.52834070305211 16.055807374000462 16.516421661779617 16.916484374585746 17.262908595847506 17.562577608383947 17.821943140115458 18.04681608289834 18.242291371716 18.412756086521956 18.56194298618969 18.69300434170798 18.808590810798304 18.910926944365613 19.00187928605607 19.08301562443344 19.15565538916098 19.22091189359721 19.27972741592878 19.332902168984916 19.381118149636865 19.424958748228278 19.464924873069005 19.501448223241248)

;; If we get no photons at all, then we've probably got a low brightness star.
;; It'll most likely be a one, might well be a two, could be a three, and so on.
;; So we know the mean *must* be bigger than one.

;; Lots of photons is telling us that we've probably got a twenty, but might
;; have a 19, and could have an 18, so the mean should be less than 20.

;; That's enough evidence of correctness that I believe my calculation, even
;; though at first I was pretty sure that I'd made an off-by-one error
;; somewhere.

;; If anyone's still reading, can they see whether I got it right or not?

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

;; This way of calculating things is called Bayesian Statistics, after the
;; Reverend Thomas Bayes.  It's very general, but was considered to be just a
;; curiosity until computers came along and made it possible to do these
;; calculations on real problems.

;; Poisson processes are a nice simple case that I happened to be thinking about
;; as a result of reading David Mackay's beautiful book Information Theory,
;; Inference, and Learning Algorithms. The nice thing is that the techniques
;; work for any distributions or sets of models at all.

;; This half-witted way of thinking about what the calculations mean is my own,
;; and I don't want to blame it on Bayes or anyone else. But it works for me.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
;; This is a more conventional way to calculate the chance of a Poisson
;; distribution with constant 9 producing 7 events:

(/
 (* (Math/exp -9)
    (reduce * (repeat 7 9)))
 (reduce * (map inc (range 7)))) ; 0.11711612445290907

;; But it produces the same answer as mine.

(nth (poisson 9) 7) ; 0.11711612445290907

;; Proof by it totally works on this example.

Favourite keys for EMACS and SLIME

2011-04-06T21:01:00.002+01:00

I made this post about two years ago. It's been really useful to me, and over the years I've added more commands. These are all the things that make emacs a joy and a horror for clojure editing for me.

Emacs is weird. I often hate it for how much it distracts me from what I'm actually trying to do. And yet when I try to use any other editor, I find myself missing emacs features so much that I end up going back to it. And the more you learn about it, the more addicted you get.

I think that there are fundamental usability mistakes in emacs, mostly to do with its habit of rearranging your windows. It was a great relief when I found out that you can store the window arrangement in a register, and then create a key that says 'put it all back! (you bastard)'. I dream of an editor as powerful as emacs but with a sanely designed UI that obeys modern conventions.

In a very real sense, an editor that is to emacs as clojure is to common lisp. What would be wonderful is if it lived in clojure's standard library like idle does for python.

On the other hand, having got the hang of emacs for clojure editing, it very rarely distracts me these days, and in that mode it's a pure joy. I wouldn't use any other tool. It mainly drives me up the wall now when I'm trying to use it to program in some other language.

(Latest: It was an even greater relief when I found winner-mode, which gives you a key for 'put them back the way they were just now, you stupid bastard'. In fact it gives you an undo history. I've been using that for a while now, and it's great.

I do wonder if the very fact that you can eventually hack the damned thing around so much that it stops being annoying is what has prevented people from starting the programmable editor project again from scratch but starting with something a bit more modern. )

the Editor for Middle Aged Computer Scientists

Is a very good way of editing Clojure, once you have slime and swank installed.

However, the key combinations don't stay in the brain very well if you don't use it all the time, so I thought I'd write down my favourite keys for Clojure editing whilst I'm using it every day, to serve as a cheat sheet for others, and for me in a few months time.

General Editing Keys

C-h m
describe the mode you're in

C-/
Incomprehensible undo

C-w    cut
M-w   copy
C-y     paste
M-y    paste the next most recent thing instead

C-x 0      kill this window
C-x 1      make this window the only window

C-x 2
C-x 3 split the window (vertically and horizontally)

C-x o         go to other window

M-g M-g    go to a particular line by number

C-x h mark whole buffer

M-q reformat paragraph. good for long comments

M-<           beginning
M->           and end of file

C-x C-+      adjust the font size (nice repeating behaviour, try C-x C-+ C-+)
C-x C- -

There's a nice list of useful keys at
http://www.math.uh.edu/~bgb/emacs_keys.html

And this is a nice guide
http://www.gnu.org/software/emacs/tour/

Keyboard Macros

C-x (          start recording
C-x )          stop

C-x e                      execute (then e to reexecute immediately)
C-u 10 C-x e          execute 10 times

C-x C-k C-p move backwards
C-x C-k C-n and forwards in the macro ring

Also see: http://www.emacswiki.org/emacs/KeyboardMacros
and: http://www.gnu.org/software/emacs/manual/html_node/emacs/Keyboard-Macros.html#Keyboard-Macros

Search

C-s forward and reverse incremental search
C-r

Afterwards RET drops you in the buffer, then C-x C-x should get you back because mark was set at the start of the search!

Slimey stuff

M-x slime                         start slime and a clojure instance
M-x slime-connect           start slime, connect to a running swank server

<F12> r find the slime REP
<F12> j    find the clojure buffers
(well, that's what it means once you've made this hack in your .emacs file)


;;bind the slime selector to f12 and add a method for finding clojure buffers
(define-key global-map (kbd "<f12>") 'slime-selector)
(def-slime-selector-method ?j
  "most recently visited clojure-mode buffer."
  (slime-recently-visited-buffer 'clojure-mode))

Completion of Symbols

C-c TAB slime-complete-symbol
M-/

Switching Namespace
C-c M-p

At the REPL

M-p M-n   previous and next in history

C-x C-o flush the last output
C-x M-o    flush the whole REPL buffer

Lisp editing and evaluating keys

Note that in Gnome, Ctrl-Alt-D is 'show desktop', Ctrl-Alt-L is'lock', and so on.
I turn all those off so that I can edit LISP codeeasily.
Use System/Preferences/Keyboard Shortcuts)

C-M-x     evaluate this top-level expression

C-x C-e evaluate the expression behind the cursor

C-u C-x C-e evaluate the expression and paste the result into the buffer!

C-c C-k save compile and load this buffer

M-.
M-,
go to the definition of the expression under the cursor
come back from it

C-c C-d C-d
describe symbol under cursor

C-M-k
kill by s-expression

C-M-h
mark this top-level s-expression

C-M-@
mark next s-expression

(these last three repeat nicely. C-M-h C-M-h marks the next two top-level expressions)

C-M-q
reindent the s-expression

C-M-t
transpose s-expressions (weirdly useful)
repeated presses move an s-expression along inside a larger one.

C-M-a
C-M-e
beginning and end of top level s-expressions (defn ....)

C-M-f
C-M-b
C-M-p
C-M-n
backwards and forwards by list, previous and next s-expression
(the difference is subtle and not terribly important)

C-M-u
C-M-d
up and down list nest levels

Desktops and Windows

a desktop is a list of buffers, places in buffers, etc.

M-: (desktop-read "~/data/dir") will read a desktop file from an arbitrary place

M-x desktop-save will put one down

window arrangements (the layout of the screen) can be saved in registers

C-x rfa store window layout ( in register a , you can use other registers too)

C-x rja restore it

Other commands

M-x htmlize-file is great for saving files with their syntax highlighting in html.

M-x kill-rectangle, yank-rectangle and string-rectangle

Not an emacs command, but I have Alt-F11 bound to full-screen windows in GNOME.

A useful keyboard macro

In clojure-mode, trying to make a keyboard macro to paste results into the buffer
always seems to go horribly wrong. But this particular sequence works to write the result of an expression into the buffer above the expression itself. Place the cursor at the start of the expression.

If you do it by hand, the result ends up under the expression. If you record it and play it back, the result is over. No idea why, but this is a lot less haywire than every other time I've tried to do this.

M-C-f ;; forward-sexp
C-u C-x C-e ;; slime-eval-last-expression
M-C-b ;; backward-sexp
RET ;; reindent-then-newline-and-indent
<up> ;; previous-line

Bookmarks

C-x r m make bookmark
C-x r b jump to bookmark

Hello Web : Dynamic Compojure Web Development at the REPL

2011-03-23T01:11:00.001Z

;; In an article which convinced many people to have a look at LISP

;; http://www.paulgraham.com/avg.html

;; Paul Graham wrote about how nice it was to write a web application from an
;; environment where everything was controllable from the REPL.

;; Here's how to create a web server at the REPL


;; First we'll need a function to generate a simple web page
(defn yo [] "Hello")

;; Now we attach that to /
(use 'compojure.core)
(defroutes main-routes (GET "/" [] (yo)))

;; We should wrap our routing table in a handler (note the #')
(use 'compojure.handler)
(def app (site #'main-routes))

;; And then we can run it with jetty (again note the #')
(use 'ring.adapter.jetty)
(def server (run-jetty #'app {:port 8080 :join? false}))

;; So now we're serving web pages. Go look at http://localhost:8080/

;; The web browser of the gods:
;; $ watch -d -n 1 curl -sv http://localhost:8080/ 

;; We can stop the server
(.stop server)

;; And restart it:
(.start server)

;; We can redefine those functions, and the change takes effect immediately:
(defn yo [] "<h1>Hello<h1/>")

;; (you will need to refresh the page, unless you use twbotg)

;; HTML is a bit of a pain. We already have a syntax for tree-structured data:

(use 'hiccup.page-helpers)

(defn yo []      
  (html5 [:head
          [:title "Hello World"]]
         [:body
          [:h1 "Lisp is the future"]]))


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

;; I'm confused: see http://stackoverflow.com/questions/5398569/can-clojure-be-made-dynamic

;; The reason for the #' above is to make sure that we can redefine main-routes
;; and app at runtime. Without the indirection, the changes don't seem to get picked up.

;; And yet they do get picked up when you redefine yo.

;; I do not understand what is going on here. Any explanation would be welcome!

To start a REPL with the necessary libraries available, you can use maven.

To install maven on ubuntu:
$ sudo apt-get install maven2

Put this xml in a file named pom.xml, and then type
$ mvn clojure:repl
in the same directory.For a swank server that you can use with emacs, use
$ mvn clojure:swank

The first time you run it it will do an incredible amount of downloading to get all the necessary bits, but after that it's fairly speedy.

<project>
  <modelVersion>4.0.0</modelVersion>
  <groupId>com.aspden</groupId>
  <artifactId>hello-compojure</artifactId>
  <version>0.1-SNAPSHOT</version>
  <name>hello-compojure</name>
  <description>maven, clojure, emacs, swank, compojure: together at last</description>
  <url>http://www.learningclojure.com</url>

  <!-- repositories -->
  <repositories>
    <repository>
      <id>clojars</id>
      <url>http://clojars.org/repo/</url>
    </repository>
    <repository>
      <id>clojure</id>
      <url>http://build.clojure.org/releases</url>
    </repository>
    <repository>
      <id>central</id>
      <url>http://repo1.maven.org/maven2</url>
    </repository>
  </repositories>

  <!-- libraries -->
  <dependencies>
    <dependency>
      <groupId>compojure</groupId>
      <artifactId>compojure</artifactId>
      <version>0.6.2</version>
    </dependency>
    <dependency>
      <groupId>hiccup</groupId>
      <artifactId>hiccup</artifactId>
      <version>0.3.4</version>
    </dependency>
    <dependency>
      <groupId>ring</groupId>
      <artifactId>ring-jetty-adapter</artifactId>
      <version>0.3.7</version>
    </dependency>
    <dependency>
      <groupId>org.clojure</groupId>
      <artifactId>clojure</artifactId>
      <version>1.2.0</version>
    </dependency>
    <dependency>
      <groupId>org.clojure</groupId>
      <artifactId>clojure-contrib</artifactId>
      <version>1.2.0</version>
    </dependency>
    <dependency>
      <groupId>jline</groupId>
      <artifactId>jline</artifactId>
      <version>0.9.94</version>
    </dependency>
    <dependency>
      <groupId>swank-clojure</groupId>
      <artifactId>swank-clojure</artifactId>
      <version>1.2.1</version>
    </dependency>
  </dependencies>

  <build>
    <plugins>
      <!-- talios' clojure-maven-plugin provides mvn clojure:swank etc -->
      <plugin>
        <groupId>com.theoryinpractise</groupId>
        <artifactId>clojure-maven-plugin</artifactId>
        <version>1.3.2</version>
      </plugin>
      <!-- The versions plugin allows you to find out what you can upgrade -->
      <!-- mvn versions:help -->
      <!-- mvn versions:display-dependency-updates -->
      <plugin>
        <groupId>org.codehaus.mojo</groupId>
        <artifactId>versions-maven-plugin</artifactId>
      </plugin>
    </plugins>
  </build>
  

</project>

Advice to a Newbie

2011-02-23T16:50:00.002Z

Someone wrote to me the other day asking me to teach him Clojure over Skype.

Well, I am sceptical of the feasibility of this scheme, but it is well known that the words "£50/hour" ring pleasantly in my ears, although in fact for this kind of one-hour every so often-type piecework it might be truer to say that the real figure is rather higher, but what the hell, this guy is a beginner interested in lisp, so I can cut him some slack....

However it occurs to me that he doesn't really need me, and as a very wise man once pointed out to me, my clients pay me for my advice, rather than for my typing.

So this is my reply:

I have removed any identifying details from it, I hope, but if the guy in question wishes to out himself here, then that is perfectly fine, and maybe he will find a study partner who lives somewhere in his area.

And if anyone reading this would like to recommend their own route in, or just flame me for being so utterly ignorant and old-fashioned, then feel free to leave long comments, at the risk that I may use them as material for further posts.

Dear .......,

I'm sorry I haven't replied for so long. I've been rushed off my feet, but I was very flattered by your offer.

Clojure's exceedingly cool, but to understand it you need to speak lisp. (Scheme is the purest and easiest lisp to understand).

I'd recommend the route I took into all that sort of thing, the excellent series of lectures given to first years at MIT, and filmed here as given to a load of HP execs by the original lecturers.

http://groups.csail.mit.edu/mac/classes/6.001/abelson-sussman-lectures/

You should get the excellent PLT Racket (which used to be called PLT Scheme),

http://racket-lang.org/

which has a superb beginners' scheme editor and environment, called DrScheme (or these days I presume DrRacket)
and set the required language to Scheme R4RS or R5RS (not modern racket, which is a similar but different language to the one used in the book).

That way you won't have to deal with the incomprehensibly weird emacs at the same time as learning lisp.

When they write anything on the board in the videos, stop the lecture and try it out in DrRacket.

There's also the excellent companion textbook:

http://mitpress.mit.edu/sicp/full-text/book/book.html

which has been used as an introductory CompSci text for many years.

I seriously recommend you to do *all* the exercises in the book. Even the ones that look boring are in fact spectacularly well chosen and interesting once you start, and each one will teach you something new.

Until you've done them all, you won't understand the material well enough to be able to go onto the next chapter.

You will of course ignore this advice. I know I would.

When you are half way through Chapter 2, and a bit bewildered and wondering if you are not clever enough, or why anyone should be interested in this weird and incomprehensible way of doing simple things, remember the advice and go back and do all the exercises in Chapter 1.

Be aware that this is a very long process. Reckon on being able to do one exercise per day, and if you're treating it as a hobby, reckon on taking a year to understand the whole book.

But it's also terrific fun, and you'll get regular rewards of small doses of enlightenment throughout.

I actually did all this. I've never been to any sort of computer training. I learned lisp from that wonderful book, those wonderful lectures, and DrScheme.

It really is the easiest way in. That's where all the smuggest lisp weenies come from, or go eventually.

And once you've learned scheme (or more accurately, once you've learned all the new ways of thinking that Abelson and Sussman demonstrate), you'll look at clojure and think:

'Oh yes, that's nice... I see...'

-----

That said, I am a complete slut, and amongst the things I am prepared to do for £50 an hour is to try to teach someone Clojure over Skype.

But I really think that you'll get on better with the video lectures.

Also, even though I am rather busy at the moment, I am not entirely without public spiritedness.

If you read the first chapter of SICP, and watch the associated lectures, and have an honest go at doing all the exercises in chapter one, then you can skype me and I'll spend a couple of hours leading you through any exercises you had difficulty with for free.

Cheers, John.

P.S. I find that I have just written a fine blog post, which I will publish. I will remove your name.

P.P.S. Hypocrite note. I never did get round to chapter 5 of SICP. I intend to one day.

P.P.P.S. In fact, you can probably get away with just reading the first three before you can understand clojure. But chapter 4 is definitely the best of the first four, so it would be a shame to stop there.
And you won't understand the why of lisp until you've read chapter 4.

Clojure Dojo 4: Symbolic Differentiation

2011-02-16T16:45:00.000Z

We've written a fairly good Newton-Raphson solver, but it's a bit hacky to do numerical differentiation when differentiation is so easy by hand.

Here I've tried to create the simplest possible symbolic differentiator.

It relies on four basic rules:

The derivative of a constant is 0.
The derivative of a variable is 1 with respect to itself, 0 with respect to a different variable.
The derivative of a sum is the sum of the derivatives of the parts.
The derivative of a product is the sum of the derivative of the first part times the second part and the first part times the derivative of the second part.

If you tell the computer this then it can differentiate polynomials like
x^5+x^2+2x+3 as long as you express them in terms of binary operators.

(+ (+ (* x (* x x)) (* 2 (* x x))) 1)

It's easily extended to deal with arbitrary operators like (* x y z 2), and with functions like sin and cos.

A routine to simplify the answers would also be nice, so that (+ (* 0 (* (* x x) x)) (* x 1)) might be more readably represented as x.

;; The simplest possible symbolic differentiator

;; Functions to create and unpack additions like (+ 1 2)
(defn make-add [ a b ] (list '+ a b))
(defn addition? [x] (and (=(count x) 3) (= (first x) '+)))
(defn add1   [x] (second x))
(defn add2   [x] (second (rest x)))


;; Similar for multiplications (* 1 2)
(defn make-mul [ a b ] (list '* a b))
(defn multiplication? [x] (and (=(count x) 3) (= (first x) '*)))
(defn mul1   [x] (second x))
(defn mul2   [x] (second (rest x)))


;; Differentiation. 
(defn deriv [exp var]
  (cond (number? exp) 0                                                              ;; d/dx c -> 0
        (symbol? exp) (if (= exp var) 1 0)                                           ;; d/dx x -> 1, d/dx y -> 0

        (addition? exp) (make-add (deriv (add1 exp) var) (deriv (add2 exp) var))     ;; d/dx a+b -> d/dx a + d/dx b
        (multiplication? exp) (make-add (make-mul (deriv (mul1 exp) var) (mul2 exp)) ;; d/dx a*b -> d/dx a * b + a * d/dx b
                                        (make-mul (mul1 exp) (deriv (mul2 exp) var)))
        :else :error))


;;an example of use: create the function x -> x^3 + 2x^2 + 1 and its derivative 
(def poly '(+ (+ (* x (* x x)) (* 2 (* x x))) 1))

(defn poly->fnform [poly] (list 'fn '[x] poly))

(def polyfn  (eval (poly->fnform poly)))
(def dpolyfn (eval (poly->fnform (deriv poly 'x))))



;;tests

(use 'clojure.test)

(deftest deriv-test
  (testing "binary operators"
    (is (= (let [m '(* a b)] [(multiplication? m) (make-mul (mul1 m) (mul2 m))]) [true  '(* a b)]))
    (is (= (let [m '(* a b)] [(addition? m)       (make-add (add1 m) (add2 m))]) [false '(+ a b)])))
  (testing "derivative function"

    (is (= (deriv '0 'x)               '0))
    (is (= (deriv '1 'x)               '0))
    (is (= (deriv 'x 'x)               '1))
    (is (= (deriv 'y 'x)               '0))
    (is (= (deriv '(+ x x) 'x)         '(+ 1 1)))
    (is (= (deriv '(* x x) 'x)         '(+ (* 1 x) (* x 1))))
    (is (= (deriv '(* x x) 'y)         '(+ (* 0 x) (* x 0))))
    (is (= (deriv '(* x (* x x)) 'x)   '(+ (* 1 (* x x)) (* x (+ (* 1 x) (* x 1)))))))
  (testing "function creation: d/dx (x^3 + 2x^2 + 1) = 3x^2 + 4x "
    (let [poly '(+ (+ (* x (* x x)) (* 2 (* x x))) 1)]
      (is (= ((eval (poly->fnform poly)) 3) 46))
      (is (= ((eval (poly->fnform (deriv poly 'x))) 3))))))

Clojure Dojo 3: From Heron to Newton-Raphson

2011-02-16T16:41:00.000Z



;; We've so far written a square root solver using Heron's method

(defn average [a b] 
  (/ (+ a b) 2))

(defn abs[x] 
  (if (< x 0) (- x) x))

(defn make-good-enough? [n]
  (fn [guess] (< (abs (- n (* guess guess))) 1e-6)))

(defn make-improver [n]
  (fn [guess] (average guess (/ n guess))))

(defn iterative-improve [x improve good?]
  (if (good? x) x
      (iterative-improve (improve x) improve good?)))

(defn square-root [n]
     (iterative-improve 1.0 (make-improver n) (make-good-enough? n)))



;; Now in fact, Heron's method turns out to be only a simple example of a more general method of root finding known as Newton-Raphson.
;; Their great idea was to use the derivative of a function to help find its roots.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Finding Derivatives (by cheating)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Suppose we've got a function, say cube

(defn cube [x] (* x x x))


;; And we want to find out what its derivative is, i.e. how much it changes when you change the argument.

(cube 4)        ; i s 64
(cube 4.000001) ; is 64.000048000012

;; In other words, if we add 0.000001 to x, then (cube x) goes up by 0.00004800..... which is about 48 times as much.

;; We say that the derivative of cube at 4 is (about) 48.

;; It's different in different places.


(cube 3)        ; is 27
(cube 3.000001) ; is 27.000027000009005

;; So the rate of change of cube at 3 is (about) 27.

;;In general, we want a function that takes a function, and gives back a function that tells us how much it changes if you add a tiny bit to its argument.

(defn deriv [f]
  (fn [x] (/ 
           (- 
            (f (+ x 0.000001)) 
            (f x))
           0.000001)))


;; Let's try that out:

((deriv cube) 4)       48.00001200067072
((deriv cube) 3)       27.00000900546229

;; Notice what we did here! This is a function which takes a function and gives back a function.

;; We could name the answer!
(def dcube (deriv cube))

(map dcube (range 10))        ;(1.0E-12 3.0000029997978572 12.000006002210739 27.00000900546229 48.00001200067072 75.00001501625775 108.00001800248538 147.00002100198617 192.00002384422987 243.0000268986987)


;; So now we know how to calculate (a good approximation to) the derivative of any function.


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Newton Raphson
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Imagine that we wanted to find the cube root of 100.

;; That's the same as finding a number which, when you cube it and subtract 100, gives you zero


;; So we should take as our equation to solve: (cube x) - 100 = 0

;; or (f x) = 0 where f is 

(defn f[x] (- (* x x x) 100))


;; Isaac Newton told us that if you have an equation like that, and you have a guess at a number
;; which will make it zero, then you can find a much better guess if you know the derivative.


(def df (deriv f))

;; Newton said: Suppose that I guess that 4 is the cube root of 100.

;; Then I try it (f 4) = -36

;; That's not so close. What's the derivative there?

(df 4) 48.00001200067072

;; What that means is that when we raise x by a tiny amount, say 4.0000001, then (f x) goes up by about 48 times that.


;; So Newton tells us:

;; The function is too low by 36.
;; If we make the guess a bit larger, then the function will go up by ~ 48 times as much.
;; We should try adding 36/48 as our next guess.

;; Let's try

(+ 4.0 (/ 36 48)) ;4.75

(f 4.75)          ;7.171875

;; ok, better. We were under by 36, now we're over by 7. 4.75 is too high!

;; Take the derivative again:

(df 4.75) ;67.68751426022845

;;Divide the amount we're off by the derivative
(/ (f 4.75) (df 4.75)) ;0.10595565634789487


;; So our next guess should be (- 4.75 0.10595565634789487) , which is 4.644044343652105

;; Let's try 4.644 :

(f 4.644) ;0.15592198400001678 ...homing in...

(- 4.644 (/ (f 4.644) (df 4.644)))  ;4.641590085793267

(f 4.64159) ;7.538717167676623E-5 ...very good...


;; Obviously we've got another guess-improving function here

(defn improve-cbrt100 [guess]
  (- guess (/ (f guess) (df guess))))


(improve-cbrt100 4.0)     4.749999812489567
(improve-cbrt100 (improve-cbrt100 4.0)) 4.644044335287887

(take 10 (iterate improve-cbrt100 4.0))   ; (4.0 4.749999812489567 4.644044335287887 4.6415901322397985 4.641588833613422 4.641588833612778 4.641588833612778 4.641588833612778 4.641588833612778 4.641588833612778)


;; Let's try 4.641588833612778, which seems to be about as good as floating point arithmetic can get us.

(f 4.641588833612778)      ; -2.8421709430404007E-14

(cube 4.641588833612778)   ; 99.99999999999997

;; Pretty good for only 5 steps!

;; What about a good-enough function to tell us when to stop iterating?


(defn good-enough-cbrt100? [x]
  (< (abs (f x)) 0.0000001))


;; Now we can plug the improver and the good enough function into iterative-improve, as before. We'll use 1.0 as our initial guess again.
(iterative-improve 1.0 improve-cbrt100 good-enough-cbrt100?) ;; 4.641588833613406

(cube 4.641588833613406)   ;; 100.00000000004056



;; But of course, nothing we did depended on (f x) being (- (cube x) 100)
;; We can make a guess-improver for any function



(defn make-improver [f]
  (fn [guess] (- guess (/ (f guess) ((deriv f) guess)))))


;; and a function to tell us whether it's good enough

(defn make-good-enough [f]
  (fn [guess] (< (abs (f guess)) 0.0000001)))

;; What if we'd like to solve the equation x^3 + x^2 + x + 1 = 0 ?


(defn solve [f guess]
  (iterative-improve guess (make-improver f) (make-good-enough f)))

(solve (fn [x] (+ (* x x x) (* x x) x 1)) 1.0) ;  -1.0000000235152005

;; Let's see how good an answer that is.

((fn [x] (+ (* x x x) (* x x) x 1)) -1.0000000235152005) ; -4.703040201725628E-8




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Newton Raphson Solver
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Here's the whole of our general function solver, but let's make the two 'magic numbers' into variables too.

(defn deriv [f dx]
  (fn [x] (/  (- (f (+ x dx)) (f x) )  dx)))

(defn make-improver [f dx]
  (fn [guess] (- guess (/ (f guess) ((deriv f dx) guess)))))

(defn make-good-enough [f tolerance]
  (fn [guess] (< (abs (f guess)) tolerance)))

(defn iterative-improve [x improve good?]
  (if (good? x) x
      (iterative-improve (improve x) improve good?)))

(defn solve [f guess dx tolerance]
  (iterative-improve guess (make-improver f dx) (make-good-enough f tolerance)))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Some applications
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;usually we'll want an initial guess of 1.0, a tolerance of 0.0000001 and a dx of 0.0000001
(defn solve-it [f] 
  (solve f 1.0 0.000001 0.000001))

;; Square roots are the solutions of (- (* x x) n) = 0


(defn sqrt [n] 
  (solve-it 
   (fn [x] (- (* x x) n))))

(sqrt 2) 

;; Cube roots are the solutions of (- (* x x x) n) = 0

(defn cubert [n] 
  (solve-it 
   (fn[x] (- (* x x x) n))))

(cubert 2)


;; General positive integer powers can be expressed recursively:

(defn pow [x m]
  (if (= m 0) 1
      (* x (pow x (- m 1)))))

(pow 2 3) ;8

;; And if we can express them, we can use our method to find mth roots.


(defn mthrt [m n] 
  (solve-it 
   (fn[x] (- (pow x m) n))))

(mthrt 4 10000) ; 10.0
(mthrt 123 234) ; 1.045350466874579

(pow 1.045350466874579 123) ;; 234.0000000243792


;; But of course, we can solve any equation which is amenable to Newton's method.
;; What number, when multiplied by itself 5 times, and added to its cube, is equal to 700?

(solve-it (fn[x] (+ (* x x x x x) (* x x x) (- 700))))         ; 3.653807138770812

((fn[x] (+ (* x x x x x) (* x x x))) 3.653807138770812) ; 700.0000003041337

;; And if we'd like a better answer, we can tighten the tolerances

(solve (fn[x] (+ (* x x x x x) (* x x x) (- 700))) 1.0 0.0001 0.0000000001) ; 3.6538071384442974


((fn[x] (+ (* x x x x x) (* x x x))) 3.6538071384442974) ;  700.0000000000817









;; In theory, we'd be able to use exact arithmetic to get arbitrarily accurate answers.
;; Unfortunately, there's a problem:


(def it (make-improver (fn[x] (+ (* x x x x x) (* x x x) (- 700))) 1/1000000)) ; let's use an exact number for the derivative approximation

;; and start with an exact guess

(it 1)                  ; 706000013000011000005000001/8000013000011000005000001
(it (it 1))             ; 701635321996143769309418549344377227910245191235491382603948489939699061403896688219959595828904782223381164841011935230760200930031084004228000441000030000001/9938320550977024134920249980653903204519665284470123223646489755187942893061019955878983167997749088761731121640594804980166030029688004228000441000030000001

(it (it (it 1)))        ; 680237917441885524336418138840317390387477847859858876595840234272149425820646224687954143292874376454521646334951672503154114880181300843385953617515081285781202235479686701822851428206060076054501437170597117295185786622448092483300621863099620331252416621986363001610323019591307968502301991515214451683404838437583955509995832608492690263029459069929976980341191436712932448855080982518511841862770021155980758185192008036141486759432138585947862423162360655529584841649322416573707879742014567586254212003318865179086757633089843219756812844011861068141391297639719428726706532258163565342026653477095733524189143027946838723155848407864793520497990564324241350635343653518356286709120797176348156180962741126448716721627967605789740379911361208362221973935469409915879343921377079168493057613428011966000155000001/12044286332140129752156903683583126447244938902908926261607579552042467607850219073475448382767207875539137993119265514088232278316554495937166615076210884321449749330293196083213705419806298445939325718843197812663753008605821508903073125871316024961234609794639656768170967982087039663880189123415565181043003755509161080647377737131911463013608843269269085395872146930064532858794789612245264182978492519500766787287053880913766446555513734299742826803443919934807497308923876636662216418140448109767596616697192628702623509472332595967524822599329261416093615390398256303064181652715052013431614546617380753560488219680455031256525883523608904103537508019124887595921116604264464741508409082509232101707795378534258208496424946224563887570504744604783840489200142804159760458106765068490963613428011966000155000001
(it (it (it (it 1))))   ; 582683064888886254627237413128914355133205070107730345306566818667169839077043977063385150567627868898954627396888803873298542915735291068231854737660749871284970969691761755237023440170586132711096078205774661890184691373273510273283202351760392701768411497903999874829223861118399238720882465581354352630000971969184626801423877016182068328118391265021597001152296290812122250720956031547021906861678543752400823983366892551426008467510890046723277161664312028285174276774217375277065253742530859064530016440400316790801769530544371284081707010989839631883556283241113998175234840943170722268328720701199347867588322060049770935643979286527454236006759544836737208537749694757872416290454093930058675287658912516799102545283953465977040744182800487171767210596289815719837983663631950117171893161135122746529146479448859227394441049923190055342857801930221659120398700765917475396756356933748228771216329690971598121024906393909199740001655437226675292089740865401850681717862593295044249444227096530596676969537355119944076825364494231692463504463385356777816358752007019352340833078969061391374565063342370784189339508167803312120387247644380450529051088582048118902223778854166538058459895429188150443034952695968773967230092326545545925098065274390529306120532407597280407465373342806740914557404279803404281875098412568052206856189504731510965158306039056809060766012400914292618679128113015139253847973084750652345634913827827043851034605686643123187076375915598639958363937697829956243245104853978610467265918061529051311071288413580756443693449482971629647374397753310557728689142976074021187362385338932147567642210185904671961035565315840293902065425141448705637965586698203041117408261444705945865047531954486514780505078832025358426644461672866572253641788067682868313500867194340509941792198538459647737009443599755837217723308620719398896201880698984796565219264280137596609642845400547125754456209756196169964659473873586279764214453671645956003619308063592648337612774637171241884690676604002029522916589191903605327335495103912144305652194435971484948567480613108531306049822517644426959900990535402937677569195079220088350105434599994075059763835418492356744894787315782423704801923419774013400657508864679447262609755156219146757585536473633170206959340454211578712959760853352548696411404380901268846622462357056396277653709161131992015278410666198586576834004214031859566710336893361077448333057320403538099709665120727817043295276197596679826722665947871368777705493013687853852488123574246109063974638016285172124627788918704893034469219280568313768497634816004319733836637800496141268248807282891107795078010343054429751145414810367778479411850027960387501665140659520089944949565082934980430105203072555206594274562501052867096198848952152330124216130716659045181152784334507386057959224529758674635647399911089583876402241119812642070687710563501331168402312038757841360479166547376873053789334438351535437225196907596309808445013472858435271020713482531473249425954713138213993853690048279402486411516169465400907740591254918376665681348518024822387460412485718766760961958943780464278973812982576607564820484862661878406223946061213224303753767486290224002094709777464746472645414893385625996401380435392614667079096134453611792569645198318837308794387253480934803497220222891193495355138449277282306273087395040974247699578508944806995861199096983093778899088019827135297023837526798957092191839173362464128894165682017815738236599955940542605196688759536725013516191327597134380086354086023426735965448551045630879574037331690144148829216502025494635660367861837832650027912461998002111883455796634997105423180569440930861025594379890711775754266620157029760796036618547810092105837604672377921726136778076482261042587195382143724289727648562704136017853509498366890447377518135820019350488687720760940855188770925017427141691614329578171795155245502836823003566045473401512439769020488822239922252814521793930372045109677392922267054560823873124298362347797278136358895228352710628394677753104291145730759681035568127201415728926581298819651031969039124011747083909428303966000780000001/12896626453600069493698610007462574074360012803192517550652653321288458704437628147812351238096058570880646809642684891850066999395344051036048084092196565019966845694603226043392692831832408470853928980193581747252243200468103361176860667342346033755639391767792595834026685821301417812313001324714923798171069125692865759809093740132786844460822717151100184874113026243342113737733017773505307764700961069595789996009548181632271209787338405203417470570243624108296656484641803537134016638223204512186999071878554046958682488174902216383193717986143532189831044292308519842126873269313853560876247133978458353574480851312759631684358889032900084642745496288725888743073969677137896558829603592262230802522695739419744389018487062275548598153935413565238867472051597918163575900094811377931457308002870398147828447466466856465628782660404829455110767346686507193439951612181709790999494520905843712219628937059390531846716829881905679762662797797772057202910931972704357372155619580235092749976889122762716133368842037661984150683629394903567236216963194674605967568392431208794043431169892910445094675188563758743114907987788746591274250207082015901133815715910101081963885658315948409184499380715077775262542275170959223244947248063594716133379982128239636514206935394270725149984811164184816315321673225284930879745768846720157626335366782084507443322775922411153536241445645457909068808146547276992161599402164331016798612763537783551540734336332823264126557187267639802943507420912300582217231591699364446110633032192667875702437386209756724502546707178905247996448449229725676664533990896268840252253987954869584921155755185596659051602784734099723461617865205317857691298850939510647703168248071840078142301962876095991821410664592847145496334778994271130316719075179553320153706951867058953414639227099786449755622188371176105537376510675813747465892564129310902573333355338988514058192477421861710533925137494638802746649014987114280003213204406857866359359280991867594433240560958642355778026603456310571266226381230212501928363154073509767737834411081511741666078389988426255700754713921887266615912226559289922110288321637513352794999588449807206908304314944589521175612052878625277657474456489493344970961101196034555416517832906272790448598835563702334203053761460216483188413922273789078424479313993102344574019369085117678979346999182040112416824764512591041037973913657403354399232242224999268943921208792520828958494399634230250208344108593016608457066759423447434130404474515620954899757483291938447043440289286090559545316946393298710577732756942628616513884579585056022530728783711122920500882377115985715201353863643271746094708630198929184352363067785421438134909914864758356500163689964428638647687495591442052889991128729671129301459039511338151733777970655351650119214267552140948065933092934877227143743916767768243644993862805442252649852133799824872693054837337674268061282267623038576543034125243858631285899198303109966963172209152874196720199801119238717086988413507996967711063340008844027377086930137204987454203538469697377749646520851213248589905149021586496254737543585289159013650128559023578482611563124559893592352427180734048684239663086639329223310790701294430611951452851943347414551549066596198403402509218874978312155181397196069302152786264584396523101878569322693955950044160266238264131771759791930583114290189937813208245285469977694328792320835576523730805116731577483465467653183154605688937756438221438317995677209102959548247744261605934633001969642392983652820479907300076279299375999438660224810860948463249531812666274785951673129729793010427319485797202888099118900435857053255162204718549364268309991788374329365868674979202327345545939405245168318555414785299112969636869024689405450765132794123797478424267092928772922594681937404357022580366085874534312247818159047202717252632019175403527618827567476810823460590807535175349230955481084676292524894761199878853286342204470938790343375932066863536221551379561495530913704293344747209001634639584668143339854848154541379814946625988808902291845280350306322621482114145676960211744291909428303966000780000001

;; after 4 iterations we've got 8000 digits in our fraction! It might take a while to calculate each step.

;; It's possible to control this crazed exponentiation of digits. 
;; I leave the problem of modifying the method to produce arbitrarily accurate answers as a problem to the reader.

Clojure Dojo 2: What about numbers that aren't 10?

2011-02-16T16:40:00.000Z



;; We've so far implemented Heron's method for finding the square root of
;; the number 10. The program looks like this:

(defn average [a b] 
  (/ (+ a b) 2))

(defn improve-guess [guess]
  (average guess (/ 10 guess)))

(defn abs[x] 
  (if (< x 0) (- x) x))

(defn good-enough? [guess] 
  (< (abs (- 10 (* guess guess))) 1e-6))

(defn good-enough-guess [x]
  (if (good-enough? x) x
      (good-enough-guess (improve-guess x))))


;; At this point, it might be a good idea to put the functions in a file. 
;; If we use a lisp-conscious IDE (EMACS, the Editor for Middle Aged Computer Scientists, is my favourite)
;; Then we can carry on our conversational style of development.

;; In fact, this program is easily generalized. Our first task is to make a more general
;; iterative-improve function.  This function should take a guess, a function to improve guesses,
;; and a function to tell when guesses are good enough, and return an answer which is good enough,
;; which it should find by repeatedly improving the guesses.

;; We almost have this function already. We just need to make the constants in good-enough-guess into arguments.


(defn iterative-improve [x improve good?]
  (if (good? x) x
      (iterative-improve (improve x) improve good?)))

;; Let's test that: (In emacs, put the cursor at the end of the expression and type C-x C-e to evaluate the expression)

(iterative-improve 1.0 improve-guess good-enough?)

;; Or in fact, we can use C-u C-x C-e to evaluate it and paste it into the buffer.


(iterative-improve 1.0 improve-guess good-enough?)
3.162277665175675

;; (* 3.162277665175675 3.162277665175675) is 10.000000031668918, so we haven't broken our method.

;; Of course to find the square roots of numbers other than 10, we need functions that will make our guesses closer and tell whether the answers are good enough.

;; It would be a terrible thing to have to hand code them every time. But we can make functions which make functions:

(defn make-improver [n]
  (fn [guess] (average guess (/ n guess))))

(def f (make-improver 25))

(f 1.0)                                 ;13.0

(f (f 1.0))                             ;7.461538461538462
(take 10 (iterate f 1.0))               ;(1.0 13.0 7.461538461538462 5.406026962727994 5.015247601944898 5.000023178253949 5.000000000053722 5.0 5.0 5.0)

(defn make-good-enough? [n]
  (fn [guess] (< (abs (- n (* guess guess))) 1e-6)))

(def g? (make-good-enough? 25))

(take 10 (map g? (iterate f 1.0)))      ;(false false false false false false true true true true)


(iterative-improve 1.0 f g?)            ;5.000000000053722

;; So now we can find arbitrary square roots

(defn square-root [n]
     (iterative-improve 1.0 (make-improver n) (make-good-enough? n)))

(square-root 2)                         ;1.4142135623746899

Clojure Dojo: The method of Heron of Alexandria

2011-02-16T16:39:00.002Z

To my further surprise, I was asked to give this talk again this year, only twice. So I've moved this back up to the top in case anyone who was there wants to read through the original draft.

Somewhat to my surprise, I've been asked to give a talk on Clojure in London in a few weeks. http://wiki.2010.dev8d.org/w/Clojure_Lab

It will be for an audience of people who don't have LISP experience, and it has to be in the form of a 'programming dojo', where people take turns to write real code.

I decided that talking about all the things that make Clojure the best new LISP for 30 years would maybe be a bit redundant with such an audience, so instead I thought I'd try and go through the example that got me hooked on LISP, the generalized iterative improver from Structure and Interpretation of Computer Programs.

Here's a draft of the first half. I don't actually expect to get any further than this in two hours,but just in case, there's also a second half that will turn this into a Newton-Raphson solver, and then go on to symbolic differentiation of functions.

;; When we write lisp, we have magic powers, and the reason is that our code is in a form that is
;; easy to manipulate programmatically.

;; Like all magic, there is a price to pay.

;; The price is that we have to rip the front end off our compiler.

;; The front end is the thing that takes 1*2+3*5 and says:

;; 1 is valid but 1* isn't, * is an operator, and so on, so that must really be 1 * 2 + 3 * 5 , and
;; * binds tighter than +, so that's really ( ( 1 * 2 ) + ( 3 * 5 ) ) Which is really a tree:


             (+ 
              (* 5 3)
              (* 1 2))

;; Meaning take 1 and 2 and multiply them, take 5 and 3 and multiply them. Add the results.

;; This is clearly seriously annoying, and like monads in Haskell, it's the first thing you hit when
;; you learn lisp, but I promise that after a month or so of using it you stop noticing it, and
;; although it's never quite as good for actual arithmetic, it's actually much nicer as a notation
;; for a generalized function call.

;; More importantly, it allows us to treat all functions uniformly, including the ones we define

;; ourselves. And this is the source of the magic.

;; If we do well today we'll be able to symbolically differentiate a function.  Using a very short
;; program. And that's hard in a language with a syntax.

;; So my aim for today is to get us used to the (lack of) syntax of lisp.

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

;; We're going to use clojure, which is a modern member of the LISP family, in the same sense that
;; Java is a member of the ALGOL family.


;; Clojure's a wonderful language, which uses Java bytecode as its machine language, and runs on the
;; JVM, and so has easy natural access to all the libraries that Java has.

;; My friends at the Chemistry department in Cambridge are using it because they've written a lot of
;; Java over the years, and while they like it for all sorts of reasons, they've found that Clojure
;; is an easier way to write Java than Java is.

;; It has pervasive laziness throughout the language, which allows us to disconnect the things that
;; need to be done from the order they need to be done in.


;; It has extraordinary pure-functional data structures, and baked-in software transactional memory,
;; which together make for a style of programming that can run in parallel on many cores.

;; It has maps and sets and vectors and regular expressions built right into the language, and so
;; natural to use that they're everyday objects like they are for PERL and Python people.

;; But I'm going to ignore all of that, and we're going to get over the pons asinorum of lisp, which
;; is all those damned brackets.

;; We're going to write a program that you could have written in LISP forty years ago, when the

;; newest member of the ALGOL family was ALGOL itself.


;; -------------------------------------------------------------------------------------------------

;; The method of Heron of Alexandria
;; ---------------------------------

;; I don't know if any of you remember not having a calculator.  Once upon a time they didn't even
;; have mathematical tables, and so they had to calculate things like square roots by hand.


;; One popular method of calculating a square root is due to Hero(n) of Alexandria, who also
;; invented steam power, the windmill, the syringe and the vending machine.

;; What Heron said to do was this:

;; Suppose you have a number that you want to find the square root of, let's say 10.

;; And suppose you have a guess for where that square root might be, let's say 1.  It's not a very
;; good guess, because the square of 1 is 1, not ten, but it will do to get started.

;; So Heron tells us, if we know the square root, then if we divide 10 by it, we'll get it back.

;; Like say the square root of 4 is 2, and if we divide 4 by 2, then we get 2.  And that's what it
;; means to *be* a square root.

;; But if we guessed too low, then the thing we get back will be too high.  And if it's too high,
;; then the thing that we get back will be too low.

;; So Heron says to take the average of what we have and what we get back when we do this division,
;; and he promises us that will be a better guess.

;; Let's try that, for our problem number 10, and our guess 1

user> (/ 10 1)
10

;; We guessed 1, we divided 10 by 1, we got back 10.
;; What's the average?
user> (/ (+ 10 1) 2)
11/2 

;; Or if we do the whole calculation at once, it looks like:
user> (/ (+ (/ 10 1) 1) 2)
11/2

;; That's getting a bit hard to read, so we should define a function to give us averages:

user> (defn average [a b] (/ (+ a b) 2))
#'user/average

;; We can test it:
user> (average 10 1)
11/2

;; So our old calculation goes like:
user> (average (/ 10 1) 1)
11/2


;; We might as well make another function which just makes our guesses better.
user> (defn improve-guess [guess] (average guess (/ 10 guess)))
#'user/improve-guess

;; Let's try that:
user> (improve-guess 1)
11/2

;; Of course, a better guess can also be improved.

user> (improve-guess (improve-guess 1))
161/44

;; And improved again.
user> (improve-guess (improve-guess (improve-guess 1)))
45281/14168

;; Now you'll notice that clojure is doing exact arithmetic, and giving us back a fraction, just
;; like a human being would do.  

;; But if we start it off with an inexact guess, a decimal rather than a ratio, 
;; say 1.0 rather than 1, then every answer we get back will be a decimal.

  
;; Because you can't make an exact answer from an inexact input.
user> (improve-guess (improve-guess (improve-guess 1.0)))
3.196005081874647

;; Inexactness is contagious.
user> (+ 1 2)
3
user> (+ 1.0 2)
3.0


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

;; Let's go on a bit of a side-track now and look at just how good Heron's idea was.  We can

;; use a magic function called iterate to make an infinite sequence of guesses.  It's best only to
;; look at the first bit of an infinite sequence.

;; Otherwise, like with the Medusa of Greek myth, the REPL turns to stone.  Let's look at the first
;; five values in our sequence.
user> (take 5 (iterate improve-guess 1))
(1 11/2 161/44 45281/14168 4057691201/1283082416)


;; We get a better idea of what's going on here if we use decimal fractions
user> (take 5 (iterate improve-guess 1.))
(1.0 5.5 3.659090909090909 3.196005081874647 3.16245562280389)

;; It only takes a few iterations before we've hit the limit of floating-point accuracy
user> (take 10 (iterate improve-guess 1.))
(1.0 5.5 3.659090909090909 3.196005081874647 3.16245562280389 3.162277665175675 3.162277660168379 3.162277660168379 3.162277660168379 3.162277660168379)


;; Let's look at the squares of the values in our sequence:
user> (map (fn[x](* x x)) (iterate improve-guess 1.))
;; ...time passes...

;; Oops. We looked at the medusa! 
;; In some environments the REPL is smart enough to print out the values one-by-one, 

;; but here, it's trying to do the whole thing before printing anything.
;; Ctrl-C will stop the calculation and restore the REPL.


;; Try again:
user> (take 10 (map (fn[x](* x x)) (iterate improve-guess 1.)))
(1.0 30.25 13.388946280991735 10.21444848336857 10.001125566203939 10.000000031668918 9.999999999999998 9.999999999999998 9.999999999999998 9.999999999999998)


;; Iterate and map are magic functions, so let's not worry too much about them yet, 
;; but the magic goes away once you know how it's done.

;; If they weren't built into the language we could write our own versions.
;; They're one liners.

;; Later I'll show you how they work.  

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

;; How do we know when to stop improving our guesses?

        
;; First off, we'd like a definition of when our guess is good enough.
;; How about, square it and if the answer's close to 10 then that goes?
;; How good a guess is 3.19?
user> (- 10 (* 3.19 3.19))
-0.17609999999999992

;; Rats, we need an absolute value function
user> (defn abs[x] (if (< x 0) -x x))

; Barf....

; Sigh. *Every* function needs the same call syntax. Unary minus is a function just like anything else.
user> (defn abs[x] (if (< x 0) (- x) x))
#'user/abs


; Let's use map to test the function.
user> (map abs (list -1 1 0))
(1 1 0)

;; So:
user> (abs (- 10 (* 3.19 3.19)))
0.17609999999999992

;; What's good enough? Let's say we're happy if the difference is less than 1/10^6, or 1e-6
user> 1e-6
1.0E-6


;; Here's our good-enough test
user> (< (abs (- 10 (* 3.19 3.19))) 1e-6)
false
;; 3.19 isn't good enough

;; Wrap the test up
user> (defn good-enough? [guess] (< (abs (- 10 (* guess guess))) 1e-6))
#'user/good-enough?


;; Let's see whether the first five values are good enough
user> (take 5 (map good-enough? (iterate improve-guess 1.)))
(false false false false false)

;; What about the first ten?
user> (take 10 (map good-enough? (iterate improve-guess 1.)))
(false false false false false true true true true true)


;; What were those answers again?
user> (take 10 (iterate improve-guess 1.))
(1.0 5.5 3.659090909090909 3.196005081874647 3.16245562280389 3.162277665175675 3.162277660168379 3.162277660168379 3.162277660168379 3.162277660168379)

;; So what if we just want a function that will give us an answer that is good enough?  We'll call
;; it good-enough-guess. 

;; We give it a guess. If that guess is good enough, then it gives us it back.

;; If it's not, then it makes it better, and tries again.
user> (defn good-enough-guess [x]
        (if (good-enough? x) x
            (good-enough-guess (improve-guess x))))
#'user/good-enough-guess

;; It doesn't really matter what our initial guess is. Anything will work.
user> (good-enough-guess 1.0)
3.162277665175675

user> (good-enough-guess 3.0)
3.162277660169842


;; How good is that guess?
user> (* (good-enough-guess 3.0)(good-enough-guess 3.0))
10.00000000000925

;; Sweet. Here endeth the method of Heron of Alexandria for finding the square root of 10.

;; Let's have a look at our program in its entirety

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

(defn average [a b] 
  (/ (+ a b) 2))

(defn improve-guess [guess]
  (average guess (/ 10 guess)))

(defn abs[x] 
  (if (< x 0) (- x) x))

(defn good-enough? [guess] 
  (< (abs (- 10 (* guess guess))) 1e-6))

(defn good-enough-guess [x]
  (if (good-enough? x) x
      (good-enough-guess (improve-guess x))))


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

RabbitMQ & Clojure: Hello World

2011-02-08T22:24:00.003Z


;; RabbitMQ is a way for different processes to talk to one another.

;; Here we're going to start a process and get it to send messages to itself via
;; rabbitMQ.

;; First we need to install rabbitmq. I'm using Ubuntu 10.10, so that's:

;; $ sudo apt-get install rabbitmq-server

;; Unfortunately Ubuntu has an old version, which the current client library is
;; unable to talk to.

;; Since the protocol has changed, it seems silly to use the old version,
;; especially since rabbit provide a debian package which seems to work fine:

;; To install the latest rabbitmq release (2.3.1 at time of writing)
;; $ wget http://www.rabbitmq.com/releases/rabbitmq-server/v2.3.1/rabbitmq-server_2.3.1-1_all.deb
;; $ sudo dpkg -i rabbitmq-server_2.3.1-1_all.deb

;; I installed this over the top of the default ubuntu installation and this
;; seems to work.  I don't know what happens if you install it
;; out-of-the-blue. Things may not get set up right, although the rabbitmq
;; website seems to indicate that it will be ok.

;; You can check that rabbitMQ is working (woc-desktop is the name of my
;; machine, don't ask..) with

;; $ sudo rabbitmqctl -n rabbit@woc-desktop status

;; Which should respond:

;; Status of node 'rabbit@woc-desktop' ...
;; [{running_applications,[{rabbit,"RabbitMQ","2.3.1"},
;;                         {mnesia,"MNESIA  CXC 138 12","4.4.12"},
;;                         {os_mon,"CPO  CXC 138 46","2.2.4"},
;;                         {sasl,"SASL  CXC 138 11","2.1.8"},
;;                         {stdlib,"ERTS  CXC 138 10","1.16.4"},
;;                         {kernel,"ERTS  CXC 138 10","2.13.4"}]},
;;  {nodes,[{disc,['rabbit@woc-desktop']}]},
;;  {running_nodes,['rabbit@woc-desktop']}]

;; The corresponding java client library is in the central maven repository so
;; you should add this snippet to your pom.xml or the corresponding incantation
;; to your project.clj if you use leiningen:

;;    <dependency>
;;      <groupId>com.rabbitmq</groupId>
;;      <artifactId>amqp-client</artifactId>
;;      <version>2.3.1</version>
;;    </dependency>


;; First we have to import the classes from the java library.
(import (com.rabbitmq.client ConnectionFactory Connection Channel QueueingConsumer))

;; And then we have to translate the equivalent java hello world program using
;; Clojure's excellent interop.

;; It feels very strange writing this sort of ceremony-oriented imperative code
;; in Clojure:

;; Make a connection factory on the local host
(def connection-factory
     (doto (ConnectionFactory.)
       (.setHost "localhost")))

;; and get it to make you a connection
(def connection (.newConnection connection-factory))

;; get that to make you a channel
(def channel (. connection createChannel))

;; And on that channel, create (or at least ensure the existence of) a queue,
;; named "hello", which is neither durable nor exclusive nor auto-deleted:
(. channel queueDeclare "hello" false false false nil)

;; Now we'll send ten messages to that queue:
(dotimes [ i 10 ]
  (. channel basicPublish "" "hello" nil (. (format "Hello World! (%d)" i) getBytes)))

;; Strictly we're sending the messages to the default exchange "", using the
;; routing key "hello"

;; Now we want to retrieve our messages from the queue.

;; We create a consumer
(def consumer (new QueueingConsumer channel))

;; And attach it to the channel
(. channel basicConsume "hello" true consumer)

;; consumer.nextDelivery() will block waiting for a message to arrive in the
;; queue, so we can just keep looping
(loop []
  (let [delivery (. consumer nextDelivery)
        str (String. (. delivery getBody))]
          (println str)
          (recur)))

;; The fun part is that we could run several copies of this program simultaneously

;; If you're using the maven-clojure-plugin, then the command to run this script
;; from the command line is:

;; $ mvn -Dclojure.script=rabbitmq.clj clojure:run

;; So you can run it on various separate terminals at once.

;; The first one you run will just send messages to itself and print them out.
;; Later copies will send messages over which earlier ones will fight. It's most entertaining!

Turning Exceptions into Return Values

2011-01-30T22:21:00.002Z


;; I often find that the default behaviour of exceptions in emacs/slime is a bit
;; annoying.

;; Up pops up a debugger window, but that moves the focus away from where
;; you're working, and the new window then needs dismissing, and has often
;; buggered up the size of other windows, etc...

;; Also when playing with experimental clojure versions, the behaviour is a bit
;; random, and often the expression just gets swallowed, with no clue as to what
;; it was:

;; Anyway I wrote this little macro, which turns an exception into a return value:
(defmacro tryc[ & e]
  `(try ~@e
       (catch Exception a# a#)))


(tryc (java.io.File ".")) ;; #<ClassCastException java.lang.ClassCastException: java.lang.Class cannot be cast to clojure.lang.IFn>

;; Obviously you should only use this at the top level!

;; I find I'm using it all the time when writing code in emacs

Finding Something in a Vector, Parsing CSV Files

2011-01-30T21:24:00.003Z


;; Indexing Into a Vector: A Sequence of Answers


;; The other day, I had a csv file to read (csv files are comma-separated
;; values, often produced by exporting spreadsheets, or a reasonably standard
;; human-readable format for grids of data)


;; There's a lovely library which takes all the hard work out of reading csv
;; files

;; The library's on clojars, so you use it by adding this to your maven pom.xml

;; <dependency>
;;   <groupId>clojure-csv</groupId>
;;   <artifactId>clojure-csv</artifactId>
;;   <version>1.2.0</version>
;; </dependency>

;; and then you can use clojure-csv.core/parse-csv to read your file in as a
;; sequence.  It seems to magically handle all the intricies of all the
;; different possible csv file arrangements.

;; Thankyou David Santiago, it works like a charm.

;; But then I hit a problem: I wanted to read the email-addresses from every
;; record.

;; I've rather simplified here, but let's imagine that the data after parsing
;; looked like this:

(def csv-data
     [["Name" "E-mail 1" "Address" "E-mail 2" "E-mail 3"]
      ["John" "john@mailinator.com" "3 Church St" "xyz@learnclojure.com" ""]
      ["Fred" "fred@mailinator.com" "7 Park Street" "fred@gmail.com" "fred@googlemail.com" ]])

(def header (first csv-data))
(def data   (rest csv-data))

;; As it happened, most of the email-addresses were stored in the second column,
;; but some records had two or even three addresses, and it looked as though the
;; generating program might be labelling its columns according to the number of
;; e-mail addresses it needed to output.

;; It seemed very likely that another data set might have an "E-mail 4" column,
;; and it seemed unwise to rely on the needed columns always being 1,3,4 and
;; possibly 5. What if the program introduced another field entirely?

;; So obviously I needed a function to look up things in the header row somehow.
;; And there didn't seem to be one, although there were the interesting
;; functions keep-indexed and map-indexed, which I hadn't noticed before.

;; And I couldn't find one. So I figured that I could either write one, or I
;; could ask on Stack Overflow.

;; And so I did both, expecting to find that I'd re-invent something in the
;; standard library, or at least in contrib, that I hadn't been able to find.

;; http://stackoverflow.com/questions/4830900/how-do-i-find-the-index-of-an-item-in-a-vector

;; So the function(s) I came up with were:

(defn indices-of [f coll]
  (keep-indexed #(if (f %2) %1 nil) coll))

(defn first-index-of [f coll]
  (first (indices-of f coll)))

(defn find-thing [value coll]
  (first-index-of #(= % value) coll))

;; And here are some examples:

(indices-of #(= "Name" %) header) ; (0)
(indices-of (partial re-matches #".*E-mail.*") header) ; (1 3 4)
(first-index-of #(= "Name" %) header) ; 0
(find-thing "Name" header) ; 0

;; But I was a bit nervous about these solutions, because I thought I must just
;; have re-invented some sort of wheel, and also because they're happy to return
;; answers for sets and maps, where the question doesn't really make much sense

;;fine
(find-thing "two" ["one" "two" "three" "two"]) ; 1
(find-thing "two" '("one" "two" "three")) ; 1

;;but these answers are a bit silly
(find-thing "two" #{"one" "two" "three"}) ; 1
(find-thing "two" {"one" "two" "two" "three"}) ; nil

;; But I went back to Stack Overflow, in order to answer my own question, and
;; found a couple of much better answers.

;; Brian Carper pointed out:

(.indexOf header "Name") ; 0

;; Which is clearly the answer for searching vectors.

;; And ponzao pointed out this lovely thing, originally due to Stuart Halloway:
(require 'clojure.contrib.seq)
(first (clojure.contrib.seq/positions #{"Name"} header)) ; 0

;; positions is like indices-of, but using a set as the predicate is really clever.

;; anyway, now I could say:
(map #( % (.indexOf header "Name")) data) ; ("John" "Fred")
(map #( % (.indexOf header "E-mail 1")) data) ; ("john@mailinator.com" "fred@mailinator.com")

;; Or even, for fans of terseness and crypticity (and forgive me Lord, for I am
;; such a fan), something like:
(map #(vector (% (.indexOf header "Name"))
              (for [i (clojure.contrib.seq/positions
                       (partial re-matches #"E-mail \d+") header)]
                (% i))) data)

;; (["John" ("john@mailinator.com" "xyz@learnclojure.com" "")]
;;  ["Fred" ("fred@mailinator.com" "fred@gmail.com" "fred@googlemail.com")])

;; But a little later, there was another answer from cgrand, who warned me
;; against using indices on general principles. And I agree with that, so I
;; asked what I should do in the particular case of csv files. And there was an
;; answer from Alex Stoddard, which I believe is the answer to the question that
;; I should have asked.

;; We can make a map from strings to indices
(def hmap (zipmap header (iterate inc 0)))

;; And use it like this:
(map #(% (hmap "Name")) data) ; ("John" "Fred")

;; or this:
(def e-mails (filter (partial re-matches #"E-mail \d+") header))

(map #(vector (% (hmap "Name")) (for [e e-mails] (% (hmap e)))) data)

;; or this:
(map #(vector
       (% (hmap "Name"))
       (for [e (filter (partial re-matches #"E-mail \d+") header)]
         (% (hmap e)))) data)

;; to get:
;; (["John" ("john@mailinator.com" "xyz@learnclojure.com" "")]
;;  ["Fred" ("fred@mailinator.com" "fred@gmail.com" "fred@googlemail.com")])

;; or even this (although now you do have to rely on the name coming before the e-mails):
(map #(for [e (filter (partial re-matches #"E-mail \d+|Name") header)]
        (% (hmap e)))
     data)

;; to get:
;; (("John" "john@mailinator.com" "xyz@learnclojure.com" "")
;;  ("Fred" "fred@mailinator.com" "fred@gmail.com" "fred@googlemail.com"))

;;If we want to abstract out a pattern, then:
(defn columns [f header data]
  (let [hmap (zipmap header (iterate inc 0))
        cols (filter f header)
    (map #(for [e cols] (% (hmap e))) data)))

;; allows:
(columns #{"Name"} header data) ; (("John") ("Fred"))
(columns (partial re-matches #"E-mail \d+") header data) ; (("john@mailinator.com" "xyz@learnclojure.com" "") ("fred@mailinator.com" "fred@gmail.com" "fred@googlemail.com"))

£500 if you can find me a job (version 1.0)

2011-01-27T20:10:00.000Z

That worked a treat, £500 awarded to lucky winner who found me a fun short contract.

Now I'd like another one, so I repeat the offer:

If, within the next six months, I take a job which lasts longer than one month, and that is not obtained through an agency, then on the day the first cheque from that job cashes, I'll give £500 to the person who provided the crucial introduction.

If there are a number of people involved somehow, then I'll apportion it fairly between them. And if the timing conditions above are not quite met, or if someone points me at a short contract which the £500 penalty makes not worth taking, then I'll do something fair and proportional anyway. (The thing Simon pointed me at only lasted three weeks and I paid him in full anyway, because it was neat.)

And this offer applies even to personal friends, and to old contacts who I have not got round to calling yet, and to people who are themselves offering work, because why wouldn't it?

And obviously if I find one through my own efforts then I'll keep the money. But my word is generally thought to be good, and I have made a public promise on my own blog to this effect, so if I cheat you you can blacken my name and ruin my reputation for honesty, which is worth much more to me than £500.

Anyhow, my CV is at http://www.aspden.com, and any advice on how it could be improved will be gratefully received.

I'll also repeat the original advert: (job-hunt.contrib 1.0)

Anyone in Cambridge need a programmer? Obviously Clojure is a speciality, and my current obsession, but I'm also pretty good with C (especially the embedded variety), microcontrollers, and Python, and I have a particular facility with mathematical concepts and algorithms of all kinds. My obsessions can be pretty quickly changed when I need them to be.

I have a (deserved) reputation for being able to produce heavily optimised but nevertheless bug-free and readable code, but I also know how to hack together sloppy, bug-ridden prototypes, and I know which style is appropriate when, and how to slide along the continuum between them.

I've worked in telecoms, commercial research, banking, university research, a chip design company, server virtualization, a couple of startups, and occasionally completely alone.

I've worked on many sizes of machine. I've written programs for tiny 8-bit microcontrollers, and once upon a time every IBM machine in one building in Imperial College was running my partial differential equation solvers in parallel in the background.

I'm smart and I get things done. I'm confident enough in my own abilities that if I can't do something I admit it and find someone who can.

I also have various ancient and rusty skills with things like Java, C++, R, OCaml, Common LISP, Scheme, FORTRAN and Pascal which can be brushed up if necessary. Like all lispers, I occasionally write toy interpreters for made-up languages for fun.

If you're a local company using Java, who might be interested in giving Clojure a try (motivation here, in Paul Graham's classic Beating the Averages), I'd love to try to show you what all the fuss is about.

CV here if you're interested: http://www.aspden.com

I've never used a CV before, having always found work through word of mouth. So I expect that it can be improved. If anyone's got any suggestions as to how it could be better written, do please leave comments or e-mail cv@aspden.com.

K-means : An Algorithm for Clustering Data

2011-01-25T23:37:00.000Z


;; K-means

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

;; K-means is a Clustering Algorithm.

;; We use it when we have some data, and we want to split the data into separate categories.

;; For instance, an early biologist, let's call him Adam, might measure all
;; sorts of things about the living objects he encounters in the world.

;; (black, feathers, flies)
;; (green, tasty)
;; (green, slithers, poisonous)
;; ... and so on

;; After he collects enough data, he might go looking for structure in it.

;; Uninformed by theory, he might nevertheless notice that many things that do
;; not move are green, and that many things that are not green move.

;; He might name these two obvious groups the Animals and the Plants.

;; Further analysis of the data might split the Plants into Trees and Flowers,
;; and the Animals into Mammals, Birds, and Fish.

;; Theoretically, this process could continue further, extracting 'natural
;; categories' from the observed structure of the data, without any theory about
;; how the various properties come about

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

;; Let's consider a very simple clustering situation. We have some numbers,
;; and we'd like to see if they form groups.

;; Suppose we want to cluster the points 2 3 5 6 10 11 100 101 102
(def data '(2 3 5 6 10 11 100 101 102))
;; You may be able to see some natural groupings in this data.

;; It's easy enough to say how far one number is from another
(defn distance[a b] = (if (< a b) (- b a) (- a b)))

;; To do K-means, we need to start with some guesses about where the clusters are.
;; They don't have to be terribly good guesses.
(def guessed-means '(0 10))

;; Given a particular point, we want to know which of our means is closest
(defn closest [point means distance]
  (first (sort-by #(distance % point) means)))

;; In our little dataset, 2 is closest to the guess of 0, and 100 is closest to the guess of 10
(closest 2   guessed-means distance) ; 0
(closest 100 guessed-means distance) ; 10

;; So we can talk about the group of all points for which 0 is the best guess
;; and the group of all points for which 10 is the best guess.
(defn point-groups [means data distance]
  (group-by #(closest % means distance) data))

(point-groups guessed-means data distance) ; {0 [2 3 5], 10 [6 10 11 100 101 102]}

;; We can take an average of a group of points
(defn average [& list] (/ (reduce + list) (count list)))

(average 6 10 11 100 101 102) ; 55

;; So we can take the average of each group, and use it as a new guess for where
;; the clusters are. If a mean doesn't have a group, then we'll leave it where
;; it is.
(defn new-means [average point-groups old-means]
  (for [o old-means]
    (if (contains? point-groups o)
      (apply average (get point-groups o)) o)))

(new-means average (point-groups guessed-means data distance) guessed-means) ; (10/3 55)

;; So if we know we've got a particular set of points, and a particular idea of
;; distance, and a way of averaging things, that gives us a way of making a new
;; list of guesses from our original list of guesses
(defn iterate-means [data distance average]
  (fn [means] (new-means average (point-groups means data distance) means)))

((iterate-means data distance average) '(0 10)) ; (10/3 55)

;; and of course we can use that as a new guess, and improve it again.
((iterate-means data distance average) '(10/3 55)) ; (37/6 101)

;; We can do this repeatedly until it settles down.
(iterate (iterate-means data distance average) '(0 10)) ; ((0 10) (10/3 55) (37/6 101) (37/6 101) .....)

;; K-means with two means seems to be trying to tell us that we've got a group
;; centered around 6 and another centred around 101

;; These groups are:
(defn groups [data distance means]
  (vals (point-groups means data distance)))

(groups data distance '(37/6 101)) ; ([2 3 5 6 10 11] [100 101 102])

;; Ideally we'd like to iterate until the groups stop changing.
;; I described a function for doing this in a previous post:
(defn take-while-unstable 
  ([sq] (lazy-seq (if-let [sq (seq sq)]
                    (cons (first sq) (take-while-unstable (rest sq) (first sq))))))
  ([sq last] (lazy-seq (if-let [sq (seq sq)]
                         (if (= (first sq) last) '() (take-while-unstable sq))))))

(take-while-unstable '(1 2 3 4 5 6 7 7 7 7)) ; (1 2 3 4 5 6 7)

(take-while-unstable (map #(groups data distance %) (iterate (iterate-means data distance average) '(0 10))))
                                        ; (([2 3 5] [6 10 11 100 101 102])
                                        ;  ([2 3 5 6 10 11] [100 101 102]))

;; Shows that our first guesses group 2,3 and 5 (nearer to 0 than 10) vs all the rest.
;; K-means modifies that instantly to separate out the large group of three.


;; We can make a function, which takes our data, notion of distance, and notion of average,
;; and gives us back a function which, for a given set of initial guesses at the means,
;; shows us how the group memberships change.
(defn k-groups [data distance average]
  (fn [guesses]
    (take-while-unstable
     (map #(groups data distance %)
          (iterate (iterate-means data distance average) guesses)))))

(def grouper (k-groups data distance average))

(grouper '(0 10))
                                        ; (([2 3 5] [6 10 11 100 101 102])
                                        ;  
                                        ;  ([2 3 5 6 10 11] [100 101 102]))


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

;; Nothing we said above limits us to only having two guesses
(grouper '(1 2 3))
                                        ; (([2] [3 5 6 10 11 100 101 102])
                                        ;  ([2 3 5 6 10 11] [100 101 102])
                                        ;  ([2 3] [5 6 10 11] [100 101 102])
                                        ;  ([2 3 5] [6 10 11] [100 101 102])
                                        ;
                                        ;  ([2 3 5 6] [10 11] [100 101 102]))


;; The more means we start with, the more detailed our clustering.
(grouper '(0 1 2 3 4))
                                        ; (([2] [3] [5 6 10 11 100 101 102])
                                        ;  ([2] [3 5 6 10 11] [100 101 102])
                                        ;  ([2 3] [5 6 10 11] [100 101 102])
                                        ;  ([2 3 5] [6 10 11] [100 101 102])
                                        ;  ([2] [3 5 6] [10 11] [100 101 102])
                                        ;
                                        ;  ([2 3] [5 6] [10 11] [100 101 102]))

;; We have to be careful not to start off with too many means, or we just get our data back:
(grouper (range 200)) ; (([2] [3] [100] [5] [101] [6] [102] [10] [11]))

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

;; Generalizing to Other Spaces

;; In fact, nothing we said above depends on our inputs being numbers

;; We can use any data where we can define a distance, and a method of averaging:

;; One of the easiest things to do this for would be vectors:

(defn vec-distance [a b] (reduce + (map #(* % %) (map - a b))))
(defn vec-average  [& l] (map #(/ % (count l)) (apply map + l)))

(vec-distance [1 2 3][5 6 7]) ; 48
(vec-average  [1 2 3][5 6 7]) ; (3 4 5)

;; Here's a little set of vectors
(def vector-data '( [1 2 3] [3 2 1] [100 200 300] [300 200 100] [50 50 50]))

;; And choosing three guesses in a fairly simple-minded manner, we can see how the algorithm
;; divides them into three different groups.
((k-groups vector-data vec-distance vec-average) '([1 1 1] [2 2 2] [3 3 3]))
                                        ; (([[1 2 3] [3 2 1]] [[100 200 300] [300 200 100] [50 50 50]])

                                        ;  ([[1 2 3] [3 2 1] [50 50 50]]
                                        ;   [[100 200 300] [300 200 100]])

                                        ;  ([[1 2 3] [3 2 1]]
                                        ;   [[100 200 300] [300 200 100]]
                                        ;   [[50 50 50]]))

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

;; Pedantic Footnote

;; Note that the algorithm as described above isn't quite the classic K-means.
;; I don't think the difference is terribly important, and I thought it would complicate the explanation to deal with it.

;; In the usual K-means, if you have two identical means, then you're only supposed to update one of them.

;; Here our two identical guesses are both getting updated
(new-means average (point-groups '(0 0) '(0 1 2 3 4) distance) '(0 0)) ; (2 2)

;; Our update function:
(defn new-means [average point-groups old-means]
  (for [o old-means]
    (if (contains? point-groups o)
      (apply average (get point-groups o)) o)))

;; Needs to be changed so that if there are two identical means only one of them will be changed:

(defn update-seq [sq f]
  (let [freqs (frequencies sq)]
    (apply concat
     (for [[k v] freqs]
       (if (= v 1) (list (f k))
           (cons (f k) (repeat (dec v) k)))))))


(defn new-means [average point-groups old-means]
  (update-seq old-means (fn[o]
                          (if (contains? point-groups o)
                            (apply average (get point-groups o)) o))))

;; Now only one will get updated at once
(new-means average (point-groups '(0 0) '(0 1 2 3 4) distance) '(0 0)) ; (2 0)


;; Now we don't lose groups when the means get aliased.
((k-groups '(0 1 2 3 4) distance average) '(0 1)) ; (([0] [1 2 3 4]) ([0 1] [2 3 4]))
((k-groups '(0 1 2 3 4) distance average) '(0 0)) ; (([0 1 2 3 4]) ([0] [1 2 3 4]) ([0 1] [2 3 4]))

((k-groups vector-data vec-distance vec-average) '([1 1 1] [1 1 1] [1 1 1])) ;
                                        ; (([[1 2 3] [3 2 1] [100 200 300] [300 200 100] [50 50 50]])
                                        ;  ([[1 2 3] [3 2 1]] [[100 200 300] [300 200 100] [50 50 50]])
                                        ;  ([[1 2 3] [3 2 1] [50 50 50]] [[100 200 300] [300 200 100]])
                                        ;  ([[1 2 3] [3 2 1]] [[100 200 300] [300 200 100]] [[50 50 50]]))

;; Although it's still possible that a mean never acquires any points, so we can still get out fewer groups than means.
((k-groups '(0 1 2 3 4) distance average) '(0 5 10)) ; (([0 1 2] [3 4]))

Cleaning Old Definitions from the REPL : shred-user

2011-01-25T15:22:00.002Z


;; I often find myself restarting the clojure repl in order to get rid of old
;; definitions which introduce subtle bugs.

;; Consider:

(defn factorial [n] (if (< n 2) 1 (* n (factorial (dec n))))) ; #'user/factorial

(factorial 10) ; 3628800

;; Suppose instead I wanted (factorial2 n) to be (* 1 2 2 2 3 2 4 .... 2 n)

(defn factorial2 [n] (if (< n 2) 1 (* 2 n (factorial (dec n))))) ; #'user/factorial2

;; There is a subtle bug introduced by my failing to rename the recursive call

(map factorial2 '(1 2 3 4)) ; (1 4 24 192)

;; But (factorial2 4) should be (* 1 2 2 2 3 2 4) = 192

;; I suspect that there is an old definition hanging around confusing matters.

;; To find it, I can run the program from the command line, or restart the repl and recompile, or:

(doseq [s (map first (ns-interns 'user))](ns-unmap 'user s)) ; nil

(defn factorial2 [n] (if (< n 2) 1 (* 2 n (factorial (dec n))))) ; Unable to resolve symbol: factorial in this context

;; Aha! 

(defn factorial2 [n] (if (< n 2) 1 (* 2 n (factorial2 (dec n))))) ; #'user/factorial2

(factorial2 4) ; 192


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

;; For convenience I'd like to make that a function:

(defn shred-user []
  (doseq [s (map first (ns-interns 'user))] (ns-unmap 'user s)))

(shred-user) ; nil

(shred-user) ; Unable to resolve symbol: shred-user in this context ; bugger!

(defn shred-user []
  (doseq [s (filter (complement #{'shred-user}) (map first (ns-interns 'user)))] (ns-unmap 'user s)))

(shred-user) ; nil

(shred-user) ; nil

take-while-unstable

2011-01-25T13:43:00.000Z


;; I often find myself wanting to iterate something until it converges.

(defn collatz [n] (if (even? n) (recur (/ n 2)) (+ 1 (* 3 n))))

(iterate collatz 100) ; (100 76 58 88 34 52 40 16 4 4 4 4 4 4 ...)

;; And so I often find myself writing:

(defn take-while-unstable
  ([sq] (if (empty? sq) '()
            (cons (first sq) (#'take-while-unstable (rest sq) (first sq)))))
  ([sq last] (cond (empty? sq) '()
                   (= (first sq) last) '()
                   :else (cons (first sq) (#'take-while-unstable (rest sq) (first sq))))))

;; Which allows me to stop the iteration once it's settled down.
(take-while-unstable (iterate collatz 100)) ; (100 76 58 88 34 52 40 16 4)

;; You can see how it works here:
(require 'clojure.contrib.trace)

(binding  [*print-length* 20] ; taking care not to look at the medusa
  (clojure.contrib.trace/dotrace (take-while-unstable) (take-while-unstable (iterate collatz 100))))

;; TRACE t15884: (take-while-unstable (100 76 58 88 34 52 40 16 4 4 4 4 4 4 4 4 4 4 4 4 ...))
;; TRACE t15885: |    (take-while-unstable (76 58 88 34 52 40 16 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 100)
;; TRACE t15886: |    |    (take-while-unstable (58 88 34 52 40 16 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 76)
;; TRACE t15887: |    |    |    (take-while-unstable (88 34 52 40 16 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 58)
;; TRACE t15888: |    |    |    |    (take-while-unstable (34 52 40 16 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 88)
;; TRACE t15889: |    |    |    |    |    (take-while-unstable (52 40 16 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 34)
;; TRACE t15890: |    |    |    |    |    |    (take-while-unstable (40 16 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 52)
;; TRACE t15891: |    |    |    |    |    |    |    (take-while-unstable (16 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 40)
;; TRACE t15892: |    |    |    |    |    |    |    |    (take-while-unstable (4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 16)
;; TRACE t15893: |    |    |    |    |    |    |    |    |    (take-while-unstable (4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...) 4)
;; TRACE t15893: |    |    |    |    |    |    |    |    |    => ()
;; TRACE t15892: |    |    |    |    |    |    |    |    => (4)
;; TRACE t15891: |    |    |    |    |    |    |    => (16 4)
;; TRACE t15890: |    |    |    |    |    |    => (40 16 4)
;; TRACE t15889: |    |    |    |    |    => (52 40 16 4)
;; TRACE t15888: |    |    |    |    => (34 52 40 16 4)
;; TRACE t15887: |    |    |    => (88 34 52 40 16 4)
;; TRACE t15886: |    |    => (58 88 34 52 40 16 4)
;; TRACE t15885: |    => (76 58 88 34 52 40 16 4)
;; TRACE t15884: => (100 76 58 88 34 52 40 16 4)
;; (100 76 58 88 34 52 40 16 4)

;; Anyway, this function comes in so handy that I thought I'd lazyize it and give it some tests

(use 'clojure.test)

(with-test
  
  (defn take-while-unstable
    ([sq] (lazy-seq (if-let [sq (seq sq)]
                      (cons (first sq) (take-while-unstable (rest sq) (first sq))))))
    ([sq last] (lazy-seq (if-let [sq (seq sq)]
                           (if (= (first sq) last) '() (take-while-unstable sq))))))

  (is (= (take-while-unstable '()) '()))
  (is (= (take-while-unstable '(1)) '(1)))
  (is (= (take-while-unstable '(1 1)) '(1)))
  (is (= (take-while-unstable '(1 2)) '(1 2)))
  (is (= (take-while-unstable '(1 1 2)) '(1)))
  (is (= (take-while-unstable '(1 1 1)) '(1)))
  (is (= (take-while-unstable '(1 1 2 1)) '(1)))
  (is (= (take-while-unstable '(1 2 3 4 5 6 7 7 7 7)) '(1 2 3 4 5 6 7)))
  (is (= (take 10000 (take-while-unstable (range))) (take 10000 (range))))
  (is (= (take-while-unstable (concat (take 10000 (range)) '(10000) (drop 10000 (range)))) (range 10001)))
  (is (= (take-while-unstable '(nil)) '(nil)))
  (is (= (take-while-unstable '(nil nil)) '(nil)))
  (is (= (take-while-unstable '[nil nil]) '(nil)))
  (is (= (take-while-unstable [ :a :b :c :d :d :e]) '(:a :b :c :d))))

(run-tests) ; {:type :summary, :test 1, :pass 14, :fail 0, :error 0}

;; Hopefully someone will find it useful. Can anyone break it?

A Very Gentle Introduction to Information Theory: Guessing the Entropy

2011-01-17T23:09:00.003Z


;; A Very Gentle Introduction to Information Theory : Part IV

;; Entropy and Huffman Coding

;; Once again, I'll keep some code from the earlier parts, without explanation

;; Probability distributions

(defn combine-keywords [& a] (keyword (apply str (mapcat #(drop 1 (str %)) a))))

(defn combine-distributions
  ([P] P)
  ([P1 P2]
     (into {}
           (for [[s1 p1] P1
                 [s2 p2] P2]
             [(combine-keywords s1 s2) (* p1 p2)])))
  ([P1 P2 & Plist] (reduce combine-distributions (cons P1 (cons P2 Plist)))))

;; Coding and Decoding

(defn decoder
  ([code-tree stream] (decoder code-tree code-tree stream))
  ([current-code-tree code-tree stream]
     (lazy-seq
        (if (keyword? current-code-tree)
          (cons current-code-tree (decoder code-tree code-tree stream))
          (if-let [stream (seq stream)]
            (if (= (first stream) 1)
              (decoder (first current-code-tree)  code-tree (rest stream))
              (decoder (second current-code-tree) code-tree (rest stream))))))))

(defn encoder [code stream] (mapcat code stream))

(defn make-encoder [code]  (fn [s] (encoder code s)))
(defn make-decoder [code-tree] (fn[s] (decoder code-tree s)))

;; Huffman encoding

(defn symbols [prefix code-tree]
  (if (keyword? code-tree) (list prefix code-tree)
      (concat (symbols (cons 1 prefix) (first code-tree))
              (symbols (cons 0 prefix) (second code-tree)))))

(defn make-code [code-tree]
  (into {} (map (fn[[c s]][s (reverse c)]) (partition 2 (symbols '() code-tree)))))

;; The original make-code-tree was very slow, because it re-sorted the list every
;; iteration. We can speed it up considerably by using a priority queue instead
;; of sorting the list every iteration. Clojure doesn't have one built in, so
;; here's a poor man's version built out of a sorted map of lists.

(defn madd [m [k p]] (assoc m p (cons k (m p))))

(defn mpop [m]
  (let [[k vlist] (first m)]
    [k (first vlist)
    (if (empty? (rest vlist))
      (dissoc m k)
      (assoc m k (rest vlist)))]))

(defn mcombine [m]
  (let [[pa a pop1] (mpop m)
        [pb b pop2] (mpop pop1)]
    (madd pop2 [[a b] (+ pa pb)])))

(defn make-code-tree [P]
  (second (mpop
          (nth (iterate mcombine (reduce madd (sorted-map) P)) (dec (count P))))))

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

(def fair-coin {:H 1 :T 1})
(def unfair-coin {:H 3 :T 1})
(def unfair-triples (combine-distributions unfair-coin unfair-coin unfair-coin))

(def triple-code-tree (make-code-tree unfair-triples))
(def triple-code (make-code triple-code-tree))


;; We don't need to estimate the cost of a code by generating a long random stream and transmitting it.

;; Given a probability distribution and a code, we can just calculate the expected cost:

;; We make a distribution over the transmitted symbols
(defn code-distribution [P code]  (for [s (keys P)] [(P s) (code s)]))

;; e.g.
(code-distribution unfair-triples triple-code)
;;([27 (0)] [9 (1 1 1)] [9 (1 1 0)] [3 (1 0 0 1 0)] [9 (1 0 1)] [3 (1 0 0 0 1)] [3 (1 0 0 0 0)] [1 (1 0 0 1 1)])

;; And from that, it's easy to calculate the expected length of a sequence
(defn expected-code-length [P code]
  (let [cd (code-distribution P code)]
    (/ (reduce + (for [[k v] cd] (* k (count v))))
       (reduce + (map first cd)))))

;; So the expected cost per symbol is:
(expected-code-length unfair-triples triple-code) ;79/32
;; or per coin-toss:
(float ( / 79/32 3)) ; 0.8229167

;; So we can get the noise out of our table:

(defn cost-for-n-code [ P n ]
     (let [Pn (apply combine-distributions (repeat n P))
           code (make-code (make-code-tree Pn))]
       (float (/ (expected-code-length Pn code) n))))

(cost-for-n-code unfair-coin 1) ; 1.0
(cost-for-n-code unfair-coin 2) ; 0.84375
(cost-for-n-code unfair-coin 3) ; 0.8229167
(cost-for-n-code unfair-coin 4) ; 0.8183594
(cost-for-n-code unfair-coin 5) ; 0.81777346
(cost-for-n-code unfair-coin 6) ; 0.8186849
(cost-for-n-code unfair-coin 7) ; 0.81685966
(cost-for-n-code unfair-coin 8) ; 0.81575775
(cost-for-n-code unfair-coin 9) ; 0.81493336
(cost-for-n-code unfair-coin 10) ; 0.8141917
(cost-for-n-code unfair-coin 11) ; 0.8137328
(cost-for-n-code unfair-coin 12) ; 0.81351095

;; It looks like something is converging, although the convergence isn't monotonic. 
;; I'm now revising my estimate of the cost of sending the results of a 1:3 process to be about 0.813

;; But we don't know whether that's the limit of the coding process, or a property of the 1:3 distribution,
;; or whether it's in some way specific to transmitting over a binary channel.

;; Let's look at some other distributions.

;; For the fair coin distribution, huffman coding triples doesn't help at all.
(cost-for-n-code fair-coin 1) ; 1.0
(cost-for-n-code fair-coin 2) ; 1.0
(cost-for-n-code fair-coin 3) ; 1.0
(cost-for-n-code fair-coin 4) ; 1.0

;; But for an even choice between three things, it does:
(def triad {:A 1 :B 1 :C 1})

(cost-for-n-code triad 1) ; 1.6666666
(cost-for-n-code triad 2) ; 1.6111112
(cost-for-n-code triad 3) ; 1.6049383
(cost-for-n-code triad 4) ; 1.6049383
(cost-for-n-code triad 5) ; 1.5893004
(cost-for-n-code triad 6) ; 1.5992227
(cost-for-n-code triad 7) ; 1.5895878
(cost-for-n-code triad 8) ; 1.5939262



;; For a choice between four things, it makes no difference
(def quad {:A 1 :B 1 :C 1 :D 1})

(cost-for-n-code quad 1) ; 2.0
(cost-for-n-code quad 2) ; 2.0
(cost-for-n-code quad 3) ; 2.0
(cost-for-n-code quad 4) ; 2.0
(cost-for-n-code quad 5) ; 2.0

;; For five it's a good thing to do
(def quint {:A 1 :B 1 :C 1 :D 1 :E 1})

(cost-for-n-code quint 1) ; 2.4
(cost-for-n-code quint 2) ; 2.36
(cost-for-n-code quint 3) ; 2.3253334
(cost-for-n-code quint 4) ; 2.3404
(cost-for-n-code quint 5) ; 2.337856
(cost-for-n-code quint 6) ; 2.3252373


;; And again, for the next power of two, no difference.
(def octet {:A 1 :B 1 :C 1 :D 1 :E 1 :F 1 :G 1 :H 1})

(cost-for-n-code octet 1) ; 3.0
(cost-for-n-code octet 2) ; 3.0
(cost-for-n-code octet 3) ; 3.0

;; I think that we might have guessed that it would take three bits to decide between
;; eight equally likely things, and two bits for four things, but what about the other numbers?

;; If 8 = 2*2*2 -> 3, and 4 = 2*2 -> 2, and 2 -> 1, what's the easiest pattern we can fit to that?

(defn bits [n] (/ (Math/log n)(Math/log 2))) ;; Also known as log2

(map bits (range 2 10)) ; (1.0 1.5849625007211563 2.0 2.321928094887362 2.584962500721156 2.807354922057604 3.0 3.1699250014423126)

;; So let's make a prediction
(def sextet {:A 1 :B 1 :C 1 :D 1 :E 1 :F 1})
(bits 6) ; 2.584962500721156

(cost-for-n-code sextet 1) ; 2.6666667
(cost-for-n-code sextet 2) ; 2.6111112
(cost-for-n-code sextet 3) ; 2.6049383
(cost-for-n-code sextet 4) ; 2.6049383
(cost-for-n-code sextet 5) ; 2.5893004
(cost-for-n-code sextet 6) ; 2.5992227


;; It looks as though the cost of coding an even distribution using huffman
;; encoding of runs is pretty close to being the logarithm (to base 2) of the number of symbols.

A Very Gentle Introduction to Information Theory: Huffman Coding

2011-01-15T18:46:00.000Z


;; A Very Gentle Introduction to Information Theory : Part III

;; Entropy and Huffman Coding

;; Once again, I'll keep some code from the first two parts, without explanation

(defn random-stream [P]
  (let [pseq (vec (mapcat (fn[[k v]](repeat v k )) P))]
    (for [i (range)] (rand-nth pseq))))

(defn cost [encoder decoder message]
  (let [coded (encoder message)]
    (if (= (decoder coded) message) (count coded) :fail)))

(defn decoder
  ([code-tree stream] (decoder code-tree code-tree stream))
  ([current-code-tree code-tree stream]
     (lazy-seq
        (if (keyword? current-code-tree)
          (cons current-code-tree (decoder code-tree code-tree stream))
          (if-let [stream (seq stream)]
            (if (= (first stream) 1)
              (decoder (first current-code-tree)  code-tree (rest stream))
              (decoder (second current-code-tree) code-tree (rest stream))))))))

(defn encoder [code stream] (mapcat code stream))

(defn make-encoder [code]  (fn [s] (encoder code s)))
(defn make-decoder [code-tree] (fn[s] (decoder code-tree s)))

(defn combine-keywords [& a] (keyword (apply str (mapcat #(drop 1 (str %)) a))))
(defn split-keyword [a] (map #(keyword (str %)) (drop 1 (str a))))

(defn make-combination-encoder [code n]
  (fn [s] (encoder code (map #(apply combine-keywords %) (partition n s)))))

(defn make-combination-decoder [code-tree]
  (fn [s] (mapcat split-keyword (decoder code-tree s))))

;; So far we've looked at three probability distributions:
(def fair-coin {:H 1 :T 1})
(def unfair-coin {:H 3 :T 1})
(def unfair-pairs {:HH 9, :HT 3, :TH 3, :TT 1})

;; And two codes:
(def fair-code-tree [:H :T])
(def fair-code {:H '(1) :T '(0)})
     
(def unfair-code-tree [ :HH [ :HT [ :TH :TT]]])
(def unfair-code {:HH '(1) :HT '(0 1) :TH '(0 0 1) :TT '(0 0 0)})

;; We should add a fourth probability distribution to represent pairs of fair coin toss results
(def fair-pairs {:HH 1 :HT 1 :TH 1 :TT 1})

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

;; We found that

(defn estimate-cost [P encoder decoder]
  (let [n 100000
        c (cost encoder decoder (take n (random-stream P)))]
    (if (number? c) (float (/ c n)) c)))

;; Using the best code we can think of for the fair coin resulted in a transmission cost of 1 (£/symbol)
(estimate-cost fair-coin   (make-encoder fair-code) (make-decoder fair-code-tree)) ; 1.0
;; And that that was also the cost for the unfair coin with this code:
(estimate-cost unfair-coin (make-encoder fair-code) (make-decoder fair-code-tree)) ; 1.0

;; But that we could come up with a code for pairs of coin tosses
;; which did substantially better for pairs of unfair coin tosses
(estimate-cost unfair-pairs (make-encoder unfair-code) (make-decoder unfair-code-tree)) ; 1.68338
;; but substantially worse for pairs of fair coin tosses
(estimate-cost fair-pairs   (make-encoder unfair-code) (make-decoder unfair-code-tree)) ; 2.24722
;; remember that that's the transmission cost per symbol, and that each symbol represents two coin tosses

;; In case you think there's any sleight of hand going on there, here's how we'd use the pairs code to transmit
;; the original unpaired streams
(estimate-cost unfair-coin  (make-combination-encoder unfair-code 2) (make-combination-decoder unfair-code-tree)) ; 0.84561
(estimate-cost fair-coin    (make-combination-encoder unfair-code 2) (make-combination-decoder unfair-code-tree)) ; 1.12407
;; Notice that the costs here are per-toss, showing that the unfair code is actually an improvement 

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

;; Now it seems that, if we can send the results of unfair-coin more efficiently
;; by considering {:HH 9, :HT 3, :TH 3, :TT 1}, the distribution of pairs of
;; tosses, then we should have a look at the distribution of triples, and work out a code for that:

(defn combine-distributions
  ([P] P)
  ([P1 P2]
     (into {}
           (for [[s1 p1] P1
                 [s2 p2] P2]
             [(combine-keywords s1 s2) (* p1 p2)])))
  ([P1 P2 & Plist] (reduce combine-distributions (cons P1 (cons P2 Plist)))))

(def unfair-triples (combine-distributions unfair-coin unfair-coin unfair-coin))

;; unfair-triples is {:HHH 27, :HHT 9, :HTH 9, :HTT 3, :THH 9, :THT 3, :TTH 3, :TTT 1}

;; Now how should we work out a code for this distribution?

;; Huffman tells us that we should combine the two lowest probability events
;; so that
{:HHH 27, :HHT 9, :HTH 9, :HTT 3, :THH 9, :THT 3, :TTH 3, :TTT 1}
;; goes to
{:HHH 27, :HHT 9, :HTH 9, :HTT 3, :THH 9, :THT 3, {:TTH 3, :TTT 1} 4}
;; and then do it again, so that
{:HHH 27, :HHT 9, :HTH 9, :HTT 3, :THH 9, :THT 3, {:TTH 3, :TTT 1} 4}
;; goes to
{:HHH 27, :HHT 9, :HTH 9, :HTT 3, :THH 9, {:THT 3, {:TTH 3, :TTT 1} 4} 7}
;; and so on .....

(defn huffman-combine [P]
  (let [plist (sort-by second P)
        newelement (into {} (take 2 plist))]
    (into {} (cons [newelement (reduce + (vals newelement))] (drop 2 plist)))))

(nth (iterate huffman-combine unfair-triples) (dec (dec (count unfair-triples))))
;; {{{:HHT 9, :HTH 9} 18, {:THH 9, {{:TTT 1, :HTT 3} 4, {:THT 3, :TTH 3} 6} 10} 19} 37, :HHH 27}

;; At the end, we get a sort of nested binary probability distribution
;; You could think of this as being a way to generate the triples by tossing strangely biased coins!

;; From that, we can generate our code tree directly by just throwing away the numbers
(require 'clojure.walk)
(defn make-code-tree [P]
  (clojure.walk/postwalk #(if (map? %) (into[] (map first  %)) %)
                         (nth (iterate huffman-combine P) (dec (dec (count P))))))

(def triple-code-tree (make-code-tree unfair-triples))
;;[[[:HHT :HTH] [:THH [[:TTT :HTT] [:THT :TTH]]]] :HHH]

;; If we have the decoder, then we can use it to generate the coder!

(defn symbols [prefix code-tree]
  (if (keyword? code-tree) (list prefix code-tree)
      (concat (symbols (cons 1 prefix) (first code-tree))
              (symbols (cons 0 prefix) (second code-tree)))))

(defn make-code [code-tree]
  (into {} (map (fn[[c s]][s (reverse c)]) (partition 2 (symbols '() code-tree)))))

(def triple-code (make-code triple-code-tree))
;; {:HHT (0 0 0), :HTH (0 0 1), :THH (0 1 0), :TTT (0 1 1 0 0), :HTT (0 1 1 0 1), :THT (0 1 1 1 0), :TTH (0 1 1 1 1), :HHH (1)}

;; Let's see how this does
(estimate-cost unfair-triples (make-encoder triple-code) (make-decoder triple-code-tree)) ; 2.4615

;; £2.46 per symbol, or 0.82p per toss

;; So while going from single tosses to pairs allowed us to go from 1->0.85, going from pairs to triples only allowed us to get from 0.85->0.82.

;; Is there an end to this process?

(defn bit-rate [P n]
  (let [Pn (apply combine-distributions (repeat n P))
        tree (make-code-tree Pn)]
    (/ (estimate-cost Pn (make-encoder (make-code tree)) (make-decoder tree)) n)))

(bit-rate unfair-coin 1) ; 1.0
(bit-rate unfair-coin 2) ; 0.844435
(bit-rate unfair-coin 3) ; 0.82466
(bit-rate unfair-coin 4) ; 0.8196275
(bit-rate unfair-coin 5) ; 0.818912
(bit-rate unfair-coin 6) ; 0.81896996
(bit-rate unfair-coin 7) ; 0.8166514
(bit-rate unfair-coin 8) ; 0.815995
(bit-rate unfair-coin 9) ; 0.81352335
(bit-rate unfair-coin 10) ; 0.81462896
(bit-rate unfair-coin 11) ; 0.8137691

;; To me, this is at least suggestive that there might be something fundamental
;; about a cost of 82p to transmit the result of a 3:1 random result over a
;; binary channel.

A Very Gentle Introduction to Information Theory : Generalizing the Encoder and Decoder

2011-01-14T00:57:00.000Z


;; A Very Gentle Introduction to Information Theory : Part II

;; Entropy and Huffman Coding

;; Here's the essential code from part I, which I'm not going to explain again:

(defn random-stream [P]
  (let [pseq (vec (mapcat (fn[[k v]](repeat v k )) P))]
    (for [i (range)] (rand-nth pseq))))

(defn cost [encoder decoder message]
  (let [coded (encoder message)]
    (if (= (decoder coded) message) (count coded) :fail)))

(def unfair-pairs {:HH 9, :HT 3, :TH 3, :TT 1})

;; We're trying to transmit the output of the random process represented by:

(def stream (random-stream unfair-pairs))

(take 20 stream) ;(:HH :HH :HH :HH :HH :HH :HH :HH :HT :HH :HH :TH :HH :HH :HH :TT :HH :HH :HT :HT)

;; And we're using the code HH -> 1, HT ->01 TH->001, TT-> 000

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

;; Our code seems to have something of a tree structure

;; 1-> :HH
;; 0-> 1 -> :HT
;;     0 -> 1 -> :TH
;;          0 -> :TT


;; Let's see if we can find some way of expressing that, so that we don't have to hand-code a decoder
;; for every different code.

(def code-tree [ :HH [ :HT [ :TH :TT]]])

(defn decoder
  ([code-tree stream] (decoder code-tree code-tree stream))
  ([current-code-tree code-tree stream]
     (lazy-seq
        (if (keyword? current-code-tree)
          (cons current-code-tree (decoder code-tree code-tree stream))
          (if-let [stream (seq stream)]
            (if (= (first stream) 1)
              (decoder (first current-code-tree)  code-tree (rest stream))
              (decoder (second current-code-tree) code-tree (rest stream))))))))

(decoder code-tree '(0 1 1 0 1 0 1 1 0 0 0)) ;(:HT :HH :HT :HT :HH :TT)
  
;; A general encoder, by comparison, is fairly straightforward:

(def code {:HH '(1) :HT '(0 1) :TH '(0 0 1) :TT '(0 0 0)})

(defn encoder [code stream]
  (mapcat code stream))

(take 20 (encoder code stream)) ;(1 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 1 1 0 0)

;; Trying the two together:

(take 20  (decoder code-tree (encoder code stream))) ;(:HH :HH :HH :HH :HH :HH :HH :HH :HT :HH :HH :TH :HH :HH :HH :TT :HH :HH :HT :HT)

;; And finally:

(defn make-encoder [code]  (fn [s] (encoder code s)))
(defn make-decoder [code-tree] (fn[s] (decoder code-tree s)))

(cost (make-encoder code) (make-decoder code-tree) (take 10000 stream)) ; £16992

;; It costs us £16992 to send 10000 symbols.

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

;; We can use our code to send the output from the original biased coin very easily

(def unfair-coin {:H 3 :T 1})

(def unfair-stream (random-stream unfair-coin))

(take 20 unfair-stream) ; (:H :H :H :H :H :H :T :H :H :H :H :H :H :H :T :H :H :H :H :T)

(defn combine-keywords [& a] (keyword (apply str (mapcat #(drop 1 (str %)) a))))
(defn split-keyword [a] (map #(keyword (str %)) (drop 1 (str a))))


(defn make-combination-encoder [code n]
  (fn [s] (encoder code (map #(apply combine-keywords %) (partition n s)))))

(defn make-combination-decoder [code-tree]
  (fn [s] (mapcat split-keyword (decoder code-tree s))))

(cost (make-combination-encoder code 2) (make-combination-decoder code-tree) (take 10000 unfair-stream)) ; £8460

;; So our method of coding for {:HH 9, :HT 3, :TH 3, :TT 1} has given us a method of coding for {:H 3, :T 1}
;; which is 16% more efficient than the obvious one.

;; What if we try it on the output from the fair coin?

(def fair-stream (random-stream {:H 1 :T 1}))

(cost (make-combination-encoder code 2) (make-combination-decoder code-tree) (take 10000 fair-stream)) ; £ 11257

;; Using this code on the output from an unbiased coin actually makes it more expensive to transmit!

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

;; To recap:

;; We can transmit the output from a series of coin tosses, or other random
;; processes, down a binary channel, if we choose a code.

;; The code can be trivial, like :H -> 0 :T -> 1,
;; or it can be complex, like :HH -> 1, :HT -> 01, :TH ->001, :TT -> 000

;; Different codes can result in different costs of transmission for the outputs
;; of different processes

;; So far, we've seen costs of £1/symbol for the fair coin with the trivial code
;; or £1.12/symbol with the more complex code

;; And we've seen costs of £1/symbol and £0.84/symbol for the unfair coin with
;; the trivial and complex code respectively.

;; It seems that choosing the right code can make transmission cheaper, and
;; choosing the wrong code can make it more expensive.