(ns clatrix.core
  (:refer-clojure :exclude [get set map-indexed map rand vector? + - * pp])
  (:use [slingshot.slingshot :only [throw+]])
  (:import [org.jblas DoubleMatrix ComplexDoubleMatrix ComplexDouble
            Decompose Decompose$LUDecomposition Eigen Solve Geometry
            Singular MatrixFunctions]))

(declare get permute size matrix matrix? row? nrows ncols vstack)

(deftype Matrix [^DoubleMatrix me ^Boolean vector? ^clojure.lang.IPersistentMap metadata]
  Object
  (toString [^Matrix mat]
    (str (list `matrix
               (vec (clojure.core/map vec (vec (.toArray2 ^DoubleMatrix (.me mat))))))))
  clojure.lang.IObj
  (withMeta [this metadata]
    (Matrix. me vector? metadata))
  (meta [this]
    metadata)
  clojure.lang.ISeq
  (equiv [this that]
    (cond
      (matrix? that) (.equals (.me this) (.me that))
      (coll? that) (and (= (count this) (count that))
                       (every? true? (clojure.core/map #(== %1 %2) this that))) 
      (number? that) (and (= [1 1] (size this)) (= that (get this 0 0)))
      :else (.equals (.me this) that)))
  (first [this]
    (let [[r c] (size this)]
      (cond
        (or (zero? r) (zero? c)) nil
        (and vector? (or (= r 1) (= c 1))) (get this 0 0)
        :else (let [out (get this 0 (range c))]
                (if (number? out)
                  (matrix (vector out))
                  out)))))
  (more [this]
    (if-let [nxt (next this)]
      nxt
      (matrix [])))
  (cons [this x]
    (cond
      (matrix? x)  (vstack this x)
      (and (coll? x) (number? (first x))) (vstack this (matrix (vector x)))
      :else (vstack this (matrix x))))
  (seq [this]
    (let [[r c] (size this)]
      (when-not (or (zero? r) (zero? c))
        (Matrix. me vector? nil))))
  (next [this]
    (let [[r c] (size this)]
      (cond
        (and vector? (= 1 c) (> r 1)) (matrix (.me (get this (range 1 r) 0)) true nil)
        (and vector? (= 1 r) (> c 1)) (matrix (.me (get this 0 (range 1 c))) true nil)
        (and (not vector?) (> r 1)) (matrix (.me (get this (range 1 r) (range c))) false nil)
        :else nil)))
  (empty [this]
    (matrix []))
  clojure.lang.Counted
  (count [this]
    (if vector?
      (clojure.core/* (nrows this) (ncols this))
      (nrows this)))
  clojure.lang.Sequential)

(defn me [^Matrix mat]
  (.me mat))

(defmacro dotom [name m & args]
  `(~name ^DoubleMatrix (me ~m) ~@args))

(defmacro promote-cfun* [defname name fname]
  (let [n (gensym)
        m (gensym)]
    `(~defname ~name
               ([^long ~n] (matrix (~fname ~n)))
               ([^long ~n ~m] (matrix (~fname ~n ~m))))))

(defmacro promote-mfun* [defname name fname]
  (let [m (gensym)]
    `(~defname ~name [^Matrix ~m] (dotom ~fname ~m))))

(defn matrix?
  [m] (isa? (class m) Matrix))

(promote-mfun* defn ncols .columns)
(promote-mfun* defn nrows .rows)
(defn size    [^Matrix m] [(nrows m) (ncols m)])
(defn vector? [^Matrix m] (.vector? m))
(defn row?    [^Matrix m] (== 1 (first (size m))))
(defn column? [^Matrix m] (== 1 (second (size m))))
(defn square? [^Matrix m] (reduce == (size m)))

(defn- int-arraytise
  [x] (if (coll? x) (int-array x) x))

(defn get [^Matrix m r c]
  (let [out (dotom .get m (int-arraytise r) (int-arraytise c))]
    (if (number? out)
      out
      (matrix out))))
(defn set [^Matrix m ^long r ^long c ^double e]
  (dotom .put m r c e))

(defn trace
  [^Matrix mat]
  (if (square? mat)
    (let [[n _] (size mat)]
      (reduce #(clojure.core/+ (get mat %2 %2) %1) 0 (range n)))
    (throw+ {:error "Cannot take trace of non-square matrix."})))

(defn dense
    [^Matrix m]
  (vec (clojure.core/map vec (vec (dotom .toArray2 m)))))

(defn as-vec
  [^Matrix m]
  (if (vector? m)
    (vec (dotom .toArray m))
    (dense m)))

(defn column  ;; TODO remove from api
  [^doubles seq]
  (matrix seq))

(derive java.util.Collection ::collection)
(derive DoubleMatrix ::double-matrix)
(defmulti matrix
  (fn [m & args] (class m)))

(defmethod matrix ::double-matrix
  ([^DoubleMatrix x]
   (matrix x (.isVector x) nil)) 
  ([^DoubleMatrix x vector? meta]
   (Matrix. x vector? meta)))

(defmethod matrix ::collection
  [seq-of-seqs]
  (if (number? (first seq-of-seqs))
    (matrix (DoubleMatrix. (into-array Double/TYPE (clojure.core/map double seq-of-seqs))))
    (let [lengths (clojure.core/map count seq-of-seqs)
          flen    (first lengths)]
      (cond
        (or (= (count lengths) 0) (some zero? lengths)) (matrix (DoubleMatrix. 0 0))
        (every? (partial = flen) lengths)
        (matrix
          (DoubleMatrix.
            ^"[[D" (into-array (clojure.core/map #(into-array Double/TYPE (clojure.core/map double %)) seq-of-seqs))))
        :else (throw+ {:error "Cannot create a ragged matrix."})))))

(defn diag
  [seq-or-matrix]
  (if (matrix? seq-or-matrix)
    (let [mat ^Matrix seq-or-matrix]
      ;; We'll extract largest diagonals from non-square matrices
      ;; since this isn't really a matrix algebraic property
      (let [n (apply min (size mat))]
        (clojure.core/map #(get mat % %) (range n))))
    (let [di (clojure.core/map double (seq seq-or-matrix))]
      (matrix (DoubleMatrix/diag (DoubleMatrix. ^doubles (into-array Double/TYPE di)))))))

(promote-cfun* defn  ones DoubleMatrix/ones) ;; TODO wrap inside matrix

(promote-cfun* defn zeros DoubleMatrix/zeros)
(defn constant
  ([^long n ^double c]
   (matrix (doto (DoubleMatrix/ones n)
             (.muli c))))
  ([^long n ^long m ^double c]
   (matrix (doto (DoubleMatrix/ones n m)
             (.muli c)))))

(defn id
  [^long n] (matrix (DoubleMatrix/eye n)))

(defn reshape
  [^Matrix A p q]
  (let [[n m] (size A)]
    (if (= (clojure.core/* n m) (clojure.core/* p q))
      (dotom A p q)
      (throw+ {:exception "Cannot change the number of elements during a reshape."
               :previous (clojure.core/* n m)
               :new (clojure.core/* p q)}))))

(defn from-sparse
  [^long n ^long m specs]
  (let [m (zeros n m)]
    (doseq [[i j v] specs]
      (set m i j v))
    m))

(defn from-indices
  [^long n ^long m fun]
  (let [ary (make-array Double/TYPE n m)]
    (doseq [i (range n)
            j (range n)]
      (aset ary i j (fun i j)))
    (matrix (DoubleMatrix. ^"[[D" ary))))

(promote-cfun* defn  rand   DoubleMatrix/rand)
(promote-cfun* defn- randn* DoubleMatrix/randn)

(defn rnorm
  ([^double mu ^double sigma ^long n ^long m]
   (matrix (doto (dotom .muli (randn* n m) sigma)
             (.addi mu))))
  ([^double mu ^double sigma ^long n]
   (matrix (doto (dotom .muli (randn* n) sigma)
             (.addi mu))))
  ([^long n ^long m] (randn* n m))
  ([^long n] (randn* n)))

(defn t
  [^Matrix mat] (matrix (dotom .transpose mat)))

(defn hstack
  [& vec-seq]
  (let [row-counts (clojure.core/map nrows vec-seq)
        rows (first row-counts)]
    (if (every? (partial == rows) row-counts)
      (matrix (reduce #(DoubleMatrix/concatHorizontally
                         ^DoubleMatrix %1
                         ^DoubleMatrix (me %2))
                      (me (first vec-seq))
                      (rest vec-seq))))))

(defn vstack
  [& vec-seq]
  (let [col-counts (clojure.core/map ncols vec-seq)
        cols (first col-counts)]
    (if (every? (partial == cols) col-counts)
      (matrix (reduce #(DoubleMatrix/concatVertically
                         ^DoubleMatrix %1
                         ^DoubleMatrix (me %2))
                      (me (first vec-seq))
                      (rest vec-seq))))))

(defn rows
  ([^Matrix mat]   ;; TODO use map and get on datastruct
   (let [[n m] (size mat)]
     (rows mat (range n))))
  ([^Matrix m ^longs idxs]
   (clojure.core/map #(matrix (dotom .getRow m %)) idxs)))

(defn cols
  ([^Matrix mat]
   (cols mat (range (ncols mat))))
  ([^Matrix m ^longs idxs]
   (clojure.core/map #(matrix (dotom .getColumn m %)) idxs)))

(defn- permute-rows [^Matrix mat rowspec]
  (if rowspec (apply vstack (rows mat rowspec)) mat))

(defn- permute-cols [^Matrix mat colspec]
  (if colspec (apply hstack (cols mat colspec)) mat))

(defn permute
    [^Matrix
                                                            mat & {:keys [r c rowspec colspec]}]
  (let [[n m] (size mat)
        r (or r rowspec)
        c (or c colspec)]
    (cond (and r (some #(> % (dec n)) r))
          (throw+ {:error "Row index out of bounds" :num-rows n :rowspec r})
          (and c (some #(> % (dec m)) c))
          (throw+ {:error "Column index out of bounds" :num-columns n :colspec c})
          :else (permute-cols (permute-rows mat r) c))))

(defn- iswild
  [sym]
  ;; This is a sort of silly way to do it, but I can't get the regex
  ;; to work for both '_ and #'user/_
  (let [name1 (second (re-find #"/(.+)" (str sym)))
        name2 (str sym)]
    (or (= name1 "_")
        (= name1 ".")
        (= name1 "*")
        (= name2 "_")
        (= name2 ".")
        (= name2 "*"))))

(defn- make-constr
  [e]
  (if (matrix? e)
    {:matrix e
     :rows (first (size e))
     :cols (second (size e))}
    {:constant e}))

(defn- update-hash-with-constr
  [hsh constr i j key]
  (let [{n key} constr
        old (hsh [i j])]
    (if (key old)
      (if (not= (key old) n) ;the constraint doesn't match the hash, uh oh
        (throw+ {:error "Block matrix diagram sizes are inconsistent."
                 :type :constraint-error
                 :location [i j]})
        hsh)
      ;; if there isn't an old key then we can fix that constraint now
      (assoc hsh [i j]
             (assoc old key n)))))

(defn block-fn
  [matrices]
  ;; We must do size-constraint propagation along the rows and columns
  ;; of the block-diagram in order to (a) ensure that the input isn't
  ;; in error and (b) find the proper sizes for the constant matrices.
  (let [n       (count matrices)
        lengths (clojure.core/map count matrices)
        m       (first lengths)]
    (if (not (every? (partial == m) lengths))
      (throw+ {:error "Block matrices cannot be ragged."})
      ;; Build the constraints map
      (let [indices (for [i (range n) j (range m)] [i j])
            ;; The initial hash map contains what we know before
            ;; constraint propagation.
            init-map (reduce (fn [hsh [i j]]
                               (assoc hsh [i j]
                                      (make-constr (nth (nth matrices i) j))))
                             (hash-map) indices)
            ;; Walk over the index set and propagate all the constraints
            ;; over each row and column
            constraint-map
            (reduce
              (fn [hash [i j]]
                (let [constr (init-map [i j])]
                  (if (or (:rows constr) (:cols constr))
                    ;; Look up and to the left for constraint violations,
                    ;; locking in constraints if they don't already exist
                    (reduce #(update-hash-with-constr %1 constr i %2 :rows)
                            (reduce #(update-hash-with-constr %1 constr %2 j :cols)
                                    hash (range n))
                            (range m))
                    hash)))
              init-map
              indices)]
        ;; Use the constraint map to build the final matrix
        (apply vstack
               (for [i (range n)]
                 (apply hstack
                        (for [j (range m)]
                          (let [constr (constraint-map [i j])]
                            (if (:matrix constr)
                              (:matrix constr)
                              ;; Constants are assumed to be 1x1
                              ;; unless otherwise constrained
                              (constant (:rows constr 1)
                                        (:cols constr 1)
                                        (:constant constr))))))))))))

(defmacro block
  [blockspec]
  `(block-fn
     ~(vec (clojure.core/map
             #(vec (clojure.core/map
                     (fn [e] (if (iswild e) 0 e)) %))
             blockspec))))

(defn- slicer
  ([^Matrix matrix rowspec colspec]
   (cond (and (iswild rowspec) (iswild colspec)) matrix
         (iswild rowspec) `(matrix (dotom .getColumn ~matrix ~colspec))
         (iswild colspec) `(matrix (dotom .getRow    ~matrix ~rowspec))
         :else            `(get ~matrix ~rowspec ~colspec)))
  ([^DoubleMatrix matrix rowspec colspec values]
   (let [m (gensym)
         form (cond (and (iswild rowspec) (iswild colspec)) `(.copy ~m ~values)
                    (iswild rowspec) `(dotom .putColumn ~m ~colspec ~values)
                    (iswild colspec) `(dotom .putRow    ~m ~rowspec ~values)
                    :else            `(set        ~m ~rowspec ~colspec ~values))]
     `(let [~m ~matrix]
        (do ~form ~m)))))

(defmacro slice
  [^Matrix matrix rowspec colspec & values?]
  (apply slicer matrix rowspec colspec values?))

(defmacro slices
  [^Matrix matrix rowspec colspec & values?]
  `(as-vec ~(apply slicer matrix rowspec colspec values?)))

(defn symmetric
  [^Matrix m] (with-meta m {:symmetric true}))

(defn positive
  [^Matrix m] (with-meta m {:symmetric true :positive true}))

(defn arbitrary
  [^Matrix m] (with-meta m {:symmetric false :positive false}))

(defn symmetric? [^Matrix m] (:symmetric (meta m)))
(defn positive?  [^Matrix m] (:positive (meta m)))
(defn arbitrary? [^Matrix m] (not (or (symmetric? m) (positive? m))))

(defn maybe-symmetric
  [^Matrix m]
  (if (or (symmetric? m) (= (t m) m))
    (symmetric m)
    (with-meta m {:symmetric false})))

(defn maybe-positive
  [^Matrix m]
  (let [m (maybe-symmetric m)]
    (if (symmetric? m)
      ;; We'll have faster access to the eigenvalues later...
      (let [vals (dotom Eigen/eigenvalues m)
            rvals (seq (.toArray (.real vals)))
            ivals (seq (.toArray (.real vals)))
            mags (clojure.core/map #(Math/sqrt
                                      (clojure.core/+ (Math/pow %1 2)
                                                      (Math/pow %2 2))) rvals ivals)]
        (if (every? #(> % 0) mags)
          (positive m)
          m))
      m)))

(defn map-indexed
  [fun ^Matrix mat]
  (let [[n m] (size mat)]
    (from-sparse n m
                 (for [i (range n)
                       j (range m)]
                   [i j (fun i j (get mat i j))]))))

(defn map
  [fun ^Matrix mat]
  (let [[n m] (size mat)]
    (from-sparse n m
                 (for [i (range n)
                       j (range m)]
                   [i j (fun (get mat i j))]))))

(defn norm
  [^Matrix mat]
  (dotom .norm2 mat))

(defn normalize
  [^Matrix mat & [flags]]
  (if (column? mat)
    (matrix (dotom Geometry/normalize mat))
    (matrix (dotom Geometry/normalizeColumns mat))))

(defn +
  ([a b] (cond (and (matrix? a) (matrix? b))
               (if (= (size a) (size b))
                 (matrix (dotom .add a ^DoubleMatrix (me b)))
                 (throw+ {:exception "Matrices of different sizes cannot be summed."
                          :asize (size a)
                          :bsize (size b)}))
               (matrix? a) (matrix
                             (dotom .add a (double b)))
               (matrix? b) (matrix
                             (dotom .add b (double a)))
               :else       (clojure.core/+ a b)))
  ([a b & as] (reduce + a (cons b as))))

(defn *
  ([a b] (cond (and (matrix? a) (matrix? b))
               (if (= (second (size a)) (first (size b)))
                 (matrix (dotom .mmul a ^DoubleMatrix (me b)))
                 (throw+ {:exception "Matrix products must have compatible sizes."
                          :a-cols (ncols a)
                          :b-rows (nrows b)}))
               (matrix? a) (matrix
                             (dotom .mmul a (double b)))
               (matrix? b) (matrix
                             (dotom .mmul b (double a)))
               :else       (clojure.core/* a b)))
  ([a b & as] (reduce * a (cons b as))))

(defn mult
  ([a b] (cond (and (matrix? a) (matrix? b))
               (matrix (dotom .mul a ^DoubleMatrix (me b)))
               (matrix? a) (matrix
                             (dotom .mul a (double b)))
               (matrix? b) (matrix
                             (dotom .mul b (double a)))
               :else       (clojure.core/* a b))))

(defn div
  ([a b] (cond (and (matrix? a) (matrix? b))
               (matrix (dotom .div a ^DoubleMatrix (me b)))
               (matrix? a) (matrix
                             (dotom .div a (double b)))
               (matrix? b) (matrix
                             (dotom .rdiv b (double a)))
               :else       (clojure.core// a b))))

(defn -
  ([a] (* -1 a))
  ([a b] (cond (and (matrix? a) (matrix? b))
               (if (= (size a) (size b))
                 (matrix (dotom .sub a ^DoubleMatrix (me b)))
                 (throw+ {:exception "Matrices of different sizes cannot be differenced."
                          :asize (size a)
                          :bsize (size b)}))
               (matrix? a) (matrix
                             (dotom .sub a (double b)))
               (matrix? b) (matrix
                             (dotom .rsub b (double a)))
               :else       (clojure.core/- a b)))
  ([a b & as] (reduce - a (cons b as))))

(defn dot
  [^Matrix m1 ^Matrix m2]
  (dotom .dot m1 ^DoubleMatrix (me m2)))

(defn rreflection
  [n]
  (let [v (Geometry/normalize (DoubleMatrix/randn n))]
    (matrix
      (.sub (DoubleMatrix/eye n) (.mmul (.mmul v (.transpose v)) (double 2))))))

(defn rspectral
  [n-or-spectrum]
  (let [[n spectrum]
        (if (sequential? n-or-spectrum)
          [(count n-or-spectrum) n-or-spectrum]
          [n-or-spectrum (repeatedly n-or-spectrum clojure.core/rand)])
        V ^DoubleMatrix
        (nth (iterate (fn [^DoubleMatrix prod]
                        (.mmul prod ^DoubleMatrix (me (rreflection n))))
                      (DoubleMatrix/eye n))
             (clojure.core/* 2 n))
        L (DoubleMatrix/diag
            (DoubleMatrix.
              ^doubles (into-array Double/TYPE spectrum)))
        A (matrix (-> V (.mmul L) (.mmul (.transpose V))))]
    (if (every? pos? spectrum)
      (positive A)
      A)))

(defn solve
  [^Matrix A ^Matrix B]
  (matrix
    (cond
      (positive? A)  (Solve/solvePositive ^DoubleMatrix (me A)
                                          ^DoubleMatrix (me B))
      (symmetric? A) (Solve/solveSymmetric ^DoubleMatrix (me A)
                                           ^DoubleMatrix (me B))
      :else          (Solve/solve ^DoubleMatrix (me A)
                                  ^DoubleMatrix (me B)))))

(defn i
  [^Matrix mat]
  (if (not (square? mat))
    (throw+ {:exception "Cannot invert a non-square matrix."})
    (let [[n _] (size mat)]
      (solve mat (id n)))))

(defn eigen
  ([^Matrix A]
   (cond
     (symmetric? A) (let [[vecs vals] (clojure.core/map #(matrix %)
                                                        (seq (dotom Eigen/symmetricEigenvectors A)))]
                      {:vectors vecs :values (diag vals)})
     :else          (let [[^ComplexDoubleMatrix vecs ^ComplexDoubleMatrix vals]
                          (seq (dotom Eigen/eigenvectors A))
                          rvecs (matrix (.real vecs))
                          ivecs (matrix (.imag vecs))
                          rvals (diag (matrix (.real vals)))
                          ivals (diag (matrix (.imag vals)))
                          out {:vectors rvecs
                               :values  rvals}]
                      (if (some (partial not= 0.0) ivals)
                        (merge out
                               {:ivectors ivecs
                                :ivalues  ivals})
                        out))))
  ([^Matrix A ^Matrix B]
   (let [A (maybe-symmetric A)
         B (maybe-symmetric B)]
     (if (and (symmetric? A) (symmetric? B))
       (let [[vecs vals]
             (clojure.core/map #(matrix %)
                               (seq (Eigen/symmetricGeneralizedEigenvectors ^DoubleMatrix (me A)
                                                                            ^DoubleMatrix (me B))))]
         {:vectors vecs :values (as-vec vals)})
       (throw+ {:error "Cannot do generalized eigensystem for non-symmetric matrices."})))))

(defn svd
  [^Matrix A & {:keys [type] :or {type :sparse}}]
  (let [[U L V] (if (= type :full)
                  (dotom Singular/fullSVD A)
                  (dotom Singular/sparseSVD A)) 
        left (matrix U)
        right (matrix V)
        values (seq (.toArray L))]
    (if (= type :values)
      {:values (seq (.toArray (dotom Singular/SVDValues A)))}
      {:left left
       :right right
       :values values
       :rank (count values)})))

(defn rank
  [^Matrix A] (:rank (svd A)))

(defn- complex-power
  [^ComplexDouble z ^double e]
  (let [r (.real z)
        i (.imag z)
        m (Math/sqrt (clojure.core/+ (Math/pow r 2) (Math/pow i 2)))
        a (Math/atan2 i r)
        m2 (Math/pow m e)
        a2 (clojure.core/* a e)]
    (ComplexDouble. (clojure.core/* m2 (Math/cos a2))
                    (clojure.core/* m2 (Math/sin a2)))))

(defn pow
  [^Matrix A e]
  (if (not (square? A))
    (throw+ {:exception "Cannot take power of non-square matrix."
             :size (size A)})
    (let [[n _] (size A)
          [^ComplexDoubleMatrix V
           ^ComplexDoubleMatrix L]
          (seq (dotom Eigen/eigenvectors A))
          new-spect (clojure.core/map #(complex-power % e)
                                      (seq (.toArray (.diag L))))]
      ;; Update L with the new spectrum
      (doseq [i (range n)]
        (.put L i i (nth new-spect i)))
      ;; Compute the product
      (-> V
        (.mmul L)
        (.mmul (.transpose V))
        .real
        (matrix)))))

(defn cholesky
  [^Matrix mat]
  (let [mat (maybe-positive mat)]
    (if (positive? mat)
      (matrix (dotom Decompose/cholesky mat))
      (throw+ {:exception "Cholesky decompositions require positivity."}))))

(defn lu
  [^Matrix mat]
  (if (not (square? mat))
    (throw+ {:exception "Cannot compute LU decomposition of a non-square matrix."
             :size (size mat)})
    (let [lu (dotom Decompose/lu mat)]
      {:p (matrix ^DoubleMatrix (.p lu))
       :l (matrix ^DoubleMatrix (.l lu))
       :u (matrix ^DoubleMatrix (.u lu))})))

(defmethod print-method Matrix [mat ^java.io.Writer w]
  (let [[nbits prec] [3 2]
        [n m] (size mat)
        small-rows (< n (* nbits 2))
        small-cols (< m (clojure.core/* nbits 2))
        rowset (if (< n (clojure.core/* nbits 2))
                 (range n)
                 (concat (range nbits) (range (clojure.core/- n nbits) n)))
        colset (if (< m (clojure.core/* nbits 2))
                 (range m)
                 (concat (range nbits) (range (clojure.core/- m nbits) m)))
        submat (apply hstack (cols (apply vstack (rows mat rowset)) colset))
        [n m] (size submat)
        fmt (str "% ." prec "e ")
        header (apply format " A %dx%d matrix\n" (size mat))]
    ;; Print the header
    (.write w header)
    (.write w " ")
    (doseq [i (range (dec (count header)))] (.write w "-"))
    (.write w "\n")
    ;; Print the matrix
    (if small-rows
      (doseq [i (range n)]
        (if small-cols
          (doseq [j (range m)]
            (.write w (format fmt (slice submat i j))))
          (do (doseq [j (range nbits)]
                (print (format fmt (slice submat i j))))
            (.write w " . ")
            (doseq [j (range nbits (clojure.core/* 2 nbits))]
              (.write w (format fmt (slice submat i j))))))
        (.write w "\n"))
      (do (doseq [i (range nbits)]
            (if small-cols
              (doseq [j (range m)]
                (.write w (format fmt (slice submat i j))))
              (do (doseq [j (range nbits)]
                    (.write w (format fmt (slice submat i j))))
                (.write w " . ")
                (doseq [j (range nbits (clojure.core/* 2 nbits))]
                  (.write w (format fmt (slice submat i j))))))
            (.write w "\n"))
        (.write w " ... \n")
        (doseq [i (range nbits (clojure.core/* 2 nbits))]
          (if small-cols
            (doseq [j (range m)]
              (.write w (format fmt (slice submat i j))))
            (do (doseq [j (range nbits)]
                  (.write w (format fmt (slice submat i j))))
              (.write w " . ")
              (doseq [j (range nbits (clojure.core/* 2 nbits))]
                (.write w (format fmt (slice submat i j))))))
          (.write w "\n"))))))

(defmacro promote-mffun* [defname name fname]
  (let [m (gensym)]
    `(~defname ~name ^Matrix [^Matrix ~m] (matrix (dotom ~fname ~m)))))

(promote-mffun* defn exp    MatrixFunctions/exp)
(promote-mffun* defn abs    MatrixFunctions/abs)
(promote-mffun* defn acos   MatrixFunctions/acos)
(promote-mffun* defn asin   MatrixFunctions/asin)
(promote-mffun* defn atan   MatrixFunctions/atan)
(promote-mffun* defn cbrt   MatrixFunctions/cbrt)
(promote-mffun* defn ceil   MatrixFunctions/ceil)
(promote-mffun* defn cos    MatrixFunctions/cos)
(promote-mffun* defn cosh   MatrixFunctions/cosh)
(promote-mffun* defn exp    MatrixFunctions/exp)
(promote-mffun* defn floor  MatrixFunctions/floor)
(promote-mffun* defn log    MatrixFunctions/log)
(promote-mffun* defn log10  MatrixFunctions/log10)
(promote-mffun* defn signum MatrixFunctions/signum)
(promote-mffun* defn sin    MatrixFunctions/sin)
(promote-mffun* defn sinh   MatrixFunctions/sinh)
(promote-mffun* defn sqrt   MatrixFunctions/sqrt)
(promote-mffun* defn tan    MatrixFunctions/tan)
(promote-mffun* defn tanh   MatrixFunctions/tanh)

(promote-mffun* defn exp!    MatrixFunctions/expi)
(promote-mffun* defn abs!    MatrixFunctions/absi)
(promote-mffun* defn acos!   MatrixFunctions/acosi)
(promote-mffun* defn asin!   MatrixFunctions/asini)
(promote-mffun* defn atan!   MatrixFunctions/atani)
(promote-mffun* defn cbrt!   MatrixFunctions/cbrti)
(promote-mffun* defn ceil!   MatrixFunctions/ceili)
(promote-mffun* defn cos!    MatrixFunctions/cosi)
(promote-mffun* defn cosh!   MatrixFunctions/coshi)
(promote-mffun* defn exp!    MatrixFunctions/expi)
(promote-mffun* defn floor!  MatrixFunctions/floori)
(promote-mffun* defn log!    MatrixFunctions/logi)
(promote-mffun* defn log10!  MatrixFunctions/log10i)
(promote-mffun* defn signum! MatrixFunctions/signumi)
(promote-mffun* defn sin!    MatrixFunctions/sini)
(promote-mffun* defn sinh!   MatrixFunctions/sinhi)
(promote-mffun* defn sqrt!   MatrixFunctions/sqrti)
(promote-mffun* defn tan!    MatrixFunctions/tani)
(promote-mffun* defn tanh!   MatrixFunctions/tanhi)