text:
module Stats where

-- import Fraction
import Misc
import Maths
import Constants
import Distribution_Tables

import List
import Matrix
import Control.Monad
import System.Random
import Tuple
import Text.Printf

import qualified Data.ByteString.Lazy as BL
import qualified Data.List
import Control.Monad.Trans.State
import System.Process
import System.IO
import Statistics.Distribution.FDistribution
import Numeric.SpecFunctions (incompleteBeta,invIncompleteBeta )
import Data.List(foldl')
------------------------- VALUES

path_data_stats = "/home/mika/Documents/Cours/Stats/data/"

------------------------- DATA TYPES

-- readDF s =

data DataFrame α = DataFrame [Col α]
  -- deriving Show

-- data DFCell = C String

-- instance Show DFCell where
--   show (C s) = s

-- instance Eq DFCell where
--   (C s) == (C s') = s==s'

-- instance Ord DFCell where
--   compare (C s) (C s') = compare s s'

-- instance Show a => Show (DataFrame a) where
--   show (DataFrame cols) = unlines . map (concatWith "|") . transpose . map prepareColumn $ cols
--     where prepareColumn (Col h c) = all_same_size . (h:) . map show $ c
--             where all_same_size l = map (fix'size k) l
--                     where k = reduce max $ map count l

--   show'dataframe :: Matrix (Col String) -> String
-- show'dataframe (Matrix r c cells) =
--   showTable . transpose
--   . map (\(Col t c) -> [t]++c)
--   $ cells

------------------------- DATA CONSTRUCTORS

-- readTable sep = map (splitWhen sep) . lines

-- data DF = DF {
--   getStrings :: DFStr,
--   getNums    :: DFNum
--   }

-- data DFSth α = DFSth {getNames :: [String], getCols :: [C α]}
-- type DFStr = DFSth String
-- type DFNum = DFSth Double

-- data C α = C [Maybe α]

-- colToStrings :: C -> [String]
-- colToStrings (CStr values) = map jstr values
--   where jstr Nothing = "NA"
--         jstr (Just s) = s
-- colToStrings (CNb values) = map jnb values
--   where
--         jnb Nothing = "NA"
--         jnb (Just n) = show n

-- -- not gonna happen
-- -- instance Functor C where
-- --   fmap f (CStr c) = CStr $ f <$> c
-- --   fmap f (CNb c) = CNb $ f <$> c

-- instance Show DF where
-- --  show (DF cols) = ""
--   show (DF names cols) = unlines . map (concatWith "|")
--                          . transpose $ t
--     where scol (title, values) = title : horizbar
--             : colToStrings values
--           horizbar = take 99 (repeat '-')
--           namesbar = map (\t -> [t] ++ [horizbar]) names
--           t = map fixCol $ zipWith (++) namesbar $ map colToStrings cols
-- --          t = map (fixCol . scol) cols -- [[String]]
--           fixCol c = map (fixedLength ' ' k) c
--             where k = foldr max 0 . map length . dropwall $ c
--                   dropwall (x:y:zs) = x:zs

-- -- instance Functor DF where
-- --   fmap f (DF df) = DF . f $ df

-- -- list to fixed length
-- -- fill with e's if too short,
-- -- cut the end if too long
-- fixedLength e n l = (take k $ repeat e) ++ take k' l
--   where k = n - length l
--         k' = min (length l) n

-- -- column-first table from a <Char>-separated string
-- parsecsv :: Char -> String -> [[String]]
-- parsecsv csv_sep str = transpose $ p [] [] [] str False
--   where p table [] [] [] q = table
--         p table line [] [] q = p (table++[line]) [] [] [] q
--         p table line cell [] q = p table (line++[cell]) [] [] q
--         p table line cell s@(x:xs) False
--           | x=='\n'    = p (table ++ [line ++ [cell]]) [] [] xs False
--           | x==csv_sep = p table (line++[cell]) [] xs False
--         p table line cell ('"':'"':xs) quoted =
--           p table line (cell++"\"") xs quoted
--         p table line cell ('"':xs) quoted =
--           p table line cell xs (not quoted)
--         p table line cell (x:xs) quoted =
--           p table line (cell++[x]) xs quoted

-- csv_to_df :: [[String]] -> DF
-- csv_to_df t = DF (flatten names) cols
--   where (names, cols_str) = unzip . map (sliceN 1) $ t
--         cols = map p cols_str
--         p values
--           | all isNumber . concat $  values = CNb
--             $ map (\c -> if length c>0
--                          then Just (readDouble c)
--                          else Nothing) values
--           | otherwise = CStr . map Just $ values

-- -- R like subsetting with boolean vectors
-- -- subdf :: [Bool] -> [Bool] -> DF -> DF
-- -- subdf r c (DF df) = DF . fc . fr $ df

-- fc b (DF names cols) = DF (f names) (f cols)
--   where f = filterBool b
-- -- fc bs df = foldr (\(b,v) -> \vals -> if b then v:vals else vals)
-- --         [] $ zip bs df

-- -- fr b (DF names cols) = DF names $ map (filterBool b) cols
-- --   where fboolcol =

-- filter values by a boolean vector
filterBool b l = foldl f [] $ zip b l
  where f l (b,e) = if b then l++[e] else l

-- fr bs df =

-- table :: C -> T
-- table (CStr c) = 1
-- table (CNb c) = 1

-- weird function (similar to "by")
-- create DataFrame from a table [row]
-- by specifying how to get the factor
dataTableFromTable kfactor kdata l =
  cols
  where factorl = map (!!kfactor) l
        datal   = map (!!kdata) l
        tuples = zip factorl datal
        groups = groupBy comp tuples
          where comp a b = compare (fst a) (fst b)
        titles = map (fst . head) groups
        vals = map (map (readDouble . snd)) groups
        cols = zipWith Col titles vals

-- takes the first element as a header
-- assumes it is a list of columns, (not a list of rows)
tableToDF :: (String -> a) ->  [[String]] -> [Col a]
tableToDF f = map (\(h:c) -> Col h (map f c))

---------------- PROBABILITY

----- PRN

-- rollDie :: (Random a,Num a) => IO a
-- rollDie = randomRIO(1,6)

-- rollNdice n = mapM (\_ -> rollDie) [1..n]

-- rollNdice :: Int -> [Int]
-- rollNdice 0 = []
-- rollNdice n = rollDie

-- rollDie :: State StdGen Int
-- rollDie = state $ randomR (1, 6)

rollDie :: State StdGen Int
rollDie = do
  generator <- get
  let (value,newGen) = randomR (1,6) generator
  put newGen
  return value

-- frd :: (Monad m) => (State ((),s))
-- frd = do
--   a <- get
--   put $ a
--   return ()

-------

----------- Probability

-- choose :: Integral α => α -> α -> α
-- choose n k
--   | n < k     = error "Choose n k"
--   | otherwise = div a b
--   where α = [1..n]
--         β = [1..k] ++ [1..(n-k)]
--         a = foldr (*) 1 α
--         b = foldr (*) 1 β

choose :: Integral α => α -> α -> α
choose n k = reduce (*) α `div` reduce (*) β
  where α = [ max k (n-k) + 1 .. n ]
        β = [1 .. min k (n-k)]

showchooseInt :: Int -> Int -> [[Int]]
showchooseInt n k = showchooseList [0..n-1] k

showchooseList :: [α] -> Int -> [[α]]
showchooseList _ 0 = []
showchooseList l 1 = map (:[]) l
showchooseList l k = flatmap f [drop x l | x <- [0..n-k]]
  where n = count l
        f (x:xs) = map (x:) (showchooseList xs (k-1))

-- binom :: (Integral α) => α -> Fraction -> α -> Fraction
-- binom :: (Fractional b, Integral b) => b -> b -> b -> b
-- binom :: (Num a,Integral b) => b -> a -> b -> a
binom n p k = fromIntegral (choose n k) * p^k * (1-p)^(n-k)

-- dbinom :: (Integral b, Fractional b) => b -> b -> [b]
dbinom n p = map (binom n p) [0..n]

-- cloppearson
--   :: (Show t, Fractional t, Integral t) => t -> Fraction -> [(t, Fraction)]

-- cloppearson :: Fraction -> Fraction -> [(Fraction, Fraction)]
cloppearson n p = let
  dbin = dbinom n p
  props = map (/fromInteger n) [0..(fromInteger n)]
  risk = 0.025
  l = zip props dbin

  f (pacc, acc) (p,b) = if acc+b>risk then (pacc,acc) else (p, acc+b)

  min = foldr f (0,0) l
  max = foldr f (0,0) . reverse $ l

  in [min,max]

-- dbinom :: (Integral α,Real β) => α -> β -> [Fraction]
-- dbinom n p = la ++ middle ++ reverse la
--   where (a,b) = (\(a,b) -> (fromInteger a,b)) $ divMod (toInteger n) 2
--         q=1-p
--         la = map f [0..a-1]
--         middle = map f $ if b == 1 then [] else [a]
--         f = binom n p
        -- something m = (choose n m *) . toFraction $  p^m * q^(n-m)

-- graph'dbinom :: (Integral α,Real β) => α -> β -> IO ()
-- graph'dbinom a b = f $ dbinom a b
--   where f l = putStrLn . ("\n"++) . unlines . map g $ l
--           where m = reduce max l
--                 -- g :: Fraction -> String
--                 g e = ((++) (fix'size 6 ' ' $ printPrecision 2 (e * 100))) . take k . repeat $ '#'
--                   where k = fromInteger . round $ k'
--                         k' =  fromFraction $ 29 * e / m

-- pbinom :: (Integral α, Real β) => α -> β -> Fraction
-- pbinom

-- -- total = number of balls
-- -- goal  = number of "right" balls
-- data Lottery = Lottery {
--   total :: Integer,
--   goal  :: Integer}
--   deriving Show

-- -- compute the probability to get
-- -- k good balls
-- -- from a lottery
-- proba_lottery :: Lottery -> Integer -> Fraction
-- proba_lottery (Lottery n b) k =
--   choose b k * choose (n-b) (b-k) / choose n b

-- all_proba_lottery l = map (proba_lottery l) [0..g]
--   where g = goal l

-- pretty_all_proba_lottery l@(Lottery n b) =
--   putStrLn.unlines.map (printPrecision 15).all_proba_lottery $ l

random01 :: (RandomGen γ) => γ -> (Double, γ)
random01 = (\(a,b) -> (a :: Double, b)) . randomR (0,1)

-- for computation
sqrt2pi = sqrt $ 2 * pi
recipsqrt2pi = 1 / sqrt2pi

{- Normal Distribution -}

-- box-muller
-- normalr mean sd g = ([f cos, f sin],gb)
--   where (ra,ga) = random01 g
--         (rb,gb) = random01 ga
--         c = sqrt $ (-2) * log ra
--         f = (+mean) . (*sd) . (*c) . ($ 2 * pi * rb)

normalr :: (RandomGen g) => Double -> Double -> State g [Double]
normalr mean sd = do
  g <- get
  let (ra,g') = random01 g
      (rb,g'') = random01 g'
      c = sqrt $ (-2) * log ra
      f = (+mean) . (*sd) . (*c) . ($ 2 * pi * rb)
  put g''
  return [f cos, f sin]

-- rnorm mean sd k g
--   | k <= 0    = ([],g)
--   | otherwise = ((take n la) ++ lb ++ lc,g''')
--   where n = min k 2
--         (k',r) = (k-2) `divMod` 2
--         k'' = k' + r
--         ( la  ,g')   = normalr  mean sd     g
--         ( lb  ,g'')  = rnorm mean sd k'  g'
--         ( lc  ,g''') = rnorm mean sd k'' g''

rnorm :: (RandomGen g) => Double -> Double -> Int -> State g [Double]
rnorm mean sd k
  | k <= 0 = return []
  | otherwise = do
      let n = min k 2
          (k',r) = (k-2) `divMod` 2
          k'' = k' + r
      ns <- normalr mean sd
      la <- rnorm mean sd k'
      lb <- rnorm mean sd k''
      return $ take n ns ++ la ++ lb

-- density function
dnorm x mean sd =
  ( recipsqrt2pi / sd )
  * e ** (negate ((x-mean)^2) / (2*sd^2))

-- approximation of the error function
-- http://www.johndcook.com/blog/python_phi/
erf x = signum x * y
  where
    --  constants
    -- todo define out of function ? efficiency? no idea
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    x' = abs(x)

    -- A & S 7.1.26
    t = 1.0/(1.0 + p*x')
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t* e ** (-x'*x')

-- Cody, W. D. (1993) Algorithm 715: SPECFUN – A portable FORTRAN package of special function routines and test drivers. ACM Transactions on Mathematical Software 19, 22–32
-- integral from -Inf
-- cumulative distribution of the normal function (cdf)
pnorm = (/2) . (+1) . erf . (/sqrt2)

-- Wichura, M. J. (1988) Algorithm AS 241: The percentage points of the normal distribution. Applied Statistics, 37, 477–484
qnorm = undefined

-- beta :: Double
-- beta x y = sqrt2pi * f x * f y / f (x+y)
--   where f e = e ^ (e-0.5)

-----------------------------
{-     TABLE GENERATOR     -}
-----------------------------

-- just like https://homes.cs.washington.edu/~jrl/normal_cdf.pdf
table_cumulative_normal_proba :: String
table_cumulative_normal_proba =
  unlines . addRowSep
  . map (concat . addColSep . intercal "\t" . map tfl)
  $ header : body
  where centi = (/100) <$> [0..9]
        range =  (/10) <$> [-34..34]
        header = ("z":) $ printPrecision 2 . toFraction <$> centi
        row x =
          (printPrecision 1 . toFraction $ x) :
          (printPrecision 4 . toFraction <$> (( pnorm . (add x))
                                              <$> centi))
        body = row <$> range
        add a b = signum a * (abs a + b)
        addColSep (x:y:xs) = x : "┃" : xs
        addRowSep (x:xs) = x : sep : xs
          where l = (`replicate` '━')
                sep = l 7 ++ ('╋' : l 80)
        -- toFixedLength 0 _      = []
        -- toFixedLength k ""     = toFixedLength (k-1) "" ++ " "
        -- toFixedLength k (x:xs) = x : toFixedLength (k-1) xs
        toFixedLength k s = (replicate (k - count s) ' ') ++ take (max k (count s)) s
        tfl = toFixedLength 7

-- Generate F-Distribution tables
-- from a call to Rscript, and R's "qf" function
table_qf_from_R_lol :: Fraction -> IO F_TABLE
table_qf_from_R_lol α =
  ((F_TABLE nums denoms α)<$>)
  $ sequence
  $ (`map`cmds)
  $ callRscript
  . combineLine

  where
        αstr = printPrecision 4 $ 1-α
        nums   = [1..10] ++ [12,15,20,24,30,40,60,120,infinint]
        denoms = [1..30] ++ [40,60,120,infinint]
        -- cmds = gen_cmd <$> nums <*> denoms
        cmds = cmd_generator nums denoms
        combineLine l = "cat(" ++ concatWith [',','"',' ','"',','] l ++ ")"

        gen_cmd num denom = "qf(" ++ αstr
                            ++ ", df1=" ++ show num
                            ++ ", df2=" ++ show denom
                            ++")"

        cmd_generator nums denoms =
          (<$>denoms) $ (\el -> fmap ($ e) (gen_cmd <$> nums))

--callRscript :: String -> IO [Double]
        callRscript s =
          ((\(_,Just hout,_,_) -> hout)
           <$> createProcess (proc "Rscript" ["-e",s]){std_out=CreatePipe})
          >>= ((convert <$>) . hGetContents)
          where convert = map readDouble . words

f_distrib_alphas = [0.1,0.05,0.025,0.01,0.001] :: [Fraction]
io_f_distrib_tables = sequence $ table_qf_from_R_lol <$> f_distrib_alphas

-- -- returns the F critical value
-- -- for a given set of α , df num and df denom
-- -- Nothing if not in pre-computed tables
-- f_stat_critical :: Fraction -> Int -> Int -> Maybe Double
-- f_stat_critical α num denom = safe_head t >>= lookup num denom
--   where t = filter ((==α) . alpha) f_distribution_tables
--         lookup num denom (F_TABLE n d α v) = do
--           ni <- safe_indexOf num   n
--           di <- safe_indexOf denom d
--           return (v !! di !! ni)

f_stat_critical :: Fraction -> Int -> Int -> Maybe Double
f_stat_critical α' n' d' =
  let α = fromFraction α'
      n = fromIntegral n'
      d = fromIntegral d'
      y = n * α
  -- in Just $ incompleteBeta (0.5 * n) (0.5 * d) (y/(d+y))
      x = invIncompleteBeta (0.5 * n) (0.5 * d) α
  in Just $ d * x / (n * (1 - x))

-----------------------------
{- DESCRIPTIVE STATISTICS  -}
-----------------------------

-- mean :: (Fractional α,Foldable t) => t α -> α
-- -- mean l = (/) (count l) (sum l)
-- mean = uncurry (/) . foldr (\a (s,c) -> (s+a,c+1) ) (0,0)

-- sumSq :: (Fractional α,Foldable t,Functor t) => t α -> α
-- sumSq = sum . fmap (^2)

-- varPop :: (Fractional α,Functor t,Foldable t) => t α -> α
-- -- varPop l = (/µ) . sumSq . fmap (µ-) $ l
-- varPop l = sumSq l - (µ^2)
--   where µ = mean l

-- var :: (Fractional α,Functor t,Foldable t) => t α -> α
-- var l = (/(n-1)) . sumSq . fmap (µ-) $ l
--   where n = count l
--         µ = (/n) $ sum l

-- std :: (Real α,Fractional α, Foldable t, Functor t,Floating β) => t α -> β
-- std = sqrt . fromFraction . var

median :: (Ord α,Fractional α) => [α] -> α
median l = (if r == 0
            then head
            else (µ.take 2))
           . drop (q-1) . Data.List.sort $ l
  where (q,r) = (`divMod`2) . count $ l

-- mean after applying function
-- offset to number of values (for n-1)
mean :: (Fractional α) => (α -> α) -> α -> [α] -> α
-- mean f offset values = (/ (count values + offset)) . sum . map f $ values
mean f offset = uncurry (/) . foldr (\ v (s,c) -> (s+f v,c+1) ) (0,offset)

-- µ :: (Fractional α) => [α] -> α
-- µ = (\(a,b)->a/b) . foldl (\(s,k) e->(s+e,k+1)) (0,0)

µ :: (Fractional α) => [α] -> α
µ = mean id 0

-- computes a list's deviations from it's mean
µdeviations :: (Fractional α) => [α] -> [α]
-- µdeviations l = map (subtract mean) l
--   where mean = µ l
µdeviations = (uncurry ((.(subtract.µ)) . (flip map))) . duplicate

-- mean deviation statistic = sum |xi-µx|
mean_deviation δ = mean abs  δ . µdeviations
variance δ       = mean (^2) δ . µdeviations

squares :: (Fractional α) => [α] -> [α]
squares = map (^2)
-- sumsquares = sum . squares

ss :: (Num α,Fractional α) => [α] -> α
ss = sum . squares . µdeviations

σ² :: Fractional c => [c] -> c
σ² = variance 0

σ²' :: Fractional c => [(c,c)] -> c
σ²' l = (sum . map (uncurry (*).applyToFst (^2)) $ l) - ((^2) . sum . map (uncurry (*)) $ l)

σ = sqrt . σ²

s² = variance (-1)
s = sqrt . s²

md_pop    = mean_deviation 0
md_sample = mean_deviation (-1)

-- anniversary :: Integer -> Fraction
-- anniversary k = 1 - reduce (*) [365-k+1..365] % (365^k)

-----------------------------
{- INFERENTIAL STATISTICS  -}
-----------------------------

-- unicode utf-8 utf8 characters
-- ∪∩∀∈

cov :: (Fractional α) => α -> [α] -> [α] -> α
-- cov x y = µ $ deviationsProducts x y
cov δ = ((mean id δ).) . deviationsProducts

-- population's covariance
σxy :: Fractional α => [α] -> [α] -> α
σxy = cov 0

-- sample's covariance
sxy :: Fractional α => [α] -> [α] -> α
sxy = cov (-1)

-- ∀i∈{1,n} (xi - µx) * (yi - µy)
deviationsProducts x y = zipWith (*) (µdeviations x) (µdeviations y)

-- Pearson's correlation coefficient
r :: (Fractional α,Real α) => [α] -> [α] -> Double
r x y =  ssxy / sqrt ( ss x * ss y )
  where ssxy = fromFraction $ sum $ deviationsProducts x y
        ss = fromFraction . sum . squares . µdeviations

r² :: (Fractional α,Real α) => [α] -> [α] -> Double
r² = ((^2).) . r

-- test = do
--   a <- getLine
--   b <- getLine
--   putStrLn $ a ++ b

data AnovaTable = AnovaTable {
                             df :: [Int],
                             ssq :: [Double],
                             ms :: [Double],
                             fStat :: Double,
                             fCritical :: Double,
                             isSignificant :: Bool,
                             significance:: Int}
--                  deriving Show

instance Show AnovaTable where
  show (AnovaTable df ssq ms fStat fCritical _ significance) =
    unlines . (++[show significance,legend]). map (concatWith "\t") $ [a,b,c,d,e]
    where f s f' l = ([s]++) $ map f' l
          f' s l = f s (printf "%.4f") l
          a = f "Df" show df
          b = f' "SS" ssq
          c = f' "MS" ms
          d = f' "F_value" [fStat] ++ [sl]
          e = f' "F_crit" [fCritical]
          sl = [" ",".","*","**","***"] !! significance
          legend = "Signif.codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
--oneWay_anova :: (Real α,Fractional α,) => Double -> [[α]] -> AnovaTable
oneWay_anova α cols =
  AnovaTable {
             df = [dfb,dfw,dft],
             ssq = [ssb,ssw,sst],
             ms = [msb,msw],
             fStat = fstat,
             fCritical = fcritical,
             isSignificant = fstat > fcritical,
             significance  = significanceLevel}
  where l = map (\(Col _ v) -> v) cols
        vs = concat l
        mean = µ vs
        n = count vs
        means = map µ l
        counts = map count l
        dfb = count l - 1
        dft = count vs - 1
        dfw = dft - dfb
        ssb = sum $ zipWith (\m c -> c*(m-mean)^2 ) means counts
        ssw = sum . map ss $ l
        sst = ssb + ssw
        msb = ssb / (fromIntegral dfb)
        msw = ssw / (fromIntegral dfw)
        (Just fcritical) = f_stat_critical α dfb dfw
        (Just criticalList) =
          sequence $ map (\a -> f_stat_critical a dfb dfw) ([0.9,0.95,0.99,0.999] :: [Fraction])
        fstat = msb / msw

        significanceLevel = count . filter (fstat>) $ criticalList

test_data_anova = [[7,4,6,8,6,6,2,9],
                   [5,5,3,4,4,7,2,2],
                   [2,4,7,1,2,1,5,5]] :: [[Double]]

-- example:
-- (oneWay_anova 0.95 . dataTableFromTable 4 3 . tail . readTable) <$> readFile (path_data_stats ++ "fisher_iris.txt")

--------------
{- SAMPLING -}
--------------

-- data Sampler = Sampler {
--   randomgen :: RandomGen,

--         --                }
-- pick1 :: (RandomGen γ) => State γ [α] -> State γ (α,[α])
-- pick1 ( (pop,γ)) = do
--   k <- random g
--   let n = count pop
--       (hd,(t:ts)) = splitAt (k `mod` n) pop
--   return (t,hd++ts)

-- pick1 :: (RandomGen γ) => State γ [α] -> State γ (α,[α])
-- pick1 pop = do
--   k <- random
--   let (hd,(t:ts)) = splitAt (k `mod` n) pop
--   return (t,hd++ts)

-- pick1 pop g = ((t,hd++ts),g')
--   where (k,g') = random g
--         n = count pop
--         (hd,(t:ts)) = splitAt (k `mod` n) pop

-- pick One element from pop at random
pick1 :: [α] -> State StdGen (α,[α])
pick1 pop = do
  g <- get
  let (k,g') = random g
      n = count pop
      (hd,(t:ts)) = splitAt (k `mod` n) pop
  put g'
  return (t,hd++ts)

-- get a sample from pop (without replacement)
-- sample :: (Integral β,RandomGen δ) => δ -> β -> [α] -> [α]
sample 0 pop = return ([],pop)
sample k pop = do
  (e,pop') <- pick1 pop
  (l, pop'') <- sample (k-1) pop'
  return (e:l, pop'')

-- get n samples of size k from pop
samples 0 _ pop = return []
samples n k pop = do
  (s, _) <- sample k pop
  l <- samples (n-1) k pop
  return $ s : l

freqs :: (Foldable f , Eq α) => f α -> [(α,Int)]
freqs = quickSortBy snd . foldl' f []
  where f [] e = [(e,1)]
        f (x@(e',k):xs) e | e' == e   = ((e,k+1):xs)
                          | otherwise = x : f xs e

max_stars_repo_name: mikaunix/hlib

max_issues_repo_name: mikaunix/hlib

max_forks_repo_name: mikaunix/hlib


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MultiParamTypeClasses #-}
module Statistics.BBVI.Observed
  ( Obs(..)
  , defaultObs
  )
where

import qualified Data.Vector                   as V
import           Statistics.BBVI.Class
import           Statistics.BBVI.Propagator     ( DistCell(..)
                                                , SampleVector
                                                , SampleDouble
                                                )
import           System.Random.MWC              ( uniformR )

-- | Dummy distribution representing a vector of observations, for use
-- with building "observation" distribution cells.  these can be used
-- to easily/selectively subsample (and resample) observations using
-- existing typeclass methods (useful with stochatic gradient
-- updates).
newtype Obs a = O (V.Vector a) deriving (Show, Eq, Ord, Read)

-- | helper function to build a distribution cell of observations
defaultObs
  :: V.Vector a -- ^ observations
  -> DistCell (Obs a)
defaultObs d = (Node 1 V.empty (O d))

instance DistUtil (Obs Double) where
  nParams _x = 0
  toParamVector _ = V.empty
  fromParamVector _ = O V.empty

instance Dist (Obs Double) SampleDouble where
  resample (O d) gen = return . (d V.!) =<< uniformR (0, (V.length d - 1)) gen
  logProb _d _x = 0
  paramGradOfLogQ _d _x = V.empty

instance DistUtil (Obs SampleVector) where
  nParams _x = 0
  toParamVector _ = V.empty
  fromParamVector _ = O V.empty

instance Dist (Obs SampleVector) SampleVector where
  resample (O d) gen = return . (d V.!) =<< uniformR (0, (V.length d - 1)) gen
  logProb _d _x = 0
  paramGradOfLogQ _d _x = V.empty

ext: hs

max_stars_repo_name: massma/propagator-bbvi

max_issues_repo_name: massma/propagator-bbvi

max_forks_repo_name: massma/propagator-bbvi


text:
-- |
-- Copyright  : (c) 2012 Daniël de Kok
-- License    : Apache 2
--
-- Maintainer : Daniël de Kok <me@danieldk.eu>
-- Stability  : experimental
--
-- Approximate randomization test (Noreen, 1989)

{-# LANGUAGE DoAndIfThenElse #-}

module Main where

import           Control.Monad (liftM, when)
import           Control.Monad.Mersenne.Random (evalRandom)
import           Data.List (intersperse)
import qualified Data.Vector.Generic as VG
import qualified Data.Vector.Unboxed as V
import           Data.Word (Word64)
import           Statistics.Test.ApproxRand
import           Statistics.Test.Types (TestType(..))
import           Statistics.Types (Sample)
import           System.Console.GetOpt
import           System.Environment (getArgs)
import           System.Exit (exitFailure)
import           System.Random.Mersenne.Pure64 (PureMT, newPureMT, pureMT)
import           Text.Printf (printf)

import           ChartHistogram
import           SampleIO
import           TextHistogram

main :: IO ()
main = do
  -- Read command-line options and arguments.
  (opts, args) <- getOptions

  -- Read score files
  let col = pred $ optColumn opts
  v1 <- liftM V.fromList $ readFileCol (head args) col
  v2 <- liftM V.fromList $ readFileCol (args !! 1) col

  let stat = optTestStatistic opts

  prng <- case optPRNGSeed opts of
    Just seed -> return $ pureMT seed
    Nothing   -> newPureMT

  if optPrintStats opts then
    printStats opts stat prng v1 v2
  else
    applyTest opts stat prng v1 v2

applyTest :: Options -> TestStatistic -> PureMT -> Sample ->
  Sample -> IO ()
applyTest opts stat prng v1 v2 = do
  putStrLn $ printf "Iterations: %d" $ optIterations opts
  putStrLn $ printf "Sample sizes: %d %d" (V.length v1) (V.length v2)

  let testOptions = TestOptions (optTestType opts)
                                stat
                                (optIterations opts)
                                (optSigP opts)

  let testType = optTestType opts

  let pTest = optSigP opts
  let pTail = case testType of
                OneTailed -> pTest
                TwoTailed -> pTest / 2

  -- Test information
  putStrLn $ "Test type: " ++ show testType
  putStrLn $ printf "Test significance: %f" pTest
  putStrLn $ printf "Tail significance: %f" pTail

  -- Approximate randomization testing.
  let test = approxRandTest testOptions v1 v2
  let result = evalRandom test prng

  -- Print test statistic for the samples.
  putStrLn $ printf "Test statistic: %f" $ trStat result

  case trSignificance result of
    Significant    p -> putStrLn $ printf "Significant: %f" p
    NotSignificant p -> putStrLn $ printf "Not significant: %f" p

  when (optPrintHistogram opts) $ do
    putStrLn ""
    printHistogram 21 result
  case (optWriteHistogram opts) of
    Just fn ->
      writeHistogram testOptions 31 result fn
    Nothing ->
      return ()

printStats :: Options -> TestStatistic -> PureMT -> Sample ->
  Sample -> IO ()
printStats opts stat prng v1 v2 =
  VG.mapM_ (putStrLn . printf "%f") $
    evalRandom (approxRandStats stat (optIterations opts) v1 v2) prng

data Options = Options {
  optColumn         :: Int,
  optPrintHistogram :: Bool,
  optIterations     :: Int,
  optPRNGSeed       :: Maybe Word64,
  optPrintStats     :: Bool,
  optSigP           :: Double,
  optTestStatistic  :: TestStatistic,
  optTestType       :: TestType,
  optWriteHistogram :: Maybe String
}

defaultOptions :: Options
defaultOptions = Options {
  optColumn         = 1,
  optIterations     = 10000,
  optPRNGSeed       = Nothing,
  optPrintHistogram = False,
  optPrintStats     = False,
  optSigP           = 0.01,
  optTestStatistic  = meanDifference,
  optTestType       = TwoTailed,
  optWriteHistogram = Nothing
}

options :: [OptDescr (Options -> Options)]
options =
  chartHistogramOption backendFormats : mandatoryOptions

mandatoryOptions :: [OptDescr (Options -> Options)]
mandatoryOptions =
  [ Option ['c'] ["column"]
      (ReqArg (\arg opt -> opt { optColumn = read arg }) "NUMBER")
      "column number (starting at 1)",
    Option ['h']    ["print-histogram"]
      (NoArg (\opt -> opt { optPrintHistogram = True }))
      "print a histogram of randomized sample statistics",
    Option ['i'] ["iterations"]
      (ReqArg (\arg opt -> opt { optIterations = read arg }) "NUMBER")
      "number of iterations",
    Option ['o'] ["one-tailed"]
      (NoArg (\opt -> opt { optTestType = OneTailed }))
      "perform a one-tailed test",
    Option ['p'] []
      (ReqArg (\arg opt -> opt {optSigP = read arg }) "NUMBER")
      "significant p-value",
    Option []    ["print-scores"]
      (NoArg (\opt -> opt { optPrintStats = True }))
      "output statistics of permuted vectors",
    Option ['s'] ["seed"]
      (ReqArg (\arg opt -> opt { optPRNGSeed = Just $ read arg}) "NUMBER")
      "pseudorandom number generator seed",
    Option ['t'] ["test-statistic"]
      (ReqArg (\arg opt -> opt { optTestStatistic = parseStatistic arg}) "NAME")
      "test statistic (mean_diff, var_ratio)"
  ]

chartHistogramOption :: [String] -> OptDescr (Options -> Options)
chartHistogramOption formats =
  Option ['w'] ["write-histogram"]
    (ReqArg (\arg opt -> opt { optWriteHistogram = Just arg}) "FILENAME")
    $ printf "write a histogram (supported file extensions: %s)" formatStr
  where
    formatStr = concat $ intersperse ", " formats

getOptions :: IO (Options, [String])
getOptions = do
  args <- getArgs
  case getOpt Permute options args of
    (actions, nonOpts, [])   -> do
      when (length nonOpts /= 2) usageExit
      let opts = foldl (flip ($)) defaultOptions actions
      return (opts, nonOpts)
    (_,    _,          _)    ->
      usageExit
  where
  usageExit = do
    putStrLn $ usageInfo header options
    exitFailure
    where
      header = "Usage: approx-rand-test [OPTION...] scores scores2"

parseStatistic :: String -> TestStatistic
parseStatistic "mean_diff" = meanDifference
parseStatistic "var_ratio" = varianceRatio
parseStatistic _           = error "Unknown test statistic"

max_stars_repo_name: danieldk/approx-rand-test

max_issues_repo_name: danieldk/approx-rand-test

max_forks_repo_name: danieldk/approx-rand-test


text:
module Algebra
    (Adjoin (..)
    ,FreeGroup (..)
    ,Group (..)
    ,Z2
    ,plugIn
    ,Algebra
    ,represent
    ,cutoff
    ,module N
    ) where

import Algebra.Adjoin
import Algebra.Group
import Algebra.FreeGroup
import Algebra.Z2
import Numeric.LinearAlgebra as N

type Algebra = Adjoin Z2 FreeGroup
{--type Algegra' = [(Z2,Vector Z,FreeGroup)]

plush :: (Z2,Vector Z,FreeGroup) -> Algebra' -> Algebra'
plush e [] = [e]
plush e1@(z1,v1,g1) (e2@(z2,v2,g2):ls) = if v1 == v2 && g1 == g2 then (if z1+z2 == 0 then (0,N.fromList $ map (\_ -> 0) $ N.tolist v1,mempty) else (z1+z2,v1,g1)):ls else e2:(plush e1 ls)

instance Num Algebra' where
    l1 + l2 = map (\e -> plush e l2) l1
    l1 * l2 = sum $ map (\(z1,v1,g1) -> map (\(z2,v2,g2) -> (z1*z2,v1+v2,g1<>g2)) l2) l1
    abs = map (\(z,v,g) -> (abs z,v,g))
    signum = map (\(z,v,g) -> (signum z,v,g))
--    fromInteger x = [(fromInteger x,
-}
cutoff :: R
cutoff = 0.0000000001

plugIn :: (Eq f,Fractional f) => (Char -> Adjoin f FreeGroup) -> Adjoin f FreeGroup -> Adjoin f FreeGroup
plugIn f = adjPlugIn (groupPlugIn f)

represent :: [Char] -> Algebra -> Maybe [Vector Z]
represent [] _ = Nothing
represent elems expr = do
                        { chain <- pairPlusChain expr
                        ; vectors <- mapM (\(_,g) -> toPseudoVector g) chain
                        ; return $ map (\l -> (length elems) |> [toEnum $ maybe 0 id $ lookup (elems !! i) l | i <- [0..]]) vectors
                        }

max_stars_repo_name: Creatorri/Legendrian-Knots-UROP

max_issues_repo_name: Creatorri/Legendrian-Knots-UROP

max_forks_repo_name: Creatorri/Legendrian-Knots-UROP


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TypeApplications #-}

module ProblemSpec where

import Commons
import Control.Applicative (liftA2)
import Control.Concurrent
import Control.Monad (forM_, replicateM, replicateM_, unless, when)
import Data.Array
import Data.Complex (Complex (..))
import qualified Data.IntMap.Strict as IM
import Data.List (intercalate, sort)
import Data.List.Split (splitOn)
import qualified Data.Map.Strict as Map
import Data.Map.Strict (fromList)
import Data.Maybe (fromJust)
import qualified Data.Set as Set
import qualified Data.Text as T
import qualified Data.Text.IO as TIO
import Debug.Trace (traceShowId)
import GHC.IO.Exception (ExitCode (..))
import HashedExpression.Derivative
import HashedExpression.Internal.CollectDifferential
import HashedExpression.Internal.Expression
import HashedExpression.Internal.Inner
import HashedExpression.Internal.Node
import HashedExpression.Internal.Normalize (normalize)
import HashedExpression.Internal.Utils
import HashedExpression.Interp
import HashedExpression.Operation
import HashedExpression.Prettify (showExp, showExpDebug)
import HashedExpression.Problem
import HashedExpression.Value
import System.Process (readProcess, readProcessWithExitCode)
import Test.HUnit
import Test.Hspec
import Test.QuickCheck
import Var
import Prelude hiding ((^))
import qualified Prelude

---- |
--isOneAfterAnother :: MemMap -> [Int] -> Bool
--isOneAfterAnother memMap nIds = all isOk xs
--  where
--    xs = zip nIds (tail nIds)
--    isOk (cur, nxt) =
--      let (offsetCur, _, shapeCur) =
--            fromJust $ IM.lookup cur (entryMap memMap)
--          (offsetNxt, _, _) = fromJust $ IM.lookup nxt (entryMap memMap)
--       in offsetCur + Prelude.product shapeCur == offsetNxt

-- |
prop_constructProblemNoConstraint :: SuiteScalarR -> Expectation
prop_constructProblemNoConstraint (Suite exp valMap) = do
  let names = Map.keys valMap
      pdMap = partialDerivativeMaps . collectDifferentials . exteriorDerivative (Set.fromList names) $ exp
      constructResult = constructProblem exp names (Constraint [])
  case constructResult of
    NoVariables -> return () -- it is possible that the random expression doesn't have any variables
    ProblemInvalid reason ->
      assertFailure $ "Can't construct problem: " ++ reason
    ProblemValid Problem {..} -> do
      let vars = map varName variables
      vars `shouldBe` Map.keys pdMap -- vars should be keys of partial differential map
      let ok variable
            | Just pId <- Map.lookup (varName variable) pdMap,
              pId == partialDerivativeId variable =
              True
            | otherwise = False
      assertBool "partial derivative ids aren't correct" $ all ok variables

--      assertBool "variables are not allocated consecutively" $
--        isOneAfterAnother memMap (map nodeId variables)

-- |
makeValidBoxConstraint :: (String, Shape) -> IO ConstraintStatement
makeValidBoxConstraint (name, shape) =
  case shape of
    [] -> do
      let x = variable name
      val1 <- VScalar <$> generate arbitrary
      val2 <- VScalar <$> generate arbitrary
      generate $
        elements [x .<= val1, x .>= val2, x `between` (val1, val2)]
    [size] -> do
      let x = variable1D @Default1D name
      val1 <- V1D . listArray (0, size - 1) <$> generate (vectorOf size arbitrary)
      val2 <- V1D . listArray (0, size - 1) <$> generate (vectorOf size arbitrary)
      generate $ elements [x .<= val1, x .>= val2, x `between` (val1, val2)]
    [size1, size2] -> do
      let x = variable2D @Default2D1 @Default2D2 name
      val1 <- V2D . listArray ((0, 0), (size1 - 1, size2 - 1)) <$> generate (vectorOf (size1 * size2) arbitrary)
      val2 <- V2D . listArray ((0, 0), (size1 - 1, size2 - 1)) <$> generate (vectorOf (size1 * size2) arbitrary)
      generate $ elements [x .<= val1, x .>= val2, x `between` (val1, val2)]

--        [size1, size2, size3] -> do TODO -- add 3D for tests
--            val1 <-
--                V3D . listArray ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1)) <$>
--                generate (vectorOf (size1 * size2 * size3) arbitrary)
--            val2 <-
--                V3D . listArray ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1)) <$>
--                generate (vectorOf (size1 * size2 * size3) arbitrary)
--            generate $
--                elements [x .<= val1, x .>= val2, x `between` (val1, val2)]

-- |
prop_constructProblemBoxConstraint :: SuiteScalarR -> Expectation
prop_constructProblemBoxConstraint (Suite exp valMap) = do
  let names = Map.keys valMap
  let df@(Expression dfN dfMp) =
        collectDifferentials . exteriorDerivative (Set.fromList names) $ exp
      pdMap = partialDerivativeMaps df
      vars = Map.keys pdMap
      varsWithShape = Map.toList $ Map.map (`retrieveShape` dfMp) pdMap
  bcs <- mapM makeValidBoxConstraint varsWithShape
  sampled <- generate $ sublistOf bcs
  let constraints = Constraint sampled
  let constructResult = constructProblem exp names constraints
  case constructResult of
    NoVariables -> return () -- it is possible that the random expression doesn't have any variables
    ProblemInvalid reason ->
      assertFailure $ "Can't construct problem: " ++ reason
    ProblemValid Problem {..} -> do
      let vars = map varName variables
      vars `shouldBe` Map.keys pdMap -- vars should be keys of partial differential map
      let ok variable
            | Just pId <- Map.lookup (varName variable) pdMap,
              pId == partialDerivativeId variable =
              True
            | otherwise = False
      assertBool "partial derivative ids aren't correct" $
        all ok variables
      --      assertBool "variables are not allocated consecutively" $
      --        isOneAfterAnother memMap (map nodeId variables)
      case (sampled, boxConstraints) of
        (_ : _, []) ->
          assertFailure
            "Valid box constraints but not appear in the problem"
        (_, bs) -> length sampled `shouldBe` length bs

-- |
makeValidScalarConstraint :: IO ConstraintStatement
makeValidScalarConstraint = do
  sc <-
    fst <$> generate (sized (genScalarR @Default1D @Default2D1 @Default2D2))
  val1 <- VScalar <$> generate arbitrary
  val2 <- VScalar <$> generate arbitrary
  generate $ elements [sc .<= val1, sc .>= val2, sc `between` (val1, val2)]

-- |
prop_constructProblemScalarConstraints :: SuiteScalarR -> Expectation
prop_constructProblemScalarConstraints (Suite exp valMap) = do
  let names = Map.keys valMap
  let df@(Expression dfN dfMp) =
        collectDifferentials . exteriorDerivative (Set.fromList names) $ exp
      pdMap = partialDerivativeMaps df
      vars = Map.keys pdMap
      varsWithShape = Map.toList $ Map.map (`retrieveShape` dfMp) pdMap
  -- box constraints
  bcs <- mapM makeValidBoxConstraint varsWithShape
  sampled <- generate $ sublistOf bcs
  -- scalar constraints
  numScalarConstraint <- generate $ elements [2 .. 4]
  scc <- replicateM numScalarConstraint makeValidScalarConstraint
  let constraints = Constraint $ sampled ++ scc
  let constructResult = constructProblem exp names constraints
  case constructResult of
    NoVariables -> return () -- it is possible that the random expression doesn't have any variables
    ProblemInvalid reason ->
      assertFailure $ "Can't construct problem: " ++ reason
    ProblemValid Problem {..} -> do
      let vars = map varName variables
      vars `shouldBe` Map.keys pdMap -- vars should be keys of partial differential map
      let ok variable
            | Just pId <- Map.lookup (varName variable) pdMap,
              pId == partialDerivativeId variable =
              True
            | otherwise = False
      assertBool "partial derivative ids aren't correct" $
        all ok variables
      --      assertBool "variables are not allocated consecutively" $
      --        isOneAfterAnother memMap (map nodeId variables)
      case (scc, scalarConstraints) of
        ([], _) -> return ()
        (_ : _, []) ->
          assertFailure
            "Having scalar constraints but not present in problem"
        (_, sConstraints) -> do
          let isOk sc =
                length (constraintPartialDerivatives sc)
                  `shouldBe` length vars
          assertBool "Empty constraint ?" $ not (null sConstraints)
          mapM_ isOk sConstraints

-- | List of hand-written problems and the expected result whether this problem is valid
problemsRepo :: [(ProblemResult, Bool)]
problemsRepo =
  [ ( let [x, y, z, t] = map (variable2D @128 @128) ["x", "y", "z", "t"]
          f = x <.> y + z <.> t
          constraints =
            Constraint
              [ x .>= VFile (TXT "x_lb.txt"),
                y .<= VFile (TXT "y_ub.txt"),
                x <.> z .>= VScalar 3
              ]
          vars = ["x", "y", "z", "t"]
       in constructProblem f vars constraints,
      True
    ),
    ( let [x, y, z, t] = map (variable2D @128 @128) ["x", "y", "z", "t"]
          f = x <.> y + z <.> t
          constraints =
            Constraint
              [ x .>= VScalar 5,
                y .<= VFile (TXT "y_ub.txt"),
                x <.> z .>= VScalar 3
              ]
          vars = ["x", "y", "z", "t"]
       in constructProblem f vars constraints,
      False
    ),
    ( let [x, y, z, t] = map (variable2D @128 @128) ["x", "y", "z", "t"]
          f = x <.> y + z <.> t
          constraints =
            Constraint [x .>= VNum 5, y .<= VNum 10, x <.> z .>= VNum 18]
          vars = ["x", "y", "z", "t"]
       in constructProblem f vars constraints,
      True
    )
  ]

spec :: Spec
spec =
  describe "Hash Solver spec " $ do
    specify "test hand-written problems" $ forM_ problemsRepo $ \(problemResult, expected) ->
      case (problemResult, expected) of
        (ProblemInvalid _, True) ->
          assertFailure "This problem is valid but fail to construct"
        (ProblemValid _, False) ->
          assertFailure "This problem is invalid but success to construct"
        _ -> return ()
    specify "valid problem should be constructed successfully" $
      property prop_constructProblemNoConstraint
    specify "valid box constrained problem should be constructed successfully" $
      property prop_constructProblemBoxConstraint
    specify "valid scalar constraints problem should be successfully successfully" $
      property prop_constructProblemScalarConstraints

max_stars_repo_name: Turboscient/HashedExpression

max_issues_repo_name: Turboscient/HashedExpression

max_forks_repo_name: Turboscient/HashedExpression


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE GADTs                    #-}
{-# LANGUAGE Strict                   #-}
module Grenade.Layers.Internal.Update (
    descendMatrix
  , descendVector
  , MatrixInputValues (..)
  , MatrixResult (..)
  , VectorInputValues (..)
  , VectorResult (..)
  , descendMatrixV
  , descendVectorV
  , MatrixInputValuesV (..)
  , MatrixResultV (..)
  , VectorInputValuesV (..)
  , VectorResultV (..)
  ) where

import           Control.Parallel.Strategies
import           Data.Maybe                   (fromJust, fromMaybe)
import qualified Data.Vector.Storable         as U (unsafeFromForeignPtr0,
                                                    unsafeToForeignPtr0)
import qualified Data.Vector.Storable         as V
import           Foreign                      (mallocForeignPtrArray, withForeignPtr)
import           Foreign.Ptr                  (Ptr)
import           GHC.TypeLits
import           Numeric.LinearAlgebra        (Vector, flatten)
import qualified Numeric.LinearAlgebra.Devel  as U
import           Numeric.LinearAlgebra.Static
import           System.IO.Unsafe             (unsafePerformIO)

import           Grenade.Core.Optimizer
import           Grenade.Layers.Internal.CUDA
import           Grenade.Types

data MatrixInputValues rows columns
  = MatrixValuesSGD
      !(L rows columns) -- ^ current weights
      !(L rows columns) -- ^ gradients
      !(L rows columns) -- ^ last update (old momentum)
  | MatrixValuesAdam
      !Int              -- ^ Step
      !(L rows columns) -- ^ current weights
      !(L rows columns) -- ^ gradients
      !(L rows columns) -- ^ current m
      !(L rows columns) -- ^ current v

data MatrixInputValuesV
  = MatrixValuesSGDV
      !(V.Vector RealNum) -- ^ current weights
      !(V.Vector RealNum) -- ^ gradients
      !(V.Vector RealNum) -- ^ old momentum
  | MatrixValuesAdamV
      !Int                -- ^ Step
      !(V.Vector RealNum) -- ^ current weights
      !(V.Vector RealNum) -- ^ gradients
      !(V.Vector RealNum) -- ^ current m
      !(V.Vector RealNum) -- ^ current v

data MatrixResult rows columns
  = MatrixResultSGD
      { matrixActivations :: !(L rows columns) -- ^ new activations (weights)
      , matrixMomentum    :: !(L rows columns) -- ^ new momentum
      }
  | MatrixResultAdam
      { matrixActivations :: !(L rows columns) -- ^ new activations (weights)
      , matrixM           :: !(L rows columns) -- ^ new m
      , matrixV           :: !(L rows columns) -- ^ new v
      }

data MatrixResultV
  = MatrixResultSGDV
      { matrixActivationsV :: !(V.Vector RealNum) -- ^ new activations (weights)
      , matrixMV           :: !(V.Vector RealNum) -- ^ new m
      }
  | MatrixResultAdamV
      { matrixActivationsV :: !(V.Vector RealNum) -- ^ new activations (weights)
      , matrixMV           :: !(V.Vector RealNum) -- ^ new m
      , matrixVV           :: !(V.Vector RealNum) -- ^ new v
      }

data VectorInputValues r
  = VectorValuesSGD
      !(R r) -- ^ current weights
      !(R r) -- ^ gradients
      !(R r) -- ^ last update (old momentum)
  | VectorValuesAdam
      !Int   -- ^ Step
      !(R r) -- ^ current weights
      !(R r) -- ^ current gradients
      !(R r) -- ^ current m
      !(R r) -- ^ current v

data VectorInputValuesV
  = VectorValuesSGDV
      !(V.Vector RealNum) -- ^ current weights
      !(V.Vector RealNum) -- ^ current gradients
      !(V.Vector RealNum) -- ^ current m
  | VectorValuesAdamV
      !Int                -- ^ Step
      !(V.Vector RealNum) -- ^ current weights
      !(V.Vector RealNum) -- ^ current gradients
      !(V.Vector RealNum) -- ^ current m
      !(V.Vector RealNum) -- ^ current v

data VectorResult r
  = VectorResultSGD
      { vectorBias     :: !(R r) -- ^ new activations (bias)
      , vectorMomentum :: !(R r) -- ^ new momentum
      }
  | VectorResultAdam
      { vectorBias :: !(R r) -- ^ new activations (bias)
      , vectorM    :: !(R r) -- ^ new m
      , vectorV    :: !(R r) -- ^ new v
      }

data VectorResultV
  = VectorResultSGDV
      { vectorBiasV :: !(V.Vector RealNum) -- ^ new activations (bias)
      , vectorMV    :: !(V.Vector RealNum) -- ^ new m
      }
  | VectorResultAdamV
      { vectorBiasV :: !(V.Vector RealNum) -- ^ new activations (bias)
      , vectorMV    :: !(V.Vector RealNum) -- ^ new m
      , vectorVV    :: !(V.Vector RealNum) -- ^ new v
      }

descendMatrix :: (KnownNat rows, KnownNat columns) => Optimizer o -> MatrixInputValues rows columns -> MatrixResult rows columns
descendMatrix (OptSGD rate momentum regulariser) (MatrixValuesSGD weights gradient lastUpdate) =
  let (rows, cols) = size weights
      len          = rows * cols
      -- Most gradients come in in ColumnMajor,
      -- so we'll transpose here before flattening them
      -- into a vector to prevent a copy.
      --
      -- This gives ~15% speed improvement for LSTMs.
      weights'     = flatten . tr . extract $ weights
      gradient'    = flatten . tr . extract $ gradient
      lastUpdate'  = flatten . tr . extract $ lastUpdate
      (vw, vm)     = descendUnsafeSGD len rate momentum regulariser weights' gradient' lastUpdate'

      -- Note that it's ColumnMajor, as we did a transpose before
      -- using the internal vectors.
      mw           = U.matrixFromVector U.ColumnMajor rows cols vw
      mm           = U.matrixFromVector U.ColumnMajor rows cols vm
  in  MatrixResultSGD (fromJust . create $ mw) (fromJust . create $ mm)
descendMatrix (OptAdam alpha beta1 beta2 epsilon lambda) (MatrixValuesAdam step weights gradient m v) =
  let (rows, cols) = size weights
      len          = rows * cols
      -- Most gradients come in in ColumnMajor,
      -- so we'll transpose here before flattening them
      -- into a vector to prevent a copy.
      --
      -- This gives ~15% speed improvement for LSTMs.
      weights'  = flatten . tr . extract $ weights
      gradient' = flatten . tr . extract $ gradient
      m'        = flatten . tr . extract $ m
      v'        = flatten . tr . extract $ v
      (vw, vm, vv)     = descendUnsafeAdamGPU len step alpha beta1 beta2 epsilon lambda weights' gradient' m' v'

      -- Note that it's ColumnMajor, as we did a transpose before
      -- using the internal vectors.
      mw           = U.matrixFromVector U.ColumnMajor rows cols vw
      mm           = U.matrixFromVector U.ColumnMajor rows cols vm
      mv           = U.matrixFromVector U.ColumnMajor rows cols vv
  in  MatrixResultAdam (fromJust . create $ mw) (fromJust . create $ mm) (fromJust . create $ mv)
descendMatrix opt _ = error $ "optimzer does not match to MatrixInputValues in implementation! Optimizer: " ++ show opt

descendMatrixV :: Optimizer o -> MatrixInputValuesV -> MatrixResultV
descendMatrixV (OptSGD rate momentum regulariser) (MatrixValuesSGDV weights gradient lastUpdate) =
  let len          = V.length weights
      (vw, vm)     = descendUnsafeSGD len rate momentum regulariser weights gradient lastUpdate
  in  MatrixResultSGDV vw vm
descendMatrixV (OptAdam alpha beta1 beta2 epsilon lambda) (MatrixValuesAdamV step weights gradient m v) =
  let len          = V.length weights
      (vw, vm, vv) = descendUnsafeAdamGPU len step alpha beta1 beta2 epsilon lambda weights gradient m v
  in  MatrixResultAdamV vw vm vv
descendMatrixV opt _ = error $ "optimzer does not match to MatrixInputValues in implementation! Optimizer: " ++ show opt

descendVector :: (KnownNat r) => Optimizer o -> VectorInputValues r -> VectorResult r
descendVector (OptSGD rate momentum regulariser) (VectorValuesSGD weights gradient lastUpdate) =
  let len          = size weights
      weights'     = extract weights
      gradient'    = extract gradient
      lastUpdate'  = extract lastUpdate
      (vw, vm)     = descendUnsafeSGD len rate momentum regulariser weights' gradient' lastUpdate'
  in  VectorResultSGD (fromJust $ create vw) (fromJust $ create vm)
descendVector (OptAdam alpha beta1 beta2 epsilon lambda) (VectorValuesAdam step weights gradient m v) =
  let len       = size weights
      weights'  = extract weights
      gradient' = extract gradient
      m'        = extract m
      v'        = extract v
      (vw, vm, vv)     = descendUnsafeAdamGPU len step alpha beta1 beta2 epsilon lambda weights' gradient' m' v'
  in  VectorResultAdam (fromJust $ create vw) (fromJust $ create vm) (fromJust $ create vv)
descendVector opt _ = error $ "optimzer does not match to VectorInputValues in implementation! Optimizer: " ++ show opt

descendVectorV :: Optimizer o -> VectorInputValuesV -> VectorResultV
descendVectorV (OptSGD rate momentum regulariser) (VectorValuesSGDV weights gradient lastUpdate) =
  let len          = V.length weights
      (vw, vm)     = descendUnsafeSGD len rate momentum regulariser weights gradient lastUpdate
  in  VectorResultSGDV vw vm
descendVectorV (OptAdam alpha beta1 beta2 epsilon lambda) (VectorValuesAdamV step weights gradient m v) =
  let len = V.length weights
      (vw, vm, vv) = descendUnsafeAdamGPU len step alpha beta1 beta2 epsilon lambda weights gradient m v
   in VectorResultAdamV vw vm vv
descendVectorV opt _ = error $ "optimzer does not match to VectorInputValues in implementation! Optimizer: " ++ show opt

-- -- | Caching of data
-- type CacheKey = (LookupType, ProxyType, StateFeatures)

-- cacheMVar :: MVar (M.Map CacheKey [Values])
-- cacheMVar = unsafePerformIO $ newMVar mempty
-- {-# NOINLINE cacheMVar #-}

-- emptyCache :: MonadIO m => m ()
-- emptyCache = liftIO $ modifyMVar_ cacheMVar (const mempty)

-- addCache :: (MonadIO m) => CacheKey -> [Values] -> m ()
-- addCache k val = liftIO $ modifyMVar_ cacheMVar (return . M.insert k val)

-- lookupCache :: (MonadIO m) => CacheKey -> m (Maybe [Values])
-- lookupCache k = liftIO $ (M.lookup k =<<) <$> tryReadMVar cacheMVar

-- -- | Get output of function f, if possible from cache according to key (st).
-- cached :: (MonadIO m) => (LookupType, ProxyType, StateFeatures) -> m [Values] -> m [Values]
-- cached st ~f = do
--   c <- lookupCache st
--   case c of
--     Nothing -> do
--       res <- f
--       res `seq` addCache st res
--       return res
--     Just res -> return res

descendUnsafeSGD :: Int -> RealNum -> RealNum -> RealNum -> Vector RealNum -> Vector RealNum -> Vector RealNum -> (Vector RealNum, Vector RealNum)
descendUnsafeSGD len rate momentum regulariser weights gradient lastUpdate =
  unsafePerformIO $ do

    outWPtr <- mallocForeignPtrArray len
    outMPtr <- mallocForeignPtrArray len
    let (wPtr, _) = U.unsafeToForeignPtr0 weights
    let (gPtr, _) = U.unsafeToForeignPtr0 gradient
    let (lPtr, _) = U.unsafeToForeignPtr0 lastUpdate

    withForeignPtr wPtr $ \wPtr' ->
      withForeignPtr gPtr $ \gPtr' ->
        withForeignPtr lPtr $ \lPtr' ->
          withForeignPtr outWPtr $ \outWPtr' ->
            withForeignPtr outMPtr $ \outMPtr' ->
              descend_sgd_cpu len rate momentum regulariser wPtr' gPtr' lPtr' outWPtr' outMPtr'

    return (U.unsafeFromForeignPtr0 outWPtr len, U.unsafeFromForeignPtr0 outMPtr len)
{-# NOINLINE descendUnsafeSGD #-}

descendUnsafeAdam ::
     Int -- Len
  -> Int -- Step
  -> RealNum -- Alpha
  -> RealNum -- Beta1
  -> RealNum -- Beta2
  -> RealNum -- Epsilon
  -> RealNum -- Lambda
  -> Vector RealNum -- Weights
  -> Vector RealNum -- Gradient
  -> Vector RealNum -- M
  -> Vector RealNum -- V
  -> (Vector RealNum, Vector RealNum, Vector RealNum)
descendUnsafeAdam len step alpha beta1 beta2 epsilon lambda weights gradient m v =
  unsafePerformIO $ do
  V.unsafeWith weights $ \wPtr' ->
    V.unsafeWith gradient $ \gPtr' ->
      V.unsafeWith m $ \mPtr' ->
        V.unsafeWith v $ \vPtr' ->
           descend_adam_cpu len step alpha beta1 beta2 epsilon lambda wPtr' gPtr' mPtr' vPtr'
  return (weights, m, v)
{-# NOINLINE descendUnsafeAdam #-}

descendUnsafeAdamGPU ::
     Int -- Len
  -> Int -- Step
  -> RealNum -- Alpha
  -> RealNum -- Beta1
  -> RealNum -- Beta2
  -> RealNum -- Epsilon
  -> RealNum -- Lambda
  -> Vector RealNum -- Weights
  -> Vector RealNum -- Gradient
  -> Vector RealNum -- M
  -> Vector RealNum -- V
  -> (Vector RealNum, Vector RealNum, Vector RealNum)
descendUnsafeAdamGPU len step alpha beta1 beta2 epsilon lambda weights gradient m v
  --  | useCuda len = -- We need some size, otherwise the overhead is too big and processing using BLAS is faster
  --   fromMaybe
  --     (descendUnsafeAdam len step alpha beta1 beta2 epsilon lambda weights gradient m v)
  --     (cudaDescendUnsafeAdamGPU len step alpha beta1 beta2 epsilon lambda weights gradient m v)
  | otherwise = descendUnsafeAdam len step alpha beta1 beta2 epsilon lambda weights gradient m v

foreign import ccall unsafe
    descend_sgd_cpu
      :: Int -> RealNum -> RealNum -> RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> IO ()

foreign import ccall unsafe
    descend_adam_cpu
      :: Int -> Int -> RealNum -> RealNum -> RealNum -> RealNum -> RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> IO ()

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE Strict #-}
{-# LANGUAGE FlexibleContexts #-}

module CPVO.IO.Reader.Ecalj.Common where
import CPVO.Numeric
import CPVO.IO

import qualified Data.Text as T
import qualified Data.Text.IO as T
import Data.List.Split
import Data.List
import Data.Maybe
import Numeric.LinearAlgebra
import Data.Char (ord)
import Control.Monad.IO.Class (MonadIO)

readCtrlAtoms :: String -> String -> IO [T.Text]
readCtrlAtoms tailer foldernya = do
    fCtrl <- T.readFile $ foldernya ++ "/ctrl." ++ tailer
    return $ catMaybes
      $ map ( T.stripPrefix "ATOM=" .  head)
      $ filter (/=[])
      $ map ( T.words . T.takeWhile (/='#') )
      $ head
      $ splitWhen (T.isPrefixOf "SPEC")
      $ last $ splitWhen (T.isPrefixOf "SITE")
      $ T.lines fCtrl

    -- uniqAtoms : [(count,noFirstAtom,atomicSymbol)]
type CntN1Atsym = (Int,Int,T.Text)
readUniqAtoms :: [T.Text] -> [CntN1Atsym]
readUniqAtoms allAtoms =
          map (\a -> (length a, snd $ head a, fst $ head a)) $
          groupBy (\a b -> fst a == fst b) $
          zip  allAtoms ([1..] :: [Int])

-- daftarCetak :: [((notCtk,((count,noFrstAtom,atomSymbol),label  ,[    AOs ])),spin)]
-- daftarCetak :: [(( 1    ,((  2  ,    13    ,  "Ni"    ),"Ni#2p",["3","4","5"])),1)]
genDaftarCetak :: [CntN1Atsym] -> String -> String -> [String]
               -> IO [((Integer,(CntN1Atsym,String,[String])),Integer)]
genDaftarCetak listAtoms _ _ aos = do
  return $ [ (i,j) | i <- daftarCetak' listAtoms , j <- [1,2] ]
  where
    daftarCetak' lAtoms = zip [1..]
                      $ map (\(a,label,b) -> (head $ filter (\(_,_,aa) -> aa == (T.pack a)) lAtoms , label, b) )
                      $ map ( (\(a:label:as) -> (a,label,as)) . splitOn ":") aos

-- ctrlAtomicAOs :: [(atomNumber,(atomSym,(label ,[intAOs])))]
-- ctrlAtomicAOs :: [(    1     ,(  "O"  ,("O#2p",[ 2,3,4])))]
-- ctrlAtomicAOs :: 14 atoms
genCtrlAtomicAOs :: [(String, String, [Int])]
                 -> [T.Text]
                 -> [(Int, (String, (String, [Int])))]
genCtrlAtomicAOs aoSet ctrlAtoms =  map (\x -> (head $ takeAOs x aoSet))
          $ concat
          $ groupBy (\(_,a:_) (_,b:_) -> (ord a) == (ord b))
          $ zip ([1..]::[Int]) $ map T.unpack ctrlAtoms

-- totalDOS :: Matrix Double [ energy, DOSspinUp, DOSspinDown ]
readTotalDOSText :: String -> String -> IO (Matrix Double)
readTotalDOSText tailer foldernya = loadMatrix $ foldernya ++ "/dos.tot." ++ tailer

-----------------------------------------------------------
getLastLLMF :: MonadIO io =>
  String -> io [T.Text]
getLastLLMF foldernya = inshell2text $ concat ["ls -laht ", foldernya,"/llmf{,_gwscend.*} | head -1|awk '{print $NF}'" ]
-----------------------------------------------------------

--readHeaderData (texFile:jd:jdHead:colAlign:xr:ymax':wTot::invS:tailer:foldernya:aos) = do
readHeaderData :: [String]
               -> IO ( Either String ( Double, Double, Double, Double, [T.Text]
                      , [CntN1Atsym]
                      , [(Int, (String, (String,[Int])))]
                      , String, [String], String, String, String, String))
readHeaderData (texFile:jd:jdHead:colAlign:xr:ymax':_:_:invS:tailer:foldernya:aos) = do
  -------------------------------reading data------------------------
    let invStat = if (invS == "flipSpin") then (-1) else 1
    let ymax = read ymax' :: Double
    let [xmin,xmax] = map (read :: String -> Double) $ splitOn ":" xr
    ctrlAtoms <- readCtrlAtoms tailer foldernya
    let jdHeads = splitOn "|" jdHead
    let uniqAtoms = readUniqAtoms ctrlAtoms
    putStrLn $ show ctrlAtoms
    putStrLn $ show uniqAtoms
    -- daftarCetak : [(nourut,,jumlah,nourut,symbol)]
    daftarCetak <- genDaftarCetak uniqAtoms tailer foldernya aos
    putStrLn $ show aos
    putStrLn $ show $ last daftarCetak
      -------------------------------generating DOS data------------------------
    totalDOS <- readTotalDOSText tailer foldernya
      -------------------------------integrating DOS data------------------------
    let intgTot = map (\i -> integrateToZero $ totalDOS ¿ [0,i]) [1,2] -- run it on spin [1,2]
    putStrLn $ show intgTot
    let aoSet = map ( (\(n:l:as) -> (n,l,map ( ((+) (-1)) . read :: String -> Int) as) ) . splitOn ":") aos
    {-
      -------------------------------generating PDOS data------------------------
              -- map ditambah -1 karena input mengikuti gnuplot
              -- input : d kolom 6-10
              -- gnuplot : d kolom 6-10
              -- hmatrix : d kolom 5-9

              -- ((namaAtom,jdAtom,[intAOs]),[(nourutAtom,namaAtom)])
              -- ((String , String,[ Int  ]),[(Int       , String )])
              -- (("O"    ,"O#2p" ,[2,3,4 ]),[(1         ,"O"     )])
      -}
    let ctrlAtomicAOs = genCtrlAtomicAOs aoSet ctrlAtoms
    let jdTable = "Table: " ++ jd
    putStrLn $ show $ head ctrlAtomicAOs
--    pdosAtomicPilihan <- readPDOS invStat tailer foldernya $ take 2 ctrlAtomicAOs
--    let integratedAtomicPDOS = integrateAtomicPDOS pdosAtomicPilihan
--    putStrLn $ show $ integratedAtomicPDOS
    putStrLn "===done:readHeaderData@CPVO/IO/Reader/Common ====================="
    return $ Right
      (invStat, ymax, xmin, xmax, ctrlAtoms, uniqAtoms, ctrlAtomicAOs,jdTable, jdHeads, foldernya, tailer, colAlign, texFile)

readHeaderData _ = return $ Left
  "===Error:readHeaderData@CPVO/IO/Reader/Common wrong args ========"

max_stars_repo_name: hasanalrasyid/cpvoh

max_issues_repo_name: hasanalrasyid/cpvoh

max_forks_repo_name: hasanalrasyid/cpvoh


text:
module Fuml.Base.Logistic where

import Numeric.LinearAlgebra
import qualified Data.Vector.Storable as VS
import Fuml.Optimisation.BFGS
import Data.List (foldl1')

-- maybe follow this: https://idontgetoutmuch.wordpress.com/2013/04/30/logistic-regression-and-automated-differentiation-3/

logit :: Floating a => a -> a
logit x = 1 / (1 + exp (negate x))

logLikelihood1 ::  Vector Double -> Vector Double -> Bool -> Double
logLikelihood1 theta x y = ind y * log (logit z) +
                           (1 - ind y) * log (1 - logit z)
  where
    z = VS.sum $ VS.zipWith (*) theta x

-- negative log likelihood
logLikelihood :: Double -> [(Vector Double, Bool)] -> Vector Double -> Double
logLikelihood delta theData theta = negate $ (a - delta*b)/l where
  l = fromIntegral $ length theData
  a = sum $ map (uncurry $ logLikelihood1 theta) theData
  b = (/2) $ VS.sum $ VS.map (^2) theta

--https://www.cs.ox.ac.uk/people/nando.defreitas/machinelearning/lecture6.pdf
gradLogLikelihood1 :: Vector Double -> Vector Double -> Bool -> Vector Double
gradLogLikelihood1 theta x y = VS.map (*(ind y - lz)) x
--                           (1 - ind y) * log (1 - logit z)
  where
    z = VS.sum $ VS.zipWith (*) theta x
    lz = logit z

gradLogLikelihood :: Double -> [(Vector Double, Bool)] -> Vector Double -> Vector Double
gradLogLikelihood delta theData theta =
  let vs = map (uncurry $ gradLogLikelihood1 theta) theData
      vsum = foldl1' vadd vs
      l = fromIntegral $ length theData
  in VS.map (negate . (/l)) vsum  `vadd` VS.map ((/l) . (*delta)) theta

logisticRegression :: [(Vector Double, Bool)] -> Vector Double
logisticRegression theData =
  let start = VS.map (const 0) $ fst $ head theData
      inisbox = VS.map (const (0.1:: Double)) $ fst $ head theData
      hessInit = ident $ size start
  in --fst $ minimizeV NMSimplex 1e-4 200 inisbox (logLikelihood theData) start

     --fst $ minimizeVD ConjugateFR 1e-10 500 0.01 0.01 (logLikelihood theData) (gradLogLikelihood theData) start -- inisbox (logLikelihood theData) start

     case bfgsWith myBOpts (logLikelihood 0 theData) (gradLogLikelihood 0 theData) start hessInit of
       Left s -> error s
       Right (p, h) -> p
      --fst $ minimizeV NMSimplex 1e-3 100 inisbox (logLikelihood theData) start

myBOpts = BFGSOpts 1e-7 2e-5 200

ind :: Bool -> Double
ind True = 1
ind False = 0

vadd :: Vector Double -> Vector Double -> Vector Double
vadd = VS.zipWith (+)

max_stars_repo_name: ekalosak/open

max_issues_repo_name: ekalosak/open

max_forks_repo_name: ekalosak/open


text:
module Lib (
    readExpr
  , eval
  , valueToPrint
) where

import Control.Monad
import Control.Monad.Error
import Control.Monad.Error.Class
import Data.Array
import Data.Complex
import Data.Ratio
import Numeric

import Text.ParserCombinators.Parsec hiding (spaces)

data LispVal = Atom String
             | Vector (Array Int LispVal)
             | List [LispVal]
             | DottedList [LispVal] LispVal
             | Number Integer
             | Real Double
             | Rational Rational
             | Complex (Complex Double)
             | String String
             | Character Char
             | Bool Bool

data LispError = NumArgs Integer [LispVal]
               | TypeMismatch String LispVal
               | Parser ParseError
               | BadSpecialForm String LispVal
               | NotFunction String String
               | UnboundVar String String
               | Default String

type ThrowsError = Either LispError

showVals :: [LispVal] -> String
showVals = unwords . map show

instance Show LispVal where
    show (Atom name) = name
    show (Vector contents) = "(" ++ (showVals $ elems contents) ++ ")"
    show (List contents) = "(" ++ showVals contents ++ ")"
    show (DottedList init last) = "(" ++ showVals init ++ " . " ++ show last ++ ")"
    show (Number n) = show n
    show (Real x) = show x
    show (Rational r) = show r
    show (Complex c) = show c
    show (String content) = content
    show (Character c) = charToString c
    show (Bool True) = "#t"
    show (Bool False) = "#f"

instance Show LispError where
    show (UnboundVar message varname) = message ++ ": " ++ varname
    show (BadSpecialForm message form) = message ++ ": " ++ show form
    show (NotFunction message func) = message ++ ": " ++ show func
    show (NumArgs expected found) = "Expected " ++ show expected ++ " arguments. Found values " ++ showVals found
    show (TypeMismatch expected found) = "Invalid type: expected " ++ expected ++ ", found " ++ show found
    show (Parser parseErr) = "Parse error at " ++ show parseErr

instance Error LispError where
     noMsg = Default "An error has occurred"
     strMsg = Default

symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

parseAtom :: Parser LispVal
parseAtom = do
    first <- letter <|> symbol
    rest <- many (letter <|> digit <|> symbol)
    let atom = first:rest
    return $ Atom atom

spaces :: Parser ()
spaces = skipMany1 space

parseVector' :: Parser LispVal
parseVector' = do
    vectorValues <- sepBy parseExpr spaces
    return $ Vector (listArray (0, (length vectorValues - 1)) vectorValues)

parseVector :: Parser LispVal
parseVector = do
    string "#("
    x <- parseVector'
    char ')'
    return x

parseDatum :: Parser LispVal
parseDatum = do
    char '.'
    spaces
    parseExpr

parseListContent :: Parser LispVal
parseListContent = do
    list <- sepEndBy parseExpr spaces
    datum <- optionMaybe parseDatum
    return $ case datum of
        Nothing -> List list
        Just datum -> DottedList list datum

parseList :: Parser LispVal
parseList = do
    char '('
    list <- parseListContent
    char ')'
    return list

parseQuoted :: Parser LispVal
parseQuoted = do
    char '\''
    x <- parseExpr
    return $ List [Atom "quote", x]

parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
    char '`'
    x <- parseExpr
    return $ List [Atom "quasiquote", x]

parseUnquoted :: Parser LispVal
parseUnquoted = do
    char ','
    x <- parseExpr
    return $ List [Atom "unquote", x]

parseDecimal :: Parser LispVal
parseDecimal = do
    x <- many1 digit
    (return . Number . read) x

parseDec :: Parser LispVal
parseDec = do
    try $ string "#d"
    x <- many1 digit
    (return . Number . read) x

hex2dig x = fst . head $ readHex x

parseHex :: Parser LispVal
parseHex = do
    try $ string "#h"
    x <- many1 hexDigit
    (return . Number . hex2dig) x

oct2dig x = fst . head $ readOct x

parseOct :: Parser LispVal
parseOct = do
    try $ string "#o"
    x <- many1 octDigit
    (return . Number . oct2dig) x

bin2dig = bin2dig' 0
bin2dig' value "" = value
bin2dig' value (x:xs) = let newValue = 2 * value + (if x == '0' then 0 else 1) in
                        bin2dig' newValue xs

parseBin :: Parser LispVal
parseBin = do
    try $ string "#b"
    let binDigit = oneOf "10"
    x <- many1 binDigit
    (return . Number . bin2dig) x

parseNumber :: Parser LispVal
parseNumber = parseDecimal
          <|> parseDec
          <|> parseHex
          <|> parseOct
          <|> parseBin

parseReal :: Parser LispVal
parseReal = do
    x <- many1 digit
    char '.'
    y <- many1 digit
    (return . Real . fst . head . readFloat) (y ++ "." ++ x)

parseRational :: Parser LispVal
parseRational = do
    x <- many1 digit
    char '/'
    y <- many1 digit
    (return . Rational) ((read x) % (read y))

toDouble :: LispVal -> Double
toDouble (Real r) = realToFrac r
toDouble (Number n) = fromIntegral n

parseComplex :: Parser LispVal
parseComplex = do
    x <- try parseReal <|> parseDecimal
    char '+'
    y <- try parseReal <|> parseDecimal
    char 'i'
    (return . Complex) (toDouble x :+ toDouble y)

escapedChars :: Parser Char
escapedChars = do
    char '\\'
    x <- oneOf "\\\"nrt"
    return $ case x of
        '\\' -> x
        '"'  -> x
        'n'  -> '\n'
        'r'  -> '\r'
        't'  -> '\t'

parseString :: Parser LispVal
parseString = do
    char '"'
    x <- many $ try escapedChars <|> noneOf "\""
    char '"'
    return $ String x

charToString :: Char -> String
charToString x = [x]

parseCharacter :: Parser LispVal
parseCharacter = do
    string "#\\"
    value <- try (string "newline" <|> string "space")
        <|> do
            x <- anyChar
            notFollowedBy alphaNum
            return $ charToString x
    return $ Character $ case value of
        "space" -> ' '
        "newline" -> '\n'
        otherwise -> value !! 0

parseBool :: Parser LispVal
parseBool = do
    char '#'
    val <- oneOf "ft"
    return $ case val of
        'f' -> Bool False
        't' -> Bool True

parseExpr :: Parser LispVal
parseExpr = try parseAtom
        <|> try parseVector
        <|> try parseList
        <|> try parseQuoted
        <|> try parseQuasiQuoted
        <|> try parseUnquoted
        <|> try parseReal
        <|> try parseRational
        <|> try parseComplex
        <|> try parseNumber
        <|> try parseString
        <|> try parseCharacter
        <|> try parseBool

primitives :: [(String, [LispVal] -> ThrowsError LispVal)]
primitives = [("+", numberOp (+)),
              ("-", numberOp (-)),
              ("*", numberOp (*)),
              ("/", numberOp div),
              ("mod", numberOp mod),
              ("quotient", numberOp quot),
              ("remainder", numberOp rem),
              ("=", numberOpBool (==)),
              ("<", numberOpBool (<)),
              (">", numberOpBool (>)),
              ("/=", numberOpBool (/=)),
              ("<=", numberOpBool (<=)),
              (">=", numberOpBool (>=)),
              ("&&", boolOpBool (&&)),
              ("||", boolOpBool (||)),
              ("string=?", stringOpBool (==)),
              ("string<?", stringOpBool (<)),
              ("string>?", stringOpBool (>)),
              ("string<=?", stringOpBool (<=)),
              ("string>=?", stringOpBool (>=)),
              ("symbol?", isAtom),
              ("vector?", isVector),
              ("list?", isList),
              ("number?", isNumber),
              ("real?", isReal),
              ("rational?", isRational),
              ("complex?", isComplex),
              ("string?", isString),
              ("character?", isCharacter),
              ("bool?", isBool),
              ("string->symbol", string2symbol),
              ("symbol->string", symbol2string)]

unpackNumber :: LispVal -> ThrowsError Integer
unpackNumber (Number n) = return n
unpackNumber notNumber = throwError $ TypeMismatch "number" notNumber

unpackString :: LispVal -> ThrowsError String
unpackString (String n) = return n
unpackString notString = throwError $ TypeMismatch "string" notString

unpackBool :: LispVal -> ThrowsError Bool
unpackBool (Bool n) = return n
unpackBool notBool = throwError $ TypeMismatch "bool" notBool

getNumberValues :: [LispVal] -> ThrowsError [Integer]
getNumberValues lispvals = mapM unpackNumber lispvals

numberOp :: (Integer -> Integer -> Integer) -> [LispVal] -> ThrowsError LispVal
numberOp op [] = throwError $ NumArgs 2 []
numberOp op arg@[x] = throwError $ NumArgs 2 arg
numberOp op args = getNumberValues args >>= return . Number . foldl1 op

opBool :: (LispVal -> ThrowsError a) -> (a -> a -> Bool) -> [LispVal] -> ThrowsError LispVal
opBool unpacker op args = if length args /= 2
                          then throwError $ NumArgs 2 args
                          else do left <- unpacker $ head args
                                  right <- unpacker $ last args
                                  return $ Bool $ left `op` right

numberOpBool = opBool unpackNumber
stringOpBool = opBool unpackString
boolOpBool = opBool unpackBool

isAtom :: [LispVal] -> ThrowsError LispVal
isAtom [Atom _] = return $ Bool True
isAtom [List (Atom "quoted":xs)] = return $ Bool True
isAtom [_] = return $ Bool False
isAtom x = throwError $ NumArgs 1 x

isVector :: [LispVal] -> ThrowsError LispVal
isVector [Vector _] = return $ Bool True
isVector [_] =return $  Bool False
isVector x = throwError $ NumArgs 1 x

isList :: [LispVal] -> ThrowsError LispVal
isList [List _] = return $ Bool True
isList [DottedList _ _] = return $ Bool True
isList [_] = return $ Bool False
isList x = throwError $ NumArgs 1 x

isNumber :: [LispVal] -> ThrowsError LispVal
isNumber [Number _] = return $ Bool True
isNumber [_] = return $ Bool False
isNumber x = throwError $ NumArgs 1 x

isReal :: [LispVal] -> ThrowsError LispVal
isReal [Number _] = return $ Bool True
isReal [_] = return $ Bool False
isReal x = throwError $ NumArgs 1 x

isRational :: [LispVal] -> ThrowsError LispVal
isRational [Number _] = return $ Bool True
isRational [_] = return $ Bool False
isRational x = throwError $ NumArgs 1 x

isComplex :: [LispVal] -> ThrowsError LispVal
isComplex [Number _] = return $ Bool True
isComplex [_] = return $ Bool False
isComplex x = throwError $ NumArgs 1 x

isString :: [LispVal] -> ThrowsError LispVal
isString [String _] = return $ Bool True
isString [_] = return $ Bool False
isString x = throwError $ NumArgs 1 x

isCharacter :: [LispVal] -> ThrowsError LispVal
isCharacter [Character _] = return $ Bool True
isCharacter [_] = return $ Bool False
isCharacter x = throwError $ NumArgs 1 x

isBool :: [LispVal] -> ThrowsError LispVal
isBool [Bool _] = return $ Bool True
isBool [_] = return $ Bool False
isBool x = throwError $ NumArgs 1 x

string2symbol :: [LispVal] -> ThrowsError LispVal
string2symbol [String s] = return $ Atom s
string2symbol [x] = throwError $ TypeMismatch "string" x
string2symbol x = throwError $ NumArgs 1 x

symbol2string :: [LispVal] -> ThrowsError LispVal
symbol2string [Atom s] = return $ String s
symbol2string [x] = throwError $ TypeMismatch "symbol" x
symbol2string x = throwError $ NumArgs 1 x

valueToPrint = extractValue . trapError

trapError action = catchError action (return . show)

extractValue :: ThrowsError a -> a
extractValue (Right val) = val

apply :: String -> [LispVal] -> ThrowsError LispVal
apply func args = maybe (throwError $ NotFunction "Unrecognized primitive function args" func)
                        ($ args)
                        (lookup func primitives)

eval :: LispVal -> ThrowsError LispVal
eval atom@(Atom _) = return atom
eval vector@(Vector _) = return vector
eval list@(List [Atom "quote", value]) = return value
eval list@(List (Atom func : args)) = mapM eval args >>= apply func
eval list@(List _) = return list
eval dotted@(DottedList _ _) = return dotted
eval number@(Number _) = return number
eval real@(Real _) = return real
eval rational@(Rational _) = return rational
eval complex@(Complex _) = return complex
eval string@(String _) = return string
eval char@(Character _) = return char
eval bool@(Bool _) = return bool

readExpr :: String -> ThrowsError LispVal
readExpr input = case parse parseExpr "lisp" input of
    Left err -> throwError $ Parser err
    Right value -> return value

max_stars_repo_name: prcastro/haskell-scheme

max_issues_repo_name: prcastro/haskell-scheme

max_forks_repo_name: prcastro/haskell-scheme


text:
{-# LANGUAGE DeriveAnyClass    #-}
{-# LANGUAGE DeriveGeneric     #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE Strict            #-}
-- Type.hs ---
--
-- Filename: Type.hs
-- Description:
-- Author: Manuel Schneckenreither
-- Maintainer:
-- Created: Wed Nov  1 17:40:44 2017 (+0100)
-- Version:
-- Package-Requires: ()
-- Last-Updated:
--           By:
--     Update #: 94
-- URL:
-- Doc URL:
-- Keywords:
-- Compatibility:
--
--

-- Commentary:
--
--
--
--

-- Change Log:
--
--
--
--
--
--
--

-- Code:

module SimSim.Order.Type where

import           ClassyPrelude
import           Control.DeepSeq
import           Data.Serialize
import           GHC.Generics
import           Statistics.Distribution
import           System.Random.MWC

import           SimSim.Block
import           SimSim.ProductType
import           SimSim.Time

type OrderId = Integer

-- | The message type is a record of user name, content and timestamp. Equality is defined by the orderId, that means
-- right after creating an order with @newOrder@ all orders are equal until fed into the simulation.
data Order = Order
  { orderId          :: !OrderId
  , productType      :: !ProductType
  , arrivalDate      :: !ArrivalDate
  , dueDate          :: !DueDate
  , released         :: !(Maybe Time)
  , prodStart        :: !(Maybe Time)
  , prodEnd          :: !(Maybe Time)
  , shipped          :: !(Maybe Time)
  , lastBlock        :: !Block
  , blockStartTime   :: !Time
  , nextBlock        :: !Block
  , orderCurrentTime :: !Time
  } deriving (Ord, Show, Generic, Serialize)

instance NFData Order where
  rnf (Order i tp ar dd rel st end sh last blSt nBl t) =
    rnf i `seq` rnf tp `seq` rnf ar `seq` rnf dd `seq` rnf1 rel `seq` rnf1 st `seq` rnf1 end `seq` rnf1 sh `seq` rnf last `seq` rnf blSt `seq` rnf nBl `seq` rnf t

instance Eq Order where
  x == y = orderId x == orderId y

-- | This function creates a new order.
newOrder :: ProductType -> ArrivalDate -> DueDate -> Order
newOrder productType arrDate dueDate =
  Order (-1) productType arrDate dueDate Nothing Nothing Nothing Nothing OrderPool arrDate OrderPool 0

orderSlackTime :: Order -> Time
orderSlackTime order = dueDate order - arrivalDate order

orderFgiTime :: Order -> Maybe Time
orderFgiTime order = (-) <$> shipped order <*> prodEnd order

orderBoTime :: Order -> Maybe Time
orderBoTime order = max 0 <$> ((-) <$> shipped order <*> Just (dueDate order))

orderWipTime :: Order -> Maybe Time
orderWipTime order = (-) <$> prodEnd order <*> prodStart order

orderFinishedProduction :: Order -> Bool
orderFinishedProduction order = isJust (prodEnd order)

orderFinishedAndTardyProduction :: Order -> Bool
orderFinishedAndTardyProduction order = maybe False (> dueDate order) (prodEnd order)

orderFinishedAndTardyShipped :: Order -> Bool
orderFinishedAndTardyShipped order = maybe False (> dueDate order) (shipped order)

orderTardinessProduction :: Order -> Maybe Time
orderTardinessProduction order = packTardy . subtract (dueDate order) =<< prodEnd order
  where packTardy x | x > 0 = return x
                    | otherwise = Nothing -- "not tardy"

orderTardinessShipped :: Order -> Maybe Time
orderTardinessShipped order = packTardy . subtract (dueDate order) =<< shipped order
  where packTardy x | x > 0 = return x
                    | otherwise = Nothing -- "not shipped"

setOrderId :: OrderId -> Order -> Order
setOrderId nr o = o {orderId = nr}

setOrderCurrentTime :: Time -> Order -> Order
setOrderCurrentTime t o = o {orderCurrentTime = t}

setOrderBlockStartTime :: Time -> Order -> Order
setOrderBlockStartTime t o = o { blockStartTime = t}

addToOrderCurrentTime :: Time -> Order -> Order
addToOrderCurrentTime t o = o {orderCurrentTime = orderCurrentTime o + t}

setReleaseTime :: Time -> Order -> Order
setReleaseTime t o = o {released = Just t}

setProdStartTime :: Time -> Order -> Order
setProdStartTime t o =
  if isNothing (prodStart o)
    then o {prodStart = Just t}
    else o

setProdEndTime :: Time -> Order -> Order
setProdEndTime t o =
  if isNothing (prodEnd o)
    then o {prodEnd = Just t}
    else o

setLastBlock :: Block -> Order -> Order
setLastBlock b o = o {lastBlock = b}

setNextBlock :: Block -> Order -> Order
setNextBlock b o = o {nextBlock = b}

setShippedTime :: Time -> Order -> Order
setShippedTime t o = o {shipped = Just t}

--
-- Type.hs ends here

max_stars_repo_name: schnecki/simsim

max_issues_repo_name: schnecki/simsim

max_forks_repo_name: schnecki/simsim


text:
{-# LANGUAGE ForeignFunctionInterface #-}
module Grenade.Layers.Internal.Update (
    descendMatrix
  , descendVector
  ) where

import           Data.Maybe ( fromJust )
import qualified Data.Vector.Storable as U ( unsafeToForeignPtr0, unsafeFromForeignPtr0 )

import           Foreign ( mallocForeignPtrArray, withForeignPtr )
import           Foreign.Ptr ( Ptr )
import           GHC.TypeLits

import           Numeric.LinearAlgebra ( Vector, flatten )
import           Numeric.LinearAlgebra.Static
import qualified Numeric.LinearAlgebra.Devel as U

import           System.IO.Unsafe ( unsafePerformIO )

descendMatrix :: (KnownNat rows, KnownNat columns) => Double -> Double -> Double -> L rows columns -> L rows columns -> L rows columns -> (L rows columns, L rows columns)
descendMatrix rate momentum regulariser weights gradient lastUpdate =
  let (rows, cols) = size weights
      len          = rows * cols
      -- Most gradients come in in ColumnMajor,
      -- so we'll transpose here before flattening them
      -- into a vector to prevent a copy.
      --
      -- This gives ~15% speed improvement for LSTMs.
      weights'     = flatten . tr . extract $ weights
      gradient'    = flatten . tr . extract $ gradient
      lastUpdate'  = flatten . tr . extract $ lastUpdate
      (vw, vm)     = descendUnsafe len rate momentum regulariser weights' gradient' lastUpdate'

      -- Note that it's ColumnMajor, as we did a transpose before
      -- using the internal vectors.
      mw           = U.matrixFromVector U.ColumnMajor rows cols vw
      mm           = U.matrixFromVector U.ColumnMajor rows cols vm
  in  (fromJust . create $ mw, fromJust . create $ mm)

descendVector :: (KnownNat r) => Double -> Double -> Double -> R r -> R r -> R r -> (R r, R r)
descendVector rate momentum regulariser weights gradient lastUpdate =
  let len          = size weights
      weights'     = extract weights
      gradient'    = extract gradient
      lastUpdate'  = extract lastUpdate
      (vw, vm)     = descendUnsafe len rate momentum regulariser weights' gradient' lastUpdate'
  in  (fromJust $ create vw, fromJust $ create vm)

descendUnsafe :: Int -> Double -> Double -> Double -> Vector Double -> Vector Double -> Vector Double -> (Vector Double, Vector Double)
descendUnsafe len rate momentum regulariser weights gradient lastUpdate =
  unsafePerformIO $ do
    outWPtr <- mallocForeignPtrArray len
    outMPtr <- mallocForeignPtrArray len
    let (wPtr, _) = U.unsafeToForeignPtr0 weights
    let (gPtr, _) = U.unsafeToForeignPtr0 gradient
    let (lPtr, _) = U.unsafeToForeignPtr0 lastUpdate

    withForeignPtr wPtr $ \wPtr' ->
      withForeignPtr gPtr $ \gPtr' ->
        withForeignPtr lPtr $ \lPtr' ->
          withForeignPtr outWPtr $ \outWPtr' ->
            withForeignPtr outMPtr $ \outMPtr' ->
              descend_cpu len rate momentum regulariser wPtr' gPtr' lPtr' outWPtr' outMPtr'

    return (U.unsafeFromForeignPtr0 outWPtr len, U.unsafeFromForeignPtr0 outMPtr len)

foreign import ccall unsafe
    descend_cpu
      :: Int -> Double -> Double -> Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> IO ()

max_issues_repo_name: Alien-Inc/grenade

max_forks_repo_name: Alien-Inc/grenade


text:
{-# LANGUAGE Safe #-}
-- | Import all modules from prelude that should be safe
module Main where

import Numeric
import Prelude
-- import Foreign

import Control.Applicative
import Control.Arrow
import Control.Category

-- import Control.ST
-- import Control.ST.Lazy
-- import Control.ST.Strict

import Control.Concurrent
import Control.Concurrent.Chan
import Control.Concurrent.MVar

import Control.Exception
import Control.Exception.Base

import Control.Monad
import Control.Monad.Fix
import Control.Monad.Zip

import Data.Bits
import Data.Bool
import Data.Char
import Data.Complex
import Data.Either
import Data.Eq
import Data.Fixed
import Data.Foldable
import Data.Function
import Data.Functor
import Data.IORef
import Data.Int
import Data.Ix
import Data.List
import Data.Maybe
import Data.Monoid
import Data.Ord
import Data.Ratio
import Data.String
import Data.Traversable
import Data.Tuple
import Data.Typeable
import Data.Unique
import Data.Version
import Data.Word

import Data.STRef
import Data.STRef.Lazy
import Data.STRef.Strict

-- import Debug.Trace

import Foreign.Concurrent
-- import Foreign.ForeignPtr
import Foreign.Ptr
import Foreign.StablePtr
import Foreign.Storable

import Foreign.C
import Foreign.C.Error
import Foreign.C.String
import Foreign.C.Types

-- import Foreign.Marshal
import Foreign.Marshal.Alloc
import Foreign.Marshal.Array
import Foreign.Marshal.Error
import Foreign.Marshal.Pool
import Foreign.Marshal.Utils

import System.CPUTime
import System.Environment
import System.Exit
import System.Info
import System.Mem
import System.Timeout

import System.Console.GetOpt

import System.IO
import System.IO.Error

import System.Mem.StableName
import System.Mem.Weak

import System.Posix.Internals
import System.Posix.Types

import Text.Printf

import Text.Read
import Text.Read.Lex

import Text.Show
import Text.Show.Functions

import Text.ParserCombinators.ReadP
import Text.ParserCombinators.ReadPrec

-- import Unsafe.Coerce

f :: Int
f = 2

main :: IO ()
main = putStrLn $ "X is: " ++ show f

max_stars_repo_name: ocharles/ghc

max_issues_repo_name: ocharles/ghc

max_forks_repo_name: ocharles/ghc


text:
{-# LANGUAGE FlexibleContexts,FlexibleInstances,IncoherentInstances #-}
{- Assignment 1 Extra Credit
 - Name: TODO add full name
 - Date: TODO add of completion
 -}
module Assign_2_ExtraCredit where

import Data.Complex

macid = "TODO: put your mac id here"

data GaussianInt a = a :@ a
  deriving Show

class GaussianIntegral g where
  gaussZero :: Integral a => g a
  gaussReal :: Integral a => g a -> a
  gaussImag :: Integral a => g a -> a
  gaussConj :: Integral a => g a -> g a
  gaussAdd :: Integral a => g a -> g a -> g a
  gaussMult :: Integral a => g a -> g a -> g a

{- TODO
 -   implement instances of GaussianIntegral
 -   implmenet instances of Eq, Ord
 -}

gaussNorm :: (Integral a, GaussianIntegral g) => g a -> a
gaussNorm g = error "TODO: implement gaussNorm"

maxGaussNorm :: (Integral a, GaussianIntegral g) => [g a] -> g a
maxGaussNorm gs = error "TODO: implement maxGaussNorm"

max_stars_repo_name: SaribK/Gaussian-Integers

max_issues_repo_name: SaribK/Gaussian-Integers

max_forks_repo_name: SaribK/Gaussian-Integers


text:
{-# language OverloadedStrings #-}
{-# language TypeApplications #-}

import Geometry
import Data.Complex

iterateN g f = zipWith g [1..] . iterate f

main = do
  let romashka =
        group [ aSegment # rotate (18 * fromIntegral i) #: label (show i) <> white
              | i <- [1..20] ]
      clock = aCircle <+> group (modularScale 24 (aCircle # rotate 90 ))
      hand h = mkTriangle @CN [0.03:+0, 0:+0.9, (-0.03):+0] #: fill "red" # rotate (-15 * h)
               <+> aCircle # scale 0.05 #: fill "red"

  writeSVG 450 "figures/1.5/26.svg" $ let
    seg = aSegment # scale (1/8) #: "1" <+> point (1/8,0)
    segs = aPoint <+> (group $ take 8 $ iterate (translate (1/8, 0)) seg)
    rays = aSegment <+>
           origin #: "O" <+>
           (group $ take 7 $ tail $ iterate (rotate 8) (aSegment #: white)) <+>
           aSegment # rotate 64 #: "1"
    in rays <+> segs # translate (0,-0.25)

  writeSVG 500 "figures/1.5/28.svg" $
    group [ aSegment # rotate (18 * fromIntegral i) #: label (show i)
          | i <- [1..20] ]

  writeSVG 500 "figures/1.5/30.svg" $
    (group [ aSegment # rotate (18 * fromIntegral i) #: label (show i) <> white
           | i <- [1..6] ]) <+>
    aSegment # rotate (18 * 7) #: "7" <+>
    aSegment <+>
    aCircle #: invisible

  writeSVG 500 "figures/1.5/31.svg" $
    (group [ aSegment # rotate (18 * fromIntegral i) #: label (show i) <> white
           | i <- [1..6] ]) <+>
    aSegment # rotate (18 * 7) #: "7" <+>
    (group [ aSegment # rotate (18 * fromIntegral i) #: label (show (i - 7)) <> white
           | i <- [8..19] ]) <+>
    aSegment # rotate (18 * 20) #: "13" <+>
    aCircle #: invisible

  writeSVG 500 "figures/1.5/33.svg" $
    romashka <+> aSegment # rotate (18 * 13)

  writeSVG 500 "figures/1.5/34.svg" $
    romashka <+> aSegment # rotate (18 * 3)

  writeSVG 300 "figures/1.5/35.svg" $
    clock <+> hand 11

  writeSVG 300 "figures/1.5/37.svg" $
    clock <+> hand 17

  writeSVG 200 "figures/1.5/38.svg" $
    let dow l i = aLabel
                  # on aCircle (-i * 1/7) # rotate 90
                  #: label l
    in group $ zipWith dow ["Пн","Вт","Ср","Чт","Пт","Сб","Вс"] [1..]

  writeSVG 300 "figures/1.5/51.svg" $
    let dec = regularPoly 10
        sc = group $ modularScale 10 dec
    in dec <+> sc

  writeSVG 450 "figures/1.5/53.svg" $
    group [ aSegment # scale 0.9 # rotate i <+>
            aLabel # at' i #: label (show i)
          | i <- asDeg <$> [10,20..360] ]

  writeSVG 450 "figures/1.5/54.svg" $
    anAngle 180 #: "#" <+> aPoint <+>
    aSegment <+> aSegment # rotate 180

  writeSVG 450 "figures/1.5/55.svg" $
    anAngle 90 #: "#" <+> aPoint <+>
    aSegment <+> aSegment # rotate 180 <+>
    aSegment # rotate 90

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
module RayLib where

import Numeric.LinearAlgebra
import qualified Graphics.Image as G
import qualified Debug.Trace as Debug

data Ray = Ray {rayOrigin :: Vector Double,
                rayDirection :: Vector Double}
    deriving (Show, Eq)

eps :: Double
eps = 0.000001

data Isect = Isect {iT :: Double,
                    iNormal :: Vector Double,
                    iObject :: Geometry,
                    iMaterial :: Material}

instance Eq Isect where
    (Isect iT1 iNormal1 _ _) == (Isect iT2 iNormal2 _ _) = (iT1 == iT2) && (iNormal1 == iNormal2)

instance Ord Isect where
    compare (Isect iT1 _ _ _) (Isect iT2 _ _ _) = compare iT1 iT2

instance Show Isect where
    show (Isect t n _ _) = "time: " ++ show t ++ ", normal: " ++ show n

class At a where
    at :: a -> Ray -> Vector Double

instance At Double where
    at t (Ray p d) = p + (fromList [t] :: Vector Double) * d

instance At Isect where
    at (Isect t _ _ _) (Ray p d) = p + (fromList [t] :: Vector Double) * d

type Color = Vector Double

data MatParam = MatParam Color
    deriving (Show, Read, Eq)

data Material = Material {ka :: MatParam,
                          ks :: MatParam,
                          kd :: MatParam,
                          ke :: MatParam,
                          kr :: MatParam,
                          kt :: MatParam,
                          shininess :: Double}
    deriving (Show, Read, Eq)

data Light =    PointLight       {lColor :: Color,
                                  lPosition :: Vector Double,
                                  constantTerm :: Double,
                                  linearTerm :: Double,
                                  quadraticTerm :: Double}
            |   DirectionalLight {lColor :: Color,
                                  lOrientation :: Vector Double}
    deriving (Show, Read, Eq)

type Lights = [Light]

lightColor :: Light -> Color
lightColor (PointLight c _ _ _ _) = c
lightColor (DirectionalLight c _) = c

getDirection :: Light -> Vector Double -> Vector Double
getDirection (PointLight _ pos _ _ _) p = normalize (pos - p)
getDirection (DirectionalLight _ o) _ = -o

distanceAttenuation :: Light -> Vector Double -> Double
distanceAttenuation (PointLight _ pos c b a) p = min 1.0 atten
    where atten = (c * dist * dist + b * dist + a) ** (-1.0)
          dist = norm_2 (p - pos)
distanceAttenuation (DirectionalLight _ _) p = 1.0

shadowAttenuation :: Scene -> Light -> Vector Double -> Double
shadowAttenuation scene (PointLight _ pos _ _ _) p =
    let shadowDir = normalize (pos - p) in
    let ray = Ray p shadowDir in
    let isect = intersect scene ray in
    case isect of
        Nothing -> 1
        Just (Isect t _ _ _) ->
            if (t < 0) || (t > norm_2 (pos - p)) then
                1
            else
                0
shadowAttenuation scene dl@(DirectionalLight _ _) p =
    let shadowDir = normalize (getDirection dl p) in
    let ray = Ray p shadowDir in
    let isect = intersect scene ray in
    case isect of
        Nothing -> 1
        Just (Isect t _ _ _) ->
            if t < 0 then
                1
            else
                0

data Geometry = Sphere {sphereCenter :: Vector Double,
                        sphereRadius :: Double,
                        sphereMaterial :: Material}
    deriving (Show, Read, Eq)

type Objects = [Geometry]

data Scene = Scene {scLights :: Lights,
                    scObjects :: Objects,
                    scAmbient :: Color}
    deriving (Show, Read, Eq)

data Camera = Camera {eye :: Vector Double,
                      viewDir :: Vector Double,
                      upDir :: Vector Double,
                      aspect :: Double,
                      fov :: Double}
    deriving (Show, Read, Eq)

rayThrough :: Camera -> Double -> Double -> Ray
rayThrough (Camera eye viewDir upDir aspect fov) xx yy =
    let z = -viewDir in
    let y = upDir in
    let x = cross y z in
    let m = fromRows [x,y,z] in
    let fov' = fov * pi / 180.0 in
    let normalizedHeight = 2 * tan (fov' / 2.0) in
    let u = (flatten (m <> (fromColumns [(fromList [1,0,0])]))) * (fromList [normalizedHeight * aspect]) in
    let v = (flatten (m <> (fromColumns [(fromList [0,1,0])]))) * (fromList [normalizedHeight]) in
    let look = (flatten (m <> (fromColumns [(fromList [0,0,-1])]))) in
    let x' = fromList [xx - 0.5] in
    let y' = fromList [yy - 0.5] in
    let dir = normalize (look + x' * u + y' * v) in
    Ray eye dir

intersectObj :: Ray -> Geometry -> Maybe Isect
intersectObj r@(Ray rayPos rayDir) s@(Sphere center radius material) =
    let v = center - rayPos in
    let b = v <.> rayDir in
    let discriminant = b * b - v <.> v + radius * radius in
    if discriminant < 0.0 then
        Nothing
    else
        let det = sqrt discriminant in
        let (t1,t2) = (b - det, b + det) in
        if t2 < 0.0 then
            Nothing
        else if t1 > 0.0 then
            pure (Isect t1 (normalize (at t1 r)) s material)
        else
            pure (Isect t2 (normalize (at t2 r)) s material)

intersect :: Scene -> Ray -> Maybe Isect
intersect scene@(Scene _ objects _) r =
    let isects = filter (/= Nothing) . map (intersectObj r) $ objects in
    if length isects == 0 then
        Nothing
    else
        minimum isects

perLightDiffuse :: Light -> Vector Double -> Vector Double -> Vector Double -> Double
perLightDiffuse light isectPoint rayDir normal =
    let lightDir = getDirection light isectPoint in
    let ddot = lightDir <.> (if (rayDir <.> normal > 0) then -normal else normal) in
    let dmax = max ddot 0.0 in
    dmax

perLightSpecular :: Light -> Vector Double -> Vector Double -> Vector Double -> Double
perLightSpecular light isectPoint rayDir normal =
    let lightDir = getDirection light isectPoint in
    let reflect = lightDir - (fromList [2 * (lightDir <.> normal)]) * normal in
    if norm_2 reflect == 0 then
        0
    else
        let reflectNorm = normalize reflect in
        let sdot = rayDir <.> reflect in
        let smax = max sdot 0.0 in
        smax

combineLightTerms :: Scene -> Vector Double -> Vector Double -> Vector Double -> Double -> Light -> Color
combineLightTerms scene isectPoint rayDir normal alpha light =
    let diffuse = perLightDiffuse light isectPoint rayDir normal in
    let specular = perLightSpecular light isectPoint rayDir normal in
    let color = lightColor light in
    let atten = (distanceAttenuation light isectPoint) * (shadowAttenuation scene light isectPoint) in
    color * (fromList [(diffuse + specular) * atten])

shade :: Scene -> Ray -> Isect -> Color
shade scene@(Scene lights objects ambient) r@(Ray p d) i@(Isect t normal _ material) =
    let (Material kaM ksM kdM keM krM ktM shiny) = material in
    let (MatParam ke) = keM in
    let emissive = ke in
    let (MatParam ka) = kaM in
    let amb = ka * ambient in
    let isectPoint = at t r in
    let alpha = shiny in
    emissive + amb + (sum . map (combineLightTerms scene isectPoint d normal alpha) $ lights)

data Screen = Screen {screenWidth :: Int,
                      screenHeight :: Int}
    deriving (Show, Read, Eq)

clamp :: Double -> Double -> Vector Double -> Vector Double
clamp a b = fromList . map (\x -> max a (min b x)) . toList

trace :: Scene -> Camera -> Double -> Double -> Color
trace scene camera x y =
    let ray = rayThrough camera x y in
    let color = traceRay scene ray 0 0 in
    clamp 0 1 color

tracePixel :: Scene -> Screen -> Camera -> Int -> Int -> Color
tracePixel scene (Screen width height) camera i j =
    let x = (fromIntegral i :: Double) / (fromIntegral width :: Double) in
    let y = (fromIntegral j :: Double) / (fromIntegral height :: Double) in
    let color = trace scene camera x y in
    (fromList [255.0]) * color

reflection :: Scene -> Ray -> Isect -> Int -> Color
reflection scene ray@(Ray p d) isect@(Isect t n _ mat) depth =
    let reflectR = d - (fromList [2 * d <.> n]) * n in
    if reflectR == 0 then
        fromList ([0,0,0])
    else
        let isectPointApprox = at (t - eps) ray in
        let reflRay = Ray isectPointApprox (normalize reflectR) in
        let (Material _ _ _ _ (MatParam kr)  _ _) = mat in
        kr * (traceRay scene reflRay (depth - 1) 0)

traceRay :: Scene -> Ray -> Int -> Double -> Color
traceRay scene ray@(Ray p d) depth t =
    let mIsect = intersect scene ray in
    case mIsect of
        Nothing -> fromList[0,0,0]
        (Just isect) -> let colorC = shade scene ray isect in
                        if depth == 0 then
                            colorC
                        else
                            let colorRefl = reflection scene ray isect depth in
                            colorC + colorRefl

type Image = ([[Color]], Int, Int)

traceImage :: Scene -> Screen -> Camera -> Image
traceImage scene screen@(Screen width height) camera =
    let (w, h) = (width - 1, height - 1) in
    let img = [[tracePixel scene screen camera i j | j <- [0..h]] | i <- [0..w]] in
    (img, width, height)

rgbify :: [Double] -> G.Pixel G.RGB Double
rgbify vec = G.PixelRGB (r / 255.0) (g / 255.0) (b / 255.0)
    where r = vec !! 0
          g = vec !! 1
          b = vec !! 2

render :: Image -> G.Image G.VU G.RGB Double
render (image, width, height) = G.fromLists . format $ image
    where format = map (map rgbify) . map (map toList)