population-processing

population-processing.hs

{-# LANGUAGE TupleSections #-}
module Lib
    ( someFunc
    ) where

import Text.Parsec
import System.IO
import Data.List
import Graphics.Gnuplot.Simple
import Control.Monad.Reader

parseColumn :: Parsec String () (String, Float)
parseColumn = do
  year <- many1 $ satisfy (/= ',')
  char ','
  char '\"'
  population <- many1 $ satisfy (/= '\"')
  skipMany space
  char '\"'
  skipMany (char ',')
  char '\n'
  --skipMany space
  return (year, read $ filter (/=',') population)

parseCsv :: Parsec String () [(String, Float)]
parseCsv = do
  column <- parseColumn
  rest <- (eof >> return []) <|> (parseCsv)
  return (column:rest)

average :: [Float] -> Float
average xs = sum xs / fromIntegral(length xs)

chunk :: Int -> [a] -> [[a]]
chunk _ [] = []
chunk n xs = (take n xs) : (chunk n (drop n xs))

lstsAdder :: [[Float]] -> [Float]
lstsAdder [] = repeat 0
lstsAdder (x:xs) = zipWith (+) x $ lstsAdder xs

population :: ReaderT [(String, Float)] IO [Float]
population = do
  csv <- (ask :: ReaderT [(String, Float)] IO [(String, Float)])
  let _popu = map snd csv
  let result = _popu ++ replicate (12 - ((length _popu)`mod`12)) (last _popu)
  return result

averagePopulation :: ReaderT [(String, Float)] IO Float
averagePopulation = do
  p <- population
  return $ average p

getTrend ::  ReaderT [(String, Float)] IO [Float]
getTrend = do
  popu <- population
  let trend = map (average . take 12) $ unfoldr eachTail popu
  return trend

  where eachTail [] = Nothing
        eachTail (a:as) = Just (a:as, as)

getApril :: Int -> ReaderT [(String, Float)] IO [Float]
getApril aprilNum = do
  popu <- population
  let april = map (!! aprilNum) $ chunk 12 popu
  return april

getAprilGraph :: Int -> ReaderT [(String, Float)] IO [Float]
getAprilGraph aprilNum  = do
  april <- getApril aprilNum
  dummy <- averagePopulation
  let aprilGraph = replicate (aprilNum ) dummy ++ (concat . map (replicate 12) $ april)
  return aprilGraph

getAverageSlope :: ReaderT [(String, Float)] IO Float
getAverageSlope = do
  popu <- population
  let average_dif = average $ zipWith (flip (-)) popu (tail popu)
  return average_dif

getAprilSlope_av  :: Int -> ReaderT [(String, Float)] IO Float
getAprilSlope_av aprilNum = do
  april <- getApril aprilNum
  let aprilSlope_av = (average $ zipWith (flip (-)) april (tail april)) / 12
  return aprilSlope_av

getPeriod :: ReaderT [(String, Float)] IO [Float]
getPeriod = do
  popu <- population
  trend <- getTrend
  return $ zipWith (-) popu trend

getAvPeriod ::  ReaderT [(String, Float)] IO [Float]
getAvPeriod = do
  period <- getPeriod
  popu <- population
  let av_period = map (/ (fromIntegral (length popu)/12)) $ lstsAdder $ chunk 12 period
  return av_period

getAvPeriodGraph ::  ReaderT [(String, Float)] IO [Float]
getAvPeriodGraph = do
  popu <- population
  av_period <- getAvPeriod
  let av_period_graph = take (length popu) $ concat . repeat $ av_period
  return av_period_graph

deviance :: [Float] -> [Float]
deviance xs = let x_av = average xs
               in map ((-) x_av) xs

variance :: [Float] -> Float
variance xs = average $ map (^2) $ deviance xs

covariance :: [Float] -> [Float] -> Float
covariance xs ys = let
  x_dev = deviance xs
  y_dev = deviance ys
  in average $ zipWith (*) x_dev y_dev

linear_regression_slope :: [Float] -> [Float] -> Float
linear_regression_slope xs ys = covariance xs ys / variance xs

linear_regression_constant  :: [Float] -> [Float] -> Float
linear_regression_constant xs ys = average ys - (linear_regression_slope xs ys) * average xs


getA0 :: Int -> ReaderT [(String, Float)] IO Float
getA0 aprilNum = do
  av_period <- getAvPeriod
  trend <- getTrend
  a1 <- getA1
  let
    xs = map fromIntegral [0..((length trend)-1)] :: [Float]
    a0 =  linear_regression_constant xs trend
    --a0 = (av_period !! aprilNum) + (trend !! 0)
  return a0

getA1 :: ReaderT [(String, Float)] IO Float
getA1 = do
  trend <- getTrend
  let xs = map fromIntegral [0..((length trend)-1)]
  return (linear_regression_slope xs trend)
  --april <- (getAprilSlope_av 3)
  --return april

getA2 :: ReaderT [(String, Float)] IO Float
getA2 = do
  av_period <- getAvPeriod
  let
    av_period_max = foldl1 max av_period
    av_period_min = foldl1 min av_period
    a2 = av_period_max - av_period_min
  return a2

getApproximation :: Int -> ReaderT [(String, Float)] IO (Float -> Float)
getApproximation aprilNum = do
  a0 <- getA0 aprilNum
  a1 <- getA1
  a2 <- getA2
  let appro_period n = a2 * (sin (pi/12 * (n-(fromIntegral aprilNum))))^2
  let app_year_av n = a1 * n
  let approximation n = app_year_av n + appro_period n + a0 :: Float
  return approximation

myStyle = PlotStyle Points $ lineSpec (defaultStyle :: PlotStyle)

output :: ReaderT [(String, Float)] IO ()
output = do
  popu <- population
  let toList f = map f [0..(fromIntegral(length popu)-1)]
  average_popu <- averagePopulation
  trend <- getTrend
  let aprilNum = 3
  april <- getApril aprilNum
  aprilGraph <- getAprilGraph aprilNum
  -- average_slope <- getAverageSlope
  aprilSlope_av <- getAprilSlope_av aprilNum
  period <- getPeriod
  av_period <- getAvPeriod
  av_period_graph <- getAvPeriodGraph

  a0 <- getA0 aprilNum
  liftIO $ putStrLn $ "length av_period: " ++ show (length av_period)
  -- let a1 = aprilSlope_av
  a1 <- getA1
  a2 <- getA2
  let appro_period n = a2 * (sin (pi/12 * (n-(fromIntegral aprilNum))))^2 + (av_period !! aprilNum):: Float
  let app_year_av n = a1 * n + (trend !! 0)
  let approximation n = app_year_av n + appro_period n :: Float

  let errorRate = map (^2) $ zipWith (-) (toList approximation)  popu

  liftIO $ putStrLn $ "a0: " ++ show a0
  liftIO $ putStrLn $ "a1: " ++ show a1
  liftIO $ putStrLn $ "a2: " ++ show a2

  liftIO $ plotListsStyle [(Title "Total")] $ map (myStyle,) $
            [popu,
            --trend,
            aprilGraph,
            toList approximation
            ]
  {-plotListsStyle [(Title "Error")] $ map (myStyle,) $
    [ errorRate
    , map ((*1000) . flip (-) (average april)) aprilGraph]-}

  liftIO $ plotListsStyle [(Title "Period")] $ map (myStyle,) $
    [period,
    map (flip (-) (average april)) aprilGraph,
    toList appro_period,
    av_period_graph]

  liftIO $ plotListsStyle [(Title "Trend")] $ map (myStyle,) $
    [trend
    , toList (\x -> a0 + a1 * x)
    ]

  liftIO $ print $ sum errorRate


  -- plotPath [(Title "hello")] $ zip xs trend


{-run :: Parsec String () [(String, Float)] -> String -> IO ()
run p input = case parse p "hoge" input of
                   Left err -> putStr "parse error at" >> print err
                   Right csv -> (`runReaderT` csv) $ output -}

getPopu1 :: Parsec String () [(String, Float)] -> String -> IO [(String, Float)]
getPopu1 p input = case parse p "hoge" input of
                   Left err -> error $ show err
                   Right csv -> return csv

showPopulation :: [(String, Float)] -> IO ()
showPopulation sfs = forM_ sfs $ \(year, popu) -> do
  putStrLn $ year ++ ": " ++ show popu


someFunc :: IO ()
someFunc = do
  _nagoya_popu <- readFile "./nagoyashi_popu2.csv"
  _nagakute_popu <- readFile "./nagakute_popu.csv"
  _sapporo_popu <- readFile "./sapporo_shi2.csv"
  nagoya_popu1 <- getPopu1 parseCsv _nagoya_popu
  nagakute_popu1 <- getPopu1 parseCsv _nagakute_popu
  sapporo_popu1 <- getPopu1 parseCsv _sapporo_popu

  f <- runReaderT (getApproximation 3) nagoya_popu1
  g <- runReaderT (getApproximation 3) sapporo_popu1

  h <- runReaderT (getApproximation 3) $ zipWith (\(str, i) (str', j) -> (str++","++str', i+j)) nagoya_popu1 sapporo_popu1

  let f_g_list = map ((\x -> f x + g x) . fromIntegral) [1..length(sapporo_popu1)]
  let h_list = map ((\x -> h x) . fromIntegral) [1..length(sapporo_popu1)]

  print (f_g_list!!5)
  print (h_list!!5)

  plotListsStyle [(Title "Sum")] $ map (myStyle,) $
    [map ((\x -> h x) . fromIntegral) [1..length(nagoya_popu1)]
    , map ((\x -> f x + g x) . fromIntegral) [1..length(nagoya_popu1)]
    ]
  runReaderT output $ sapporo_popu1