Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Linear equations method.
module Main where
import Control.Monad.Fix (fix)
import Data.Bool (bool)
import Data.Tuple.Extra (first, second, dupe, (***))
import Data.Maybe (maybe)
import Graphics.Rendering.Chart.Backend.Cairo
import Graphics.Rendering.Chart.Easy
import System.IO (IOMode (..), withFile, hGetLine, hIsEOF)
import System.Exit (exitFailure)
readDta :: String -> IO (Maybe [(Double, Double)])
readDta fname = withFile fname ReadMode $ \h -> ($ []) . fix $ \f s ->
(>>=) (hIsEOF h) $ flip bool (return $ Just s) $ do
line <- words <$> hGetLine h
if length line == 2 then
f $ (uncurry (,) . first head . second last . dupe . map (read :: String -> Double)) line : s
else return Nothing
linearEquations :: [(Double, Double)] -> (Double, Double)
linearEquations p = ((len * sxy - sx * sy) / (len * sx2 - sx^^2), (sx2 * sy - sxy * sx) / (len * sx2 - sx^^2))
where
len = fromIntegral $ length p
(sxy, sx2) = foldr1 (uncurry (***) . both (+)) $ map (uncurry (,) . first (uncurry (*)) . second ((^^2) . fst) . dupe) p
(sx, sy) = foldr1 (uncurry (***) . both (+)) p
main :: IO ()
main = (>>=) (readDta "./mile.dta") $ maybe exitFailure $ \xy -> toFile def "linearEquations.png" $ do -- Reading this dataset: http://math.arizona.edu/~dsl/bmile.htm
let (a, b) = linearEquations xy
layout_title .= "The record of Mile run"
setColors [opaque red, opaque blue]
plot $ points "When the men's world record for the mile was broken this century" xy
plot $ line "Line fitted by least squares method" [[(x, a * x + b) | x <- [fst $ last xy .. fst $ head xy]]]
module Main where
import Control.Monad.Fix (fix)
import Data.Bool (bool)
import Data.Tuple.Extra (first, second, dupe, both, (***))
import Data.Maybe (maybe)
import qualified Graphics.Rendering.Chart.Backend.Cairo as GR
import qualified Graphics.Rendering.Chart.Easy as GR
import System.IO (IOMode (..), withFile, hGetLine, hIsEOF)
import System.Exit (exitFailure)
readDta :: String -> IO (Maybe [(Double, Double)])
readDta fname = withFile fname ReadMode $ \h -> ($ []) . fix $ \f s -> (>>=) (hIsEOF h) $ flip bool (return $ Just s) $ do
line <- words <$> hGetLine h
if length line == 2 then f $ (uncurry (,) . first head . second last . dupe . map (read :: String -> Double)) line : s else return Nothing
setMaxAbs :: (Num a, Ord a) => Int -> [[a]] -> [[a]]
setMaxAbs i = uncurry (++) . first (take i) . second (setMaxAbs' i . drop i) . dupe
where
swap i j row | i > j = swap j i row | otherwise = take i row ++ [row !! j] ++ take (j - i - 1) (drop (i + 1) row) ++ [row !! i] ++ drop (j + 1) row
setMaxAbs' ii row
| succ ii < length row = setMaxAbs' (succ ii) $ if abs ((row !! ii) !! ii) < abs ((row !! succ ii) !! ii) then swap ii (succ ii) row else row
| otherwise = row
forward :: (Fractional a, Ord a) => [[a]] -> [[a]]
forward = forward' 0
where
forward' n xs | length xs - 1 > n = forward' (succ n) $ toZero n $ setMaxAbs n xs | otherwise = xs
toZero i xs = take (i + 1) xs ++ [let c = (r !! i) / (ur !! i) in zipWith ((+) . negate . (*c)) ur r | r <- drop (i + 1) xs]
where
ur = xs !! i
backward :: (Fractional a, Ord a) => [[a]] -> [[a]]
backward = uncurry backward' . first ((+(-1)) . length) . dupe
where
backward' n xs | n >= 0 = backward' (pred n) $ toZero n xs | otherwise = xs
toZero i xs = [let ur = xs !! i; c = (r !! i) / (ur !! i) in zipWith ((+) . negate . (*c)) ur r | r <- take i xs] ++ drop i xs
cofactorExpand :: Num a => Int -> [[a]] -> a
cofactorExpand _ [[x]] = x
cofactorExpand _ [[x1, y1], [x2, y2]] = x1 * y2 - y1 * x2
cofactorExpand i mx = expand i mx
where
expand i mx = sum $ map (\(e, j) -> (-1)^(i+j) * e * cofactorExpand (succ i) (takeCofactor j mx)) $ zip (head mx) $ take (length mx) (iterate (+1) 0)
takeCofactor j mx = map (uncurry (++) . first (take j) . second (drop (j + 1)) . dupe) (drop 1 mx)
gaussianElim :: (Fractional a, Ord a) => [[a]] -> Maybe [[a]]
gaussianElim m | (length m > 1) && all ((>= length m) . length) m = Just $ backward $ forward m | otherwise = Nothing
determinant :: Num a => [[a]] -> Maybe a
determinant m | all ((==length m) . length) m = Just $ cofactorExpand 1 m | otherwise = Nothing
resolveLinearEq :: (Fractional a, Ord a) => [[a]] -> [a] -> Maybe [a]
resolveLinearEq = (.) (fmap (map (uncurry (/) . first last . second (sum . init) . dupe))) . (.) gaussianElim . zipWith (flip (.) (:[]) . (++))
linearEquations :: (Fractional a, Ord a) => [(a, a)] -> Maybe [a]
linearEquations p = let m = [[fromIntegral $ length p, sx, sx2], [sx, sx2, sx3], [sx2, sx3, sx4]] in
($ determinant m) $ maybe Nothing $ bool Nothing (resolveLinearEq m [sy, sxy, sx2y]) . (/=0)
where
sumPair = foldr1 (uncurry (***) . both (+))
(sxy, sx2y) = sumPair $ map (uncurry (,) . first (uncurry (*)) . second (uncurry (*) . first (^^2)) . dupe) p
(sx2, sx3) = sumPair $ map (first ((^^2). fst) . second ((^^3) . fst) . dupe) p
(sx, sy) = sumPair p
sx4 = sum $ map ((^^4) . fst) p
main :: IO ()
main = (>>=) (readDta "./range.dta") $ maybe exitFailure $ \xy -> flip (maybe exitFailure) (linearEquations xy) $ \[b, a0, a1] -> GR.toFile GR.def "linearEquations1.png" $ do
GR.layout_title GR..= "A study of the trajectories of projectiles"
GR.setColors [GR.opaque GR.red, GR.opaque GR.blue]
GR.plot $ GR.points "the trajectories" xy
GR.plot $ GR.line "Line fitted by least squares method" [[(x, a0 * x + a1 * x^^2 + b) | x <- [fst $ last xy .. fst $ head xy]]]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment