dynamics/src/math.clj

   1 (ns src.math)
   2 ;(alias 'cl 'clojure.core)
   3
   4 (defn neg-index [cnt x]
   5   (if (neg? x) (+ cnt x) x))
   6
   7 (defn swap [i j s]
   8   (let [[i j] (map (partial neg-index (count s)) [i j])]
   9     (assoc (assoc (vec s) i (nth s j))
  10            j (nth s i)
  11            )))
  12
  13 (defn search-non-zero [mat pivot]
  14   (let [width (count (first mat))
  15         height (count mat)]
  16     (first (for [j (range pivot (dec width))
  17                  i (range pivot height)
  18                  :when (not= 0 (nth (nth mat i) j))]
  19                 [i j]))))
  20
  21 (defn sweep1 [mat pivot]
  22   (let [denom (nth (nth mat pivot) pivot)
  23         p-row (nth mat pivot)]
  24     (map (fn [i row]
  25            (if (= i pivot)
  26              (map #(/ % denom) row)
  27              (let [multiplier (/ (nth row pivot) denom)]
  28                (map #(- %1 %2) row (map #(* % multiplier) p-row))
  29                )))
  30          (iterate inc 0) mat)))
  31
  32 (defn lin-solve [mat]
  33   (let [width  (count (first mat))
  34         height (count mat)]
  35     (loop [pivot 0 m (vec (map vec mat)) swap-acc ()]
  36       (if (>= pivot (max (dec width) height))
  37         [m swap-acc]
  38         (let [[col row] (search-non-zero mat pivot)]
  39           (if (nil? col)
  40             m
  41             (recur (inc pivot)
  42                    (sweep1 (swap pivot row
  43                                  (map (partial swap pivot col) m))
  44                            pivot)
  45                    (cons [pivot col] swap-acc)
  46                    )))))))
  47
  48 ;
  49 ;  __ 0      |\  2        3
  50 ; |__>-------| >o--------+
  51 ;           1|/          |
  52 ;                        |
  53 ;                        |          __
  54 ;                        +---------<__|
  55 ;                       4         5
  56 ; nodes:
  57 ;  [0 in] [1 mid] [2 mid] [3 mid] [4 mid] [5 out]
  58 ; edges:
  59 ;  [hor 0 1] [not1 1 2] [hor 2 3] [ver 3 4] [hor 4 5]
  60 ;
  61 ; [x0 y0 x1 y1 x2 y2 x3 y3 x4 y4 x5 y5        -const]
  62 ; [ 1  0  0  0  0  0  0  0  0  0  0  0 (-  inport-x)] ; [0 in]
  63 ; [ 0  0  0  0  0  0  0  0  0  0  1  0 (- outport-x)] ; [5 out]
  64 ; [ 0  1  0 -1  0  0  0  0  0  0  0  0             0] ; [hor 0 1]
  65 ; [ 0  0 -1  0  1  0  0  0  0  0  0  0    not1-width] ; [not1 1 2]
  66 ; [ 0  0  0  1  0 -1  0  0  0  0  0  0             0] ; [not1 1 2]
  67 ; [ 0  0  0  0  0  1  0 -1  0  0  0  0             0] ; [hor 2 3]
  68 ; [ 0  0  0  0  0  0  1  0 -1  0  0  0             0] ; [ver 3 4]
  69 ; [ 0  0  0  0  0  0  0  0  0  1  0 -1             0] ; [hor 4 5]
  70
  71 (defn one-hot [tuple idx]
  72   (cons (first tuple)
  73         (assoc (vec (take (count (rest tuple)) (repeat 0)))
  74                idx 1)))
  75
  76 (defn all-one [tuple]
  77   (cons (first tuple)
  78         (map (fn [x] (if (seq? x) (recur x) 1))
  79              (rest tuple))))
  80
  81 (defn trans-t [tuple] ; transposes tuple
  82   (cons ({'up 'down, 'down 'up} (first tuple))
  83         (rest tuple)))
  84
  85 (def f-table
  86   {'+ +, '- -, '* *, '/ /,
  87    'ref (fn [tuple idx]
  88           (nth tuple (inc idx)))
  89    'one-hot one-hot
  90    'all-one all-one
  91    'trans-t trans-t
  92    })
  93
  94 (defn neval [expr arg env]
  95   (cond (isa? (class expr) Number) expr
  96         (= expr 'arg) arg
  97         (symbol? expr) (env expr)
  98         :else (apply (f-table (first expr))
  99                      (map #(neval % arg env) (rest expr))
 100                      )))
 101
 102 (defn m*v [m v]
 103   (map #(apply + (map * % v)) m))
 104
 105 (defn s*m [s m]
 106   (map (fn [v] (map #(* s %) v))
 107        m))
 108
 109 (defn m-m [m0 m1]
 110   (map #(map - %1 %2)
 111        m0 m1))
 112
 113 (defn i-mat [size]
 114   (loop [i size acc []]
 115     (if (<= i 0)
 116       acc
 117       (recur (dec i)
 118              (cons (cons 1 (take (count acc) (repeat 0)))
 119                    (map #(cons 0 %) acc)
 120                    )))))
 121
 122 (def *dxi* 0.001)
 123
 124 (defn jaco [phys xis env]
 125   (let [xis (vec xis)
 126         xis+ (map #(assoc xis %1 (+ (nth xis %1) (/ *dxi* 2)))
 127                   (range (count xis)))
 128         xis- (map #(assoc xis %1 (- (nth xis %1) (/ *dxi* 2)))
 129                   (range (count xis)))]
 130     (map (fn [phy] (map #(/ (- (neval phy (cons '_ %1) env)
 131                                (neval phy (cons '_ %2) env))
 132                             *dxi*)
 133                          xis+ xis-))
 134          phys)))
 135
 136 ; dx_i/dt = f_i(x)
 137 ; x^i(t+h) = x^i(t) + f^i(x(t+h))h ; implicit Eular method
 138 ; x^i(t+h) = x^i(t) + f^i(x(t) + x(t+h) - x(t))h
 139 ; x^i(t+h) = x^i(t) + f^i(x(t))h + df^i(x(t))/dx^j * (x^j(t+h) - x^j(t)) * h
 140 ; x~i(t+h) - h*df^i(x(t))/dx^j*x^j(t+h)
 141 ;  = x^i(t) + f^i(x(t))h - h*df^i(x(t))/dx^j*x^j(t)
 142 (def *h* 0.02)
 143
 144 (defn next-iter [phys xis env]
 145   (let [j (jaco phys xis env)
 146         lhs (m-m (i-mat (count xis))
 147                  (s*m *h* j))
 148         rhs (map #(+ %1 %2 (- %3))
 149                  xis
 150                  (map #(* *h* (neval % (cons '_ xis) env)) phys)
 151                  (m*v (s*m *h* j) xis))
 152         [solved varswap] (lin-solve (map #(conj %1 %2)
 153                                          (map vec lhs) rhs))]
 154     (reduce #(swap (nth %2 0) (nth %2 1) %1)
 155             (map #(nth % (count xis)) solved)
 156             varswap)))
 157
 158 ; symbol type
 159 (defn stype [x]
 160   (cond (isa? (class x) Number) 'num
 161         (symbol? x)             'symb
 162         :else                   (first x)
 163         ))
 164
 165 (defn stype-order [x]
 166   (let [order {'num  0
 167                'symb 1
 168                'up   2
 169                'down 3
 170                '+    4}]
 171     (get order (stype x))
 172     ))
 173
 174 ; symbolic addtion 2 terms
 175 (defn s+2 [x0 x1]
 176   (let [[t1 t2] [(stype x0) (stype x1)]]
 177     (cond (= [t1 t2] '[+    +   ]) (concat ['+] (rest x0) (rest x1))
 178           (= [t1 t2] '[symb +   ]) (concat ['+ x0] (rest x1))
 179           (=     t2         +    ) ['+ (s+2 x0 (first x1)) (rest x1)]
 180           (= [t1 t2] '[num  num ]) (+ x0 x1)
 181           (= [t1 t2] '[up   up  ]) (cons 'up (map + (rest x0) (rest x1)))
 182           (= [t1 t2] '[down down]) (cons 'down (map + (rest x0) (rest x1)))
 183           )))
 184
 185 ; symbolic add
 186 (defn s+ [& args]
 187   (letfn [(rec [x & xs]
 188             (if (empty? xs)
 189               x
 190               (s+2 x (apply rec xs))
 191               ))]
 192     (apply rec (sort-by stype-order args))
 193     ))
 194
 195 (defn n*2 [x0 x1]
 196   (let [[t0 t1] (map stype [x0 x1])]
 197     (cond (= [t0 t1] '[num num]) (* x0 x1)
 198
 199           (or (= [t0 t1] '[up down])
 200               (= [t0 t1] '[down up]))
 201           (apply + (map * (rest x0) (rest x1)))
 202
 203           (= t1 'down) (cons 'down (map #(n*2 x0 %) (rest x1)))
 204           (= t1 'up)   (cons 'up (map #(n*2 x0 %) (rest x1)))
 205           )))
 206
 207 (defn onehot [idxs tuple]
 208   (letfn [(rec [tuple pos]
 209             (if (coll? tuple)
 210               (cons (first tuple)
 211                     (map #(rec %1 (cons %2 pos))
 212                          (rest tuple)
 213                          (iterate inc 0)))
 214               (if (= pos idxs) 1 0)
 215               ))]
 216     (rec tuple [])))
 217
 218 (defn diff-tuple [tuple & [sel]]
 219   (letfn [(rec [pos sub-tuple]
 220             (if (coll? sub-tuple)
 221               (cons ({'up 'down, 'down 'up} (first sub-tuple))
 222                     (map #(rec (cons %2 pos) %1)
 223                          (rest sub-tuple) (iterate inc 0)
 224                          ))
 225               (reduce (fn [acc x] (nth acc (inc x)))
 226                       (onehot pos tuple)
 227                       (reverse sel))))]
 228     (rec [] tuple)))
 229
 230 ; symbolic differentiation
 231 (defn sdiff [expr]
 232   (cond (= expr 'arg) '[D arg]
 233         (not (coll? expr)) 0
 234         (= (first expr) 'arg) (cons 'D expr)
 235
 236         (= (first expr) 'expt)
 237         ['* (if (= (nth expr 2) 2)
 238               (nth expr 1)
 239               ['expt (nth expr 1) (dec (nth expr 2))])
 240             (sdiff (nth expr 1))
 241             ]))
 242