comparison sicm/utils.clj @ 2:b4de894a1e2e

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 (in-ns 'sicm.utils)
5
6 ;; Let some objects have spin
7
8 (defprotocol Spinning
9 (up? [this])
10 (down? [this]))
11
12 (defn spin
13 "Returns the spin of the Spinning s, either :up or :down"
14 [#^Spinning s]
15 (cond (up? s) :up (down? s) :down))
16
17
18 ;; DEFINITION: A tuple is a sequence with spin
19
20 (deftype Tuple
21 [spin coll]
22
23 clojure.lang.Seqable
24 (seq [this] (seq (.coll this)))
25
26 clojure.lang.Counted
27 (count [this] (count (.coll this)))
28
29 Spinning
30 (up? [this] (= ::up (.spin this)))
31 (down? [this] (= ::down (.spin this))))
32
33 (defmethod print-method Tuple
34 [o w]
35 (print-simple
36 (if (up? o)
37 (str "u" (.coll o))
38 (str "d" (vec(.coll o))))
39 w))
40
41 (def tuple? #(= (type %) Tuple))
42
43 ;; CONSTRUCTORS
44
45 (defn up
46 "Create a new up-tuple containing the contents of coll."
47 [coll]
48 (Tuple. ::up coll))
49
50 (defn down
51 "Create a new down-tuple containing the contents of coll."
52 [coll]
53 (Tuple. ::down coll))
54
55 (defn same-spin
56 "Creates a tuple which has the same spin as tuple and which contains
57 the contents of coll."
58 [tuple coll]
59 (if (up? tuple)
60 (up coll)
61 (down coll)))
62
63 (defn opposite-spin
64 "Create a tuple which has opposite spin to tuple and which contains
65 the contents of coll."
66 [tuple coll]
67 (if (up? tuple)
68 (down coll)
69 (up coll)))
70 (in-ns 'sicm.utils)
71 (require 'incanter.core) ;; use incanter's fast matrices
72
73
74 (defn all-equal? [coll]
75 (if (empty? (rest coll)) true
76 (and (= (first coll) (second coll))
77 (recur (rest coll)))))
78
79
80 (defprotocol Matrix
81 (rows [matrix])
82 (cols [matrix])
83 (diagonal [matrix])
84 (trace [matrix])
85 (determinant [matrix])
86 (transpose [matrix])
87 (conjugate [matrix])
88 )
89
90 (extend-protocol Matrix incanter.Matrix
91 (rows [rs] (map down (apply map vector (apply map vector rs))))
92 (cols [rs] (map up (apply map vector rs)))
93 (diagonal [matrix] (incanter.core/diag matrix) )
94 (determinant [matrix] (incanter.core/det matrix))
95 (trace [matrix] (incanter.core/trace matrix))
96 (transpose [matrix] (incanter.core/trans matrix)))
97
98 (defn count-rows [matrix]
99 ((comp count rows) matrix))
100
101 (defn count-cols [matrix]
102 ((comp count cols) matrix))
103
104 (defn square? [matrix]
105 (= (count-rows matrix) (count-cols matrix)))
106
107 (defn identity-matrix
108 "Define a square matrix of size n-by-n with 1s along the diagonal and
109 0s everywhere else."
110 [n]
111 (incanter.core/identity-matrix n))
112
113
114 (defn matrix-by-rows
115 "Define a matrix by giving its rows."
116 [& rows]
117 (if
118 (not (all-equal? (map count rows)))
119 (throw (Exception. "All rows in a matrix must have the same number of elements."))
120 (incanter.core/matrix (vec rows))))
121
122 (defn matrix-by-cols
123 "Define a matrix by giving its columns"
124 [& cols]
125 (if (not (all-equal? (map count cols)))
126 (throw (Exception. "All columns in a matrix must have the same number of elements."))
127 (incanter.core/matrix (vec (apply map vector cols)))))
128
129 (defn identity-matrix
130 "Define a square matrix of size n-by-n with 1s along the diagonal and
131 0s everywhere else."
132 [n]
133 (incanter.core/identity-matrix n))
134
135 (in-ns 'sicm.utils)
136 (use 'clojure.contrib.generic.arithmetic
137 'clojure.contrib.generic.collection
138 'clojure.contrib.generic.functor
139 'clojure.contrib.generic.math-functions)
140
141 (defn numbers?
142 "Returns true if all arguments are numbers, else false."
143 [& xs]
144 (every? number? xs))
145
146 (defn tuple-surgery
147 "Applies the function f to the items of tuple and the additional
148 arguments, if any. Returns a Tuple of the same type as tuple."
149 [tuple f & xs]
150 ((if (up? tuple) up down)
151 (apply f (seq tuple) xs)))
152
153
154
155 ;;; CONTRACTION collapses two compatible tuples into a number.
156
157 (defn contractible?
158 "Returns true if the tuples a and b are compatible for contraction,
159 else false. Tuples are compatible if they have the same number of
160 components, they have opposite spins, and their elements are
161 pairwise-compatible."
162 [a b]
163 (and
164 (isa? (type a) Tuple)
165 (isa? (type b) Tuple)
166 (= (count a) (count b))
167 (not= (spin a) (spin b))
168
169 (not-any? false?
170 (map #(or
171 (numbers? %1 %2)
172 (contractible? %1 %2))
173 a b))))
174
175 (defn contract
176 "Contracts two tuples, returning the sum of the
177 products of the corresponding items. Contraction is recursive on
178 nested tuples."
179 [a b]
180 (if (not (contractible? a b))
181 (throw
182 (Exception. "Not compatible for contraction."))
183 (reduce +
184 (map
185 (fn [x y]
186 (if (numbers? x y)
187 (* x y)
188 (contract x y)))
189 a b))))
190
191
192
193
194
195 (defmethod conj Tuple
196 [tuple & xs]
197 (tuple-surgery tuple #(apply conj % xs)))
198
199 (defmethod fmap Tuple
200 [f tuple]
201 (tuple-surgery tuple (partial map f)))
202
203
204
205 ;; TODO: define Scalar, and add it to the hierarchy above Number and Complex
206
207
208 (defmethod * [Tuple Tuple] ; tuple*tuple
209 [a b]
210 (if (contractible? a b)
211 (contract a b)
212 (map (partial * a) b)))
213
214
215 (defmethod * [java.lang.Number Tuple] ;; scalar * tuple
216 [a x] (fmap (partial * a) x))
217
218 (defmethod * [Tuple java.lang.Number]
219 [x a] (* a x))
220
221 (defmethod * [java.lang.Number incanter.Matrix] ;; scalar * matrix
222 [x M] (incanter.core/mult x M))
223
224 (defmethod * [incanter.Matrix java.lang.Number]
225 [M x] (* x M))
226
227 (defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
228 [M1 M2]
229 (incanter.core/mmult M1 M2))
230
231 (defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
232 [M v]
233 (if (and (apply numbers? v) (up? v))
234 (* M (matrix-by-cols v))
235 (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
236 ))
237
238 (defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
239 [v M]
240 (if (and (apply numbers? v) (down? v))
241 (* (matrix-by-rows v) M)
242 (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
243 ))
244
245
246 (defmethod exp incanter.Matrix
247 [M]
248 (incanter.core/exp M))
249
250
251 (in-ns 'sicm.utils)
252 (use 'clojure.contrib.seq
253 'clojure.contrib.generic.arithmetic
254 'clojure.contrib.generic.collection
255 'clojure.contrib.generic.math-functions)
256
257 ;;∂
258
259 ;; DEFINITION : Differential Term
260
261 ;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
262 ;; coefficient of the quantity is returned by the 'coefficient' method,
263 ;; while the sequence of differential parameters is returned by the
264 ;; method 'partials'.
265
266 ;; Instead of using (potentially ambiguous) letters to denote
267 ;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].
268
269 ;; The coefficient can be any arithmetic object; the
270 ;; partials must be a nonrepeating sorted sequence of nonnegative
271 ;; integers.
272
273 (deftype DifferentialTerm [coefficient partials])
274
275 (defn differential-term
276 "Make a differential term from a coefficient and list of partials."
277 [coefficient partials]
278 (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials))
279 (DifferentialTerm. coefficient (set partials))
280 (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))
281
282
283 ;; DEFINITION : Differential Sequence
284 ;; A differential sequence is a sequence of differential terms, all with different partials.
285 ;; Internally, it is a map from the partials of each term to their coefficients.
286
287 (deftype DifferentialSeq
288 [terms]
289 ;;clojure.lang.IPersistentMap
290 clojure.lang.Associative
291 (assoc [this key val]
292 (DifferentialSeq.
293 (cons (differential-term val key) terms)))
294 (cons [this x]
295 (DifferentialSeq. (cons x terms)))
296 (containsKey [this key]
297 (not(nil? (find-first #(= (.partials %) key) terms))))
298 (count [this] (count (.terms this)))
299 (empty [this] (DifferentialSeq. []))
300 (entryAt [this key]
301 ((juxt #(.partials %) #(.coefficient %))
302 (find-first #(= (.partials %) key) terms)))
303 (seq [this] (seq (.terms this))))
304
305 (def differential? #(= (type %) DifferentialSeq))
306
307 (defn zeroth-order?
308 "Returns true if the differential sequence has at most a constant term."
309 [dseq]
310 (and
311 (differential? dseq)
312 (every?
313 #(= #{} %)
314 (keys (.terms dseq)))))
315
316 (defmethod fmap DifferentialSeq
317 [f dseq]
318 (DifferentialSeq.
319 (fmap f (.terms dseq))))
320
321
322
323
324 ;; BUILDING DIFFERENTIAL OBJECTS
325
326 (defn differential-seq
327 "Define a differential sequence by specifying an alternating
328 sequence of coefficients and lists of partials."
329 ([coefficient partials]
330 (DifferentialSeq. {(set partials) coefficient}))
331 ([coefficient partials & cps]
332 (if (odd? (count cps))
333 (throw (Exception. "differential-seq requires an even number of terms."))
334 (DifferentialSeq.
335 (reduce
336 #(assoc %1 (set (second %2)) (first %2))
337 {(set partials) coefficient}
338 (partition 2 cps))))))
339
340
341
342 (defn big-part
343 "Returns the part of the differential sequence that is finite,
344 i.e. not infinitely small. If the sequence is zeroth-order, returns
345 the coefficient of the zeroth-order term instead. "
346 [dseq]
347 (if (zeroth-order? dseq) (get (.terms dseq) #{})
348 (let [m (.terms dseq)
349 keys (sort-by count (keys m))
350 smallest-var (last (last keys))]
351 (DifferentialSeq.
352 (reduce
353 #(assoc %1 (first %2) (second %2))
354 {}
355 (remove #((first %) smallest-var) m))))))
356
357
358 (defn small-part
359 "Returns the part of the differential sequence that infinitely
360 small. If the sequence is zeroth-order, returns zero."
361 [dseq]
362 (if (zeroth-order? dseq) 0
363 (let [m (.terms dseq)
364 keys (sort-by count (keys m))
365 smallest-var (last (last keys))]
366 (DifferentialSeq.
367 (reduce
368 #(assoc %1 (first %2) (second %2)) {}
369 (filter #((first %) smallest-var) m))))))
370
371
372
373 (defn cartesian-product [set1 set2]
374 (reduce concat
375 (for [x set1]
376 (for [y set2]
377 [x y]))))
378
379 (defn nth-subset [n]
380 (if (zero? n) []
381 (let [lg2 #(/ (log %) (log 2))
382 k (int(java.lang.Math/floor (lg2 n)))
383 ]
384 (cons k
385 (nth-subset (- n (pow 2 k)))))))
386
387 (def all-partials
388 (lazy-seq (map nth-subset (range))))
389
390
391 (defn differential-multiply
392 "Multiply two differential sequences. The square of any differential
393 variable is zero since differential variables are infinitesimally
394 small."
395 [dseq1 dseq2]
396 (DifferentialSeq.
397 (reduce
398 (fn [m [[vars1 coeff1] [vars2 coeff2]]]
399 (if (not (empty? (clojure.set/intersection vars1 vars2)))
400 m
401 (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
402 {}
403 (cartesian-product (.terms dseq1) (.terms dseq2)))))
404
405
406
407 (defmethod * [DifferentialSeq DifferentialSeq]
408 [dseq1 dseq2]
409 (differential-multiply dseq1 dseq2))
410
411 (defmethod + [DifferentialSeq DifferentialSeq]
412 [dseq1 dseq2]
413 (DifferentialSeq.
414 (merge-with + (.terms dseq1) (.terms dseq2))))
415
416 (defmethod * [java.lang.Number DifferentialSeq]
417 [x dseq]
418 (fmap (partial * x) dseq))
419
420 (defmethod * [DifferentialSeq java.lang.Number]
421 [dseq x]
422 (fmap (partial * x) dseq))
423
424 (defmethod + [java.lang.Number DifferentialSeq]
425 [x dseq]
426 (+ (differential-seq x []) dseq))
427 (defmethod + [DifferentialSeq java.lang.Number]
428 [dseq x]
429 (+ dseq (differential-seq x [])))
430
431 (defmethod - DifferentialSeq
432 [x]
433 (fmap - x))
434
435
436 ;; DERIVATIVES
437
438
439
440 (defn linear-approximator
441 "Returns an operator that linearly approximates the given function."
442 ([f df|dx]
443 (fn [x]
444 (let [big-part (big-part x)
445 small-part (small-part x)]
446 ;; f(x+dx) ~= f(x) + f'(x)dx
447 (+ (f big-part)
448 (* (df|dx big-part) small-part)
449 ))))
450
451 ([f df|dx df|dy]
452 (fn [x y]
453 (let [X (big-part x)
454 Y (big-part y)
455 DX (small-part x)
456 DY (small-part y)]
457 (+ (f X Y)
458 (* DX (f df|dx X Y))
459 (* DY (f df|dy X Y)))))))
460
461
462
463
464
465 (defn D[f]
466 (fn[x] (f (+ x (differential-seq 1 [0] 1 [1] 1 [2])))))
467
468 (defn d[partials f]
469 (fn [x]
470 (get
471 (.terms ((D f)x))
472 (set partials)
473 0
474 )))
475
476 (defmethod exp DifferentialSeq [x]
477 ((linear-approximator exp exp) x))
478
479 (defmethod sin DifferentialSeq
480 [x]
481 ((linear-approximator sin cos) x))
482
483 (defmethod cos DifferentialSeq
484 [x]
485 ((linear-approximator cos #(- (sin %))) x))
486
487 (defmethod log DifferentialSeq
488 [x]
489 ((linear-approximator log (fn [x] (/ x)) ) x))
490
491 (defmethod / [DifferentialSeq DifferentialSeq]
492 [x y]
493 ((linear-approximator /
494 (fn [x y] (/ 1 y))
495 (fn [x y] (- (/ x (* y y)))))
496 x y))