#+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