haskell-numerics / hmatrix

Linear algebra and numerical computation
381 stars 104 forks source link

Odd behaviour with 'minIndex' and 'maxIndex' #213

Closed fedeinthemix closed 7 years ago

fedeinthemix commented 7 years ago

I defined the following function to provide a default value when referring to an out of range index of a vector:

atIndexOrDefault :: (Container Vector e) => e -> Vector e -> IndexOf Vector -> e
atIndexOrDefault d v i
  | i < 0 || i >= (size v) = d
  | otherwise = v `atIndex` i

Then I tried the following variant:

atIndexOrDefault' :: (Container Vector e) => e -> Vector e -> IndexOf Vector -> e
atIndexOrDefault' d v i
  | i < minIndex v || i > maxIndex v = d
  | otherwise = v `atIndex` i

Both work with trivial examples in GHCi. However, in a bigger application, while the former gives the expected result, the latter give a different and unexpected result. Here is a snippet of the code using the above function:

import Numeric.LinearAlgebra
import qualified Data.Vector as V

...

adaptiveFirAdapt
  :: (Field a)
  => (Int -> Vector a -> a) 
  -> ((Vector a -> a) -> Vector a -> Vector a -> a -> AdaptiveFIRVT a)
  -> Vector a -> Vector a -> Vector a
  -> AdaptiveFIRVT a
adaptiveFirAdapt muOfkx step input desired w0 = do
  c <- ask
  let kEnd = adaptiveFirEndTime c
      o = adaptiveFirOrder c
      firIn = firInput input o
      ks = V.fromList [0..kEnd-1]
      atIndex' v i = atIndexOrDefault' 0 v i
  V.foldM' (\w k -> step (muOfkx k) w (firIn k) (desired `atIndex'` k)) w0 ks

Here the unexpected result with the second variant wrong.pdf and the correct, expected one with the first one correct.pdf

albertoruiz commented 7 years ago

minIndex gives the index of the minimum element, not the minimum index, which is always 0 (if the vector is not empty). In the same way, maxIndex is the maximum element, not the maximum index.

fedeinthemix commented 7 years ago

Oh, that explains the difference. Thanks!