1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/sicm/utils.org	Fri Oct 28 00:03:05 2011 -0700
     1.3 @@ -0,0 +1,690 @@
     1.4 +#+TITLE: Building a Classical Mechanics Library in Clojure
     1.5 +#+author: Dylan Holmes
     1.6 +#+EMAIL:     rlm@mit.edu
     1.7 +#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js"
     1.8 +#+STYLE: <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />
     1.9 +#+OPTIONS:   H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
    1.10 +#+SETUPFILE: ../templates/level-0.org
    1.11 +#+INCLUDE: ../templates/level-0.org
    1.12 +#+BABEL: :noweb yes :results silent
    1.13 +
    1.14 +* Useful data types
    1.15 +
    1.16 +** Complex numbers
    1.17 +
    1.18 +** Power series
    1.19 +
    1.20 +** Tuples and tensors
    1.21 +
    1.22 +*** Tuples are \ldquo{}sequences with spin\rdquo{}
    1.23 +
    1.24 +#+srcname: tuples
    1.25 +#+begin_src clojure
    1.26 +(in-ns 'sicm.utils)
    1.27 +
    1.28 +;; Let some objects have spin
    1.29 +
    1.30 +(defprotocol Spinning
    1.31 +  (up? [this])
    1.32 +  (down? [this]))
    1.33 +
    1.34 +(defn spin
    1.35 +  "Returns the spin of the Spinning s, either :up or :down"
    1.36 +  [#^Spinning s]
    1.37 +  (cond (up? s) :up (down? s) :down))
    1.38 +
    1.39 +
    1.40 +;; DEFINITION: A tuple is a sequence with spin
    1.41 +
    1.42 +(deftype Tuple
    1.43 +  [spin coll]
    1.44 +
    1.45 +  clojure.lang.Seqable
    1.46 +  (seq [this] (seq (.coll this)))
    1.47 +
    1.48 +  clojure.lang.Counted
    1.49 +  (count [this] (count (.coll this)))
    1.50 +
    1.51 +  Spinning
    1.52 +  (up? [this] (= ::up (.spin this)))
    1.53 +  (down? [this] (= ::down (.spin this))))
    1.54 +
    1.55 +(defmethod print-method Tuple
    1.56 +  [o w]
    1.57 +  (print-simple
    1.58 +   (if (up? o)
    1.59 +     (str "u" (.coll o))
    1.60 +     (str "d" (vec(.coll o))))
    1.61 +   w))
    1.62 +
    1.63 +(def tuple? #(= (type %) Tuple))
    1.64 +
    1.65 +;; CONSTRUCTORS
    1.66 +
    1.67 +(defn up
    1.68 +  "Create a new up-tuple containing the contents of coll."
    1.69 +  [coll]
    1.70 +  (Tuple. ::up coll))       
    1.71 +
    1.72 +(defn down
    1.73 +  "Create a new down-tuple containing the contents of coll."
    1.74 +  [coll]
    1.75 +  (Tuple. ::down coll))
    1.76 +
    1.77 +(defn same-spin
    1.78 +  "Creates a tuple which has the same spin as tuple and which contains
    1.79 +the contents of coll."
    1.80 +  [tuple coll]
    1.81 +  (if (up? tuple)
    1.82 +    (up coll)
    1.83 +    (down coll)))
    1.84 +
    1.85 +(defn opposite-spin
    1.86 +  "Create a tuple which has opposite spin to tuple and which contains
    1.87 +the contents of coll."
    1.88 +  [tuple coll]
    1.89 +  (if (up? tuple)
    1.90 +    (down coll)
    1.91 +    (up coll)))
    1.92 +#+end_src
    1.93 +
    1.94 +
    1.95 +
    1.96 +*** Matrices
    1.97 +#+srcname:matrices
    1.98 +#+begin_src clojure
    1.99 +(in-ns 'sicm.utils)
   1.100 +(require 'incanter.core) ;; use incanter's fast matrices
   1.101 +
   1.102 +
   1.103 +(defn all-equal? [coll]
   1.104 +  (if (empty? (rest coll)) true
   1.105 +      (and (= (first coll) (second coll))
   1.106 +	   (recur (rest coll)))))
   1.107 +
   1.108 +
   1.109 +(defprotocol Matrix
   1.110 +  (rows [matrix])
   1.111 +  (cols [matrix])
   1.112 +  (diagonal [matrix])
   1.113 +  (trace [matrix])
   1.114 +  (determinant [matrix])
   1.115 +  (transpose [matrix])
   1.116 +  (conjugate [matrix])
   1.117 +)
   1.118 +
   1.119 +(extend-protocol Matrix incanter.Matrix
   1.120 +  (rows [rs] (map down (apply map vector (apply map vector rs))))
   1.121 +  (cols [rs] (map up (apply map vector rs)))
   1.122 +  (diagonal [matrix] (incanter.core/diag matrix) )
   1.123 +  (determinant [matrix] (incanter.core/det matrix))
   1.124 +  (trace [matrix] (incanter.core/trace matrix))
   1.125 +  (transpose [matrix] (incanter.core/trans matrix)))
   1.126 +
   1.127 +(defn count-rows [matrix]
   1.128 +  ((comp count rows) matrix))
   1.129 +
   1.130 +(defn count-cols [matrix]
   1.131 +  ((comp count cols) matrix))
   1.132 +
   1.133 +(defn square? [matrix]
   1.134 +  (= (count-rows matrix) (count-cols matrix)))
   1.135 +
   1.136 +(defn identity-matrix
   1.137 +  "Define a square matrix of size n-by-n with 1s along the diagonal and
   1.138 +  0s everywhere else."
   1.139 +  [n]
   1.140 +  (incanter.core/identity-matrix n))
   1.141 +
   1.142 +
   1.143 +(defn matrix-by-rows
   1.144 +  "Define a matrix by giving its rows."
   1.145 +  [& rows]
   1.146 +  (if
   1.147 +   (not (all-equal? (map count rows)))
   1.148 +   (throw (Exception. "All rows in a matrix must have the same number of elements."))
   1.149 +   (incanter.core/matrix (vec rows))))
   1.150 +
   1.151 +(defn matrix-by-cols
   1.152 +  "Define a matrix by giving its columns"
   1.153 +  [& cols]
   1.154 +  (if (not (all-equal? (map count cols)))
   1.155 +   (throw (Exception. "All columns in a matrix must have the same number of elements."))
   1.156 +   (incanter.core/matrix (vec (apply map vector cols)))))
   1.157 +
   1.158 +(defn identity-matrix
   1.159 +  "Define a square matrix of size n-by-n with 1s along the diagonal and
   1.160 +  0s everywhere else."
   1.161 +  [n]
   1.162 +  (incanter.core/identity-matrix n))
   1.163 +
   1.164 +#+end_src
   1.165 +
   1.166 +** Generic Operations
   1.167 +#+srcname:arith-tuple
   1.168 +#+begin_src clojure
   1.169 +(in-ns 'sicm.utils)
   1.170 +(use 'clojure.contrib.generic.arithmetic
   1.171 +     'clojure.contrib.generic.collection
   1.172 +     'clojure.contrib.generic.functor
   1.173 +     'clojure.contrib.generic.math-functions)
   1.174 +
   1.175 +(defn numbers?
   1.176 +  "Returns true if all arguments are numbers, else false."
   1.177 +  [& xs]
   1.178 +  (every? number? xs))
   1.179 +
   1.180 +(defn tuple-surgery
   1.181 +  "Applies the function f to the items of tuple and the additional
   1.182 +  arguments, if any. Returns a Tuple of the same type as tuple."
   1.183 +  [tuple f & xs]
   1.184 +  ((if (up? tuple) up down)
   1.185 +   (apply f (seq tuple) xs)))
   1.186 +
   1.187 +
   1.188 +
   1.189 +;;; CONTRACTION collapses two compatible tuples into a number.
   1.190 +
   1.191 +(defn contractible?
   1.192 +  "Returns true if the tuples a and b are compatible for contraction,
   1.193 +  else false. Tuples are compatible if they have the same number of
   1.194 +  components, they have opposite spins, and their elements are
   1.195 +  pairwise-compatible."
   1.196 +  [a b]
   1.197 +  (and
   1.198 +   (isa? (type a) Tuple)
   1.199 +   (isa? (type b) Tuple)
   1.200 +   (= (count a) (count b))
   1.201 +   (not= (spin a) (spin b))
   1.202 +   
   1.203 +   (not-any? false?
   1.204 +	     (map #(or
   1.205 +		    (numbers? %1 %2)
   1.206 +		    (contractible? %1 %2))
   1.207 +		  a b))))
   1.208 +
   1.209 +(defn contract
   1.210 +  "Contracts two tuples, returning the sum of the
   1.211 +  products of the corresponding items. Contraction is recursive on
   1.212 +  nested tuples."
   1.213 +  [a b]
   1.214 +  (if (not (contractible? a b))
   1.215 +    (throw
   1.216 +     (Exception. "Not compatible for contraction."))
   1.217 +    (reduce +
   1.218 +	    (map
   1.219 +	     (fn [x y]
   1.220 +	       (if (numbers? x y)
   1.221 +		 (* x y)
   1.222 +		 (contract x y)))
   1.223 +	     a b))))
   1.224 +
   1.225 +
   1.226 +
   1.227 +
   1.228 +
   1.229 +(defmethod conj Tuple
   1.230 +  [tuple & xs]
   1.231 +  (tuple-surgery tuple #(apply conj % xs)))
   1.232 +
   1.233 +(defmethod fmap Tuple
   1.234 +  [f tuple]
   1.235 +  (tuple-surgery tuple (partial map f)))
   1.236 +
   1.237 +
   1.238 +
   1.239 +;; TODO: define Scalar, and add it to the hierarchy above Number and Complex
   1.240 +
   1.241 +					
   1.242 +(defmethod * [Tuple Tuple]             ; tuple*tuple
   1.243 +  [a b]
   1.244 +  (if (contractible? a b)
   1.245 +    (contract a b)
   1.246 +    (map (partial * a) b)))
   1.247 +
   1.248 +
   1.249 +(defmethod * [java.lang.Number Tuple]  ;; scalar *  tuple
   1.250 +  [a x] (fmap (partial * a) x))
   1.251 +
   1.252 +(defmethod * [Tuple java.lang.Number]
   1.253 +  [x a] (* a x))
   1.254 +
   1.255 +(defmethod * [java.lang.Number incanter.Matrix] ;; scalar *  matrix
   1.256 +  [x M] (incanter.core/mult x M))
   1.257 +
   1.258 +(defmethod * [incanter.Matrix java.lang.Number]
   1.259 +  [M x] (* x M))
   1.260 +
   1.261 +(defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
   1.262 +  [M1 M2]
   1.263 +  (incanter.core/mmult M1 M2))
   1.264 +
   1.265 +(defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
   1.266 +  [M v]
   1.267 +  (if (and (apply numbers? v) (up? v)) 
   1.268 +    (* M (matrix-by-cols v))
   1.269 +    (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
   1.270 +    ))
   1.271 +
   1.272 +(defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
   1.273 +  [v M]
   1.274 +  (if (and (apply numbers? v) (down? v))
   1.275 +    (* (matrix-by-rows v) M)
   1.276 +    (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
   1.277 +    ))
   1.278 +
   1.279 +
   1.280 +(defmethod exp incanter.Matrix
   1.281 +  [M]
   1.282 +  (incanter.core/exp M))
   1.283 +
   1.284 +#+end_src
   1.285 +
   1.286 +* Operators and Differentiation
   1.287 +** Operators
   1.288 +#+scrname: operators
   1.289 +#+begin_src clojure
   1.290 +(in-ns 'sicm.utils)
   1.291 +(use 'clojure.contrib.seq
   1.292 +     'clojure.contrib.generic.arithmetic
   1.293 +     'clojure.contrib.generic.collection
   1.294 +     'clojure.contrib.generic.math-functions)
   1.295 +
   1.296 +(defmethod + [clojure.lang.IFn clojure.lang.IFn]
   1.297 +	   [f g]
   1.298 +	   (fn [& args]
   1.299 +	     (+ (apply f args) (apply g args))))
   1.300 +
   1.301 +(defmethod * [clojure.lang.IFn clojure.lang.IFn]
   1.302 +	   [f g]
   1.303 +	   (fn [& args]
   1.304 +	     (* (apply f args) (apply g args))))
   1.305 +
   1.306 +(defmethod / [clojure.lang.IFn java.lang.Number]
   1.307 +  [f x]
   1.308 +  (fn [& args]
   1.309 +    (/ (apply f args) x)))
   1.310 +
   1.311 +
   1.312 +(defmethod - [clojure.lang.IFn]
   1.313 +  [f]
   1.314 +  (fn [& args]
   1.315 +    (- (apply f args))))
   1.316 +
   1.317 +(defmethod - [clojure.lang.IFn clojure.lang.IFn]
   1.318 +  [f g]
   1.319 +  (fn [& args]
   1.320 +    (- (apply f args) (apply g args))))
   1.321 +
   1.322 +(defmethod pow [clojure.lang.IFn java.lang.Number]
   1.323 +  [f x]
   1.324 +  (fn [& args]
   1.325 +    (pow (apply f args) x)))
   1.326 +
   1.327 +
   1.328 +(defmethod + [java.lang.Number clojure.lang.IFn]
   1.329 +	   [x f]
   1.330 +	   (fn [& args]
   1.331 +	     (+ x (apply f args))))
   1.332 +
   1.333 +(defmethod * [java.lang.Number clojure.lang.IFn]
   1.334 +	   [x f]
   1.335 +	   (fn [& args]
   1.336 +	     (* x (apply f args))))
   1.337 +
   1.338 +(defmethod * [clojure.lang.IFn java.lang.Number]
   1.339 +	   [f x]
   1.340 +	   (* x f))
   1.341 +(defmethod + [clojure.lang.IFn java.lang.Number]
   1.342 +	   [f x]
   1.343 +	   (+ x f))
   1.344 +
   1.345 +#+end_src
   1.346 +
   1.347 +** Differential Terms and Sequences
   1.348 +#+srcname: differential
   1.349 +#+begin_src clojure
   1.350 +(in-ns 'sicm.utils)
   1.351 +(use 'clojure.contrib.seq
   1.352 +     'clojure.contrib.generic.arithmetic
   1.353 +     'clojure.contrib.generic.collection
   1.354 +     'clojure.contrib.generic.math-functions)
   1.355 +
   1.356 +;;∂
   1.357 +
   1.358 +;; DEFINITION : Differential Term
   1.359 +
   1.360 +;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
   1.361 +;; coefficient of the quantity is returned by the 'coefficient' method,
   1.362 +;; while the sequence of differential parameters is returned by the
   1.363 +;; method 'partials'.
   1.364 +
   1.365 +;; Instead of using (potentially ambiguous) letters to denote
   1.366 +;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].
   1.367 +
   1.368 +;; The coefficient can be any arithmetic object; the
   1.369 +;; partials must be a nonrepeating sorted sequence of nonnegative
   1.370 +;; integers.
   1.371 +
   1.372 +;; (deftype DifferentialTerm [coefficient partials])
   1.373 +
   1.374 +;; (defn differential-term
   1.375 +;;   "Make a differential term from a  coefficient and list of partials."
   1.376 +;;   [coefficient partials]
   1.377 +;;   (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials)) 
   1.378 +;;     (DifferentialTerm. coefficient (set partials))
   1.379 +;;     (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))
   1.380 +
   1.381 +
   1.382 +;; DEFINITION : Differential Sequence
   1.383 +;; A differential sequence is a sequence of differential terms, all with different partials.
   1.384 +;; Internally, it is a map from the partials of each term to their coefficients.
   1.385 +
   1.386 +(deftype DifferentialSeq
   1.387 +  [terms]
   1.388 +  ;;clojure.lang.IPersistentMap
   1.389 +  clojure.lang.Associative
   1.390 +  (assoc [this key val]
   1.391 +    (DifferentialSeq.
   1.392 +     (cons (differential-term val key) terms)))
   1.393 +  (cons [this x]
   1.394 +	(DifferentialSeq. (cons x terms)))
   1.395 +  (containsKey [this key]
   1.396 +	       (not(nil? (find-first #(= (.partials %) key) terms))))
   1.397 +  (count [this] (count (.terms this)))
   1.398 +  (empty [this] (DifferentialSeq. []))
   1.399 +  (entryAt [this key]
   1.400 +	   ((juxt #(.partials %) #(.coefficient %))
   1.401 +	    (find-first #(= (.partials %) key) terms)))
   1.402 +  (seq [this] (seq (.terms this))))
   1.403 +
   1.404 +(def differential? #(= (type %) DifferentialSeq))
   1.405 +
   1.406 +(defn zeroth-order?
   1.407 +  "Returns true if the differential sequence has at most a constant term."
   1.408 +  [dseq]
   1.409 +  (and
   1.410 +   (differential? dseq)
   1.411 +   (every?
   1.412 +    #(= #{} %)
   1.413 +    (keys (.terms dseq)))))
   1.414 +
   1.415 +(defmethod fmap DifferentialSeq
   1.416 +  [f dseq]
   1.417 +  (DifferentialSeq.
   1.418 +   (fmap f (.terms dseq))))
   1.419 +
   1.420 +
   1.421 +
   1.422 +
   1.423 +;; BUILDING DIFFERENTIAL OBJECTS
   1.424 +
   1.425 +(defn differential-seq
   1.426 +    "Define a differential sequence by specifying an alternating
   1.427 +sequence of coefficients and lists of partials."
   1.428 +  ([coefficient partials]
   1.429 +     (DifferentialSeq. {(set partials) coefficient}))
   1.430 +  ([coefficient partials & cps]
   1.431 +     (if (odd? (count cps))
   1.432 +       (throw (Exception. "differential-seq requires an even number of terms."))
   1.433 +       (DifferentialSeq.
   1.434 +	(reduce
   1.435 +	 #(assoc %1 (set (second %2)) (first %2))
   1.436 +	 {(set partials) coefficient}
   1.437 +	 (partition 2 cps))))))
   1.438 +  
   1.439 +
   1.440 +
   1.441 +(defn big-part
   1.442 +  "Returns the part of the differential sequence that is finite,
   1.443 +  i.e. not infinitely small. If the sequence is zeroth-order, returns
   1.444 +  the coefficient of the zeroth-order term instead. "
   1.445 +  [dseq]
   1.446 +  (if (zeroth-order? dseq) (get (.terms dseq) #{})
   1.447 +      (let [m (.terms dseq)
   1.448 +	    keys (sort-by count (keys m))
   1.449 +	    smallest-var (last (last keys))]
   1.450 +	(DifferentialSeq.
   1.451 +	 (reduce
   1.452 +	  #(assoc %1 (first %2) (second %2))
   1.453 +	  {}
   1.454 +	  (remove #((first %) smallest-var) m))))))
   1.455 +
   1.456 +
   1.457 +(defn small-part
   1.458 +  "Returns the part of the differential sequence that infinitely
   1.459 +  small. If the sequence is zeroth-order, returns zero."
   1.460 +  [dseq]
   1.461 +  (if (zeroth-order? dseq) 0
   1.462 +      (let [m (.terms dseq)
   1.463 +	    keys (sort-by count (keys m))
   1.464 +	    smallest-var (last (last keys))]
   1.465 +	(DifferentialSeq.
   1.466 +	 (reduce
   1.467 +	  #(assoc %1 (first %2) (second %2)) {}
   1.468 +	  (filter #((first %) smallest-var) m))))))
   1.469 +
   1.470 +
   1.471 +
   1.472 +(defn cartesian-product [set1 set2]
   1.473 +  (reduce concat
   1.474 +	  (for [x set1]
   1.475 +	    (for [y set2]
   1.476 +	      [x y]))))
   1.477 +
   1.478 +(defn nth-subset [n]
   1.479 +  (if (zero? n) []
   1.480 +      (let [lg2 #(/ (log %) (log 2))
   1.481 +	    k (int(java.lang.Math/floor (lg2 n)))
   1.482 +	    ]
   1.483 +	(cons k
   1.484 +	 (nth-subset  (- n (pow 2 k)))))))
   1.485 +    
   1.486 +(def all-partials
   1.487 +     (lazy-seq (map nth-subset (range))))
   1.488 +
   1.489 +
   1.490 +(defn differential-multiply
   1.491 +  "Multiply two differential sequences. The square of any differential
   1.492 +  variable is zero since differential variables are infinitesimally
   1.493 +  small."
   1.494 +  [dseq1 dseq2]
   1.495 +  (DifferentialSeq.
   1.496 +   (reduce
   1.497 +    (fn [m [[vars1 coeff1] [vars2 coeff2]]]
   1.498 +      (if (not (empty? (clojure.set/intersection vars1 vars2)))
   1.499 +	m
   1.500 +	(assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
   1.501 +    {}
   1.502 +    (cartesian-product (.terms dseq1) (.terms dseq2)))))
   1.503 +
   1.504 +
   1.505 +
   1.506 +(defmethod * [DifferentialSeq DifferentialSeq]
   1.507 +  [dseq1 dseq2]
   1.508 +   (differential-multiply dseq1 dseq2))
   1.509 +
   1.510 +(defmethod + [DifferentialSeq DifferentialSeq]
   1.511 +  [dseq1 dseq2]
   1.512 +  (DifferentialSeq.
   1.513 +   (merge-with + (.terms dseq1) (.terms dseq2))))
   1.514 +
   1.515 +(defmethod * [java.lang.Number DifferentialSeq]
   1.516 +  [x dseq]
   1.517 +  (fmap (partial * x) dseq))
   1.518 +
   1.519 +(defmethod * [DifferentialSeq java.lang.Number]
   1.520 +  [dseq x]
   1.521 +  (fmap (partial * x) dseq))
   1.522 +
   1.523 +(defmethod + [java.lang.Number DifferentialSeq]
   1.524 +  [x dseq]
   1.525 +  (+ (differential-seq x []) dseq))
   1.526 +(defmethod + [DifferentialSeq java.lang.Number]
   1.527 +  [dseq x]
   1.528 +  (+ dseq (differential-seq x [])))
   1.529 +
   1.530 +(defmethod - DifferentialSeq
   1.531 +  [x]
   1.532 +  (fmap - x))
   1.533 +
   1.534 +
   1.535 +;; DERIVATIVES
   1.536 +
   1.537 +	      
   1.538 +
   1.539 +(defn linear-approximator
   1.540 +  "Returns an operator that linearly approximates the given function."
   1.541 +  ([f df|dx]
   1.542 +      (fn [x]
   1.543 +	(let [big-part (big-part x)
   1.544 +	      small-part (small-part x)]
   1.545 +	  ;; f(x+dx) ~= f(x) + f'(x)dx
   1.546 +	  (+ (f big-part)
   1.547 +	     (* (df|dx big-part) small-part)
   1.548 +	  ))))
   1.549 +     
   1.550 +  ([f df|dx df|dy]
   1.551 +     (fn [x y]
   1.552 +       (let [X (big-part x)
   1.553 +	     Y (big-part y)
   1.554 +	     DX (small-part x)
   1.555 +	     DY (small-part y)]
   1.556 +	 (+ (f X Y)
   1.557 +	    (* DX (f df|dx X Y))
   1.558 +	    (* DY (f df|dy X Y)))))))
   1.559 +
   1.560 +
   1.561 +
   1.562 +
   1.563 +
   1.564 +(defn D[f]
   1.565 +  (fn[x] (f (+ x (differential-seq  1 [0] 1 [1] 1 [2])))))
   1.566 +  
   1.567 +(defn d[partials f]
   1.568 +  (fn [x]
   1.569 +    (get 
   1.570 +     (.terms ((D f)x))
   1.571 +     (set partials)
   1.572 +     0
   1.573 +    )))
   1.574 +
   1.575 +(defmethod exp DifferentialSeq [x]
   1.576 +	   ((linear-approximator exp exp) x))
   1.577 +
   1.578 +(defmethod sin DifferentialSeq
   1.579 +  [x]
   1.580 +  ((linear-approximator sin cos) x))
   1.581 +
   1.582 +(defmethod cos DifferentialSeq
   1.583 +  [x]
   1.584 +  ((linear-approximator cos #(- (sin %))) x))
   1.585 +
   1.586 +(defmethod log DifferentialSeq
   1.587 +  [x]
   1.588 +  ((linear-approximator log (fn [x] (/ x)) ) x))
   1.589 +
   1.590 +(defmethod / [DifferentialSeq DifferentialSeq]
   1.591 +  [x y]
   1.592 +  ((linear-approximator /
   1.593 +			(fn [x y] (/ 1 y))
   1.594 +			(fn [x y] (- (/ x (* y y)))))
   1.595 +   x y))
   1.596 +
   1.597 +#+end_src
   1.598 +
   1.599 +
   1.600 +
   1.601 +* Derivatives Revisited
   1.602 +#+begin_src clojure
   1.603 +(in-ns 'sicm.utils)
   1.604 +(use 'clojure.contrib.seq
   1.605 +     'clojure.contrib.generic.arithmetic
   1.606 +     'clojure.contrib.generic.collection
   1.607 +     'clojure.contrib.generic.math-functions)
   1.608 +
   1.609 +(defn replace-in
   1.610 +  "Replaces the nth item in coll with the given item. If n is larger
   1.611 +  than the size of coll, adds n to the end of the collection."
   1.612 +  [coll n item]
   1.613 +  (concat
   1.614 +   (take n coll)
   1.615 +   [item]
   1.616 +   (drop (inc n) coll)))
   1.617 +
   1.618 +(defn euclidean-structure [f partials]
   1.619 +  (fn sd[g v]
   1.620 +    (cond
   1.621 +     (tuple? v)
   1.622 +     (opposite-spin
   1.623 +      v
   1.624 +      (map
   1.625 +       (fn [n]
   1.626 +	 (sd (fn [xn]
   1.627 +	       (g
   1.628 +		(same-spin v
   1.629 +			   (replace-in v n xn))))
   1.630 +	     (nth v n)))
   1.631 +       (range (count v)))))))
   1.632 +   
   1.633 +      
   1.634 +  
   1.635 +
   1.636 +#+end_src
   1.637 +
   1.638 +
   1.639 +* Symbolic Quantities
   1.640 +
   1.641 +#+srcname: symbolic
   1.642 +#+begin_src clojure
   1.643 +(in-ns 'sicm.utils)
   1.644 +(use 'clojure.contrib.generic.arithmetic
   1.645 +     'clojure.contrib.generic.math-functions)
   1.646 +
   1.647 +(deftype Symbolic
   1.648 +  [type
   1.649 +   expression])
   1.650 +
   1.651 +(defn print-expression [s]
   1.652 +  (print (.expression  s)))
   1.653 +
   1.654 +(defn symbolic-number
   1.655 +  [symbol]
   1.656 +  (Symbolic. java.lang.Number symbol))
   1.657 +    
   1.658 +(defn simplify-expression [s]
   1.659 +  (let [e (.expression s)]
   1.660 +    (cond
   1.661 +     (not(list? e)) e
   1.662 +     (= (first e) '+) 
   1.663 +)
   1.664 +
   1.665 +
   1.666 +
   1.667 +(defmethod + [Symbolic Symbolic]
   1.668 +  [x y]
   1.669 +  (Symbolic. (.type x) (list '+ (.expression x) (.expression y))))
   1.670 +
   1.671 +(defmethod + [Symbolic java.lang.Number]
   1.672 +  [s x]
   1.673 +  (if (zero? x)
   1.674 +    s
   1.675 +    (Symbolic. (.type s) (list '+ (.expression s) x))))
   1.676 +
   1.677 +(defmethod sin Symbolic
   1.678 +  [s]
   1.679 +  (Symbolic. (.type s) (list 'sin (.expression s))))
   1.680 +
   1.681 +#+end_src
   1.682 +
   1.683 +* COMMENT To-be-tangled Source
   1.684 +#+begin_src clojure :tangle utils.clj
   1.685 +(ns sicm.utils)
   1.686 +
   1.687 +<<tuples>>
   1.688 +<<matrices>>
   1.689 +<<arith-tuple>>
   1.690 +
   1.691 +<<differential>>
   1.692 +#+end_src
   1.693 +