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