ttaaa

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,
  5.         isThreaded,
  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
  26. isThreaded StraightLine = True
  27. isThreaded _ = False
  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