Hi,
I've taken a look at your code and refactored a bit to use Haskell's
list functions. The functionality should be identical.
One of the advantages of using the list functions is that you don't
have to worry about an out-of-bounds exception - which is what your
limIndex function corrects for. Also instead of grabbing list indices
I've created a data structure for points.
I'll just present the code and let you ask me any questions.
-deech
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
type Point = ((Double,Double),(Double,Double))
getPnts :: Double -> [Point] -> Point
getPnts x [] = error "empty list"
getPnts x xs = case (break (\((i,_),_) -> i > x) xs) of
(_,[]) -> last xs
(_,xs) -> head xs
byTwos :: [a] -> [(a,a)]
byTwos [] = []
byTwos (x:[]) = []
byTwos (x:y:xs) = (x,y) : byTwos (y:xs)
interp :: [Double] -> [Double] -> Double -> Double
interp xi yi x = let ((x1,x2),(y1,y2)) = getPnts x $ zip (byTwos xi)
(byTwos yi)
in
(y2-y1)/(x2-x1) * (x-x1) + x2
main = print $ interp xi yi 11
On Tue, Jan 25, 2011 at 3:05 PM, gutti
Hi,
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
Cheers Phil
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
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)]
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
-- Calc y = interp xi yi x
main = do -- Output Data print (y)
-- View this message in context: http://haskell.1045720.n5.nabble.com/HMatrix-Vector-Matrix-interpolation-ala... Sent from the Haskell - Haskell-Cafe mailing list archive at Nabble.com.
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe