From 26a1088f0890ec3348edabc22bf4168c82f7e91a Mon Sep 17 00:00:00 2001 From: John Lawrence Aspden Date: Mon, 20 Dec 2010 15:00:13 +0000 Subject: [PATCH] hamming codes --- itila1.4.clj | 178 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 178 insertions(+) create mode 100644 itila1.4.clj diff --git a/itila1.4.clj b/itila1.4.clj new file mode 100644 index 00000000..107bca0f --- /dev/null +++ b/itila1.4.clj @@ -0,0 +1,178 @@ +;; Information Theory, Inference, and Learning Algorithms +;; Exercise 1.4 + +;; A binary symmetric channel of noise f flips each bit with probability f + +(defn flip[bit] (if (= bit 0) 1 0)) + +(defn flip-rand[f bit] (if (< (rand) f) (flip bit) bit)) + +(defn corrupt[f sq] (map #(flip-rand f %) sq)) + +(defn corrupt[f sq] + (lazy-seq + (when-let [sq (seq sq)] + (cons (if (< (rand) f) (flip (first sq)) (first sq)) + (corrupt f (rest sq)))))) + + +(defn binary-symmetric-channel[f] (fn[sq] (corrupt f sq))) + +(defn estimate-bit-error [coder decoder n] +; (with-precision 15 +; (bigdec + (let [channel (binary-symmetric-channel 0.1) + message-in (repeat n (rand-int 2)) + tx (coder message-in) + rx (channel tx) + message-out (decoder rx) + ferrs (frequencies (map bit-xor message-in message-out))] + (/ (get ferrs 1 0) (+ (get ferrs 0 0) (get ferrs 1 0))))) + + + + +;; repetition codes just repeat each bit n times +(defn repetition-code [n sq] + (lazy-seq + (when-let [sq (seq sq)] + (concat (repeat n (first sq)) (repetition-code n (rest sq)))))) + +(defn repetition-code [n sq] (mapcat #(repeat n %) sq)) + + +(defn repetition-decode [n sq] + (lazy-seq + (when-let [sq (seq sq)] + (let [ssq (take n sq) + d (first (last (sort-by second (frequencies ssq))))] + (cons d (repetition-decode n (drop n sq))))))) + +(defn repetition-decode [n sq] + (map ffirst + (map #(sort-by second > %) + (map frequencies + (partition n sq))))) + +(defn repetition-coder [n] (fn [sq] (repetition-code n sq))) +(defn repetition-decoder [n] (fn [sq] (repetition-decode n sq))) + +;; these keep coming in handy. uncomment if needed +;; (defn factorial [n] (if (< n 2) 1N (* n (factorial (dec n))))) +;; (defn choose [ n r ] (/ (factorial n) (factorial r) (factorial (- n r)))) + + +;; Hamming's 7 4 code +(defn parity [ & sq] (reduce bit-xor sq)) + + +(defn hamming-7-4 [[a b c d]] + (let [e (parity a b c) + f (parity b c d) + g (parity c d a)] + [a b c d e f g])) + +(defn hamming-7-4-syndrome [[a b c d e f g]] + [(parity a b c e) (parity b c d f) (parity c d a g)]) + +(defn syndrome-to-bit [a b c] + (case [a b c] + [1 0 1] 0 + [1 1 0] 1 + [1 1 1] 2 + [0 1 1] 3 + 'ok)) +(defn hamming-7-4-correct [[a b c d e f g :as v]] + (def-let [syn (apply syndrome-to-bit (hamming-7-4-syndrome v))] + (if (= syn 'ok) [a b c d] + (vec (concat (take syn v) (cons (flip (nth v syn)) (drop (inc syn) (list a b c d)))))))) + + + +(defn hamming-7-4-encoder [sq] + (mapcat hamming-7-4 (partition 4 sq))) + +(defn hamming-7-4-decoder [sq] + (mapcat hamming-7-4-correct (partition 7 sq))) + + +;; Tests and sanity checks + +(use 'clojure.test) +(deftest corruption + (is (= ((binary-symmetric-channel 1) (repeat 10 1)) (repeat 10 0))) + (is (= ((binary-symmetric-channel 0) (repeat 10 1)) (repeat 10 1))) + (is (= (count (frequencies ((binary-symmetric-channel 0.5) (repeat 10 1)))) 2)) + + (is (= ((repetition-coder 3) (repeat 10 1)) (repeat 30 1))) + (is (= ((repetition-decoder 3) (repeat 30 0)) (repeat 10 0))) + (is (= ((repetition-coder 7) (repeat 10 1)) (repeat 70 1))) + (is (= ((repetition-decoder 7) (repeat 70 0)) (repeat 10 0))) + + (is (< 0.01 (estimate-bit-error (repetition-coder 3) (repetition-decoder 3) 1000) 0.1)) + (is (< (estimate-bit-error (repetition-coder 9) (repetition-decoder 9) 1000) 0.003)) + + (is (= (hamming-7-4 '(0 0 0 0)) [0 0 0 0 0 0 0])) + (is (= (hamming-7-4 '(0 0 0 1)) [0 0 0 1 0 1 1])) + (is (= (hamming-7-4 '(0 0 1 0)) [0 0 1 0 1 1 1])) + (is (= (hamming-7-4 '(0 0 1 1)) [0 0 1 1 1 0 0])) + (is (= (hamming-7-4 '(0 1 0 0)) [0 1 0 0 1 1 0])) + (is (= (hamming-7-4 '(0 1 0 1)) [0 1 0 1 1 0 1])) + (is (= (hamming-7-4 '(0 1 1 0)) [0 1 1 0 0 0 1])) + (is (= (hamming-7-4 '(0 1 1 1)) [0 1 1 1 0 1 0])) + (is (= (hamming-7-4 '(1 0 0 0)) [1 0 0 0 1 0 1])) + (is (= (hamming-7-4 '(1 0 0 1)) [1 0 0 1 1 1 0])) + (is (= (hamming-7-4 '(1 0 1 0)) [1 0 1 0 0 1 0])) + + + (is (= (hamming-7-4-correct [0 0 0 0 0 0 0]) [0 0 0 0])) + (is (= (hamming-7-4-correct [0 0 0 0 0 0 1]) [0 0 0 0])) + (is (= (hamming-7-4-correct [0 0 0 0 0 1 0]) [0 0 0 0])) + (is (= (hamming-7-4-correct [0 0 0 0 1 0 0]) [0 0 0 0])) + (is (= (hamming-7-4-correct [0 0 0 1 0 0 0]) [0 0 0 0])) + (is (= (hamming-7-4-correct [0 0 1 0 0 0 0]) [0 0 0 0])) + (is (= (hamming-7-4-correct [0 1 0 0 0 0 0]) [0 0 0 0])) + (is (= (hamming-7-4-correct [1 0 0 0 0 0 0]) [0 0 0 0])) + + (is (= (hamming-7-4-syndrome [0 0 0 0 0 0 0]) [0 0 0])) + (is (= (hamming-7-4-syndrome [1 0 0 0 0 0 0]) [1 0 1])) + (is (= (hamming-7-4-syndrome [0 1 0 0 0 0 0]) [1 1 0])) + (is (= (hamming-7-4-syndrome [0 0 1 0 0 0 0]) [1 1 1])) + (is (= (hamming-7-4-syndrome [0 0 0 1 0 0 0]) [0 1 1])) + (is (= (hamming-7-4-syndrome [0 0 0 0 1 0 0]) [1 0 0])) + (is (= (hamming-7-4-syndrome [0 0 0 0 0 1 0]) [0 1 0])) + (is (= (hamming-7-4-syndrome [0 0 0 0 0 0 1]) [0 0 1])) + + (is (= (syndrome-to-bit 0 0 0) 'ok)) + (is (= (syndrome-to-bit 1 1 1) 2)) + (is (= (syndrome-to-bit 0 1 1) 3)) + + (is (= (parity 1 1 1) 1)) + + (is (= (apply syndrome-to-bit (hamming-7-4-syndrome [0 0 0 0 0 0 0])) 'ok)) + (is (= (apply syndrome-to-bit (hamming-7-4-syndrome [0 0 0 1 0 0 0])) 3)) + + (is (= (hamming-7-4-encoder '(0 0 0 0 0 1 0 1 1 1 1 1)) + '(0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 1 1 1 1 1 1))) + (is (= (hamming-7-4-decoder (hamming-7-4-encoder '(0 1 0 1 0 0 1 0 1 0 1 1 ))) + '(0 1 0 1 0 0 1 0 1 0 1 1 ))) + + + + + + + ) +(run-tests) + + +(def errrates + (map + #(estimate-bit-error (repetition-coder %) (repetition-decoder %) 10000) + '(3 5 7 9 11))) ; (0.0272M 0.0088M 0.0026M 0.001M 0.0004M) + +;; Theory predicts 0.36 = 4*0.1*0.9 asymptotically +(map float (map / (drop 1 errrates) errrates)) ;; (0.372 0.290 0.37 0.5) + + +(float (estimate-bit-error hamming-7-4-encoder hamming-7-4-decoder 10000))