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

Henning Thielemann lemming at henning-thielemann.de
Tue Jan 25 22:27:10 CET 2011

```On Tue, 25 Jan 2011, gutti wrote:

> I created some code from scratch - probably "ugly" beginners style - so
> I'm keen to get tips how to make it more pretty and faster

> import Data.List
>
> -- Input Data
> xi :: [Double]
> xi = [0 .. 10]
> yi :: [Double]
> yi = [2, 3, 5, 6, 7, 8, 9, 10 , 9, 8, 7]
> x = 11 :: Double

> -- Functions
> limIndex xi idx
>    | idx < 0 = 0
>    | idx > (length xi)-2 = (length xi)-2
>    | otherwise = idx

limIndex :: [a] -> Int -> Int
limIndex xi idx = max 0 (min (length xi - 2) idx)

> getIndex xi x = limIndex xi (maybe (length xi) id (findIndex (>x) xi)-1)
>
> getPnts xi yi idx = [xi !! idx, xi !! (idx+1), yi !! idx, yi !! (idx+1)]

Since this list has always four elements, I suggest a quadruple:

getPnts xi yi idx = (xi !! idx, xi !! (idx+1), yi !! idx, yi !! (idx+1))

(!!) is not very efficient, but for now I imagine that's an access to a
HMatrix-Vector.

> interp xi yi x =
> 	let pts = getPnts xi yi (getIndex xi x)
> 	in (pts!!3-pts!!2)/(pts!!1-pts!!0)*(x-pts!!0)+pts!!2

let (x0,x1,y0,y1) = getPnts xi yi (getIndex xi x)
in  (y1-y0)/(x1-x0)*(x-x0)+y0

For more clarity you might define a function for linear interpolation
between two nodes. I use the following implementation that is more
symmetric. I hope it is more robust with respect to cancelations:

interpolateLinear :: Fractional a => (a,a) -> (a,a) -> a -> a
interpolateLinear (x0,y0) (x1,y1) x =
(y0*(x1-x) + y1*(x-x0))/(x1-x0)

> -- Calc
> y = interp xi yi x
>
> main = do
> 	-- Output Data
>    	print (y)

print y   is just fine, or    print (interp xi yi x)

```