Mercurial > dylan

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