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