# Interpolation

Jan 17th, 2022 (edited)
255
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. module Interpolation
2.     (
3.         InterpolationType,
4.         typeFromString,
6.         Interpolation,
7.         getL,
8.         getR,
9.         interpolate,
10.         createInterpolation
11.     ) where
12.
13. import Point
14. import Data.Text as Text
15.
16. data InterpolationType = StraightLine | LeastSquares deriving (Enum, Show, Eq)
17.
18. typeFromString :: String -> InterpolationType
19. typeFromString s =
20.     case Text.unpack . Text.toLower . Text.pack \$ s of
21.         "straight" -> StraightLine
22.         "leastsquares" -> LeastSquares
23.         _ -> error "Unknown method"
24.
25. isThreaded :: InterpolationType -> Bool
28.
29. data Interpolation = Func Double Double (Double -> Double)
30.
31. instance Show Interpolation where
32.     show (Func l r f) = "Interpolation from " ++ show l ++ " to " ++ show r
33.
34.
35. getL :: Interpolation -> Double
36. getL (Func l _ _) = l
37.
38. getR :: Interpolation -> Double
39. getR (Func _ r _) = r
40.
41. interpolate :: Interpolation -> Double -> Point
42. interpolate (Func l r f) x =
43.     if (l <= x) && (x <= r) then
44.         createPoint x (f \$ x)
45.     else
46.         error ("Interpolation defined only on range from " ++ show l ++ " to " ++ show r)
47.
48. merge :: Interpolation -> Interpolation -> Interpolation
49. merge (Func l1 r1 f1) (Func l2 r2 f2) =
50.     if (r1 /= l2 && r2 /= l1) then
51.         error ("Can't merge to Interpolation 1: " ++ show (Func l1 r1 f1) ++ " 2: " ++ show (Func l2 r2 f2))
52.     else let
53.         l = min l1 l2
54.         m = min r1 r2
55.         r = max r1 r2
56.         lf = if l == l1 then f1
57.         else f2
58.         mf = if l == l1 then f2
59.         else f1
60.         f = \x -> if (x <= m) then lf \$ x
61.                   else mf \$ x
62.     in Func l r f
63.
64. createInterpolation :: InterpolationType -> [Point] -> Interpolation
65. createInterpolation t p =
66.     case t of
67.         StraightLine -> createStraightLineInterpolation p
68.         LeastSquares -> createLeastSquaresInterpolation p
69.         _ -> error "Unimplimented method"
70.
71. createStraightLineInterpolation :: [Point] -> Interpolation
72. createStraightLineInterpolation (p1:p2:[]) =
73.     let x1 = getX p1
74.         x2 = getX p2
75.         y1 = getY p1
76.         y2 = getY p2
77.     in Func x1 x2 (\x -> (y1 + (x - x1) * (y2 - y1) / (x2 - x1)))
78. createStraightLineInterpolation (p1:tl@(p2:_)) =
79.     merge (createStraightLineInterpolation (p1:p2:[])) (createStraightLineInterpolation tl)
80. createStraightLineInterpolation _ = error "To interpolate function we need 2 points at least"
81.
82. createLeastSquaresInterpolation :: [Point] -> Interpolation
83. createLeastSquaresInterpolation p =
84.     let sumx = sum (Prelude.map getX p)
85.         sumy = sum (Prelude.map getY p)
86.         sumx2 = sum (Prelude.map (pow2 . getX) p)
87.         sumxy = sum (Prelude.map mulxy p)
88.
89.         n = fromIntegral (Prelude.length p) :: Double
90.         b = (sumxy * n - sumx * sumy) / (sumx2 * n - sumx * sumx)
91.         a = (sumy - b * sumx) / n
92.
93.         left = getX (Prelude.head p)
94.         right = getX (Prelude.head (Prelude.reverse p))
95.     in Func left right (linear a b)
96.
97. linear :: Double -> Double -> Double -> Double
98. linear a b x = a + b * x
99.
100. mulxy :: Point -> Double
101. mulxy p = (getX p) * (getY p)
102.
103. pow2 :: Double -> Double
104. pow2 x = x * x
105.
RAW Paste Data Copied