Mercurial > dylan

annotate sicm/bk/utils.org @ 3:44d3dc936f6a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
moved backup files
author Robert McIntyre <rlm@mit.edu>
date Fri, 28 Oct 2011 00:06:37 -0700
parents b4de894a1e2e
children
rev   line source
rlm@2 1 #+TITLE:Building a Classical Mechanics Library in Clojure
rlm@2 2 #+author: Robert McIntyre & Dylan Holmes
rlm@2 3 #+EMAIL: rlm@mit.edu
rlm@2 4 #+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js"
rlm@2 5 #+STYLE: <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />
rlm@2 6 #+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
rlm@2 7 #+SETUPFILE: ../templates/level-0.org
rlm@2 8 #+INCLUDE: ../templates/level-0.org
rlm@2 9 #+BABEL: :noweb yes :results silent
rlm@2 10
rlm@2 11 * Generic Arithmetic
rlm@2 12
rlm@2 13 #+srcname: generic-arithmetic
rlm@2 14 #+begin_src clojure
rlm@2 15 (ns sicm.utils)
rlm@2 16 (in-ns 'sicm.utils)
rlm@2 17
rlm@2 18
rlm@2 19
rlm@2 20 (defn all-equal? [coll]
rlm@2 21 (if (empty? (rest coll)) true
rlm@2 22 (and (= (first coll) (second coll))
rlm@2 23 (recur (rest coll)))))
rlm@2 24
rlm@2 25
rlm@2 26 (defprotocol Arithmetic
rlm@2 27 (zero [this])
rlm@2 28 (one [this])
rlm@2 29 (negate [this])
rlm@2 30 (invert [this])
rlm@2 31 (conjugate [this]))
rlm@2 32
rlm@2 33 (defn zero? [x]
rlm@2 34 (= x (zero x)))
rlm@2 35 (defn one? [x]
rlm@2 36 (= x (one x)))
rlm@2 37
rlm@2 38
rlm@2 39
rlm@2 40 (extend-protocol Arithmetic
rlm@2 41 java.lang.Number
rlm@2 42 (zero [x] 0)
rlm@2 43 (one [x] 1)
rlm@2 44 (negate [x] (- x))
rlm@2 45 (invert [x] (/ x))
rlm@2 46 )
rlm@2 47
rlm@2 48 (extend-protocol Arithmetic
rlm@2 49 clojure.lang.IFn
rlm@2 50 (one [f] identity)
rlm@2 51 (negate [f] (comp negate f))
rlm@2 52 (invert [f] (comp invert f)))
rlm@2 53
rlm@2 54
rlm@2 55 (extend-protocol Arithmetic
rlm@2 56 clojure.lang.Seqable
rlm@2 57 (zero [this] (map zero this))
rlm@2 58 (one [this] (map one this))
rlm@2 59 (invert [this] (map invert this))
rlm@2 60 (negate [this] (map negate this)))
rlm@2 61
rlm@2 62
rlm@2 63 (defn ordered-like
rlm@2 64 "Create a comparator using the sorted collection as an
rlm@2 65 example. Elements not in the sorted collection are sorted to the
rlm@2 66 end."
rlm@2 67 [sorted-coll]
rlm@2 68 (let [
rlm@2 69 sorted-coll? (set sorted-coll)
rlm@2 70 ascending-pair?
rlm@2 71 (set(reduce concat
rlm@2 72 (map-indexed
rlm@2 73 (fn [n x]
rlm@2 74 (map #(vector x %) (nthnext sorted-coll n)))
rlm@2 75 sorted-coll)))]
rlm@2 76 (fn [x y]
rlm@2 77 (cond
rlm@2 78 (= x y) 0
rlm@2 79 (ascending-pair? [x y]) -1
rlm@2 80 (ascending-pair? [y x]) 1
rlm@2 81 (sorted-coll? x) -1
rlm@2 82 (sorted-coll? y) 1))))
rlm@2 83
rlm@2 84
rlm@2 85
rlm@2 86 (def type-precedence
rlm@2 87 (ordered-like [incanter.Matrix]))
rlm@2 88
rlm@2 89 (defmulti add
rlm@2 90 (fn [x y]
rlm@2 91 (sort type-precedence [(type x)(type y)])))
rlm@2 92
rlm@2 93 (defmulti multiply
rlm@2 94 (fn [x y]
rlm@2 95 (sort type-precedence [(type x) (type y)])))
rlm@2 96
rlm@2 97 (defmethod add [java.lang.Number java.lang.Number] [x y] (+ x y))
rlm@2 98 (defmethod multiply [java.lang.Number java.lang.Number] [x y] (* x y))
rlm@2 99
rlm@2 100 (defmethod multiply [incanter.Matrix java.lang.Integer] [x y]
rlm@2 101 (let [args (sort #(type-precedence (type %1)(type %2)) [x y])
rlm@2 102 matrix (first args)
rlm@2 103 scalar (second args)]
rlm@2 104 (incanter.core/matrix (map (partial map (partial multiply scalar)) matrix))))
rlm@2 105
rlm@2 106 #+end_src
rlm@2 107
rlm@2 108
rlm@2 109 * Useful Data Types
rlm@2 110
rlm@2 111 ** Complex Numbers
rlm@2 112 #+srcname: complex-numbers
rlm@2 113 #+begin_src clojure
rlm@2 114 (in-ns 'sicm.utils)
rlm@2 115
rlm@2 116 (defprotocol Complex
rlm@2 117 (real-part [z])
rlm@2 118 (imaginary-part [z])
rlm@2 119 (magnitude-squared [z])
rlm@2 120 (angle [z])
rlm@2 121 (conjugate [z])
rlm@2 122 (norm [z]))
rlm@2 123
rlm@2 124 (defn complex-rectangular
rlm@2 125 "Define a complex number with the given real and imaginary
rlm@2 126 components."
rlm@2 127 [re im]
rlm@2 128 (reify Complex
rlm@2 129 (real-part [z] re)
rlm@2 130 (imaginary-part [z] im)
rlm@2 131 (magnitude-squared [z] (+ (* re re) (* im im)))
rlm@2 132 (angle [z] (java.lang.Math/atan2 im re))
rlm@2 133 (conjugate [z] (complex-rectangular re (- im)))
rlm@2 134
rlm@2 135 Arithmetic
rlm@2 136 (zero [z] (complex-rectangular 0 0))
rlm@2 137 (one [z] (complex-rectangular 1 0))
rlm@2 138 (negate [z] (complex-rectangular (- re) (- im)))
rlm@2 139 (invert [z] (complex-rectangular
rlm@2 140 (/ re (magnitude-squared z))
rlm@2 141 (/ (- im) (magnitude-squared z))))
rlm@2 142
rlm@2 143 Object
rlm@2 144 (toString [_]
rlm@2 145 (if (and (zero? re) (zero? im)) (str 0)
rlm@2 146 (str
rlm@2 147 (if (not(zero? re))
rlm@2 148 re)
rlm@2 149 (if ((comp not zero?) im)
rlm@2 150 (str
rlm@2 151 (if (neg? im) "-" "+")
rlm@2 152 (if ((comp not one?) (java.lang.Math/abs im))
rlm@2 153 (java.lang.Math/abs im))
rlm@2 154 "i")))))))
rlm@2 155
rlm@2 156 (defn complex-polar
rlm@2 157 "Define a complex number with the given magnitude and angle."
rlm@2 158 [mag ang]
rlm@2 159 (reify Complex
rlm@2 160 (magnitude-squared [z] (* mag mag))
rlm@2 161 (angle [z] angle)
rlm@2 162 (real-part [z] (* mag (java.lang.Math/cos ang)))
rlm@2 163 (imaginary-part [z] (* mag (java.lang.Math/sin ang)))
rlm@2 164 (conjugate [z] (complex-polar mag (- ang)))
rlm@2 165
rlm@2 166 Arithmetic
rlm@2 167 (zero [z] (complex-polar 0 0))
rlm@2 168 (one [z] (complex-polar 1 0))
rlm@2 169 (negate [z] (complex-polar (- mag) ang))
rlm@2 170 (invert [z] (complex-polar (/ mag) (- ang)))
rlm@2 171
rlm@2 172 Object
rlm@2 173 (toString [_] (str mag " * e^(i" ang")"))
rlm@2 174 ))
rlm@2 175
rlm@2 176
rlm@2 177 ;; Numbers are complex quantities
rlm@2 178
rlm@2 179 (extend-protocol Complex
rlm@2 180 java.lang.Number
rlm@2 181 (real-part [x] x)
rlm@2 182 (imaginary-part [x] 0)
rlm@2 183 (magnitude [x] x)
rlm@2 184 (angle [x] 0)
rlm@2 185 (conjugate [x] x))
rlm@2 186
rlm@2 187
rlm@2 188 #+end_src
rlm@2 189
rlm@2 190
rlm@2 191
rlm@2 192 ** Tuples and Tensors
rlm@2 193
rlm@2 194 A tuple is a vector which is spinable\mdash{}it can be either /spin
rlm@2 195 up/ or /spin down/. (Covariant, contravariant; dual vectors)
rlm@2 196 #+srcname: tuples
rlm@2 197 #+begin_src clojure
rlm@2 198 (in-ns 'sicm.utils)
rlm@2 199
rlm@2 200 (defprotocol Spinning
rlm@2 201 (up? [this])
rlm@2 202 (down? [this]))
rlm@2 203
rlm@2 204 (defn spin
rlm@2 205 "Returns the spin of the Spinning s, either :up or :down"
rlm@2 206 [#^Spinning s]
rlm@2 207 (cond (up? s) :up (down? s) :down))
rlm@2 208
rlm@2 209
rlm@2 210 (deftype Tuple
rlm@2 211 [spin coll]
rlm@2 212 clojure.lang.Seqable
rlm@2 213 (seq [this] (seq (.coll this)))
rlm@2 214 clojure.lang.Counted
rlm@2 215 (count [this] (count (.coll this))))
rlm@2 216
rlm@2 217 (extend-type Tuple
rlm@2 218 Spinning
rlm@2 219 (up? [this] (= ::up (.spin this)))
rlm@2 220 (down? [this] (= ::down (.spin this))))
rlm@2 221
rlm@2 222 (defmethod print-method Tuple
rlm@2 223 [o w]
rlm@2 224 (print-simple (str (if (up? o) 'u 'd) (.coll o)) w))
rlm@2 225
rlm@2 226
rlm@2 227
rlm@2 228 (defn up
rlm@2 229 "Create a new up-tuple containing the contents of coll."
rlm@2 230 [coll]
rlm@2 231 (Tuple. ::up coll))
rlm@2 232
rlm@2 233 (defn down
rlm@2 234 "Create a new down-tuple containing the contents of coll."
rlm@2 235 [coll]
rlm@2 236 (Tuple. ::down coll))
rlm@2 237
rlm@2 238
rlm@2 239 #+end_src
rlm@2 240
rlm@2 241 *** Contraction
rlm@2 242 Contraction is a binary operation that you can apply to compatible
rlm@2 243 tuples. Tuples are compatible for contraction if they have the same
rlm@2 244 length and opposite spins, and if the corresponding items in each
rlm@2 245 tuple are both numbers or both compatible tuples.
rlm@2 246
rlm@2 247 #+srcname: tuples-2
rlm@2 248 #+begin_src clojure
rlm@2 249 (in-ns 'sicm.utils)
rlm@2 250
rlm@2 251 (defn numbers?
rlm@2 252 "Returns true if all arguments are numbers, else false."
rlm@2 253 [& xs]
rlm@2 254 (every? number? xs))
rlm@2 255
rlm@2 256 (defn contractible?
rlm@2 257 "Returns true if the tuples a and b are compatible for contraction,
rlm@2 258 else false. Tuples are compatible if they have the same number of
rlm@2 259 components, they have opposite spins, and their elements are
rlm@2 260 pairwise-compatible."
rlm@2 261 [a b]
rlm@2 262 (and
rlm@2 263 (isa? (type a) Tuple)
rlm@2 264 (isa? (type b) Tuple)
rlm@2 265 (= (count a) (count b))
rlm@2 266 (not= (spin a) (spin b))
rlm@2 267
rlm@2 268 (not-any? false?
rlm@2 269 (map #(or
rlm@2 270 (numbers? %1 %2)
rlm@2 271 (contractible? %1 %2))
rlm@2 272 a b))))
rlm@2 273
rlm@2 274
rlm@2 275
rlm@2 276 (defn contract
rlm@2 277 "Contracts two tuples, returning the sum of the
rlm@2 278 products of the corresponding items. Contraction is recursive on
rlm@2 279 nested tuples."
rlm@2 280 [a b]
rlm@2 281 (if (not (contractible? a b))
rlm@2 282 (throw
rlm@2 283 (Exception. "Not compatible for contraction."))
rlm@2 284 (reduce +
rlm@2 285 (map
rlm@2 286 (fn [x y]
rlm@2 287 (if (numbers? x y)
rlm@2 288 (* x y)
rlm@2 289 (contract x y)))
rlm@2 290 a b))))
rlm@2 291
rlm@2 292 #+end_src
rlm@2 293
rlm@2 294 *** Matrices
rlm@2 295 #+srcname: matrices
rlm@2 296 #+begin_src clojure
rlm@2 297 (in-ns 'sicm.utils)
rlm@2 298 (require 'incanter.core) ;; use incanter's fast matrices
rlm@2 299
rlm@2 300 (defprotocol Matrix
rlm@2 301 (rows [matrix])
rlm@2 302 (cols [matrix])
rlm@2 303 (diagonal [matrix])
rlm@2 304 (trace [matrix])
rlm@2 305 (determinant [matrix])
rlm@2 306 (transpose [matrix])
rlm@2 307 (conjugate [matrix])
rlm@2 308 )
rlm@2 309
rlm@2 310 (extend-protocol Matrix
rlm@2 311 incanter.Matrix
rlm@2 312 (rows [rs] (map down (apply map vector (apply map vector rs))))
rlm@2 313 (cols [rs] (map up (apply map vector rs)))
rlm@2 314 (diagonal [matrix] (incanter.core/diag matrix) )
rlm@2 315 (determinant [matrix] (incanter.core/det matrix))
rlm@2 316 (trace [matrix] (incanter.core/trace matrix))
rlm@2 317 (transpose [matrix] (incanter.core/trans matrix))
rlm@2 318 )
rlm@2 319
rlm@2 320 (defn count-rows [matrix]
rlm@2 321 ((comp count rows) matrix))
rlm@2 322
rlm@2 323 (defn count-cols [matrix]
rlm@2 324 ((comp count cols) matrix))
rlm@2 325
rlm@2 326 (defn square? [matrix]
rlm@2 327 (= (count-rows matrix) (count-cols matrix)))
rlm@2 328
rlm@2 329 (defn identity-matrix
rlm@2 330 "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2 331 0s everywhere else."
rlm@2 332 [n]
rlm@2 333 (incanter.core/identity-matrix n))
rlm@2 334
rlm@2 335
rlm@2 336
rlm@2 337
rlm@2 338
rlm@2 339 (defn matrix-by-rows
rlm@2 340 "Define a matrix by giving its rows."
rlm@2 341 [& rows]
rlm@2 342 (if
rlm@2 343 (not (all-equal? (map count rows)))
rlm@2 344 (throw (Exception. "All rows in a matrix must have the same number of elements."))
rlm@2 345 (incanter.core/matrix (vec rows))))
rlm@2 346
rlm@2 347 (defn matrix-by-cols
rlm@2 348 "Define a matrix by giving its columns"
rlm@2 349 [& cols]
rlm@2 350 (if (not (all-equal? (map count cols)))
rlm@2 351 (throw (Exception. "All columns in a matrix must have the same number of elements."))
rlm@2 352 (incanter.core/matrix (vec (apply map vector cols)))))
rlm@2 353
rlm@2 354 (defn identity-matrix
rlm@2 355 "Define a square matrix of size n-by-n with 1s along the diagonal and
rlm@2 356 0s everywhere else."
rlm@2 357 [n]
rlm@2 358 (incanter.core/identity-matrix n))
rlm@2 359
rlm@2 360
rlm@2 361
rlm@2 362 (extend-protocol Arithmetic
rlm@2 363 incanter.Matrix
rlm@2 364 (one [matrix]
rlm@2 365 (if (square? matrix)
rlm@2 366 (identity-matrix (count-rows matrix))
rlm@2 367 (throw (Exception. "Non-square matrices have no multiplicative unit."))))
rlm@2 368 (zero [matrix]
rlm@2 369 (apply matrix-by-rows (map zero (rows matrix))))
rlm@2 370 (negate [matrix]
rlm@2 371 (apply matrix-by-rows (map negate (rows matrix))))
rlm@2 372 (invert [matrix]
rlm@2 373 (incanter.core/solve matrix)))
rlm@2 374
rlm@2 375
rlm@2 376
rlm@2 377 (defmulti coerce-to-matrix
rlm@2 378 "Converts x into a matrix, if possible."
rlm@2 379 type)
rlm@2 380
rlm@2 381 (defmethod coerce-to-matrix incanter.Matrix [x] x)
rlm@2 382 (defmethod coerce-to-matrix Tuple [x]
rlm@2 383 (if (apply numbers? (seq x))
rlm@2 384 (if (up? x)
rlm@2 385 (matrix-by-cols (seq x))
rlm@2 386 (matrix-by-rows (seq x)))
rlm@2 387 (throw (Exception. "Non-numerical tuple cannot be converted into a matrix."))))
rlm@2 388
rlm@2 389
rlm@2 390
rlm@2 391
rlm@2 392
rlm@2 393 ;; (defn matrix-by-cols
rlm@2 394 ;; "Define a matrix by giving its columns."
rlm@2 395 ;; [& cols]
rlm@2 396 ;; (cond
rlm@2 397 ;; (not (all-equal? (map count cols)))
rlm@2 398 ;; (throw (Exception. "All columns in a matrix must have the same number of elements."))
rlm@2 399 ;; :else
rlm@2 400 ;; (reify Matrix
rlm@2 401 ;; (cols [this] (map up cols))
rlm@2 402 ;; (rows [this] (map down (apply map vector cols)))
rlm@2 403 ;; (diagonal [this] (map-indexed (fn [i col] (nth col i) cols)))
rlm@2 404 ;; (trace [this]
rlm@2 405 ;; (if (not= (count-cols this) (count-rows this))
rlm@2 406 ;; (throw (Exception.
rlm@2 407 ;; "Cannot take the trace of a non-square matrix."))
rlm@2 408 ;; (reduce + (diagonal this))))
rlm@2 409
rlm@2 410 ;; (determinant [this]
rlm@2 411 ;; (if (not= (count-cols this) (count-rows this))
rlm@2 412 ;; (throw (Exception.
rlm@2 413 ;; "Cannot take the determinant of a non-square matrix."))
rlm@2 414 ;; (reduce * (map-indexed (fn [i col] (nth col i)) cols))))
rlm@2 415 ;; )))
rlm@2 416
rlm@2 417 (extend-protocol Matrix Tuple
rlm@2 418 (rows [this] (if (down? this)
rlm@2 419 (list this)
rlm@2 420 (map (comp up vector) this)))
rlm@2 421
rlm@2 422 (cols [this] (if (up? this)
rlm@2 423 (list this)
rlm@2 424 (map (comp down vector) this))
rlm@2 425 ))
rlm@2 426
rlm@2 427 (defn matrix-multiply
rlm@2 428 "Returns the matrix resulting from the matrix multiplication of the given arguments."
rlm@2 429 ([A] (coerce-to-matrix A))
rlm@2 430 ([A B] (incanter.core/mmult (coerce-to-matrix A) (coerce-to-matrix B)))
rlm@2 431 ([M1 M2 & Ms] (reduce matrix-multiply (matrix-multiply M1 M2) Ms)))
rlm@2 432
rlm@2 433 #+end_src
rlm@2 434
rlm@2 435
rlm@2 436 ** Power Series
rlm@2 437 #+srcname power-series
rlm@2 438 #+begin_src clojure
rlm@2 439 (in-ns 'sicm.utils)
rlm@2 440 (use 'clojure.contrib.def)
rlm@2 441
rlm@2 442
rlm@2 443
rlm@2 444 (defn series-fn
rlm@2 445 "The function corresponding to the given power series."
rlm@2 446 [series]
rlm@2 447 (fn [x]
rlm@2 448 (reduce +
rlm@2 449 (map-indexed (fn[n x] (* (float (nth series n)) (float(java.lang.Math/pow (float x) n)) ))
rlm@2 450 (range 20)))))
rlm@2 451
rlm@2 452 (deftype PowerSeries
rlm@2 453 [coll]
rlm@2 454 clojure.lang.Seqable
rlm@2 455 (seq [this] (seq (.coll this)))
rlm@2 456
rlm@2 457 clojure.lang.Indexed
rlm@2 458 (nth [this n] (nth (.coll this) n 0))
rlm@2 459 (nth [this n not-found] (nth (.coll this) n not-found))
rlm@2 460
rlm@2 461 ;; clojure.lang.IFn
rlm@2 462 ;; (call [this] (throw(Exception.)))
rlm@2 463 ;; (invoke [this & args] args
rlm@2 464 ;; (let [f
rlm@2 465 ;; )
rlm@2 466 ;; (run [this] (throw(Exception.)))
rlm@2 467 )
rlm@2 468
rlm@2 469 (defn power-series
rlm@2 470 "Returns a power series with the items of the coll as its
rlm@2 471 coefficients. Trailing zeros are added to the end of coll."
rlm@2 472 [coeffs]
rlm@2 473 (PowerSeries. coeffs))
rlm@2 474
rlm@2 475 (defn power-series-indexed
rlm@2 476 "Returns a power series consisting of the result of mapping f to the non-negative integers."
rlm@2 477 [f]
rlm@2 478 (PowerSeries. (map f (range))))
rlm@2 479
rlm@2 480
rlm@2 481 (defn-memo nth-partial-sum
rlm@2 482 ([series n]
rlm@2 483 (if (zero? n) (first series)
rlm@2 484 (+ (nth series n)
rlm@2 485 (nth-partial-sum series (dec n))))))
rlm@2 486
rlm@2 487 (defn partial-sums [series]
rlm@2 488 (lazy-seq (map nth-partial-sum (range))))
rlm@2 489
rlm@2 490
rlm@2 491
rlm@2 492
rlm@2 493 (def cos-series
rlm@2 494 (power-series-indexed
rlm@2 495 (fn[n]
rlm@2 496 (if (odd? n) 0
rlm@2 497 (/
rlm@2 498 (reduce *
rlm@2 499 (reduce * (repeat (/ n 2) -1))
rlm@2 500 (range 1 (inc n)))
rlm@2 501 )))))
rlm@2 502
rlm@2 503 (def sin-series
rlm@2 504 (power-series-indexed
rlm@2 505 (fn[n]
rlm@2 506 (if (even? n) 0
rlm@2 507 (/
rlm@2 508 (reduce *
rlm@2 509 (reduce * (repeat (/ (dec n) 2) -1))
rlm@2 510 (range 1 (inc n)))
rlm@2 511 )))))
rlm@2 512
rlm@2 513 #+end_src
rlm@2 514
rlm@2 515
rlm@2 516 * Basic Utilities
rlm@2 517
rlm@2 518 ** Sequence manipulation
rlm@2 519
rlm@2 520 #+srcname: seq-manipulation
rlm@2 521 #+begin_src clojure
rlm@2 522 (ns sicm.utils)
rlm@2 523
rlm@2 524 (defn do-up
rlm@2 525 "Apply f to each number from low to high, presumably for
rlm@2 526 side-effects."
rlm@2 527 [f low high]
rlm@2 528 (doseq [i (range low high)] (f i)))
rlm@2 529
rlm@2 530 (defn do-down
rlm@2 531 "Apply f to each number from high to low, presumably for
rlm@2 532 side-effects."
rlm@2 533 [f high low]
rlm@2 534 (doseq [i (range high low -1)] (f i)))
rlm@2 535
rlm@2 536
rlm@2 537 (defn all-equal? [coll]
rlm@2 538 (if (empty? (rest coll)) true
rlm@2 539 (and (= (first coll) (second coll))
rlm@2 540 (recur (rest coll))))))
rlm@2 541
rlm@2 542 (defn multiplier
rlm@2 543 "Returns a function that 'multiplies' the members of a collection,
rlm@2 544 returning unit if called on an empty collection."
rlm@2 545 [multiply unit]
rlm@2 546 (fn [coll] ((partial reduce multiply unit) coll)))
rlm@2 547
rlm@2 548 (defn divider
rlm@2 549 "Returns a function that 'divides' the first element of a collection
rlm@2 550 by the 'product' of the rest of the collection."
rlm@2 551 [divide multiply invert unit]
rlm@2 552 (fn [coll]
rlm@2 553 (apply
rlm@2 554 (fn
rlm@2 555 ([] unit)
rlm@2 556 ([x] (invert x))
rlm@2 557 ([x y] (divide x y))
rlm@2 558 ([x y & zs] (divide x (reduce multiply y zs))))
rlm@2 559 coll)))
rlm@2 560
rlm@2 561 (defn left-circular-shift
rlm@2 562 "Remove the first element of coll, adding it to the end of coll."
rlm@2 563 [coll]
rlm@2 564 (concat (rest coll) (take 1 coll)))
rlm@2 565
rlm@2 566 (defn right-circular-shift
rlm@2 567 "Remove the last element of coll, adding it to the front of coll."
rlm@2 568 [coll]
rlm@2 569 (cons (last coll) (butlast coll)))
rlm@2 570 #+end_src
rlm@2 571
rlm@2 572
rlm@2 573
rlm@2 574
rlm@2 575 ** Ranges, Airity and Function Composition
rlm@2 576 #+srcname: arity
rlm@2 577 #+begin_src clojure
rlm@2 578 (in-ns 'sicm.utils)
rlm@2 579 (def infinity Double/POSITIVE_INFINITY)
rlm@2 580 (defn infinite? [x] (Double/isInfinite x))
rlm@2 581 (def finite? (comp not infinite?))
rlm@2 582
rlm@2 583 (defn arity-min
rlm@2 584 "Returns the smallest number of arguments f can take."
rlm@2 585 [f]
rlm@2 586 (apply
rlm@2 587 min
rlm@2 588 (map (comp alength #(.getParameterTypes %))
rlm@2 589 (filter (comp (partial = "invoke") #(.getName %))
rlm@2 590 (.getDeclaredMethods (class f))))))
rlm@2 591
rlm@2 592 (defn arity-max
rlm@2 593 "Returns the largest number of arguments f can take, possibly
rlm@2 594 Infinity."
rlm@2 595 [f]
rlm@2 596 (let [methods (.getDeclaredMethods (class f))]
rlm@2 597 (if (not-any? (partial = "doInvoke") (map #(.getName %) methods))
rlm@2 598 (apply max
rlm@2 599 (map (comp alength #(.getParameterTypes %))
rlm@2 600 (filter (comp (partial = "invoke") #(.getName %)) methods)))
rlm@2 601 infinity)))
rlm@2 602
rlm@2 603
rlm@2 604 (def ^{:arglists '([f])
rlm@2 605 :doc "Returns a two-element list containing the minimum and
rlm@2 606 maximum number of args that f can take."}
rlm@2 607 arity-interval
rlm@2 608 (juxt arity-min arity-max))
rlm@2 609
rlm@2 610
rlm@2 611
rlm@2 612 ;; --- intervals
rlm@2 613
rlm@2 614 (defn intersect
rlm@2 615 "Returns the interval of overlap between interval-1 and interval-2"
rlm@2 616 [interval-1 interval-2]
rlm@2 617 (if (or (empty? interval-1) (empty? interval-2)) []
rlm@2 618 (let [left (max (first interval-1) (first interval-2))
rlm@2 619 right (min (second interval-1) (second interval-2))]
rlm@2 620 (if (> left right) []
rlm@2 621 [left right]))))
rlm@2 622
rlm@2 623 (defn same-endpoints?
rlm@2 624 "Returns true if the left endpoint is the same as the right
rlm@2 625 endpoint."
rlm@2 626 [interval]
rlm@2 627 (= (first interval) (second interval)))
rlm@2 628
rlm@2 629 (defn naturals?
rlm@2 630 "Returns true if the left endpoint is 0 and the right endpoint is
rlm@2 631 infinite."
rlm@2 632 [interval]
rlm@2 633 (and (zero? (first interval))
rlm@2 634 (infinite? (second interval))))
rlm@2 635
rlm@2 636
rlm@2 637 (defn fan-in
rlm@2 638 "Returns a function that pipes its input to each of the gs, then
rlm@2 639 applies f to the list of results. Consequently, f must be able to
rlm@2 640 take a number of arguments equal to the number of gs."
rlm@2 641 [f & gs]
rlm@2 642 (fn [& args]
rlm@2 643 (apply f (apply (apply juxt gs) args))))
rlm@2 644
rlm@2 645 (defn fan-in
rlm@2 646 "Returns a function that pipes its input to each of the gs, then applies f to the list of results. The resulting function takes any number of arguments, but will fail if given arguments that are incompatible with any of the gs."
rlm@2 647 [f & gs]
rlm@2 648 (comp (partial apply f) (apply juxt gs)))
rlm@2 649
rlm@2 650
rlm@2 651
rlm@2 652 (defmacro airty-blah-sad [f n more?]
rlm@2 653 (let [syms (vec (map (comp gensym (partial str "x")) (range n)))
rlm@2 654 optional (gensym "xs")]
rlm@2 655 (if more?
rlm@2 656 `(fn ~(conj syms '& optional)
rlm@2 657 (apply ~f ~@syms ~optional))
rlm@2 658 `(fn ~syms (~f ~@syms)))))
rlm@2 659
rlm@2 660 (defmacro airt-whaa* [f n more?]
rlm@2 661 `(airty-blah-sad ~f ~n ~more?))
rlm@2 662
rlm@2 663
rlm@2 664
rlm@2 665
rlm@2 666 (defn fan-in*
rlm@2 667 "Returns a function that pipes its input to each of the gs, then
rlm@2 668 applies f to the list of results. Unlike fan-in, fan-in* strictly
rlm@2 669 enforces arity: it will fail if the gs do not have compatible
rlm@2 670 arities."
rlm@2 671 [f & gs]
rlm@2 672 (let [arity-in (reduce intersect (map arity-interval gs))
rlm@2 673 left (first arity-in)
rlm@2 674 right (second arity-in)
rlm@2 675 composite (fan-in f gs)
rlm@2 676 ]
rlm@2 677 (cond
rlm@2 678 (empty? arity-in)
rlm@2 679 (throw (Exception. "Cannot compose functions with incompatible arities."))
rlm@2 680
rlm@2 681 (not
rlm@2 682 (or (= left right)
rlm@2 683 (and (finite? left)
rlm@2 684 (= right infinity))))
rlm@2 685
rlm@2 686 (throw (Exception.
rlm@2 687 "Compose can only handle arities of the form [n n] or [n infinity]"))
rlm@2 688 :else
rlm@2 689 (airty-blah-sad composite left (= right infinity)))))
rlm@2 690
rlm@2 691
rlm@2 692
rlm@2 693 (defn compose-n "Compose any number of functions together."
rlm@2 694 ([] identity)
rlm@2 695 ([f] f)
rlm@2 696 ([f & fs]
rlm@2 697 (let [fns (cons f fs)]
rlm@2 698 (compose-bin (reduce fan-in (butlast fs)) (last fs))))
rlm@2 699 )
rlm@2 700
rlm@2 701
rlm@2 702
rlm@2 703
rlm@2 704
rlm@2 705
rlm@2 706 (defn iterated
rlm@2 707 ([f n id] (reduce comp id (repeat n f)))
rlm@2 708 ([f n] (reduce comp identity (repeat n f))))
rlm@2 709
rlm@2 710 (defn iterate-until-stable
rlm@2 711 "Repeatedly applies f to x, returning the first result that is close
rlm@2 712 enough to its predecessor."
rlm@2 713 [f close-enough? x]
rlm@2 714 (second (swank.util/find-first
rlm@2 715 (partial apply close-enough?)
rlm@2 716 (partition 2 1 (iterate f x)))))
rlm@2 717
rlm@2 718 (defn lexical< [x y]
rlm@2 719 (neg? (compare (str x) (str y))))
rlm@2 720
rlm@2 721
rlm@2 722 ;; do-up
rlm@2 723 ;; do-down
rlm@2 724 (def make-pairwise-test comparator)
rlm@2 725 ;;all-equal?
rlm@2 726 (def accumulation multiplier)
rlm@2 727 (def inverse-accumulation divider)
rlm@2 728 ;;left-circular-shift
rlm@2 729 ;;right-circular-shift
rlm@2 730 (def exactly-n? same-endpoints?)
rlm@2 731 (def any-number? naturals?)
rlm@2 732 ;; TODO compose
rlm@2 733 ;; TODO compose-n
rlm@2 734 ;; identity
rlm@2 735 (def compose-2 fan-in)
rlm@2 736 (def compose-bin fan-in*)
rlm@2 737 (def any? (constantly true))
rlm@2 738 (def none? (constantly false))
rlm@2 739 (def constant constantly)
rlm@2 740 (def joint-arity intersect)
rlm@2 741 (def a-reduce reduce)
rlm@2 742 ;; filter
rlm@2 743 (def make-map (partial partial map) )
rlm@2 744 (def bracket juxt)
rlm@2 745 ;; TODO apply-to-all
rlm@2 746 ;; TODO nary-combine
rlm@2 747 ;; TODO binary-combine
rlm@2 748 ;; TODO unary-combine
rlm@2 749 ;; iterated
rlm@2 750 ;; iterate-until-stable
rlm@2 751 (def make-function-of-vector (partial partial map))
rlm@2 752 (def make-function-of-arguments (fn [f] (fn [& args] (f args))))
rlm@2 753 (def alphaless lexical<)
rlm@2 754
rlm@2 755 #+end_src
rlm@2 756
rlm@2 757
rlm@2 758
rlm@2 759
rlm@2 760
rlm@2 761
rlm@2 762
rlm@2 763 :
rlm@2 764
rlm@2 765 * Numerical Methods
rlm@2 766 #+srcname: numerical-methods
rlm@2 767 #+begin_src clojure
rlm@2 768 (in-ns 'sicm.utils)
rlm@2 769 (import java.lang.Math)
rlm@2 770 (use 'clojure.contrib.def)
rlm@2 771
rlm@2 772 ;; ---- USEFUL CONSTANTS
rlm@2 773
rlm@2 774 (defn machine-epsilon
rlm@2 775 "Smallest value representable on your machine, as determined by
rlm@2 776 successively dividing a number in half until consecutive results are
rlm@2 777 indistinguishable."
rlm@2 778 []
rlm@2 779 (ffirst
rlm@2 780 (drop-while
rlm@2 781 (comp not zero? second)
rlm@2 782 (partition 2 1
rlm@2 783 (iterate (partial * 0.5) 1)))))
rlm@2 784
rlm@2 785
rlm@2 786 (def pi (Math/PI))
rlm@2 787 (def two-pi (* 2 pi))
rlm@2 788
rlm@2 789 (def eulers-gamma 0.5772156649015328606065)
rlm@2 790
rlm@2 791 (def phi (/ (inc (Math/sqrt 5)) 2))
rlm@2 792
rlm@2 793 (def ln2 (Math/log 2))
rlm@2 794 (def ln10 (Math/log 10))
rlm@2 795 (def exp10 #(Math/pow 10 %))
rlm@2 796 (def exp2 #(Math/pow 2 %))
rlm@2 797
rlm@2 798
rlm@2 799 ;;
rlm@2 800
rlm@2 801 ;; ---- ANGLES AND TRIGONOMETRY
rlm@2 802
rlm@2 803 (defn angle-restrictor
rlm@2 804 "Returns a function that ensures that angles lie in the specified interval of length two-pi."
rlm@2 805 [max-angle]
rlm@2 806 (let [min-angle (- max-angle two-pi)]
rlm@2 807 (fn [x]
rlm@2 808 (if (and
rlm@2 809 (<= min-angle x)
rlm@2 810 (< x max-angle))
rlm@2 811 x
rlm@2 812 (let [corrected-x (- x (* two-pi (Math/floor (/ x two-pi))))]
rlm@2 813 (if (< corrected-x max-angle)
rlm@2 814 corrected-x
rlm@2 815 (- corrected-x two-pi)))))))
rlm@2 816
rlm@2 817 (defn angle-restrict-pi
rlm@2 818 "Coerces angles to lie in the interval from -pi to pi."
rlm@2 819 [angle]
rlm@2 820 ((angle-restrictor pi) angle))
rlm@2 821
rlm@2 822 (defn angle-restrict-two-pi
rlm@2 823 "Coerces angles to lie in the interval from zero to two-pi"
rlm@2 824 [angle]
rlm@2 825 ((angle-restrictor two-pi) angle))
rlm@2 826
rlm@2 827
rlm@2 828
rlm@2 829
rlm@2 830 (defn invert [x] (/ x))
rlm@2 831 (defn negate [x] (- x))
rlm@2 832
rlm@2 833 (defn exp [x] (Math/exp x))
rlm@2 834
rlm@2 835 (defn sin [x] (Math/sin x))
rlm@2 836 (defn cos [x] (Math/cos x))
rlm@2 837 (defn tan [x] (Math/tan x))
rlm@2 838
rlm@2 839 (def sec (comp invert cos))
rlm@2 840 (def csc (comp invert sin))
rlm@2 841
rlm@2 842 (defn sinh [x] (Math/sinh x))
rlm@2 843 (defn cosh [x] (Math/cosh x))
rlm@2 844 (defn tanh [x] (Math/tanh x))
rlm@2 845
rlm@2 846 (def sech (comp invert cosh))
rlm@2 847 (def csch (comp invert sinh))
rlm@2 848
rlm@2 849
rlm@2 850 ;; ------------
rlm@2 851
rlm@2 852 (defn factorial
rlm@2 853 "Computes the factorial of the nonnegative integer n."
rlm@2 854 [n]
rlm@2 855 (if (neg? n)
rlm@2 856 (throw (Exception. "Cannot compute the factorial of a negative number."))
rlm@2 857 (reduce * 1 (range 1 (inc n)))))
rlm@2 858
rlm@2 859 (defn exact-quotient [n d] (/ n d))
rlm@2 860
rlm@2 861 (defn binomial-coefficient
rlm@2 862 "Computes the number of different ways to choose m elements from n."
rlm@2 863 [n m]
rlm@2 864 (assert (<= 0 m n))
rlm@2 865 (let [difference (- n m)]
rlm@2 866 (exact-quotient
rlm@2 867 (reduce * (range n (max difference m) -1 ))
rlm@2 868 (factorial (min difference m)))))
rlm@2 869
rlm@2 870 (defn-memo stirling-1
rlm@2 871 "Stirling numbers of the first kind: the number of permutations of n
rlm@2 872 elements with exactly m permutation cycles. "
rlm@2 873 [n k]
rlm@2 874 ;(assert (<= 1 k n))
rlm@2 875 (if (zero? n)
rlm@2 876 (if (zero? k) 1 0)
rlm@2 877 (+ (stirling-1 (dec n) (dec k))
rlm@2 878 (* (dec n) (stirling-1 (dec n) k)))))
rlm@2 879
rlm@2 880 (defn-memo stirling-2 ;;count-partitions
rlm@2 881 "Stirling numbers of the second kind: the number of ways to partition a set of n elements into k subsets."
rlm@2 882 [n k]
rlm@2 883 (cond
rlm@2 884 (= k 1) 1
rlm@2 885 (= k n) 1
rlm@2 886 :else (+ (stirling-2 (dec n) (dec k))
rlm@2 887 (* k (stirling-2 (dec n) k)))))
rlm@2 888
rlm@2 889 (defn harmonic-number [n]
rlm@2 890 (/ (stirling-1 (inc n) 2)
rlm@2 891 (factorial n)))
rlm@2 892
rlm@2 893
rlm@2 894 (defn sum
rlm@2 895 [f low high]
rlm@2 896 (reduce + (map f (range low (inc high)))))
rlm@2 897
rlm@2 898 #+end_src
rlm@2 899
rlm@2 900
rlm@2 901
rlm@2 902
rlm@2 903
rlm@2 904
rlm@2 905
rlm@2 906
rlm@2 907
rlm@2 908
rlm@2 909
rlm@2 910
rlm@2 911 * Differentiation
rlm@2 912
rlm@2 913 We compute derivatives by passing special *differential objects* $[x,
rlm@2 914 dx]$ through functions. Roughly speaking, applying a function $f$ to a
rlm@2 915 differential object \([x, dx]\) should produce a new differential
rlm@2 916 object $[f(x),\,Df(x)\cdot dx]$.
rlm@2 917
rlm@2 918 \([x,\,dx]\xrightarrow{\quad f \quad}[f(x),\,Df(x)\cdot dx]\)
rlm@2 919 Notice that you can obtain the derivative of $f$ from this
rlm@2 920 differential object, as it is the coefficient of the $dx$ term. Also,
rlm@2 921 as you apply successive functions using this rule, you get the
rlm@2 922 chain-rule answer you expect:
rlm@2 923
rlm@2 924 \([f(x),\,Df(x)\cdot dx]\xrightarrow{\quad g\quad} [gf(x),\,
rlm@2 925 Dgf(x)\cdot Df(x) \cdot dx ]\)
rlm@2 926
rlm@2 927 In order to generalize to multiple variables and multiple derivatives,
rlm@2 928 we use a *power series of differentials*, a sortred infinite sequence which
rlm@2 929 contains all terms like $dx\cdot dy$, $dx^2\cdot dy$, etc.
rlm@2 930
rlm@2 931
rlm@2 932 #+srcname:differential
rlm@2 933 #+begin_src clojure
rlm@2 934 (in-ns 'sicm.utils)
rlm@2 935 (use 'clojure.contrib.combinatorics)
rlm@2 936 (use 'clojure.contrib.generic.arithmetic)
rlm@2 937
rlm@2 938 (defprotocol DifferentialTerm
rlm@2 939 "Protocol for an infinitesimal quantity."
rlm@2 940 (coefficient [this])
rlm@2 941 (partials [this]))
rlm@2 942
rlm@2 943 (extend-protocol DifferentialTerm
rlm@2 944 java.lang.Number
rlm@2 945 (coefficient [this] this)
rlm@2 946 (partials [this] []))
rlm@2 947
rlm@2 948 (deftype DifferentialSeq
rlm@2 949 [terms]
rlm@2 950 clojure.lang.IPersistentCollection
rlm@2 951 (cons [this x]
rlm@2 952 (throw (Exception. x))
rlm@2 953 (DifferentialSeq. (cons x terms)))
rlm@2 954 ;; (seq [this] (seq terms))
rlm@2 955 (count [this] (count terms))
rlm@2 956 (empty [this] (empty? terms)))
rlm@2 957
rlm@2 958
rlm@2 959
rlm@2 960
rlm@2 961
rlm@2 962
rlm@2 963
rlm@2 964
rlm@2 965
rlm@2 966
rlm@2 967
rlm@2 968 (defn coerce-to-differential-seq [x]
rlm@2 969 (cond
rlm@2 970 (= (type x) DifferentialSeq) x
rlm@2 971 (satisfies? DifferentialTerm x) (DifferentialSeq. x)))
rlm@2 972
rlm@2 973
rlm@2 974 (defn differential-term
rlm@2 975 "Returns a differential term with the given coefficient and
rlm@2 976 partials. Coefficient may be any arithmetic object; partials must
rlm@2 977 be a list of non-negative integers."
rlm@2 978 [coefficient partials]
rlm@2 979 (reify DifferentialTerm
rlm@2 980 (partials [_] (set partials))
rlm@2 981 (coefficient [_] coefficient)))
rlm@2 982
rlm@2 983
rlm@2 984
rlm@2 985 (defn differential-seq*
rlm@2 986 ([coefficient partials]
rlm@2 987 (DifferentialSeq. [(differential-term coefficient partials)]))
rlm@2 988 ([coefficient partials & cps]
rlm@2 989 (if cps
rlm@2 990
rlm@2 991
rlm@2 992
rlm@2 993
rlm@2 994 (defn differential-seq
rlm@2 995 "Constructs a sequence of differential terms from a numerical
rlm@2 996 coefficient and a list of keys for variables. If no coefficient is supplied, uses 1."
rlm@2 997 ([variables] (differential-seq 1 variables))
rlm@2 998 ([coefficient variables & cvs]
rlm@2 999 (if (number? coefficient)
rlm@2 1000 (conj (assoc {} (apply sorted-set variables) coefficient)
rlm@2 1001 (if (empty? cvs)
rlm@2 1002 nil
rlm@2 1003 (apply differential-seq cvs)))
rlm@2 1004 (apply differential-seq 1 coefficient 1 variables cvs)
rlm@2 1005 )))
rlm@2 1006
rlm@2 1007
rlm@2 1008 (defn differential-add
rlm@2 1009 "Add two differential sequences by combining like terms."
rlm@2 1010 [dseq1 dseq2]
rlm@2 1011 (merge-with + dseq1 dseq2))
rlm@2 1012
rlm@2 1013 (defn differential-multiply
rlm@2 1014 "Multiply two differential sequences. The square of any differential variable is zero since differential variables are infinitesimally small."
rlm@2 1015 [dseq1 dseq2]
rlm@2 1016 (reduce
rlm@2 1017 (fn [m [[vars1 coeff1] [vars2 coeff2]]]
rlm@2 1018 (if (empty? (clojure.set/intersection vars1 vars2))
rlm@2 1019 (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))
rlm@2 1020 m))
rlm@2 1021 {}
rlm@2 1022 (clojure.contrib.combinatorics/cartesian-product
rlm@2 1023 dseq1
rlm@2 1024 dseq2)))
rlm@2 1025
rlm@2 1026
rlm@2 1027 (defn big-part
rlm@2 1028 "Returns the part of the differential sequence that is finite,
rlm@2 1029 i.e. not infinitely small."
rlm@2 1030 [dseq]
rlm@2 1031 (let
rlm@2 1032 [keys (sort-by count (keys dseq))
rlm@2 1033 smallest-var (last(last keys))]
rlm@2 1034 (apply hash-map
rlm@2 1035 (reduce concat
rlm@2 1036 (remove (comp smallest-var first) dseq)))))
rlm@2 1037
rlm@2 1038 (defn small-part
rlm@2 1039 "Returns the part of the differential sequence that is
rlm@2 1040 infinitesimal."
rlm@2 1041 [dseq]
rlm@2 1042 (let
rlm@2 1043 [keys (sort-by count (keys dseq))
rlm@2 1044 smallest-var (last(last keys))]
rlm@2 1045 (apply hash-map
rlm@2 1046 (reduce concat
rlm@2 1047 (filter (comp smallest-var first) dseq)))))
rlm@2 1048
rlm@2 1049 (defn big-part
rlm@2 1050 "Returns the 'finite' part of the differential sequence."
rlm@2 1051 [dseq]
rlm@2 1052 (let
rlm@2 1053 [keys (sort-by count (keys dseq))
rlm@2 1054 smallest-var (last (last keys))]
rlm@2 1055 (apply hash-map
rlm@2 1056 (reduce concat
rlm@2 1057 (filter (remove smallest-var first) dseq)
rlm@2 1058 ))))
rlm@2 1059
rlm@2 1060
rlm@2 1061 (defn small-part
rlm@2 1062 "Returns the 'infinitesimal' part of the differential sequence."
rlm@2 1063 [dseq]
rlm@2 1064 (let
rlm@2 1065 [keys (sort-by count (keys dseq))
rlm@2 1066 smallest-var (last (last keys))]
rlm@2 1067
rlm@2 1068 (apply hash-map
rlm@2 1069 (reduce concat
rlm@2 1070 (filter (comp smallest-var first) dseq)
rlm@2 1071 ))))
rlm@2 1072
rlm@2 1073 (defn linear-approximator
rlm@2 1074 "Returns the function that corresponds to unary-fn but which accepts and returns differential objects."
rlm@2 1075 ([unary-f dfdx]
rlm@2 1076 (fn F [dseq]
rlm@2 1077 (differential-add
rlm@2 1078 (unary-f (big-part dseq))
rlm@2 1079 (differential-multiply
rlm@2 1080 (dfdx (big-part dseq))
rlm@2 1081 (small-part dseq))))))
rlm@2 1082
rlm@2 1083
rlm@2 1084
rlm@2 1085
rlm@2 1086
rlm@2 1087
rlm@2 1088 ;; ;; A differential term consists of a numerical coefficient and a
rlm@2 1089 ;; ;; sorted
rlm@2 1090 ;; (defrecord DifferentialTerm [coefficient variables])
rlm@2 1091 ;; (defmethod print-method DifferentialTerm
rlm@2 1092 ;; [o w]
rlm@2 1093 ;; (print-simple
rlm@2 1094 ;; (apply str (.coefficient o)(map (comp (partial str "d") name) (.variables o)))
rlm@2 1095 ;; w))
rlm@2 1096
rlm@2 1097
rlm@2 1098 ;; (defn differential-seq
rlm@2 1099 ;; "Constructs a sequence of differential terms from a numerical
rlm@2 1100 ;; coefficient and a list of keywords for variables. If no coefficient is
rlm@2 1101 ;; supplied, uses 1."
rlm@2 1102 ;; ([variables] (differential-seq 1 variables))
rlm@2 1103 ;; ([coefficient variables]
rlm@2 1104 ;; (list
rlm@2 1105 ;; (DifferentialTerm. coefficient (apply sorted-set variables))))
rlm@2 1106 ;; ([coefficient variables & cvs]
rlm@2 1107 ;; (sort-by
rlm@2 1108 ;; #(vec(.variables %))
rlm@2 1109 ;; (concat (differential-seq coefficient variables) (apply differential-seq cvs)))))
rlm@2 1110
rlm@2 1111 #+end_src
rlm@2 1112
rlm@2 1113
rlm@2 1114
rlm@2 1115
rlm@2 1116
rlm@2 1117
rlm@2 1118
rlm@2 1119
rlm@2 1120
rlm@2 1121
rlm@2 1122
rlm@2 1123
rlm@2 1124
rlm@2 1125 * Symbolic Manipulation
rlm@2 1126
rlm@2 1127 #+srcname:symbolic
rlm@2 1128 #+begin_src clojure
rlm@2 1129 (in-ns 'sicm.utils)
rlm@2 1130
rlm@2 1131 (deftype Symbolic [type expression]
rlm@2 1132 Object
rlm@2 1133 (equals [this that]
rlm@2 1134 (cond
rlm@2 1135 (= (.expression this) (.expression that)) true
rlm@2 1136 :else
rlm@2 1137 (Symbolic.
rlm@2 1138 java.lang.Boolean
rlm@2 1139 (list '= (.expression this) (.expression that)))
rlm@2 1140 )))
rlm@2 1141
rlm@2 1142
rlm@2 1143
rlm@2 1144
rlm@2 1145 (deftype AbstractSet [glyph membership-test])
rlm@2 1146 (defn member? [abstract-set x]
rlm@2 1147 ((.membership-test abstract-set) x))
rlm@2 1148
rlm@2 1149 ;; ------------ Some important AbstractSets
rlm@2 1150
rlm@2 1151
rlm@2 1152 (def Real
rlm@2 1153 (AbstractSet.
rlm@2 1154 'R
rlm@2 1155 (fn[x](number? x))))
rlm@2 1156
rlm@2 1157
rlm@2 1158 ;; ------------ Create new AbstractSets from existing ones
rlm@2 1159
rlm@2 1160 (defn abstract-product
rlm@2 1161 "Gives the cartesian product of abstract sets."
rlm@2 1162 ([sets]
rlm@2 1163 (if (= (count sets) 1) (first sets)
rlm@2 1164 (AbstractSet.
rlm@2 1165 (symbol
rlm@2 1166 (apply str
rlm@2 1167 (interpose 'x (map #(.glyph %) sets))))
rlm@2 1168 (fn [x]
rlm@2 1169 (and
rlm@2 1170 (coll? x)
rlm@2 1171 (= (count sets) (count x))
rlm@2 1172 (reduce #(and %1 %2)
rlm@2 1173 true
rlm@2 1174 (map #(member? %1 %2) sets x)))))))
rlm@2 1175 ([abstract-set n]
rlm@2 1176 (abstract-product (repeat n abstract-set))))
rlm@2 1177
rlm@2 1178
rlm@2 1179
rlm@2 1180 (defn abstract-up
rlm@2 1181 "Returns the abstract set of all up-tuples whose items belong to the
rlm@2 1182 corresponding abstract sets in coll."
rlm@2 1183 ([coll]
rlm@2 1184 (AbstractSet.
rlm@2 1185 (symbol (str "u["
rlm@2 1186 (apply str
rlm@2 1187 (interpose " "
rlm@2 1188 (map #(.glyph %) coll)))
rlm@2 1189 "]"))
rlm@2 1190 (fn [x]
rlm@2 1191 (and
rlm@2 1192 (satisfies? Spinning x)
rlm@2 1193 (up? x)
rlm@2 1194 (= (count coll) (count x))
rlm@2 1195 (reduce
rlm@2 1196 #(and %1 %2)
rlm@2 1197 true
rlm@2 1198 (map #(member? %1 %2) coll x))))))
rlm@2 1199 ([abstract-set n]
rlm@2 1200 (abstract-up (repeat n abstract-set))))
rlm@2 1201
rlm@2 1202
rlm@2 1203 (defn abstract-down
rlm@2 1204 "Returns the abstract set of all down-tuples whose items belong to the
rlm@2 1205 corresponding abstract sets in coll."
rlm@2 1206 ([coll]
rlm@2 1207 (AbstractSet.
rlm@2 1208 (symbol (str "d["
rlm@2 1209 (apply str
rlm@2 1210 (interpose " "
rlm@2 1211 (map #(.glyph %) coll)))
rlm@2 1212 "]"))
rlm@2 1213 (fn [x]
rlm@2 1214 (and
rlm@2 1215 (satisfies? Spinning x)
rlm@2 1216 (down? x)
rlm@2 1217 (= (count coll) (count x))
rlm@2 1218 (reduce
rlm@2 1219 #(and %1 %2)
rlm@2 1220 true
rlm@2 1221 (map #(member? %1 %2) coll x))))))
rlm@2 1222 ([abstract-set n]
rlm@2 1223 (abstract-down (repeat n abstract-set))))
rlm@2 1224
rlm@2 1225
rlm@2 1226
rlm@2 1227
rlm@2 1228
rlm@2 1229 ;-------ABSTRACT FUNCTIONS
rlm@2 1230 (defrecord AbstractFn
rlm@2 1231 [#^AbstractSet domain #^AbstractSet codomain])
rlm@2 1232
rlm@2 1233
rlm@2 1234 (defmethod print-method AbstractFn
rlm@2 1235 [o w]
rlm@2 1236 (print-simple
rlm@2 1237 (str
rlm@2 1238 "f:"
rlm@2 1239 (.glyph (:domain o))
rlm@2 1240 "-->"
rlm@2 1241 (.glyph (:codomain o))) w))
rlm@2 1242 #+end_src
rlm@2 1243
rlm@2 1244
rlm@2 1245 * COMMENT
rlm@2 1246 #+begin_src clojure :tangle utils.clj
rlm@2 1247 (ns sicm.utils)
rlm@2 1248
rlm@2 1249 ;***** GENERIC ARITHMETIC
rlm@2 1250 <<generic-arithmetic>>
rlm@2 1251
rlm@2 1252
rlm@2 1253 ;***** TUPLES AND MATRICES
rlm@2 1254 <<tuples>>
rlm@2 1255 <<tuples-2>>
rlm@2 1256 ;***** MATRICES
rlm@2 1257 <<matrices>>
rlm@2 1258
rlm@2 1259 #+end_src
rlm@2 1260
rlm@2 1261
rlm@2 1262