HMatrix Vector/Matrix interpolation ala Matlab interp/interp2 ??
Dear Haskellers, I'm looking for Vector and especially Matric interpolation ala: z = interp2 (xMatrix, yMatrix, zMatrix, x, y) - x and y can be single values or also vectors or matrices - indeally with the options nearest, linear, quadratic, qubic.. Any hope that there is something similar especially using the HMatrix matrices (Matrix Double). - If I have to code it any suggestions ? Cheers Phil -- 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.
Hi Phil, On 22/01/11 14:07, gutti wrote:
Dear Haskellers,
I'm looking for Vector and especially Matric interpolation ala:
z = interp2 (xMatrix, yMatrix, zMatrix, x, y)
- x and y can be single values or also vectors or matrices - indeally with the options nearest, linear, quadratic, qubic..
Any hope that there is something similar especially using the HMatrix matrices (Matrix Double). - If I have to code it any suggestions ?
Cheers Phil
I'm not sure if this is what you mean, as I'm not familiar with matlab, so apologies if this is totally a waste of time.. But I had a simple cubic interpolation implemented already, it wasn't too hard to lift it into matrices (there are many elevation permutations, I imagine), I used (fromLists ...stuff... toLists), but maybe there's a better way: {-# LANGUAGE NoMonomorphismRestriction #-} module Interpolate where import Numeric.LinearAlgebra import Data.List (zipWith5) -- cubic interpolation cubic t a b c d = let a1 = 0.5 * (c - a) a2 = a - 2.5 * b + 2.0 * c - 0.5 * d a3 = 0.5 * (d - a) + 1.5 * (b - c) in ((a3 * t + a2) * t + a1) * t + b -- boring manual lifting liftMatrix5 f a b c d e = let la = toLists a lb = toLists b lc = toLists c ld = toLists d le = toLists e in fromLists (zipWith5 (zipWith5 f) la lb lc ld le) -- test mt = (3><3) [0, 0.1 ..] ma = (3><3) [0, 1 ..] mb = (3><3) [0, 2 ..] mc = (3><3) [0, 3 ..] md = (3><3) [0, 4 ..] test = liftMatrix5 cubic mt ma mb mc md {- test output (3><3) [ 0.0, 2.1, 4.4 , 6.9, 9.6, 12.5 , 15.600000000000001, 18.9, 22.4 ] -} -- http://claudiusmaximus.goto10.org
Hi Claude, thanks a lot. - Your code looks interesting. I just don't fully get how it works: - are t a b c d points or curve parameters ? - how does lifting to matrix create a 1d spline to a 2d spline ? -- I don't see how it works Cheers Phil -- 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.
Hi Claude, thanks a lot. - Your code looks interesting. I just don't fully get how it works: - are t a b c d points or curve parameters ? - how does lifting to matrix create a 1d spline to a 2d spline ? -- I don't see how it works Cheers Phil -- 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.
On Sat, 22 Jan 2011, gutti wrote:
I just don't fully get how it works:
- are t a b c d points or curve parameters ? - how does lifting to matrix create a 1d spline to a 2d spline ? -- I don't see how it works
I think it interpolates cubically between four equidistant nodes, then it lifts the interpolation from scalar values to matrices. However, I would avoid interim lists, but just perform a zipMatrix5. I hope there is one ...
Hi Phil, On 22/01/11 23:13, gutti wrote:
- are t a b c d points or curve parameters ?
a b c d are points, t is the interpolation coefficient (between 0 and 1)
- how does lifting to matrix create a 1d spline to a 2d spline ? -- I don't see how it works
essentially, it creates a matrix of 1d splines, but now I see that this isn't what you wanted... for interpolated 2d matrix lookup, something like this, perhaps: -- using the cubic interpolator from earlier in the thread m @@>+ (x, y) = let (i, j) = (floor x, floor y) (s, t) = (x - fromIntegral i, y - fromIntegral j) cx j' = cubic s (m@@>(i-1,j')) (m@@>(i,j')) (m@@>(i+1,j')) (m@@>(i+2,j')) in cubic t (cx (j-1)) (cx j) (cx (j+1)) (cx (j+2)) test = let m = (16><16) [0 ..] n = 36 r = 5 (x0, y0) = (8, 8) in [ m @@>+ (x, y) | a <- [0 .. n - 1] , let a' = 2 * pi * fromIntegral a / fromIntegral n , let x = x0 + r * cos a' , let y = y0 + r * sin a' ] Claude -- http://claudiusmaximus.goto10.org
On Sun, 23 Jan 2011, Claude Heiland-Allen wrote:
essentially, it creates a matrix of 1d splines, but now I see that this isn't what you wanted...
for interpolated 2d matrix lookup, something like this, perhaps:
Interpolated matrix or vector lookup can of course be written as interpolation of sub-matrices or sub-vectors. For lazy matrices and vectors this would be almost as efficient. interp i v = vectorIndex (floor i) $ interpolateVectorSpace (fraction i) (Vector.take (n-3) $ Vector.drop 0 v) (Vector.take (n-3) $ Vector.drop 1 v) (Vector.take (n-3) $ Vector.drop 2 v) (Vector.take (n-3) $ Vector.drop 3 v) (Sorry for the many fictional functions.)
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.
On Tuesday 25 January 2011 22:05:34, gutti wrote:
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
limIndex xi idx = max 0 (min idx (length xi - 2)) not necessarily better, but different
getIndex xi x = limIndex xi (maybe (length xi) id (findIndex (>x) xi)-1)
maybe val id is the same as fromMaybe val
getPnts xi yi idx = [xi !! idx, xi !! (idx+1), yi !! idx, yi !! (idx+1)]
getPnts xi yi idx = case drop idx $ zip xi yi of ((x1,y1):(x2,y2):_) -> [x1,x2,y1,y2] _ -> error "Not enough points"
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 (x1:x2:y1:y2:_) = getPnts xi yi (getIndex xi x) in (y2 - y1)/(x2 - x1)*(x-x1) + y1
-- Calc y = interp xi yi x
main = do -- Output Data print (y)
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...
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)
Hi Henning, thanks for the code review -- reason that I don't use the type declaration a lot -- It causes trouble , because I don't yet fully understand it. When I declare what I think is right is fails - see Message at bottom -- so what's wrong ? By the way I just used lists so far - no arrays -- Why is !! inefficient ? - what is better - Cheers Phil 1 import Data.List 2 3 -- Input Data 4 xi :: [Double] 5 xi = [0 .. 10] 6 yi :: [Double] 7 yi = [2, 3, 5, 6, 7, 8, 9, 10 , 9, 8, 7] 8 x = 11 :: Double 9 10 -- Functions 11 limIndex :: [a] -> Int -> Int 12 limIndex xi idx 13 | idx < 0 = 0 14 | idx > (length xi)-2 = (length xi)-2 15 | otherwise = idx 16 17 getIndex :: [a] -> Double -> Int 18 getIndex xi x = limIndex xi (maybe (length xi) id (findIndex (>x) xi)-1) 19 20 getPnts :: [a] -> [a] -> Int -> [a] 21 getPnts xi yi idx = [xi !! idx, xi !! (idx+1), yi !! idx, yi !! (idx+1)] 22 23 interp :: [a] -> [a] -> Double -> Double 24 interp xi yi x = 25 let pts = getPnts xi yi (getIndex xi x) 26 in (pts!!3-pts!!2)/(pts!!1-pts!!0)*(x-pts!!0)+pts!!2 27 28 -- Calc 29 y = interp xi yi x 30 31 main = do 32 -- Output Data 33 print (y) === Compiler Error Message === interp_v4.hs:18:66: Couldn't match expected type `Double' against inferred type `a' `a' is a rigid type variable bound by the type signature for `getIndex' at interp_v4.hs:17:13 Expected type: [Double] Inferred type: [a] In the second argument of `findIndex', namely `xi' In the third argument of `maybe', namely `(findIndex (> x) xi)' interp_v4.hs:26:39: Couldn't match expected type `Double' against inferred type `a' `a' is a rigid type variable bound by the type signature for `interp' at interp_v4.hs:23:11 In the second argument of `(-)', namely `pts !! 0' In the second argument of `(*)', namely `(x - pts !! 0)' In the first argument of `(+)', namely `(pts !! 3 - pts !! 2) / (pts !! 1 - pts !! 0) * (x - pts !! 0)' -- 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.
On Tuesday 25 January 2011 23:16:49, gutti wrote:
Hi Henning,
thanks for the code review -- reason that I don't use the type declaration a lot -- It causes trouble , because I don't yet fully understand it.
When I declare what I think is right is fails - see Message at bottom -- so what's wrong ?
See inline.
By the way I just used lists so far - no arrays -- Why is !! inefficient ? - what is better - Cheers Phil
(!!) is inefficient because it has to traverse the list from the start to the desired index. That's fine if you do it once, but if you do it a lot, you're probably doing something wrong.
1 import Data.List 2 3 -- Input Data 4 xi :: [Double] 5 xi = [0 .. 10] 6 yi :: [Double] 7 yi = [2, 3, 5, 6, 7, 8, 9, 10 , 9, 8, 7] 8 x = 11 :: Double 9 10 -- Functions 11 limIndex :: [a] -> Int -> Int 12 limIndex xi idx 13 | idx < 0 = 0 14 | idx > (length xi)-2 = (length xi)-2 15 | otherwise = idx 16 17 getIndex :: [a] -> Double -> Int 18 getIndex xi x = limIndex xi (maybe (length xi) id (findIndex (>x) xi)-1)
The type of (>) is (>) :: Ord a => a -> a -> Bool , so both arguments to (>) must have the same type and that type must be an instance of the Ord class. Since the type of x is stated to be Double, (> x) :: Double -> Bool, so the list needs to have the type [Double]. You can ask ghci what the (most general) type of your function is, here it's getIndex :: Ord a => [a] -> a -> Int
19 20 getPnts :: [a] -> [a] -> Int -> [a] 21 getPnts xi yi idx = [xi !! idx, xi !! (idx+1), yi !! idx, yi !! (idx+1)] 22 23 interp :: [a] -> [a] -> Double -> Double 24 interp xi yi x = 25 let pts = getPnts xi yi (getIndex xi x) 26 in (pts!!3-pts!!2)/(pts!!1-pts!!0)*(x-pts!!0)+pts!!2
(-) :: Num a => a -> a -> a The two arguments of (-) must have the same type and the result has the same type too. x is stated to be a Double by the signature, so we must have pts :: [Double] and hence xi :: [Double], yi :: [Double]. For the most general type, the use of getIndex implies an Ord constraint, (-) and (*) give a Num constraint and (/) :: Fractional a => a -> a -> a , so the use of (/) adds a Fractional constraint. Fractional implies Num, hence interp :: (Ord a, Fractional a) => [a] -> [a] -> a -> a
27 28 -- Calc 29 y = interp xi yi x 30 31 main = do 32 -- Output Data 33 print (y)
=== Compiler Error Message ===
interp_v4.hs:18:66: Couldn't match expected type `Double' against inferred type `a' `a' is a rigid type variable bound by the type signature for `getIndex' at interp_v4.hs:17:13 Expected type: [Double] Inferred type: [a] In the second argument of `findIndex', namely `xi' In the third argument of `maybe', namely `(findIndex (> x) xi)'
interp_v4.hs:26:39: Couldn't match expected type `Double' against inferred type `a' `a' is a rigid type variable bound by the type signature for `interp' at interp_v4.hs:23:11 In the second argument of `(-)', namely `pts !! 0' In the second argument of `(*)', namely `(x - pts !! 0)' In the first argument of `(+)', namely `(pts !! 3 - pts !! 2) / (pts !! 1 - pts !! 0) * (x - pts !! 0)'
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
On Sat, 22 Jan 2011, gutti wrote:
I'm looking for Vector and especially Matric interpolation ala:
z = interp2 (xMatrix, yMatrix, zMatrix, x, y)
- x and y can be single values or also vectors or matrices - indeally with the options nearest, linear, quadratic, qubic..
Any hope that there is something similar especially using the HMatrix matrices (Matrix Double). - If I have to code it any suggestions ?
Here is the power of Haskell - I have implemented some interpolations in terms of Vector spaces, that is you can interpolate Vectors, Matrices, Functions, or whatever you like: http://code.haskell.org/synthesizer/core/src/Synthesizer/Interpolation/Modul... However this code is specialised to interpolate between adjacent objects in a stream, and needs numeric-prelude, and you have to declare an instance for Algebra.Module.C Double (HMatrix Double) and this would be an orphan instance that should better go to an official separate package ...
Hi Henning, thanks - I don't claim I understand that yet with my limited haskell knowledge, but it looks rather advanced passing modules around. If there isn't anything "standard" of the shelf (which is important feedback), than its probably the best for me to pull it up from scratch and go through the learning curve. Many thanks, Philipp -- 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.
participants (6)
-
aditya siram -
Claude Heiland-Allen -
Daniel Fischer -
Ferenc Wagner -
gutti -
Henning Thielemann