rlm@2: 
rlm@2: (ns sicm.utils)
rlm@2: 
rlm@2: 					;***** GENERIC ARITHMETIC
rlm@2: (ns sicm.utils)
rlm@2: (in-ns 'sicm.utils)
rlm@2: 
rlm@2: (defprotocol Arithmetic
rlm@2:   (zero [this])
rlm@2:   (one [this]))
rlm@2: 
rlm@2: 
rlm@2: (extend-protocol Arithmetic
rlm@2:   java.lang.Number
rlm@2:   (zero [this] 0)
rlm@2:   (one [this] 1))
rlm@2: 
rlm@2: (extend-protocol Arithmetic
rlm@2:   clojure.lang.Seqable
rlm@2:   (zero [this] (map zero this))
rlm@2:   (one [this] (map one this)))
rlm@2: 
rlm@2: 
rlm@2: 					;***** TUPLES AND MATRICES
rlm@2: (in-ns 'sicm.utils)
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: (deftype Tuple
rlm@2:   [spin coll]
rlm@2:   clojure.lang.Seqable
rlm@2:   (seq [this] (seq (.coll this)))
rlm@2:   clojure.lang.Counted
rlm@2:   (count [this] (count (.coll this))))
rlm@2: 
rlm@2: (extend-type Tuple
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 (str (if (up? o) 'u 'd) (.coll o))  w))
rlm@2: 
rlm@2: 
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: 
rlm@2: (in-ns 'sicm.utils)
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 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: 
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: 					;***** MATRICES
rlm@2: (in-ns 'sicm.utils)
rlm@2: (require 'incanter.core) ;; use incanter's fast matrices
rlm@2: 
rlm@2: (defprotocol Matrix
rlm@2:   (rows [this])
rlm@2:   (cols [this])
rlm@2:   (diagonal [this])
rlm@2:   (trace [this])
rlm@2:   (determinant [this]))
rlm@2: 
rlm@2: (extend-protocol Matrix
rlm@2:   incanter.Matrix
rlm@2:   (rows [this] (map down this)))
rlm@2: 
rlm@2: 
rlm@2: 
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: 
rlm@2: (defn matrix-by-rows
rlm@2:   "Define a matrix by giving its rows."
rlm@2:   [& rows]
rlm@2:   (cond
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:    :else
rlm@2:    (reify Matrix
rlm@2: 	  (rows [this] (map down rows))
rlm@2: 	  (cols [this] (map up (apply map vector rows)))
rlm@2: 	  (diagonal [this] (map-indexed (fn [i row] (nth row i) rows)))
rlm@2: 	  (trace [this]
rlm@2: 		 (if (not= (count-rows this) (count-cols this))
rlm@2: 		   (throw (Exception.
rlm@2: 			   "Cannot take the trace of a non-square matrix."))
rlm@2: 		   (reduce + (diagonal this))))
rlm@2: 	  
rlm@2: 	  (determinant [this]
rlm@2: 		       (if (not= (count-rows this) (count-cols this))
rlm@2: 			 (throw (Exception.
rlm@2: 				 "Cannot take the determinant of a non-square matrix."))
rlm@2: 			 (reduce * (diagonal this))))
rlm@2:    )))
rlm@2: 
rlm@2: 
rlm@2: (defn matrix-by-cols
rlm@2:   "Define a matrix by giving its columns."
rlm@2:   [& cols]
rlm@2:   (cond
rlm@2:    (not (all-equal? (map count cols)))
rlm@2:    (throw (Exception. "All columns in a matrix must have the same number of elements."))
rlm@2:    :else
rlm@2:    (reify Matrix
rlm@2: 	  (cols [this] (map up cols))
rlm@2: 	  (rows [this] (map down (apply map vector cols)))
rlm@2: 	  (diagonal [this] (map-indexed (fn [i col] (nth col i) cols)))
rlm@2: 	  (trace [this]
rlm@2: 		 (if (not= (count-cols this) (count-rows this))
rlm@2: 		   (throw (Exception.
rlm@2: 			   "Cannot take the trace of a non-square matrix."))
rlm@2: 		   (reduce + (diagonal this))))
rlm@2: 	  
rlm@2: 	  (determinant [this]
rlm@2: 		       (if (not= (count-cols this) (count-rows this))
rlm@2: 			 (throw (Exception.
rlm@2: 				 "Cannot take the determinant of a non-square matrix."))
rlm@2: 			 (reduce * (map-indexed (fn [i col] (nth col i)) cols))))
rlm@2: 	  )))
rlm@2: 
rlm@2: (extend-protocol Matrix Tuple
rlm@2: 		 (rows [this] (if (down? this)
rlm@2: 				(list this)
rlm@2: 				(map (comp up vector) this)))
rlm@2: 
rlm@2: 		 (cols [this] (if (up? this)
rlm@2: 				(list this)
rlm@2: 				(map (comp down vector) this))))
rlm@2: 
rlm@2: (defn matrix-multiply [A B]
rlm@2:   (apply matrix-by-rows
rlm@2: 	 (for [a (rows A)]
rlm@2: 	   (for [b (cols B)]
rlm@2: 	     (contract a b)))))