text:
-- |
-- Module      : Cartesian.Plane.Core
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, year
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created date year

-- TODO | - Which constraints are appropriate (Num is probably too generic, should be Real, maybe RealFrac)
--        - Strictness, performance
--        - Rename (?)

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Cartesian.Plane.Core where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
-- import Data.List (sort, minimumBy)
-- import Data.Ord  (comparing)
-- import Data.Complex hiding (magnitude)

-- import           Control.Monad (when)
import           Control.Applicative

import           Control.Lens ((^.))

import Cartesian.Internal.Types
import Cartesian.Internal.Lenses (begin, end)
import Cartesian.Internal.Core

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------

-- | Determines if a point lies within a polygon using the odd-even method.
--
-- TODO: Use epsilon (?)
-- TODO: How to treat points that lie on an edge
-- inside :: Num n => Polygon n -> Vector2D n -> Bool
-- inside polygon (Vector2D px py) = error "Cartesian.Plane.inside is still a work in progress"
--   where
--     edges = polygon ++ [head polygon] -- Close the loop
    -- between (Line (Vector ax ay) (Vector bx by)) = _

-- |
-- instance Convertible (Vector2D f, Vector3D f) where
  -- _

-- TODO: Use type families for this stuff (?)

-- |
-- to3D :: Num f => Vector2D f -> Vector3D f
-- to3D (Vector2D x' y') = Vector3D x' y' 0

-- -- |
-- from3D :: Num f => Vector3D f -> Vector2D f
-- from3D (Vector3D x' y' _) = Vector2D x' y'

-- -- | Perform some unary operation on a 2D vector as a 3D vector, converting the result back to 2D by discarding the z component.
-- -- TODO: Rename (?)
-- -- TODO: Loosen Num restriction (eg. to anything with a 'zero' value) (?)
-- in3D :: (Num f, Num f') => (Vector3D f -> Vector3D f') -> Vector2D f -> Vector2D f'
-- in3D f = from3D . f . to3D

-- | Same as in3D, but for binary operations.
-- _ :: _
-- _ f = from3D . (dotmap f) . to3D

-- Vector math -----------------------------------------------------------------------------------------------------------------------------

-- | Angle (in radians) between the positive X-axis and the vector
-- argument :: (Floating a, Eq a) => Vector a -> a
-- argument (Vector 0 0) = 0
-- argument (Vector x y) = atan $ y/x
--
--
-- arg :: (Floating a, Eq a) => Vector a -> a
-- arg = argument
--
--
-- -- | Vector -> (magnitude, argument)
-- polar :: (Floating a, Eq a) => Vector a -> (a, a)
-- polar v@(Vector x y) = (magnitude v, argument v)

-- Linear functions ------------------------------------------------------------------------------------------------------------------------

-- | Yields the intersection point of two finite lines. The lines are defined inclusively by
--   their endpoints. The result is wrapped in a Maybe value to account for non-intersecting
--   lines.
-- TODO: Refactor
-- TODO: Move equation solving to separate function (two linear functions)
-- TODO: Invariants, check corner cases
-- TODO: Deal with vertical lines
-- TODO: Factor out infinite-line logic
-- TODO: Decide how to deal with identical lines
-- TODO: Factor out domain logic (eg. write restrict or domain function)
-- TODO: Return Either instead of Maybe (eg. Left "parallel") (?)
-- TODO: Visual debugging functions
-- TODO: Math notes, MathJax or LaTex
-- TODO: Intersect for curves (functions) and single points (?)
-- TODO: Polymorphic, typeclass (lines, shapes, ranges, etc.) (?)
-- TODO: Intersect Rectangles
intersect :: RealFloat f => Line (Vector2D f) -> Line (Vector2D f) -> Maybe (Vector2D f)
intersect f' g' = do
  p <- mp
  indomain f' p
  indomain g' p
  where
    -- indomain :: RealFloat f => Line (Vector2D f) -> Vector2D f -> Maybe (Vector2D f)
    indomain h' = restrict (h'^.begin) (h'^.end) -- TODO: Rename

    -- mp :: Maybe (Vector2D f)
    mp = case (linear f', linear g') of
      (Just f, Nothing) -> let x' = g'^.begin.x in Just . Vector2D x' $ plotpoint f x'
      (Nothing, Just g) -> let x' = f'^.begin.x in Just . Vector2D x' $ plotpoint g x'
      (Just f,  Just g) -> linearIntersect f g
      _                 -> Nothing

-- | Gives the linear function overlapping the given segment, or Nothing if there is no such function
linear :: RealFloat f => Line (Vector2D f) -> Maybe (Linear f)
linear line = Linear <$> intercept line <*> slope line

-- | Applies a linear function to the given value
-- TODO: Rename (?)
plotpoint :: RealFloat f => Linear f -> f -> f
plotpoint f x' = slopeOf f*x' + interceptOf f

-- | Finds the intersection (if any) of two linear functions
-- TODO: Use Epsilon (?)
-- TODO: Rename (eg. 'solve') (?)
linearIntersect :: RealFloat f => Linear f -> Linear f -> Maybe (Vector2D f)
linearIntersect f g
  | slopeOf f == slopeOf g = Nothing
  | otherwise = let x' = (β-b)/(a-α) in Just $ Vector2D x' (a*x' + b)
  where
    (a, α) = (slopeOf     f, slopeOf     g)
    (b, β) = (interceptOf f, interceptOf g)

-- |
slope :: RealFloat f => Line (Vector2D f) -> Maybe f
slope (Line fr to)
  | dx == 0   = Nothing
  | otherwise = Just $ dy/dx
  where
    (Vector2D dx dy) = liftA2 (-) to fr

-- |
intercept :: RealFloat f => Line (Vector2D f) -> Maybe f
intercept line = do
  slope' <- slope line
  return $ y' - slope'*x'
  where
    (Vector2D x' y') = line^.begin

--------------------------------------------------------------------------------------------------------------------------------------------

-- | Ensures that a given point lies within the domain and codomain
-- TODO: Make polymorphic
-- TODO: Let this function work on scalars, write another function for domain and codomain (?)
-- restrict domain codomain p = _
-- restrict :: (Applicative v, Traversable v, Num n, Ord n) => v n -> v n -> v n -> Maybe (v n)
-- restrict a b p@(Vector2D x' y')
--   | indomain && incodomain = Just p
--   | otherwise              = Nothing
--   where
--     (Vector2D lowx lowy)   = dotwise min a b
--     (Vector2D highx highy) = dotwise max a b
--     indomain   = between lowx highx x'
--     incodomain = between lowy highy y'

max_stars_repo_name: jordanemedlock/Cartesian

max_issues_repo_name: jordanemedlock/Cartesian

max_forks_repo_name: jordanemedlock/Cartesian


text:
{-|
Module      : ExponentialSum
Description : Compute the exponential sum of the day

Compute the exponential sum of the day, to be plotted later.
 -}

module ExponentialSum
  (
    -- * Main function
    expSums
    -- * Exported for convenience
  , expVal
  , realArg
  ) where

import Data.Complex

-- | Compute the value of the real argument of the exponential, at point n,
-- given the initial coefficients.
--
-- For example, here is the third value computed where the coefficients are
-- current date:
--
-- >>> realArg [12, 5, 2018] 3
-- 2.063379583746283
-- >>> realArg [12, 5, 18] 3
-- 3.55
--
-- There are some properties which might be useful:
--
-- 1. If the coefficients are missing, we want to return a 0 for any index:
--
--     prop> realArg [] x == 0
--
-- 2. If the index is 0, we output 0 regardless of the coefficients, provided
-- at least one is non-zero
--
--     prop> a == 0 || realArg (a:xs) 0 == 0
--
-- There might be a few more properties but we can ignore them for now, as
-- some might fail due to numerical inaccuracies. We are also sidestepping
-- NaNs and infinities at the moment.
realArg :: [Double] -> Double -> Double
realArg xs n = sum $ zipWith power xs [1..]
  where
    power :: Double -> Int -> Double
    power x i = n ^ i / x

-- | Computes the value of the complex exponential for a term at an index
-- Since we need to multiply the result of `realArg` with i, we just put it in
-- the imaginary part of the newly minted complex number.
--
-- Examples:
--
-- >>> expVal [12, 5, 2018] 3
-- 0.9217505006686714 :+ 0.3877834634393964
-- >>> expVal [12, 5, 18] 3
-- (-0.9510565162951544) :+ (-0.3090169943749449)
--
-- Also, we have the following properties:
--
-- prop> x == 0 || expVal [] x == 1
-- prop> a >= 0 || expVal (a:xs) 0 == 1
expVal :: [Double] -> Double -> Complex Double
expVal xs n = exp $ 2 * pi * (0 :+ realArg xs n)

-- | Computes the partial sums of the exponential sum of the day
--
-- >>> take 3 $ expSums [12, 5, 18]
-- [1.0 :+ 0.0,0.4700807357667952 :+ 0.8480480961564261,(-0.3779673603896302) :+ 1.3779673603896319]
-- >>> take 3 $ expSums [12, 5, 2018]
-- [1.0 :+ 0.0,0.7890437903507289 :+ 0.9774955128339019,1.7720662136902803 :+ 0.794010049617778]
-- >>> take 3 $ expSums []
-- [1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0]
expSums :: [Double] -> [Complex Double]
expSums xs = scanl1 (+) $ map (expVal xs) [0..]

ext: hs

max_stars_repo_name: mihaimaruseac/esotd

max_issues_repo_name: mihaimaruseac/esotd

max_forks_repo_name: mihaimaruseac/esotd


text:
{-# LANGUAGE  OverloadedStrings #-}
{-# LANGUAGE  DeriveGeneric #-}
{-# LANGUAGE  DeriveAnyClass #-}
{-# LANGUAGE  NamedFieldPuns #-}
{-# LANGUAGE  TypeSynonymInstances #-}
{-# LANGUAGE  FlexibleInstances #-}
{-# LANGUAGE  QuasiQuotes #-}
{-# LANGUAGE  StandaloneDeriving #-}
{-# LANGUAGE  StrictData #-}
{-# LANGUAGE  TupleSections #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE QuantifiedConstraints #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE LambdaCase #-}

module ImdbStats where

import qualified Data.Map.Strict               as Map
import qualified Data.Text                     as Text
import qualified Data.Text.Lazy                 as Text.Lazy
import qualified Data.Text.Read                as Text.Read
import           Text.HTML.Scalpel
import           Control.Monad
import           Data.Char                      ( isDigit )
import           Data.List.Extra                ( enumerate )
import           Data.Function                  ( (&) )
import           Debug.Trace
import           Data.List.Extra                ( intercalate )
import           Data.Maybe                     ( catMaybes
                                                , fromMaybe
                                                )
import qualified Network.URL                   as URL
import           Data.String.Interpolate
import           GHC.Generics                   ( Generic )
import qualified Data.Aeson                    as Aeson
import qualified Data.Aeson.Text                 as Aeson.Text
import           Data.Aeson.Types               ( Parser )
import           Control.Monad.Extra            ( fromMaybeM )
import           Control.Monad.Fail             ( MonadFail )
import qualified Data.Vector                   as Vector
import           Data.Vector                    ( Vector )
import qualified Data.ByteString               as BS
import qualified Data.ByteString.Lazy          as LBS
import qualified Data.ByteString.Char8         as BS.Char8
import qualified Data.Attoparsec.ByteString.Char8
                                               as AP.BS
import           Control.DeepSeq                ( NFData )
import qualified Streamly as S
import qualified Streamly.Prelude as S
import qualified Streamly.Data.Fold as FL
import qualified Streamly.Memory.Array as A
import qualified Streamly.FileSystem.Handle as FH
import qualified Streamly.Data.String as DS
import qualified Data.Word as Word
import System.IO
import qualified Statistics.Test.StudentT as Stat
import qualified Statistics.Test.Types as Stat
import qualified Statistics.Sample as Stat
import Control.Applicative
import qualified Graphics.Rendering.Chart.Easy as Chart
import Graphics.Rendering.Chart.Easy ((.=))
import qualified Graphics.Rendering.Chart.Backend.Diagrams as Chart
import qualified Network.HTTP.Req as Req
import Network.HTTP.Req ((/~), (/:), (=:))
import Data.Default
import qualified Data.List.Split as Split
import qualified Control.Retry as Retry

data Age = AllAges | LT18 | A18to29 | A30to44 | GT45 deriving (Show, Eq, Ord, Bounded, Enum, Generic, NFData, Aeson.ToJSON, Aeson.FromJSON)
data Gender = AllGenders | Male | Female deriving (Show, Eq, Ord, Bounded, Enum, Generic, NFData, Aeson.ToJSON, Aeson.FromJSON)

type Demographic = (Age, Gender)
-- TODO: rational or scientific instead of double
newtype Rating = Rating Double deriving (Show, Generic, NFData, Aeson.ToJSON, Eq, Aeson.FromJSON)
unRating (Rating x) = x -- TODO: so as not to mess with the avaliable JSON
newtype NumUsers = NumUsers Int deriving (Show, Generic, NFData, Aeson.ToJSON, Eq, Aeson.FromJSON)
data RatingStats = RatingStats {rating :: Rating, numUsers :: NumUsers} deriving (Generic, NFData, Aeson.ToJSON, Eq, Aeson.FromJSON)

instance Show RatingStats where
    show (RatingStats (Rating rating) (NumUsers numUsers)) =
        (show rating) ++ "@" ++ (show numUsers)

type RatingTable = Map.Map Demographic RatingStats

instance {-# OVERLAPS #-} Show RatingTable where
    show table = Map.toList table & map show & intercalate "\n"

a @. b = a @: [hasClass b]

tableToMap :: [[RatingStats]] -> RatingTable
tableToMap table = Map.fromList $ zip
    [ (age, gender)
    | gender <- enumerate :: [Gender]
    , age <- enumerate :: [Age]
    ]
    (mconcat table)

-- TODO: operator (//) implies arbitrarly deep nesting, i need direct nesting (children)

scrapeRatingTable :: Scraper BS.ByteString RatingTable
scrapeRatingTable =
    chroots ("td" @. "ratingTable") cell
        & chroots "tr"
        & fmap (tableToMap . tail) -- skip the first row, it's just headings (it will be an empty list)
        & chroots ("div" @: ["class" @= "title-ratings-sub-page"] // "table")
        >>= \case -- the second table is the relevant one
                _ : x : _ : [] -> return x
                _ -> fail "table isn't on the second spot"
  where
    cell :: Scraper BS.ByteString RatingStats
    cell = do
        rating <- text $ "div" @. "bigcell"
        numUsers <- text $ "div" @. "smallcell" // "a"

        -- TODO: maybe a monad transformer for this?
        case
                ( AP.BS.parseOnly AP.BS.double rating
                , BS.Char8.readInt $ BS.filter AP.BS.isDigit_w8 numUsers
                )
            of
                (Right parsedRating, Just (parsedNumUsers, _)) ->
                    return $ RatingStats (Rating parsedRating)
                                         (NumUsers parsedNumUsers)
                _ -> fail "can't parse cell"

scrapeMoviesFromListPage :: Scraper BS.ByteString (Vector URL.URL)
scrapeMoviesFromListPage =
    ("div" @. "lister-list" // "h3" @. "lister-item-header" // "a")
        & attrs "href"
        & fmap
              (Vector.fromList . catMaybes . fmap
                  (URL.importURL . BS.Char8.unpack)
              )

-- TODO: scrape the number of movies in a list, calculate the number of pages, scrape them in parallel
scrapeWholeList :: String -> IO (Vector URL.URL)
scrapeWholeList url = rec url 1 Vector.empty
  where
    rec :: String -> Int -> Vector URL.URL -> IO (Vector URL.URL)
    rec url page collected = do
        Just links <- scrapeURL
            (url ++ "?sort=list_order,asc&st_dt=&mode=detail&page=" ++ show page)
            scrapeMoviesFromListPage
        if Vector.null links
            then -- stop iterating when scraping fails (on 404)
                 return collected
            else rec url (page + 1) (collected <> links)

loadLinks :: FilePath -> IO [URL.URL]
loadLinks path = do
    Just links <- Aeson.decode <$> LBS.readFile path -- TODO mtl, handling errors?
    return links

newtype MovieId = MovieId { getMovieId :: Text.Text } deriving (Show)

extractId :: URL.URL -> MovieId
extractId = MovieId . Text.pack . (!! 1) . Split.splitOn "/" . URL.url_path

downloadAllMovies :: FilePath -> [MovieId] -> IO ()
downloadAllMovies folder =
    S.drain
    . S.mapM (print . fst)
    . S.indexed
    . S.asyncly
    . S.mapM saveMovie
    . S.maxRate 10
    . S.mapM downloadMovie
    . S.fromList
    where
        downloadMovie :: MovieId -> IO (MovieId, LBS.ByteString)
        downloadMovie movieId =
            Req.runReq
                def { Req.httpConfigRetryPolicy = Retry.retryPolicy $ const $ Just 1000 }
             $  (movieId,)
            <$> Req.responseBody
            <$> Req.req
                Req.GET
                (Req.https "www.imdb.com" /: "title" /: (getMovieId movieId) /: "ratings")
                -- TODO unite with proper url making
                Req.NoReqBody
                Req.lbsResponse
                ("ref_" =: ("tt_ov_rt" :: Text.Text))

        saveMovie :: (MovieId, LBS.ByteString) -> IO ()
        saveMovie (movieId, body) =
            LBS.writeFile [i|#{folder}/#{getMovieId movieId}.html|] body -- TODO: path package

badRequestTest :: IO LBS.ByteString
badRequestTest =
    Req.runReq
        def -- { Req.httpConfigRetryPolicy = Retry.retryPolicy $ const $ Just 1000 }
    $ Req.responseBody
    <$> Req.req
        Req.GET
        (Req.https "httpbin.org" /: "status" /: (Text.pack $ show 504))
        -- TODO unite with proper url making
        Req.NoReqBody
        Req.lbsResponse
        mempty

data Movie = Movie { movieURL :: URL.URL, ratingTable :: RatingTable }
    deriving (Show, Generic, NFData, Aeson.ToJSON, Aeson.FromJSON)

moviePageToRatingsPage :: URL.URL -> URL.URL
moviePageToRatingsPage (URL.URL url_type url_path url_params) =
    URL.URL url_type (url_path ++ "ratings") []

addImdbHost :: URL.URL -> URL.URL
addImdbHost url = url
    { URL.url_type = URL.Absolute
                         $ URL.Host (URL.HTTP True) "www.imdb.com" Nothing
    }

getMovie :: URL.URL -> IO (Maybe Movie)
getMovie url = do
    let ratingsURL = addImdbHost $ moviePageToRatingsPage url
    maybeMovie <- scrapeURL (URL.exportURL ratingsURL) scrapeRatingTable
    print url
    return (fmap (Movie ratingsURL) maybeMovie)

getAllMovies :: [URL.URL] -> S.AsyncT IO Movie
getAllMovies = S.mapMaybeM getMovie . S.fromList

getAllMoviesIO :: [URL.URL] -> IO [Movie]
getAllMoviesIO = S.toList . S.asyncly . getAllMovies

streamEncodeArray :: (S.MonadAsync m, Aeson.ToJSON a) => S.SerialT m a -> S.SerialT m Char
streamEncodeArray =
    S.cons '['
    . (flip (<>)) (S.yield ']')
    . S.concatMap (S.fromList . Text.Lazy.unpack)
    . S.intersperse ","
    . S.map Aeson.Text.encodeToLazyText

encodeMoviesStr :: [URL.URL] -> IO LBS.ByteString
encodeMoviesStr =
    (fmap LBS.pack) .S.toList . DS.encodeChar8 . streamEncodeArray . S.asyncly . getAllMovies

writeStream :: FilePath -> S.Serial Char -> IO ()
writeStream file stream = do
    fh <- openFile file WriteMode
    S.fold (FH.write fh) $ DS.encodeChar8 stream

instance Aeson.ToJSON URL.URL where
    toJSON = Aeson.String . Text.pack . URL.exportURL

instance Aeson.FromJSON URL.URL where
    parseJSON (Aeson.String str) = lift $ URL.importURL $ Text.unpack str
      where
-- TODO: there's gotta be a function that does this. mtl?
        lift :: (MonadFail m) => Maybe a -> m a
        lift (Just x) = return x
        lift Nothing = fail "got Nothing while trying to lift"
    parseJSON _ = fail "URL has to be a string"

deriving instance Generic URL.URL
instance NFData URL.URL
deriving instance Generic URL.URLType
instance NFData URL.URLType
deriving instance Generic URL.Host
instance NFData URL.Host
deriving instance Generic URL.Protocol
instance NFData URL.Protocol

loadMovies :: IO (Maybe (Vector Movie))
loadMovies = Aeson.decode <$> LBS.readFile "data/ratings_all.json"

mean v = Vector.sum v / (fromIntegral $ Vector.length v)

paired = Vector.mapMaybe $ \movie ->
    liftA2 (,) (genderedRating Male movie) (genderedRating Female movie)
    where genderedRating gender = ((fmap (unRating . rating)) . Map.lookup (AllAges, gender) . ratingTable)

genderDiff :: IO ()
genderDiff = do
    Just movies <- loadMovies
    print $ Stat.pairedTTest Stat.BGreater $ paired movies

plotGenderDiff :: IO ()
plotGenderDiff = do
    Just movies <- loadMovies
    Chart.toFile Chart.def "plots/scatter.png" $ do
        Chart.layout_title .= "Male x Female ratings"
        Chart.setColors [Chart.opaque Chart.blue]
        Chart.plot (Chart.points "points" $ Vector.toList $ paired movies)

data Range r
    = FromTo r r
    | LessThan r
    | MoreThan r

instance Show r => Show (Range r) where
    show (FromTo a b) = [i|#{a} - #{b}|]
    show (LessThan x) = [i|< #{x}|]
    show (MoreThan x) = [i|> #{x}|]

genderDiffHistogram :: IO ()
genderDiffHistogram = do
    Just movies <- loadMovies
    let diff = Vector.map (uncurry (-)) $ paired movies
        sd = Stat.stdDev diff
        m = Stat.mean diff

    putStrLn "hey"

max_stars_repo_name: SandiCat/imdb-stats

max_issues_repo_name: SandiCat/imdb-stats

max_forks_repo_name: SandiCat/imdb-stats


text:
-- |
-- Module      : Occlusion.Core
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created September 20 2015

-- TODO | - Angles utilities, units, normalise angles
--        -

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------
-- {-# LANGUAGE ScopedTypeVariables #-}

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Occlusion.Core where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.Complex
import Data.Function
import Data.Fixed
import Data.List (minimumBy, maximumBy)
import Data.Ord  (comparing)

import Southpaw.Math.Constants
import Southpaw.Utilities.Utilities (pairwise)

import Occlusion.Types
import Occlusion.Lenses
import Occlusion.Vector

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------
-- |
angles :: RealFloat f => Complex f -> Polygon f -> [f]
angles = map . angle

-- |
angle :: RealFloat f => Complex f -> Complex f -> f
angle a = snd . polar . subtract a

-- | Finds the minimum and maximum value in a list, by some arbitrary criterion.
-- TODO: Strictness, performance, move to library
minmaxBy :: (a -> a -> Ordering) -> [a] -> Maybe ((Int, a), (Int, a))
minmaxBy _ []     = Nothing
minmaxBy f (x:xs) = Just . foldr (\n (mini, maxi) -> (minBy f n mini, maxBy f n maxi)) ((0, x), (0, x)) $ zip [1..] xs
  where
    minBy f n mini = maximumBy (f `on` snd) [n, mini]
    maxBy f n maxi = minimumBy (f `on` snd) [n, maxi]

-- |
-- TODO: Rename (?)
-- TODO: Maybe it would be a good idea if a function called 'anglespan' actually returned some angles.
-- TODO: Think radar
anglespan :: RealFloat f => Complex f -> Polygon f -> Maybe ((Int, Complex f), (Int, Complex f))
anglespan p shape = minmaxBy (comparing $ normalise . angle p) shape

-- |
normalise :: RealFloat f => f -> f
normalise = flip mod' $ 2*π

-- |
-- TODO: Rename (eg. occlusionEdge, frontEdge, near/far, cover, etc.) (?)
nearestEdge :: RealFloat f => Complex f -> Polygon f -> Edge f
nearestEdge p poly = error "Function 'nearestEdge' is currently on holiday"

-- |
-- TODO: Rename (eg. distant, etc.) (?)
-- TODO: Rigorous algorithm
-- TODO: I could solve this with intersect testing...
-- TODO: Or comparing (ai < bi) to (α < β)
distantEdge :: (RealFloat f, Ord f) => Complex f -> Polygon f -> Maybe (Edge f)
distantEdge p shape = case span' of
  -- Just ((ai, α), (bi, β)) -> Just $ slice (min ai bi) (max ai bi) shape -- TODO: This line needs some love and attention
  -- Just ((ai, α), (bi, β)) -> Just $ slice (max ai bi) (max ai bi+length shape - min ai bi) $ cycle shape -- TODO: This line needs some love and attention
  Just ((ai, a), (bi, b)) -> Just $ slice (max ai bi) (length shape + max ai bi - 1) $ cycle shape
  -- Just ((ai, fr), (bi, to)) -> Just $ if (ai < bi) == (normalise (angle p fr) < normalise (angle p to))
  --                                       then slice (max ai bi) (min ai bi + length shape) (cycle shape)
  --                                       else slice (min ai bi) (max ai bi) (shape)
  Nothing                 -> Nothing
  where
    span' = anglespan p shape
    -- gap   =
    -- from = _
    -- to   = _

-- |
-- TODO: Cyclic slice
slice :: Int -> Int -> [a] -> [a]
slice fr to = take (to-fr) . drop fr

--------------------------------------------------------------------------------------------------------------------------------------------

-- |
-- TODO: Cyclic (?)
-- TODO: Caching, performance, drop path segments as we walk past them
-- TODO: Arbitrary paths (not just straight lines)
-- TODO: Refactor, simplify
-- TODO: Invariants, tests, QuickCheck
walkalong :: RealFloat f => [Complex f] -> f -> Maybe (Complex f)
walkalong []   _         = Nothing
walkalong path progress' = case which of
  []     -> Nothing                    -- We've walked off the beaten path and now we're lost.
  (pl:_) -> Just $ uncurry walkline pl --
  where
    progress = mod' progress' (sum lengths) -- TODO: Remove this
    lengths  = pairwise distance path       -- The length of each segment, in order
    segments = pairwise (,) path            -- A list of consecutive endpoints
    which    = dropWhile ((<progress) . snd) $ zip segments (scanl1 (+) lengths)           -- Drops segments that have already been passed (leaving us with the line we're on as the first value)
    distance = (realPart . abs) .: (-)                                                     -- Distance (as a real number) between two points
    walkline (fr, to) l = fr + (to-fr) * (((progress-l+distance fr to)/distance fr to):+0) --
    (.:) f g = (f .) . g

max_stars_repo_name: SwiftsNamesake/Occlusion

max_issues_repo_name: SwiftsNamesake/Occlusion

max_forks_repo_name: SwiftsNamesake/Occlusion


text:
{-|
Module: Math.Ftensor.Algebra
Copyright: (c) 2015 Michael Benfield
License: ISC

Classes supporting arithmetic operations. Haskell's numeric classes are too
coarse for our needs. Types which are instances of more than one of these
classes should satisfy the obvious compatibility relations between the various
operations.

For instance, the equalities @x -. x = zero@ and @x *. (y +. z) = x *. y +. x
*. z@ and @one *. x = x@ should hold. If the type in question implements
inexact arithmetic, they should hold approximately.

For an instance of @Num@, these should hold whenever the type is also an
instance of the classes defining these functions:

* @(+.) = (+)@

* @(*.) = (*)@

* @(-.) = (-)@

* @neg = negate@

* @one = fromInteger 1@

* @zero = fromInteger 0@

Similarly, for an instance of @Fractional@

* @(\/.) = (\/)@

* @inv = recip@
-}

{-# LANGUAGE CPP #-}

module Math.Ftensor.Algebra (
    Additive(..),
    WithZero(..),
    WithNegatives(..),
    Multiplicative(..),
    WithOne(..),
    WithReciprocals(..),
    WithScalars(..),
    Rg,
    Rng,
    Rig,
    Ring,
    Field,
    Module,
    VectorSpace,
    Algebra,
    UnitalAlgebra,
) where

import Data.Complex (Complex)
import Data.Int
import Data.Ratio (Ratio)
import Data.Word
import Numeric.Natural

infixl 6 +.

class Additive a where
    (+.) :: a -> a -> a

-- | Types that have an additive identity. Should satisfy:
--
-- * @x +. zero = zero +. x = x@
class Additive a => WithZero a where
    zero :: a

infixl 6 -.

-- | Types that have additive inverses. Should satisfy:
-- @neg x = zero -. x@ and @x -. x = zero@.
class WithZero a => WithNegatives a where
    {-# MINIMAL neg | (-.) #-}
    neg :: a -> a
    neg x = zero -. x
    (-.) :: a -> a -> a
    lhs -. rhs = lhs +. neg rhs

infixl 7 *.

class Multiplicative a where
    (*.) :: a -> a -> a

-- | Types with a multiplicative identity. Should satisfy:
--
-- * @one *. x = x *. one = x@.
class Multiplicative a => WithOne a where
    one :: a

infixl 7 /.

-- | Types with multiplicative inverse.
-- @inv x@ and @y /. x@
-- may throw an exception or behave in undefined ways if
-- @invertible x = False@.
--
-- * @inv x = one ./ x@
--
-- * @y *. x /. y = x@
--
-- * @x /. y *. x = y@
class WithOne a => WithReciprocals a where

    {-# MINIMAL invertible, (inv | (/.)) #-}
    invertible :: a -> Bool

    inv :: a -> a
    inv x = one /. x

    (/.) :: a -> a -> a
    lhs /. rhs = lhs *. inv rhs

infixl 7 *:

-- | Types like mathematical vectors that can be multiplied by another type.
class WithScalars a where
    type Scalar a
    (*:) :: Scalar a -> a -> a

-- | A @Rg@ is a Ring without one or negatives.
type Rg a = (WithZero a, Multiplicative a)

-- | A @Rng@ is a Ring without one.
type Rng a = (Rg a, WithNegatives a)

-- | A @Rig@ is a Ring without negatives.
type Rig a = (Rg a, WithOne a)

type Ring a = (Rng a, Rig a)

type Field a = (Ring a, WithReciprocals a)

type Module a = (WithNegatives a, WithScalars a, Ring (Scalar a))

type VectorSpace a = (Module a, Field (Scalar a))

type Algebra a = (VectorSpace a, Rng a)

type UnitalAlgebra a = (Algebra a, WithOne a)

#define INSTANCES_NUM(typ, ctxt) \
instance ctxt Additive (typ) where { \
; {-# INLINE (+.) #-} \
; (+.) = (+) \
} ; instance ctxt WithZero (typ) where { \
; {-# INLINE zero #-} \
; zero = 0 \
} ; instance ctxt WithNegatives (typ) where { \
; {-# INLINE neg #-} \
; neg = negate \
; {-# INLINE (-.) #-} \
; (-.) = (-) \
} ; instance ctxt Multiplicative (typ) where { \
; {-# INLINE (*.) #-} \
; (*.) = (*) \
} ; instance ctxt WithOne (typ) where { \
; {-# INLINE one #-} \
; one = 1 \
} ; instance ctxt WithScalars (typ) where { \
; type Scalar (typ) = typ \
; {-# INLINE (*:) #-} \
; (*:) = (*) \
}

#define INSTANCES_FRACTIONAL(typ, ctxt) \
INSTANCES_NUM(typ, ctxt) ; \
instance ctxt WithReciprocals (typ) where { \
; {-# INLINE invertible #-} \
; invertible = (/= 0) \
; {-# INLINE inv #-} \
; inv = recip \
; {-# INLINE (/.) #-} \
; (/.) = (/) \
}

INSTANCES_NUM(Natural, )
INSTANCES_NUM(Integer, )
INSTANCES_NUM(Int, )
INSTANCES_NUM(Int8, )
INSTANCES_NUM(Int16, )
INSTANCES_NUM(Int32, )
INSTANCES_NUM(Int64, )
INSTANCES_NUM(Word, )
INSTANCES_NUM(Word8, )
INSTANCES_NUM(Word16, )
INSTANCES_NUM(Word32, )
INSTANCES_NUM(Word64, )

INSTANCES_FRACTIONAL(Float, )
INSTANCES_FRACTIONAL(Double, )
INSTANCES_FRACTIONAL(Complex a, RealFloat a =>)
INSTANCES_FRACTIONAL(Ratio a, Integral a =>)

max_stars_repo_name: mikebenfield/ftensor

max_issues_repo_name: mikebenfield/ftensor

max_forks_repo_name: mikebenfield/ftensor


text:
{-# Language DeriveGeneric     #-}
{-# Language OverloadedStrings #-}
{-# Language DeriveFunctor     #-}
{-# Language DeriveFoldable    #-}
{-# Language RecordWildCards   #-}

module Lib where

import           Control.Arrow                  ( (&&&) )
import           Data.Csv
import qualified Data.ByteString.Lazy          as BL
import qualified Data.Vector                   as V
import           Data.Either                    ( fromRight )
import           GHC.Generics
import           Statistics.Sample
import           System.IO.Unsafe
import           Data.Foldable                  ( toList
                                                , foldl'
                                                )
import           Data.List                      ( transpose )

-- The definition of a spambase entry.
-- It is SpamEntry a as that allows me to be able to map over the data easily.
type SpamEntry = SpamE Double -- Alias for nicer type signatures
data SpamE a =
  S {
    _make        ::  !a
  , _address     ::  !a
  , _all         ::  !a
  , _3d          ::  !a
  , _our         ::  !a
  , _over        ::  !a
  , _remove      ::  !a
  , _internet    ::  !a
  , _order       ::  !a
  , _mail        ::  !a
  , _receive     ::  !a
  , _will        ::  !a
  , _people      ::  !a
  , _report      ::  !a
  , _addresses   ::  !a
  , _free        ::  !a
  , _business    ::  !a
  , _email       ::  !a
  , _you         ::  !a
  , _credit      ::  !a
  , _your        ::  !a
  , _font        ::  !a
  , _000         ::  !a
  , _money       ::  !a
  , _hp          ::  !a
  , _hpl         ::  !a
  , _george      ::  !a
  , _650         ::  !a
  , _lab         ::  !a
  , _labs        ::  !a
  , _telnet      ::  !a
  , _857         ::  !a
  , _data        ::  !a
  , _415         ::  !a
  , _85          ::  !a
  , _technology  ::  !a
  , _1999        ::  !a
  , _parts       ::  !a
  , _pm          ::  !a
  , _direct      ::  !a
  , _cs          ::  !a
  , _meeting     ::  !a
  , _original    ::  !a
  , _project     ::  !a
  , _re          ::  !a
  , _edu         ::  !a
  , _table       ::  !a
  , _conference  ::  !a
  , _semicolon   ::  !a
  , _lparen      ::  !a
  , _pipe        ::  !a
  , _bang        ::  !a
  , _dollarsign  ::  !a
  , _hash        ::  !a
  , _caps_avg    ::  !a
  , _caps_max    ::  !a
  , _caps_total  ::  !a
  , _spam        ::  !Int
} deriving (Generic, Show, Functor, Foldable)
  -- Gives me mapping, toList, folding, decoding, etc., for free.

-- CSV decoding instanes for free because of Generic
instance FromField a => FromRecord (SpamE a)
instance ToField   a => ToRecord   (SpamE a)

--
-- This is written in a semi unidiomatic haskell style where I don't really
-- have a main function that assembles everything together.  Instead, I use
-- unsafePerformIO to grab the data and treat it as a pure value so I can just
-- hardcode the assignment here. It allows for faster iteratino time.
--

-- Decodes the file into a vector of spambase entries
{-# NOINLINE spamEntries #-}
spamEntries :: V.Vector SpamEntry
spamEntries = fromRight undefined . decode NoHeader $ mydata
 where
  -- Load the data from a file. Bit of a hacky way to do it, but it makes the
  -- rest of the code easier to incrementally write.
  mydata :: BL.ByteString
  mydata = unsafePerformIO $ BL.readFile "spambase.data"
  {-# NOINLINE mydata #-} -- if you inline this function, bad shit happens

-- This works by indexing the vector of SpamEntry entries
-- Then the vector is split into two vectors, one with all the even indices,
-- one with all odd (this effectively splits the data sets equally)
--
-- Then the extra indices are stripped from both vectors
test, train :: V.Vector SpamEntry
(test, train) = splitData spamEntries
 where
  splitData
    :: V.Vector SpamEntry
    -> (V.Vector SpamEntry, V.Vector SpamEntry)
  splitData e = both (snd <$>) $ V.partition (even . fst) (V.indexed e)

--
-- Probabilities of the prior
--
probSpam :: V.Vector SpamEntry -> Double
probSpam s = V.length (fst $ splitSpam s) // V.length s

pSpamTrain, pSpamTest :: Double
(pSpamTrain, pSpamTest) = both probSpam (train, test)

-- Takes the spam entries and collects the (mean, standard deviation) of all 57
-- features
stats :: Foldable t => V.Vector (t Double) -> [(Double, Double)]
stats entries = (mean &&& stdDev) <$> collectedStats
 where
  --  This is just data shuffling to get all 57 features from the different
  --  spamEntries into a list of feature-vectors
  collectedStats :: [V.Vector Double]
  collectedStats = V.fromList . normalize <$$> transpose $ toList <$> V.toList entries
  -- Normalize is a safety precaution to prevent zeroes
  normalize       = fmap (\x -> if x < 0.0001 then 0.0001 else x)

-- Splits into a tuple of (spam, non spam)
splitSpam :: V.Vector SpamEntry -> (V.Vector SpamEntry, V.Vector SpamEntry)
splitSpam = V.partition ((== 1) . _spam)

-- Just some globals to use the functions just defined above
trainSpamStat, trainNonStat, testSpamStat, testNonStat :: [(Double, Double)]
(trainSpamStat, trainNonStat) = both stats $ splitSpam train
(testSpamStat, testNonStat) = both stats $ splitSpam test

--  This is the N(xi;μ,σ) function from the text
normal :: Floating a => (a, a) -> a -> a
normal (μ, σ) x = 1 / root * exp (-(x - μ) ^ 2 / (2 * σ ^ 2))
  where root = sqrt (2 * pi) * σ

--
-- The meat of the library. This classifies things as spam or nonSpam
--
data Class = Spam | NonSpam deriving (Show)

classify :: SpamEntry -> Class
classify e = if testSpamLogs > testNonLogs then Spam else NonSpam
 where
  testSpamLogs = log pSpamTrain + sum (logify testSpamStat)
  testNonLogs  = log (1 - pSpamTrain) + sum (logify testNonStat)

  logify :: [(Double, Double)] -> [Double]
  logify stat = log <$> zipWith ($) (normal <$> stat) (toList e)

classify' :: SpamEntry -> Accuracy
classify' e = case (guess, actual) of
  (Spam   , Spam   ) -> TruePos
  (NonSpam, NonSpam) -> TrueNeg
  (Spam   , NonSpam) -> FalsePos
  (NonSpam, Spam   ) -> FalseNeg
 where
  guess  = classify e
  actual = if _spam e == 1 then Spam else NonSpam

--
-- Accuracy and analysis section
--
data Accuracy = TruePos | TrueNeg | FalsePos | FalseNeg deriving (Show)

-- The higher-ordered data structure is, again, so I can map over the data
-- easily.  foldl' (+) accuracyData 0 will sum up all of the items in this data
-- structure, for example.
type AccuracyData = ACount Int
data ACount a =
  ACount {
    _truePos  :: !a
  , _trueNeg  :: !a
  , _falsePos :: !a
  , _falseNeg :: !a
} deriving (Show, Functor, Foldable)

-- Default starting point
accuracyData :: AccuracyData
accuracyData = ACount 0 0 0 0

-- Updater functions: succ a <=> a+1
truePos :: AccuracyData -> AccuracyData
truePos a@ACount {..} = a { _truePos = succ _truePos }

trueNeg :: AccuracyData -> AccuracyData
trueNeg a@ACount {..} = a { _trueNeg = succ _trueNeg }

falsePos :: AccuracyData -> AccuracyData
falsePos a@ACount {..} = a { _falsePos = succ _falsePos }

falseNeg :: AccuracyData -> AccuracyData
falseNeg a@ACount {..} = a { _falseNeg = succ _falseNeg }

-- Go over all of the classified data; categorize it into the confusion matrix,
-- then tally up all of the results.
tally :: AccuracyData
tally = foldl' accum accuracyData classes
 where
  classes = classify' <$> test
  accum a e = case e of
    TruePos  -> truePos a
    TrueNeg  -> trueNeg a
    FalsePos -> falsePos a
    FalseNeg -> falseNeg a

--
-- The statistics:
--
accuracy :: Double
accuracy = allPos // total
 where
  allPos = _truePos tally + _trueNeg tally
  total  = foldl' (+) 0 tally

recall :: Double
recall = tp // (tp + _falseNeg tally) where tp = _truePos tally

precision :: Double
precision = tp // allPos
 where
  tp     = _truePos tally
  allPos = tp + _falsePos tally

--
-- Main
--
main = do
  putStrLn $ "The accuracy is: "         ++ show accuracy
  putStrLn $ "The recall is: "           ++ show recall
  putStrLn $ "The precision is: "        ++ show precision
  putStrLn $ "The confusion matrix is: " ++ show tally

--
-- Helper functions
--

-- Apply a function to both sides of a tuple at once
both :: (t -> b) -> (t, t) -> (b, b)
both f (a, b) = (f a, f b)

-- Nested mapping
infixl 5 <$$>
(<$$>) :: (Functor f2, Functor f1) => (a -> b) -> f1 (f2 a) -> f1 (f2 b)
(<$$>) = fmap fmap fmap

-- Division with implicit number conversion
(//) :: (Integral a1, Integral a, Fractional a2) => a -> a1 -> a2
a // b = fromIntegral a / fromIntegral b

max_stars_repo_name: jared-w/Haskell-Spambase

max_issues_repo_name: jared-w/Haskell-Spambase

max_forks_repo_name: jared-w/Haskell-Spambase


text:
{-|
Module: MachineLearning.Optimization.GradientDescent
Description: Gradient Descent
Copyright: (c) Alexander Ignatyev, 2016
License: BSD-3
Stability: experimental
Portability: POSIX

-}

module MachineLearning.Optimization.GradientDescent
(
  gradientDescent
)

where

import MachineLearning.Types (R, Vector, Matrix)
import MachineLearning.Regularization (Regularization)
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA

import qualified MachineLearning.Model as Model

-- | Gradient Descent method implementation. See "MachineLearning.Regression" for usage details.
gradientDescent :: Model.Model a => R-> a -> R -> Int -> Regularization -> Matrix -> Vector -> Vector -> (Vector, Matrix)
gradientDescent alpha model eps maxIters lambda x y theta = helper theta maxIters []
  where gradient = Model.gradient model lambda
        cost = Model.cost model lambda
        helper theta nIters optPath =
          let theta' = theta - (LA.scale alpha (gradient x y theta))
              j = cost x y theta'
              gradientTest = LA.norm_2 (theta' - theta) < eps
              optPathRow = V.concat [LA.vector [(fromIntegral $ maxIters - nIters), j], theta']
              optPath' = optPathRow : optPath
          in if gradientTest || nIters <= 1
             then (theta, LA.fromRows $ reverse optPath')
             else helper theta' (nIters - 1) optPath'

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Lib
    ( Mandelbrot(..)
    , makeMandelbrot
    , drawMandelbrot
    ) where

import Graphics.Image hiding (realPart, imagPart, scale)
import Data.Complex

maxDepth :: Int
maxDepth = 200

-- The function to test on
f :: Complex Double -> Complex Double -> Complex Double
f z c = z ** 2 + c

-- Return infinite list of all numbers in sequence
nextValues :: Complex Double -> Complex Double -> [Complex Double]
nextValues z c = z:nextValues (f z c) c

series :: Complex Double -> [Complex Double]
series pos = nextValues (0 :+ 0) pos

-- See if number is within radius of convergence (using the pythagorean theorem)
mightConverge :: Complex Double -> Bool
mightConverge z = (realPart z) ** 2 + (imagPart z) ** 2 <= 2 ** 2

-- Count how quickly the infinite list becomes guaranteed to never converge
speedOfBlowup :: Complex Double -> Int
speedOfBlowup pos = length $ takeWhile mightConverge $ take maxDepth $ series pos

-- A color value between 0 and 1, based on how quick the blowup is
valueOf :: Complex Double -> Double
valueOf pos = (fromIntegral $ speedOfBlowup pos) / (fromIntegral maxDepth)

-- Mandelbrot settings
data Mandelbrot = Mandelbrot
  { width :: Int
  , height :: Int
  , scale :: Double
  , x_offset :: Double
  , y_offset :: Double
  }

makeMandelbrot :: Int -> Int -> Mandelbrot
makeMandelbrot w h = Mandelbrot
  { width = w
  , height = h
  , scale = 1
  , x_offset = 0
  , y_offset = 0
  }

-- Pixel coordinates to complex numer
coordinate :: Mandelbrot -> (Int, Int) -> Complex Double
coordinate m (x, y) =
  let
    x_proportional = fromIntegral x / fromIntegral (width m)
    y_proportional = fromIntegral y / fromIntegral (height m)

    -- Instead of going from 0 to 1, go from -2 to 2
    x_normalised = (x_proportional - 0.5) * (2 / 0.5)
    y_normalised = (y_proportional - 0.5) * (2 / 0.5)
  in
    ((x_normalised - x_offset m) / scale m) :+ ((y_normalised - y_offset m) / scale m)

-- Make a mandelbrot image of the given size and zoom
drawMandelbrot :: Mandelbrot -> Image VS Y Double
drawMandelbrot m = makeImage (width m, height m) $ \(y, x) -> PixelY $ valueOf $ coordinate m (x, y)


text:
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE DataKinds #-}
{-# OPTIONS_GHC -Wno-simplifiable-class-constraints #-}
module Data.Matrix.Static.LinearAlgebra.Types
    ( Numeric(..)
    , Matrix
    , MMatrix
    , SparseMatrix
    ) where

import Data.Vector.Storable (Vector, Storable)
import Data.Vector.Storable.Mutable (MVector)
import Data.Complex (Complex)
import Foreign.C.Types

import qualified Data.Matrix.Static.Dense as D
import qualified Data.Matrix.Static.Dense.Mutable as DM
import qualified Data.Matrix.Static.Sparse as S

class (S.Zero a, Storable a, Num a) => Numeric a where
    foreignType :: a -> CInt

instance Numeric Float where foreignType _ = 0
instance Numeric CFloat where foreignType _ = 0
instance Numeric Double where foreignType _ =1
instance Numeric (Complex Float) where foreignType _ = 2
instance Numeric (Complex Double) where foreignType _ = 3

type Matrix r c a = D.Matrix r c Vector a
type MMatrix r c s a = DM.MMatrix r c MVector s a

type SparseMatrix r c a = S.SparseMatrix r c Vector a

max_stars_repo_name: kaizhang/matrix-sized

max_issues_repo_name: kaizhang/matrix-sized

max_forks_repo_name: kaizhang/matrix-sized


text:
#if __GLASGOW_HASKELL__ >= 703
{-# LANGUAGE Safe #-}
#endif
-----------------------------------------------------------------------------
-- |
-- Module      :  Control.DeepSeq
-- Copyright   :  (c) The University of Glasgow 2001-2009
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  libraries@haskell.org
-- Stability   :  stable
-- Portability :  portable
--
-- This module provides an overloaded function, 'deepseq', for fully
-- evaluating data structures (that is, evaluating to \"Normal Form\").
--
-- A typical use is to prevent resource leaks in lazy IO programs, by
-- forcing all characters from a file to be read. For example:
--
-- > import System.IO
-- > import Control.DeepSeq
-- >
-- > main = do
-- >     h <- openFile "f" ReadMode
-- >     s <- hGetContents h
-- >     s `deepseq` hClose h
-- >     return s
--
-- 'deepseq' differs from 'seq' as it traverses data structures deeply,
-- for example, 'seq' will evaluate only to the first constructor in
-- the list:
--
-- > > [1,2,undefined] `seq` 3
-- > 3
--
-- While 'deepseq' will force evaluation of all the list elements:
--
-- > > [1,2,undefined] `deepseq` 3
-- > *** Exception: Prelude.undefined
--
-- Another common use is to ensure any exceptions hidden within lazy
-- fields of a data structure do not leak outside the scope of the
-- exception handler, or to force evaluation of a data structure in one
-- thread, before passing to another thread (preventing work moving to
-- the wrong threads).
--
module Control.DeepSeq (
     deepseq, ($!!), force,
     NFData(..),
  ) where

import Data.Int
import Data.Word
import Data.Ratio
import Data.Complex
import Data.Array
import Data.Fixed
import Data.Version

infixr 0 $!!

-- | 'deepseq': fully evaluates the first argument, before returning the
-- second.
--
-- The name 'deepseq' is used to illustrate the relationship to 'seq':
-- where 'seq' is shallow in the sense that it only evaluates the top
-- level of its argument, 'deepseq' traverses the entire data structure
-- evaluating it completely.
--
-- 'deepseq' can be useful for forcing pending exceptions,
-- eradicating space leaks, or forcing lazy I/O to happen.  It is
-- also useful in conjunction with parallel Strategies (see the
-- @parallel@ package).
--
-- There is no guarantee about the ordering of evaluation.  The
-- implementation may evaluate the components of the structure in
-- any order or in parallel.  To impose an actual order on
-- evaluation, use 'pseq' from "Control.Parallel" in the
-- @parallel@ package.
--
deepseq :: NFData a => a -> b -> b
deepseq a b = rnf a `seq` b

-- | the deep analogue of '$!'.  In the expression @f $!! x@, @x@ is
-- fully evaluated before the function @f@ is applied to it.
($!!) :: (NFData a) => (a -> b) -> a -> b
f $!! x = x `deepseq` f x

-- | a variant of 'deepseq' that is useful in some circumstances:
--
-- > force x = x `deepseq` x
--
-- @force x@ fully evaluates @x@, and then returns it.  Note that
-- @force x@ only performs evaluation when the value of @force x@
-- itself is demanded, so essentially it turns shallow evaluation into
-- deep evaluation.

force :: (NFData a) => a -> a
force x = x `deepseq` x

-- | A class of types that can be fully evaluated.
class NFData a where
    -- | rnf should reduce its argument to normal form (that is, fully
    -- evaluate all sub-components), and then return '()'.
    --
    -- The default implementation of 'rnf' is
    --
    -- > rnf a = a `seq` ()
    --
    -- which may be convenient when defining instances for data types with
    -- no unevaluated fields (e.g. enumerations).
    rnf :: a -> ()
    rnf a = a `seq` ()

instance NFData Int
instance NFData Word
instance NFData Integer
instance NFData Float
instance NFData Double

instance NFData Char
instance NFData Bool
instance NFData ()

instance NFData Int8
instance NFData Int16
instance NFData Int32
instance NFData Int64

instance NFData Word8
instance NFData Word16
instance NFData Word32
instance NFData Word64

instance NFData (Fixed a)

-- |This instance is for convenience and consistency with 'seq'.
-- This assumes that WHNF is equivalent to NF for functions.
instance NFData (a -> b)

--Rational and complex numbers.

instance (Integral a, NFData a) => NFData (Ratio a) where
  rnf x = rnf (numerator x, denominator x)

instance (RealFloat a, NFData a) => NFData (Complex a) where
  rnf (x:+y) = rnf x `seq`
               rnf y `seq`
               ()

instance NFData a => NFData (Maybe a) where
    rnf Nothing  = ()
    rnf (Just x) = rnf x

instance (NFData a, NFData b) => NFData (Either a b) where
    rnf (Left x)  = rnf x
    rnf (Right y) = rnf y

instance NFData Data.Version.Version where
    rnf (Data.Version.Version branch tags) = rnf branch `seq` rnf tags

instance NFData a => NFData [a] where
    rnf [] = ()
    rnf (x:xs) = rnf x `seq` rnf xs

instance (Ix a, NFData a, NFData b) => NFData (Array a b) where
    rnf x = rnf (bounds x, Data.Array.elems x)

instance (NFData a, NFData b) => NFData (a,b) where
  rnf (x,y) = rnf x `seq` rnf y

instance (NFData a, NFData b, NFData c) => NFData (a,b,c) where
  rnf (x,y,z) = rnf x `seq` rnf y `seq` rnf z

instance (NFData a, NFData b, NFData c, NFData d) => NFData (a,b,c,d) where
  rnf (x1,x2,x3,x4) = rnf x1 `seq`
		        rnf x2 `seq`
		        rnf x3 `seq`
		        rnf x4

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5) =>
         NFData (a1, a2, a3, a4, a5) where
  rnf (x1, x2, x3, x4, x5) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6) =>
         NFData (a1, a2, a3, a4, a5, a6) where
  rnf (x1, x2, x3, x4, x5, x6) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7) =>
         NFData (a1, a2, a3, a4, a5, a6, a7) where
  rnf (x1, x2, x3, x4, x5, x6, x7) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7, NFData a8) =>
         NFData (a1, a2, a3, a4, a5, a6, a7, a8) where
  rnf (x1, x2, x3, x4, x5, x6, x7, x8) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7 `seq`
                  rnf x8

instance (NFData a1, NFData a2, NFData a3, NFData a4, NFData a5, NFData a6, NFData a7, NFData a8, NFData a9) =>
         NFData (a1, a2, a3, a4, a5, a6, a7, a8, a9) where
  rnf (x1, x2, x3, x4, x5, x6, x7, x8, x9) =
                  rnf x1 `seq`
                  rnf x2 `seq`
                  rnf x3 `seq`
                  rnf x4 `seq`
                  rnf x5 `seq`
                  rnf x6 `seq`
                  rnf x7 `seq`
                  rnf x8 `seq`
                  rnf x9

max_stars_repo_name: emorins/nahc

max_issues_repo_name: emorins/nahc

max_forks_repo_name: emorins/nahc


text:
module Hasql.Private.Prelude
(
  module Exports,
  LazyByteString,
  ByteStringBuilder,
  LazyText,
  TextBuilder,
  forMToZero_,
  forMFromZero_,
  strictCons,
  mapLeft,
)
where

-- base
-------------------------
import Control.Applicative as Exports hiding (WrappedArrow(..))
import Control.Arrow as Exports hiding (first, second)
import Control.Category as Exports
import Control.Concurrent as Exports
import Control.Exception as Exports
import Control.Monad as Exports hiding (fail, mapM_, sequence_, forM_, msum, mapM, sequence, forM)
import Control.Monad.IO.Class as Exports
import Control.Monad.Fail as Exports
import Control.Monad.Fix as Exports hiding (fix)
import Control.Monad.ST as Exports
import Data.Bifunctor as Exports
import Data.Bits as Exports
import Data.Bool as Exports
import Data.Char as Exports
import Data.Coerce as Exports
import Data.Complex as Exports
import Data.Data as Exports
import Data.Dynamic as Exports
import Data.Either as Exports
import Data.Fixed as Exports
import Data.Foldable as Exports hiding (toList)
import Data.Function as Exports hiding (id, (.))
import Data.Functor as Exports
import Data.Functor.Compose as Exports
import Data.Functor.Contravariant as Exports
import Data.Int as Exports
import Data.IORef as Exports
import Data.Ix as Exports
import Data.List as Exports hiding (sortOn, isSubsequenceOf, uncons, concat, foldr, foldl1, maximum, minimum, product, sum, all, and, any, concatMap, elem, foldl, foldr1, notElem, or, find, maximumBy, minimumBy, mapAccumL, mapAccumR, foldl')
import Data.List.NonEmpty as Exports (NonEmpty(..))
import Data.Maybe as Exports
import Data.Monoid as Exports hiding (Alt, (<>))
import Data.Ord as Exports
import Data.Proxy as Exports
import Data.Ratio as Exports
import Data.Semigroup as Exports hiding (First(..), Last(..))
import Data.STRef as Exports
import Data.String as Exports
import Data.Traversable as Exports
import Data.Tuple as Exports
import Data.Unique as Exports
import Data.Version as Exports
import Data.Void as Exports
import Data.Word as Exports
import Debug.Trace as Exports
import Foreign.ForeignPtr as Exports
import Foreign.Ptr as Exports
import Foreign.StablePtr as Exports
import Foreign.Storable as Exports
import GHC.Conc as Exports hiding (orElse, withMVar, threadWaitWriteSTM, threadWaitWrite, threadWaitReadSTM, threadWaitRead)
import GHC.Exts as Exports (IsList(..), lazy, inline, sortWith, groupWith)
import GHC.Generics as Exports (Generic)
import GHC.IO.Exception as Exports
import GHC.OverloadedLabels as Exports
import Numeric as Exports
import Prelude as Exports hiding (Read, fail, concat, foldr, mapM_, sequence_, foldl1, maximum, minimum, product, sum, all, and, any, concatMap, elem, foldl, foldr1, notElem, or, mapM, sequence, id, (.))
import System.Environment as Exports
import System.Exit as Exports
import System.IO as Exports (Handle, hClose)
import System.IO.Error as Exports
import System.IO.Unsafe as Exports
import System.Mem as Exports
import System.Mem.StableName as Exports
import System.Timeout as Exports
import Text.ParserCombinators.ReadP as Exports (ReadP, readP_to_S, readS_to_P)
import Text.ParserCombinators.ReadPrec as Exports (ReadPrec, readPrec_to_P, readP_to_Prec, readPrec_to_S, readS_to_Prec)
import Text.Printf as Exports (printf, hPrintf)
import Unsafe.Coerce as Exports

-- transformers
-------------------------
import Control.Monad.IO.Class as Exports
import Control.Monad.Trans.Class as Exports
import Control.Monad.Trans.Cont as Exports hiding (shift, callCC)
import Control.Monad.Trans.Except as Exports (ExceptT(ExceptT), Except, except, runExcept, runExceptT, mapExcept, mapExceptT, withExcept, withExceptT, throwE, catchE)
import Control.Monad.Trans.Maybe as Exports
import Control.Monad.Trans.Reader as Exports (Reader, ask, runReader, mapReader, withReader, ReaderT(ReaderT), runReaderT, mapReaderT, withReaderT)
import Control.Monad.Trans.State.Strict as Exports (State, runState, evalState, execState, mapState, withState, StateT(StateT), runStateT, evalStateT, execStateT, mapStateT, withStateT)
import Control.Monad.Trans.Writer.Strict as Exports (Writer, runWriter, execWriter, mapWriter, WriterT(..), execWriterT, mapWriterT)
import Data.Functor.Compose as Exports
import Data.Functor.Identity as Exports

-- mtl
-------------------------
import Control.Monad.Error.Class as Exports (MonadError (..))

-- profunctors
-------------------------
import Data.Profunctor.Unsafe as Exports

-- contravariant
-------------------------
import Data.Functor.Contravariant.Divisible as Exports

-- hashable
-------------------------
import Data.Hashable as Exports (Hashable(..))

-- text
-------------------------
import Data.Text as Exports (Text)

-- bytestring
-------------------------
import Data.ByteString as Exports (ByteString)

-- vector
-------------------------
import Data.Vector as Exports (Vector)

-- dlist
-------------------------
import Data.DList as Exports (DList)

-- postgresql-binary
-------------------------
import PostgreSQL.Binary.Data as Exports (UUID)

-- custom
-------------------------
import qualified Data.Text.Lazy
import qualified Data.Text.Lazy.Builder
import qualified Data.ByteString.Lazy
import qualified Data.ByteString.Builder

type LazyByteString =
  Data.ByteString.Lazy.ByteString

type ByteStringBuilder =
  Data.ByteString.Builder.Builder

type LazyText =
  Data.Text.Lazy.Text

type TextBuilder =
  Data.Text.Lazy.Builder.Builder

{-# INLINE forMToZero_ #-}
forMToZero_ :: Applicative m => Int -> (Int -> m a) -> m ()
forMToZero_ !startN f =
  ($ pred startN) $ fix $ \loop !n -> if n >= 0 then f n *> loop (pred n) else pure ()

{-# INLINE forMFromZero_ #-}
forMFromZero_ :: Applicative m => Int -> (Int -> m a) -> m ()
forMFromZero_ !endN f =
  ($ 0) $ fix $ \loop !n -> if n < endN then f n *> loop (succ n) else pure ()

{-# INLINE strictCons #-}
strictCons :: a -> [a] -> [a]
strictCons !a b =
  let !c = a : b in c

{-# INLINE mapLeft #-}
mapLeft :: (a -> c) -> Either a b -> Either c b
mapLeft f =
  either (Left . f) Right

max_stars_repo_name: guaraqe/hasql

max_issues_repo_name: guaraqe/hasql

max_forks_repo_name: guaraqe/hasql


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}

module MLP.Config (
   readConf,
   mkDataSet
) where

import Data.Aeson (eitherDecode)
import Numeric.LinearAlgebra.Data (loadMatrix, (><), Matrix)
import Numeric.LinearAlgebra.Static (create, L)
import Control.Lens ((^.), foldlOf, (%~), (&))
import MLP.DataSet (toPat, normPat)
import MLP.Network (AllCon, Net, Learn, MLP)
import MLP.Types (trainFile,
                  testFile,
                  topology,
                  normalise,
                  topoFold,
                  DataSet(..),
                  Parameters(..),
                  State(..),
                  Topology,
                  MonadFileSystem(..))
import Control.Monad ((=<<))
import System.FilePath (splitFileName, takeFileName, (</>))
import GHC.TypeLits (KnownNat, Nat)
import Control.Monad.Random.Class (MonadRandom)

liftEither :: (MonadFail m) => Either String b -> m b
liftEither = either fail pure

readConf :: MonadFileSystem m => FilePath -> m Parameters
readConf fl = do
   s <- liftEither . eitherDecode =<< readFileM fl
   pure $ s & trainFile %~ absPath
            & testFile %~ absPath

   where
      (base,_) = splitFileName fl
      absPath x = base </> takeFileName x

mkDataSet :: (KnownNat i, KnownNat o, MonadFileSystem m) => Parameters -> m (DataSet i o)
mkDataSet p = do
      trainSet <- mkPat (p ^. trainFile)
      testSet <- mkPat (p ^. testFile)

      return $
         if norm then
            DataSet (normPat trainSet) (normPat testSet)
         else
            DataSet trainSet testSet
   where
      mkPat fl = readMatrixM fl >>= liftEither . toPat
      norm = p ^. normalise

max_stars_repo_name: ovanr/Typed-MLP

max_issues_repo_name: ovanr/Typed-MLP

max_forks_repo_name: ovanr/Typed-MLP


text:
{- |
Module      :  Numeric.GSL.Fitting
Copyright   :  (c) Alberto Ruiz 2010
License     :  GPL

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses ffi

Nonlinear Least-Squares Fitting

<http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html>

The example program in the GSL manual (see examples/fitting.hs):

@dat = [
 ([0.0],([6.0133918608118675],0.1)),
 ([1.0],([5.5153769909966535],0.1)),
 ([2.0],([5.261094606015287],0.1)),
 ...
 ([39.0],([1.0619821710802808],0.1))]

expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]

expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]

(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]

\> path
(6><5)
 [ 1.0,  76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
 , 2.0, 37.683816318260355,  2.858760367632973,  8.092094813253975e-2, 1.4479636296208662
 , 3.0,    9.5807893736187,  4.948995119561291,   0.11942927999921617, 1.0945766509238248
 , 4.0,  5.630494933603935,  5.021755718065913,   0.10287787128056883, 1.0338835440862608
 , 5.0,  5.443976278682909,  5.045204331329302,   0.10405523433131504,  1.019416067207375
 , 6.0, 5.4439736648994685,  5.045357818922331,   0.10404905846029407, 1.0192487112786812 ]
\> sol
[(5.045357818922331,6.027976702418132e-2),
(0.10404905846029407,3.157045047172834e-3),
(1.0192487112786812,3.782067731353722e-2)]@

-}
-----------------------------------------------------------------------------

module Numeric.GSL.Fitting (
    -- * Levenberg-Marquardt
    nlFitting, FittingMethod(..),
    -- * Utilities
    fitModelScaled, fitModel
) where

import Data.Packed.Internal
import Numeric.LinearAlgebra
import Numeric.GSL.Internal

import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
import Foreign.C.Types
import System.IO.Unsafe(unsafePerformIO)

-------------------------------------------------------------------------

data FittingMethod = LevenbergMarquardtScaled -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom.
                   | LevenbergMarquardt -- ^ This is an unscaled version of the lmder algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of lmder converges too slowly, or the function is already scaled appropriately.
        deriving (Enum,Eq,Show,Bounded)

-- | Nonlinear multidimensional least-squares fitting.
nlFitting :: FittingMethod
      -> Double                     -- ^ absolute tolerance
      -> Double                     -- ^ relative tolerance
      -> Int                        -- ^ maximum number of iterations allowed
      -> (Vector Double -> Vector Double)     -- ^ function to be minimized
      -> (Vector Double -> Matrix Double)   -- ^ Jacobian
      -> Vector Double                   -- ^ starting point
      -> (Vector Double, Matrix Double)  -- ^ solution vector and optimization path

nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit

nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
    let p   = dim xiv
        n   = dim (f xiv)
    fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
    jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
    rawpath <- createMatrix RowMajor maxit (2+p)
    app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
    let it = round (rawpath @@> (maxit-1,0))
        path = takeRows it rawpath
        [sol] = toRows $ dropRows (it-1) path
    freeHaskellFunPtr fp
    freeHaskellFunPtr jp
    return (subVector 2 p sol, path)

foreign import ccall safe "nlfit"
    c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM

-------------------------------------------------------

checkdim1 n _p v
    | dim v == n = v
    | otherwise = error $ "Error: "++ show n
                        ++ " components expected in the result of the function supplied to nlFitting"

checkdim2 n p m
    | rows m == n && cols m == p = m
    | otherwise = error $ "Error: "++ show n ++ "x" ++ show p
                        ++ " Jacobian expected in nlFitting"

------------------------------------------------------------

err (model,deriv) dat vsol = zip sol errs where
    sol = toList vsol
    c = max 1 (chi/sqrt (fromIntegral dof))
    dof = length dat - (rows cov)
    chi = norm2 (fromList $ cost (resMs model) dat sol)
    js = fromLists $ jacobian (resDs deriv) dat sol
    cov = inv $ trans js <> js
    errs = toList $ scalar c * sqrt (takeDiag cov)

-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and
-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
-- instances (x, (y,sigma)) to be fitted.

fitModelScaled
         :: Double -- ^ absolute tolerance
         -> Double -- ^ relative tolerance
         -> Int    -- ^ maximum number of iterations allowed
         -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
         -> [(x, ([Double], Double))] -- ^ instances
         -> [Double] -- ^ starting point
         -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where
    (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit
                (fromList . cost (resMs model) dt . toList)
                (fromLists . jacobian (resDs deriv) dt . toList)
                (fromList xin)

-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
-- instances (x,y) to be fitted.

fitModel :: Double -- ^ absolute tolerance
         -> Double -- ^ relative tolerance
         -> Int    -- ^ maximum number of iterations allowed
         -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
         -> [(x, [Double])]  -- ^ instances
         -> [Double] -- ^ starting point
         -> ([Double], Matrix Double) -- ^ solution and optimization path
fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where
    (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
                (fromList . cost (resM model) dt . toList)
                (fromLists . jacobian (resD deriv) dt . toList)
                (fromList xin)

cost model ds vs = concatMap (model vs) ds

jacobian modelDer ds vs = concatMap (modelDer vs) ds

-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s

-- | Associated derivative for 'resMs'.
resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x)

-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent
-- to 'resMs' with all sigmas = 1.
resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b

-- | Associated derivative for 'resM'.
resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
resD m v = \(x,_) -> m v x

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE DeriveGeneric, DeriveAnyClass, BangPatterns #-}

module Report where

import           Control.DeepSeq
import           Control.Monad
import           Criterion.IO
import           Criterion.Types
import           Data.List
import qualified Data.Map                      as M
import           Data.Maybe
import           GHC.Generics
import           Statistics.Types
import           System.Directory
import           System.FilePath
import           Text.Printf
import           Data.Bifunctor
import           Data.Char

-- | Map test names to test measures
type Measures = M.Map String Measure

data Measure = Measure
  { mTest :: String -- ^Test name, with the format <test-kind>/<test-obj>-<test-pkg>
  , mValue :: Double -- ^Execution time (in ms or bytes, depending on the test)
  } deriving (Show, Read, Eq, Generic, NFData, Ord)

{- |
>>> import Data.List
>>> let m = Measure "deserialization (time)/BinTree Direction-binary" 1989
>>> let m2 = Measure "serialization (time)/BinTree Direction-binary" 1004
>>> let m3 = Measure "size (bytes)/BinTree Direction-binary" 6291455
>>> let m4 = Measure "size (bytes)/Cars-binary" 301455
>>> let m5 = Measure "size (bytes)/Cars-flat" 300000
>>> let ms = measures [m,m2,m3,m4,m5]

>>> mKindPkg m
("deserialization (time)/BinTree Direction","binary")

>>> mKindPkgVal m
("deserialization (time)/BinTree Direction",("binary",1989.0))

>>> mKind m
"deserialization (time)/BinTree Direction"

>>> mObj m
"BinTree Direction"

>>> mObjSub m
"BinTree Direction-binary"

>>> mSub m
"binary"

>>> allOf mObj ms
["BinTree Direction","Cars"]

>>> allOf mType ms
["deserialization (time)","serialization (time)","size (bytes)"]

>>> summaryTable ms
[("BinTree Direction",[("deserialization (time)",[(1989.0,[Measure {mTest = "deserialization (time)/BinTree Direction-binary", mValue = 1989.0}])]),("serialization (time)",[(1004.0,[Measure {mTest = "serialization (time)/BinTree Direction-binary", mValue = 1004.0}])]),("size (bytes)",[(6291455.0,[Measure {mTest = "size (bytes)/BinTree Direction-binary", mValue = 6291455.0}])])]),("Cars",[("size (bytes)",[(300000.0,[Measure {mTest = "size (bytes)/Cars-flat", mValue = 300000.0}]),(301455.0,[Measure {mTest = "size (bytes)/Cars-binary", mValue = 301455.0}])])])]

>>> renderTable ms
"|Dataset\\Measure|deserialization|serialization|size|\n| ---| ---| ---| ---|\n|BinTree Direction|[binary](https://hackage.haskell.org/package/binary)|[binary](https://hackage.haskell.org/package/binary)|[binary](https://hackage.haskell.org/package/binary)|\n|Cars|||[flat](https://hackage.haskell.org/package/flat),[binary](https://hackage.haskell.org/package/binary)|\n"

>>> addTransfers_ ms
fromList [("deserialization (time)/BinTree Direction-binary",Measure {mTest = "deserialization (time)/BinTree Direction-binary", mValue = 1989.0}),("serialization (time)/BinTree Direction-binary",Measure {mTest = "serialization (time)/BinTree Direction-binary", mValue = 1004.0}),("size (bytes)/BinTree Direction-binary",Measure {mTest = "size (bytes)/BinTree Direction-binary", mValue = 6291455.0}),("size (bytes)/Cars-binary",Measure {mTest = "size (bytes)/Cars-binary", mValue = 301455.0}),("size (bytes)/Cars-flat",Measure {mTest = "size (bytes)/Cars-flat", mValue = 300000.0}),("transfer [10 MBits] (time)/BinTree Direction-binary",Measure {mTest = "transfer [10 MBits] (time)/BinTree Direction-binary", mValue = 8026.164}),("transfer [100 MBits] (time)/BinTree Direction-binary",Measure {mTest = "transfer [100 MBits] (time)/BinTree Direction-binary", mValue = 3496.3164}),("transfer [1000 MBits] (time)/BinTree Direction-binary",Measure {mTest = "transfer [1000 MBits] (time)/BinTree Direction-binary", mValue = 3043.33164})]
-}
renderTable :: Measures -> String
renderTable ms =
  let
    tests = allOf mObj ms
    kinds = allOf mType ms
    vals  = allOf mSub ms
    lines s = s : map (\t -> showPkgs $ map (\v -> tos t s v ms) vals) kinds
  in
    unlines
    . map mdRow
    $ ("Dataset\\Measure" : map short kinds)
    : replicate (length kinds + 1) " ---"
    : map lines tests
 where
  pkgVals = map snd . tops . sort . catMaybes . map
    ((\m -> (mValue m, pkgRef $ mSub m)) <$>)
  tops [] = []
  tops hs =
    let limit = fst (head hs) * 1.3 in takeWhile (\e -> fst e <= limit) hs
  showPkgs = intercalate "," . pkgVals
  mdRow vs = concat ["|", intercalate "|" vs, "|"]
  pkgRef name =
    concat ["[", name, "](https://hackage.haskell.org/package/", name, ")"]
  short = unwords . init . words

tos :: String -> String -> String -> Measures -> Maybe Measure
tos t o s ms = M.lookup (concat [t, "/", o, "-", s]) ms

summaryTable :: M.Map k Measure -> [(String, [(String, [(Double, [Measure])])])]
summaryTable =
  map (second $ (map (second $ by mValue) . by mType)) . by mObj . M.elems

addTransfers :: FilePath -> IO ()
addTransfers workDir = addMeasures__ workDir addTransfers_ >> return ()

addTransfers_ :: Measures -> Measures
addTransfers_ ms =
  let
    tests = allOf mObjSub ms
    addT :: Measures -> String -> Measures
    addT ms o =
      fromMaybe ms
        $   (\ser des siz ->
              let
                add1 :: Double -> Measures -> Measures
                add1 megaBits ms =
                  let
                    trans =
                      concat
                        [ "transfer ["
                        , show (round megaBits :: Integer)
                        , " MBits] (time)/"
                        , o
                        ]
                    time =
                      mValue ser + mValue des + mValue siz * 8 / (megaBits * 1000)
                  in
                    add (Measure trans time) ms
              in
                add1 1000 (add1 100 (add1 10 ms))
            )
        <$> M.lookup ("serialization (time)/" ++ o) ms
        <*> M.lookup ("deserialization (time)/" ++ o) ms
        <*> M.lookup ("size (bytes)/" ++ o) ms
  in
    foldl addT ms tests

add :: Measure -> Measures -> Measures
add m = M.insert (mTest m) m

mKindPkg :: Measure -> (String, String)
mKindPkg m = let (k, _ : p) = break (== '-') $ mTest m in (k, p)

mKind :: Measure -> String
mKind = fst . mKindPkg

mType :: Measure -> String
mType = brk1 '/' . mKind

mObj :: Measure -> String
mObj = brk2 '/' . mKind

mObjSub :: Measure -> String
mObjSub = brk2 '/' . mTest

mSub :: Measure -> String
mSub = brk2 '-' . mObjSub

brk1 :: Eq a => a -> [a] -> [a]
brk1 sep = fst . break (== sep)

brk2 :: Eq a => a -> [a] -> [a]
brk2 sep = tail . snd . break (== sep)

mKindPkgVal :: Measure -> (String, (String, Double))
mKindPkgVal m = let (k, p) = mKindPkg m in (k, (p, mValue m))

byTestKind :: Measures -> [(String, (String, Double))]
byTestKind = map mKindPkgVal . M.elems

allKinds :: Measures -> [String]
allKinds = allOf mKind

allOf :: Eq a1 => (a2 -> a1) -> M.Map k a2 -> [a1]
allOf f = nub . map f . M.elems

allTests :: [(String, (String, Double))] -> [(String, [(String, Double)])]
allTests = -- by fst
           sort . map (\g -> (fst . head $ g, map snd g)) . groupBy
  (\a b -> fst a == fst b)

-- order by dimension
by :: (Ord a, Ord t) => (t -> a) -> [t] -> [(a, [t])]
by f = sort . map (\g -> (f $ head g, g)) . groupBy (\a b -> f a == f b)

toMeasures :: [Report] -> Measures
toMeasures = measures . map toMeasure

measures :: [Measure] -> Measures
measures = M.fromList . map (\m -> (mTest m, m))

toMeasure :: Report -> Measure
toMeasure r =
  Measure (reportName r) ((1000 *) . estPoint . anMean . reportAnalysis $ r)

deleteMeasures :: FilePath -> IO ()
deleteMeasures = removeFile . measuresFile

updateMeasures :: FilePath -> IO ()
updateMeasures = void . updateMeasures_

updateMeasures_ :: FilePath -> IO (Measures, Measures, Measures)
updateMeasures_ dir = do
  m' <- readCriterionMeasures dir
  addMeasures_ dir m'

addMeasures_ :: FilePath -> Measures -> IO (Measures, Measures, Measures)
addMeasures_ dir m' = addMeasures__ dir (const m')

addMeasures__
  :: FilePath -> (Measures -> Measures) -> IO (Measures, Measures, Measures)
addMeasures__ dir f = do
  !m <- readMeasures dir
  let m'  = f m
  let m'' = M.union m' m
  writeMeasures dir m''
  return (m, m', m'')

addMeasures :: FilePath -> String -> [(String, Double)] -> IO ()
addMeasures dir name ms = void $ addMeasures_
  dir
  (M.fromList $ map
    (\(pkg, val) -> let n = concat [name, "-", pkg] in (n, Measure n val))
    ms
  )

readMeasures :: FilePath -> IO Measures
readMeasures dir = do
  let f = measuresFile dir
  fe <- doesFileExist f
  if not fe then return M.empty else (force . read) <$> readFile f

writeMeasures :: FilePath -> Measures -> IO ()
writeMeasures dir = writeFile (measuresFile dir) . show

readCriterionMeasures :: FilePath -> IO Measures
readCriterionMeasures dir = toMeasures <$> readReports (reportsFile dir)

measuresFile :: FilePath -> FilePath
measuresFile dir = dir </> "measures"

reportsFile :: FilePath -> FilePath
reportsFile dir = dir </> "report.json"

reportsMDFile :: FilePath -> FilePath
reportsMDFile dir = dir </> "report.md"

printMeasures :: FilePath -> IO ()
printMeasures dir = reportMeasures_ dir >>= putStrLn

printSummary :: FilePath -> (String -> Bool) -> IO ()
printSummary dir f =
  readMeasures dir >>= putStrLn . renderTable . M.filter (f . mTest)

reportMeasures :: FilePath -> IO ()
reportMeasures dir =
  reportMeasures_ dir
    >>= writeFile (reportsMDFile dir)
    . ("Results that are within 30% of the best result are displayed in **bold**.\n\n" ++
      )

reportMeasures_ :: FilePath -> IO String
reportMeasures_ dir = reportMeasures__ (const True) <$> readMeasures dir

printMeasuresCurrent :: (Measures, Measures, Measures) -> IO ()
printMeasuresCurrent (_, m', m'') =
  let currentKinds = allKinds m'
  in  printMeasuresAll_ (\ts -> fst ts `elem` currentKinds) m''

printMeasuresAll_ :: ((String, [(String, Double)]) -> Bool) -> Measures -> IO ()
printMeasuresAll_ f = putStrLn . reportMeasures__ f

reportMeasures__ :: ((String, [(String, Double)]) -> Bool) -> Measures -> String
reportMeasures__ f =
  concatMap (uncurry report) . filter f . allTests . byTestKind

printMeasuresDiff :: (Measures, Measures, Measures) -> IO ()
printMeasuresDiff (m, m', _) =
  mapM_ (\(n, d) -> putStrLn (unwords [n, show d, "%"]))
    . M.toList
    $ M.intersectionWith
        (\a b -> round ((mValue b / mValue a - 1) * 100) :: Int)
        m
        m'

readReports :: FilePath -> IO [Report]
readReports jsonReportFile = do
  fe <- doesFileExist jsonReportFile
  if not fe
    then return [] -- M.empty
    else do
      er <- readJSONReports jsonReportFile
      return $ case er of
        Right (_, _, reports) -> reports
        _                     -> [] -- M.empty

report :: String -> [(String, Double)] -> String
report _ [] = []
report name rs =
  let (_, rss) = report_ rs
      width    = maximum . map (length . fst) $ rs
      out =
        [ "| ---| ---|"
        , "| package | performance |"
        , ""
        , unwords ["####", name, "(best first)"]
        ]
          -- out = ["| ---| ---| ---|","| package | measure | relative measure |","",unwords ["####",name, "(best first)"]]
  in  unlines
      . reverse
      $ ""
      : foldl'
          (\out (n, r, a) ->
            let bold n = concat ["**", n, "**"]
                marked n = -- fix prob with github bold display
                           if r <= 1.3
                  then if isSpace (head n)
                    then let (n1, n2) = span isSpace n in n1 ++ bold n2
                    else
                      let (n1, n2) = span (not . isSpace) n in bold n1 ++ n2
                  else n
            in  unwords
                    [ "|"
                    , marked (printString width n) -- pkg name
                    , "|"
                  -- , printf "%11.1f" a -- absolute measure
                  -- ,"|"
                    , marked (printDouble r) -- relative measure
                    , "|"
                    ]
                  : out
          )
          out
          rss

report_ :: (Fractional c, Ord c) => [(a, c)] -> ((a, c), [(a, c, c)])
report_ rs =
  let rss  = sortOn snd rs
      best = snd . head $ rss
  in  (head rss, map (\(n, v) -> (n, v / best, v)) rss)

printDouble :: Double -> String
printDouble = printf "%7.1f"

printInt :: Double -> String
printInt = printf "%.0f"

printString :: Int -> String -> String
printString width = printf ("%-" ++ show width ++ "s")

max_stars_repo_name: effectfully/flat

max_issues_repo_name: tittoassini/flat

max_forks_repo_name: tittoassini/flat


text:
{-# LANGUAGE DataKinds, GADTs, LambdaCase, NoImplicitPrelude, RankNTypes #-}
{-# LANGUAGE TypeApplications, TypeOperators                             #-}
{-# OPTIONS_GHC -fno-warn-unused-imports #-}
{-# OPTIONS_GHC -Wno-type-defaults #-}
{-# OPTIONS_GHC -fplugin Data.Singletons.TypeNats.Presburger #-}
module Algebra.Algorithms.Groebner.ZeroDimSpec where
import           Algebra.Algorithms.Groebner
import           Algebra.Algorithms.ZeroDim
import           Algebra.Internal
import           Algebra.Internal
import           Algebra.Prelude.Core             hiding ((%))
import           Algebra.Ring.Ideal
import           Algebra.Ring.Polynomial
import           Algebra.Ring.Polynomial.Quotient
import qualified Prelude                          as P

import           Control.Monad
import           Control.Monad.Random
import           Data.Complex
import           Data.Convertible       (convert)
import qualified Data.Foldable          as F
import qualified Data.Matrix            as M
import           Data.Maybe
import qualified Data.Sized.Builtin     as SV
import qualified Data.Vector            as V
import           Numeric.Field.Fraction (Fraction, (%))
import           Numeric.Natural
import           Test.Hspec
import           Test.Hspec.QuickCheck
import           Test.HUnit
import           Test.QuickCheck        hiding (promote)
import           Utils

default (Fraction Integer, Natural)

asGenListOf :: Gen [a] -> a -> Gen [a]
asGenListOf = const

spec :: Spec
spec = parallel $ do
  describe "solveLinear" $ do
    it "solves data set correctly" $
      forM_ linSet $ \set ->
      solveLinear (M.fromLists $ inputMat set) (V.fromList $ inputVec set)
      `shouldBe` Just (V.fromList $ answer set)
    prop "solves any solvable cases" $
      forAll (resize 10 arbitrary) $ \(MatrixCase ms) ->
      let mat = M.fromLists ms :: M.Matrix (Fraction Integer)
      in rank mat == M.ncols mat ==>
           forAll (vector (length $ head ms)) $ \v ->
             let ans = M.getCol 1 $ mat P.* M.colVector (V.fromList v)
             in solveLinear mat ans == Just (V.fromList v)
    it "cannot solve unsolvable cases" $
      pendingWith "need example"
  describe "univPoly" $ modifyMaxSuccess (const 25) $ modifyMaxSize (const 4) $
    prop "produces monic generators of the elimination ideal" $
      checkForTypeNat [2..4] prop_univPoly
  describe "radical" $
    it "really computes radical" $
      pendingWith "We can verify correctness by comparing with singular, but it's not quite smart way..."

  describe "fglm" $ modifyMaxSuccess (const 25) $ modifyMaxSize (const 3) $ do
    prop "computes monomial basis" $
      checkForTypeNat [2..4] $ \sdim ->
        case zeroOrSucc sdim of
          IsZero -> error "impossible"
          IsSucc k ->
            withWitness (lneqSucc k) $ withWitness (lneqZero k) $
            withKnownNat k $
            forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
            let base =
                  reifyQuotient (mapIdeal (changeOrder Lex) ideal) $ \ii ->
                  map quotRepr $ fromJust $ standardMonomials' ii
            in stdReduced (snd $ fglm ideal) == stdReduced base
    prop "computes lex base" $
      checkForTypeNat [2..4] $ \sdim ->
        case zeroOrSucc sdim of
          IsZero -> error "impossible"
          IsSucc k -> withKnownNat k $
            withWitness (lneqSucc k) $
            withWitness (lneqZero k) $
            forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
            stdReduced (fst $ fglm ideal)
              == stdReduced (calcGroebnerBasisWith Lex ideal)
    prop "returns lex base in descending order" $
      checkForTypeNat [2..4] $ \sdim ->
      case zeroOrSucc sdim of
        IsZero -> error "impossible"
        IsSucc k -> withKnownNat k $
          case (lneqSucc k, lneqZero k) of
            (Witness, Witness) ->
              forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
              isDescending (map leadingMonomial $ fst $ fglm ideal)
  describe "solve'" $ modifyMaxSuccess (const 50) $ modifyMaxSize (const 4) $ do
    it "solves equation with admissible error" $
      checkForTypeNat [2..4] $ prop_isApproximateZero 1e-5 (pure . solve' 1e-5)
    describe "solves regression cases correctly" $
      forM_ companionRegressions $
      approxZeroTestCase 1e-5 (pure . solve' 1e-5)

  describe "solveViaCompanion" $
    modifyMaxSuccess (const 50) $ modifyMaxSize (const 4) $ do
    it "solves equation with admissible error" $
      checkForTypeNat [2..4] $ prop_isApproximateZero 1e-5
        (pure . solveViaCompanion 1e-5)
    describe "solves regression cases correctly" $
      forM_ companionRegressions $
      approxZeroTestCase 1e-5 (pure . solveViaCompanion 1e-5)

  describe "solveM" $ modifyMaxSuccess (const 50) $ modifyMaxSize (const 4) $ do
    prop "solves equation with admissible error" $
      checkForTypeNat [2..4] $ prop_isApproximateZero 1e-9 solveM
    describe "solves regressions correctly" $
      forM_ companionRegressions $
      approxZeroTestCase 1e-9 solveM

isDescending :: Ord a => [a] -> Bool
isDescending xs = and $ zipWith (>=) xs (drop 1 xs)

prop_isApproximateZero
  :: (KnownNat n)
  => Double
  -> (forall m. ((0 < m) ~ 'True, KnownNat m) =>
      Ideal (Polynomial (Fraction Integer) m) -> IO [Sized m (Complex Double)])
  -> SNat n -> Property
prop_isApproximateZero err solver sn =
  case zeroOrSucc sn of
    IsSucc _ -> withKnownNat sn $
      forAll (zeroDimOf sn) $ ioProperty . checkSolverApproxZero err solver
    IsZero -> error "prop_isApproximateZero must be called with non-zero typenats!"

companionRegressions
  :: [SolverTestCase]
companionRegressions =
  [ SolverTestCase @3 $
    let [x,y,z] = vars
    in ZeroDimIdeal $ toIdeal
        [ 2*x^2 + (0.5 :: Fraction Integer) .*. x * y
        , -(0.5 :: Fraction Integer) .*. y ^ 2 + y * z
        , z^2 - z]
  , SolverTestCase @3 $
    let [x,y,z] = vars
    in ZeroDimIdeal $ toIdeal
        [ 2*x + 1, 2*y^2,-z^2 - 1]
  , SolverTestCase @2 $
    let [x,y] = vars
    in ZeroDimIdeal $ toIdeal
        [ x^2 - y, -2*y^2]
  ]

approxZeroTestCase
  :: Double
  -> (forall n. ((0 < n) ~ 'True, KnownNat n)
        => Ideal (Polynomial (Fraction Integer) n) -> IO [Sized n (Complex Double)]
      )
  -> SolverTestCase
  -> SpecWith ()
approxZeroTestCase err calc (SolverTestCase f@(ZeroDimIdeal i)) =
  it (show $ getIdeal f) $ do
    isZeroDimensional (F.toList i) @? "Not zero-dimensional!"
    checkSolverApproxZero err calc f
      @? "Solved correctly"

checkSolverApproxZero
  :: ((0 < n) ~ 'True, KnownNat n)
  => Double
  -> (Ideal (Polynomial (Fraction Integer) n) -> IO [Sized n (Complex Double)])
  -> ZeroDimIdeal n
  -> IO Bool
checkSolverApproxZero err solver (ZeroDimIdeal ideal) = do
  anss <- solver ideal
  let mul r d = convert r * d
  pure
    $ all (\as -> all ((<err) . magnitude . substWith mul as)
    $ generators ideal) anss

data SolverTestCase where
  SolverTestCase
    :: ((0 < n) ~ 'True, KnownNat n)
    => ZeroDimIdeal n -> SolverTestCase

prop_univPoly :: KnownNat n => SNat n -> Property
prop_univPoly sdim =
  forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
  let ods = enumOrdinal sdim
  in conjoin $ flip map ods $ \nth ->
  let gen = univPoly nth ideal
  in forAll (unaryPoly sdim nth) $ \f ->
  (f `modPolynomial` [gen] == 0) == (f `isIdealMember` ideal)

rank :: (Ord r, Fractional r) => M.Matrix r -> Int
rank mat =
  let Just (u, _, _, _,_, _) = M.luDecomp' mat
  in V.foldr (\a acc -> if a /= 0 then acc + 1 else acc) (0 :: Int) $ M.getDiag u

data TestSet = TestSet { inputMat :: [[Fraction Integer]]
                       , inputVec :: [Fraction Integer]
                       , answer   :: [Fraction Integer]
                       } deriving (Show, Eq, Ord)

linSet :: [TestSet]
linSet =
  [TestSet
   [[1 ,0 ,0 ,0 ,0 ]
   ,[0 ,(-2) ,(-2) ,(-2) ,(-2) ]
   ,[0 ,0 ,3 % 2,0 ,(-1) % 2]
   ,[0 ,0 ,0 ,0 ,(-5) % 2]
   ,[0 ,1 ,1 ,1 ,1 ]
   ,[0 ,0 ,(-2) ,1 ,(-1) ]
   ]
   [0 ,(-2) ,19 % 5,14 % 5,1 ,0 ]
   [0 ,(-81) % 25,54 % 25,16 % 5,(-28) % 25]
  ]

max_stars_repo_name: MxmUrw/computational-algebra

max_issues_repo_name: MxmUrw/computational-algebra

max_forks_repo_name: MxmUrw/computational-algebra


text:
{-# LANGUAGE DefaultSignatures      #-}
{-# LANGUAGE DeriveFoldable         #-}
{-# LANGUAGE DeriveFunctor          #-}
{-# LANGUAGE DeriveTraversable      #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE RankNTypes             #-}
{-# LANGUAGE RecordWildCards        #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TemplateHaskell        #-}
{-# LANGUAGE TypeFamilies           #-}
{-# LANGUAGE TypeOperators          #-}
module Circuits.Analysis.DC where

import           Control.Lens
import           Control.Monad.Primitive
import           Control.Monad.ST
import           Data.Primitive.MutVar
import qualified Data.Vector                  as V
import qualified Data.Vector.Storable         as VS
import qualified Data.Vector.Storable.Mutable as VSM
import           Generics.Deriving.Lens
import           GHC.Generics
import qualified Numeric.LinearAlgebra        as HMatrix

import           Debug.Trace

import           Circuits.Circuit

-- | Describes the time scale of a DC simulation. We can either perform a DC steady-state or a transient analysis.
data TimeScale
  = SteadyState
  -- ^ solving for steady state
  | Step
    { _currentTime :: Double -- ^ the current simulation time in seconds since t0
    , _timeStep    :: Double -- ^ the size of current time step to be taken
    }
    -- ^ solving for next time step

makeLenses ''TimeScale

-- | Type class providing an interface for the component stamps used in DC analysis.
class SimulateDC c where
  -- | Beginning to solve the next time step.
  beginSolve :: Simulation m Double sim => c -> TimeScale -> sim -> m ()
  default beginSolve :: (Generic c, GSimulateDC (Rep c), Simulation m Double sim) => c -> TimeScale -> sim -> m ()
  beginSolve c t s = forOf_ generic c $ \gc -> gbeginSolve gc t s
  -- | Beginning the next iteration of solving the non-linear parts.
  -- This should set up the right hand side based on the current element state.
  beginIteration :: Simulation m Double sim => c -> TimeScale -> sim -> m ()
  default beginIteration :: (Generic c, GSimulateDC (Rep c), Simulation m Double sim) => c -> TimeScale -> sim -> m ()
  beginIteration c t s = forOf_ generic c $ \gc -> gbeginIteration gc t s
  -- | End of the current iteration, after the referenced variables have been updated.
  -- It should update the state of non-linear parts, but not the state of linear parts.
  endIteration :: c -> TimeScale -> c
  default endIteration :: (Generic c, GSimulateDC (Rep c)) => c -> TimeScale -> c
  endIteration c t = over generic (`gendIteration` t) c
  -- | Finished solving the current time step. Here, the state of linear parts should be updated.
  endSolve :: c -> TimeScale -> c
  default endSolve :: (Generic c, GSimulateDC (Rep c)) => c -> TimeScale -> c
  endSolve c t = over generic (`gendSolve` t) c

instance SimulateDC Int where
  beginSolve _ _ _ = pure ()
  beginIteration _ _ _ = pure ()
  endIteration x _ = x
  endSolve x _ = x

instance SimulateDC Double where
  beginSolve _ _ _ = pure ()
  beginIteration _ _ _ = pure ()
  endIteration x _ = x
  endSolve x _ = x

instance SimulateDC Variable
instance SimulateDC Node
instance SimulateDC Branch

instance SimulateDC c => SimulateDC (Circuit c)

-- * Generic deriving of SimulateDC instances

class GSimulateDC f where
  gbeginSolve :: Simulation m Double sim => f p -> TimeScale -> sim -> m ()
  gbeginIteration :: Simulation m Double sim => f p -> TimeScale -> sim -> m ()
  gendIteration :: f p -> TimeScale -> f p
  gendSolve :: f p -> TimeScale -> f p

instance GSimulateDC V1 where
  gbeginSolve _ _ _ = return ()
  gbeginIteration _ _ _ = return ()
  gendIteration x _ = x
  gendSolve x _ = x

instance GSimulateDC U1 where
  gbeginSolve _ _ _ = return ()
  gbeginIteration _ _ _ = return ()
  gendIteration x _ = x
  gendSolve x _ = x

instance SimulateDC a => GSimulateDC (K1 i a) where
  gbeginSolve (K1 x) = beginSolve x
  gbeginIteration (K1 x) = beginIteration x
  gendIteration (K1 x) = K1 . endIteration x
  gendSolve (K1 x) = K1 . endSolve x

instance GSimulateDC f => GSimulateDC (M1 i t f) where
  gbeginSolve (M1 x) = gbeginSolve x
  gbeginIteration (M1 x) = gbeginIteration x
  gendIteration (M1 x) = M1 . gendIteration x
  gendSolve (M1 x) = M1 . gendSolve x

instance (GSimulateDC f, GSimulateDC g) => GSimulateDC (f :+: g) where
  gbeginSolve (L1 x) = gbeginSolve x
  gbeginSolve (R1 x) = gbeginSolve x
  gbeginIteration (L1 x) = gbeginIteration x
  gbeginIteration (R1 x) = gbeginIteration x
  gendIteration (L1 x) = L1 . gendIteration x
  gendIteration (R1 x) = R1 . gendIteration x
  gendSolve (L1 x) = L1 . gendSolve x
  gendSolve (R1 x) = R1 . gendSolve x

instance (GSimulateDC f, GSimulateDC g) => GSimulateDC (f :*: g) where
  gbeginSolve (x :*: y) t s = gbeginSolve x t s >> gbeginSolve y t s
  gbeginIteration (x :*: y) t s = gbeginIteration x t s >> gbeginIteration y t s
  gendIteration (x :*: y) t = gendIteration x t :*: gendIteration y t
  gendSolve (x :*: y) t = gendSolve x t :*: gendSolve y t

-- * DC Simulation

-- | A mutable system of linear equations suitable for representing linear circuits.
data SimulationState s v = SimulationState
  { size                 :: Int -- ^ the number of variables in this system
  , coefficients         :: VSM.MVector s v -- ^ the underlying storage of the equations coefficients in row-major order
  , coefficientsModified :: MutVar s Bool -- ^ flag indicating whether the coefficient matrix has been modified
  , rhs                  :: VSM.MVector s v -- ^ the right hand side of the equations
  , rhsModified          :: MutVar s Bool -- ^ flag indicating whether the RHS has been modified
  }

-- | Creates a new mutable linear circuit equation system with the given number of variables.
newSimulation :: (PrimMonad m, VSM.Storable v, Num v) => Int -> m (SimulationState (PrimState m) v)
newSimulation n =
  SimulationState n <$> VSM.replicate (n * n) 0 <*> newMutVar False <*> VSM.replicate n 0 <*> newMutVar False

cleanCoefficients :: (PrimMonad m) => SimulationState (PrimState m) v -> m ()
cleanCoefficients s = writeMutVar (coefficientsModified s) False

cleanRhs :: (PrimMonad m) => SimulationState (PrimState m) v -> m ()
cleanRhs s = writeMutVar (rhsModified s) False

-- | Computes the index into the flat buffer of matrix variables from a row and column.
flatIndex :: SimulationState s v -- ^ the equation system
          -> Int -- ^ matrix row
          -> Int -- ^ matrix column
          -> Int -- ^ flat buffer index
flatIndex sim row col = row * size sim + col

class Monad m => Simulation m v sys | sys -> v  where
  stampMatrix :: sys -> Variable -> Variable -> v -> m ()
  stampRhs :: sys -> Variable -> v -> m ()

instance (PrimMonad m, PrimState m ~ s, Num v, VS.Storable v) => Simulation m v (SimulationState s v) where
  stampMatrix sys (Index row _) (Index col _) value = do
    VSM.modify (coefficients sys) (+ value) (flatIndex sys row col)
    writeMutVar (coefficientsModified sys) True
  stampMatrix _ _ _ _ = return () -- no stamping when ground is involved

  stampRhs sys (Index row _) value = do
    VSM.modify (rhs sys) (+ value) row
    writeMutVar (rhsModified sys) True
  stampRhs _ Ground _ = return () -- no stamping for ground

initialSolution :: HasVariables c => Circuit c -> VS.Vector Double
initialSolution c = runST $ do
  vec <- VSM.new (c ^. variableCount)
  forOf_ (model . variables . _Index) c $ uncurry (VSM.write vec)
  VS.unsafeFreeze vec

step :: forall c. (HasVariables c, SimulateDC c) => TimeScale -> Circuit c -> Circuit c
step time c = runST simulate where
  simulate :: ST s (Circuit c)
  simulate = do
    sim <- newSimulation (c ^. variableCount)
    -- initialize matrix
    beginSolve c time sim
    mat <- HMatrix.reshape (c ^. variableCount) <$> VS.unsafeFreeze (coefficients sim)
    cleanCoefficients sim
    -- iterate
    writeMutVar (rhsModified sim) True
    c' <- solveNonLinear sim (HMatrix.luPacked mat) (initialSolution c) c
    -- finalize
    return $ endSolve c' time

  solveNonLinear :: SimulationState s Double -> HMatrix.LU Double -> VS.Vector Double -> Circuit c -> ST s (Circuit c)
  solveNonLinear sim lu prevSolution prevCircuit = do
    -- prepare RHS
    beginIteration c time sim
    rhsUpdated <- readMutVar (rhsModified sim)
    cleanRhs sim
    if rhsUpdated
      then do
        -- solve and update components
        rhs' <- VS.unsafeFreeze $ rhs sim
        let solution = HMatrix.flatten $ HMatrix.luSolve lu (HMatrix.asColumn rhs')
            c' = endIteration (updateValues solution c) time
            maxChange = VS.maximum $ VS.zipWith (\x y -> abs (x - y)) solution prevSolution
        if maxChange < 1e-6
          then return c'
          else solveNonLinear sim lu solution c'
      else return prevCircuit

-- | Update all variable references inside a circuit description with its new values.
updateValues :: HasVariables a => VS.Vector Double -> a -> a
updateValues values = over variables replace where
  replace Ground      = Ground
  replace (Index i _) = vars V.! i
  vars = V.imap Index $ VS.convert values

max_stars_repo_name: fatho/circuits

max_issues_repo_name: fatho/circuits

max_forks_repo_name: fatho/circuits


text:
-- |
-- Module      :  Test.FFTW
-- Copyright   :  (c) 2016 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Test.FFTW (
    fft
  ) where

import Data.Complex
import qualified Data.Vector.Storable as V
import Numeric.FFT.Vector.Unnormalized as FFTW

fft :: Int -> V.Vector (Complex Double) -> V.Vector (Complex Double)
fft n = FFTW.execute $ FFTW.plan FFTW.dft n

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral