Feynman' Arrows: OK, so What are the Complex Numbers?

;; Feynman's Arrows II : OK, so what are the complex numbers?
;; requires [simple-plotter "0.1.2"] 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Here's some code from the previous post

;; use pomegranate to get the library, if it's not already on your classpath
(require 'cemerick.pomegranate) 
(cemerick.pomegranate/add-dependencies 
 :coordinates '[[simple-plotter "0.1.2"]] 
 :repositories {"clojars" "http://clojars.org/repo"})

(use 'simple-plotter)

;; Make blackboards to draw arrows on
(defn make-blackboard [title size]
  (create-window title 400 400 white black (- size) size (- size) size)
  (axes)
  (ink yellow))

;; Draw an arrow shape from (a,b) to (c,d)
(defn draw-offset-arrow [[a b][c d]]
  (let [headx (+ a c) 
        heady (+ b d)]
    (line a b headx heady)
    (line headx heady (+ headx (* -0.086 c) (* -0.05 d)) (+ heady (*  0.05 c) (* -0.086 d)))
    (line headx heady (+ headx (* -0.086 c) (*  0.05 d)) (+ heady (* -0.05 c) (* -0.086 d)))))

;; Draw one of our arrows, which always have their tails at 0
(defn draw-arrow [[a b]] (draw-offset-arrow [0 0] [a b]))

;; Here's everything we know about the arrows so far:
(defn add-arrows[[a b][c d]]   
  [(+ a c) (+ b d)])

(defn multiply-arrows[[a b][c d]]
  [(- (* a c) (* d b)) (+ (* a d) (* c b))])

;; That's easier to read in the standard prefix notation
;; (a,b)+(c,d) -> (a+b, c+d)
;; (a,b)*(c,d) -> (ac-db, ad+cb)

;; The addition rule says: treat arrows as if they represented vector displacements, and add them nose to tail
;; The multiplication rule says: treat arrows as if they were zooms and rotations, and use one to zoom and rotate the other.

;; We also had a few favourite arrows that we'd played with:
(def arrow1 [3,4])
(def arrow2 [4,-3])
(def arrow3 [1,1/10])

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

(do (make-blackboard "Favourite Arrows" 6)
    (doseq [i [arrow1 arrow2 arrow3]] (draw-arrow i)))

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

;; The first thing I want you to notice is that there's a subset of
;; the arrows that works exactly like the numbers that we all got
;; comfortable with in primary school. They're the ones where the 
;; arrows point due east or west.

(add-arrows [10 0] [5 0]) ;-> [15 0]
;; (10,0) + (5,0) -> (15,0) is like 10 + 5 -> 15


(do (make-blackboard "Adding positive real numbers" 16)
    (draw-arrow [5 0])
    (draw-arrow [10 0])
    (ink red)
    (draw-arrow (add-arrows [10 0] [5 0])))
    


(multiply-arrows [4 0] [3 0]) ;-> [12 0]
;; (4,0) * (3,0) -> (12,0) is like 4 * 3 -> 12

(do (make-blackboard "Multiplying positive real numbers" 16)
    (draw-arrow [4 0])
    (draw-arrow [3 0])
    (ink red)
    (draw-arrow (multiply-arrows [4 0] [3 0])))



(multiply-arrows [-1 0] [-2 0]) ;-> [2 0]
;; (-1,0) * (-2,0) -> (2,0) is like   -1 * -2 -> 2

(do (make-blackboard "Multiplying negative real numbers" 3)
    (draw-arrow [-1 0])
    (draw-arrow [-2 0])
    (ink red)
    (draw-arrow (multiply-arrows [-1 0] [-2 0])))



;; You probably learned 'minus times minus is plus' as an arbitrary
;; rule, but it's obvious if you think of -1 as meaning 'leave it the
;; same size but turn it 180 degrees.'

;; We say that the 'real numbers' are 'embedded' in the arrows. 

;; What we mean is that the horizontal arrows behave just like the
;; real numbers, and so in any place where we were going to use real
;; numbers we can just use horizontal arrows instead and everything
;; will work out exactly the same.

;; And so since it doesn't matter whether we think about (3,0) or 3,
;; we'll usually just forget about the difference, and sometimes write (3,0),
;; and sometimes write 3, depending on convenience.

;; So now I can say -3 * 4 -> -12, and that's a statement about arrows!

;; There's another subset of the arrows, that point due north and due south.

;; Under addition, they're just like the real numbers too.

;; (0, 7)+ (0, 3) -> (0, 10)

(add-arrows [0 7] [0 3]) ;-> [0 10]

(do (make-blackboard "Adding Vertical Arrows" 12)
    (draw-arrow [0 7])
    (draw-arrow [0 3])
    (ink red)
    (draw-arrow (add-arrows [0 7] [0 3])))



;; But when you multiply them, they end up turning each other into horizontal arrows
;; (literally 'turning')

;; (0, 2) * (0, 3) -> (0x0-2x2, 0x3+2x0) = (-6, 0)

(multiply-arrows [0 2] [0 3]) ;-> [-6 0]

(do (make-blackboard "Multiplying Vertical Arrows" 12)
    (draw-arrow [0 2])
    (draw-arrow [0 3])
    (ink red)
    (draw-arrow (multiply-arrows [0 2] [0 3])))


;; You can see why pretty easily. 

;; (0, 10) represents 'turn 90 degrees clockwise and magnify by 10',
;; and (0, 4) means 'turn 90 degrees clockwise and magnify by 4'

;; And the product (-40, 0) means turn 180 degrees and multiply by 40. 

;; In this view, (-1, 0), or just -1 means 'turn 180 degrees (no zooming!)'

;; And the (0, 1) * (0, 1) -> (-1, 0) is just the fact 'if you turn 90
;; degrees anticlockwise, and then you turn another 90
;; degrees anticlockwise, then that's the same as if you'd
;; turned 180 degrees'.

;; So as long as we're talking about arrows, the thing we've called
;; -1, or the pair (-1,0), or the arrow length 1 that points east, or
;; the idea of turning through 180 degrees, does have a 'square root'.

;; There is a thing, the arrow length 1 that points straight north, or
;; the pair (0,1), or the idea of turning 90 degrees, that if you
;; multiply it by itself you get -1.

;; That's important, but it's also trivial.

;; Two quarter-turns clockwise make a 180 turn.

(multiply-arrows [0,1] [0,1]) ;-> [-1 0]

;; Two quarter-turns make a half-turn
(do (make-blackboard "Something Whose Square is (-1,0)" 2)
    (draw-arrow [0,1])
    (draw-arrow [0,1])
    (ink red)
    (draw-arrow (multiply-arrows [0,1] [0,1]))))

;; We call that upwards pointing, length 1 arrow i, for historical reasons. 
(do (make-blackboard "The Mysterious and Magical i" 2)
    (draw-arrow [0,1])))

;; And those vertical numbers, again, for historical reasons, are
;; called Imaginary Numbers.

;; Why on earth would anyone call sideways-pointing arrows 'real' and
;; upwards-pointing arrows 'imaginary'?

;; Well, it's hard to say exactly. The mathematicians of the 17th
;; century weren't thinking about rotations and scalings. They were
;; thinking about roots of equations.

;; Despite what's usually said, it never bothered them in the
;; slightest that x^2+1 didn't have a solution.  It was completely
;; obvious that it couldn't have, and that was fine.

;; What did worry them is that they had a formula for the cubic
;; equation, and when they used it to solve cubic equations that
;; really obviously did have three roots, like x^3=x, which is true
;; when x is -1, and when x is 0, and when x is 1, their formula kept
;; insisting that they do weird things, like taking the cube root of
;; the square root of -1.

;; And they had no idea what that meant, but some brave ones just ran
;; with the idea and found that if you followed the absurdity through
;; far enough then at the end the formula would give you the right
;; answers.

;; But they were, quite reasonably, very suspicious of the whole
;; procedure, and they called these weird things 'fictitious', or
;; 'imaginary' numbers.

;; I'm not sure why it was so important to solve cubic equations in
;; the 17th century. But they used to do it a lot, apparently. And
;; often the formula was the only way to get the answers.

;; And of course because they didn't really know what they were doing,
;; they got very confused and kept making mistakes, but nevertheless,
;; they did keep getting the right answers to their problems.

;; Eventually the breakthrough came in 1799, when a Danish
;; cartographer, Caspar Wessel, published 'On the Analytical
;; Representation of Direction', where he thought about using numbers
;; to represent directions and distances.

;; In this ground-breaking paper, Wessel both invented the idea of the
;; vector, and realized that the vectors that he was using to
;; represent directions and rotations and scalings were the same thing
;; as the 'imaginary numbers' that had been scaring and confusing
;; people for nearly two hundred years.

;; Unfortunately, "Om directionens analytiske betegning", was
;; published in Danish by the Royal Danish Academy of Sciences and
;; Letters, and so no-one could understand it, and so it vanished
;; without trace.

;; Luckily, seven years later, a bookstore manager in Paris,
;; Jean-Robert Argand, made the same discovery, and immediately
;; noticed that if you thought about arrows, and thought about
;; polynomial equations as being about arrows, then it was blindingly
;; obvious what was going on.

;; Indeed Argand realised that the question 
;; 'What are the roots of x^2 + 1 =0?'
;; is the same question as 'What do you have to do twice to turn round?'.

;; And he went quickly from that to showing that any such question has
;; an answer, which is called the Fundamental Theorem of Algebra.

;; And Argand published his idea in French, which is why today the
;; pictures of arrows that people use to reason about directions and
;; rotations and scalings are called 'Argand Diagrams'.

(do (make-blackboard "Wessel^h^h^h^h^h^h Argand Diagram" 10)
    (dotimes [i 100] (draw-arrow [(- (rand-int 20) 10) (- (rand-int 20) 10)])))



;; Anyway, even though things that point upwards aren't really any
;; more imaginary than things that point sideways, the name Imaginary
;; Numbers has stuck, and it has the twin virtues of:

;; (a) Making the Complex Numbers look Really Cool and Mysterious,
;; adding to the Aura of Mathematics.

;; (b) Scaring and Confusing Scientists and Engineers, and indeed
;; anyone who needs to think about things that rotate and get bigger
;; or smaller. Mathematicians love doing this. 

;; Engineers are practical people, and they use complex numbers all
;; the time when thinking about electricity and whether buildings will
;; fall over, and whether bridges will stay up, and eventually they
;; just accept that it's a really useful incomprehensible mystery and
;; get used to it. They say things like 'a capacitor is an imaginary
;; resistor', without even appearing to notice how strange that is.

;; But I think everyone else is puzzled and a bit frightened. I once
;; met a man with a PhD in Quantum Mechanics, who said that he
;; couldn't really believe that Quantum Mechanics was the true theory
;; of the world, because how could a thing that dealt with imaginary
;; quantities describe reality?

;; There is no need for any of this. Just Arrows. Just Rotations. Not
;; Scary. Children can do this.

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

;; Anyway, if 1 is another name for (1,0), and i is another name for
;; (0,1), then notice that 5+3i is (5,0)+(3,0)*(0,1) -> (5,3)

;; So we have yet another way of describing our pairs:

(defn print-arrow [[a b]]
   (str "the pair (" a "," b "), "
        "also known as the complex number " a "+"b"i, "
        "also known as the arrow " (cond (> a 0) (str a " north") (= a 0) "" :else (str (- a) "south"))
        " and " (cond (> b 0) (str b " east") (= b 0) "" :else (str (- b) " west"))))

(print-arrow arrow1) ;-> "the pair (3,4), also known as the complex number 3+4i, also known as the arrow 3 north and 4 east"
(print-arrow arrow2) ;-> "the pair (4,-3), also known as the complex number 4+-3i, also known as the arrow 4 north and 3 west"
(print-arrow arrow3) ;-> "the pair (1,1/10), also known as the complex number 1+1/10i, also known as the arrow 1 north and 1/10 east"
(print-arrow [1/11 0.25]) ;-> "the pair (1/11,0.25), also known as the complex number 1/11+0.25i, also known as the arrow 1/11 north and 0.25 east"

Feynman's Arrows : What are the Complex Numbers?

;; Feynman's Arrows : What are the Complex Numbers?

;; I've lost count of the number of times I've met physicists and even
;; professional mathematicians who think there's something a bit
;; spooky and unreal about the complex numbers.

;; There isn't.

;; Imagine that you're sitting in front of a photo-manipulation
;; program, and you're able to enlarge a photograph by a factor of
;; 2, and rotate it 30 degrees clockwise.

;; I don't think anyone is going to find that spooky. If you do, go
;; find a photo manipulation program and do it a couple of times.

;; If you do it twice, then you'll find that that gives you the same
;; effect as enlarging the photo by a factor of four and rotating it
;; 60 degrees clockwise.

;; In fact, how would we get our original photo back? My first guess
;; would be that I should shrink it by a factor of four and rotate it
;; 60 degrees anticlockwise. And that's the right answer. Go and try
;; it if it's not obvious.

;; And if you have your head round that, then you understand the
;; complex numbers.

;; If the mathematicians of the 16th century had been thinking about
;; how to use photoshop, instead of worrying about how to solve cubic
;; equations, then they'd have come up with the complex numbers in
;; about fifteen minutes flat, and they'd have thought they were the
;; most obvious thing in the world, and it would never have occurred
;; to them to talk about 'imaginary numbers', and an awful lot of
;; terror and confusion would have been avoided over the years.

;; And if the complex numbers had been found that way, then I think
;; we'd teach them to little children at about the same time we teach
;; them about fractions, and well before we teach them about the
;; really weird stuff like the square roots of two.

;; And the little children would have no problem at all with them, 
;; and they would think that they were fun, and easy.


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

;; In his wonderful little book QED, Richard Feynman managed to
;; explain the soul of Quantum Field Theory without using the complex
;; numbers at all.

;; Instead he talks about little arrows, which rotate like the hands
;; on a stopwatch, and occasionally he needs to add them up and rotate
;; them to work out what light and electrons are going to do.

;; Before we talk about the complex numbers, let's think about
;; Feynman's arrows.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; GRAPHICS CODE

;; I'd like to draw some pictures.  And I'd like to use my
;; 'simple-plotter' library, which is on clojars:
;; https://clojars.org/simple-plotter

;; This is probably the most complicated bit of this post, but it's
;; just graphics, and it doesn't matter in the slightest, so feel free
;; to skip over it unless you're interested in how to draw little
;; arrows on a computer.

;; You can add this:
;; [simple-plotter "0.1.2"] 

;; to your project file and restart everything, or if you are wise and
;; have the incomparable pomegranate on your classpath you can do:

(require 'cemerick.pomegranate)
(cemerick.pomegranate/add-dependencies 
 :coordinates '[[simple-plotter "0.1.2"]] 
 :repositories {"clojars" "http://clojars.org/repo"})

;; Either way, once the library is added to classpath
(use 'simple-plotter)

;; We can use it to draw arrows with tiny heads on them, so:
(defn draw-offset-arrow [[a b][c d]]
  (let [headx (+ a c) 
        heady (+ b d)]
    (line a b headx heady)
    (line headx heady (+ headx (* -0.086 c) (* -0.05 d)) (+ heady (*  0.05 c) (* -0.086 d)))
    (line headx heady (+ headx (* -0.086 c) (*  0.05 d)) (+ heady (* -0.05 c) (* -0.086 d)))))

;; And we'll usually think of our arrows as having their tails at [0,0]
(defn draw-arrow [[a b]] (draw-offset-arrow [0 0] [a b]))

(defn make-blackboard [title size]
  (create-window title 400 400 white black (- size) size (- size) size)
  (axes)
  (ink yellow))

;; END OF GRAPHICS CODE
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; So, with a few bits of graphics code so that the computer can do
;; what Feynman could do on a blackboard:

;; Let us make a blackboard:
(make-blackboard "Arrows!" 10)

;; This arrow points roughly northeast (3 miles east and 4 miles north, in fact)
(def arrow1 [3,4])

;; put it on the blackboard
(draw-arrow arrow1)

;; And this one points more southeast (4 miles east and -3 north, or 3 miles south)
(def arrow2 [4,-3])

;; put it on the blackboard
(draw-arrow arrow2) 

;; And we can add these arrows in what I'd hope is a really obvious way.
(defn add-arrows [arw1 arw2] (mapv + arw1 arw2))

(add-arrows arrow1 arrow2) ;-> [7 1]

;; Draw that too, but in red chalk
(do (ink red) 
    (draw-arrow (add-arrows arrow1 arrow2)))

;; This addition rule looks very simple in terms of coordinates, and geometrically it's 
;; very simple too. 

;; You draw one arrow (it's already there actually)
(do (ink yellow) (draw-arrow arrow1))

;; And then you draw the second one, but starting at the head of the first one:
(do (ink green) (draw-offset-arrow arrow1 arrow2))

;; And the sum is the arrow (in red) that points at the head of the
;; second green arrow.



;; It works the same whichever order you do the drawing in.
(make-blackboard "sworrA!" 10)
(draw-arrow arrow1)
(draw-arrow arrow2)
(ink green)
(draw-offset-arrow arrow2 arrow1)
(ink red)
(draw-arrow (add-arrows arrow2 arrow1))

;; It shouldn't be too hard to see how that's going to work for any
;; two arrows to produce a new one. And it should come as no surprise
;; that it works the same way whichever order you add the arrows in.

(add-arrows arrow2 arrow1) ;-> (7 1)
(add-arrows arrow1 arrow2) ;-> (7 1)

;; So the addition of arrows is fairly straightforward, and it gives a
;; way of talking about moving around in the plane. It works exactly
;; as it works for vectors.

;; Now we're going to do something slightly weirder, but in the same spirit.

;; Adding two arrows went like (a,b)+(c,d) -> (a+b, c+d)

;; But we'll say that to multiply two arrows, 

;; (a,b)*(c,d) -> (ac-db, ad+cb)

;; This rule looks a bit odd, but it turns out to be the way to think
;; about composing zooms and twists that I was talking about above,
;; that the ancient philosophers might have come up with if they'd
;; spent their days messing about with photoshop rather than dying of
;; plague or being set on fire by religious people.

;; At any rate it's a very easy rule to explain to a computer:
(defn multiply-arrows[[a b][c d]]
  [(- (* a c) (* d b)) (+ (* a d) (* c b))])

;; And, like arrow-addition, it's also symmetric.
;; Do a few by hand to see how it works!

;; Here's what we get if we multiply our two favourite arrows.
(multiply-arrows arrow1 arrow2) ;-> [24 7]
(multiply-arrows arrow2 arrow1) ;-> [24 7]

;; Here's another way of saying the same thing
(defn multiply-arrows[arrow1 arrow2]
  [(reduce - (map * arrow1 arrow2))
   (reduce + (map * arrow1 (reverse arrow2)))])

(multiply-arrows arrow1 arrow2) ;-> [24 7]
(multiply-arrows arrow2 arrow1) ;-> [24 7]

;; You might be able to see that there's a sense here of multiplying things, 
;; and of exchanging x coordinates and y coordinates, and taking pluses to minuses.

;; You might notice above that arrow2 [4, -3] is arrow1 [3,4] rotated by 90 degrees,
;; and that to get one from the other (3,4) -> (-4, 3), we exchange x and y and
;; swap +4 for -4.

;; In fact, we've encoded the idea of 'rotating and zooming' into our multiplication law.

;; Let's draw our example arrows and what they multiply to

(do 
  (make-blackboard "Arrow1 * Arrow1" 30) ;; <- Much bigger blackboard needed!
  (ink yellow)
  (draw-arrow arrow1)
  (draw-arrow arrow2)
  (ink red)
  (draw-arrow (multiply-arrows arrow1 arrow2)))

;; It turns out that the product of [3,4] and [4,-3] is a gigantic arrow.

(multiply-arrows arrow1 arrow2) ;-> [24 7]


;; The two arrows we started with were five long, but the result of multiplying them is 25 long.
;; Of course, we could have guessed that. 500% zoom followed by 500% zoom is 2500% zoom!


;; What happens, if, instead of taking arrow1 and arrow2 we multiply
;; arrow1 by a much smaller arrow, so that it will not grow so much?

;; Arrow 3 will be a tiny arrow, just a bit longer than 1 mile long, and 
;; just a bit to the north of due east.
(def arrow3 [1, 1/10])

(do 
  (make-blackboard "Arrow3 * Arrow1" 10)
  (ink yellow)
  (draw-arrow arrow1)
  (ink cyan)
  (draw-arrow arrow3)
  (ink red)
  (draw-arrow (multiply-arrows arrow3 arrow1)))




;; Arrow 3 represents a very small turn anticlockwise, and a very small zoom




;; Look at what happens if we keep multiplying arrow1 by arrow3

(ink red)
(draw-arrow (multiply-arrows arrow3 arrow1))
(draw-arrow (multiply-arrows arrow3 (multiply-arrows arrow3 arrow1)))
(draw-arrow (multiply-arrows arrow3 (multiply-arrows arrow3 (multiply-arrows arrow3 arrow1))))
;; and so on
(def arrows (iterate (fn[x] (multiply-arrows arrow3 x)) arrow1))
(draw-arrow (nth arrows 3))
(draw-arrow (nth arrows 4))
(draw-arrow (nth arrows 5))
(doseq [ i (take 100 arrows)]
  (draw-arrow i))

;; By repeatedly applying very small turns and very small zooms, we've
;; made a beautiful expanding spiral.


;; And I think, I genuinely think, that you might be able to get small
;; children to play with these arrows and to make this pretty spiral, just
;; after you've taught them about fractions.

The Triumph of the Median

;; The Triumph of the Median

;; http://lesswrong.com/lw/isk/a_voting_puzzle_some_political_science_and_a_nerd/
;; And originally from Sideways Stories from Wayside School by Louis Sachar

;; The students have Mrs. Jewl's class have been given the privilege
;; of voting on the height of the school's new flagpole. She has each
;; of them write down what they think would be the best height for the
;; flagpole. The votes are distributed as follows:

;;     1 student votes for 6 feet.
;;     1 student votes for 10 feet.
;;     7 students vote for 25 feet.
;;     1 student votes for 30 feet.
;;     2 students vote for 50 feet.
;;     2 students vote for 60 feet.
;;     1 student votes for 65 feet.
;;     3 students vote for 75 feet.
;;     1 student votes for 80 feet, 6 inches.
;;     4 students vote for 85 feet.
;;     1 student votes for 91 feet.
;;     5 students vote for 100 feet.

;; At first, Mrs. Jewls declares 25 feet the winning answer, but one
;; of the students who voted for 100 feet convinces her there should
;; be a runoff between 25 feet and 100 feet. In the runoff, each
;; student votes for the height closest to their original answer. But
;; after that round of voting, one of the students who voted for 85
;; feet wants their turn, so 85 feet goes up against the winner of the
;; previous round of voting, and the students vote the same way, with
;; each student voting for the height closest to their original
;; answer. Then the same thing happens again with the 50 foot
;; option. And so on, with each number, again and again, "very much
;; like a game of tether ball."

;; Question: if this process continues until it settles on an answer
;; that can't be beaten by any other answer, how tall will the new
;; flagpole be?

;; The most interesting part of this problem is to work out an incantation in emacs
;; that can turn

;;     1 student votes for 6 feet.
;;     1 student votes for 10 feet.
;;     7 students vote for 25 feet.
;;     1 student votes for 30 feet.
;;     2 students vote for 50 feet.
;;     2 students vote for 60 feet.
;;     1 student votes for 65 feet.
;;     3 students vote for 75 feet.
;;     1 student votes for 80 feet, 6 inches.
;;     4 students vote for 85 feet.
;;     1 student votes for 91 feet.
;;     5 students vote for 100 feet.

;; into

(def votes (concat
            (repeat 1,6)
            (repeat 1,10)
            (repeat 7,25)
            (repeat 1,30)
            (repeat 2,50)
            (repeat 2,60)
            (repeat 1,65)
            (repeat 3,75)
            (repeat 1,80.5)
            (repeat 4,85)
            (repeat 1,91)
            (repeat 5,100)))

;; which I leave as an exercise to the reader.

;; the lucky winner is 25 feet, with seven votes, and 100 feet getting 5 votes
(reverse (sort-by second (frequencies votes))) ;-> ([25 7] [100 5] [85 4] [75 3] [60 2] [50 2] [30 1] [91 1] [10 1] [6 1] [65 1] [80.5 1])

;; how do we do a run off between two options?

;; firstly let's ask whether a single voter would prefer option 1
(defn closer [vote option1 option2]
  (< (Math/abs (- option1 vote))
     (Math/abs (- option2 vote))))

(closer 26 25 100) ;-> true (one who likes 26 prefers 25 to 100)
(closer 50 25 100) ;-> true 
(closer 75 25 100) ;-> false (one who prefers 75 prefers 100 over 25)

;; given this choice between 100 and 25, how many prefer 25?
(frequencies (map (fn[x] (closer x 25 100)) votes)) ;-> {true 14, false 15}



;; So 100 is (just) stable against challenges from 25

;; But although 100 survived the challenge from 25, it still has to face 
;; a second challenge from partisans of 85

(defn winner [incumbent challenger votes]
  (let [v (frequencies (map (fn[x] (closer x incumbent challenger)) votes))]
    (if (< (v false 0) (v true 0))
      incumbent
      challenger)))


(winner 25 100 votes) ;-> 100
(winner 85 100 votes) ;-> 85

;; Now 50 gets its turn....

(winner 50 85 votes) ;-> 50

;; as do all the other attempts
(distinct votes) ;-> (6 10 25 30 50 60 65 75 80.5 85 91 100)

;; fifty can survive a lot of challenges but eventually 65 can beat it.
(reductions (fn [a b] (winner a b votes)) 50 (distinct votes))
 ;-> (50 50 50 50 50 50 60 65 65 65 65 65 65)

;; but 65 can survive any challenge
(reductions (fn [a b] (winner a b votes)) 65 (distinct votes))
 ;-> (65 65 65 65 65 65 65 65 65 65 65 65 65)

;; But it's not at all clear that this outcome is inevitable

(reduce (fn [a b] (winner a b votes)) (rand-nth votes) (shuffle votes)) ;-> 65

;; Careful examination of various possible run-offs shows that 
(for [a (distinct votes)]
  (for [b (distinct votes)]
  (winner a b votes)))

((6 10 25 30 50 60 65 75 80.5 85 91 100) 
 (10 10 25 30 50 60 65 75 80.5 85 91 100) 
 (25 25 25 30 50 60 65 75 80.5 85 91 100) 
 (30 30 30 30 50 60 65 75 80.5 85 91 100) 
 (50 50 50 50 50 60 65 75 50 50 50 50) 
 (60 60 60 60 60 60 65 60 60 60 60 60) 
 (65 65 65 65 65 65 65 65 65 65 65 65) 
 (75 75 75 75 75 60 65 75 75 75 75 75) 
 (80.5 80.5 80.5 80.5 50 60 65 75 80.5 80.5 80.5 80.5) 
 (85 85 85 85 50 60 65 75 80.5 85 85 85) 
 (91 91 91 91 50 60 65 75 80.5 85 91 91) 
 (100 100 100 30 50 60 65 75 80.5 85 91 100))


(count votes) ;-> 29

;; There's an easier way to get this answer: 
(def m (quot 29 2))
(nth (sort votes) m)

;; The votes will settle on the median. By symmetry I imagine that if
;; there were an even number of people there might be two stable
;; answers, and whichever got picked first would win forever.

;; This sort-of-explains why:

;; (a) Two ice-cream stalls will set up right next to each other in
;; the middle of a beach, rather than far apart

;; (b) Political parties in a democracy end up having very similar
;; views, which means that voters get no effective choice in an
;; election, whilst nevertheless getting a government that is
;; reasonably in tune with the popular mood.

2-SAT : Generalizing a Bit

;; 2-SAT : Generalizing a Bit

;; Let's generalize our problem a bit:

;; The instance of 2-SAT from last time:

;; Either x1 must be true or x2 must be true
;; Either x3 must be true or x4 must be true
;; Either x1 must be false or x3 must be false
;; Either x2 must be false or x4 must be false

;; in a more general form might look like:

(def problem [{:x1 true  :x2 true}
              {:x3 true  :x4 true}
              {:x1 false :x3 false}
              {:x2 false :x4 false}])

;; take an example clause 

(def example-clause {:x1 true :x2 true})

;; and a sample variable assignment

(def sample-assignment {:x1 true :x2 false :x3 true :x4 false})

;; and we can see whether the assignment satisfies either of the 
;; conditions in the clause

(for [[k v] example-clause] 
  (= (sample-assignment k) v)) ;-> (false true)

;; In this case, the clause says that either variable 1 must be true
;; (which is true) or variable 2 must be true (which is false).  Since
;; one of the clauses' conditions is satisfied, this clause overall is happy.


;; That allows us to make a function that tells us whether an
;; assignment satisfies a particular clause.
(defn satisfies [assignment clause]
  (some identity (for [[k v] clause] 
                   (= (assignment k) v))))


(satisfies sample-assignment example-clause) ;-> true

;; We can try each of the clauses in our example problem on our sample assignment
(satisfies sample-assignment {:x1 true  :x2 true}) ;-> true
(satisfies sample-assignment {:x3 true  :x4 true}) ;-> true
(satisfies sample-assignment {:x1 false :x3 false}) ;-> nil
(satisfies sample-assignment {:x2 false :x4 false}) ;-> true

;; And we find as before that it satisfies the first, second and
;; fourth clauses but not the third one.

;; Another way to do that would be:
(map (partial satisfies sample-assignment) problem ) ;-> (true true nil true)

;; And so a way to check whether they're all satisfied at once is:
(every? identity (map (partial satisfies sample-assignment) problem)) ;-> false

;; So here's a general mechanism to check whether an assignment is a
;; solution to a problem
(defn solves [assignment problem]
  (every? identity 
          (map (partial satisfies assignment) problem)))

(solves sample-assignment problem) ;-> false

;; Just to check, here's an assignment that isn't a solution
(solves {:x1 true :x2 false :x3 true :x4 false} problem) ;-> false
;; And here's one that is
(solves {:x1 true :x2 false :x3 false :x4 true} problem) ;-> true

;; How should we do our exhaustive search now?

;; Firstly find all the variables mentioned in the problem
(distinct (mapcat keys problem)) ;-> (:x2 :x1 :x4 :x3)

;; Then find a way of generating all possible assignments of true and false to 
;; those variables
(defn all-assignments [vars]
  (if (empty? vars) '({})
      (apply concat (for [a (all-assignments (rest vars))]
                      (list (assoc a (first vars) true)
                            (assoc a (first vars) false))))))
        
      

(all-assignments '()) ;-> ({})
(all-assignments '(:x1)) ;-> ({:x1 true} {:x1 false})
(all-assignments '(:x1 :x2)) ;-> ({:x1 true, :x2 true} {:x1 false, :x2 true} {:x1 true, :x2 false} {:x1 false, :x2 false})
(all-assignments '(:x1 :x2 :x3)) ;-> ({:x1 true, :x2 true, :x3 true} {:x1 false, :x2 true, :x3 true} {:x1 true, :x2 false, :x3 true} {:x1 false, :x2 false, :x3 true} {:x1 true, :x2 true, :x3 false} {:x1 false, :x2 true, :x3 false} {:x1 true, :x2 false, :x3 false} {:x1 false, :x2 false, :x3 false})

;; This gets back the two possible satisfying assignments that we worked out earlier:
(filter #(solves % problem) (all-assignments (distinct (mapcat keys problem))))
;; {:x2 true, :x1 false, :x4 false, :x3 true} 
;; {:x2 false, :x1 true, :x4 true, :x3 false}


(defn exhaustive-search-solve [problem]
  (filter #(solves % problem) 
          (all-assignments (distinct (mapcat keys problem)))))


;; We can now generate some test cases for future use:

;; If we just say that either x1 or x2 needs to be true
(exhaustive-search-solve
 [{:x1 true  :x2 true}]) 

;; Then there are three ways of doing that
;-> ({:x2 true, :x1 true} {:x2 false, :x1 true} {:x2 true, :x1 false})

;; But if we want one of them true and one of them false
(exhaustive-search-solve
 [{:x1 true  :x2 true} {:x1 false  :x2 false}]) 
;-> ({:x2 false, :x1 true} {:x2 true, :x1 false})

;; Then there are only two ways

;; Or we could specify that we want one with :x1 true
(exhaustive-search-solve
 [{:x1 true  :x2 true} {:x1 false  :x2 false} {:x1 true}]) 
;-> ({:x2 false, :x1 true})

;; Or even demand that x1 is true and x2 is true
(exhaustive-search-solve
 [{:x1 true  :x2 true} {:x1 false  :x2 false} {:x1 true} {:x2 true}]) 
;-> ()

;; That leads to an unsatifiable problem.

;; I've just realized that I've accidentally over-generalised this to
;; be able to represent the general problem of boolean satisfiability.

(exhaustive-search-solve
 [{:x1 false :x2 false} 
  {:x2 false :x3 false :x4 false :x5 true}
  {:x1 true} 
  {:x4 false :x1 false}
  {:x3 true :x2 false}
  {:x3 true :x4 true}
  {:x4 true :x5 true}
  ]) 

;-> ({:x2 false, :x1 true, :x4 false, :x5 true, :x3 true})

;; Good Old Clojure!

2SAT and P = NP

;; 2-SAT and P = NP

;; The archetypal NP-complete problem is 3-SAT. Given a set of logical
;; triples that must be enforced, is there a satisfying allocation of
;; variables?

;; To warm up for this, let's try its tractable cousin 2-SAT

;; In 2-SAT, the clauses to be satisfied are all of form either:
;; (a must be true/false) or (b must be true/false)

;; So here's an example problem, with four variables and four clauses

;; Either x1 must be true or x2 must be true
;; Either x3 must be true or x4 must be true
;; Either x1 must be false or x3 must be false
;; Either x2 must be false or x4 must be false

;; A solution to the problem is either a satisfying assignment, or a
;; proof that no such assignment exists.

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

;; To make this trial problem more precise:

;; Represent an assignment as a vector of true, false values, then 
(def test-assignment [ true false true false ])

;; We can represent the first clause as:
(defn satisfies-1 [v] (or (nth v (dec 1) (nth v (dec 2)))))

(satisfies-1 test-assignment) ;-> true  ;; So far so good

;; The second clause as:
(defn satisfies-2 [v] (or (nth v (dec 3)) (nth v (dec 4))))

(satisfies-2 test-assignment) ;-> true  ;; Hooray

;; The third clause as:
(defn satisfies-3 [v] (or (not (nth v (dec 1))) (not (nth v (dec 3)))))

(satisfies-3 test-assignment) ;-> false  ;; Darn

;; And finally the fourth clause as:
(defn satisfies-4 [v] (or (not (nth v (dec 2))) (not (nth v (dec 4)))))

(satisfies-4 test-assignment) ;-> true   ;; Three out of four ain't bad

;; So actually we only failed one of the clauses with our guess.

;; What it means for the problem to be in NP is that given a proposed
;; satisfying assignment, we can check all these conditions in
;; polynomial time. (polynomial in n, the length of the input, here
;; take n to be the greater of the number of clauses and variables)

;; In fact it's obvious that it's linear in n. We need to check each
;; clause exactly once on a proposed solution.

(defn satisfies? [v]
  (and (satisfies-1 v)
       (satisfies-2 v)
       (satisfies-3 v)
       (satisfies-4 v)))

;; So 2-SAT is definitely in NP.

(satisfies? test-assignment) ;-> false

;; Given speedy checking of proposed solutions, we can find (or rule out) a satisfying solution by brute-force search:

(def possible-assignments 
  (for [i [true false]
        j [true false]
        k [true false]
        l [true false]] [i j k l]))

possible-assignments ;-> ([true true true true] [true true true false] [true true false true] [true true false false] [true false true true] [true false true false] [true false false true] [true false false false] [false true true true] [false true true false] [false true false true] [false true false false] [false false true true] [false false true false] [false false false true] [false false false false])

(count possible-assignments) ;-> 16  ;; Not too bad, we can just try them all:

(for [a possible-assignments :when (satisfies? a)] a) ;-> ([true false false true] [false true true false])

;; So we've found two possible satisfying assignments to the variables

(satisfies? [true false false true]) ;-> true
(satisfies? [false true true false]) ;-> true

;; One would have done. And if we'd found none, then that would have
;; shown the problem to be unsatisfiable, which is also a valid answer
;; to this 'decision problem'.

;; Because the space of possible solutions is exponential in the
;; number of variables, which is constrained by the size of the input,
;; the fact that 2-SAT is in NP and so each trial solution can be
;; assessed in polynomial time, guarantees us an exponential-time
;; solution to the problem.

;; The interesting question is 'can we do better than that?'

;; In the case of 2-SAT we can do much better. In fact it will turn
;; out that 2-SAT can be solved in (small) polynomial time, which
;; means that we can solve very large problems of this kind.

;; In the case of 3-SAT, almost certainly not. There's a vast class of
;; problems in NP which are all easier than 3 SAT, in the sense that
;; if we had a polynomial time 3-SAT algorithm then we'd automatically
;; also have polynomial time algorithms for this vast class, and
;; people have been trying to solve them for many years. And nobody
;; has ever come up with a polynomial time algorithm for any of them. 

;; This is known as the P=NP problem. It's thought not, but nobody knows.
;; It's probably the most important open maths problem there is.

Efficiency and Progress V: Recap and Best Code so Far

;; Efficiency and Progress V:

;; Recap

;; This is slow:
(time (reduce + (map + (range 1000000) (range 1000000))))

;; C can perform roughly equivalent task in 8.6 ms, and Java in around 16 ms

;; So do we have to drop into C or Java when we want to make algorithms fast?
;; I hope not! 

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

;; I'm starting clojure like this, so that I know it's running as fast as it can
;; rlwrap java -server -classpath ~/.m2/repository/org/clojure/clojure/1.5.1/clojure-1.5.1.jar:. clojure.main

(clojure-version) ;-> "1.5.1"
((into{} (System/getProperties)) "java.version") ;-> "1.7.0_25"
((into{} (System/getProperties)) "java.vm.name") ;-> "OpenJDK Server VM"
((into{} (System/getProperties)) "sun.management.compiler") ;-> "HotSpot Tiered Compilers"

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

;; Repeated runs through this little benchmark now show 
;; Hotspot doing its thing and speeding up the calculation
(time (reduce + (map + (range 1000000) (range 1000000))))
999999000000
"Elapsed time: 2852.570643 msecs"
"Elapsed time: 2756.410014 msecs"
"Elapsed time: 1907.089513 msecs"
"Elapsed time: 1872.189534 msecs"
"Elapsed time: 1870.054495 msecs"
"Elapsed time: 1907.627285 msecs"

;; Setting these two variables is a good thing when trying to achieve C/Java like speeds
(set! *warn-on-reflection* true)
(set! *unchecked-math* true)

;; And doesn't seem to make any difference to this code
(time (reduce + (map + (range 1000000) (range 1000000))))
999999000000
"Elapsed time: 1883.237194 msecs"

;; So I think it's safe to conclude that Clojure written idiomatically
;; is around 100x slower than Java And around 200x slower than C

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

;; I have a little microbenchmark written in C which represents the
;; sort of things I am trying to do by adding length 1000000 vectors
;; and then adding up all the numbers in the vectors repeatedly.

;; Appropriately compiled, it runs in 8.6 seconds

;; Approximately the same program, translated into Java, runs in
;; around 16 seconds with java7 -server

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

;; In idiomatic Clojure I reckon this program would take about 1600
;; seconds, or about half an hour. And that's actually reflective of
;; the difference between C and Clojure that I've observed while
;; trying to write algorithms code as part of the Coursera/Stanford
;; Algorithms II class (which I very highly recommend!)

;; My best shot in non-idiomatic Clojure in the last post was:

(def N 1000000)
(def a (int-array (range N)))
(def b (int-array N))

(time 
 (let [a ^ints a b ^ints b N ^int N]
   (loop [count 0 sum 0]
     (if (= count 1000) sum
         (do 
           (println count sum)
           (dotimes [i N] (aset b i (+ (aget a i)(aget b i))))
           (recur (inc count) (+ sum (loop [i 0 ret 0] 
                                       (if (= i N) ret
                                           (recur (unchecked-inc i) (+ ret (aget b i))))))))))))

"Elapsed time: 177181.749304 msecs"
250249749750000000

;; Which is unreadable and still about 11x slower than Java, but gets
;; the right answer

;; I am so unused to using mutation in Clojure that I keep forgetting
;; to reset the variables and then being surprised when the answers
;; are wrong. It really screws up the REPL way of programming. No
;; wonder LISP developed the functional style that's now becoming so
;; fashionable.



;; After an awful lot of essentially random screwing around, I managed
;; to concoct another version, which uses the areduce macro instead of
;; the explicit inner loop above.
(def N 1000000)
(def a (int-array (range N)))
(def b (int-array N))


(time 
 (let [a ^ints a b ^ints b]
   (loop [count 0 sum 0]
     (if (= count 1000) sum
         (do 
           (println count sum)
           (dotimes [i N] (aset b i (+ (aget a i)(aget b i))))
           (recur (inc count) (+ sum (areduce b i ret 0 (+ ret (aget b i))))))))))

"Elapsed time: 63657.893856 msecs"
250249749750000000

;; I have no idea why this runs so much faster. A profiler would be
;; very useful here, but jvisualvm, which I used to find really useful
;; when tuning clojure, is giving me no real information and wasting
;; an awful lot of time in return. It only seems to profile at the
;; level of Java classes, and so this loop code isn't visible to it.

;; Still, down to about 4x slower than Java, maybe 7x slower than C,
;; and (I think) quite a lot more readable.



;; Various commenters have been kind enough to suggest improvements over my pitiful effort.

;; Dmitry Groshev's:

(def N 1000000)
(def a (int-array (range N)))
(def b (int-array N))

(defn test3 []
  (let [^ints a a
        ^ints b b
        N (int N)]
    (loop [count (int 0) sum (long 0)]
      (if (== count 1000) sum
          (do
            (loop [i (int 0)]
              (when (< i N)
                (aset b i (+ (aget a i) (aget b i)))
                (recur (inc i))))
            (recur (inc count)
                   (+ sum (long (loop [i (int 0) ret (long 0)]
                                  (if (== i N) ret
                                      (recur (inc i)
                                             (+ ret (aget b i)))))))))))))


(time (test3))

"Elapsed time: 45213.470027 msecs"
250249749750000000

;; Gets down to 3x slower than Java, 5x slower than C




;; And James Reeves, by private e-mail after he couldn't comment on this blog

(defn asum [^ints xs]
  (areduce xs i s 0 (unchecked-add s (aget xs i))))

(defn amap-add [^ints xs ^ints ys]
  (dotimes [i (alength xs)]
    (aset xs i (unchecked-add (aget xs i) (aget ys i)))))

(defn test-low-level []
  (let [a (int-array 1000000)
        b (int-array 1000000)]
    (dotimes [i (alength a)]
      (aset a i i))
    (loop [count 0, sum 0]
      (amap-add b a)
      (if (< count 1000)
        (recur (unchecked-inc count) (unchecked-add sum (asum b)))
        sum))))

(time (test-low-level))
"Elapsed time: 46496.528226 msecs"
250249749750000000



;; Has managed to get down to the same speed, whilst splitting the
;; inner loops off into nice little functions!

;; I think this version wins hands down for general comprehensibility.

;; So well done James, this is this blog's collective best shot at
;; this problem so far.

;; We're down to 3x slower than Java, 5x slower than C, and readable
;; if not quite as readable as it would be in a language designed for
;; imperative loops over arrays.

;; If I can learn to produce code like this reliably and without too
;; much buggering about, then I may be able to stay in Clojure for the
;; tight loops parts of my algorithms code. A factor of 5 I can
;; tolerate for the convenience of doing everything else in Clojure.

;; As Dmitry pointed out, it may be possible to use macros to make a
;; little language to make this easier

;; And James shows us another possible route with functions.

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

;; I still don't understand why these similar looking bits of code
;; (things like areduce are only macros that write code much like the
;; explicit loops) run at such different speeds, and that worries me,
;; because that's going to translate to a lot of screwing around while
;; trying to write such things.

;; And also, why even after all this are we not as fast as Java?

;; Bernard said something about 'removing synchronization', which
;; sounds scary.

;; Another approach, of course, would be to write pure Java classes to
;; do the heavy lifting and call them from Clojure code which does the
;; pre-processing, but that sounds like a fairly nasty and fragile
;; approach itself, and I'd rather avoid doing that if I can. But it
;; should be a bit easier than using C within python, at least.

;; I should also investigate core.matrix and hiphip, two projects
;; targeted at this sort of thing.

;; If anyone knows how to use them to solve this problem, please chip in.