haskell-numerics / hmatrix

Linear algebra and numerical computation
381 stars 104 forks source link

change unwrap to extract in overMatL' and overMatM' #364

Open jjdosa opened 1 year ago

jjdosa commented 1 year ago

@lsg6140 and I've found this bug in the overMatL' and overMatM' functions.

A reduced example case we found is that

{-# LANGUAGE DataKinds #-}

module Main where

import qualified Numeric.LinearAlgebra.Static as H
rmat = H.matrix [1,2,3,4] :: H.L 2 2
cmat = H.matrix [1,2,3,4] :: H.M 2 2

invDiagR = H.inv $ H.diag (H.vector [1,2] :: H.R 2)
invDiagC = H.inv $ H.diag (H.vector [1,2] :: H.C 2)

After loading this module in ghci, executing invDiagR and invDiagC produces these errors:

*Main> invDiagR
*** Exception: inv of nonsquare (1x3) matrix
CallStack (from HasCallStack):
  error, called at src/Internal/Algorithms.hs:706:21 in hmatrix-0.20.2-inplace:Internal.Algorithms
*Main> invDiagC
*** Exception: inv of nonsquare (1x3) matrix
CallStack (from HasCallStack):
  error, called at src/Internal/Algorithms.hs:706:21 in hmatrix-0.20.2-inplace:Internal.Algorithms

By replacing unwrap with extract in the definition of overMatL and overMatM functions, we've checked that we have the expected results:

*Main> invDiagR 
(matrix
 [ 1.0, 0.0
 , 0.0, 0.5 ] :: L 2 2)
*Main> invDiagC
(matrix
 [ 1.0 :+ 0.0, 0.0 :+ 0.0
 , 0.0 :+ 0.0, 0.5 :+ 0.0 ] :: M 2 2)
idontgetoutmuch commented 1 year ago

I am struggling with the linker on my m1 Mac atm - I hope to fix this tomorrow and then I can test your PR - thanks for contributing :-)

jjdosa commented 1 year ago

Thank you @idontgetoutmuch for the feedback!

You can ignore the rmat and camt matrices as they are unused.

rmat = H.matrix [1,2,3,4] :: H.L 2 2
cmat = H.matrix [1,2,3,4] :: H.M 2 2

They are there because I forgot to remove them after reducing these examples

invDiagR = H.inv $ H.diag (H.takeDiag rmat)
invDiagC = H.inv $ H.diag (H.takeDiag cmat)

into the ones in the PR message.