Mercurial > dylan

annotate sicm/utils.clj @ 2:b4de894a1e2e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
initial import
author Robert McIntyre <rlm@mit.edu>
date Fri, 28 Oct 2011 00:03:05 -0700
parents
children
rev   line source
rlm@2 1
rlm@2 2 (ns sicm.utils)
rlm@2 3
rlm@2 4 (in-ns 'sicm.utils)
rlm@2 5
rlm@2 6 ;; Let some objects have spin
rlm@2 7
rlm@2 8 (defprotocol Spinning
rlm@2 9 (up? [this])
rlm@2 10 (down? [this]))
rlm@2 11
rlm@2 12 (defn spin
rlm@2 13 "Returns the spin of the Spinning s, either :up or :down"
rlm@2 14 [#^Spinning s]
rlm@2 15 (cond (up? s) :up (down? s) :down))
rlm@2 16
rlm@2 17
rlm@2 18 ;; DEFINITION: A tuple is a sequence with spin
rlm@2 19
rlm@2 20 (deftype Tuple
rlm@2 21 [spin coll]
rlm@2 22
rlm@2 23 clojure.lang.Seqable
rlm@2 24 (seq [this] (seq (.coll this)))
rlm@2 25
rlm@2 26 clojure.lang.Counted
rlm@2 27 (count [this] (count (.coll this)))
rlm@2 28
rlm@2 29 Spinning
rlm@2 30 (up? [this] (= ::up (.spin this)))
rlm@2 31 (down? [this] (= ::down (.spin this))))
rlm@2 32
rlm@2 33 (defmethod print-method Tuple
rlm@2 34 [o w]
rlm@2 35 (print-simple
rlm@2 36 (if (up? o)
rlm@2 37 (str "u" (.coll o))
rlm@2 38 (str "d" (vec(.coll o))))
rlm@2 39 w))
rlm@2 40
rlm@2 41 (def tuple? #(= (type %) Tuple))
rlm@2 42
rlm@2 43 ;; CONSTRUCTORS
rlm@2 44
rlm@2 45 (defn up
rlm@2 46 "Create a new up-tuple containing the contents of coll."
rlm@2 47 [coll]
rlm@2 48 (Tuple. ::up coll))
rlm@2 49
rlm@2 50 (defn down
rlm@2 51 "Create a new down-tuple containing the contents of coll."
rlm@2 52 [coll]
rlm@2 53 (Tuple. ::down coll))
rlm@2 54
rlm@2 55 (defn same-spin
rlm@2 56 "Creates a tuple which has the same spin as tuple and which contains
rlm@2 57 the contents of coll."
rlm@2 58 [tuple coll]
rlm@2 59 (if (up? tuple)
rlm@2 60 (up coll)
rlm@2 61 (down coll)))
rlm@2 62
rlm@2 63 (defn opposite-spin
rlm@2 64 "Create a tuple which has opposite spin to tuple and which contains
rlm@2 65 the contents of coll."
rlm@2 66 [tuple coll]
rlm@2 67 (if (up? tuple)
rlm@2 68 (down coll)
rlm@2 69 (up coll)))
rlm@2 70 (in-ns 'sicm.utils)
rlm@2 71 (require 'incanter.core) ;; use incanter's fast matrices
rlm@2 72
rlm@2 73
rlm@2 74 (defn all-equal? [coll]
rlm@2 75 (if (empty? (rest coll)) true
rlm@2 76 (and (= (first coll) (second coll))
rlm@2 77 (recur (rest coll)))))
rlm@2 78
rlm@2 79
rlm@2 80 (defprotocol Matrix
rlm@2 81 (rows [matrix])
rlm@2 82 (cols [matrix])
rlm@2 83 (diagonal [matrix])
rlm@2 84 (trace [matrix])
rlm@2 85 (determinant [matrix])
rlm@2 86 (transpose [matrix])
rlm@2 87 (conjugate [matrix])
rlm@2 88 )
rlm@2 89
rlm@2 90 (extend-protocol Matrix incanter.Matrix
rlm@2 91 (rows [rs] (map down (apply map vector (apply map vector rs))))
rlm@2 92 (cols [rs] (map up (apply map vector rs)))
rlm@2 93 (diagonal [matrix] (incanter.core/diag matrix) )
rlm@2 94 (determinant [matrix] (incanter.core/det matrix))
rlm@2 95 (trace [matrix] (incanter.core/trace matrix))
rlm@2 96 (transpose [matrix] (incanter.core/trans matrix)))
rlm@2 97
rlm@2 98 (defn count-rows [matrix]
rlm@2 99 ((comp count rows) matrix))
rlm@2 100
rlm@2 101 (defn count-cols [matrix]
rlm@2 102 ((comp count cols) matrix))
rlm@2 103
rlm@2 104 (defn square? [matrix]
rlm@2 105 (= (count-rows matrix) (count-cols matrix)))
rlm@2 106
rlm@2 107 (defn identity-matrix
rlm@2 108 "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2 109 0s everywhere else."
rlm@2 110 [n]
rlm@2 111 (incanter.core/identity-matrix n))
rlm@2 112
rlm@2 113
rlm@2 114 (defn matrix-by-rows
rlm@2 115 "Define a matrix by giving its rows."
rlm@2 116 [& rows]
rlm@2 117 (if
rlm@2 118 (not (all-equal? (map count rows)))
rlm@2 119 (throw (Exception. "All rows in a matrix must have the same number of elements."))
rlm@2 120 (incanter.core/matrix (vec rows))))
rlm@2 121
rlm@2 122 (defn matrix-by-cols
rlm@2 123 "Define a matrix by giving its columns"
rlm@2 124 [& cols]
rlm@2 125 (if (not (all-equal? (map count cols)))
rlm@2 126 (throw (Exception. "All columns in a matrix must have the same number of elements."))
rlm@2 127 (incanter.core/matrix (vec (apply map vector cols)))))
rlm@2 128
rlm@2 129 (defn identity-matrix
rlm@2 130 "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2 131 0s everywhere else."
rlm@2 132 [n]
rlm@2 133 (incanter.core/identity-matrix n))
rlm@2 134
rlm@2 135 (in-ns 'sicm.utils)
rlm@2 136 (use 'clojure.contrib.generic.arithmetic
rlm@2 137 'clojure.contrib.generic.collection
rlm@2 138 'clojure.contrib.generic.functor
rlm@2 139 'clojure.contrib.generic.math-functions)
rlm@2 140
rlm@2 141 (defn numbers?
rlm@2 142 "Returns true if all arguments are numbers, else false."
rlm@2 143 [& xs]
rlm@2 144 (every? number? xs))
rlm@2 145
rlm@2 146 (defn tuple-surgery
rlm@2 147 "Applies the function f to the items of tuple and the additional
rlm@2 148 arguments, if any. Returns a Tuple of the same type as tuple."
rlm@2 149 [tuple f & xs]
rlm@2 150 ((if (up? tuple) up down)
rlm@2 151 (apply f (seq tuple) xs)))
rlm@2 152
rlm@2 153
rlm@2 154
rlm@2 155 ;;; CONTRACTION collapses two compatible tuples into a number.
rlm@2 156
rlm@2 157 (defn contractible?
rlm@2 158 "Returns true if the tuples a and b are compatible for contraction,
rlm@2 159 else false. Tuples are compatible if they have the same number of
rlm@2 160 components, they have opposite spins, and their elements are
rlm@2 161 pairwise-compatible."
rlm@2 162 [a b]
rlm@2 163 (and
rlm@2 164 (isa? (type a) Tuple)
rlm@2 165 (isa? (type b) Tuple)
rlm@2 166 (= (count a) (count b))
rlm@2 167 (not= (spin a) (spin b))
rlm@2 168
rlm@2 169 (not-any? false?
rlm@2 170 (map #(or
rlm@2 171 (numbers? %1 %2)
rlm@2 172 (contractible? %1 %2))
rlm@2 173 a b))))
rlm@2 174
rlm@2 175 (defn contract
rlm@2 176 "Contracts two tuples, returning the sum of the
rlm@2 177 products of the corresponding items. Contraction is recursive on
rlm@2 178 nested tuples."
rlm@2 179 [a b]
rlm@2 180 (if (not (contractible? a b))
rlm@2 181 (throw
rlm@2 182 (Exception. "Not compatible for contraction."))
rlm@2 183 (reduce +
rlm@2 184 (map
rlm@2 185 (fn [x y]
rlm@2 186 (if (numbers? x y)
rlm@2 187 (* x y)
rlm@2 188 (contract x y)))
rlm@2 189 a b))))
rlm@2 190
rlm@2 191
rlm@2 192
rlm@2 193
rlm@2 194
rlm@2 195 (defmethod conj Tuple
rlm@2 196 [tuple & xs]
rlm@2 197 (tuple-surgery tuple #(apply conj % xs)))
rlm@2 198
rlm@2 199 (defmethod fmap Tuple
rlm@2 200 [f tuple]
rlm@2 201 (tuple-surgery tuple (partial map f)))
rlm@2 202
rlm@2 203
rlm@2 204
rlm@2 205 ;; TODO: define Scalar, and add it to the hierarchy above Number and Complex
rlm@2 206
rlm@2 207
rlm@2 208 (defmethod * [Tuple Tuple] ; tuple*tuple
rlm@2 209 [a b]
rlm@2 210 (if (contractible? a b)
rlm@2 211 (contract a b)
rlm@2 212 (map (partial * a) b)))
rlm@2 213
rlm@2 214
rlm@2 215 (defmethod * [java.lang.Number Tuple] ;; scalar * tuple
rlm@2 216 [a x] (fmap (partial * a) x))
rlm@2 217
rlm@2 218 (defmethod * [Tuple java.lang.Number]
rlm@2 219 [x a] (* a x))
rlm@2 220
rlm@2 221 (defmethod * [java.lang.Number incanter.Matrix] ;; scalar * matrix
rlm@2 222 [x M] (incanter.core/mult x M))
rlm@2 223
rlm@2 224 (defmethod * [incanter.Matrix java.lang.Number]
rlm@2 225 [M x] (* x M))
rlm@2 226
rlm@2 227 (defmethod * [incanter.Matrix incanter.Matrix] ;; matrix * matrix
rlm@2 228 [M1 M2]
rlm@2 229 (incanter.core/mmult M1 M2))
rlm@2 230
rlm@2 231 (defmethod * [incanter.Matrix Tuple] ;; matrix * tuple
rlm@2 232 [M v]
rlm@2 233 (if (and (apply numbers? v) (up? v))
rlm@2 234 (* M (matrix-by-cols v))
rlm@2 235 (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
rlm@2 236 ))
rlm@2 237
rlm@2 238 (defmethod * [Tuple incanter.Matrix] ;; tuple * Matrix
rlm@2 239 [v M]
rlm@2 240 (if (and (apply numbers? v) (down? v))
rlm@2 241 (* (matrix-by-rows v) M)
rlm@2 242 (throw (Exception. "Currently, you can only multiply a matrix by a tuple of *numbers*"))
rlm@2 243 ))
rlm@2 244
rlm@2 245
rlm@2 246 (defmethod exp incanter.Matrix
rlm@2 247 [M]
rlm@2 248 (incanter.core/exp M))
rlm@2 249
rlm@2 250
rlm@2 251 (in-ns 'sicm.utils)
rlm@2 252 (use 'clojure.contrib.seq
rlm@2 253 'clojure.contrib.generic.arithmetic
rlm@2 254 'clojure.contrib.generic.collection
rlm@2 255 'clojure.contrib.generic.math-functions)
rlm@2 256
rlm@2 257 ;;∂
rlm@2 258
rlm@2 259 ;; DEFINITION : Differential Term
rlm@2 260
rlm@2 261 ;; A quantity with infinitesimal components, e.g. x, dxdy, 4dydz. The
rlm@2 262 ;; coefficient of the quantity is returned by the 'coefficient' method,
rlm@2 263 ;; while the sequence of differential parameters is returned by the
rlm@2 264 ;; method 'partials'.
rlm@2 265
rlm@2 266 ;; Instead of using (potentially ambiguous) letters to denote
rlm@2 267 ;; differential parameters (dx,dy,dz), we use integers. So, dxdz becomes [0 2].
rlm@2 268
rlm@2 269 ;; The coefficient can be any arithmetic object; the
rlm@2 270 ;; partials must be a nonrepeating sorted sequence of nonnegative
rlm@2 271 ;; integers.
rlm@2 272
rlm@2 273 (deftype DifferentialTerm [coefficient partials])
rlm@2 274
rlm@2 275 (defn differential-term
rlm@2 276 "Make a differential term from a coefficient and list of partials."
rlm@2 277 [coefficient partials]
rlm@2 278 (if (and (coll? partials) (every? #(and (integer? %) (not(neg? %))) partials))
rlm@2 279 (DifferentialTerm. coefficient (set partials))
rlm@2 280 (throw (java.lang.IllegalArgumentException. "Partials must be a collection of integers."))))
rlm@2 281
rlm@2 282
rlm@2 283 ;; DEFINITION : Differential Sequence
rlm@2 284 ;; A differential sequence is a sequence of differential terms, all with different partials.
rlm@2 285 ;; Internally, it is a map from the partials of each term to their coefficients.
rlm@2 286
rlm@2 287 (deftype DifferentialSeq
rlm@2 288 [terms]
rlm@2 289 ;;clojure.lang.IPersistentMap
rlm@2 290 clojure.lang.Associative
rlm@2 291 (assoc [this key val]
rlm@2 292 (DifferentialSeq.
rlm@2 293 (cons (differential-term val key) terms)))
rlm@2 294 (cons [this x]
rlm@2 295 (DifferentialSeq. (cons x terms)))
rlm@2 296 (containsKey [this key]
rlm@2 297 (not(nil? (find-first #(= (.partials %) key) terms))))
rlm@2 298 (count [this] (count (.terms this)))
rlm@2 299 (empty [this] (DifferentialSeq. []))
rlm@2 300 (entryAt [this key]
rlm@2 301 ((juxt #(.partials %) #(.coefficient %))
rlm@2 302 (find-first #(= (.partials %) key) terms)))
rlm@2 303 (seq [this] (seq (.terms this))))
rlm@2 304
rlm@2 305 (def differential? #(= (type %) DifferentialSeq))
rlm@2 306
rlm@2 307 (defn zeroth-order?
rlm@2 308 "Returns true if the differential sequence has at most a constant term."
rlm@2 309 [dseq]
rlm@2 310 (and
rlm@2 311 (differential? dseq)
rlm@2 312 (every?
rlm@2 313 #(= #{} %)
rlm@2 314 (keys (.terms dseq)))))
rlm@2 315
rlm@2 316 (defmethod fmap DifferentialSeq
rlm@2 317 [f dseq]
rlm@2 318 (DifferentialSeq.
rlm@2 319 (fmap f (.terms dseq))))
rlm@2 320
rlm@2 321
rlm@2 322
rlm@2 323
rlm@2 324 ;; BUILDING DIFFERENTIAL OBJECTS
rlm@2 325
rlm@2 326 (defn differential-seq
rlm@2 327 "Define a differential sequence by specifying an alternating
rlm@2 328 sequence of coefficients and lists of partials."
rlm@2 329 ([coefficient partials]
rlm@2 330 (DifferentialSeq. {(set partials) coefficient}))
rlm@2 331 ([coefficient partials & cps]
rlm@2 332 (if (odd? (count cps))
rlm@2 333 (throw (Exception. "differential-seq requires an even number of terms."))
rlm@2 334 (DifferentialSeq.
rlm@2 335 (reduce
rlm@2 336 #(assoc %1 (set (second %2)) (first %2))
rlm@2 337 {(set partials) coefficient}
rlm@2 338 (partition 2 cps))))))
rlm@2 339
rlm@2 340
rlm@2 341
rlm@2 342 (defn big-part
rlm@2 343 "Returns the part of the differential sequence that is finite,
rlm@2 344 i.e. not infinitely small. If the sequence is zeroth-order, returns
rlm@2 345 the coefficient of the zeroth-order term instead. "
rlm@2 346 [dseq]
rlm@2 347 (if (zeroth-order? dseq) (get (.terms dseq) #{})
rlm@2 348 (let [m (.terms dseq)
rlm@2 349 keys (sort-by count (keys m))
rlm@2 350 smallest-var (last (last keys))]
rlm@2 351 (DifferentialSeq.
rlm@2 352 (reduce
rlm@2 353 #(assoc %1 (first %2) (second %2))
rlm@2 354 {}
rlm@2 355 (remove #((first %) smallest-var) m))))))
rlm@2 356
rlm@2 357
rlm@2 358 (defn small-part
rlm@2 359 "Returns the part of the differential sequence that infinitely
rlm@2 360 small. If the sequence is zeroth-order, returns zero."
rlm@2 361 [dseq]
rlm@2 362 (if (zeroth-order? dseq) 0
rlm@2 363 (let [m (.terms dseq)
rlm@2 364 keys (sort-by count (keys m))
rlm@2 365 smallest-var (last (last keys))]
rlm@2 366 (DifferentialSeq.
rlm@2 367 (reduce
rlm@2 368 #(assoc %1 (first %2) (second %2)) {}
rlm@2 369 (filter #((first %) smallest-var) m))))))
rlm@2 370
rlm@2 371
rlm@2 372
rlm@2 373 (defn cartesian-product [set1 set2]
rlm@2 374 (reduce concat
rlm@2 375 (for [x set1]
rlm@2 376 (for [y set2]
rlm@2 377 [x y]))))
rlm@2 378
rlm@2 379 (defn nth-subset [n]
rlm@2 380 (if (zero? n) []
rlm@2 381 (let [lg2 #(/ (log %) (log 2))
rlm@2 382 k (int(java.lang.Math/floor (lg2 n)))
rlm@2 383 ]
rlm@2 384 (cons k
rlm@2 385 (nth-subset (- n (pow 2 k)))))))
rlm@2 386
rlm@2 387 (def all-partials
rlm@2 388 (lazy-seq (map nth-subset (range))))
rlm@2 389
rlm@2 390
rlm@2 391 (defn differential-multiply
rlm@2 392 "Multiply two differential sequences. The square of any differential
rlm@2 393 variable is zero since differential variables are infinitesimally
rlm@2 394 small."
rlm@2 395 [dseq1 dseq2]
rlm@2 396 (DifferentialSeq.
rlm@2 397 (reduce
rlm@2 398 (fn [m [[vars1 coeff1] [vars2 coeff2]]]
rlm@2 399 (if (not (empty? (clojure.set/intersection vars1 vars2)))
rlm@2 400 m
rlm@2 401 (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))))
rlm@2 402 {}
rlm@2 403 (cartesian-product (.terms dseq1) (.terms dseq2)))))
rlm@2 404
rlm@2 405
rlm@2 406
rlm@2 407 (defmethod * [DifferentialSeq DifferentialSeq]
rlm@2 408 [dseq1 dseq2]
rlm@2 409 (differential-multiply dseq1 dseq2))
rlm@2 410
rlm@2 411 (defmethod + [DifferentialSeq DifferentialSeq]
rlm@2 412 [dseq1 dseq2]
rlm@2 413 (DifferentialSeq.
rlm@2 414 (merge-with + (.terms dseq1) (.terms dseq2))))
rlm@2 415
rlm@2 416 (defmethod * [java.lang.Number DifferentialSeq]
rlm@2 417 [x dseq]
rlm@2 418 (fmap (partial * x) dseq))
rlm@2 419
rlm@2 420 (defmethod * [DifferentialSeq java.lang.Number]
rlm@2 421 [dseq x]
rlm@2 422 (fmap (partial * x) dseq))
rlm@2 423
rlm@2 424 (defmethod + [java.lang.Number DifferentialSeq]
rlm@2 425 [x dseq]
rlm@2 426 (+ (differential-seq x []) dseq))
rlm@2 427 (defmethod + [DifferentialSeq java.lang.Number]
rlm@2 428 [dseq x]
rlm@2 429 (+ dseq (differential-seq x [])))
rlm@2 430
rlm@2 431 (defmethod - DifferentialSeq
rlm@2 432 [x]
rlm@2 433 (fmap - x))
rlm@2 434
rlm@2 435
rlm@2 436 ;; DERIVATIVES
rlm@2 437
rlm@2 438
rlm@2 439
rlm@2 440 (defn linear-approximator
rlm@2 441 "Returns an operator that linearly approximates the given function."
rlm@2 442 ([f df|dx]
rlm@2 443 (fn [x]
rlm@2 444 (let [big-part (big-part x)
rlm@2 445 small-part (small-part x)]
rlm@2 446 ;; f(x+dx) ~= f(x) + f'(x)dx
rlm@2 447 (+ (f big-part)
rlm@2 448 (* (df|dx big-part) small-part)
rlm@2 449 ))))
rlm@2 450
rlm@2 451 ([f df|dx df|dy]
rlm@2 452 (fn [x y]
rlm@2 453 (let [X (big-part x)
rlm@2 454 Y (big-part y)
rlm@2 455 DX (small-part x)
rlm@2 456 DY (small-part y)]
rlm@2 457 (+ (f X Y)
rlm@2 458 (* DX (f df|dx X Y))
rlm@2 459 (* DY (f df|dy X Y)))))))
rlm@2 460
rlm@2 461
rlm@2 462
rlm@2 463
rlm@2 464
rlm@2 465 (defn D[f]
rlm@2 466 (fn[x] (f (+ x (differential-seq 1 [0] 1 [1] 1 [2])))))
rlm@2 467
rlm@2 468 (defn d[partials f]
rlm@2 469 (fn [x]
rlm@2 470 (get
rlm@2 471 (.terms ((D f)x))
rlm@2 472 (set partials)
rlm@2 473 0
rlm@2 474 )))
rlm@2 475
rlm@2 476 (defmethod exp DifferentialSeq [x]
rlm@2 477 ((linear-approximator exp exp) x))
rlm@2 478
rlm@2 479 (defmethod sin DifferentialSeq
rlm@2 480 [x]
rlm@2 481 ((linear-approximator sin cos) x))
rlm@2 482
rlm@2 483 (defmethod cos DifferentialSeq
rlm@2 484 [x]
rlm@2 485 ((linear-approximator cos #(- (sin %))) x))
rlm@2 486
rlm@2 487 (defmethod log DifferentialSeq
rlm@2 488 [x]
rlm@2 489 ((linear-approximator log (fn [x] (/ x)) ) x))
rlm@2 490
rlm@2 491 (defmethod / [DifferentialSeq DifferentialSeq]
rlm@2 492 [x y]
rlm@2 493 ((linear-approximator /
rlm@2 494 (fn [x y] (/ 1 y))
rlm@2 495 (fn [x y] (- (/ x (* y y)))))
rlm@2 496 x y))