Search This Blog

Wednesday, February 23, 2011

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.

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

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:

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.

Wednesday, February 16, 2011

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


(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))

;; 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] 
   (fn [x] (- (* x x) n))))

(sqrt 2) 

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

(defn cubert [n] 
   (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] 
   (fn[x] (- (pow x m) n))))

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

(pow 1.045350466874579 123) ;; 234.0000000243792

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

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

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

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

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

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

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

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

;; and start with an exact guess

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(def f (make-improver 25))

(f 1.0)                                 ;13.0

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

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

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

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

(iterative-improve 1.0 f g?)            ;5.000000000053722

;; So now we can find arbitrary square roots

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

(square-root 2)                         ;1.4142135623746899

Clojure Dojo: The method of Heron of Alexandria

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

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

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

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

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

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

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

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

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

              (* 5 3)
              (* 1 2))

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

;; What Heron said to do was this:

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

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

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

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

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

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

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

user> (/ 10 1)

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

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

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

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

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

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

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

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

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

user> (improve-guess (improve-guess 1))

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

; Barf....

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

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

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

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

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

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

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

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

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

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

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

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

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

user> (good-enough-guess 3.0)

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

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

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


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

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

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

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

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


Tuesday, February 8, 2011

RabbitMQ & Clojure: Hello World

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

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

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

;; $ sudo apt-get install rabbitmq-server

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

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

;; To install the latest rabbitmq release (2.3.1 at time of writing)
;; $ wget
;; $ sudo dpkg -i rabbitmq-server_2.3.1-1_all.deb

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

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

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

;; Which should respond:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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