[Haskell-cafe] HMatrix Vector/Matrix interpolation ala Matlab interp/interp2 ??

[Daniel Fischer]

On Tuesday 25 January 2011 23:16:49, gutti wrote:
> Hi Henning,
>
> thanks for the code review -- reason that I don't use the type
> declaration a lot -- It causes trouble , because I don't yet fully
> understand it.
>
> When I declare what I think is right is fails - see Message at bottom --
> so what's wrong ?

See inline.

>
> By the way I just used lists so far - no arrays -- Why is !! inefficient
> ? - what is better - Cheers Phil

(!!) is inefficient because it has to traverse the list from the start to 
the desired index. That's fine if you do it once, but if you do it a lot, 
you're probably doing something wrong.

>
> 1 import Data.List
> 2
> 3 -- Input Data
> 4 xi :: [Double]
> 5 xi = [0 .. 10]
> 6 yi :: [Double]
> 7 yi = [2, 3, 5, 6, 7, 8, 9, 10 , 9, 8, 7]
> 8 x = 11 :: Double
> 9
> 10 -- Functions
> 11 limIndex :: [a] -> Int -> Int
> 12 limIndex xi idx
> 13    | idx < 0 = 0
> 14     | idx > (length xi)-2 = (length xi)-2
> 15    | otherwise = idx
> 16
> 17 getIndex :: [a] -> Double -> Int
> 18 getIndex xi x = limIndex xi (maybe (length xi) id
>                                          (findIndex (>x) xi)-1)

The type of (>) is

(>) :: Ord a => a -> a -> Bool

, so both arguments to (>) must have the same type and that type must be an 
instance of the Ord class.
Since the type of x is stated to be Double, (> x) :: Double -> Bool, so the 
list needs to have the type [Double].

You can ask ghci what the (most general) type of your function is, here 
it's

getIndex :: Ord a => [a] -> a -> Int

> 19
> 20 getPnts :: [a] -> [a] -> Int -> [a]
> 21 getPnts xi yi idx = [xi !! idx, xi !! (idx+1), yi !! idx,
>                                                           yi !! (idx+1)]
> 22
> 23 interp :: [a] -> [a] -> Double -> Double
> 24 interp xi yi x =
> 25 	let pts = getPnts xi yi (getIndex xi x)
> 26 	in (pts!!3-pts!!2)/(pts!!1-pts!!0)*(x-pts!!0)+pts!!2

(-) :: Num a => a -> a -> a

The two arguments of (-) must have the same type and the result has the 
same type too. x is stated to be a Double by the signature, so we must have 
pts :: [Double] and hence xi :: [Double], yi :: [Double].

For the most general type, the use of getIndex implies an Ord constraint, 
(-) and (*) give a Num constraint and

(/) :: Fractional a => a -> a -> a

, so the use of (/) adds a Fractional constraint. Fractional implies Num, 
hence

interp :: (Ord a, Fractional a) => [a] -> [a] -> a -> a

> 27
> 28 -- Calc
> 29 y = interp xi yi x
> 30
> 31 main = do
> 32	-- Output Data
> 33   	print (y)
>
>  === Compiler Error Message ===
>
> interp_v4.hs:18:66:
>     Couldn't match expected type `Double' against inferred type `a'
>       `a' is a rigid type variable bound by
>           the type signature for `getIndex' at interp_v4.hs:17:13
>       Expected type: [Double]
>       Inferred type: [a]
>     In the second argument of `findIndex', namely `xi'
>     In the third argument of `maybe', namely `(findIndex (> x) xi)'
>
> interp_v4.hs:26:39:
>     Couldn't match expected type `Double' against inferred type `a'
>       `a' is a rigid type variable bound by
>           the type signature for `interp' at interp_v4.hs:23:11
>     In the second argument of `(-)', namely `pts !! 0'
>     In the second argument of `(*)', namely `(x - pts !! 0)'
>     In the first argument of `(+)', namely
>         `(pts !! 3 - pts !! 2) / (pts !! 1 - pts !! 0) * (x - pts !! 0)'