2010年6月13日日曜日

ICFP 2009 orbit main

Here's the code that glues it all together and starts a simulation of the first Hohmann scenario:


;; Problem specification at http://www.ittc.ku.edu/icfp-contest/

(set! *warn-on-reflection* true)

(ns orbit-main
(:use orbit-vm orbit-binary-reader orbit-plotter orbit-controllers (clojure.contrib repl-utils))
(:import (java.io FileInputStream) (java.net URL)))

(def obf (obf-seq (byte-seq
;; (.. (URL. "http://www.ittc.ku.edu/icfp-contest/binaries/bin1.obf") openConnection getInputStream)
(FileInputStream. "bin1.obf") ; in case you downloaded the file, use this instead
)))

(comment do
(println)
(println "DATA OPCODE PARAM1 PARAM2")
(println "--------------------------------")
(let [fn-to-name { > ">" < "<" = "=" >= ">=" <= "<=" }]
(doseq [[data [opcode p1 p2]] (take 10 obf)]
(let [p1 (if (fn? p1) (fn-to-name p1) p1)]
(printf "%+08.2e %-6s %6s %6s%n" data opcode p1 p2)))))

(def vm (create-vm obf))

(def objects (atom {}))
(def stats (atom []))

(defn simulate
[vm controller]
(loop [vm vm countdown (long 3e6)] ; countdown can be used to set an upper bound for the number of iterations
;; (Thread/sleep 1) ; no sleep necessary on MY machine, running THIS code
(if (and (> countdown 0) (zero? (or (get-in vm [:outputs 0]) 0)))
(let [[new-inputs new-objects new-stats] (controller (:inputs vm) (:outputs vm))
vm (assoc-in vm [:inputs] new-inputs)]
(when (zero? (mod countdown 80)) ; throttle gui updates
(reset! stats new-stats)
(reset! objects new-objects))
(recur (one-simulation-step vm) (dec countdown)))
vm)))

(def orbit-panel (create-orbit-panel objects))
(def stats-panel (create-stats-panel stats))

(def frame (create-frame orbit-panel stats-panel))
(.setSize frame 500 500)

;; spawn another thread so I can continue using the repl while it spins
(.start (Thread. #(do
(def new-vm (simulate vm (create-hohmann-controller 1001)))
(:outputs new-vm))))

ICFP 2009 orbit controller

After setting up the infrastructure (binary reader, vm, visualizer), the controllers will be the part where the actual problems are solved.

For now, I won't solve any of the actual problems, I will just introduce the controller framework to show how it can be done. Here's the code:



(ns orbit-controllers
(:use orbit-plotter orbit-util)
(:import (java.awt Color)))

(def *earth-mass* 6e24)
(def *earth-radius* 6.357e6)

(def *the-earth* (struct orbit-object 0 0 *earth-radius* Color/BLUE))

(defn example-controller
"A controller that only initiliazes the scenario to 1003, returns some random stats and space objects,
but otherwise does nothing and is meant to act as an example of a minimal controller.

A controller is a function that gets called by the vm. The parameters are the inputs from last round
(inputs are what you communicate to the vm, like thrust values) and the outputs (they allow you to
read your x/y coordinates, fuel level, score etc.).
The controller has to return a vector
[new-inputs objects stats]
where new-inputs is a map of the vm input ports you want to write, objects is a name/orbit-object map
of space objects to be plotted by the visualizer, and stats is a list of key-value pairs to be displayed
in the stats table of the visualizer."

[inputs outputs]
(let [new-inputs (if (inputs 0x3e80) inputs (assoc inputs 0x3e80 1003))
objects { :earth *the-earth*
:random-satellite (struct orbit-object *earth-radius* *earth-radius* (/ *earth-radius* 32) Color/YELLOW) }
stats [[ "r" 1.234] ["x" 2.345] ["y" 3.456]]]
[new-inputs objects stats]))

(defn create-hohmann-controller
"Creates a controller for one of the Hohmann scenarios - this is just the framework needed to solve the problem.
So what's missing now is calculation of velocity/position/acceleration to calculate the (at least) 2 thrusts needed."

[scenario]
; close over state used inside the controller
(let [history (atom ())
memory (atom {})
last-update (atom nil)]
(fn [inputs outputs]
(let [delta-x (outputs 2 0) delta-y (outputs 3 0)
new-inputs (if-not (inputs 0x3e80) ; not yet initialized?
(assoc inputs 0x3e80 scenario) ; initialize scenario
(if (zero? (inputs 3 0)) ; currently we are not thrusting (in the y-axis)
(if-not (@memory :last-thrust-x) ; we haven't thrusted yet, so let's do it
(let [thrust-x 1000 thrust-y -1000]
(swap! memory assoc :last-thrust-x thrust-x :last-thrust-y thrust-y)
(assoc (assoc inputs 3 thrust-y) 2 thrust-x))
inputs) ; we have thrusted before - don't change
(assoc (assoc inputs 3 0) 2 0))) ; stop thrusting
now (System/nanoTime)
last @last-update
ups (if last (/ 1e9 (- now last)) -1)
objects { :earth *the-earth*, :satellite (struct orbit-object delta-x delta-y (/ *earth-radius* 64) Color/YELLOW) }
stats (partition 2 [ "r" (hypot delta-x delta-y)
"x" delta-x
"y" delta-y
"ups" ups])]
(reset! last-update now)
[new-inputs objects stats]))))

ICFP 2009 orbit visualizer

Visualization helps to better understand each of the problems. This is the visualization of the simplest problem, the Hohmann transfer:




And this screen shot shows the part around the earth in the "Operation Clear Skies" problem:



First some utilities we'll use:


(ns orbit-util)

(defmacro pipi [& body]
"quick hack debug prn macro that can be prepended before any fn/macro call without having to add any parens"
`(let [result# (~@body)]
(prn '~body result#)
result#))

(defn abs [n]
"Math/abs has caused me some trouble before - don't remember exactly, but this is what I always use now"
(if (< n 0) (* -1 n) n))

(defn hypot [a b]
"For some reason faster than Math/hypot"
(Math/sqrt (+ (* a a) (* b b))))


And here's the visualization code:


(ns orbit-plotter
(:use orbit-util)
(:import (java.awt BorderLayout Color Graphics Polygon RenderingHints Toolkit)
(java.awt.geom Ellipse2D$Double)
(javax.swing JComponent JFrame JLabel JPanel)
(java.awt.event ActionListener MouseAdapter MouseListener MouseMotionListener MouseWheelListener)))

(defstruct orbit-object :x :y :radius :color)

(def border-layout-order (list BorderLayout/CENTER BorderLayout/EAST BorderLayout/SOUTH BorderLayout/WEST BorderLayout/NORTH))

(defn create-frame [& components]
"Create a JFrame and add the components to it, the first will be the center, the others fill the border clockwise
starting from EAST"

(let [frame (JFrame. "ORBIT")]
(doseq [[component position] (map list components border-layout-order)]
(.add frame #^java.awt.Component component position))
(.setDefaultCloseOperation frame JFrame/DISPOSE_ON_CLOSE)
(.show frame)
frame))

(defn create-orbit-panel
"Create the main panel, showing the objects in space"
[objects-ref]
(let [scale (atom Long/MAX_VALUE)
panel
(proxy [JPanel] [(BorderLayout.)]
(paintComponent [#^java.awt.Graphics2D g2]
(let [panel #^javax.swing.JPanel this
panel-w (.getWidth panel) panel-h (.getHeight panel)]
(.setColor g2 Color/BLACK)
(.fillRect g2 0 0 panel-w panel-h)
(.setRenderingHint g2 RenderingHints/KEY_ANTIALIASING RenderingHints/VALUE_ANTIALIAS_ON)
(let [[max-rx max-ry]
(reduce (fn [[x y] [name obj]]
(let [r (:radius obj)]
[(max x (+ (abs (:x obj)) r)) (max y (+ (abs (:y obj)) r))])) [0 0] @objects-ref)
x-scale (if (zero? max-rx) 1.0 (/ panel-w (+ max-rx max-rx)))
y-scale (if (zero? max-ry) 1.0 (/ panel-h (+ max-ry max-ry)))
s (swap! scale #(min % x-scale y-scale))]
(.translate g2 #^java.lang.Double (+ (/ panel-w 2)) #^java.lang.Double (+ (/ panel-h 2)))
(.scale g2 s s)
)
(doseq [[name object] @objects-ref]
(let [r (:radius object)
dia (+ r r)]
(.setColor g2 (:color object))
(.fill g2 (Ellipse2D$Double. (- (:x object) r) (- (:y object) r) dia dia)))))))]
(add-watch objects-ref "orbit-panel"
(fn [key ref old-state new-state]
(if (not= old-state new-state)
(.repaint panel))))
panel))


(defn- generate-html-table
"Generate an HTML table for the statistics"
[stats]
(str (apply str
(reduce (fn [strs [k v]]
(let [v (if (= Double (class v)) (format "%+08.2e" v) v)]
(conj strs (format " <tr><td>%s</td><td align='right'>%s</td></tr>%n" k v)))) [ "<html><table>\n" ] stats))
"</table></html>"))

(defn create-stats-panel
"Create the stats panel, showing a table of the stats"
[stats-ref]
(let [label (doto
(javax.swing.JLabel.)
(.setOpaque true)
(.setBackground java.awt.Color/BLACK)
(.setForeground java.awt.Color/WHITE))]
(add-watch stats-ref "stats-panel"
(fn [key ref old-state new-state]
(if (not= old-state new-state)
(.setText label (generate-html-table new-state)))))
label))

2010年6月12日土曜日

ICFP 2009 orbit vm

The orbit virtual machine - pretty simple, as the binary reader took care of all the idiosyncracies of the binary format already.

The input and output ports are implemented as maps, as the 14-bit addressable memory is only sparsely used, allocating vectors/arrays for that seemed like a bit of a waste.


(ns orbit-vm)

(defstruct vm-struct :data :instructions :program-counter :inputs :outputs :status)

(defn create-vm
"Creates an orbit vm structure, which is used to hold the vm state and an instruction list"
[obf-seq]
(let [[data instructions] (apply map vector obf-seq)]
(struct vm-struct data instructions 0 {} {} 0)))

(defn process-instruction
"Processes the given (single) instruction on the vm. Returns a vm with the changed state"
[vm [opcode p1 p2]]
(let [r (fn [addr] (get-in vm [:data addr]))
w (fn [val] (assoc-in vm [:data (:program-counter vm)] val))
binary-op (fn [op] (w (op (r p1) (r p2))))
vm (condp = opcode
; d-instructions: p1: r1, p2: r2
:add (binary-op +)
:sub (binary-op -)
:mult (binary-op *)
:div (binary-op #(if (zero? %2) 0.0 (/ %1 %2)))
:output (do (assoc-in vm [:outputs p1] (r p2)))
:phi (w (r (if (zero? (:status vm)) p2 p1)))

; s-instructions: p1: "immediate", p2: r1
:noop vm
:cmpz (assoc-in vm [:status] (if (p1 (r p2) 0.0) 1 0))
:sqrt (w (Math/sqrt (r p2)))
:copy (w (r p2))
:input (let [v (or (get-in vm [:inputs p2]) 0)] (do (w v))))]
(assoc-in vm [:program-counter] (inc (:program-counter vm )))))

(defn one-simulation-step
"Processes one full round of instructions. Returns a vm with the changed state"
[vm]
(let [vm (reduce process-instruction vm (:instructions vm))]
(assoc-in vm [:program-counter] 0)))

ICFP 2009 orbit binary file reader

In preparation for next weekend's ICFP contest, I decided to do some clojure (<= 1.1) exercises this weekend, using last year's problem.

I will not do the actual problems, but I will implement the complete infrastructure, i.e a binary reader, the vm, a binary writer for the output file, and a graphical visualizer, as that is what will probably help me most in terms of improving my clojure skills.

This post takes care of the orbit binary file (obf) reader. I tried to implement the reader the clojure way, i.e. using sequences and operations on sequences rather than using nio and working on arrays of binary data directly to squeeze out performance - optimizations later.

Usage of the reader is straightforward: You provide the function obf-parse with a sequence of bytes (the data type is not Bytes, but Integers that can take values from 0-255, as returned by Java's InputStreams) and you get a sequence of vectors each containing a data value (type Double) and instructions (that are another vector consisting of the opcode and 2 parameters).

Anyone who has done ICFP 2009 knows that the obf format has its quirks - it look contrived, but in professional reality, you would often wish you could work with formats as nice as this one. The obf reader tries to take care of all these quirks, so that the vm gets a nice and clean stream of data and instructions.

Here's a sample output of the 10 first parsed frames of a obf file:



DATA OPCODE PARAM1 PARAM2
--------------------------------
1.00e+00 :noop 0 0
0.00e+00 :copy 0 265
3.00e+01 :noop 0 0
0.00e+00 :noop 0 0
0.00e+00 :copy 0 248
0.00e+00 :sub 4 3
0.00e+00 :cmpz = 5
0.00e+00 :phi 2 1
0.00e+00 :sub 7 0
0.00e+00 :copy 0 263



I'm posting this in the hope that I might get some feedback and that it can be used as motivation and inspiration for anyone who considers participating in the contest next week. It would be great if we got some successful clojure contributions!

If there's anything that could improve readability or usability of the code, or performance improvements that do not turn it into a mess, I'd like to know about it. And of course if you find any bugs, shoot a comment!


(ns orbit-binary-reader)

(defn byte-seq
[input-stream]
"Turns an InputStream into a sequence of bytes, lazily read, closing the stream when the end is reached"
(let [a-byte (. #^java.io.InputStream input-stream read)]
(if (= a-byte -1)
(do (. #^java.io.InputStream input-stream close) nil)
(lazy-seq (cons a-byte (byte-seq input-stream))))))

(defn- bytes-to-integer
"Turns big-endian ordered bytes into a integer number (of type int, long, etc.)"
[bytes]
(reduce #(+ (bit-shift-left %1 8) %2) 0 bytes))

(defn- bytes-to-double
"Turns big-endian ordered bytes into a double"
[bytes]
(Double/longBitsToDouble (bytes-to-integer bytes)))

(defn- obf-byte-seq-normalize
"Turns the input byte-seq into a byte-seq where data-bytes and instruction-bytes are alternating,
thus removing the even/odd oddity; also changes it to big-endian order"

[byte-seq]
(let [step (fn step [byte-seq even]
(when (not (empty? byte-seq))
(let [properly-ordered
(if even
(concat (reverse (take 8 byte-seq)) (reverse (take 4 (drop 8 byte-seq))))
(reverse (take 12 byte-seq)))]
(if (not= 12 (count properly-ordered))
(throw (Exception. (format "Expected 12 bytes, but only %d bytes left in seq" (count properly-ordered)))))
(concat properly-ordered (lazy-seq (step (drop 12 byte-seq) (not even)))))))]
(step byte-seq true)))

(def opcode-map
{
2r10001 :add
2r10010 :sub
2r10011 :mult
2r10100 :div
2r10101 :output
2r10110 :phi

2r00000 :noop
2r00001 :cmpz
2r00010 :sqrt
2r00011 :copy
2r00100 :input
})

(def cmpz-operator-map
{
2r0000000000 < ; LTZ
2r0010000000 <= ; LEZ
2r0100000000 = ; EQZ
2r0110000000 >= ; GEZ
2r1000000000 > ; GTZ
})

(defn- obf-parse-instruction
"Parses an instruction of type long into a vector [opcode arg*]"
[instruction-long]
(let [d-opcode (bit-and (bit-shift-right instruction-long 28) 2r1111)
s? (zero? d-opcode)
; we unify d-opcodes and s-opcodes into one format, d-opcodes have the fifth bit set, s-opcodes do not
; this makes it easier to look up the opcodes in the map
opcode (opcode-map (if s? (bit-and (bit-shift-right instruction-long 24) 2r1111) (bit-or d-opcode 2r10000)))
param1 (bit-and (bit-shift-right instruction-long 14) (if s? 2r1111111111 2r11111111111111))
param1 (if (= :cmpz opcode) (cmpz-operator-map (bit-and 2r1110000000 param1)) param1)
param2 (bit-and instruction-long 2r11111111111111)]
[opcode param1 param2]))


(defn- obf-normalized-parse
"Parses a normalized obf byte sequence, i.e. big-endian, and data / instruction alternating into
a sequence of [data [opcode param1 param2]] vectors, where data is a double and instruction an integer (int or long)"

[normalized-byte-seq]
(map
#(let [[data instruction] (split-at 8 %)] [(bytes-to-double data) (obf-parse-instruction (bytes-to-integer instruction))])
(partition 12 normalized-byte-seq)))

(defn obf-seq
"Parses a raw obf byte sequence (like from an obf file) into a sequence of [data [opcode param1 param2]] vectors,
where data is a double and instruction an integer (int or long)"

[byte-seq]
(obf-normalized-parse (obf-byte-seq-normalize byte-seq)))

2010年6月8日火曜日

Multilayer perceptrons in Clojure - Part II (backpropagation of errors)


(defn feed-forward-for-backprop
"Feeds an input vector through the network layer by layer and returns all raw and activated output vectors for all layers"
[network input activation activation-prime]
(reduce
(fn [all-previous layer]
(let [[[previous] before-previous] (split-at 1 all-previous)
previous-outputs (cons 1 (seq (first previous))) ; (cons 1) adds the bias
all-previous (cons (cons previous-outputs (rest previous)) before-previous) ; we'll also store the bias in the outputs we return
unactivated-out
(map
(fn [unit] (reduce + (map * previous-outputs unit)))
layer)]
(cons (list (map activation unactivated-out) unactivated-out (map activation-prime unactivated-out)) all-previous)))
(list [input nil nil])
network))

(defn propagate-back
"Trains the network with the input vector / target vector pair, using the backpropagation algorithm."
[network input target activation activation-prime eta]
(let [[[last-act last-unact last-prime] :as outputs] (feed-forward-for-backprop network input activation activation-prime)
reversed-network (reverse network)
delta-o (map (fn [out prime target] (* prime (- out target))) last-act last-prime target)]
(second (reduce
(fn [[last-deltas new-network rest-count] [old-layer [acts unacts primes]]]
(let [old-layer-trans (apply map list old-layer)
deltas (and (pos? rest-count) (map (fn [prime weights] (* prime (reduce + (map * weights last-deltas)))) primes old-layer-trans))
new-weights (map (fn [unit last-delta] (map (fn [w act] (- w (* eta act last-delta))) unit acts)) old-layer last-deltas)]
[deltas (cons new-weights new-network) (dec rest-count)]))
[delta-o () (dec (count (rest outputs)))]
(map list reversed-network (rest outputs))))))

(defn shuffle
[coll]
(let [l (java.util.ArrayList. coll)]
(java.util.Collections/shuffle l)
l))

(defn train-network
[network inputs targets activation activation-prime eta]
(reduce (fn [network [input target]]
(propagate-back network input target activation activation-prime eta))
network
(map list inputs targets)))

(defn train-network-shuffle-data
[network inputs targets activation activation-prime eta]
(let [[inputs targets] (apply map list (shuffle (map list inputs targets)))]
(train-network network inputs targets activation activation-prime eta)))

(last (take 10 (iterate #(train-network-shuffle-data % binary-inputs xor-targets (create-sigmoid 1) (create-sigmoid-prime 1) 1) a-zero-2-1-1)))
(let [network *1] (doseq [m (map #(feed-forward network %1 [(create-sigmoid 1) step-function]) xor-targets)] (prn m)) network)
(last (take 1000 (iterate #(train-network-shuffle-data % binary-inputs xor-targets (create-sigmoid 1) (create-sigmoid-prime 1) 1) *1)))

2010年6月7日月曜日

Multilayer perceptrons in Clojure - Part I (feed-forward)


(use 'clojure.contrib.pprint)

(defn rand-weight
"Return a random number r, where -1/sqrt(n) < r < 1/sqrt(n)"
[n]
(let [rand-weight-max (/ 1 (Math/sqrt n))]
(- (rand (* 2 rand-weight-max)) rand-weight-max)))

(defn create-layer
"Creates a single network layer with the specified number of units and inputs (+ one for the bias) to these units"
[inputs units]
(map
(fn [_] (map (fn [_] (rand-weight inputs)) (range (inc inputs))))
(range units)))

(defn create-network
"Creates a fully connected feed-forward neural network"
[& units-per-layer]
(map
(fn [[inputs units]] (create-layer inputs units))
(partition 2 1 units-per-layer)))

(defn feed-forward
"Feeds an input vector through the network layer by layer and returns the output vector"
[network input activations]
(reduce
(fn [previous-outputs [layer activation]]
(map
(fn [unit] (activation (reduce + (map * (cons 1 previous-outputs) unit)))) ; (cons 1) adds the bias
layer))
input
(map list network activations)))

(def binary-inputs [[0 0] [0 1] [1 0] [1 1]])

(def or-targets [[0] [1] [1] [1]])
(def and-targets [[0] [0] [0] [1]])
(def xor-targets [[0] [1] [1] [0]])

(defn create-sigmoid
[beta]
(fn [x] (/ 1 (+ 1 (Math/exp (* (- beta) x))))))

(defn create-sigmoid-prime
[beta]
(let [sigmoid (create-sigmoid beta)]
(fn [x] (let [sigmoid-x (sigmoid x)] (* beta (- 1 sigmoid-x) sigmoid-x)))))

(defn step-function
[x]
(if (> x 0) 1 0))

(def an-and-perceptron '(((-3 2 2))))
(def an-or-perceptron '(((-1 2 2))))
(def a-zero-perceptron '(((0 0 0))))

(feed-forward (create-network 2 1) [0 1] [step-function])

(do (prn) (pprint (map #(let [[t] %2 [y] (feed-forward an-and-perceptron %1 [step-function])] {:output y :target t :diff (- y t)}) binary-inputs and-targets)))

フォロワー