Numeric types are unexpectedly coerced
Hello,
I'm trying to do some very basic matrix operation using the default implementation (persistent vector) of core.matrix 0.62.0 and I'd love to keep my rational elements as clojure.lang.Ratio as long as I can during my matrix transformations.
I noticed that when I sum different diagonal matrices using rational elements they get coerced to clojure.lang.Double.
I tracked down the issue to clojure.core.matrix.impl.defaults:
(extend-protocol mp/PSpecialisedConstructors
#?(:clj Object :cljs object)
(identity-matrix [m dims]
(mp/diagonal-matrix m (repeat dims 1.0)))
(diagonal-matrix [m diagonal-values]
(let [dims (count diagonal-values)
diagonal-values (mp/convert-to-nested-vectors diagonal-values)
zs (vec (repeat dims 0.0))
dm (vec (for [i (range dims)]
(assoc zs i (nth diagonal-values i))))]
(mp/coerce-param m dm))))
Replacing zs (vec (repeat dims 0.0)) with zs (vec (repeat dims 0)) in the let bindings makes it work as expected (for me).
I didn't check if this breaks some test, but if you whish I can submit a pull request.
Here you can find a test case:
(require '[clojure.core.matrix :as matrix])
(require '[clojure.math.numeric-tower :as math])
(defn wilkinson [n]
(let [side (repeat n 1)
center (matrix/emap math/abs (range (- (/ (dec n) 2))
(+ (/ n 2))))]
(matrix/add
(matrix/shift (matrix/diagonal-matrix side) [-1])
(matrix/diagonal-matrix center)
(matrix/shift (matrix/diagonal-matrix side) [1]))))
;; without the fix
(wilkinson 6)
[[2.5 1.0 0.0 0.0 0.0 0.0]
[1.0 1.5 1.0 0.0 0.0 0.0]
[0.0 1.0 0.5 1.0 0.0 0.0]
[0.0 0.0 1.0 0.5 1.0 0.0]
[0.0 0.0 0.0 1.0 1.5 1.0]
[0.0 0.0 0.0 0.0 1.0 2.5]]
;; with the fix
(wilkinson 6)
[[5/2 1 0 0 0 0]
[1 3/2 1 0 0 0]
[0 1 1/2 1 0 0]
[0 0 1 1/2 1 0]
[0 0 0 1 3/2 1]
[0 0 0 0 1 5/2]]
As you can see I also had to map over my range with math/abs in the let binding of my example because matrix/abs coerced too my items.
Thank you.
Hi, yes this does look like unexpected behaviour. Arrays should preserve the type of elements as far as possible, although the extent to which they do this is somewhat up to the implementation.
The correct fix would be to make sure arrays are initialised with correct zeros and ones for the implementation. See c.c.m.generic/zero and c.c.m.generic/one c.c.m.generic/zero https://cljdoc.org/d/net.mikera/core.matrix/0.62.0/api/clojure.core.matrix.generic#zero
PR would be welcome!