view sicm/utils.org @ 11:1f112b4f9e8f tip

Fixed what was baroque.
author Dylan Holmes <ocsenave@gmail.com>
date Tue, 01 Nov 2011 02:30:49 -0500
parents b4de894a1e2e
children
line wrap: on
line source
1 #+TITLE: Building a Classical Mechanics Library in Clojure
2 #+author: Dylan Holmes
3 #+EMAIL: rlm@mit.edu
4 #+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js"
5 #+STYLE: <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />
6 #+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
7 #+SETUPFILE: ../templates/level-0.org
8 #+INCLUDE: ../templates/level-0.org
9 #+BABEL: :noweb yes :results silent
11 * Useful data types
13 ** Complex numbers
15 ** Power series
17 ** Tuples and tensors
19 *** Tuples are \ldquo{}sequences with spin\rdquo{}
21 #+srcname: tuples
22 #+begin_src clojure
23 (in-ns 'sicm.utils)
25 ;; Let some objects have spin
27 (defprotocol Spinning
28 (up? [this])
29 (down? [this]))
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))
37 ;; DEFINITION: A tuple is a sequence with spin
39 (deftype Tuple
40 [spin coll]
42 clojure.lang.Seqable
43 (seq [this] (seq (.coll this)))
45 clojure.lang.Counted
46 (count [this] (count (.coll this)))
48 Spinning
49 (up? [this] (= ::up (.spin this)))
50 (down? [this] (= ::down (.spin this))))
52 (defmethod print-method Tuple
53 [o w]
54 (print-simple
55 (if (up? o)
56 (str "u" (.coll o))
57 (str "d" (vec(.coll o))))
58 w))
60 (def tuple? #(= (type %) Tuple))
62 ;; CONSTRUCTORS
64 (defn up
65 "Create a new up-tuple containing the contents of coll."
66 [coll]
67 (Tuple. ::up coll))
69 (defn down
70 "Create a new down-tuple containing the contents of coll."
71 [coll]
72 (Tuple. ::down coll))
74 (defn same-spin
75 "Creates a tuple which has the same spin as tuple and which contains
76 the contents of coll."
77 [tuple coll]
78 (if (up? tuple)
79 (up coll)
80 (down coll)))
82 (defn opposite-spin
83 "Create a tuple which has opposite spin to tuple and which contains
84 the contents of coll."
85 [tuple coll]
86 (if (up? tuple)
87 (down coll)
88 (up coll)))
89 #+end_src
93 *** Matrices
94 #+srcname:matrices
95 #+begin_src clojure
96 (in-ns 'sicm.utils)
97 (require 'incanter.core) ;; use incanter's fast matrices
100 (defn all-equal? [coll]
101 (if (empty? (rest coll)) true
102 (and (= (first coll) (second coll))
103 (recur (rest coll)))))
106 (defprotocol Matrix
107 (rows [matrix])
108 (cols [matrix])
109 (diagonal [matrix])
110 (trace [matrix])
111 (determinant [matrix])
112 (transpose [matrix])
113 (conjugate [matrix])
114 )
116 (extend-protocol Matrix incanter.Matrix
117 (rows [rs] (map down (apply map vector (apply map vector rs))))
118 (cols [rs] (map up (apply map vector rs)))
119 (diagonal [matrix] (incanter.core/diag matrix) )
120 (determinant [matrix] (incanter.core/det matrix))
121 (trace [matrix] (incanter.core/trace matrix))
122 (transpose [matrix] (incanter.core/trans matrix)))
124 (defn count-rows [matrix]
125 ((comp count rows) matrix))
127 (defn count-cols [matrix]
128 ((comp count cols) matrix))
130 (defn square? [matrix]
131 (= (count-rows matrix) (count-cols matrix)))
133 (defn identity-matrix
134 "Define a square matrix of size n-by-n with 1s along the diagonal and
135 0s everywhere else."
136 [n]
137 (incanter.core/identity-matrix n))
140 (defn matrix-by-rows
141 "Define a matrix by giving its rows."
142 [& rows]
143 (if
144 (not (all-equal? (map count rows)))
145 (throw (Exception. "All rows in a matrix must have the same number of elements."))
146 (incanter.core/matrix (vec rows))))
148 (defn matrix-by-cols
149 "Define a matrix by giving its columns"
150 [& cols]
151 (if (not (all-equal? (map count cols)))
152 (throw (Exception. "All columns in a matrix must have the same number of elements."))
153 (incanter.core/matrix (vec (apply map vector cols)))))
155 (defn identity-matrix
156 "Define a square matrix of size n-by-n with 1s along the diagonal and
157 0s everywhere else."
158 [n]
159 (incanter.core/identity-matrix n))
161 #+end_src
163 ** Generic Operations
164 #+srcname:arith-tuple
165 #+begin_src clojure
166 (in-ns 'sicm.utils)
167 (use 'clojure.contrib.generic.arithmetic
168 'clojure.contrib.generic.collection
169 'clojure.contrib.generic.functor
170 'clojure.contrib.generic.math-functions)
172 (defn numbers?
173 "Returns true if all arguments are numbers, else false."
174 [& xs]
175 (every? number? xs))
177 (defn tuple-surgery
178 "Applies the function f to the items of tuple and the additional
179 arguments, if any. Returns a Tuple of the same type as tuple."
180 [tuple f & xs]
181 ((if (up? tuple) up down)
182 (apply f (seq tuple) xs)))
186 ;;; CONTRACTION collapses two compatible tuples into a number.
188 (defn contractible?
189 "Returns true if the tuples a and b are compatible for contraction,
190 else false. Tuples are compatible if they have the same number of
191 components, they have opposite spins, and their elements are
192 pairwise-compatible."
193 [a b]
194 (and
195 (isa? (type a) Tuple)
196 (isa? (type b) Tuple)
197 (= (count a) (count b))
198 (not= (spin a) (spin b))
200 (not-any? false?
201 (map #(or
202 (numbers? %1 %2)
203 (contractible? %1 %2))
204 a b))))
206 (defn contract
207 "Contracts two tuples, returning the sum of the
208 products of the corresponding items. Contraction is recursive on
209 nested tuples."
210 [a b]
211 (if (not (contractible? a b))
212 (throw
213 (Exception. "Not compatible for contraction."))
214 (reduce +
215 (map
216 (fn [x y]
217 (if (numbers? x y)
218 (* x y)
219 (contract x y)))
220 a b))))
226 (defmethod conj Tuple
227 [tuple & xs]
228 (tuple-surgery tuple #(apply conj % xs)))
230 (defmethod fmap Tuple
231 [f tuple]
232 (tuple-surgery tuple (partial map f)))
236 ;; TODO: define Scalar, and add it to the hierarchy above Number and Complex
239 (defmethod * [Tuple Tuple] ; tuple*tuple
240 [a b]
241 (if (contractible? a b)
242 (contract a b)
243 (map (partial * a) b)))
246 (defmethod * [java.lang.Number Tuple] ;; scalar * tuple
247 [a x] (fmap (partial * a) x))
249 (defmethod * [Tuple java.lang.Number]
250 [x a] (* a x))
252 (defmethod * [java.lang.Number incanter.Matrix] ;; scalar * matrix
253 [x M] (incanter.core/mult x M))
255 (defmethod * [incanter.Matrix java.lang.Number]
256 [M x] (* x M))
258 (defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
259 [M1 M2]
260 (incanter.core/mmult M1 M2))
262 (defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
263 [M v]
264 (if (and (apply numbers? v) (up? v))
265 (* M (matrix-by-cols v))
266 (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
267 ))
269 (defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
270 [v M]
271 (if (and (apply numbers? v) (down? v))
272 (* (matrix-by-rows v) M)
273 (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
274 ))
277 (defmethod exp incanter.Matrix
278 [M]
279 (incanter.core/exp M))
281 #+end_src
283 * Operators and Differentiation
284 ** Operators
285 #+scrname: operators
286 #+begin_src clojure
287 (in-ns 'sicm.utils)
288 (use 'clojure.contrib.seq
289 'clojure.contrib.generic.arithmetic
290 'clojure.contrib.generic.collection
291 'clojure.contrib.generic.math-functions)
293 (defmethod + [clojure.lang.IFn clojure.lang.IFn]
294 [f g]
295 (fn [& args]
296 (+ (apply f args) (apply g args))))
298 (defmethod * [clojure.lang.IFn clojure.lang.IFn]
299 [f g]
300 (fn [& args]
301 (* (apply f args) (apply g args))))
303 (defmethod / [clojure.lang.IFn java.lang.Number]
304 [f x]
305 (fn [& args]
306 (/ (apply f args) x)))
309 (defmethod - [clojure.lang.IFn]
310 [f]
311 (fn [& args]
312 (- (apply f args))))
314 (defmethod - [clojure.lang.IFn clojure.lang.IFn]
315 [f g]
316 (fn [& args]
317 (- (apply f args) (apply g args))))
319 (defmethod pow [clojure.lang.IFn java.lang.Number]
320 [f x]
321 (fn [& args]
322 (pow (apply f args) x)))
325 (defmethod + [java.lang.Number clojure.lang.IFn]
326 [x f]
327 (fn [& args]
328 (+ x (apply f args))))
330 (defmethod * [java.lang.Number clojure.lang.IFn]
331 [x f]
332 (fn [& args]
333 (* x (apply f args))))
335 (defmethod * [clojure.lang.IFn java.lang.Number]
336 [f x]
337 (* x f))
338 (defmethod + [clojure.lang.IFn java.lang.Number]
339 [f x]
340 (+ x f))
342 #+end_src
344 ** Differential Terms and Sequences
345 #+srcname: differential
346 #+begin_src clojure
347 (in-ns 'sicm.utils)
348 (use 'clojure.contrib.seq
349 'clojure.contrib.generic.arithmetic
350 'clojure.contrib.generic.collection
351 'clojure.contrib.generic.math-functions)
353 ;;∂
355 ;; DEFINITION : Differential Term
357 ;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
358 ;; coefficient of the quantity is returned by the 'coefficient' method,
359 ;; while the sequence of differential parameters is returned by the
360 ;; method 'partials'.
362 ;; Instead of using (potentially ambiguous) letters to denote
363 ;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].
365 ;; The coefficient can be any arithmetic object; the
366 ;; partials must be a nonrepeating sorted sequence of nonnegative
367 ;; integers.
369 ;; (deftype DifferentialTerm [coefficient partials])
371 ;; (defn differential-term
372 ;; "Make a differential term from a coefficient and list of partials."
373 ;; [coefficient partials]
374 ;; (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials))
375 ;; (DifferentialTerm. coefficient (set partials))
376 ;; (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))
379 ;; DEFINITION : Differential Sequence
380 ;; A differential sequence is a sequence of differential terms, all with different partials.
381 ;; Internally, it is a map from the partials of each term to their coefficients.
383 (deftype DifferentialSeq
384 [terms]
385 ;;clojure.lang.IPersistentMap
386 clojure.lang.Associative
387 (assoc [this key val]
388 (DifferentialSeq.
389 (cons (differential-term val key) terms)))
390 (cons [this x]
391 (DifferentialSeq. (cons x terms)))
392 (containsKey [this key]
393 (not(nil? (find-first #(= (.partials %) key) terms))))
394 (count [this] (count (.terms this)))
395 (empty [this] (DifferentialSeq. []))
396 (entryAt [this key]
397 ((juxt #(.partials %) #(.coefficient %))
398 (find-first #(= (.partials %) key) terms)))
399 (seq [this] (seq (.terms this))))
401 (def differential? #(= (type %) DifferentialSeq))
403 (defn zeroth-order?
404 "Returns true if the differential sequence has at most a constant term."
405 [dseq]
406 (and
407 (differential? dseq)
408 (every?
409 #(= #{} %)
410 (keys (.terms dseq)))))
412 (defmethod fmap DifferentialSeq
413 [f dseq]
414 (DifferentialSeq.
415 (fmap f (.terms dseq))))
420 ;; BUILDING DIFFERENTIAL OBJECTS
422 (defn differential-seq
423 "Define a differential sequence by specifying an alternating
424 sequence of coefficients and lists of partials."
425 ([coefficient partials]
426 (DifferentialSeq. {(set partials) coefficient}))
427 ([coefficient partials & cps]
428 (if (odd? (count cps))
429 (throw (Exception. "differential-seq requires an even number of terms."))
430 (DifferentialSeq.
431 (reduce
432 #(assoc %1 (set (second %2)) (first %2))
433 {(set partials) coefficient}
434 (partition 2 cps))))))
438 (defn big-part
439 "Returns the part of the differential sequence that is finite,
440 i.e. not infinitely small. If the sequence is zeroth-order, returns
441 the coefficient of the zeroth-order term instead. "
442 [dseq]
443 (if (zeroth-order? dseq) (get (.terms dseq) #{})
444 (let [m (.terms dseq)
445 keys (sort-by count (keys m))
446 smallest-var (last (last keys))]
447 (DifferentialSeq.
448 (reduce
449 #(assoc %1 (first %2) (second %2))
450 {}
451 (remove #((first %) smallest-var) m))))))
454 (defn small-part
455 "Returns the part of the differential sequence that infinitely
456 small. If the sequence is zeroth-order, returns zero."
457 [dseq]
458 (if (zeroth-order? dseq) 0
459 (let [m (.terms dseq)
460 keys (sort-by count (keys m))
461 smallest-var (last (last keys))]
462 (DifferentialSeq.
463 (reduce
464 #(assoc %1 (first %2) (second %2)) {}
465 (filter #((first %) smallest-var) m))))))
469 (defn cartesian-product [set1 set2]
470 (reduce concat
471 (for [x set1]
472 (for [y set2]
473 [x y]))))
475 (defn nth-subset [n]
476 (if (zero? n) []
477 (let [lg2 #(/ (log %) (log 2))
478 k (int(java.lang.Math/floor (lg2 n)))
479 ]
480 (cons k
481 (nth-subset (- n (pow 2 k)))))))
483 (def all-partials
484 (lazy-seq (map nth-subset (range))))
487 (defn differential-multiply
488 "Multiply two differential sequences. The square of any differential
489 variable is zero since differential variables are infinitesimally
490 small."
491 [dseq1 dseq2]
492 (DifferentialSeq.
493 (reduce
494 (fn [m [[vars1 coeff1] [vars2 coeff2]]]
495 (if (not (empty? (clojure.set/intersection vars1 vars2)))
496 m
497 (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
498 {}
499 (cartesian-product (.terms dseq1) (.terms dseq2)))))
503 (defmethod * [DifferentialSeq DifferentialSeq]
504 [dseq1 dseq2]
505 (differential-multiply dseq1 dseq2))
507 (defmethod + [DifferentialSeq DifferentialSeq]
508 [dseq1 dseq2]
509 (DifferentialSeq.
510 (merge-with + (.terms dseq1) (.terms dseq2))))
512 (defmethod * [java.lang.Number DifferentialSeq]
513 [x dseq]
514 (fmap (partial * x) dseq))
516 (defmethod * [DifferentialSeq java.lang.Number]
517 [dseq x]
518 (fmap (partial * x) dseq))
520 (defmethod + [java.lang.Number DifferentialSeq]
521 [x dseq]
522 (+ (differential-seq x []) dseq))
523 (defmethod + [DifferentialSeq java.lang.Number]
524 [dseq x]
525 (+ dseq (differential-seq x [])))
527 (defmethod - DifferentialSeq
528 [x]
529 (fmap - x))
532 ;; DERIVATIVES
536 (defn linear-approximator
537 "Returns an operator that linearly approximates the given function."
538 ([f df|dx]
539 (fn [x]
540 (let [big-part (big-part x)
541 small-part (small-part x)]
542 ;; f(x+dx) ~= f(x) + f'(x)dx
543 (+ (f big-part)
544 (* (df|dx big-part) small-part)
545 ))))
547 ([f df|dx df|dy]
548 (fn [x y]
549 (let [X (big-part x)
550 Y (big-part y)
551 DX (small-part x)
552 DY (small-part y)]
553 (+ (f X Y)
554 (* DX (f df|dx X Y))
555 (* DY (f df|dy X Y)))))))
561 (defn D[f]
562 (fn[x] (f (+ x (differential-seq 1 [0] 1 [1] 1 [2])))))
564 (defn d[partials f]
565 (fn [x]
566 (get
567 (.terms ((D f)x))
568 (set partials)
569 0
570 )))
572 (defmethod exp DifferentialSeq [x]
573 ((linear-approximator exp exp) x))
575 (defmethod sin DifferentialSeq
576 [x]
577 ((linear-approximator sin cos) x))
579 (defmethod cos DifferentialSeq
580 [x]
581 ((linear-approximator cos #(- (sin %))) x))
583 (defmethod log DifferentialSeq
584 [x]
585 ((linear-approximator log (fn [x] (/ x)) ) x))
587 (defmethod / [DifferentialSeq DifferentialSeq]
588 [x y]
589 ((linear-approximator /
590 (fn [x y] (/ 1 y))
591 (fn [x y] (- (/ x (* y y)))))
592 x y))
594 #+end_src
598 * Derivatives Revisited
599 #+begin_src clojure
600 (in-ns 'sicm.utils)
601 (use 'clojure.contrib.seq
602 'clojure.contrib.generic.arithmetic
603 'clojure.contrib.generic.collection
604 'clojure.contrib.generic.math-functions)
606 (defn replace-in
607 "Replaces the nth item in coll with the given item. If n is larger
608 than the size of coll, adds n to the end of the collection."
609 [coll n item]
610 (concat
611 (take n coll)
612 [item]
613 (drop (inc n) coll)))
615 (defn euclidean-structure [f partials]
616 (fn sd[g v]
617 (cond
618 (tuple? v)
619 (opposite-spin
620 v
621 (map
622 (fn [n]
623 (sd (fn [xn]
624 (g
625 (same-spin v
626 (replace-in v n xn))))
627 (nth v n)))
628 (range (count v)))))))
633 #+end_src
636 * Symbolic Quantities
638 #+srcname: symbolic
639 #+begin_src clojure
640 (in-ns 'sicm.utils)
641 (use 'clojure.contrib.generic.arithmetic
642 'clojure.contrib.generic.math-functions)
644 (deftype Symbolic
645 [type
646 expression])
648 (defn print-expression [s]
649 (print (.expression s)))
651 (defn symbolic-number
652 [symbol]
653 (Symbolic. java.lang.Number symbol))
655 (defn simplify-expression [s]
656 (let [e (.expression s)]
657 (cond
658 (not(list? e)) e
659 (= (first e) '+)
660 )
664 (defmethod + [Symbolic Symbolic]
665 [x y]
666 (Symbolic. (.type x) (list '+ (.expression x) (.expression y))))
668 (defmethod + [Symbolic java.lang.Number]
669 [s x]
670 (if (zero? x)
671 s
672 (Symbolic. (.type s) (list '+ (.expression s) x))))
674 (defmethod sin Symbolic
675 [s]
676 (Symbolic. (.type s) (list 'sin (.expression s))))
678 #+end_src
680 * COMMENT To-be-tangled Source
681 #+begin_src clojure :tangle utils.clj
682 (ns sicm.utils)
684 <<tuples>>
685 <<matrices>>
686 <<arith-tuple>>
688 <<differential>>
689 #+end_src