(ns sicm.utils)

 

 					;***** GENERIC ARITHMETIC

 (ns sicm.utils)

 (in-ns 'sicm.utils)

 

 (defprotocol Arithmetic

   (zero [this])

   (one [this]))

 

 

 (extend-protocol Arithmetic

   java.lang.Number

   (zero [this] 0)

   (one [this] 1))

 

 (extend-protocol Arithmetic

   clojure.lang.Seqable

   (zero [this] (map zero this))

   (one [this] (map one this)))

 

 

 					;***** TUPLES AND MATRICES

 (in-ns 'sicm.utils)

 

 (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))

 

 

 (deftype Tuple

   [spin coll]

   clojure.lang.Seqable

   (seq [this] (seq (.coll this)))

   clojure.lang.Counted

   (count [this] (count (.coll this))))

 

 (extend-type Tuple

   Spinning

   (up? [this] (= ::up (.spin this)))

   (down? [this] (= ::down (.spin this))))

 

 (defmethod print-method Tuple

   [o w]

   (print-simple (str (if (up? o) 'u 'd) (.coll o))  w))

 

 

 

 (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))

 

 

 (in-ns 'sicm.utils)

 

 (defn numbers?

   "Returns true if all arguments are numbers, else false."

   [& xs]

   (every? number? xs))

 

 (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))))

 

 					;***** MATRICES

 (in-ns 'sicm.utils)

 (require 'incanter.core) ;; use incanter's fast matrices

 

 (defprotocol Matrix

   (rows [this])

   (cols [this])

   (diagonal [this])

   (trace [this])

   (determinant [this]))

 

 (extend-protocol Matrix

   incanter.Matrix

   (rows [this] (map down this)))

 

 

 

 

 (defn count-rows [matrix]

   ((comp count rows) matrix))

 

 (defn count-cols [matrix]

   ((comp count cols) matrix))

 

 

 (defn matrix-by-rows

   "Define a matrix by giving its rows."

   [& rows]

   (cond

    (not (all-equal? (map count rows)))

    (throw (Exception. "All rows in a matrix must have the same number of elements."))

    :else

    (reify Matrix

 	  (rows [this] (map down rows))

 	  (cols [this] (map up (apply map vector rows)))

 	  (diagonal [this] (map-indexed (fn [i row] (nth row i) rows)))

 	  (trace [this]

 		 (if (not= (count-rows this) (count-cols this))

 		   (throw (Exception.

 			   "Cannot take the trace of a non-square matrix."))

 		   (reduce + (diagonal this))))

 	  

 	  (determinant [this]

 		       (if (not= (count-rows this) (count-cols this))

 			 (throw (Exception.

 				 "Cannot take the determinant of a non-square matrix."))

 			 (reduce * (diagonal this))))

    )))

 

 

 (defn matrix-by-cols

   "Define a matrix by giving its columns."

   [& cols]

   (cond

    (not (all-equal? (map count cols)))

    (throw (Exception. "All columns in a matrix must have the same number of elements."))

    :else

    (reify Matrix

 	  (cols [this] (map up cols))

 	  (rows [this] (map down (apply map vector cols)))

 	  (diagonal [this] (map-indexed (fn [i col] (nth col i) cols)))

 	  (trace [this]

 		 (if (not= (count-cols this) (count-rows this))

 		   (throw (Exception.

 			   "Cannot take the trace of a non-square matrix."))

 		   (reduce + (diagonal this))))

 	  

 	  (determinant [this]

 		       (if (not= (count-cols this) (count-rows this))

 			 (throw (Exception.

 				 "Cannot take the determinant of a non-square matrix."))

 			 (reduce * (map-indexed (fn [i col] (nth col i)) cols))))

 	  )))

 

 (extend-protocol Matrix Tuple

 		 (rows [this] (if (down? this)

 				(list this)

 				(map (comp up vector) this)))

 

 		 (cols [this] (if (up? this)

 				(list this)

 				(map (comp down vector) this))))

 

 (defn matrix-multiply [A B]

   (apply matrix-by-rows

 	 (for [a (rows A)]

 	   (for [b (cols B)]

 	     (contract a b)))))