view sicm/bk/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: Robert McIntyre & 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 * Generic Arithmetic
13 #+srcname: generic-arithmetic
14 #+begin_src clojure
15 (ns sicm.utils)
16 (in-ns 'sicm.utils)
20 (defn all-equal? [coll]
21 (if (empty? (rest coll)) true
22 (and (= (first coll) (second coll))
23 (recur (rest coll)))))
26 (defprotocol Arithmetic
27 (zero [this])
28 (one [this])
29 (negate [this])
30 (invert [this])
31 (conjugate [this]))
33 (defn zero? [x]
34 (= x (zero x)))
35 (defn one? [x]
36 (= x (one x)))
40 (extend-protocol Arithmetic
41 java.lang.Number
42 (zero [x] 0)
43 (one [x] 1)
44 (negate [x] (- x))
45 (invert [x] (/ x))
46 )
48 (extend-protocol Arithmetic
49 clojure.lang.IFn
50 (one [f] identity)
51 (negate [f] (comp negate f))
52 (invert [f] (comp invert f)))
55 (extend-protocol Arithmetic
56 clojure.lang.Seqable
57 (zero [this] (map zero this))
58 (one [this] (map one this))
59 (invert [this] (map invert this))
60 (negate [this] (map negate this)))
63 (defn ordered-like
64 "Create a comparator using the sorted collection as an
65 example. Elements not in the sorted collection are sorted to the
66 end."
67 [sorted-coll]
68 (let [
69 sorted-coll? (set sorted-coll)
70 ascending-pair?
71 (set(reduce concat
72 (map-indexed
73 (fn [n x]
74 (map #(vector x %) (nthnext sorted-coll n)))
75 sorted-coll)))]
76 (fn [x y]
77 (cond
78 (= x y) 0
79 (ascending-pair? [x y]) -1
80 (ascending-pair? [y x]) 1
81 (sorted-coll? x) -1
82 (sorted-coll? y) 1))))
86 (def type-precedence
87 (ordered-like [incanter.Matrix]))
89 (defmulti add
90 (fn [x y]
91 (sort type-precedence [(type x)(type y)])))
93 (defmulti multiply
94 (fn [x y]
95 (sort type-precedence [(type x) (type y)])))
97 (defmethod add [java.lang.Number java.lang.Number] [x y] (+ x y))
98 (defmethod multiply [java.lang.Number java.lang.Number] [x y] (* x y))
100 (defmethod multiply [incanter.Matrix java.lang.Integer] [x y]
101 (let [args (sort #(type-precedence (type %1)(type %2)) [x y])
102 matrix (first args)
103 scalar (second args)]
104 (incanter.core/matrix (map (partial map (partial multiply scalar)) matrix))))
106 #+end_src
109 * Useful Data Types
111 ** Complex Numbers
112 #+srcname: complex-numbers
113 #+begin_src clojure
114 (in-ns 'sicm.utils)
116 (defprotocol Complex
117 (real-part [z])
118 (imaginary-part [z])
119 (magnitude-squared [z])
120 (angle [z])
121 (conjugate [z])
122 (norm [z]))
124 (defn complex-rectangular
125 "Define a complex number with the given real and imaginary
126 components."
127 [re im]
128 (reify Complex
129 (real-part [z] re)
130 (imaginary-part [z] im)
131 (magnitude-squared [z] (+ (* re re) (* im im)))
132 (angle [z] (java.lang.Math/atan2 im re))
133 (conjugate [z] (complex-rectangular re (- im)))
135 Arithmetic
136 (zero [z] (complex-rectangular 0 0))
137 (one [z] (complex-rectangular 1 0))
138 (negate [z] (complex-rectangular (- re) (- im)))
139 (invert [z] (complex-rectangular
140 (/ re (magnitude-squared z))
141 (/ (- im) (magnitude-squared z))))
143 Object
144 (toString [_]
145 (if (and (zero? re) (zero? im)) (str 0)
146 (str
147 (if (not(zero? re))
148 re)
149 (if ((comp not zero?) im)
150 (str
151 (if (neg? im) "-" "+")
152 (if ((comp not one?) (java.lang.Math/abs im))
153 (java.lang.Math/abs im))
154 "i")))))))
156 (defn complex-polar
157 "Define a complex number with the given magnitude and angle."
158 [mag ang]
159 (reify Complex
160 (magnitude-squared [z] (* mag mag))
161 (angle [z] angle)
162 (real-part [z] (* mag (java.lang.Math/cos ang)))
163 (imaginary-part [z] (* mag (java.lang.Math/sin ang)))
164 (conjugate [z] (complex-polar mag (- ang)))
166 Arithmetic
167 (zero [z] (complex-polar 0 0))
168 (one [z] (complex-polar 1 0))
169 (negate [z] (complex-polar (- mag) ang))
170 (invert [z] (complex-polar (/ mag) (- ang)))
172 Object
173 (toString [_] (str mag " * e^(i" ang")"))
174 ))
177 ;; Numbers are complex quantities
179 (extend-protocol Complex
180 java.lang.Number
181 (real-part [x] x)
182 (imaginary-part [x] 0)
183 (magnitude [x] x)
184 (angle [x] 0)
185 (conjugate [x] x))
188 #+end_src
192 ** Tuples and Tensors
194 A tuple is a vector which is spinable\mdash{}it can be either /spin
195 up/ or /spin down/. (Covariant, contravariant; dual vectors)
196 #+srcname: tuples
197 #+begin_src clojure
198 (in-ns 'sicm.utils)
200 (defprotocol Spinning
201 (up? [this])
202 (down? [this]))
204 (defn spin
205 "Returns the spin of the Spinning s, either :up or :down"
206 [#^Spinning s]
207 (cond (up? s) :up (down? s) :down))
210 (deftype Tuple
211 [spin coll]
212 clojure.lang.Seqable
213 (seq [this] (seq (.coll this)))
214 clojure.lang.Counted
215 (count [this] (count (.coll this))))
217 (extend-type Tuple
218 Spinning
219 (up? [this] (= ::up (.spin this)))
220 (down? [this] (= ::down (.spin this))))
222 (defmethod print-method Tuple
223 [o w]
224 (print-simple (str (if (up? o) 'u 'd) (.coll o)) w))
228 (defn up
229 "Create a new up-tuple containing the contents of coll."
230 [coll]
231 (Tuple. ::up coll))
233 (defn down
234 "Create a new down-tuple containing the contents of coll."
235 [coll]
236 (Tuple. ::down coll))
239 #+end_src
241 *** Contraction
242 Contraction is a binary operation that you can apply to compatible
243 tuples. Tuples are compatible for contraction if they have the same
244 length and opposite spins, and if the corresponding items in each
245 tuple are both numbers or both compatible tuples.
247 #+srcname: tuples-2
248 #+begin_src clojure
249 (in-ns 'sicm.utils)
251 (defn numbers?
252 "Returns true if all arguments are numbers, else false."
253 [& xs]
254 (every? number? xs))
256 (defn contractible?
257 "Returns true if the tuples a and b are compatible for contraction,
258 else false. Tuples are compatible if they have the same number of
259 components, they have opposite spins, and their elements are
260 pairwise-compatible."
261 [a b]
262 (and
263 (isa? (type a) Tuple)
264 (isa? (type b) Tuple)
265 (= (count a) (count b))
266 (not= (spin a) (spin b))
268 (not-any? false?
269 (map #(or
270 (numbers? %1 %2)
271 (contractible? %1 %2))
272 a b))))
276 (defn contract
277 "Contracts two tuples, returning the sum of the
278 products of the corresponding items. Contraction is recursive on
279 nested tuples."
280 [a b]
281 (if (not (contractible? a b))
282 (throw
283 (Exception. "Not compatible for contraction."))
284 (reduce +
285 (map
286 (fn [x y]
287 (if (numbers? x y)
288 (* x y)
289 (contract x y)))
290 a b))))
292 #+end_src
294 *** Matrices
295 #+srcname: matrices
296 #+begin_src clojure
297 (in-ns 'sicm.utils)
298 (require 'incanter.core) ;; use incanter's fast matrices
300 (defprotocol Matrix
301 (rows [matrix])
302 (cols [matrix])
303 (diagonal [matrix])
304 (trace [matrix])
305 (determinant [matrix])
306 (transpose [matrix])
307 (conjugate [matrix])
308 )
310 (extend-protocol Matrix
311 incanter.Matrix
312 (rows [rs] (map down (apply map vector (apply map vector rs))))
313 (cols [rs] (map up (apply map vector rs)))
314 (diagonal [matrix] (incanter.core/diag matrix) )
315 (determinant [matrix] (incanter.core/det matrix))
316 (trace [matrix] (incanter.core/trace matrix))
317 (transpose [matrix] (incanter.core/trans matrix))
318 )
320 (defn count-rows [matrix]
321 ((comp count rows) matrix))
323 (defn count-cols [matrix]
324 ((comp count cols) matrix))
326 (defn square? [matrix]
327 (= (count-rows matrix) (count-cols matrix)))
329 (defn identity-matrix
330 "Define a square matrix of size n-by-n with 1s along the diagonal and
331 0s everywhere else."
332 [n]
333 (incanter.core/identity-matrix n))
339 (defn matrix-by-rows
340 "Define a matrix by giving its rows."
341 [& rows]
342 (if
343 (not (all-equal? (map count rows)))
344 (throw (Exception. "All rows in a matrix must have the same number of elements."))
345 (incanter.core/matrix (vec rows))))
347 (defn matrix-by-cols
348 "Define a matrix by giving its columns"
349 [& cols]
350 (if (not (all-equal? (map count cols)))
351 (throw (Exception. "All columns in a matrix must have the same number of elements."))
352 (incanter.core/matrix (vec (apply map vector cols)))))
354 (defn identity-matrix
355 "Define a square matrix of size n-by-n with 1s along the diagonal and
356 0s everywhere else."
357 [n]
358 (incanter.core/identity-matrix n))
362 (extend-protocol Arithmetic
363 incanter.Matrix
364 (one [matrix]
365 (if (square? matrix)
366 (identity-matrix (count-rows matrix))
367 (throw (Exception. "Non-square matrices have no multiplicative unit."))))
368 (zero [matrix]
369 (apply matrix-by-rows (map zero (rows matrix))))
370 (negate [matrix]
371 (apply matrix-by-rows (map negate (rows matrix))))
372 (invert [matrix]
373 (incanter.core/solve matrix)))
377 (defmulti coerce-to-matrix
378 "Converts x into a matrix, if possible."
379 type)
381 (defmethod coerce-to-matrix incanter.Matrix [x] x)
382 (defmethod coerce-to-matrix Tuple [x]
383 (if (apply numbers? (seq x))
384 (if (up? x)
385 (matrix-by-cols (seq x))
386 (matrix-by-rows (seq x)))
387 (throw (Exception. "Non-numerical tuple cannot be converted into a matrix."))))
393 ;; (defn matrix-by-cols
394 ;; "Define a matrix by giving its columns."
395 ;; [& cols]
396 ;; (cond
397 ;; (not (all-equal? (map count cols)))
398 ;; (throw (Exception. "All columns in a matrix must have the same number of elements."))
399 ;; :else
400 ;; (reify Matrix
401 ;; (cols [this] (map up cols))
402 ;; (rows [this] (map down (apply map vector cols)))
403 ;; (diagonal [this] (map-indexed (fn [i col] (nth col i) cols)))
404 ;; (trace [this]
405 ;; (if (not= (count-cols this) (count-rows this))
406 ;; (throw (Exception.
407 ;; "Cannot take the trace of a non-square matrix."))
408 ;; (reduce + (diagonal this))))
410 ;; (determinant [this]
411 ;; (if (not= (count-cols this) (count-rows this))
412 ;; (throw (Exception.
413 ;; "Cannot take the determinant of a non-square matrix."))
414 ;; (reduce * (map-indexed (fn [i col] (nth col i)) cols))))
415 ;; )))
417 (extend-protocol Matrix Tuple
418 (rows [this] (if (down? this)
419 (list this)
420 (map (comp up vector) this)))
422 (cols [this] (if (up? this)
423 (list this)
424 (map (comp down vector) this))
425 ))
427 (defn matrix-multiply
428 "Returns the matrix resulting from the matrix multiplication of the given arguments."
429 ([A] (coerce-to-matrix A))
430 ([A B] (incanter.core/mmult (coerce-to-matrix A) (coerce-to-matrix B)))
431 ([M1 M2 & Ms] (reduce matrix-multiply (matrix-multiply M1 M2) Ms)))
433 #+end_src
436 ** Power Series
437 #+srcname power-series
438 #+begin_src clojure
439 (in-ns 'sicm.utils)
440 (use 'clojure.contrib.def)
444 (defn series-fn
445 "The function corresponding to the given power series."
446 [series]
447 (fn [x]
448 (reduce +
449 (map-indexed (fn[n x] (* (float (nth series n)) (float(java.lang.Math/pow (float x) n)) ))
450 (range 20)))))
452 (deftype PowerSeries
453 [coll]
454 clojure.lang.Seqable
455 (seq [this] (seq (.coll this)))
457 clojure.lang.Indexed
458 (nth [this n] (nth (.coll this) n 0))
459 (nth [this n not-found] (nth (.coll this) n not-found))
461 ;; clojure.lang.IFn
462 ;; (call [this] (throw(Exception.)))
463 ;; (invoke [this & args] args
464 ;; (let [f
465 ;; )
466 ;; (run [this] (throw(Exception.)))
467 )
469 (defn power-series
470 "Returns a power series with the items of the coll as its
471 coefficients. Trailing zeros are added to the end of coll."
472 [coeffs]
473 (PowerSeries. coeffs))
475 (defn power-series-indexed
476 "Returns a power series consisting of the result of mapping f to the non-negative integers."
477 [f]
478 (PowerSeries. (map f (range))))
481 (defn-memo nth-partial-sum
482 ([series n]
483 (if (zero? n) (first series)
484 (+ (nth series n)
485 (nth-partial-sum series (dec n))))))
487 (defn partial-sums [series]
488 (lazy-seq (map nth-partial-sum (range))))
493 (def cos-series
494 (power-series-indexed
495 (fn[n]
496 (if (odd? n) 0
497 (/
498 (reduce *
499 (reduce * (repeat (/ n 2) -1))
500 (range 1 (inc n)))
501 )))))
503 (def sin-series
504 (power-series-indexed
505 (fn[n]
506 (if (even? n) 0
507 (/
508 (reduce *
509 (reduce * (repeat (/ (dec n) 2) -1))
510 (range 1 (inc n)))
511 )))))
513 #+end_src
516 * Basic Utilities
518 ** Sequence manipulation
520 #+srcname: seq-manipulation
521 #+begin_src clojure
522 (ns sicm.utils)
524 (defn do-up
525 "Apply f to each number from low to high, presumably for
526 side-effects."
527 [f low high]
528 (doseq [i (range low high)] (f i)))
530 (defn do-down
531 "Apply f to each number from high to low, presumably for
532 side-effects."
533 [f high low]
534 (doseq [i (range high low -1)] (f i)))
537 (defn all-equal? [coll]
538 (if (empty? (rest coll)) true
539 (and (= (first coll) (second coll))
540 (recur (rest coll))))))
542 (defn multiplier
543 "Returns a function that 'multiplies' the members of a collection,
544 returning unit if called on an empty collection."
545 [multiply unit]
546 (fn [coll] ((partial reduce multiply unit) coll)))
548 (defn divider
549 "Returns a function that 'divides' the first element of a collection
550 by the 'product' of the rest of the collection."
551 [divide multiply invert unit]
552 (fn [coll]
553 (apply
554 (fn
555 ([] unit)
556 ([x] (invert x))
557 ([x y] (divide x y))
558 ([x y & zs] (divide x (reduce multiply y zs))))
559 coll)))
561 (defn left-circular-shift
562 "Remove the first element of coll, adding it to the end of coll."
563 [coll]
564 (concat (rest coll) (take 1 coll)))
566 (defn right-circular-shift
567 "Remove the last element of coll, adding it to the front of coll."
568 [coll]
569 (cons (last coll) (butlast coll)))
570 #+end_src
575 ** Ranges, Airity and Function Composition
576 #+srcname: arity
577 #+begin_src clojure
578 (in-ns 'sicm.utils)
579 (def infinity Double/POSITIVE_INFINITY)
580 (defn infinite? [x] (Double/isInfinite x))
581 (def finite? (comp not infinite?))
583 (defn arity-min
584 "Returns the smallest number of arguments f can take."
585 [f]
586 (apply
587 min
588 (map (comp alength #(.getParameterTypes %))
589 (filter (comp (partial = "invoke") #(.getName %))
590 (.getDeclaredMethods (class f))))))
592 (defn arity-max
593 "Returns the largest number of arguments f can take, possibly
594 Infinity."
595 [f]
596 (let [methods (.getDeclaredMethods (class f))]
597 (if (not-any? (partial = "doInvoke") (map #(.getName %) methods))
598 (apply max
599 (map (comp alength #(.getParameterTypes %))
600 (filter (comp (partial = "invoke") #(.getName %)) methods)))
601 infinity)))
604 (def ^{:arglists '([f])
605 :doc "Returns a two-element list containing the minimum and
606 maximum number of args that f can take."}
607 arity-interval
608 (juxt arity-min arity-max))
612 ;; --- intervals
614 (defn intersect
615 "Returns the interval of overlap between interval-1 and interval-2"
616 [interval-1 interval-2]
617 (if (or (empty? interval-1) (empty? interval-2)) []
618 (let [left (max (first interval-1) (first interval-2))
619 right (min (second interval-1) (second interval-2))]
620 (if (> left right) []
621 [left right]))))
623 (defn same-endpoints?
624 "Returns true if the left endpoint is the same as the right
625 endpoint."
626 [interval]
627 (= (first interval) (second interval)))
629 (defn naturals?
630 "Returns true if the left endpoint is 0 and the right endpoint is
631 infinite."
632 [interval]
633 (and (zero? (first interval))
634 (infinite? (second interval))))
637 (defn fan-in
638 "Returns a function that pipes its input to each of the gs, then
639 applies f to the list of results. Consequently, f must be able to
640 take a number of arguments equal to the number of gs."
641 [f & gs]
642 (fn [& args]
643 (apply f (apply (apply juxt gs) args))))
645 (defn fan-in
646 "Returns a function that pipes its input to each of the gs, then applies f to the list of results. The resulting function takes any number of arguments, but will fail if given arguments that are incompatible with any of the gs."
647 [f & gs]
648 (comp (partial apply f) (apply juxt gs)))
652 (defmacro airty-blah-sad [f n more?]
653 (let [syms (vec (map (comp gensym (partial str "x")) (range n)))
654 optional (gensym "xs")]
655 (if more?
656 `(fn ~(conj syms '& optional)
657 (apply ~f ~@syms ~optional))
658 `(fn ~syms (~f ~@syms)))))
660 (defmacro airt-whaa* [f n more?]
661 `(airty-blah-sad ~f ~n ~more?))
666 (defn fan-in*
667 "Returns a function that pipes its input to each of the gs, then
668 applies f to the list of results. Unlike fan-in, fan-in* strictly
669 enforces arity: it will fail if the gs do not have compatible
670 arities."
671 [f & gs]
672 (let [arity-in (reduce intersect (map arity-interval gs))
673 left (first arity-in)
674 right (second arity-in)
675 composite (fan-in f gs)
676 ]
677 (cond
678 (empty? arity-in)
679 (throw (Exception. "Cannot compose functions with incompatible arities."))
681 (not
682 (or (= left right)
683 (and (finite? left)
684 (= right infinity))))
686 (throw (Exception.
687 "Compose can only handle arities of the form [n n] or [n infinity]"))
688 :else
689 (airty-blah-sad composite left (= right infinity)))))
693 (defn compose-n "Compose any number of functions together."
694 ([] identity)
695 ([f] f)
696 ([f & fs]
697 (let [fns (cons f fs)]
698 (compose-bin (reduce fan-in (butlast fs)) (last fs))))
699 )
706 (defn iterated
707 ([f n id] (reduce comp id (repeat n f)))
708 ([f n] (reduce comp identity (repeat n f))))
710 (defn iterate-until-stable
711 "Repeatedly applies f to x, returning the first result that is close
712 enough to its predecessor."
713 [f close-enough? x]
714 (second (swank.util/find-first
715 (partial apply close-enough?)
716 (partition 2 1 (iterate f x)))))
718 (defn lexical< [x y]
719 (neg? (compare (str x) (str y))))
722 ;; do-up
723 ;; do-down
724 (def make-pairwise-test comparator)
725 ;;all-equal?
726 (def accumulation multiplier)
727 (def inverse-accumulation divider)
728 ;;left-circular-shift
729 ;;right-circular-shift
730 (def exactly-n? same-endpoints?)
731 (def any-number? naturals?)
732 ;; TODO compose
733 ;; TODO compose-n
734 ;; identity
735 (def compose-2 fan-in)
736 (def compose-bin fan-in*)
737 (def any? (constantly true))
738 (def none? (constantly false))
739 (def constant constantly)
740 (def joint-arity intersect)
741 (def a-reduce reduce)
742 ;; filter
743 (def make-map (partial partial map) )
744 (def bracket juxt)
745 ;; TODO apply-to-all
746 ;; TODO nary-combine
747 ;; TODO binary-combine
748 ;; TODO unary-combine
749 ;; iterated
750 ;; iterate-until-stable
751 (def make-function-of-vector (partial partial map))
752 (def make-function-of-arguments (fn [f] (fn [& args] (f args))))
753 (def alphaless lexical<)
755 #+end_src
763 :
765 * Numerical Methods
766 #+srcname: numerical-methods
767 #+begin_src clojure
768 (in-ns 'sicm.utils)
769 (import java.lang.Math)
770 (use 'clojure.contrib.def)
772 ;; ---- USEFUL CONSTANTS
774 (defn machine-epsilon
775 "Smallest value representable on your machine, as determined by
776 successively dividing a number in half until consecutive results are
777 indistinguishable."
778 []
779 (ffirst
780 (drop-while
781 (comp not zero? second)
782 (partition 2 1
783 (iterate (partial * 0.5) 1)))))
786 (def pi (Math/PI))
787 (def two-pi (* 2 pi))
789 (def eulers-gamma 0.5772156649015328606065)
791 (def phi (/ (inc (Math/sqrt 5)) 2))
793 (def ln2 (Math/log 2))
794 (def ln10 (Math/log 10))
795 (def exp10 #(Math/pow 10 %))
796 (def exp2 #(Math/pow 2 %))
799 ;;
801 ;; ---- ANGLES AND TRIGONOMETRY
803 (defn angle-restrictor
804 "Returns a function that ensures that angles lie in the specified interval of length two-pi."
805 [max-angle]
806 (let [min-angle (- max-angle two-pi)]
807 (fn [x]
808 (if (and
809 (<= min-angle x)
810 (< x max-angle))
811 x
812 (let [corrected-x (- x (* two-pi (Math/floor (/ x two-pi))))]
813 (if (< corrected-x max-angle)
814 corrected-x
815 (- corrected-x two-pi)))))))
817 (defn angle-restrict-pi
818 "Coerces angles to lie in the interval from -pi to pi."
819 [angle]
820 ((angle-restrictor pi) angle))
822 (defn angle-restrict-two-pi
823 "Coerces angles to lie in the interval from zero to two-pi"
824 [angle]
825 ((angle-restrictor two-pi) angle))
830 (defn invert [x] (/ x))
831 (defn negate [x] (- x))
833 (defn exp [x] (Math/exp x))
835 (defn sin [x] (Math/sin x))
836 (defn cos [x] (Math/cos x))
837 (defn tan [x] (Math/tan x))
839 (def sec (comp invert cos))
840 (def csc (comp invert sin))
842 (defn sinh [x] (Math/sinh x))
843 (defn cosh [x] (Math/cosh x))
844 (defn tanh [x] (Math/tanh x))
846 (def sech (comp invert cosh))
847 (def csch (comp invert sinh))
850 ;; ------------
852 (defn factorial
853 "Computes the factorial of the nonnegative integer n."
854 [n]
855 (if (neg? n)
856 (throw (Exception. "Cannot compute the factorial of a negative number."))
857 (reduce * 1 (range 1 (inc n)))))
859 (defn exact-quotient [n d] (/ n d))
861 (defn binomial-coefficient
862 "Computes the number of different ways to choose m elements from n."
863 [n m]
864 (assert (<= 0 m n))
865 (let [difference (- n m)]
866 (exact-quotient
867 (reduce * (range n (max difference m) -1 ))
868 (factorial (min difference m)))))
870 (defn-memo stirling-1
871 "Stirling numbers of the first kind: the number of permutations of n
872 elements with exactly m permutation cycles. "
873 [n k]
874 ;(assert (<= 1 k n))
875 (if (zero? n)
876 (if (zero? k) 1 0)
877 (+ (stirling-1 (dec n) (dec k))
878 (* (dec n) (stirling-1 (dec n) k)))))
880 (defn-memo stirling-2 ;;count-partitions
881 "Stirling numbers of the second kind: the number of ways to partition a set of n elements into k subsets."
882 [n k]
883 (cond
884 (= k 1) 1
885 (= k n) 1
886 :else (+ (stirling-2 (dec n) (dec k))
887 (* k (stirling-2 (dec n) k)))))
889 (defn harmonic-number [n]
890 (/ (stirling-1 (inc n) 2)
891 (factorial n)))
894 (defn sum
895 [f low high]
896 (reduce + (map f (range low (inc high)))))
898 #+end_src
911 * Differentiation
913 We compute derivatives by passing special *differential objects* $[x,
914 dx]$ through functions. Roughly speaking, applying a function $f$ to a
915 differential object \([x, dx]\) should produce a new differential
916 object $[f(x),\,Df(x)\cdot dx]$.
918 \([x,\,dx]\xrightarrow{\quad f \quad}[f(x),\,Df(x)\cdot dx]\)
919 Notice that you can obtain the derivative of $f$ from this
920 differential object, as it is the coefficient of the $dx$ term. Also,
921 as you apply successive functions using this rule, you get the
922 chain-rule answer you expect:
924 \([f(x),\,Df(x)\cdot dx]\xrightarrow{\quad g\quad} [gf(x),\,
925 Dgf(x)\cdot Df(x) \cdot dx ]\)
927 In order to generalize to multiple variables and multiple derivatives,
928 we use a *power series of differentials*, a sortred infinite sequence which
929 contains all terms like $dx\cdot dy$, $dx^2\cdot dy$, etc.
932 #+srcname:differential
933 #+begin_src clojure
934 (in-ns 'sicm.utils)
935 (use 'clojure.contrib.combinatorics)
936 (use 'clojure.contrib.generic.arithmetic)
938 (defprotocol DifferentialTerm
939 "Protocol for an infinitesimal quantity."
940 (coefficient [this])
941 (partials [this]))
943 (extend-protocol DifferentialTerm
944 java.lang.Number
945 (coefficient [this] this)
946 (partials [this] []))
948 (deftype DifferentialSeq
949 [terms]
950 clojure.lang.IPersistentCollection
951 (cons [this x]
952 (throw (Exception. x))
953 (DifferentialSeq. (cons x terms)))
954 ;; (seq [this] (seq terms))
955 (count [this] (count terms))
956 (empty [this] (empty? terms)))
968 (defn coerce-to-differential-seq [x]
969 (cond
970 (= (type x) DifferentialSeq) x
971 (satisfies? DifferentialTerm x) (DifferentialSeq. x)))
974 (defn differential-term
975 "Returns a differential term with the given coefficient and
976 partials. Coefficient may be any arithmetic object; partials must
977 be a list of non-negative integers."
978 [coefficient partials]
979 (reify DifferentialTerm
980 (partials [_] (set partials))
981 (coefficient [_] coefficient)))
985 (defn differential-seq*
986 ([coefficient partials]
987 (DifferentialSeq. [(differential-term coefficient partials)]))
988 ([coefficient partials & cps]
989 (if cps
994 (defn differential-seq
995 "Constructs a sequence of differential terms from a numerical
996 coefficient and a list of keys for variables. If no coefficient is supplied, uses 1."
997 ([variables] (differential-seq 1 variables))
998 ([coefficient variables & cvs]
999 (if (number? coefficient)
1000 (conj (assoc {} (apply sorted-set variables) coefficient)
1001 (if (empty? cvs)
1002 nil
1003 (apply differential-seq cvs)))
1004 (apply differential-seq 1 coefficient 1 variables cvs)
1005 )))
1008 (defn differential-add
1009 "Add two differential sequences by combining like terms."
1010 [dseq1 dseq2]
1011 (merge-with + dseq1 dseq2))
1013 (defn differential-multiply
1014 "Multiply two differential sequences. The square of any differential variable is zero since differential variables are infinitesimally small."
1015 [dseq1 dseq2]
1016 (reduce
1017 (fn [m [[vars1 coeff1] [vars2 coeff2]]]
1018 (if (empty? (clojure.set/intersection vars1 vars2))
1019 (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))
1020 m))
1021 {}
1022 (clojure.contrib.combinatorics/cartesian-product
1023 dseq1
1024 dseq2)))
1027 (defn big-part
1028 "Returns the part of the differential sequence that is finite,
1029 i.e. not infinitely small."
1030 [dseq]
1031 (let
1032 [keys (sort-by count (keys dseq))
1033 smallest-var (last(last keys))]
1034 (apply hash-map
1035 (reduce concat
1036 (remove (comp smallest-var first) dseq)))))
1038 (defn small-part
1039 "Returns the part of the differential sequence that is
1040 infinitesimal."
1041 [dseq]
1042 (let
1043 [keys (sort-by count (keys dseq))
1044 smallest-var (last(last keys))]
1045 (apply hash-map
1046 (reduce concat
1047 (filter (comp smallest-var first) dseq)))))
1049 (defn big-part
1050 "Returns the 'finite' part of the differential sequence."
1051 [dseq]
1052 (let
1053 [keys (sort-by count (keys dseq))
1054 smallest-var (last (last keys))]
1055 (apply hash-map
1056 (reduce concat
1057 (filter (remove smallest-var first) dseq)
1058 ))))
1061 (defn small-part
1062 "Returns the 'infinitesimal' part of the differential sequence."
1063 [dseq]
1064 (let
1065 [keys (sort-by count (keys dseq))
1066 smallest-var (last (last keys))]
1068 (apply hash-map
1069 (reduce concat
1070 (filter (comp smallest-var first) dseq)
1071 ))))
1073 (defn linear-approximator
1074 "Returns the function that corresponds to unary-fn but which accepts and returns differential objects."
1075 ([unary-f dfdx]
1076 (fn F [dseq]
1077 (differential-add
1078 (unary-f (big-part dseq))
1079 (differential-multiply
1080 (dfdx (big-part dseq))
1081 (small-part dseq))))))
1088 ;; ;; A differential term consists of a numerical coefficient and a
1089 ;; ;; sorted
1090 ;; (defrecord DifferentialTerm [coefficient variables])
1091 ;; (defmethod print-method DifferentialTerm
1092 ;; [o w]
1093 ;; (print-simple
1094 ;; (apply str (.coefficient o)(map (comp (partial str "d") name) (.variables o)))
1095 ;; w))
1098 ;; (defn differential-seq
1099 ;; "Constructs a sequence of differential terms from a numerical
1100 ;; coefficient and a list of keywords for variables. If no coefficient is
1101 ;; supplied, uses 1."
1102 ;; ([variables] (differential-seq 1 variables))
1103 ;; ([coefficient variables]
1104 ;; (list
1105 ;; (DifferentialTerm. coefficient (apply sorted-set variables))))
1106 ;; ([coefficient variables & cvs]
1107 ;; (sort-by
1108 ;; #(vec(.variables %))
1109 ;; (concat (differential-seq coefficient variables) (apply differential-seq cvs)))))
1111 #+end_src
1125 * Symbolic Manipulation
1127 #+srcname:symbolic
1128 #+begin_src clojure
1129 (in-ns 'sicm.utils)
1131 (deftype Symbolic [type expression]
1132 Object
1133 (equals [this that]
1134 (cond
1135 (= (.expression this) (.expression that)) true
1136 :else
1137 (Symbolic.
1138 java.lang.Boolean
1139 (list '= (.expression this) (.expression that)))
1140 )))
1145 (deftype AbstractSet [glyph membership-test])
1146 (defn member? [abstract-set x]
1147 ((.membership-test abstract-set) x))
1149 ;; ------------ Some important AbstractSets
1152 (def Real
1153 (AbstractSet.
1154 'R
1155 (fn[x](number? x))))
1158 ;; ------------ Create new AbstractSets from existing ones
1160 (defn abstract-product
1161 "Gives the cartesian product of abstract sets."
1162 ([sets]
1163 (if (= (count sets) 1) (first sets)
1164 (AbstractSet.
1165 (symbol
1166 (apply str
1167 (interpose 'x (map #(.glyph %) sets))))
1168 (fn [x]
1169 (and
1170 (coll? x)
1171 (= (count sets) (count x))
1172 (reduce #(and %1 %2)
1173 true
1174 (map #(member? %1 %2) sets x)))))))
1175 ([abstract-set n]
1176 (abstract-product (repeat n abstract-set))))
1180 (defn abstract-up
1181 "Returns the abstract set of all up-tuples whose items belong to the
1182 corresponding abstract sets in coll."
1183 ([coll]
1184 (AbstractSet.
1185 (symbol (str "u["
1186 (apply str
1187 (interpose " "
1188 (map #(.glyph %) coll)))
1189 "]"))
1190 (fn [x]
1191 (and
1192 (satisfies? Spinning x)
1193 (up? x)
1194 (= (count coll) (count x))
1195 (reduce
1196 #(and %1 %2)
1197 true
1198 (map #(member? %1 %2) coll x))))))
1199 ([abstract-set n]
1200 (abstract-up (repeat n abstract-set))))
1203 (defn abstract-down
1204 "Returns the abstract set of all down-tuples whose items belong to the
1205 corresponding abstract sets in coll."
1206 ([coll]
1207 (AbstractSet.
1208 (symbol (str "d["
1209 (apply str
1210 (interpose " "
1211 (map #(.glyph %) coll)))
1212 "]"))
1213 (fn [x]
1214 (and
1215 (satisfies? Spinning x)
1216 (down? x)
1217 (= (count coll) (count x))
1218 (reduce
1219 #(and %1 %2)
1220 true
1221 (map #(member? %1 %2) coll x))))))
1222 ([abstract-set n]
1223 (abstract-down (repeat n abstract-set))))
1229 ;-------ABSTRACT FUNCTIONS
1230 (defrecord AbstractFn
1231 [#^AbstractSet domain #^AbstractSet codomain])
1234 (defmethod print-method AbstractFn
1235 [o w]
1236 (print-simple
1237 (str
1238 "f:"
1239 (.glyph (:domain o))
1240 "-->"
1241 (.glyph (:codomain o))) w))
1242 #+end_src
1245 * COMMENT
1246 #+begin_src clojure :tangle utils.clj
1247 (ns sicm.utils)
1249 ;***** GENERIC ARITHMETIC
1250 <<generic-arithmetic>>
1253 ;***** TUPLES AND MATRICES
1254 <<tuples>>
1255 <<tuples-2>>
1256 ;***** MATRICES
1257 <<matrices>>
1259 #+end_src