Numerical Integration: Better Refinements?

;; 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?

;; 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

;; 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

;; 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?

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)

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

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 code easily.
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