Interpolation

module Interpolation
    (
        InterpolationType,
        typeFromString,
        isThreaded,
        Interpolation,
        getL,
        getR,
        interpolate,
        createInterpolation
    ) where

import Point
import Data.Text as Text

data InterpolationType = StraightLine | LeastSquares deriving (Enum, Show, Eq)

typeFromString :: String -> InterpolationType
typeFromString s =
    case Text.unpack . Text.toLower . Text.pack $ s of
        "straight" -> StraightLine
        "leastsquares" -> LeastSquares
        _ -> error "Unknown method"

isThreaded :: InterpolationType -> Bool
isThreaded StraightLine = True
isThreaded _ = False

data Interpolation = Func Double Double (Double -> Double)

instance Show Interpolation where
    show (Func l r f) = "Interpolation from " ++ show l ++ " to " ++ show r


getL :: Interpolation -> Double
getL (Func l _ _) = l

getR :: Interpolation -> Double
getR (Func _ r _) = r

interpolate :: Interpolation -> Double -> Point
interpolate (Func l r f) x =
    if (l <= x) && (x <= r) then
        createPoint x (f $ x)
    else
        error ("Interpolation defined only on range from " ++ show l ++ " to " ++ show r)

merge :: Interpolation -> Interpolation -> Interpolation
merge (Func l1 r1 f1) (Func l2 r2 f2) =
    if (r1 /= l2 && r2 /= l1) then
        error ("Can't merge to Interpolation 1: " ++ show (Func l1 r1 f1) ++ " 2: " ++ show (Func l2 r2 f2))
    else let
        l = min l1 l2
        m = min r1 r2
        r = max r1 r2
        lf = if l == l1 then f1
        else f2
        mf = if l == l1 then f2
        else f1
        f = \x -> if (x <= m) then lf $ x
                  else mf $ x
    in Func l r f

createInterpolation :: InterpolationType -> [Point] -> Interpolation
createInterpolation t p =
    case t of
        StraightLine -> createStraightLineInterpolation p
        LeastSquares -> createLeastSquaresInterpolation p
        _ -> error "Unimplimented method"

createStraightLineInterpolation :: [Point] -> Interpolation
createStraightLineInterpolation (p1:p2:[]) =
    let x1 = getX p1
        x2 = getX p2
        y1 = getY p1
        y2 = getY p2
    in Func x1 x2 (\x -> (y1 + (x - x1) * (y2 - y1) / (x2 - x1)))
createStraightLineInterpolation (p1:tl@(p2:_)) =
    merge (createStraightLineInterpolation (p1:p2:[])) (createStraightLineInterpolation tl)
createStraightLineInterpolation _ = error "To interpolate function we need 2 points at least"

createLeastSquaresInterpolation :: [Point] -> Interpolation
createLeastSquaresInterpolation p =
    let sumx = sum (Prelude.map getX p)
        sumy = sum (Prelude.map getY p)
        sumx2 = sum (Prelude.map (pow2 . getX) p)
        sumxy = sum (Prelude.map mulxy p)

        n = fromIntegral (Prelude.length p) :: Double
        b = (sumxy * n - sumx * sumy) / (sumx2 * n - sumx * sumx)
        a = (sumy - b * sumx) / n

        left = getX (Prelude.head p)
        right = getX (Prelude.head (Prelude.reverse p))
    in Func left right (linear a b)

linear :: Double -> Double -> Double -> Double
linear a b x = a + b * x

mulxy :: Point -> Double
mulxy p = (getX p) * (getY p)

pow2 :: Double -> Double
pow2 x = x * x