1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/sicm/bk/utils.org	Fri Oct 28 00:03:05 2011 -0700
     1.3 @@ -0,0 +1,1262 @@
     1.4 +#+TITLE:Building a Classical Mechanics Library in Clojure
     1.5 +#+author: Robert McIntyre & Dylan Holmes
     1.6 +#+EMAIL:     rlm@mit.edu
     1.7 +#+MATHJAX: align:"left" mathml:t path:"../MathJax/MathJax.js"
     1.8 +#+STYLE: <link rel="stylesheet" type="text/css" href="../css/aurellem.css" />
     1.9 +#+OPTIONS:   H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
    1.10 +#+SETUPFILE: ../templates/level-0.org
    1.11 +#+INCLUDE: ../templates/level-0.org
    1.12 +#+BABEL: :noweb yes :results silent
    1.13 +
    1.14 +* Generic Arithmetic
    1.15 +
    1.16 +#+srcname: generic-arithmetic
    1.17 +#+begin_src clojure
    1.18 +(ns sicm.utils)
    1.19 +(in-ns 'sicm.utils)
    1.20 +
    1.21 +
    1.22 +
    1.23 +(defn all-equal? [coll]
    1.24 +  (if (empty? (rest coll)) true
    1.25 +      (and (= (first coll) (second coll))
    1.26 +	   (recur (rest coll)))))
    1.27 +
    1.28 +
    1.29 +(defprotocol Arithmetic
    1.30 +  (zero [this])
    1.31 +  (one [this])
    1.32 +  (negate [this])
    1.33 +  (invert [this])
    1.34 +  (conjugate [this]))
    1.35 +
    1.36 +(defn zero? [x]
    1.37 +  (= x (zero x)))
    1.38 +(defn one? [x]
    1.39 +  (= x (one x)))
    1.40 +
    1.41 +
    1.42 +
    1.43 +(extend-protocol Arithmetic
    1.44 +  java.lang.Number
    1.45 +  (zero [x] 0)
    1.46 +  (one [x] 1)
    1.47 +  (negate [x] (- x))
    1.48 +  (invert [x] (/ x))
    1.49 +  )
    1.50 +
    1.51 +(extend-protocol Arithmetic
    1.52 +  clojure.lang.IFn
    1.53 +  (one [f] identity)
    1.54 +  (negate [f] (comp negate f))
    1.55 +  (invert [f] (comp invert f)))
    1.56 +
    1.57 +    
    1.58 +(extend-protocol Arithmetic
    1.59 +  clojure.lang.Seqable
    1.60 +  (zero [this] (map zero this))
    1.61 +  (one [this] (map one this))
    1.62 +  (invert [this] (map invert this))
    1.63 +  (negate [this] (map negate this)))
    1.64 +
    1.65 +
    1.66 +(defn ordered-like
    1.67 +  "Create a comparator using the sorted collection as an
    1.68 +  example. Elements not in the sorted collection are sorted to the
    1.69 +  end."
    1.70 +  [sorted-coll]
    1.71 +  (let [
    1.72 +	sorted-coll? (set sorted-coll) 
    1.73 +	ascending-pair?
    1.74 +	(set(reduce concat
    1.75 +		    (map-indexed
    1.76 +		     (fn [n x]
    1.77 +		       (map #(vector x %) (nthnext sorted-coll n)))
    1.78 +		     sorted-coll)))]
    1.79 +    (fn [x y]
    1.80 +      (cond
    1.81 +       (= x y) 0
    1.82 +       (ascending-pair? [x y]) -1
    1.83 +       (ascending-pair? [y x]) 1
    1.84 +       (sorted-coll? x) -1
    1.85 +       (sorted-coll? y) 1))))
    1.86 +  
    1.87 +
    1.88 +
    1.89 +(def type-precedence
    1.90 +  (ordered-like [incanter.Matrix]))
    1.91 +
    1.92 +(defmulti add
    1.93 +    (fn [x y]
    1.94 +      (sort type-precedence [(type x)(type y)]))) 
    1.95 +
    1.96 +(defmulti multiply
    1.97 +  (fn [x y]
    1.98 +    (sort type-precedence [(type x) (type y)])))
    1.99 +
   1.100 +(defmethod add [java.lang.Number java.lang.Number] [x y] (+ x y))
   1.101 +(defmethod multiply [java.lang.Number java.lang.Number] [x y] (* x y))
   1.102 +
   1.103 +(defmethod multiply [incanter.Matrix java.lang.Integer] [x y]
   1.104 +	   (let [args (sort #(type-precedence (type %1)(type %2)) [x y])
   1.105 +		 matrix (first args)
   1.106 +		 scalar (second args)]
   1.107 +	     (incanter.core/matrix (map (partial map (partial multiply scalar)) matrix))))
   1.108 +	   
   1.109 +#+end_src
   1.110 +
   1.111 +
   1.112 +* Useful Data Types
   1.113 +
   1.114 +** Complex Numbers
   1.115 +#+srcname: complex-numbers
   1.116 +#+begin_src clojure
   1.117 +(in-ns 'sicm.utils)
   1.118 +
   1.119 +(defprotocol Complex
   1.120 +  (real-part [z])
   1.121 +  (imaginary-part [z])
   1.122 +  (magnitude-squared [z])
   1.123 +  (angle [z])
   1.124 +  (conjugate [z])
   1.125 +  (norm [z]))
   1.126 +
   1.127 +(defn complex-rectangular
   1.128 +  "Define a complex number with the given real and imaginary
   1.129 +  components."
   1.130 +  [re im]
   1.131 +  (reify Complex
   1.132 +	 (real-part [z] re)
   1.133 +	 (imaginary-part [z] im)
   1.134 +	 (magnitude-squared [z] (+ (* re re) (* im im)))
   1.135 +	 (angle [z] (java.lang.Math/atan2 im re))
   1.136 +	 (conjugate [z] (complex-rectangular re (- im)))
   1.137 +
   1.138 +	 Arithmetic
   1.139 +	 (zero [z] (complex-rectangular 0 0))
   1.140 +	 (one [z] (complex-rectangular 1 0))
   1.141 +	 (negate [z] (complex-rectangular (- re) (- im)))
   1.142 +	 (invert [z] (complex-rectangular
   1.143 +		      (/ re (magnitude-squared z))
   1.144 +		      (/ (- im) (magnitude-squared z))))
   1.145 +
   1.146 +	 Object
   1.147 +	 (toString [_]
   1.148 +		   (if (and (zero? re) (zero? im)) (str 0)
   1.149 +		       (str
   1.150 +			(if (not(zero? re))
   1.151 +			  re)
   1.152 +			(if ((comp not zero?) im)
   1.153 +			  (str
   1.154 +			   (if (neg? im) "-" "+")
   1.155 +			   (if ((comp not one?) (java.lang.Math/abs im))
   1.156 +			   (java.lang.Math/abs im))
   1.157 +			   "i")))))))
   1.158 +
   1.159 +(defn complex-polar
   1.160 +  "Define a complex number with the given magnitude and angle."
   1.161 +  [mag ang]
   1.162 +  (reify Complex
   1.163 +	 (magnitude-squared [z] (* mag mag))
   1.164 +	 (angle [z] angle)
   1.165 +	 (real-part [z] (* mag (java.lang.Math/cos ang)))
   1.166 +	 (imaginary-part [z] (* mag (java.lang.Math/sin ang)))
   1.167 +	 (conjugate [z] (complex-polar mag (- ang)))
   1.168 +
   1.169 +	 Arithmetic
   1.170 +	 (zero [z] (complex-polar 0 0))
   1.171 +	 (one [z] (complex-polar 1 0))
   1.172 +	 (negate [z] (complex-polar (- mag) ang))
   1.173 +	 (invert [z] (complex-polar (/ mag) (- ang)))
   1.174 +
   1.175 +	 Object
   1.176 +	 (toString [_] (str mag " * e^(i" ang")"))
   1.177 +	 ))
   1.178 +
   1.179 +
   1.180 +;; Numbers are complex quantities
   1.181 +
   1.182 +(extend-protocol Complex
   1.183 +  java.lang.Number
   1.184 +  (real-part [x] x)
   1.185 +  (imaginary-part [x] 0)
   1.186 +  (magnitude [x] x)
   1.187 +  (angle [x] 0)
   1.188 +  (conjugate [x] x))
   1.189 +
   1.190 +
   1.191 +#+end_src
   1.192 +
   1.193 +
   1.194 +
   1.195 +** Tuples and Tensors
   1.196 +
   1.197 +A tuple is a vector which is spinable\mdash{}it can be either /spin
   1.198 +up/ or /spin down/. (Covariant, contravariant; dual vectors)
   1.199 +#+srcname: tuples
   1.200 +#+begin_src clojure
   1.201 +(in-ns 'sicm.utils)
   1.202 +
   1.203 +(defprotocol Spinning
   1.204 +  (up? [this])
   1.205 +  (down? [this]))
   1.206 +
   1.207 +(defn spin
   1.208 +  "Returns the spin of the Spinning s, either :up or :down"
   1.209 +  [#^Spinning s]
   1.210 +  (cond (up? s) :up (down? s) :down))
   1.211 +
   1.212 +
   1.213 +(deftype Tuple
   1.214 +  [spin coll]
   1.215 +  clojure.lang.Seqable
   1.216 +  (seq [this] (seq (.coll this)))
   1.217 +  clojure.lang.Counted
   1.218 +  (count [this] (count (.coll this))))
   1.219 +
   1.220 +(extend-type Tuple
   1.221 +  Spinning
   1.222 +  (up? [this] (= ::up (.spin this)))
   1.223 +  (down? [this] (= ::down (.spin this))))
   1.224 +
   1.225 +(defmethod print-method Tuple
   1.226 +  [o w]
   1.227 +  (print-simple (str (if (up? o) 'u 'd) (.coll o))  w))
   1.228 +
   1.229 +
   1.230 +
   1.231 +(defn up
   1.232 +  "Create a new up-tuple containing the contents of coll."
   1.233 +  [coll]
   1.234 +  (Tuple. ::up coll))       
   1.235 +
   1.236 +(defn down
   1.237 +  "Create a new down-tuple containing the contents of coll."
   1.238 +  [coll]
   1.239 +  (Tuple. ::down coll))
   1.240 +
   1.241 +
   1.242 +#+end_src
   1.243 +
   1.244 +*** Contraction
   1.245 +Contraction is a binary operation that you can apply to compatible
   1.246 + tuples. Tuples are compatible for contraction if they have the same
   1.247 + length and opposite spins, and if the corresponding items in each
   1.248 + tuple are both numbers or both compatible tuples.
   1.249 +
   1.250 +#+srcname: tuples-2
   1.251 +#+begin_src clojure
   1.252 +(in-ns 'sicm.utils)
   1.253 +
   1.254 +(defn numbers?
   1.255 +  "Returns true if all arguments are numbers, else false."
   1.256 +  [& xs]
   1.257 +  (every? number? xs))
   1.258 +
   1.259 +(defn contractible?
   1.260 +  "Returns true if the tuples a and b are compatible for contraction,
   1.261 +  else false. Tuples are compatible if they have the same number of
   1.262 +  components, they have opposite spins, and their elements are
   1.263 +  pairwise-compatible."
   1.264 +  [a b]
   1.265 +  (and
   1.266 +   (isa? (type a) Tuple)
   1.267 +   (isa? (type b) Tuple)
   1.268 +   (= (count a) (count b))
   1.269 +   (not= (spin a) (spin b))
   1.270 +   
   1.271 +   (not-any? false?
   1.272 +	     (map #(or
   1.273 +		    (numbers? %1 %2)
   1.274 +		    (contractible? %1 %2))
   1.275 +		  a b))))
   1.276 +
   1.277 +
   1.278 +
   1.279 +(defn contract
   1.280 +  "Contracts two tuples, returning the sum of the
   1.281 +  products of the corresponding items. Contraction is recursive on
   1.282 +  nested tuples."
   1.283 +  [a b]
   1.284 +  (if (not (contractible? a b))
   1.285 +    (throw
   1.286 +     (Exception. "Not compatible for contraction."))
   1.287 +    (reduce +
   1.288 +	    (map
   1.289 +	     (fn [x y]
   1.290 +	       (if (numbers? x y)
   1.291 +		 (* x y)
   1.292 +		 (contract x y)))
   1.293 +	     a b))))
   1.294 +
   1.295 +#+end_src
   1.296 +
   1.297 +*** Matrices
   1.298 +#+srcname: matrices
   1.299 +#+begin_src clojure
   1.300 +(in-ns 'sicm.utils)
   1.301 +(require 'incanter.core) ;; use incanter's fast matrices
   1.302 +
   1.303 +(defprotocol Matrix
   1.304 +  (rows [matrix])
   1.305 +  (cols [matrix])
   1.306 +  (diagonal [matrix])
   1.307 +  (trace [matrix])
   1.308 +  (determinant [matrix])
   1.309 +  (transpose [matrix])
   1.310 +  (conjugate [matrix])
   1.311 +)
   1.312 +
   1.313 +(extend-protocol Matrix
   1.314 +  incanter.Matrix
   1.315 +  (rows [rs] (map down (apply map vector (apply map vector rs))))
   1.316 +  (cols [rs] (map up (apply map vector rs)))
   1.317 +  (diagonal [matrix] (incanter.core/diag matrix) )
   1.318 +  (determinant [matrix] (incanter.core/det matrix))
   1.319 +  (trace [matrix] (incanter.core/trace matrix))
   1.320 +  (transpose [matrix] (incanter.core/trans matrix))
   1.321 +  )
   1.322 +
   1.323 +(defn count-rows [matrix]
   1.324 +  ((comp count rows) matrix))
   1.325 +
   1.326 +(defn count-cols [matrix]
   1.327 +  ((comp count cols) matrix))
   1.328 +
   1.329 +(defn square? [matrix]
   1.330 +  (= (count-rows matrix) (count-cols matrix)))
   1.331 +
   1.332 +(defn identity-matrix
   1.333 +  "Define a square matrix of size n-by-n with 1s along the diagonal and
   1.334 +  0s everywhere else."
   1.335 +  [n]
   1.336 +  (incanter.core/identity-matrix n))
   1.337 +
   1.338 +
   1.339 +
   1.340 +
   1.341 +
   1.342 +(defn matrix-by-rows
   1.343 +  "Define a matrix by giving its rows."
   1.344 +  [& rows]
   1.345 +  (if
   1.346 +   (not (all-equal? (map count rows)))
   1.347 +   (throw (Exception. "All rows in a matrix must have the same number of elements."))
   1.348 +   (incanter.core/matrix (vec rows))))
   1.349 +
   1.350 +(defn matrix-by-cols
   1.351 +  "Define a matrix by giving its columns"
   1.352 +  [& cols]
   1.353 +  (if (not (all-equal? (map count cols)))
   1.354 +   (throw (Exception. "All columns in a matrix must have the same number of elements."))
   1.355 +   (incanter.core/matrix (vec (apply map vector cols)))))
   1.356 +
   1.357 +(defn identity-matrix
   1.358 +  "Define a square matrix of size n-by-n with 1s along the diagonal and
   1.359 +  0s everywhere else."
   1.360 +  [n]
   1.361 +  (incanter.core/identity-matrix n))
   1.362 +
   1.363 +
   1.364 +
   1.365 +(extend-protocol Arithmetic
   1.366 +  incanter.Matrix
   1.367 +  (one [matrix]
   1.368 +       (if (square? matrix)
   1.369 +	 (identity-matrix (count-rows matrix))
   1.370 +	 (throw (Exception. "Non-square matrices have no multiplicative unit."))))
   1.371 +  (zero [matrix]
   1.372 +	(apply matrix-by-rows (map zero (rows matrix))))
   1.373 +  (negate [matrix]
   1.374 +	  (apply matrix-by-rows (map negate (rows matrix))))
   1.375 +  (invert [matrix]
   1.376 +	  (incanter.core/solve matrix)))
   1.377 +
   1.378 +
   1.379 +
   1.380 +(defmulti coerce-to-matrix
   1.381 + "Converts x into a matrix, if possible."
   1.382 + type)
   1.383 +
   1.384 +(defmethod coerce-to-matrix incanter.Matrix [x] x)
   1.385 +(defmethod coerce-to-matrix Tuple [x]
   1.386 +	   (if (apply numbers? (seq x))
   1.387 +	     (if (up? x)
   1.388 +	     (matrix-by-cols (seq x))
   1.389 +	     (matrix-by-rows (seq x)))
   1.390 +	     (throw (Exception. "Non-numerical tuple cannot be converted into a matrix.")))) 
   1.391 +	     
   1.392 +  
   1.393 +
   1.394 +
   1.395 +
   1.396 +;; (defn matrix-by-cols
   1.397 +;;   "Define a matrix by giving its columns."
   1.398 +;;   [& cols]
   1.399 +;;   (cond
   1.400 +;;    (not (all-equal? (map count cols)))
   1.401 +;;    (throw (Exception. "All columns in a matrix must have the same number of elements."))
   1.402 +;;    :else
   1.403 +;;    (reify Matrix
   1.404 +;; 	  (cols [this] (map up cols))
   1.405 +;; 	  (rows [this] (map down (apply map vector cols)))
   1.406 +;; 	  (diagonal [this] (map-indexed (fn [i col] (nth col i) cols)))
   1.407 +;; 	  (trace [this]
   1.408 +;; 		 (if (not= (count-cols this) (count-rows this))
   1.409 +;; 		   (throw (Exception.
   1.410 +;; 			   "Cannot take the trace of a non-square matrix."))
   1.411 +;; 		   (reduce + (diagonal this))))
   1.412 +	  
   1.413 +;; 	  (determinant [this]
   1.414 +;; 		       (if (not= (count-cols this) (count-rows this))
   1.415 +;; 			 (throw (Exception.
   1.416 +;; 				 "Cannot take the determinant of a non-square matrix."))
   1.417 +;; 			 (reduce * (map-indexed (fn [i col] (nth col i)) cols))))
   1.418 +;; 	  )))
   1.419 +
   1.420 +(extend-protocol Matrix  Tuple
   1.421 +		 (rows [this] (if (down? this)
   1.422 +				(list this)
   1.423 +				(map (comp up vector) this)))
   1.424 +
   1.425 +		 (cols [this] (if (up? this)
   1.426 +				(list this)
   1.427 +				(map (comp down vector) this))
   1.428 +		 ))
   1.429 +  
   1.430 +(defn matrix-multiply
   1.431 +  "Returns the matrix resulting from the matrix multiplication of the given arguments."
   1.432 +  ([A] (coerce-to-matrix A))
   1.433 +  ([A B] (incanter.core/mmult (coerce-to-matrix A) (coerce-to-matrix B)))
   1.434 +  ([M1 M2 & Ms] (reduce matrix-multiply (matrix-multiply M1 M2) Ms)))
   1.435 +
   1.436 +#+end_src
   1.437 +
   1.438 +
   1.439 +** Power Series
   1.440 +#+srcname power-series
   1.441 +#+begin_src clojure
   1.442 +(in-ns 'sicm.utils)
   1.443 +(use 'clojure.contrib.def)
   1.444 +
   1.445 +
   1.446 +
   1.447 +(defn series-fn
   1.448 +  "The function corresponding to the given power series."
   1.449 +  [series]
   1.450 +  (fn [x]
   1.451 +    (reduce +
   1.452 +	    (map-indexed  (fn[n x] (* (float (nth series n)) (float(java.lang.Math/pow (float x) n)) ))
   1.453 +			  (range 20)))))
   1.454 +
   1.455 +(deftype PowerSeries
   1.456 +  [coll]
   1.457 +  clojure.lang.Seqable
   1.458 +  (seq [this] (seq (.coll this)))
   1.459 +  
   1.460 +  clojure.lang.Indexed
   1.461 +  (nth [this n] (nth (.coll this) n 0))
   1.462 +  (nth [this n not-found] (nth (.coll this) n not-found))
   1.463 +
   1.464 +  ;; clojure.lang.IFn
   1.465 +  ;; (call [this] (throw(Exception.)))
   1.466 +  ;; (invoke [this & args] args
   1.467 +  ;; 	  (let [f 
   1.468 +  ;; )
   1.469 +  ;; (run [this] (throw(Exception.)))
   1.470 +  )
   1.471 +
   1.472 +(defn power-series
   1.473 +  "Returns a power series with the items of the coll as its
   1.474 +  coefficients. Trailing zeros are added to the end of coll."
   1.475 +  [coeffs]
   1.476 +  (PowerSeries. coeffs))
   1.477 +
   1.478 +(defn power-series-indexed
   1.479 +  "Returns a power series consisting of the result of mapping f to the non-negative integers."
   1.480 +  [f]
   1.481 +  (PowerSeries. (map f (range))))
   1.482 +
   1.483 +
   1.484 +(defn-memo nth-partial-sum
   1.485 +  ([series n]
   1.486 +     (if (zero? n) (first series)
   1.487 +	 (+ (nth series n)
   1.488 +	    (nth-partial-sum series (dec n))))))
   1.489 +
   1.490 +(defn partial-sums [series]
   1.491 +  (lazy-seq (map nth-partial-sum (range))))
   1.492 +
   1.493 +
   1.494 +
   1.495 +
   1.496 +(def cos-series
   1.497 +     (power-series-indexed
   1.498 +      (fn[n]
   1.499 +	(if (odd? n) 0
   1.500 +	    (/
   1.501 +	     (reduce *
   1.502 +		     (reduce * (repeat (/ n 2) -1))
   1.503 +		     (range 1 (inc n)))
   1.504 +	     )))))
   1.505 +
   1.506 +(def sin-series
   1.507 +     (power-series-indexed
   1.508 +      (fn[n]
   1.509 +	(if (even? n) 0
   1.510 +	    (/
   1.511 +	     (reduce *
   1.512 +		     (reduce * (repeat (/ (dec n) 2) -1))
   1.513 +		     (range 1 (inc n)))
   1.514 +	     )))))
   1.515 +
   1.516 +#+end_src
   1.517 +
   1.518 +
   1.519 +* Basic Utilities
   1.520 +
   1.521 +** Sequence manipulation
   1.522 +
   1.523 +#+srcname: seq-manipulation
   1.524 +#+begin_src clojure
   1.525 +(ns sicm.utils)
   1.526 +
   1.527 +(defn do-up
   1.528 +  "Apply f to each number from low to high, presumably for
   1.529 +  side-effects."
   1.530 +  [f low high]
   1.531 +  (doseq [i (range low high)] (f i)))
   1.532 +   
   1.533 +(defn do-down
   1.534 +  "Apply f to each number from high to low, presumably for
   1.535 +   side-effects."
   1.536 +  [f high low]
   1.537 +  (doseq [i (range high low -1)] (f i)))
   1.538 +
   1.539 +
   1.540 +(defn all-equal? [coll]
   1.541 +  (if (empty? (rest coll)) true
   1.542 +      (and (= (first coll) (second coll))
   1.543 +	   (recur (rest coll))))))
   1.544 +
   1.545 +(defn multiplier
   1.546 +  "Returns a function that 'multiplies' the members of a collection,
   1.547 +returning unit if called on an empty collection."
   1.548 +  [multiply unit]
   1.549 +  (fn [coll] ((partial reduce multiply unit) coll)))
   1.550 +
   1.551 +(defn divider
   1.552 +  "Returns a function that 'divides' the first element of a collection
   1.553 +by the 'product' of the rest of the collection."
   1.554 +  [divide multiply invert unit]
   1.555 +  (fn [coll]
   1.556 +    (apply
   1.557 +     (fn
   1.558 +       ([] unit)
   1.559 +       ([x] (invert x))
   1.560 +       ([x y] (divide x y))
   1.561 +       ([x y & zs] (divide x (reduce multiply y zs))))
   1.562 +     coll)))
   1.563 +
   1.564 +(defn left-circular-shift
   1.565 +  "Remove the first element of coll, adding it to the end of coll."
   1.566 +  [coll]
   1.567 +  (concat (rest coll) (take 1 coll)))
   1.568 +
   1.569 +(defn right-circular-shift
   1.570 +  "Remove the last element of coll, adding it to the front of coll."
   1.571 +  [coll]
   1.572 +  (cons (last coll) (butlast coll)))
   1.573 +#+end_src
   1.574 +
   1.575 +
   1.576 +
   1.577 +
   1.578 +** Ranges, Airity and Function Composition
   1.579 +#+srcname: arity
   1.580 +#+begin_src clojure
   1.581 +(in-ns 'sicm.utils)
   1.582 +(def infinity Double/POSITIVE_INFINITY)
   1.583 +(defn infinite? [x] (Double/isInfinite x))
   1.584 +(def finite? (comp not infinite?)) 
   1.585 +
   1.586 +(defn arity-min
   1.587 +  "Returns the smallest number of arguments f can take."
   1.588 +  [f]
   1.589 +  (apply
   1.590 +   min
   1.591 +   (map (comp alength #(.getParameterTypes %))
   1.592 +	(filter (comp (partial = "invoke") #(.getName %))
   1.593 +		(.getDeclaredMethods (class f))))))
   1.594 +  
   1.595 +(defn arity-max
   1.596 +  "Returns the largest number of arguments f can take, possibly
   1.597 +  Infinity."
   1.598 +  [f]
   1.599 +  (let [methods (.getDeclaredMethods (class f))]
   1.600 +    (if (not-any? (partial = "doInvoke") (map #(.getName %) methods))
   1.601 +      (apply max
   1.602 +	     (map (comp alength #(.getParameterTypes %))
   1.603 +		  (filter (comp (partial = "invoke") #(.getName %))  methods)))
   1.604 +      infinity)))
   1.605 +
   1.606 +
   1.607 +(def ^{:arglists '([f])
   1.608 +       :doc "Returns a two-element list containing the minimum and
   1.609 +     maximum number of args that f can take."}
   1.610 +     arity-interval
   1.611 +     (juxt arity-min arity-max))
   1.612 +
   1.613 +
   1.614 +
   1.615 +;; --- intervals
   1.616 +
   1.617 +(defn intersect
   1.618 +  "Returns the interval of overlap between interval-1 and interval-2"
   1.619 +  [interval-1 interval-2]
   1.620 +  (if (or (empty? interval-1) (empty? interval-2)) []
   1.621 +      (let [left (max (first interval-1) (first interval-2))
   1.622 +	    right (min (second interval-1) (second interval-2))]
   1.623 +	(if (> left right) []
   1.624 +	    [left right]))))
   1.625 +
   1.626 +(defn same-endpoints?
   1.627 +  "Returns true if the left endpoint is the same as the right
   1.628 +  endpoint."
   1.629 +  [interval]
   1.630 +  (= (first interval) (second interval)))
   1.631 +
   1.632 +(defn naturals?
   1.633 +  "Returns true if the left endpoint is 0 and the right endpoint is
   1.634 +infinite."
   1.635 +  [interval]
   1.636 +  (and (zero? (first interval))
   1.637 +       (infinite? (second interval))))
   1.638 + 
   1.639 +
   1.640 +(defn fan-in
   1.641 +  "Returns a function that pipes its input to each of the gs, then
   1.642 +  applies f to the list of results. Consequently, f must be able to
   1.643 +  take a number of arguments equal to the number of gs."
   1.644 +  [f & gs]
   1.645 +  (fn [& args]
   1.646 +    (apply f (apply (apply juxt gs) args))))
   1.647 +
   1.648 +(defn fan-in
   1.649 +  "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."
   1.650 +  [f & gs]
   1.651 +  (comp (partial apply f) (apply juxt gs)))
   1.652 +
   1.653 +
   1.654 +
   1.655 +(defmacro airty-blah-sad [f n more?]
   1.656 +  (let [syms (vec (map (comp gensym (partial str "x")) (range n)))
   1.657 +	 optional (gensym "xs")]
   1.658 +    (if more?
   1.659 +      `(fn ~(conj syms '& optional)
   1.660 +	 (apply ~f ~@syms ~optional))
   1.661 +      `(fn ~syms (~f ~@syms)))))
   1.662 +
   1.663 +(defmacro airt-whaa* [f n more?]
   1.664 +  `(airty-blah-sad ~f ~n ~more?))
   1.665 +
   1.666 +
   1.667 +
   1.668 +
   1.669 +(defn fan-in*
   1.670 +  "Returns a function that pipes its input to each of the gs, then
   1.671 +  applies f to the list of results. Unlike fan-in, fan-in* strictly
   1.672 +  enforces arity: it will fail if the gs do not have compatible
   1.673 +  arities."
   1.674 +  [f & gs]
   1.675 +  (let [arity-in (reduce intersect (map arity-interval gs))
   1.676 +	left (first arity-in)
   1.677 +	right (second arity-in)
   1.678 +	composite (fan-in f gs)
   1.679 +	]
   1.680 +    (cond
   1.681 +     (empty? arity-in)
   1.682 +     (throw (Exception. "Cannot compose functions with incompatible arities."))
   1.683 +
   1.684 +     (not
   1.685 +      (or (= left right)
   1.686 +	  (and (finite? left)
   1.687 +	       (= right infinity))))
   1.688 +	 
   1.689 +     (throw (Exception.
   1.690 +	     "Compose can only handle arities of the form [n n] or [n infinity]"))
   1.691 +     :else
   1.692 +     (airty-blah-sad composite left (= right infinity)))))
   1.693 +    
   1.694 +
   1.695 +
   1.696 +(defn compose-n "Compose any number of functions together."
   1.697 +  ([] identity)
   1.698 +  ([f] f)
   1.699 +  ([f & fs]
   1.700 +  (let [fns (cons f fs)]
   1.701 +    (compose-bin (reduce fan-in (butlast fs)) (last fs))))
   1.702 +)
   1.703 +
   1.704 +
   1.705 +
   1.706 +
   1.707 +	   
   1.708 +
   1.709 +(defn iterated
   1.710 +  ([f n id] (reduce comp id (repeat n f)))
   1.711 +  ([f n] (reduce comp identity (repeat n f))))
   1.712 +
   1.713 +(defn iterate-until-stable
   1.714 +  "Repeatedly applies f to x, returning the first result that is close
   1.715 +enough to its predecessor."
   1.716 +  [f close-enough? x]
   1.717 +  (second (swank.util/find-first
   1.718 +	   (partial apply close-enough?)
   1.719 +	   (partition 2 1 (iterate f x)))))
   1.720 +
   1.721 +(defn lexical< [x y]
   1.722 +  (neg? (compare (str x) (str y))))
   1.723 +
   1.724 +
   1.725 +;; do-up
   1.726 +;; do-down
   1.727 +(def make-pairwise-test comparator)
   1.728 +;;all-equal?
   1.729 +(def accumulation multiplier)
   1.730 +(def inverse-accumulation divider)
   1.731 +;;left-circular-shift
   1.732 +;;right-circular-shift
   1.733 +(def exactly-n? same-endpoints?)
   1.734 +(def any-number? naturals?)
   1.735 +;; TODO compose
   1.736 +;; TODO compose-n
   1.737 +;; identity
   1.738 +(def compose-2 fan-in)
   1.739 +(def compose-bin fan-in*)
   1.740 +(def any? (constantly true))
   1.741 +(def none? (constantly false))
   1.742 +(def constant constantly)
   1.743 +(def joint-arity intersect)
   1.744 +(def a-reduce reduce)
   1.745 +;; filter
   1.746 +(def make-map (partial partial map) )
   1.747 +(def bracket juxt)
   1.748 +;; TODO apply-to-all
   1.749 +;; TODO nary-combine
   1.750 +;; TODO binary-combine
   1.751 +;; TODO unary-combine
   1.752 +;; iterated
   1.753 +;; iterate-until-stable
   1.754 +(def make-function-of-vector (partial partial map))
   1.755 +(def make-function-of-arguments (fn [f] (fn [& args] (f args))))
   1.756 +(def alphaless lexical<)
   1.757 +
   1.758 +#+end_src
   1.759 +
   1.760 +
   1.761 +
   1.762 +
   1.763 +
   1.764 +
   1.765 +
   1.766 +: 
   1.767 +
   1.768 +*  Numerical Methods
   1.769 +#+srcname: numerical-methods
   1.770 +#+begin_src clojure
   1.771 +(in-ns 'sicm.utils)
   1.772 +(import java.lang.Math)
   1.773 +(use 'clojure.contrib.def)
   1.774 +
   1.775 +;; ---- USEFUL CONSTANTS
   1.776 +
   1.777 +(defn machine-epsilon
   1.778 +  "Smallest value representable on your machine, as determined by
   1.779 +successively dividing a number in half until consecutive results are
   1.780 +indistinguishable."
   1.781 +  []
   1.782 +  (ffirst
   1.783 +   (drop-while
   1.784 +    (comp not zero? second)
   1.785 +    (partition 2 1
   1.786 +	       (iterate (partial * 0.5) 1)))))
   1.787 +
   1.788 +
   1.789 +(def pi (Math/PI))
   1.790 +(def two-pi (* 2 pi))
   1.791 +
   1.792 +(def eulers-gamma 0.5772156649015328606065)
   1.793 +
   1.794 +(def phi (/ (inc (Math/sqrt 5)) 2))
   1.795 +
   1.796 +(def ln2 (Math/log 2))
   1.797 +(def ln10 (Math/log 10))
   1.798 +(def exp10 #(Math/pow 10 %))
   1.799 +(def exp2 #(Math/pow 2 %))
   1.800 +
   1.801 +
   1.802 +;;
   1.803 +
   1.804 +;; ---- ANGLES AND TRIGONOMETRY
   1.805 +
   1.806 +(defn angle-restrictor
   1.807 +  "Returns a function that ensures that angles lie in the specified interval of length two-pi."
   1.808 +  [max-angle]
   1.809 +  (let [min-angle (- max-angle two-pi)]
   1.810 +    (fn [x]
   1.811 +      (if (and
   1.812 +	   (<= min-angle x)
   1.813 +	   (< x max-angle))
   1.814 +	x
   1.815 +	(let [corrected-x (- x (* two-pi (Math/floor (/ x two-pi))))]
   1.816 +	  (if (< corrected-x max-angle)
   1.817 +	    corrected-x
   1.818 +	    (- corrected-x two-pi)))))))
   1.819 +
   1.820 +(defn angle-restrict-pi
   1.821 +  "Coerces angles to lie in the interval from -pi to pi."
   1.822 +  [angle]
   1.823 +  ((angle-restrictor pi) angle))
   1.824 +
   1.825 +(defn angle-restrict-two-pi
   1.826 +  "Coerces angles to lie in the interval from zero to two-pi"
   1.827 +  [angle]
   1.828 +  ((angle-restrictor two-pi) angle))
   1.829 +
   1.830 +
   1.831 +
   1.832 +
   1.833 +(defn invert [x] (/ x))
   1.834 +(defn negate [x] (- x))
   1.835 +
   1.836 +(defn exp [x] (Math/exp x))
   1.837 +
   1.838 +(defn sin [x] (Math/sin x))
   1.839 +(defn cos [x] (Math/cos x))
   1.840 +(defn tan [x] (Math/tan x))
   1.841 +
   1.842 +(def sec (comp invert cos))
   1.843 +(def csc (comp invert sin))
   1.844 +
   1.845 +(defn sinh [x] (Math/sinh x))
   1.846 +(defn cosh [x] (Math/cosh x))
   1.847 +(defn tanh [x] (Math/tanh x))
   1.848 +
   1.849 +(def sech (comp invert cosh))
   1.850 +(def csch (comp invert sinh))
   1.851 +
   1.852 +
   1.853 +;; ------------
   1.854 +
   1.855 +(defn factorial
   1.856 +  "Computes the factorial of the nonnegative integer n."
   1.857 +  [n]
   1.858 +  (if (neg? n)
   1.859 +    (throw (Exception. "Cannot compute the factorial of a negative number."))
   1.860 +    (reduce * 1 (range 1 (inc n)))))
   1.861 +
   1.862 +(defn exact-quotient [n d] (/ n d))
   1.863 +
   1.864 +(defn binomial-coefficient
   1.865 +  "Computes the number of different ways to choose m elements from n."
   1.866 +  [n m]
   1.867 +  (assert (<= 0 m n))
   1.868 +  (let [difference (- n m)]
   1.869 +    (exact-quotient
   1.870 +     (reduce * (range n (max difference m) -1 ))
   1.871 +     (factorial (min difference m)))))
   1.872 +
   1.873 +(defn-memo stirling-1
   1.874 +  "Stirling numbers of the first kind: the number of permutations of n
   1.875 +  elements with exactly m permutation cycles. "
   1.876 +  [n k]
   1.877 +  ;(assert (<= 1 k n))
   1.878 +  (if (zero? n)
   1.879 +    (if (zero? k) 1 0)
   1.880 +    (+ (stirling-1 (dec n) (dec k))
   1.881 +       (* (dec n) (stirling-1 (dec n) k)))))
   1.882 +
   1.883 +(defn-memo stirling-2 ;;count-partitions
   1.884 +  "Stirling numbers of the second kind: the number of ways to partition a set of n elements into k subsets."
   1.885 +  [n k]
   1.886 +  (cond
   1.887 +   (= k 1) 1
   1.888 +   (= k n) 1
   1.889 +   :else (+ (stirling-2 (dec n) (dec k))
   1.890 +	    (* k (stirling-2 (dec n) k)))))
   1.891 +
   1.892 +(defn harmonic-number [n]
   1.893 +  (/ (stirling-1 (inc n) 2)
   1.894 +     (factorial n)))
   1.895 +
   1.896 +
   1.897 +(defn sum
   1.898 +  [f low high]
   1.899 +  (reduce + (map f (range low (inc high)))))
   1.900 +
   1.901 +#+end_src
   1.902 +
   1.903 +
   1.904 +
   1.905 +
   1.906 +
   1.907 +
   1.908 +
   1.909 +
   1.910 +
   1.911 +
   1.912 +
   1.913 +
   1.914 +* Differentiation
   1.915 +
   1.916 +We compute derivatives by passing special *differential objects* $[x,
   1.917 +dx]$ through functions. Roughly speaking, applying a function $f$ to a
   1.918 +differential object \([x, dx]\) should produce a new differential
   1.919 +object $[f(x),\,Df(x)\cdot dx]$.
   1.920 +
   1.921 +\([x,\,dx]\xrightarrow{\quad f \quad}[f(x),\,Df(x)\cdot dx]\)
   1.922 +Notice that you can obtain the derivative of $f$ from this
   1.923 +differential object, as it is the coefficient of the $dx$ term. Also,
   1.924 +as you apply successive functions using this rule, you get the
   1.925 +chain-rule answer you expect:
   1.926 +
   1.927 +\([f(x),\,Df(x)\cdot dx]\xrightarrow{\quad g\quad} [gf(x),\,
   1.928 +Dgf(x)\cdot Df(x) \cdot dx ]\)
   1.929 +
   1.930 +In order to generalize to multiple variables and multiple derivatives,
   1.931 +we use a *power series of differentials*, a sortred infinite sequence which
   1.932 +contains all terms like $dx\cdot dy$, $dx^2\cdot dy$, etc.
   1.933 +
   1.934 +
   1.935 +#+srcname:differential
   1.936 +#+begin_src clojure
   1.937 +(in-ns 'sicm.utils)
   1.938 +(use 'clojure.contrib.combinatorics)
   1.939 +(use 'clojure.contrib.generic.arithmetic)
   1.940 +
   1.941 +(defprotocol DifferentialTerm
   1.942 +  "Protocol for an infinitesimal quantity."
   1.943 +  (coefficient [this])
   1.944 +  (partials [this]))
   1.945 +
   1.946 +(extend-protocol DifferentialTerm
   1.947 +  java.lang.Number
   1.948 +  (coefficient [this] this)
   1.949 +  (partials [this] []))
   1.950 +
   1.951 +(deftype DifferentialSeq
   1.952 +  [terms]
   1.953 +  clojure.lang.IPersistentCollection
   1.954 +  (cons [this x]
   1.955 +	(throw (Exception. x))
   1.956 +	(DifferentialSeq. (cons x terms)))
   1.957 + ;; (seq [this] (seq terms))
   1.958 +  (count [this] (count terms))
   1.959 +  (empty [this] (empty? terms)))
   1.960 +
   1.961 +
   1.962 +
   1.963 +
   1.964 +
   1.965 +
   1.966 +
   1.967 +
   1.968 +
   1.969 +
   1.970 +
   1.971 +(defn coerce-to-differential-seq [x]
   1.972 +  (cond
   1.973 +   (= (type x) DifferentialSeq) x
   1.974 +   (satisfies? DifferentialTerm x) (DifferentialSeq. x)))
   1.975 +
   1.976 +
   1.977 +(defn differential-term
   1.978 +  "Returns a differential term with the given coefficient and
   1.979 +  partials. Coefficient may be any arithmetic object; partials must
   1.980 +  be a list of non-negative integers."
   1.981 +  [coefficient partials]
   1.982 +  (reify DifferentialTerm
   1.983 +	 (partials [_] (set partials))
   1.984 +	 (coefficient [_] coefficient)))
   1.985 +
   1.986 +
   1.987 +
   1.988 +(defn differential-seq*
   1.989 +  ([coefficient partials]
   1.990 +     (DifferentialSeq. [(differential-term coefficient partials)]))
   1.991 +  ([coefficient partials & cps]
   1.992 +     (if cps
   1.993 +       
   1.994 +	
   1.995 +
   1.996 +
   1.997 +(defn differential-seq
   1.998 +  "Constructs a sequence of differential terms from a numerical
   1.999 +coefficient and a list of keys for variables. If no coefficient is supplied, uses 1."
  1.1000 +  ([variables] (differential-seq 1 variables))
  1.1001 +  ([coefficient variables & cvs]
  1.1002 +     (if (number? coefficient)
  1.1003 +       (conj (assoc {} (apply sorted-set variables) coefficient)
  1.1004 +	     (if (empty? cvs)
  1.1005 +	       nil
  1.1006 +	       (apply differential-seq cvs)))
  1.1007 +       (apply differential-seq 1 coefficient 1 variables cvs)  
  1.1008 +       )))
  1.1009 +
  1.1010 +
  1.1011 +(defn differential-add
  1.1012 +  "Add two differential sequences by combining like terms."
  1.1013 +  [dseq1 dseq2]
  1.1014 +  (merge-with + dseq1 dseq2))
  1.1015 +
  1.1016 +(defn differential-multiply
  1.1017 +  "Multiply two differential sequences. The square of any differential variable is zero since differential variables are infinitesimally small."
  1.1018 +  [dseq1 dseq2]
  1.1019 +  (reduce
  1.1020 +   (fn [m [[vars1 coeff1] [vars2 coeff2]]]
  1.1021 +     (if (empty? (clojure.set/intersection vars1 vars2))
  1.1022 +       (assoc m (clojure.set/union vars1 vars2) (* coeff1 coeff2))
  1.1023 +       m))
  1.1024 +   {}
  1.1025 +  (clojure.contrib.combinatorics/cartesian-product
  1.1026 +   dseq1
  1.1027 +   dseq2)))
  1.1028 +
  1.1029 +
  1.1030 +(defn big-part
  1.1031 +  "Returns the part of the differential sequence that is finite,
  1.1032 +  i.e. not infinitely small."
  1.1033 +  [dseq]
  1.1034 +  (let
  1.1035 +      [keys (sort-by count (keys dseq))
  1.1036 +       smallest-var (last(last keys))]
  1.1037 +    (apply hash-map
  1.1038 +	   (reduce concat
  1.1039 +		   (remove (comp smallest-var first) dseq)))))
  1.1040 +
  1.1041 +(defn small-part
  1.1042 +  "Returns the part of the differential sequence that is
  1.1043 +  infinitesimal."
  1.1044 +  [dseq]
  1.1045 +  (let
  1.1046 +      [keys (sort-by count (keys dseq))
  1.1047 +       smallest-var (last(last keys))]
  1.1048 +    (apply hash-map
  1.1049 +	   (reduce concat
  1.1050 +		   (filter (comp smallest-var first) dseq)))))
  1.1051 +
  1.1052 +(defn big-part
  1.1053 +  "Returns the 'finite' part of the differential sequence."
  1.1054 +  [dseq]
  1.1055 +  (let
  1.1056 +      [keys (sort-by count (keys dseq))
  1.1057 +       smallest-var (last (last keys))]
  1.1058 +    (apply hash-map
  1.1059 +    (reduce concat
  1.1060 +           (filter (remove smallest-var first) dseq)
  1.1061 +	   ))))
  1.1062 +
  1.1063 +
  1.1064 +(defn small-part
  1.1065 +  "Returns the 'infinitesimal' part of the differential sequence."
  1.1066 +  [dseq]
  1.1067 +    (let
  1.1068 +      [keys (sort-by count (keys dseq))
  1.1069 +       smallest-var (last (last keys))]
  1.1070 +
  1.1071 +    (apply hash-map
  1.1072 +    (reduce concat
  1.1073 +	    (filter (comp smallest-var first) dseq)
  1.1074 +	    ))))
  1.1075 +
  1.1076 +(defn linear-approximator
  1.1077 +  "Returns the function that corresponds to unary-fn but which accepts and returns differential objects."
  1.1078 +  ([unary-f dfdx]
  1.1079 +     (fn F [dseq]
  1.1080 +       (differential-add
  1.1081 +	(unary-f (big-part dseq))
  1.1082 +	(differential-multiply
  1.1083 +	 (dfdx (big-part dseq))
  1.1084 +	 (small-part dseq))))))
  1.1085 +
  1.1086 +
  1.1087 +
  1.1088 +
  1.1089 +
  1.1090 +
  1.1091 +;; ;; A differential term consists of a numerical coefficient and a
  1.1092 +;; ;; sorted  
  1.1093 +;; (defrecord DifferentialTerm [coefficient variables])
  1.1094 +;; (defmethod print-method DifferentialTerm
  1.1095 +;;    [o w]
  1.1096 +;;     (print-simple
  1.1097 +;;      (apply str (.coefficient o)(map (comp (partial str "d") name) (.variables o)))
  1.1098 +;;      w))
  1.1099 +
  1.1100 +
  1.1101 +;; (defn differential-seq
  1.1102 +;;   "Constructs a sequence of differential terms from a numerical
  1.1103 +;; coefficient and a list of keywords for variables. If no coefficient is
  1.1104 +;; supplied, uses 1."
  1.1105 +;;   ([variables] (differential-seq 1 variables))
  1.1106 +;;   ([coefficient variables]
  1.1107 +;;      (list
  1.1108 +;;       (DifferentialTerm. coefficient (apply sorted-set variables))))
  1.1109 +;;   ([coefficient variables & cvs]
  1.1110 +;;      (sort-by
  1.1111 +;;       #(vec(.variables %))
  1.1112 +;;       (concat (differential-seq coefficient variables) (apply differential-seq cvs)))))
  1.1113 +
  1.1114 +#+end_src
  1.1115 +
  1.1116 +
  1.1117 +
  1.1118 +
  1.1119 +
  1.1120 +
  1.1121 +
  1.1122 +
  1.1123 +
  1.1124 +
  1.1125 +
  1.1126 +
  1.1127 +
  1.1128 +* Symbolic Manipulation
  1.1129 +
  1.1130 +#+srcname:symbolic
  1.1131 +#+begin_src clojure
  1.1132 +(in-ns 'sicm.utils)
  1.1133 +
  1.1134 +(deftype Symbolic [type expression]
  1.1135 +  Object
  1.1136 +  (equals [this that]
  1.1137 +	  (cond
  1.1138 +	   (= (.expression this) (.expression that)) true
  1.1139 +	   :else
  1.1140 +	   (Symbolic.
  1.1141 +	    java.lang.Boolean
  1.1142 +	    (list '= (.expression this) (.expression that)))
  1.1143 +	  )))
  1.1144 +
  1.1145 +
  1.1146 +
  1.1147 +
  1.1148 +(deftype AbstractSet [glyph membership-test])
  1.1149 +(defn member? [abstract-set x]
  1.1150 +  ((.membership-test abstract-set) x))
  1.1151 +
  1.1152 +;; ------------ Some important AbstractSets
  1.1153 +
  1.1154 +
  1.1155 +(def Real
  1.1156 +     (AbstractSet.
  1.1157 +      'R
  1.1158 +      (fn[x](number? x))))
  1.1159 +
  1.1160 +
  1.1161 +;; ------------ Create new AbstractSets from existing ones
  1.1162 +
  1.1163 +(defn abstract-product
  1.1164 +  "Gives the cartesian product of abstract sets."
  1.1165 +  ([sets]
  1.1166 +       (if (= (count sets) 1) (first sets)
  1.1167 +	   (AbstractSet.
  1.1168 +	    (symbol
  1.1169 +	     (apply str
  1.1170 +		    (interpose 'x (map #(.glyph %) sets))))
  1.1171 +	   (fn [x]
  1.1172 +	     (and
  1.1173 +	      (coll? x)
  1.1174 +	      (= (count sets) (count x))
  1.1175 +	      (reduce #(and %1 %2)
  1.1176 +		      true
  1.1177 +		      (map #(member? %1 %2) sets x)))))))
  1.1178 +  ([abstract-set n]
  1.1179 +     (abstract-product (repeat n abstract-set))))
  1.1180 +
  1.1181 +
  1.1182 +	      
  1.1183 +(defn abstract-up
  1.1184 +  "Returns the abstract set of all up-tuples whose items belong to the
  1.1185 +  corresponding abstract sets in coll."
  1.1186 +  ([coll]
  1.1187 +     (AbstractSet.
  1.1188 +      (symbol (str "u["
  1.1189 +		   (apply str
  1.1190 +			  (interpose " "
  1.1191 +				     (map #(.glyph %) coll)))
  1.1192 +		   "]"))
  1.1193 +      (fn [x]
  1.1194 +	(and
  1.1195 +	 (satisfies? Spinning x)
  1.1196 +	 (up? x)
  1.1197 +	 (= (count coll) (count x))
  1.1198 +	 (reduce
  1.1199 +	  #(and %1 %2)
  1.1200 +	  true
  1.1201 +	  (map #(member? %1 %2) coll x))))))
  1.1202 +  ([abstract-set n]
  1.1203 +     (abstract-up (repeat n abstract-set))))
  1.1204 +
  1.1205 +
  1.1206 +(defn abstract-down
  1.1207 +  "Returns the abstract set of all down-tuples whose items belong to the
  1.1208 +  corresponding abstract sets in coll."
  1.1209 +  ([coll]
  1.1210 +     (AbstractSet.
  1.1211 +      (symbol (str "d["
  1.1212 +		   (apply str
  1.1213 +			  (interpose " "
  1.1214 +				     (map #(.glyph %) coll)))
  1.1215 +		   "]"))
  1.1216 +      (fn [x]
  1.1217 +	(and
  1.1218 +	 (satisfies? Spinning x)
  1.1219 +	 (down? x)
  1.1220 +	 (= (count coll) (count x))
  1.1221 +	 (reduce
  1.1222 +	  #(and %1 %2)
  1.1223 +	  true
  1.1224 +	  (map #(member? %1 %2) coll x))))))
  1.1225 +  ([abstract-set n]
  1.1226 +     (abstract-down (repeat n abstract-set))))
  1.1227 +
  1.1228 +	      
  1.1229 +
  1.1230 +
  1.1231 +
  1.1232 +					;-------ABSTRACT FUNCTIONS
  1.1233 +(defrecord AbstractFn
  1.1234 +  [#^AbstractSet domain #^AbstractSet codomain])
  1.1235 +
  1.1236 +
  1.1237 +(defmethod print-method AbstractFn
  1.1238 +  [o w]
  1.1239 +  (print-simple
  1.1240 +   (str
  1.1241 +    "f:"
  1.1242 +   (.glyph (:domain o))
  1.1243 +   "-->"
  1.1244 +   (.glyph (:codomain o))) w))
  1.1245 +#+end_src
  1.1246 +
  1.1247 +
  1.1248 +* COMMENT 
  1.1249 +#+begin_src clojure :tangle utils.clj
  1.1250 +(ns sicm.utils)
  1.1251 +
  1.1252 +					;***** GENERIC ARITHMETIC
  1.1253 +<<generic-arithmetic>>
  1.1254 +
  1.1255 +
  1.1256 +					;***** TUPLES AND MATRICES
  1.1257 +<<tuples>>
  1.1258 +<<tuples-2>>
  1.1259 +					;***** MATRICES
  1.1260 +<<matrices>>
  1.1261 +
  1.1262 +#+end_src
  1.1263 +
  1.1264 +
  1.1265 +