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