dylan: sicm/utils.clj comparison

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

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

Mercurial > dylan

comparison sicm/utils.clj @ 2:b4de894a1e2e