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

[Henning Thielemann]

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

Can you please add type signatures? This would help me understanding.


> 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)

see also utility-ht:Data.Ord.HT.limit
   http://hackage.haskell.org/packages/archive/utility-ht/0.0.5.1/doc/html/Data-Ord-HT.html


> 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)

(Taken from http://code.haskell.org/~thielema/htam/src/Numerics/Interpolation/Linear.hs)


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

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