[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
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)
More information about the Haskell-Cafe
mailing list