rlm@2: 
rlm@2: (ns sicm.utils)
rlm@2: 
rlm@2: (in-ns 'sicm.utils)
rlm@2: 
rlm@2: ;; Let some objects have spin
rlm@2: 
rlm@2: (defprotocol Spinning
rlm@2:   (up? [this])
rlm@2:   (down? [this]))
rlm@2: 
rlm@2: (defn spin
rlm@2:   "Returns the spin of the Spinning s, either :up or :down"
rlm@2:   [#^Spinning s]
rlm@2:   (cond (up? s) :up (down? s) :down))
rlm@2: 
rlm@2: 
rlm@2: ;; DEFINITION: A tuple is a sequence with spin
rlm@2: 
rlm@2: (deftype Tuple
rlm@2:   [spin coll]
rlm@2: 
rlm@2:   clojure.lang.Seqable
rlm@2:   (seq [this] (seq (.coll this)))
rlm@2: 
rlm@2:   clojure.lang.Counted
rlm@2:   (count [this] (count (.coll this)))
rlm@2: 
rlm@2:   Spinning
rlm@2:   (up? [this] (= ::up (.spin this)))
rlm@2:   (down? [this] (= ::down (.spin this))))
rlm@2: 
rlm@2: (defmethod print-method Tuple
rlm@2:   [o w]
rlm@2:   (print-simple
rlm@2:    (if (up? o)
rlm@2:      (str "u" (.coll o))
rlm@2:      (str "d" (vec(.coll o))))
rlm@2:    w))
rlm@2: 
rlm@2: (def tuple? #(= (type %) Tuple))
rlm@2: 
rlm@2: ;; CONSTRUCTORS
rlm@2: 
rlm@2: (defn up
rlm@2:   "Create a new up-tuple containing the contents of coll."
rlm@2:   [coll]
rlm@2:   (Tuple. ::up coll))       
rlm@2: 
rlm@2: (defn down
rlm@2:   "Create a new down-tuple containing the contents of coll."
rlm@2:   [coll]
rlm@2:   (Tuple. ::down coll))
rlm@2: 
rlm@2: (defn same-spin
rlm@2:   "Creates a tuple which has the same spin as tuple and which contains
rlm@2: the contents of coll."
rlm@2:   [tuple coll]
rlm@2:   (if (up? tuple)
rlm@2:     (up coll)
rlm@2:     (down coll)))
rlm@2: 
rlm@2: (defn opposite-spin
rlm@2:   "Create a tuple which has opposite spin to tuple and which contains
rlm@2: the contents of coll."
rlm@2:   [tuple coll]
rlm@2:   (if (up? tuple)
rlm@2:     (down coll)
rlm@2:     (up coll)))
rlm@2: (in-ns 'sicm.utils)
rlm@2: (require 'incanter.core) ;; use incanter's fast matrices
rlm@2: 
rlm@2: 
rlm@2: (defn all-equal? [coll]
rlm@2:   (if (empty? (rest coll)) true
rlm@2:       (and (= (first coll) (second coll))
rlm@2: 	   (recur (rest coll)))))
rlm@2: 
rlm@2: 
rlm@2: (defprotocol Matrix
rlm@2:   (rows [matrix])
rlm@2:   (cols [matrix])
rlm@2:   (diagonal [matrix])
rlm@2:   (trace [matrix])
rlm@2:   (determinant [matrix])
rlm@2:   (transpose [matrix])
rlm@2:   (conjugate [matrix])
rlm@2: )
rlm@2: 
rlm@2: (extend-protocol Matrix incanter.Matrix
rlm@2:   (rows [rs] (map down (apply map vector (apply map vector rs))))
rlm@2:   (cols [rs] (map up (apply map vector rs)))
rlm@2:   (diagonal [matrix] (incanter.core/diag matrix) )
rlm@2:   (determinant [matrix] (incanter.core/det matrix))
rlm@2:   (trace [matrix] (incanter.core/trace matrix))
rlm@2:   (transpose [matrix] (incanter.core/trans matrix)))
rlm@2: 
rlm@2: (defn count-rows [matrix]
rlm@2:   ((comp count rows) matrix))
rlm@2: 
rlm@2: (defn count-cols [matrix]
rlm@2:   ((comp count cols) matrix))
rlm@2: 
rlm@2: (defn square? [matrix]
rlm@2:   (= (count-rows matrix) (count-cols matrix)))
rlm@2: 
rlm@2: (defn identity-matrix
rlm@2:   "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2:   0s everywhere else."
rlm@2:   [n]
rlm@2:   (incanter.core/identity-matrix n))
rlm@2: 
rlm@2: 
rlm@2: (defn matrix-by-rows
rlm@2:   "Define a matrix by giving its rows."
rlm@2:   [& rows]
rlm@2:   (if
rlm@2:    (not (all-equal? (map count rows)))
rlm@2:    (throw (Exception. "All rows in a matrix must have the same number of elements."))
rlm@2:    (incanter.core/matrix (vec rows))))
rlm@2: 
rlm@2: (defn matrix-by-cols
rlm@2:   "Define a matrix by giving its columns"
rlm@2:   [& cols]
rlm@2:   (if (not (all-equal? (map count cols)))
rlm@2:    (throw (Exception. "All columns in a matrix must have the same number of elements."))
rlm@2:    (incanter.core/matrix (vec (apply map vector cols)))))
rlm@2: 
rlm@2: (defn identity-matrix
rlm@2:   "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2:   0s everywhere else."
rlm@2:   [n]
rlm@2:   (incanter.core/identity-matrix n))
rlm@2: 
rlm@2: (in-ns 'sicm.utils)
rlm@2: (use 'clojure.contrib.generic.arithmetic
rlm@2:      'clojure.contrib.generic.collection
rlm@2:      'clojure.contrib.generic.functor
rlm@2:      'clojure.contrib.generic.math-functions)
rlm@2: 
rlm@2: (defn numbers?
rlm@2:   "Returns true if all arguments are numbers, else false."
rlm@2:   [& xs]
rlm@2:   (every? number? xs))
rlm@2: 
rlm@2: (defn tuple-surgery
rlm@2:   "Applies the function f to the items of tuple and the additional
rlm@2:   arguments, if any. Returns a Tuple of the same type as tuple."
rlm@2:   [tuple f & xs]
rlm@2:   ((if (up? tuple) up down)
rlm@2:    (apply f (seq tuple) xs)))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: ;;; CONTRACTION collapses two compatible tuples into a number.
rlm@2: 
rlm@2: (defn contractible?
rlm@2:   "Returns true if the tuples a and b are compatible for contraction,
rlm@2:   else false. Tuples are compatible if they have the same number of
rlm@2:   components, they have opposite spins, and their elements are
rlm@2:   pairwise-compatible."
rlm@2:   [a b]
rlm@2:   (and
rlm@2:    (isa? (type a) Tuple)
rlm@2:    (isa? (type b) Tuple)
rlm@2:    (= (count a) (count b))
rlm@2:    (not= (spin a) (spin b))
rlm@2:    
rlm@2:    (not-any? false?
rlm@2: 	     (map #(or
rlm@2: 		    (numbers? %1 %2)
rlm@2: 		    (contractible? %1 %2))
rlm@2: 		  a b))))
rlm@2: 
rlm@2: (defn contract
rlm@2:   "Contracts two tuples, returning the sum of the
rlm@2:   products of the corresponding items. Contraction is recursive on
rlm@2:   nested tuples."
rlm@2:   [a b]
rlm@2:   (if (not (contractible? a b))
rlm@2:     (throw
rlm@2:      (Exception. "Not compatible for contraction."))
rlm@2:     (reduce +
rlm@2: 	    (map
rlm@2: 	     (fn [x y]
rlm@2: 	       (if (numbers? x y)
rlm@2: 		 (* x y)
rlm@2: 		 (contract x y)))
rlm@2: 	     a b))))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defmethod conj Tuple
rlm@2:   [tuple & xs]
rlm@2:   (tuple-surgery tuple #(apply conj % xs)))
rlm@2: 
rlm@2: (defmethod fmap Tuple
rlm@2:   [f tuple]
rlm@2:   (tuple-surgery tuple (partial map f)))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: ;; TODO: define Scalar, and add it to the hierarchy above Number and Complex
rlm@2: 
rlm@2: 					
rlm@2: (defmethod * [Tuple Tuple]             ; tuple*tuple
rlm@2:   [a b]
rlm@2:   (if (contractible? a b)
rlm@2:     (contract a b)
rlm@2:     (map (partial * a) b)))
rlm@2: 
rlm@2: 
rlm@2: (defmethod * [java.lang.Number Tuple]  ;; scalar *  tuple
rlm@2:   [a x] (fmap (partial * a) x))
rlm@2: 
rlm@2: (defmethod * [Tuple java.lang.Number]
rlm@2:   [x a] (* a x))
rlm@2: 
rlm@2: (defmethod * [java.lang.Number incanter.Matrix] ;; scalar *  matrix
rlm@2:   [x M] (incanter.core/mult x M))
rlm@2: 
rlm@2: (defmethod * [incanter.Matrix java.lang.Number]
rlm@2:   [M x] (* x M))
rlm@2: 
rlm@2: (defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
rlm@2:   [M1 M2]
rlm@2:   (incanter.core/mmult M1 M2))
rlm@2: 
rlm@2: (defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
rlm@2:   [M v]
rlm@2:   (if (and (apply numbers? v) (up? v)) 
rlm@2:     (* M (matrix-by-cols v))
rlm@2:     (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
rlm@2:     ))
rlm@2: 
rlm@2: (defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
rlm@2:   [v M]
rlm@2:   (if (and (apply numbers? v) (down? v))
rlm@2:     (* (matrix-by-rows v) M)
rlm@2:     (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
rlm@2:     ))
rlm@2: 
rlm@2: 
rlm@2: (defmethod exp incanter.Matrix
rlm@2:   [M]
rlm@2:   (incanter.core/exp M))
rlm@2: 
rlm@2: 
rlm@2: (in-ns 'sicm.utils)
rlm@2: (use 'clojure.contrib.seq
rlm@2:      'clojure.contrib.generic.arithmetic
rlm@2:      'clojure.contrib.generic.collection
rlm@2:      'clojure.contrib.generic.math-functions)
rlm@2: 
rlm@2: ;;∂
rlm@2: 
rlm@2: ;; DEFINITION : Differential Term
rlm@2: 
rlm@2: ;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
rlm@2: ;; coefficient of the quantity is returned by the 'coefficient' method,
rlm@2: ;; while the sequence of differential parameters is returned by the
rlm@2: ;; method 'partials'.
rlm@2: 
rlm@2: ;; Instead of using (potentially ambiguous) letters to denote
rlm@2: ;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].
rlm@2: 
rlm@2: ;; The coefficient can be any arithmetic object; the
rlm@2: ;; partials must be a nonrepeating sorted sequence of nonnegative
rlm@2: ;; integers.
rlm@2: 
rlm@2: (deftype DifferentialTerm [coefficient partials])
rlm@2: 
rlm@2: (defn differential-term
rlm@2:   "Make a differential term from a  coefficient and list of partials."
rlm@2:   [coefficient partials]
rlm@2:   (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials)) 
rlm@2:     (DifferentialTerm. coefficient (set partials))
rlm@2:     (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))
rlm@2: 
rlm@2: 
rlm@2: ;; DEFINITION : Differential Sequence
rlm@2: ;; A differential sequence is a sequence of differential terms, all with different partials.
rlm@2: ;; Internally, it is a map from the partials of each term to their coefficients.
rlm@2: 
rlm@2: (deftype DifferentialSeq
rlm@2:   [terms]
rlm@2:   ;;clojure.lang.IPersistentMap
rlm@2:   clojure.lang.Associative
rlm@2:   (assoc [this key val]
rlm@2:     (DifferentialSeq.
rlm@2:      (cons (differential-term val key) terms)))
rlm@2:   (cons [this x]
rlm@2: 	(DifferentialSeq. (cons x terms)))
rlm@2:   (containsKey [this key]
rlm@2: 	       (not(nil? (find-first #(= (.partials %) key) terms))))
rlm@2:   (count [this] (count (.terms this)))
rlm@2:   (empty [this] (DifferentialSeq. []))
rlm@2:   (entryAt [this key]
rlm@2: 	   ((juxt #(.partials %) #(.coefficient %))
rlm@2: 	    (find-first #(= (.partials %) key) terms)))
rlm@2:   (seq [this] (seq (.terms this))))
rlm@2: 
rlm@2: (def differential? #(= (type %) DifferentialSeq))
rlm@2: 
rlm@2: (defn zeroth-order?
rlm@2:   "Returns true if the differential sequence has at most a constant term."
rlm@2:   [dseq]
rlm@2:   (and
rlm@2:    (differential? dseq)
rlm@2:    (every?
rlm@2:     #(= #{} %)
rlm@2:     (keys (.terms dseq)))))
rlm@2: 
rlm@2: (defmethod fmap DifferentialSeq
rlm@2:   [f dseq]
rlm@2:   (DifferentialSeq.
rlm@2:    (fmap f (.terms dseq))))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: ;; BUILDING DIFFERENTIAL OBJECTS
rlm@2: 
rlm@2: (defn differential-seq
rlm@2:     "Define a differential sequence by specifying an alternating
rlm@2: sequence of coefficients and lists of partials."
rlm@2:   ([coefficient partials]
rlm@2:      (DifferentialSeq. {(set partials) coefficient}))
rlm@2:   ([coefficient partials & cps]
rlm@2:      (if (odd? (count cps))
rlm@2:        (throw (Exception. "differential-seq requires an even number of terms."))
rlm@2:        (DifferentialSeq.
rlm@2: 	(reduce
rlm@2: 	 #(assoc %1 (set (second %2)) (first %2))
rlm@2: 	 {(set partials) coefficient}
rlm@2: 	 (partition 2 cps))))))
rlm@2:   
rlm@2: 
rlm@2: 
rlm@2: (defn big-part
rlm@2:   "Returns the part of the differential sequence that is finite,
rlm@2:   i.e. not infinitely small. If the sequence is zeroth-order, returns
rlm@2:   the coefficient of the zeroth-order term instead. "
rlm@2:   [dseq]
rlm@2:   (if (zeroth-order? dseq) (get (.terms dseq) #{})
rlm@2:       (let [m (.terms dseq)
rlm@2: 	    keys (sort-by count (keys m))
rlm@2: 	    smallest-var (last (last keys))]
rlm@2: 	(DifferentialSeq.
rlm@2: 	 (reduce
rlm@2: 	  #(assoc %1 (first %2) (second %2))
rlm@2: 	  {}
rlm@2: 	  (remove #((first %) smallest-var) m))))))
rlm@2: 
rlm@2: 
rlm@2: (defn small-part
rlm@2:   "Returns the part of the differential sequence that infinitely
rlm@2:   small. If the sequence is zeroth-order, returns zero."
rlm@2:   [dseq]
rlm@2:   (if (zeroth-order? dseq) 0
rlm@2:       (let [m (.terms dseq)
rlm@2: 	    keys (sort-by count (keys m))
rlm@2: 	    smallest-var (last (last keys))]
rlm@2: 	(DifferentialSeq.
rlm@2: 	 (reduce
rlm@2: 	  #(assoc %1 (first %2) (second %2)) {}
rlm@2: 	  (filter #((first %) smallest-var) m))))))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defn cartesian-product [set1 set2]
rlm@2:   (reduce concat
rlm@2: 	  (for [x set1]
rlm@2: 	    (for [y set2]
rlm@2: 	      [x y]))))
rlm@2: 
rlm@2: (defn nth-subset [n]
rlm@2:   (if (zero? n) []
rlm@2:       (let [lg2 #(/ (log %) (log 2))
rlm@2: 	    k (int(java.lang.Math/floor (lg2 n)))
rlm@2: 	    ]
rlm@2: 	(cons k
rlm@2: 	 (nth-subset  (- n (pow 2 k)))))))
rlm@2:     
rlm@2: (def all-partials
rlm@2:      (lazy-seq (map nth-subset (range))))
rlm@2: 
rlm@2: 
rlm@2: (defn differential-multiply
rlm@2:   "Multiply two differential sequences. The square of any differential
rlm@2:   variable is zero since differential variables are infinitesimally
rlm@2:   small."
rlm@2:   [dseq1 dseq2]
rlm@2:   (DifferentialSeq.
rlm@2:    (reduce
rlm@2:     (fn [m [[vars1 coeff1] [vars2 coeff2]]]
rlm@2:       (if (not (empty? (clojure.set/intersection vars1 vars2)))
rlm@2: 	m
rlm@2: 	(assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
rlm@2:     {}
rlm@2:     (cartesian-product (.terms dseq1) (.terms dseq2)))))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defmethod * [DifferentialSeq DifferentialSeq]
rlm@2:   [dseq1 dseq2]
rlm@2:    (differential-multiply dseq1 dseq2))
rlm@2: 
rlm@2: (defmethod + [DifferentialSeq DifferentialSeq]
rlm@2:   [dseq1 dseq2]
rlm@2:   (DifferentialSeq.
rlm@2:    (merge-with + (.terms dseq1) (.terms dseq2))))
rlm@2: 
rlm@2: (defmethod * [java.lang.Number DifferentialSeq]
rlm@2:   [x dseq]
rlm@2:   (fmap (partial * x) dseq))
rlm@2: 
rlm@2: (defmethod * [DifferentialSeq java.lang.Number]
rlm@2:   [dseq x]
rlm@2:   (fmap (partial * x) dseq))
rlm@2: 
rlm@2: (defmethod + [java.lang.Number DifferentialSeq]
rlm@2:   [x dseq]
rlm@2:   (+ (differential-seq x []) dseq))
rlm@2: (defmethod + [DifferentialSeq java.lang.Number]
rlm@2:   [dseq x]
rlm@2:   (+ dseq (differential-seq x [])))
rlm@2: 
rlm@2: (defmethod - DifferentialSeq
rlm@2:   [x]
rlm@2:   (fmap - x))
rlm@2: 
rlm@2: 
rlm@2: ;; DERIVATIVES
rlm@2: 
rlm@2: 	      
rlm@2: 
rlm@2: (defn linear-approximator
rlm@2:   "Returns an operator that linearly approximates the given function."
rlm@2:   ([f df|dx]
rlm@2:       (fn [x]
rlm@2: 	(let [big-part (big-part x)
rlm@2: 	      small-part (small-part x)]
rlm@2: 	  ;; f(x+dx) ~= f(x) + f'(x)dx
rlm@2: 	  (+ (f big-part)
rlm@2: 	     (* (df|dx big-part) small-part)
rlm@2: 	  ))))
rlm@2:      
rlm@2:   ([f df|dx df|dy]
rlm@2:      (fn [x y]
rlm@2:        (let [X (big-part x)
rlm@2: 	     Y (big-part y)
rlm@2: 	     DX (small-part x)
rlm@2: 	     DY (small-part y)]
rlm@2: 	 (+ (f X Y)
rlm@2: 	    (* DX (f df|dx X Y))
rlm@2: 	    (* DY (f df|dy X Y)))))))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defn D[f]
rlm@2:   (fn[x] (f (+ x (differential-seq  1 [0] 1 [1] 1 [2])))))
rlm@2:   
rlm@2: (defn d[partials f]
rlm@2:   (fn [x]
rlm@2:     (get 
rlm@2:      (.terms ((D f)x))
rlm@2:      (set partials)
rlm@2:      0
rlm@2:     )))
rlm@2: 
rlm@2: (defmethod exp DifferentialSeq [x]
rlm@2: 	   ((linear-approximator exp exp) x))
rlm@2: 
rlm@2: (defmethod sin DifferentialSeq
rlm@2:   [x]
rlm@2:   ((linear-approximator sin cos) x))
rlm@2: 
rlm@2: (defmethod cos DifferentialSeq
rlm@2:   [x]
rlm@2:   ((linear-approximator cos #(- (sin %))) x))
rlm@2: 
rlm@2: (defmethod log DifferentialSeq
rlm@2:   [x]
rlm@2:   ((linear-approximator log (fn [x] (/ x)) ) x))
rlm@2: 
rlm@2: (defmethod / [DifferentialSeq DifferentialSeq]
rlm@2:   [x y]
rlm@2:   ((linear-approximator /
rlm@2: 			(fn [x y] (/ 1 y))
rlm@2: 			(fn [x y] (- (/ x (* y y)))))
rlm@2:    x y))