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

Hello Web : Dynamic Compojure Web Development at the REPL

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

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

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

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