dylan: sicm/bk/utils.org comparison

comparison sicm/bk/utils.org @ 2:b4de894a1e2e

initial import

author	Robert McIntyre <rlm@mit.edu>
date	Fri, 28 Oct 2011 00:03:05 -0700
parents
children

comparison

equal deleted inserted replaced

-:8d8278e09888
+:b4de894a1e2e
+#+TITLE:Building a Classical Mechanics Library in Clojure
+#+author: Robert McIntyre & Dylan Holmes
+#+EMAIL:     rlm@mit.edu
+#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js"
+#+STYLE: <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />
+#+OPTIONS:   H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
+#+SETUPFILE: ../templates/level-0.org
+#+INCLUDE: ../templates/level-0.org
+#+BABEL: :noweb yes :results silent
+* Generic Arithmetic
+#+srcname: generic-arithmetic
+#+begin_src clojure
+(ns sicm.utils)
+(in-ns 'sicm.utils)
+(defn all-equal? [coll]
+(if (empty? (rest coll)) true
+(and (= (first coll) (second coll))
+	   (recur (rest coll)))))
+(defprotocol Arithmetic
+(zero [this])
+(one [this])
+(negate [this])
+(invert [this])
+(conjugate [this]))
+(defn zero? [x]
+(= x (zero x)))
+(defn one? [x]
+(= x (one x)))
+(extend-protocol Arithmetic
+java.lang.Number
+(zero [x] 0)
+(one [x] 1)
+(negate [x] (- x))
+(invert [x] (/ x))
+)
+(extend-protocol Arithmetic
+clojure.lang.IFn
+(one [f] identity)
+(negate [f] (comp negate f))
+(invert [f] (comp invert f)))
+(extend-protocol Arithmetic
+clojure.lang.Seqable
+(zero [this] (map zero this))
+(one [this] (map one this))
+(invert [this] (map invert this))
+(negate [this] (map negate this)))
+(defn ordered-like
+"Create a comparator using the sorted collection as an
+example. Elements not in the sorted collection are sorted to the
+end."
+[sorted-coll]
+(let [
+	sorted-coll? (set sorted-coll)
+	ascending-pair?
+	(set(reduce concat
+		    (map-indexed
+		     (fn [n x]
+		       (map #(vector x %) (nthnext sorted-coll n)))
+		     sorted-coll)))]
+(fn [x y]
+(cond
+(= x y) 0
+(ascending-pair? [x y]) -1
+(ascending-pair? [y x]) 1
+(sorted-coll? x) -1
+(sorted-coll? y) 1))))
+(def type-precedence
+(ordered-like [incanter.Matrix]))
+(defmulti add
+(fn [x y]
+(sort type-precedence [(type x)(type y)])))
+(defmulti multiply
+(fn [x y]
+(sort type-precedence [(type x) (type y)])))
+(defmethod add [java.lang.Number java.lang.Number] [x y] (+ x y))
+(defmethod multiply [java.lang.Number java.lang.Number] [x y] (* x y))
+(defmethod multiply [incanter.Matrix java.lang.Integer] [x y]
+	   (let [args (sort #(type-precedence (type %1)(type %2)) [x y])
+		 matrix (first args)
+		 scalar (second args)]
+	     (incanter.core/matrix (map (partial map (partial multiply scalar)) matrix))))
+#+end_src
+* Useful Data Types
+** Complex Numbers
+#+srcname: complex-numbers
+#+begin_src clojure
+(in-ns 'sicm.utils)
+(defprotocol Complex
+(real-part [z])
+(imaginary-part [z])
+(magnitude-squared [z])
+(angle [z])
+(conjugate [z])
+(norm [z]))
+(defn complex-rectangular
+"Define a complex number with the given real and imaginary
+components."
+[re im]
+(reify Complex
+	 (real-part [z] re)
+	 (imaginary-part [z] im)
+	 (magnitude-squared [z] (+ (* re re) (* im im)))
+	 (angle [z] (java.lang.Math/atan2 im re))
+	 (conjugate [z] (complex-rectangular re (- im)))
+	 Arithmetic
+	 (zero [z] (complex-rectangular 0 0))
+	 (one [z] (complex-rectangular 1 0))
+	 (negate [z] (complex-rectangular (- re) (- im)))
+	 (invert [z] (complex-rectangular
+		      (/ re (magnitude-squared z))
+		      (/ (- im) (magnitude-squared z))))
+	 Object
+	 (toString [_]
+		   (if (and (zero? re) (zero? im)) (str 0)
+		       (str
+			(if (not(zero? re))
+			  re)
+			(if ((comp not zero?) im)
+			  (str
+			   (if (neg? im) "-" "+")
+			   (if ((comp not one?) (java.lang.Math/abs im))
+			   (java.lang.Math/abs im))
+			   "i")))))))
+(defn complex-polar
+"Define a complex number with the given magnitude and angle."
+[mag ang]
+(reify Complex
+	 (magnitude-squared [z] (* mag mag))
+	 (angle [z] angle)
+	 (real-part [z] (* mag (java.lang.Math/cos ang)))
+	 (imaginary-part [z] (* mag (java.lang.Math/sin ang)))
+	 (conjugate [z] (complex-polar mag (- ang)))
+	 Arithmetic
+	 (zero [z] (complex-polar 0 0))
+	 (one [z] (complex-polar 1 0))
+	 (negate [z] (complex-polar (- mag) ang))
+	 (invert [z] (complex-polar (/ mag) (- ang)))
+	 Object
+	 (toString [_] (str mag " * e^(i" ang")"))
+	 ))
+;; Numbers are complex quantities
+(extend-protocol Complex
+java.lang.Number
+(real-part [x] x)
+(imaginary-part [x] 0)
+(magnitude [x] x)
+(angle [x] 0)
+(conjugate [x] x))
+#+end_src
+** Tuples and Tensors
+A tuple is a vector which is spinable\mdash{}it can be either /spin
+up/ or /spin down/. (Covariant, contravariant; dual vectors)
+#+srcname: tuples
+#+begin_src clojure
+(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))
+#+end_src
+*** Contraction
+Contraction is a binary operation that you can apply to compatible
+tuples. Tuples are compatible for contraction if they have the same
+length and opposite spins, and if the corresponding items in each
+tuple are both numbers or both compatible tuples.
+#+srcname: tuples-2
+#+begin_src clojure
+(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))))
+#+end_src
+*** Matrices
+#+srcname: matrices
+#+begin_src clojure
+(in-ns 'sicm.utils)
+(require 'incanter.core) ;; use incanter's fast matrices
+(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))
+(extend-protocol Arithmetic
+incanter.Matrix
+(one [matrix]
+(if (square? matrix)
+	 (identity-matrix (count-rows matrix))
+	 (throw (Exception. "Non-square matrices have no multiplicative unit."))))
+(zero [matrix]
+	(apply matrix-by-rows (map zero (rows matrix))))
+(negate [matrix]
+	  (apply matrix-by-rows (map negate (rows matrix))))
+(invert [matrix]
+	  (incanter.core/solve matrix)))
+(defmulti coerce-to-matrix
+"Converts x into a matrix, if possible."
+type)
+(defmethod coerce-to-matrix incanter.Matrix [x] x)
+(defmethod coerce-to-matrix Tuple [x]
+	   (if (apply numbers? (seq x))
+	     (if (up? x)
+	     (matrix-by-cols (seq x))
+	     (matrix-by-rows (seq x)))
+	     (throw (Exception. "Non-numerical tuple cannot be converted into a matrix."))))
+;; (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
+"Returns the matrix resulting from the matrix multiplication of the given arguments."
+([A] (coerce-to-matrix A))
+([A B] (incanter.core/mmult (coerce-to-matrix A) (coerce-to-matrix B)))
+([M1 M2 & Ms] (reduce matrix-multiply (matrix-multiply M1 M2) Ms)))
+#+end_src
+** Power Series
+#+srcname power-series
+#+begin_src clojure
+(in-ns 'sicm.utils)
+(use 'clojure.contrib.def)
+(defn series-fn
+"The function corresponding to the given power series."
+[series]
+(fn [x]
+(reduce +
+	    (map-indexed  (fn[n x] (* (float (nth series n)) (float(java.lang.Math/pow (float x) n)) ))
+			  (range 20)))))
+(deftype PowerSeries
+[coll]
+clojure.lang.Seqable
+(seq [this] (seq (.coll this)))
+clojure.lang.Indexed
+(nth [this n] (nth (.coll this) n 0))
+(nth [this n not-found] (nth (.coll this) n not-found))
+;; clojure.lang.IFn
+;; (call [this] (throw(Exception.)))
+;; (invoke [this & args] args
+;; 	  (let [f
+;; )
+;; (run [this] (throw(Exception.)))
+)
+(defn power-series
+"Returns a power series with the items of the coll as its
+coefficients. Trailing zeros are added to the end of coll."
+[coeffs]
+(PowerSeries. coeffs))
+(defn power-series-indexed
+"Returns a power series consisting of the result of mapping f to the non-negative integers."
+[f]
+(PowerSeries. (map f (range))))
+(defn-memo nth-partial-sum
+([series n]
+(if (zero? n) (first series)
+	 (+ (nth series n)
+	    (nth-partial-sum series (dec n))))))
+(defn partial-sums [series]
+(lazy-seq (map nth-partial-sum (range))))
+(def cos-series
+(power-series-indexed
+(fn[n]
+	(if (odd? n) 0
+	    (/
+	     (reduce *
+		     (reduce * (repeat (/ n 2) -1))
+		     (range 1 (inc n)))
+	     )))))
+(def sin-series
+(power-series-indexed
+(fn[n]
+	(if (even? n) 0
+	    (/
+	     (reduce *
+		     (reduce * (repeat (/ (dec n) 2) -1))
+		     (range 1 (inc n)))
+	     )))))
+#+end_src
+* Basic Utilities
+** Sequence manipulation
+#+srcname: seq-manipulation
+#+begin_src clojure
+(ns sicm.utils)
+(defn do-up
+"Apply f to each number from low to high, presumably for
+side-effects."
+[f low high]
+(doseq [i (range low high)] (f i)))
+(defn do-down
+"Apply f to each number from high to low, presumably for
+side-effects."
+[f high low]
+(doseq [i (range high low -1)] (f i)))
+(defn all-equal? [coll]
+(if (empty? (rest coll)) true
+(and (= (first coll) (second coll))
+	   (recur (rest coll))))))
+(defn multiplier
+"Returns a function that 'multiplies' the members of a collection,
+returning unit if called on an empty collection."
+[multiply unit]
+(fn [coll] ((partial reduce multiply unit) coll)))
+(defn divider
+"Returns a function that 'divides' the first element of a collection
+by the 'product' of the rest of the collection."
+[divide multiply invert unit]
+(fn [coll]
+(apply
+(fn
+([] unit)
+([x] (invert x))
+([x y] (divide x y))
+([x y & zs] (divide x (reduce multiply y zs))))
+coll)))
+(defn left-circular-shift
+"Remove the first element of coll, adding it to the end of coll."
+[coll]
+(concat (rest coll) (take 1 coll)))
+(defn right-circular-shift
+"Remove the last element of coll, adding it to the front of coll."
+[coll]
+(cons (last coll) (butlast coll)))
+#+end_src
+** Ranges, Airity and Function Composition
+#+srcname: arity
+#+begin_src clojure
+(in-ns 'sicm.utils)
+(def infinity Double/POSITIVE_INFINITY)
+(defn infinite? [x] (Double/isInfinite x))
+(def finite? (comp not infinite?))
+(defn arity-min
+"Returns the smallest number of arguments f can take."
+[f]
+(apply
+min
+(map (comp alength #(.getParameterTypes %))
+	(filter (comp (partial = "invoke") #(.getName %))
+		(.getDeclaredMethods (class f))))))
+(defn arity-max
+"Returns the largest number of arguments f can take, possibly
+Infinity."
+[f]
+(let [methods (.getDeclaredMethods (class f))]
+(if (not-any? (partial = "doInvoke") (map #(.getName %) methods))
+(apply max
+	     (map (comp alength #(.getParameterTypes %))
+		  (filter (comp (partial = "invoke") #(.getName %))  methods)))
+infinity)))
+(def ^{:arglists '([f])
+:doc "Returns a two-element list containing the minimum and
+maximum number of args that f can take."}
+arity-interval
+(juxt arity-min arity-max))
+;; --- intervals
+(defn intersect
+"Returns the interval of overlap between interval-1 and interval-2"
+[interval-1 interval-2]
+(if (or (empty? interval-1) (empty? interval-2)) []
+(let [left (max (first interval-1) (first interval-2))
+	    right (min (second interval-1) (second interval-2))]
+	(if (> left right) []
+	    [left right]))))
+(defn same-endpoints?
+"Returns true if the left endpoint is the same as the right
+endpoint."
+[interval]
+(= (first interval) (second interval)))
+(defn naturals?
+"Returns true if the left endpoint is 0 and the right endpoint is
+infinite."
+[interval]
+(and (zero? (first interval))
+(infinite? (second interval))))
+(defn fan-in
+"Returns a function that pipes its input to each of the gs, then
+applies f to the list of results. Consequently, f must be able to
+take a number of arguments equal to the number of gs."
+[f & gs]
+(fn [& args]
+(apply f (apply (apply juxt gs) args))))
+(defn fan-in
+"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."
+[f & gs]
+(comp (partial apply f) (apply juxt gs)))
+(defmacro airty-blah-sad [f n more?]
+(let [syms (vec (map (comp gensym (partial str "x")) (range n)))
+	 optional (gensym "xs")]
+(if more?
+`(fn ~(conj syms '& optional)
+	 (apply ~f ~@syms ~optional))
+`(fn ~syms (~f ~@syms)))))
+(defmacro airt-whaa* [f n more?]
+`(airty-blah-sad ~f ~n ~more?))
+(defn fan-in*
+"Returns a function that pipes its input to each of the gs, then
+applies f to the list of results. Unlike fan-in, fan-in* strictly
+enforces arity: it will fail if the gs do not have compatible
+arities."
+[f & gs]
+(let [arity-in (reduce intersect (map arity-interval gs))
+	left (first arity-in)
+	right (second arity-in)
+	composite (fan-in f gs)
+	]
+(cond
+(empty? arity-in)
+(throw (Exception. "Cannot compose functions with incompatible arities."))
+(not
+(or (= left right)
+	  (and (finite? left)
+	       (= right infinity))))
+(throw (Exception.
+	     "Compose can only handle arities of the form [n n] or [n infinity]"))
+:else
+(airty-blah-sad composite left (= right infinity)))))
+(defn compose-n "Compose any number of functions together."
+([] identity)
+([f] f)
+([f & fs]
+(let [fns (cons f fs)]
+(compose-bin (reduce fan-in (butlast fs)) (last fs))))
+)
+(defn iterated
+([f n id] (reduce comp id (repeat n f)))
+([f n] (reduce comp identity (repeat n f))))
+(defn iterate-until-stable
+"Repeatedly applies f to x, returning the first result that is close
+enough to its predecessor."
+[f close-enough? x]
+(second (swank.util/find-first
+	   (partial apply close-enough?)
+	   (partition 2 1 (iterate f x)))))
+(defn lexical< [x y]
+(neg? (compare (str x) (str y))))
+;; do-up
+;; do-down
+(def make-pairwise-test comparator)
+;;all-equal?
+(def accumulation multiplier)
+(def inverse-accumulation divider)
+;;left-circular-shift
+;;right-circular-shift
+(def exactly-n? same-endpoints?)
+(def any-number? naturals?)
+;; TODO compose
+;; TODO compose-n
+;; identity
+(def compose-2 fan-in)
+(def compose-bin fan-in*)
+(def any? (constantly true))
+(def none? (constantly false))
+(def constant constantly)
+(def joint-arity intersect)
+(def a-reduce reduce)
+;; filter
+(def make-map (partial partial map) )
+(def bracket juxt)
+;; TODO apply-to-all
+;; TODO nary-combine
+;; TODO binary-combine
+;; TODO unary-combine
+;; iterated
+;; iterate-until-stable
+(def make-function-of-vector (partial partial map))
+(def make-function-of-arguments (fn [f] (fn [& args] (f args))))
+(def alphaless lexical<)
+#+end_src
+:
+*  Numerical Methods
+#+srcname: numerical-methods
+#+begin_src clojure
+(in-ns 'sicm.utils)
+(import java.lang.Math)
+(use 'clojure.contrib.def)
+;; ---- USEFUL CONSTANTS
+(defn machine-epsilon
+"Smallest value representable on your machine, as determined by
+successively dividing a number in half until consecutive results are
+indistinguishable."
+[]
+(ffirst
+(drop-while
+(comp not zero? second)
+(partition 2 1
+	       (iterate (partial * 0.5) 1)))))
+(def pi (Math/PI))
+(def two-pi (* 2 pi))
+(def eulers-gamma 0.5772156649015328606065)
+(def phi (/ (inc (Math/sqrt 5)) 2))
+(def ln2 (Math/log 2))
+(def ln10 (Math/log 10))
+(def exp10 #(Math/pow 10 %))
+(def exp2 #(Math/pow 2 %))
+;;
+;; ---- ANGLES AND TRIGONOMETRY
+(defn angle-restrictor
+"Returns a function that ensures that angles lie in the specified interval of length two-pi."
+[max-angle]
+(let [min-angle (- max-angle two-pi)]
+(fn [x]
+(if (and
+	   (<= min-angle x)
+	   (< x max-angle))
+	x
+	(let [corrected-x (- x (* two-pi (Math/floor (/ x two-pi))))]
+	  (if (< corrected-x max-angle)
+	    corrected-x
+	    (- corrected-x two-pi)))))))
+(defn angle-restrict-pi
+"Coerces angles to lie in the interval from -pi to pi."
+[angle]
+((angle-restrictor pi) angle))
+(defn angle-restrict-two-pi
+"Coerces angles to lie in the interval from zero to two-pi"
+[angle]
+((angle-restrictor two-pi) angle))
+(defn invert [x] (/ x))
+(defn negate [x] (- x))
+(defn exp [x] (Math/exp x))
+(defn sin [x] (Math/sin x))
+(defn cos [x] (Math/cos x))
+(defn tan [x] (Math/tan x))
+(def sec (comp invert cos))
+(def csc (comp invert sin))
+(defn sinh [x] (Math/sinh x))
+(defn cosh [x] (Math/cosh x))
+(defn tanh [x] (Math/tanh x))
+(def sech (comp invert cosh))
+(def csch (comp invert sinh))
+;; ------------
+(defn factorial
+"Computes the factorial of the nonnegative integer n."
+[n]
+(if (neg? n)
+(throw (Exception. "Cannot compute the factorial of a negative number."))
+(reduce * 1 (range 1 (inc n)))))
+(defn exact-quotient [n d] (/ n d))
+(defn binomial-coefficient
+"Computes the number of different ways to choose m elements from n."
+[n m]
+(assert (<= 0 m n))
+(let [difference (- n m)]
+(exact-quotient
+(reduce * (range n (max difference m) -1 ))
+(factorial (min difference m)))))
+(defn-memo stirling-1
+"Stirling numbers of the first kind: the number of permutations of n
+elements with exactly m permutation cycles. "
+[n k]
+;(assert (<= 1 k n))
+(if (zero? n)
+(if (zero? k) 1 0)
+(+ (stirling-1 (dec n) (dec k))
+(* (dec n) (stirling-1 (dec n) k)))))
+(defn-memo stirling-2 ;;count-partitions
+"Stirling numbers of the second kind: the number of ways to partition a set of n elements into k subsets."
+[n k]
+(cond
+(= k 1) 1
+(= k n) 1
+:else (+ (stirling-2 (dec n) (dec k))
+	    (* k (stirling-2 (dec n) k)))))
+(defn harmonic-number [n]
+(/ (stirling-1 (inc n) 2)
+(factorial n)))
+(defn sum
+[f low high]
+(reduce + (map f (range low (inc high)))))
+#+end_src
+* Differentiation
+We compute derivatives by passing special *differential objects* $[x,
+dx]$ through functions. Roughly speaking, applying a function $f$ to a
+differential object \([x, dx]\) should produce a new differential
+object $[f(x),\,Df(x)\cdot dx]$.
+\([x,\,dx]\xrightarrow{\quad f \quad}[f(x),\,Df(x)\cdot dx]\)
+Notice that you can obtain the derivative of $f$ from this
+differential object, as it is the coefficient of the $dx$ term. Also,
+as you apply successive functions using this rule, you get the
+chain-rule answer you expect:
+\([f(x),\,Df(x)\cdot dx]\xrightarrow{\quad g\quad} [gf(x),\,
+Dgf(x)\cdot Df(x) \cdot dx ]\)
+In order to generalize to multiple variables and multiple derivatives,
+we use a *power series of differentials*, a sortred infinite sequence which
+contains all terms like $dx\cdot dy$, $dx^2\cdot dy$, etc.
+#+srcname:differential
+#+begin_src clojure
+(in-ns 'sicm.utils)
+(use 'clojure.contrib.combinatorics)
+(use 'clojure.contrib.generic.arithmetic)
+(defprotocol DifferentialTerm
+"Protocol for an infinitesimal quantity."
+(coefficient [this])
+(partials [this]))
+(extend-protocol DifferentialTerm
+java.lang.Number
+(coefficient [this] this)
+(partials [this] []))
+(deftype DifferentialSeq
+[terms]
+clojure.lang.IPersistentCollection
+(cons [this x]
+	(throw (Exception. x))
+	(DifferentialSeq. (cons x terms)))
+;; (seq [this] (seq terms))
+(count [this] (count terms))
+(empty [this] (empty? terms)))
+(defn coerce-to-differential-seq [x]
+(cond
+(= (type x) DifferentialSeq) x
+(satisfies? DifferentialTerm x) (DifferentialSeq. x)))
+(defn differential-term
+"Returns a differential term with the given coefficient and
+partials. Coefficient may be any arithmetic object; partials must
+be a list of non-negative integers."
+[coefficient partials]
+(reify DifferentialTerm
+	 (partials [_] (set partials))
+	 (coefficient [_] coefficient)))
+(defn differential-seq*
+([coefficient partials]
+(DifferentialSeq. [(differential-term coefficient partials)]))
+([coefficient partials & cps]
+(if cps
+(defn differential-seq
+"Constructs a sequence of differential terms from a numerical
+coefficient and a list of keys for variables. If no coefficient is supplied, uses 1."
+([variables] (differential-seq 1 variables))
+([coefficient variables & cvs]
+(if (number? coefficient)
+(conj (assoc {} (apply sorted-set variables) coefficient)
+	     (if (empty? cvs)
+	       nil
+	       (apply differential-seq cvs)))
+(apply differential-seq 1 coefficient 1 variables cvs)
+)))
+(defn differential-add
+"Add two differential sequences by combining like terms."
+[dseq1 dseq2]
+(merge-with + dseq1 dseq2))
+(defn differential-multiply
+"Multiply two differential sequences. The square of any differential variable is zero since differential variables are infinitesimally small."
+[dseq1 dseq2]
+(reduce
+(fn [m [[vars1 coeff1] [vars2 coeff2]]]
+(if (empty? (clojure.set/intersection vars1 vars2))
+(assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))
+m))
+{}
+(clojure.contrib.combinatorics/cartesian-product
+dseq1
+dseq2)))
+(defn big-part
+"Returns the part of the differential sequence that is finite,
+i.e. not infinitely small."
+[dseq]
+(let
+[keys (sort-by count (keys dseq))
+smallest-var (last(last keys))]
+(apply hash-map
+	   (reduce concat
+		   (remove (comp smallest-var first) dseq)))))
+(defn small-part
+"Returns the part of the differential sequence that is
+infinitesimal."
+[dseq]
+(let
+[keys (sort-by count (keys dseq))
+smallest-var (last(last keys))]
+(apply hash-map
+	   (reduce concat
+		   (filter (comp smallest-var first) dseq)))))
+(defn big-part
+"Returns the 'finite' part of the differential sequence."
+[dseq]
+(let
+[keys (sort-by count (keys dseq))
+smallest-var (last (last keys))]
+(apply hash-map
+(reduce concat
+(filter (remove smallest-var first) dseq)
+	   ))))
+(defn small-part
+"Returns the 'infinitesimal' part of the differential sequence."
+[dseq]
+(let
+[keys (sort-by count (keys dseq))
+smallest-var (last (last keys))]
+(apply hash-map
+(reduce concat
+	    (filter (comp smallest-var first) dseq)
+	    ))))
+(defn linear-approximator
+"Returns the function that corresponds to unary-fn but which accepts and returns differential objects."
+([unary-f dfdx]
+(fn F [dseq]
+(differential-add
+	(unary-f (big-part dseq))
+	(differential-multiply
+	 (dfdx (big-part dseq))
+	 (small-part dseq))))))
+;; ;; A differential term consists of a numerical coefficient and a
+;; ;; sorted
+;; (defrecord DifferentialTerm [coefficient variables])
+;; (defmethod print-method DifferentialTerm
+;;    [o w]
+;;     (print-simple
+;;      (apply str (.coefficient o)(map (comp (partial str "d") name) (.variables o)))
+;;      w))
+;; (defn differential-seq
+;;   "Constructs a sequence of differential terms from a numerical
+;; coefficient and a list of keywords for variables. If no coefficient is
+;; supplied, uses 1."
+;;   ([variables] (differential-seq 1 variables))
+;;   ([coefficient variables]
+;;      (list
+;;       (DifferentialTerm. coefficient (apply sorted-set variables))))
+;;   ([coefficient variables & cvs]
+;;      (sort-by
+;;       #(vec(.variables %))
+;;       (concat (differential-seq coefficient variables) (apply differential-seq cvs)))))
+#+end_src
+* Symbolic Manipulation
+#+srcname:symbolic
+#+begin_src clojure
+(in-ns 'sicm.utils)
+(deftype Symbolic [type expression]
+Object
+(equals [this that]
+	  (cond
+	   (= (.expression this) (.expression that)) true
+	   :else
+	   (Symbolic.
+	    java.lang.Boolean
+	    (list '= (.expression this) (.expression that)))
+	  )))
+(deftype AbstractSet [glyph membership-test])
+(defn member? [abstract-set x]
+((.membership-test abstract-set) x))
+;; ------------ Some important AbstractSets
+(def Real
+(AbstractSet.
+'R
+(fn[x](number? x))))
+;; ------------ Create new AbstractSets from existing ones
+(defn abstract-product
+"Gives the cartesian product of abstract sets."
+([sets]
+(if (= (count sets) 1) (first sets)
+	   (AbstractSet.
+	    (symbol
+	     (apply str
+		    (interpose 'x (map #(.glyph %) sets))))
+	   (fn [x]
+	     (and
+	      (coll? x)
+	      (= (count sets) (count x))
+	      (reduce #(and %1 %2)
+		      true
+		      (map #(member? %1 %2) sets x)))))))
+([abstract-set n]
+(abstract-product (repeat n abstract-set))))
+(defn abstract-up
+"Returns the abstract set of all up-tuples whose items belong to the
+corresponding abstract sets in coll."
+([coll]
+(AbstractSet.
+(symbol (str "u["
+		   (apply str
+			  (interpose " "
+				     (map #(.glyph %) coll)))
+		   "]"))
+(fn [x]
+	(and
+	 (satisfies? Spinning x)
+	 (up? x)
+	 (= (count coll) (count x))
+	 (reduce
+	  #(and %1 %2)
+	  true
+	  (map #(member? %1 %2) coll x))))))
+([abstract-set n]
+(abstract-up (repeat n abstract-set))))
+(defn abstract-down
+"Returns the abstract set of all down-tuples whose items belong to the
+corresponding abstract sets in coll."
+([coll]
+(AbstractSet.
+(symbol (str "d["
+		   (apply str
+			  (interpose " "
+				     (map #(.glyph %) coll)))
+		   "]"))
+(fn [x]
+	(and
+	 (satisfies? Spinning x)
+	 (down? x)
+	 (= (count coll) (count x))
+	 (reduce
+	  #(and %1 %2)
+	  true
+	  (map #(member? %1 %2) coll x))))))
+([abstract-set n]
+(abstract-down (repeat n abstract-set))))
+					;-------ABSTRACT FUNCTIONS
+(defrecord AbstractFn
+[#^AbstractSet domain #^AbstractSet codomain])
+(defmethod print-method AbstractFn
+[o w]
+(print-simple
+(str
+"f:"
+(.glyph (:domain o))
+"-->"
+(.glyph (:codomain o))) w))
+#+end_src
+* COMMENT
+#+begin_src clojure :tangle utils.clj
+(ns sicm.utils)
+					;***** GENERIC ARITHMETIC
+<<generic-arithmetic>>
+					;***** TUPLES AND MATRICES
+<<tuples>>
+<<tuples-2>>
+					;***** MATRICES
+<<matrices>>
+#+end_src

Mercurial > dylan

comparison sicm/bk/utils.org @ 2:b4de894a1e2e