aurellem

Building a Classical Mechanics Library in Clojure

Table of Contents

1 Generic Arithmetic

(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))))
           

2 Useful Data Types

2.1 Complex Numbers

(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))


2.2 Tuples and Tensors

A tuple is a vector which is spinable—it can be either spin up or spin down. (Covariant, contravariant; dual vectors)

(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))


2.2.1 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.

(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))))

2.2.2 Matrices

(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)))

2.3 Power Series

(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)))
             )))))

3 Basic Utilities

3.1 Sequence manipulation

(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)))

3.2 Ranges, Airity and Function Composition

(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<)

4 Numerical Methods

(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)))))

5 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.

(in-ns 'sicm.utils)

(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 'finite' part of the differential sequence."
  [dseq]
  (let
      [keys (sort-by count (keys dseq))
       pivot-key (first (last keys))]

    (apply hash-map
    (reduce concat
            (filter (comp pivot-key first) dseq)
            ))))


(defn small-part
  "Returns the 'infinitesimal' part of the differential sequence."
  [dseq]
    (let
      [keys (sort-by count (keys dseq))
       pivot-key (first (last keys))]

    (apply hash-map
    (reduce concat
            (remove (comp pivot-key first) 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)))))

6 Symbolic Manipulation

(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))

Date: 2011-08-09 18:41:37 EDT

Author: Robert McIntyre & Dylan Holmes

Org version 7.6 with Emacs version 23

Validate XHTML 1.0