rlm@2: #+TITLE:Building a Classical Mechanics Library in Clojure
rlm@2: #+author: Robert McIntyre & Dylan Holmes
rlm@2: #+EMAIL:     rlm@mit.edu
rlm@2: #+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js"
rlm@2: #+STYLE: <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />
rlm@2: #+OPTIONS:   H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
rlm@2: #+SETUPFILE: ../templates/level-0.org
rlm@2: #+INCLUDE: ../templates/level-0.org
rlm@2: #+BABEL: :noweb yes :results silent
rlm@2: 
rlm@2: * Generic Arithmetic
rlm@2: 
rlm@2: #+srcname: generic-arithmetic
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: 
rlm@2: * Useful Data Types
rlm@2: 
rlm@2: ** Complex Numbers
rlm@2: #+srcname: complex-numbers
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: ** Tuples and Tensors
rlm@2: 
rlm@2: A tuple is a vector which is spinable\mdash{}it can be either /spin
rlm@2: up/ or /spin down/. (Covariant, contravariant; dual vectors)
rlm@2: #+srcname: tuples
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: *** 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: #+srcname: tuples-2
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: *** Matrices
rlm@2: #+srcname: matrices
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: 
rlm@2: ** Power Series
rlm@2: #+srcname power-series
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: 
rlm@2: * Basic Utilities
rlm@2: 
rlm@2: ** Sequence manipulation
rlm@2: 
rlm@2: #+srcname: seq-manipulation
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: ** Ranges, Airity and Function Composition
rlm@2: #+srcname: arity
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: : 
rlm@2: 
rlm@2: *  Numerical Methods
rlm@2: #+srcname: numerical-methods
rlm@2: #+begin_src clojure
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: #+end_src
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: * 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: \([x,\,dx]\xrightarrow{\quad f \quad}[f(x),\,Df(x)\cdot dx]\)
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: \([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: 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: #+srcname:differential
rlm@2: #+begin_src clojure
rlm@2: (in-ns 'sicm.utils)
rlm@2: (use 'clojure.contrib.combinatorics)
rlm@2: (use 'clojure.contrib.generic.arithmetic)
rlm@2: 
rlm@2: (defprotocol DifferentialTerm
rlm@2:   "Protocol for an infinitesimal quantity."
rlm@2:   (coefficient [this])
rlm@2:   (partials [this]))
rlm@2: 
rlm@2: (extend-protocol DifferentialTerm
rlm@2:   java.lang.Number
rlm@2:   (coefficient [this] this)
rlm@2:   (partials [this] []))
rlm@2: 
rlm@2: (deftype DifferentialSeq
rlm@2:   [terms]
rlm@2:   clojure.lang.IPersistentCollection
rlm@2:   (cons [this x]
rlm@2: 	(throw (Exception. x))
rlm@2: 	(DifferentialSeq. (cons x terms)))
rlm@2:  ;; (seq [this] (seq terms))
rlm@2:   (count [this] (count terms))
rlm@2:   (empty [this] (empty? terms)))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defn coerce-to-differential-seq [x]
rlm@2:   (cond
rlm@2:    (= (type x) DifferentialSeq) x
rlm@2:    (satisfies? DifferentialTerm x) (DifferentialSeq. x)))
rlm@2: 
rlm@2: 
rlm@2: (defn differential-term
rlm@2:   "Returns a differential term with the given coefficient and
rlm@2:   partials. Coefficient may be any arithmetic object; partials must
rlm@2:   be a list of non-negative integers."
rlm@2:   [coefficient partials]
rlm@2:   (reify DifferentialTerm
rlm@2: 	 (partials [_] (set partials))
rlm@2: 	 (coefficient [_] coefficient)))
rlm@2: 
rlm@2: 
rlm@2: 
rlm@2: (defn differential-seq*
rlm@2:   ([coefficient partials]
rlm@2:      (DifferentialSeq. [(differential-term coefficient partials)]))
rlm@2:   ([coefficient partials & cps]
rlm@2:      (if cps
rlm@2:        
rlm@2: 	
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 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: (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: (defn big-part
rlm@2:   "Returns the part of the differential sequence that is finite,
rlm@2:   i.e. not infinitely small."
rlm@2:   [dseq]
rlm@2:   (let
rlm@2:       [keys (sort-by count (keys dseq))
rlm@2:        smallest-var (last(last keys))]
rlm@2:     (apply hash-map
rlm@2: 	   (reduce concat
rlm@2: 		   (remove (comp smallest-var first) dseq)))))
rlm@2: 
rlm@2: (defn small-part
rlm@2:   "Returns the part of the differential sequence that is
rlm@2:   infinitesimal."
rlm@2:   [dseq]
rlm@2:   (let
rlm@2:       [keys (sort-by count (keys dseq))
rlm@2:        smallest-var (last(last keys))]
rlm@2:     (apply hash-map
rlm@2: 	   (reduce concat
rlm@2: 		   (filter (comp smallest-var first) dseq)))))
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:        smallest-var (last (last keys))]
rlm@2:     (apply hash-map
rlm@2:     (reduce concat
rlm@2:            (filter (remove smallest-var 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:        smallest-var (last (last keys))]
rlm@2: 
rlm@2:     (apply hash-map
rlm@2:     (reduce concat
rlm@2: 	    (filter (comp smallest-var first) dseq)
rlm@2: 	    ))))
rlm@2: 
rlm@2: (defn linear-approximator
rlm@2:   "Returns the function that corresponds to unary-fn but which accepts and returns differential objects."
rlm@2:   ([unary-f dfdx]
rlm@2:      (fn F [dseq]
rlm@2:        (differential-add
rlm@2: 	(unary-f (big-part dseq))
rlm@2: 	(differential-multiply
rlm@2: 	 (dfdx (big-part dseq))
rlm@2: 	 (small-part dseq))))))
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: #+end_src
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: * Symbolic Manipulation
rlm@2: 
rlm@2: #+srcname:symbolic
rlm@2: #+begin_src clojure
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: #+end_src
rlm@2: 
rlm@2: 
rlm@2: * COMMENT 
rlm@2: #+begin_src clojure :tangle utils.clj
rlm@2: (ns sicm.utils)
rlm@2: 
rlm@2: 					;***** GENERIC ARITHMETIC
rlm@2: <<generic-arithmetic>>
rlm@2: 
rlm@2: 
rlm@2: 					;***** TUPLES AND MATRICES
rlm@2: <<tuples>>
rlm@2: <<tuples-2>>
rlm@2: 					;***** MATRICES
rlm@2: <<matrices>>
rlm@2: 
rlm@2: #+end_src
rlm@2: 
rlm@2: 
rlm@2: