Mercurial > dylan

annotate sicm/bk/utils.clj @ 2:b4de894a1e2e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
initial import
author Robert McIntyre <rlm@mit.edu>
date Fri, 28 Oct 2011 00:03:05 -0700
parents
children
rev   line source
rlm@2 1
rlm@2 2 (ns sicm.utils)
rlm@2 3
rlm@2 4 ;***** GENERIC ARITHMETIC
rlm@2 5 (ns sicm.utils)
rlm@2 6 (in-ns 'sicm.utils)
rlm@2 7
rlm@2 8 (defprotocol Arithmetic
rlm@2 9 (zero [this])
rlm@2 10 (one [this]))
rlm@2 11
rlm@2 12
rlm@2 13 (extend-protocol Arithmetic
rlm@2 14 java.lang.Number
rlm@2 15 (zero [this] 0)
rlm@2 16 (one [this] 1))
rlm@2 17
rlm@2 18 (extend-protocol Arithmetic
rlm@2 19 clojure.lang.Seqable
rlm@2 20 (zero [this] (map zero this))
rlm@2 21 (one [this] (map one this)))
rlm@2 22
rlm@2 23
rlm@2 24 ;***** TUPLES AND MATRICES
rlm@2 25 (in-ns 'sicm.utils)
rlm@2 26
rlm@2 27 (defprotocol Spinning
rlm@2 28 (up? [this])
rlm@2 29 (down? [this]))
rlm@2 30
rlm@2 31 (defn spin
rlm@2 32 "Returns the spin of the Spinning s, either :up or :down"
rlm@2 33 [#^Spinning s]
rlm@2 34 (cond (up? s) :up (down? s) :down))
rlm@2 35
rlm@2 36
rlm@2 37 (deftype Tuple
rlm@2 38 [spin coll]
rlm@2 39 clojure.lang.Seqable
rlm@2 40 (seq [this] (seq (.coll this)))
rlm@2 41 clojure.lang.Counted
rlm@2 42 (count [this] (count (.coll this))))
rlm@2 43
rlm@2 44 (extend-type Tuple
rlm@2 45 Spinning
rlm@2 46 (up? [this] (= ::up (.spin this)))
rlm@2 47 (down? [this] (= ::down (.spin this))))
rlm@2 48
rlm@2 49 (defmethod print-method Tuple
rlm@2 50 [o w]
rlm@2 51 (print-simple (str (if (up? o) 'u 'd) (.coll o)) w))
rlm@2 52
rlm@2 53
rlm@2 54
rlm@2 55 (defn up
rlm@2 56 "Create a new up-tuple containing the contents of coll."
rlm@2 57 [coll]
rlm@2 58 (Tuple. ::up coll))
rlm@2 59
rlm@2 60 (defn down
rlm@2 61 "Create a new down-tuple containing the contents of coll."
rlm@2 62 [coll]
rlm@2 63 (Tuple. ::down coll))
rlm@2 64
rlm@2 65
rlm@2 66 (in-ns 'sicm.utils)
rlm@2 67
rlm@2 68 (defn numbers?
rlm@2 69 "Returns true if all arguments are numbers, else false."
rlm@2 70 [& xs]
rlm@2 71 (every? number? xs))
rlm@2 72
rlm@2 73 (defn contractible?
rlm@2 74 "Returns true if the tuples a and b are compatible for contraction,
rlm@2 75 else false. Tuples are compatible if they have the same number of
rlm@2 76 components, they have opposite spins, and their elements are
rlm@2 77 pairwise-compatible."
rlm@2 78 [a b]
rlm@2 79 (and
rlm@2 80 (isa? (type a) Tuple)
rlm@2 81 (isa? (type b) Tuple)
rlm@2 82 (= (count a) (count b))
rlm@2 83 (not= (spin a) (spin b))
rlm@2 84
rlm@2 85 (not-any? false?
rlm@2 86 (map #(or
rlm@2 87 (numbers? %1 %2)
rlm@2 88 (contractible? %1 %2))
rlm@2 89 a b))))
rlm@2 90
rlm@2 91
rlm@2 92
rlm@2 93 (defn contract
rlm@2 94 "Contracts two tuples, returning the sum of the
rlm@2 95 products of the corresponding items. Contraction is recursive on
rlm@2 96 nested tuples."
rlm@2 97 [a b]
rlm@2 98 (if (not (contractible? a b))
rlm@2 99 (throw
rlm@2 100 (Exception. "Not compatible for contraction."))
rlm@2 101 (reduce +
rlm@2 102 (map
rlm@2 103 (fn [x y]
rlm@2 104 (if (numbers? x y)
rlm@2 105 (* x y)
rlm@2 106 (contract x y)))
rlm@2 107 a b))))
rlm@2 108
rlm@2 109 ;***** MATRICES
rlm@2 110 (in-ns 'sicm.utils)
rlm@2 111 (require 'incanter.core) ;; use incanter's fast matrices
rlm@2 112
rlm@2 113 (defprotocol Matrix
rlm@2 114 (rows [this])
rlm@2 115 (cols [this])
rlm@2 116 (diagonal [this])
rlm@2 117 (trace [this])
rlm@2 118 (determinant [this]))
rlm@2 119
rlm@2 120 (extend-protocol Matrix
rlm@2 121 incanter.Matrix
rlm@2 122 (rows [this] (map down this)))
rlm@2 123
rlm@2 124
rlm@2 125
rlm@2 126
rlm@2 127 (defn count-rows [matrix]
rlm@2 128 ((comp count rows) matrix))
rlm@2 129
rlm@2 130 (defn count-cols [matrix]
rlm@2 131 ((comp count cols) matrix))
rlm@2 132
rlm@2 133
rlm@2 134 (defn matrix-by-rows
rlm@2 135 "Define a matrix by giving its rows."
rlm@2 136 [& rows]
rlm@2 137 (cond
rlm@2 138 (not (all-equal? (map count rows)))
rlm@2 139 (throw (Exception. "All rows in a matrix must have the same number of elements."))
rlm@2 140 :else
rlm@2 141 (reify Matrix
rlm@2 142 (rows [this] (map down rows))
rlm@2 143 (cols [this] (map up (apply map vector rows)))
rlm@2 144 (diagonal [this] (map-indexed (fn [i row] (nth row i) rows)))
rlm@2 145 (trace [this]
rlm@2 146 (if (not= (count-rows this) (count-cols this))
rlm@2 147 (throw (Exception.
rlm@2 148 "Cannot take the trace of a non-square matrix."))
rlm@2 149 (reduce + (diagonal this))))
rlm@2 150
rlm@2 151 (determinant [this]
rlm@2 152 (if (not= (count-rows this) (count-cols this))
rlm@2 153 (throw (Exception.
rlm@2 154 "Cannot take the determinant of a non-square matrix."))
rlm@2 155 (reduce * (diagonal this))))
rlm@2 156 )))
rlm@2 157
rlm@2 158
rlm@2 159 (defn matrix-by-cols
rlm@2 160 "Define a matrix by giving its columns."
rlm@2 161 [& cols]
rlm@2 162 (cond
rlm@2 163 (not (all-equal? (map count cols)))
rlm@2 164 (throw (Exception. "All columns in a matrix must have the same number of elements."))
rlm@2 165 :else
rlm@2 166 (reify Matrix
rlm@2 167 (cols [this] (map up cols))
rlm@2 168 (rows [this] (map down (apply map vector cols)))
rlm@2 169 (diagonal [this] (map-indexed (fn [i col] (nth col i) cols)))
rlm@2 170 (trace [this]
rlm@2 171 (if (not= (count-cols this) (count-rows this))
rlm@2 172 (throw (Exception.
rlm@2 173 "Cannot take the trace of a non-square matrix."))
rlm@2 174 (reduce + (diagonal this))))
rlm@2 175
rlm@2 176 (determinant [this]
rlm@2 177 (if (not= (count-cols this) (count-rows this))
rlm@2 178 (throw (Exception.
rlm@2 179 "Cannot take the determinant of a non-square matrix."))
rlm@2 180 (reduce * (map-indexed (fn [i col] (nth col i)) cols))))
rlm@2 181 )))
rlm@2 182
rlm@2 183 (extend-protocol Matrix Tuple
rlm@2 184 (rows [this] (if (down? this)
rlm@2 185 (list this)
rlm@2 186 (map (comp up vector) this)))
rlm@2 187
rlm@2 188 (cols [this] (if (up? this)
rlm@2 189 (list this)
rlm@2 190 (map (comp down vector) this))))
rlm@2 191
rlm@2 192 (defn matrix-multiply [A B]
rlm@2 193 (apply matrix-by-rows
rlm@2 194 (for [a (rows A)]
rlm@2 195 (for [b (cols B)]
rlm@2 196 (contract a b)))))