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)
_______________________________________________
Haskell-Cafe mailing list
[email protected]
http://www.haskell.org/mailman/listinfo/haskell-cafe
